[R] make an model object (e.g. nlme) available in a user defined function (xyplot related)
Jun Shen
jun.shen.ut at gmail.com
Mon Jul 12 17:16:34 CEST 2010
Dear Deepayan,
Thank you for taking the time to look into this issue.
I have a data object called "Data", please find it at the end of the
message. Then I can run the code below separately in the console.
#Construct the nlme object
mod.nlme<-nlme(RESP~E0+(Emax-E0)*CP**gamma/(EC50**gamma+CP**gamma),data=Data,method='REML',
fixed=E0+Emax+gamma+EC50~1,
random=EC50~1,
groups=~ID,
start=list(fixed=c(E0=1,Emax=100,gamma=1,EC50=50))
)
#Plotting
xyplot(RESP~CP,Data,
groups=ID,
panel=panel.superpose,
panel.groups=function(x,y,subscripts,...){
panel.xyplot(x,y,...)
subjectData<-Data[subscripts,]
ind.pred<-predict(mod.nlme,newdata=subjectData)
panel.xyplot(x,ind.pred,type='l',lty=2)
}
)
##################################################################
Then I constructed a test function to put the two tasks together and
it seems OK. Strangely I don't even need to print() the xyplot, it is
just automatically shown on the screen.
test.function<-function(Data=Data){
mod.nlme<-nlme(RESP~E0+(Emax-E0)*CP**gamma/(EC50**gamma+CP**gamma),data=Data,method='REML',
fixed=E0+Emax+gamma+EC50~1,
random=EC50~1,
groups=~ID,
start=list(fixed=c(E0=1,Emax=100,gamma=1,EC50=50))
)
xyplot(RESP~CP,data=Data,
groups=ID,
panel=panel.superpose,
panel.groups=function(x,y,subscripts,...){
panel.xyplot(x,y,...)
subjectData<-Data[subscripts,]
ind.pred<-predict(mod.nlme,newdata=subjectData)
panel.xyplot(x,ind.pred,type='l',lty=2)
}
)
}
####################################################################
Then I have my real function as follows: If I run the code as,
>compare.curves(Data=Data)
The analytical part is working but not the plotting part (Error using
packet 1, object 'model' not found)
==============================================================
compare.curves<-function(curve='ascending',Data=stop('A data object
must be specified'),parameter='EC50',random.pdDiag=FALSE,
start.values=c(Emax=100,E0=1,EC50=50,gamma=2),...){
if (curve=='ascending') model=as.formula('RESP ~
E0+(Emax-E0)*CP**gamma/(EC50**gamma+CP**gamma)')
if (curve=='descending') model=as.formula('RESP ~
E0+(Emax-E0)*(1-CP**gamma/(EC50**gamma+CP**gamma))')
mod.nlme<-nlme(model=model,data=Data,method='REML',
fixed=Emax+E0+EC50+gamma~1,
random= if (length(parameter)==1)
eval(substitute(variable~1,list(variable=as.name(parameter))))
else {
variable<-as.name(parameter[1])
for (i in 2:length(parameter)) variable<-
paste(variable,'+',as.name(parameter[i]))
formula<-as.formula(paste(variable,'~1'))
if (random.pdDiag) list(pdDiag(formula))
else formula
},
groups=~ID,
start=list(fixed=start.values)
)
mod.nlme.RSS<-sum(resid(mod.nlme)^2)
df.mod.nlme<-dim(Data)[1]-(4+length(parameter)) # 4 fixed effects
and plus the number of random effects
constrained.fit.parameters<-coef(mod.nlme)
mod.nls.ind<-lapply(split(Data,Data$ID),function(x){
nls(formula=model,data=x,start=start.values)
})
mod.nls.ind.RSS<-do.call(sum,lapply(mod.nls.ind,function(x)resid(x)^2))
df.mod.nls.ind<-dim(Data)[1]-4*length(unique(Data$ID))
ind.fit.parameters<-do.call(rbind,lapply(mod.nls.ind,coef))
F.statistic<-mod.nlme.RSS/mod.nls.ind.RSS
F.test.p.value<-pf(F.statistic,df.mod.nlme,df.mod.nls.ind,lower.tail=FALSE)
print(
xyplot(RESP~CP,data=Data,
groups=ID,
panel=panel.superpose,
panel.groups=function(x,y,subscripts,...){
panel.xyplot(x,y,...)
subjectData<-Data[subscripts,]
ind.pred<-predict(mod.nlme,newdata=subjectData)
panel.xyplot(x,ind.pred,type='l',lty=2)
}
)
)
return(list(F_test_statistic=F.statistic,F_test_p_value=F.test.p.value,
Individual_fit=ind.fit.parameters,Constrained_fit=constrained.fit.parameters))
}
=============================================================
The data object "Data"
structure(list(ID = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L), CP = c(1, 2, 3, 4.5, 5, 7.5, 11.25, 12, 18, 30,
45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5, 11.25, 12, 18,
30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5, 11.25, 12,
18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5, 11.25,
12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5, 7.5,
11.25, 12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5, 5,
7.5, 11.25, 12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3, 4.5,
5, 7.5, 11.25, 12, 18, 30, 45, 60, 90, 120, 150, 200, 1, 2, 3,
4.5, 5, 7.5, 11.25, 12, 18, 30, 45, 60, 90), RESP = c(3.19, 2.52,
2.89, 3.28, 3.82, 7.15, 11.2, 16.25, 30.32, 55.25, 73.56, 82.07,
89.08, 95.86, 97.97, 99.03, 3.49, 4.4, 3.54, 4.99, 3.81, 10.12,
21.59, 24.93, 40.18, 61.01, 78.65, 88.81, 93.1, 94.61, 98.83,
97.86, 0.42, 0, 2.58, 5.67, 3.64, 8.01, 12.75, 13.27, 24.65,
46.1, 65.16, 77.74, 87.99, 94.4, 96.05, 100.4, 2.43, 0, 6.32,
5.59, 8.48, 12.32, 26.4, 28.36, 43.38, 69.56, 82.53, 91.36, 95.37,
98.36, 98.66, 98.8, 5.16, 2, 5.65, 3.48, 5.78, 5.5, 11.55, 8.53,
18.02, 38.11, 58.93, 70.93, 85.62, 89.53, 96.19, 96.19, 2.76,
2.99, 3.75, 3.02, 5.44, 3.08, 8.31, 10.85, 13.79, 32.06, 50.22,
63.7, 81.34, 89.59, 93.06, 92.47, 3.32, 1.14, 2.43, 2.75, 3.02,
5.4, 8.49, 7.91, 15.17, 35.01, 53.91, 68.51, 83.12, 86.85, 92.17,
95.72, 3.58, 0.02, 3.69, 4.34, 6.32, 5.15, 9.7, 11.39, 23.38,
42.9, 61.91, 71.82, 87.83)), .Names = c("ID", "CP", "RESP"), class =
"data.frame", row.names = c(NA,
-125L))
==============================================================
On Mon, Jul 12, 2010 at 1:05 AM, Deepayan Sarkar
<deepayan.sarkar at gmail.com> wrote:
> On Mon, Jul 12, 2010 at 2:51 AM, Jun Shen <jun.shen.ut at gmail.com> wrote:
>> Dear all,
>>
>> When I construct an nlme model object by calling nlme(...)->mod.nlme,
>> this object can be used in xyplot(). Something like
>>
>> xyplot(x,y,......
>> ......
>> ind.predict<-predict(mod.nlme)
>> ......
>> ) is working fine in console environment.
>>
>> But the same structure is not working in a user defined function. It
>> seems the "mod.nlme" created in a user defined function can not be
>> called in xyplot(). Why is that? Appreciate any comment. (The error
>> message says " Error in using packet 1, object "model" not found")
>
> Quoting from the footer:
>
> PLEASE [...] provide commented, minimal, self-contained, reproducible code.
>
> -Deepayan
>
>
>> Thanks.
>>
>> Jun Shen
>>
>> ______________________________________________
>> 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.
>>
>
More information about the R-help
mailing list