[R] scope, lme, ns, nlme, splines

[Jacob Wegelin]

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? 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.
>
>