[R] lme( y ~ ns(x, df=splineDF)) error

[Jacob Wegelin]

I would like to fit regression models 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 two different errors, depending on how I handle splineDF
inside the wrapper function.

A workaround is to turn the lme() command, along with the appropriate
value of splineDF, into a text string and then use
eval(parse(text=mystring)).  Is there not a more elegant solution?

The function below demonstrates this.

> sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods
base

other attached packages:
[1] nlme_3.1-104  ggplot2_0.9.1

loaded via a namespace (and not attached):
  [1] colorspace_1.1-1   dichromat_1.2-4    digest_0.5.2
  grid_2.15.1        labeling_0.1
   [6] lattice_0.20-6     MASS_7.3-18        memoise_0.1
   munsell_0.3        plyr_1.7.1
   [11] proto_0.3-9.2      RColorBrewer_1.0-5 reshape2_1.2.1
   scales_0.2.1       stringr_0.6
   [16] tools_2.15.1

The following function can also be downloaded from
http://jacobwegelin.net/tmp/r-help/.

wrapper <-function(WhichApproach=1, splineDF=4 ){
#	Create toy dataset
nID<-25
IDnames<-LETTERS[1:nID]

longdat<-data.frame( x=rep(-10:10, nID) / 3 , ID= rep( IDnames , each=21)
)
set.seed(5)
longdat$x<-longdat$x + rnorm(nrow(longdat))/10
IDeffect<-rnorm(nID) * 20
names(IDeffect) <- IDnames
longdat$IDeffect<-IDeffect[as.character(longdat$ID)]
longdat$y<- (longdat$x
 			+ longdat$x^3
 			+ (longdat$x-1)^4  / 5
 			+ 1/(abs(longdat$x/50) + 0.02)
 			+ longdat$IDeffect
 			+ rnorm(1:nrow(longdat)) * 2
 			)
longdat<-longdat[order(longdat$x),]

library(splines)
#	Calling ns within lm works fine:
mylm<- lm( y ~ ns(x,df=splineDF), data=longdat)
longdat$lmfit<-predict(mylm)
library(ggplot2)
print(
 	ggplot(longdat, aes(x, y))
 	+ geom_point(shape=1)
 	+ geom_line(aes(x=x, y=lmfit), color="red")
)


cat("Enter to attempt lme.")
readline()
library(nlme)

if(WhichApproach==1) {
# returns: "Error in eval(expr, envir, enclos) : object 'splineDF' not found"
#
# First make sure splineDF does not exist in .GlobalEnv, else we would be in WhichApproach==2.
if(exists("splineDF", where=".GlobalEnv")) remove(list="splineDF", pos=".GlobalEnv")
mylme<-lme(fixed= y ~ ns(x, df=splineDF)
 			, random= ~ 1 | ID
 			, correlation = corAR1()
 			, data=longdat
 		)
} else if(WhichApproach==2) {
# returns: "Error in model.frame.default(formula = ~ID + y + x + splineDF, data = list( : 
# variable lengths differ (found for 'splineDF')"
assign(x="splineDF", value=splineDF, pos=".GlobalEnv")
mylme<-lme(fixed= y ~ ns(x, df=splineDF)
 			, random= ~ 1 | ID
 			, correlation = corAR1()
 			, data=longdat
 		)
} else if (WhichApproach==3) {

thecommandstring<- paste("mylme<-lme(fixed= y ~ ns(x, df=", splineDF, ") , random= ~ 1 | ID , correlation = corAR1() , data=longdat)")
thecommandstring
eval(parse(text=thecommandstring))

} else {
 	stop(paste("WhichApproach=", WhichApproach, " not valid."))
}

mylme
longdat$fullfit<-predict(mylme)

library(ggplot2)
print(
 	ggplot( longdat, aes(x,y))
 	+ geom_point(shape=1)
 	+ facet_wrap( ~ ID )
 	+ geom_line( aes(x, fullfit), color="red")
)

invisible(mylme)
}

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.