Just two observations (sorry, no solution):

1.  If one plugs in the covariance estimates into SAS as starting values, 
PROC MIXED says
    that these values have a lower likelihood (and therefore higher -2 Log 
L)
    than the values that PROC MIXED ends up with.

    SAS Code:

      proc mixed data=sasdata method=REML itdetails; 
      class Subject Day; 
      model Response = Day / solution; 
      parms (0.2782600) (0.06330649);
      repeated / subject = Subject type=ar(1); 
 
    (Note that 0.06330649 above is the square of 0.2603770 (standard 
deviation from lme's Random Effects section.)

    Part of SAS output:


                                   Parameter Search

               CovP1       CovP2    Variance            Res Log Like    -2 
Res Log Like

              0.2783     0.06331     0.06672                 -5.8696  
11.7392


                                         Iteration History

           CovP1       CovP2    Iteration    Evaluations    -2 Res Log 
Like       Criterion

          0.3411     0.06871            1              2 11.63183287 
0.00000544
          0.3389     0.06862            2              1 11.63168429 
0.00000000

    So PROC MIXED ends up with the same solution even with starting with 
the solution from lme.
    This suggests (to me) that maybe there's more to it than just PROC 
MIXED or lme converging 
    too soon as the log likelihood values have a different order between 
lme and PROC MIXED.

2.  I wonder why the distances in time are treated as if they were at 
equal intervals
    when they are 1, 2, 4, and 6.  The AR(1) covariance model in PROC 
MIXED essentially assumes that
    the distances are equivalent to the ranks of the observations (getting 
you the distances
    1, 2, 3, and 4).  While corAR1 in lme allows unequal (but integer) 
distances, your code has

        as.numeric(Time)

    which gives you the values 1, 2, 3, and 4 for the "time" intervals.

    However, if Day is used in place of as.numeric(Time) in lme and if 
"type=SP(POW)(d)" is used in
    place of "type=AR(1)" in PROC MIXED (and d is a duplicate of Day - 
which is necessary
    to create because the PROC MIXED statements have already declared Day 
as a
    CLASS variable), then the time values become 1, 2, 4, and 6 and the 
PROC MIXED and lme
    estimates for all of the parameters are nearly identical.  (And 
strangely enough, using SP(POW)(t)
    where t takes on the values 1, 2, 3, and 4 in place of 1, 2, 4, and 6, 
respectively, results in 
    the same values as using AR(1) in PROC MIXED.)



Jim

Jim Baldwin
Station Statistician
Pacific Southwest Research Station
USDA Forest Service
Albany, California
510-559-6332



Ben Bolker <bbolker@gmail.com> 
Sent by: r-sig-mixed-models-bounces@r-project.org
03/27/2011 09:25 AM

To
r-sig-mixed-models@r-project.org
cc

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






Simon Bate <simontbate@...> 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

_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


	[[alternative HTML version deleted]]


