[R] mixed model for Splus and R
Douglas Bates
bates at stat.wisc.edu
Tue Nov 18 15:00:17 CET 2003
Lei Liu <liulei at l.imap.itd.umich.edu> writes:
> I try to compare the mixed model package "lme" by Splus and R. I used the
> dataset "Ovary" and the following code assuming AR(1) model for the error term:
>
> lme(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), data=Ovary, random =
> pdDiag(~sin(2*pi*Time) ) , correlation=corAR1() )
>
> But I got different results! And then I used a simpler model:
>
> lme(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), data=Ovary, random =
> pdDiag(~sin(2*pi*Time)) )
>
> This time I got the same output. I wonder the reason for the inconsistence
> between R and Splus. I attached the results as follows. Can someone tell me
> how to do it correctly in R? Thanks!
>
> Lei
>
> Result from R:
>
> ******************************************************************
> Linear mixed-effects model fit by REML
> Data: Ovary
> Log-restricted-likelihood: -775.2224
> Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time)
> (Intercept) sin(2 * pi * Time) cos(2 * pi * Time)
> 12.1895808 -2.9473432 -0.8807113
>
> Random effects:
> Formula: ~sin(2 * pi * Time) | Mare
> Structure: Diagonal
> (Intercept) sin(2 * pi * Time) Residual
> StdDev: 2.807293 0.03630784 3.665217
>
> Correlation Structure: AR(1)
> Formula: ~1 | Mare
> Parameter estimate(s):
> Phi
> 0.6073908
> Number of Observations: 308
> Number of Groups: 11
>
> *****************************************************************************************
>
> Result from Splus:
>
> *****************************************************************************************
> Linear mixed-effects model fit by REML
> Data: Ovary
> Log-restricted-likelihood: -774.724
> Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time)
> (Intercept) sin(2 * pi * Time) cos(2 * pi * Time)
> 12.18809 -2.985282 -0.8777629
>
> Random effects:
> Formula: ~ sin(2 * pi * Time) | Mare
> Structure: Diagonal
> (Intercept) sin(2 * pi * Time) Residual
> StdDev: 2.858478 1.25802 3.507097
>
> Correlation Structure: AR(1)
> Formula: ~ 1 | Mare
> Parameter estimate(s):
> Phi
> 0.5722019
> Number of Observations: 308
> Number of Groups: 11
> ******************************************************************************************
S-PLUS and R use different optimizer code so the results will not be
identical. You may want to check with intervals() to see how
well-defined the optima are. Also try setting
control=list(msVerbose=TRUE) in the call to lme to see how the actual
optimization is behaving.
The log-likelihood in complex lme and nlme models can be a difficult
function to optimize. I am aware that the results from R's optimizer(s)
are not always as good as those from the optimizer used in S-PLUS and
we are working on rectifying that problem.
More information about the R-help
mailing list