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. According to the value of
WhichApproach (1, 2, or 3), it attempts one of the approaches that
does not work or the workaround using the text string.
> 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
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.
