[R-sig-ME] Repeated measures mixed model with AR(1) correlation structure in nlme vs SAS Proc Mixed

[Ben Bolker]

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