[R] difference between splus and R

[Prof Brian D Ripley]

On Thu, 6 Apr 2000, Faheem Mitha wrote:

> I'm running splus 5 on a solaris platform remotely, and running R on linux
> on my home machine.
> 
> On Splus I have the two random effects models:
> 
> ran1 <- lme(fixed = earning ~ edu + job.pres + age,
>             random = ~ 1 + job.pres + age,cluster = ~ clus, data =
> labor.df )
>  
> ran2 <- lme(fixed = earning ~ edu + job.pres + age,
>             random = ~ 1 + job.pres, cluster = ~ clus, data = labor.df )

Might I suggest you install nlme 3.x on Splus5 too? (nlme.stat.wisc.edu)
Then you won't have to use two different syntaxes.

> anova(ran2,ran1) gives me 
> 
>      Model Df    AIC    BIC  Loglik    Test Lik.Ratio P value
> ran2     1  8 3371.5 3406.8 -1677.7
> ran1     2 11 3378.7 3427.2 -1678.3 1 vs. 2     1.189 0.75565
> 
> On R I have the random models 
> 
> ran1 <- lme(fixed = earning ~ edu + job.pres + age,
>      random = ~ 1 + job.pres + age | clus, data = labor.df )
> 
> ran2 <- lme(fixed = earning ~ edu + job.pres + age,
>             random = ~ 1 + job.pres | clus, data = labor.df )

> anova(ran2,ran1) gives me 
> 
>     Model df      AIC      BIC    logLik   Test    L.Ratio p-value
> ran2     1  8 3371.489 3406.705 -1677.745                          
> ran1     2 11 3377.515 3425.936 -1677.757 1 vs 2 0.02544750  0.9989
> 
> The models are supposed to be identical, and my understanding of the L
> ratio and the p value is that it the values corresponding to the null
> hypothesis that the smaller model is true ie. that the random effect due
> to age is zero. So a large p value in both cases corresponds to strong
> evidence in favour of the p value. 
> 
> My understanding is that both R and Splus are doing exactly this. So why
> are they returning different value. Are the models somehow different?
> 
> Another possibility is that one is using ordinary likelihood and the other
> is using REML. I see from the R documentation that REML is indeed used
> here and I thought the same was true of Splus.

(The fits for ran2 give the same statistics, so look both to be REML.) 
You should not be using anova on lme models fitted with REML. Although in
this case they are using the same fixed-effects model and so are on
comparable data, the supporting theory is for ML fits only, AFAIK.

> Also, I did summary(ran1) for both splus and R and the results which
> similar, did differ, for example in the estimation of the fixed effects.
> 
> Can anyone shed light on this? Am I missing something obvious?

The R results look like the problem, and suggest a possible optimization
failure. Please do try nlme3.x on S-PLUS to give an additional data point.
When nlme_3.1.5 is available for R, you might like to try that.


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._