Doug,
I thought perhaps that you might be interested in the comparison of
lme to the results for the same models fitted by Richard Jones' carma
(I just wrote the R interface to his Fortran code). The code to run
the example from the lme help and for the equivalent with carma is in
the file below.
The two main differences in results are
1. the random coefficients covariance matrix is quite different
because lme takes the time origin as zero age whereas carma centres
times at the mean (here age). This arbitrariness of results is one of
the fundamental problems with random time coefficient models.
2. carma gives the true full likelihood whereas lme gives the so-called
restricted likelihood. In my opinion, substituting the REML estimate
obtained from the marginal or conditional likelihood back into the
full likelihood makes no sense.
Best wishes, Jim
------------------------------------------------------------------------
library(nlme)
library(growth)
# example from lme
data(Orthodont)
summary(fm1 <- lme(distance ~ age, data = Orthodont)) # random is ~ age
summary(fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1))
# set up data object
distance <- matrix(Orthodont[1],ncol=4,byrow=T)
age <- matrix(Orthodont[2],ncol=4,byrow=T)
sex <- Orthodont[4][seq(4,108,by=4)]
data <- rmna(restovec(distance,times=age),ccov=tcctomat(sex))
rm(Orthodont,age,sex,distance)
#model fm1
carma(data,torder=1,pre=rep(1,3),pos=rbind(c(1,1),c(2,2),c(1,2)))
# uncorrelated random coefficients
carma(data,torder=1,pre=rep(1,2),pos=rbind(c(1,1),c(2,2)))
#model fm2
carma(data,ccov=~sex,torder=1,pre=1,pos=c(1,1))
# random coefficient for age
carma(data,ccov=~sex,torder=1,pre=rep(1,2),pos=rbind(c(1,1),c(2,2)))
# quadratic in age
carma(data,ccov=~sex,torder=2,pre=1,pos=c(1,1))
# uncorrelated random coefficients for age
carma(data,ccov=~sex,torder=2,pre=rep(1,3),pos=rbind(c(1,1),c(2,2),c(3,3)))
# random intercept and AR(1)
carma(data,ccov=~sex,torder=1,pre=1,pos=c(1,1),arma=c(1,0,0),par=0.5)
# interaction between sex and time
carma(data,ccov=~sex,torder=1,inter=1,pre=1,pos=c(1,1))
# interaction between sex and quadratic time
carma(data,ccov=~sex,torder=2,inter=2,pre=1,pos=c(1,1))
# interaction between sex and linear time only
carma(data,ccov=~sex,torder=2,inter=1,pre=1,pos=c(1,1))
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
