Dear friends of lme,
After so many year with lme, I feel ashamed that I cannot get this to work.
Maybe it's a syntax problem, but possibly a lack of understanding.
We have growth curves of new dental bone that can well be modeled by a
linear growth curve, for two different treatments and several subjects as
random parameter. By definition, newbone is zero at t=0, so I tried to force
the curve through 0. In Pinheiro/Bates, this is done by including the -1
term, and it works well when treatment is not included (newbone~t-1), but
seems to have no effect in (newone ~ t*treat-1).
What's wrong? Some bracket missing? I tried a few variants.
Dieter Menne
#--------------------------------
library(nlme)
library(lattice)
# Generated data
set.seed(4711)
subject = as.factor(letters[1:5])
varslope = rnorm(length(subject),0,0.02)
cslope = c (0.1,0.15)
grd = expand.grid(t=seq(5,15,by=5),
subject=subject,treat=c("contr","test"))
grd$slope = varslope[grd$subject] + cslope[grd$treat]
grd$newbone = grd$slope*grd$t+rnorm(nrow(grd),0,0.2)
xyplot(newbone~t|treat,groups=subject,data=grd,
type="l",xlim=c(0,20),ylim=c(0,3))
# With intercept
grd.lme1 = lme(newbone~t*treat,data=grd,random=~1|subject)
grd$pred1 = predict(grd.lme1,level=0)
summary(grd.lme1)
# How go force intercept = 0 ???
grd.lme0 = lme(newbone~t*treat-1,data=grd,random=~1|subject)
grd$pred0 = predict(grd.lme0,level=0)
summary(grd.lme0)
# Gives true,
all.equal(grd$pred1,grd$pred0)
# Everything as expected without treat
grd.lme2 = lme(newbone~t,data=grd,random=~1|subject)
grd$pred2 = predict(grd.lme2,level=0)
summary(grd.lme2)
# Forced intercept = 0
grd.lme3 = lme(newbone~t-1,data=grd,random=~1|subject)
grd$pred3 = predict(grd.lme3,level=0)
summary(grd.lme3)
# As expected: not equal
all.equal(grd$pred2,grd$pred3)
#-------------------------------------------------------------------
R version 2.9.0 Under development (unstable) (2009-03-13 r48127)
i386-pc-mingw32
locale:
LC_COLLATE=German_Germany.1252;LC_CTYPE=German_Germany.1252;
LC_MONETARY=German_Germany.1252;LC_NUMERIC=C;LC_TIME=German_Germany.1252
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] lattice_0.17-20 nlme_3.1-90
loaded via a namespace (and not attached):
[1] grid_2.9.0 tools_2.9.0
>
______________________________________________
R-help@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.
