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

[Duncan Murdoch]

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