[R] scope, lme, ns, nlme, splines
Duncan Murdoch
murdoch.duncan at gmail.com
Wed Dec 19 23:10:35 CET 2012
On 12-12-19 5:02 PM, Jacob Wegelin wrote:
>
> Your solution works beautifully for predict.lme in the same data that
> were used to compute the model. But what if I want to compute the
> population fitted values in newdata?
I didn't offer a solution, I offered a workaround for a bug. The
solution is to fix the bug. I don't know how to do that.
Duncan Murdoch
In the expanded example below,
> predict.lm is able to find the object called "fixed". But predict.lme is
> unable to find it and returns an error.
>
> library("nlme")
> library(splines)
> library(ggplot2)
> set.seed(5)
> junk<-data.frame(x=rnorm(100))
> junk$y<-junk$x + rnorm(nrow(junk))
> junk$idvar<- factor(sample(LETTERS[1:10], size=nrow(junk), replace=TRUE))
> PRETTYNEWDATA<-data.frame(x=seq(-2, 2, length=20))
>
> fitfunction<-function(splineDF=1) {
> fixed <- eval(
> substitute(
> expr= y ~ ns(x, df=splineDF)
> , env= list(splineDF=splineDF)
> )
> )
> thislm<- lm( fixed , data= junk)
> junk$fitlm<-predict(thislm)
> thislme<-lme( fixed= fixed , random= ~ 1 | idvar , data= junk)
> junk$fitlme<-predict(thislme,level=0)
> print(
> ggplot( junk, aes(x, y))
> + geom_point(shape=1)
> + geom_line( aes(x, fitlme), color="red", lwd=2)
> + geom_line( aes(x, fitlm), color="blue", lty=3)
> + ggtitle(deparse(fixed))
> )
> PRETTYNEWDATA$fitlm<-predict(thislm,newdata=PRETTYNEWDATA)
> print(head(PRETTYNEWDATA))
> cat("about to use newdata with lme \n")
> PRETTYNEWDATA$fitlme<-predict(thislme,level=0, newdata=PRETTYNEWDATA)
> thislme
> }
>
> Jake Wegelin
>
> On Thu, 6 Dec 2012, Duncan Murdoch wrote:
>
>> On 06/12/2012 2:28 PM, Jacob Wegelin wrote:
>>> I want to fit a series of lme() regression models that differ only in the
>>> degrees of freedom of a ns() spline. I want to use a wrapper function to do
>>> this. The models will be of the form
>>>
>>> y ~ ns(x, df=splineDF)
>>>
>>> where splineDF is passed as an argument to a wrapper function.
>>>
>>> This works fine if the regression function is lm(). But with lme(),
>>> I get an error. fitfunction() below demonstrates this.
>>>
>>> Why?
>>
>> Presumably this is a bug in nlme. Formulas have associated environments. In
>> your example, that would be the local frame of the fitfunction, and splineDF
>> lives there. It looks as though lme is not looking in the formula
>> environment for the splineDF value.
>>
>> Another workaround besides the one you found is to use substitute() to put
>> the value of splineDF into the formula, i.e.
>>
>> fitfunction<-function(splineDF=1) {
>> junk$splineDF <- splineDF
>> thislm<-lm( y ~ ns(x, df=splineDF), data= junk)
>> thislm
>> cat("finished thislm\n")
>> fixed <- eval(substitute( y ~ ns(x, df=splineDF),
>> list(splineDF=splineDF)))
>> thislme<-lme(
>> fixed= fixed
>> , random= ~ 1 | idvar
>> , data= junk)
>> thislme
>> }
>>
>> It's not much better than your kludge in this example, but it might be more
>> convenient if you want to work with different formulas.
>>
>> Duncan Murdoch
>>
>>>
>>> KLUDGEfit() below provides a clumsy solution. It turns the lme()
>>> command, along with the appropriate value of splineDF, into a text
>>> string and then uses eval(parse(text=mystring)). But is there not a
>>> more elegant solution?
>>>
>>> set.seed(5)
>>> junk<-data.frame(x=rnorm(100))
>>> junk$y<-junk$x + rnorm(nrow(junk))
>>> junk$idvar<- factor(sample(LETTERS[1:10], size=nrow(junk), replace=TRUE))
>>> library("nlme")
>>> library(splines)
>>>
>>> fitfunction<-function(splineDF=1) {
>>> thislm<-lm( y ~ ns(x, df=splineDF), data= junk)
>>> thislm
>>> cat("finished thislm\n")
>>> thislme<-lme(
>>> fixed= y ~ ns(x, df=splineDF)
>>> , random= ~ 1 | idvar
>>> , data= junk)
>>> thislme
>>> }
>>> fitfunction()
>>>
>>> KLUDGEfit<-function(splineDF=2) {
>>> thislm<-lm( y ~ ns(x, df=splineDF), data= junk)
>>> thislm
>>> cat("finished thislm\n")
>>> mystring<-paste(
>>> "thislme<-lme( fixed= y ~ ns(x, df="
>>> , splineDF
>>> , "), random= ~ 1 | idvar , data= junk)"
>>> , sep="")
>>> eval(parse(text=mystring))
>>> thislme
>>> }
>>> KLUDGEfit()
>>>
>>> Thanks for any insight
>>>
>>> Jacob A. Wegelin
>>> Assistant Professor
>>> Department of Biostatistics
>>> Virginia Commonwealth University
>>> 830 E. Main St., Seventh Floor
>>> P. O. Box 980032
>>> Richmond VA 23298-0032
>>> U.S.A.
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
More information about the R-help
mailing list