[R] mixed model for Splus and R

[Douglas Bates]

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.