Lei Liu <[EMAIL PROTECTED]> 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. ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
