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.
