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
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 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]]