[R-sig-ME] Repeated measures mixed model with AR(1) correlation structure in nlme vs SAS Proc Mixed
Ben Bolker
bbolker at gmail.com
Sun Mar 27 18:24:07 CEST 2011
Simon Bate <simontbate at ...> writes:
> Hi all, I don't know if anyone has any thoughts on this, I suspect
> it is a local convergence issue but am not sure. I have been trying
> to move from SAS Proc Mixed to R nlme and have an unusual result. I
> have several subjects measured at four time points. I want to model
> the within-subject correlation using an autoregressive
> structure. I've attached the R and SAS code I'm using along with the
> results from SAS. With R lme I get an estimate of the
> autoregressive parameter phi = 0.2782601, whereas SAS gives me an
> estimate of 0.3389 Intriguingly if I include a between subject
> factor, a covariate or delete one of the observations, then the
> results appear to agree. I'm surprised the seemingly simpler model
> if different between the two implementations whereas the more
> complex models agree. Any ideas would be most welcome! Regards,
> Simon
Hmmm. I do agree it could be a convergence problem of some sort.
I agree that it's worth trying to figure out what the difference
is, but I think it's going to be hard. I went through some basic
sanity checking with R (see below) and don't see an obvious problem.
library(nlme)
sasdata <- data.frame(Response=c(0.55,0.86,0.21,0.36,0.46,
0.32,0.11,0.24,0.36,0.29,0.48,0.93,
0.56,0.67,0.36,0.55,0.51,0.4,0.34,0.51,
1,0.61,0.65,0.41,0.99,0.86,0.64,0.86,0.31,
0.19,0.21,0.36,0.41,0.47,0.16,0.81,0.9,0.72,0.87,0.02),
Subject=c(1,1,1,1,2,2,2,2,3,3,3,3,
4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,
8,8,8,8,9,9,9,9,10,10,10,10),
## or: rep(1:10,each=4)
Day=c(1,2,4,6,1,2,4,6,1,2,4,6,1,2,4,6,
1,2,4,6,1,2,4,6,1,2,4,6,1,2,4,6,1,2,4,6,1,2,4,6))
## or: rep(c(1,2,4,6),4)
sasdata$Time<-factor(sasdata$Day)
fit0<-lme(Response~Time, random=~1|Subject,
data=sasdata, na.action = (na.omit), method = "REML")
plot(ACF(fit0),alpha=0.05,grid=TRUE)
## it actually looks like AR1 is not a very good model here
## anyway?
AR1<- update(fit0,
correlation=corAR1(form=~as.numeric(Time)|Subject,
fixed =FALSE))
plot(ACF(AR1),alpha=0.05,grid=TRUE)
## ACF(AR1) shows that the raw rho(1) is about 0.2
## identical (based on order of observations within groups)
AR1B<- update(fit0,
correlation=corAR1(form=~1|Subject,
fixed =FALSE))
## construct likelihood profile on rho
rhovec <- seq(-0.99,0.99,by=0.01)
Lvec <- sapply(rhovec,
function(r) {
logLik(update(fit0,
correlation=corAR1(form=~1|Subject,
value=r,
fixed =TRUE)))
})
plot(rhovec,-2*Lvec,type="b")
abline(h=-2*logLik(AR1)+3.84,lty=2) ## LRT cutoff
abline(v=c(-0.53,0.82),lty=2)
## CIs on rho from intervals() -- reasonable ballpark
## zoom in
plot(rhovec,-2*Lvec,type="b",ylim=c(11.5,15),xlim=c(-0.5,0.7))
abline(v=c(0.278,0.339)) ## SAS/R answers
## zoom in further
plot(rhovec,-2*Lvec,type="b",ylim=c(11.6,11.8),xlim=c(0.2,0.4),
las=1,bty="l")
abline(v=c(0.278,0.339),lty=2)
text(c(0.278,0.339),rep(par("usr")[4],2),c("lme","SAS"),
pos=3,xpd=NA)
abline(h=11.63168,lty=2)
text(par("usr")[2],11.63168,"SAS",pos=1,xpd=NA)
So ... I don't see anything obviously wrong or weird about
what R is doing. You can try
update(AR1,control=lmeControl(msVerbose=TRUE))
to see more detail on the optimization, but you will have
to dig into Pinheiro & Bates 2000 and/or the code and
documentation for details on the internal parameterizations
of rho in order to interpret the results ...
###
R Code:
SAS Code:
proc mixed;
class Subject Day;
model Response = Day / outp=pout;
repeated Day / subject = Subject type=AR(1);
run;
SAS Results:
Model Information
Data Set WORK.ALLDATA
Dependent Variable Response
Covariance Structure Autoregressive
Subject Effect Subject
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Between-Within
Class Level Information
Class Levels Values
Subject 10 1 10 2 3 4 5 6 7 8 9
Day 4 1 2 3 4
Dimensions
Covariance Parameters 2
Columns in X 5
Columns in Z 0
Subjects 10
Max Obs Per Subject 4
Number of Observations
Number of Observations Read 40
Number of Observations Used 40
Number of Observations Not Used 0
Iteration History
Iteration Evaluations -2 Res Log Like Criterion
0 1 14.67045653
1 2 11.63168913 0.00000018
2 1 11.63168429 0.00000000
Convergence criteria met.
Covariance Parameter Estimates
Cov Parm Subject Estimate
AR(1) Animal1 0.3389
Residual 0.06862
More information about the R-sig-mixed-models
mailing list