[R] lme in R (WinXP) vs. Splus (HP UNIX)

[Greg Tarpinian]

R2.2 for WinXP, Splus 6.2.1 for HP 9000 Series, HP-UX 11.0.


I am trying to get a handle on why the same lme( ) code gives
such different answers.  My output makes me wonder if the 
fact that the UNIX box is 64 bits is the reason.  The estimated
random effects are identical, but the fixed effects are very
different.  Here is my R code and output, with some columns 
and rows deleted for space considerations:

FOO.lme1 <- lme(fixed = Y ~ X1*X2, 
random = ~ X2 | X3, data = FOO,
method = "REML",
control = list(maxIter = 200, msMaxIter = 200, 
tolerance = 1e-8, niterEM = 200,
msTol = 1e-8, msVerbose = TRUE, 
optimMethod = "Nelder-Mead"))

  0      203.991: 0.924323 -0.0884338 0.493897
  1      203.991: 0.924323 -0.0884338 0.493897

> summary(FOO.lme1)
Linear mixed-effects model fit by REML
 Data: FOO 
      AIC      BIC    logLik
  357.484 373.6868 -170.7420

Random effects:
 Formula: ~ X2 | X3
 Structure: General positive-definite, 
 Log-Cholesky parametrization
               StdDev   Corr  
(Intercept)    2.025868 (Intr)
X2             4.908517 -0.475
Residual       4.493171       

Fixed effects: Y ~ X1*X2
                              Value  p-value
(Intercept)                39.08067  0.0000
X1                         12.08700  0.0000
X2                          5.85417  0.1370
X1:X2                       0.80833  0.8299 

Number of Observations: 60
Number of Groups: 3


Here is my Splus code and output, similarly edited for
space considerations:

FOO.lme1 <- lme(fixed = Y ~ X1*X2, 
random = ~ X2 | X3, data = FOO,
method = "REML",
control = list(maxIter = 200, msMaxIter = 200, 
tolerance = 1e-8, niterEM = 200,
msTol = 1e-8, msVerbose = TRUE, 
optimMethod = "Nelder-Mead"))

Iteration:  0 ,  1  function calls, F=  205.3776 
Parameters: [1]  0.89525077  0.22549223 -0.05936197

Iteration:  1 ,  2  function calls, F=  205.3776 
Parameters: [1]  0.89525078  0.22549224 -0.05936187

> summary(FOO.lme1)
Linear mixed-effects model fit by REML
 Data: FOO 
       AIC      BIC    logLik 
  360.2566 376.4594 -172.1283

Random effects:
 Formula:  ~ X2 | X3
 Structure: General positive-definite
                 StdDev   Corr 
   (Intercept) 2.025861 (Inter
            X2 4.908618 -0.475
      Residual 4.493172       

Fixed effects: Y ~ X1* X2
                        Value  p-value 
         (Intercept) 45.12417  <.0001
               X1     6.04350  <.0001
               X2     6.25833  0.0709
            X1:X2     0.40417  0.8299

Number of Observations: 60
Number of Groups: 3


I have examined the help.start() HTML pages and the lmeScale( )
function looks helpful, but I don't know how to use it to tell
R to use a particular set of starting values for the params.
It would have been instructive to see what happens if the same
starting values are specified for both Splus and R....

Below is a sample dataset that is fairly close to the data
that I have actually been using (not specifically the data
used to obtain the output above). 

Any help would be greatly appreciated.


Kind Regards,

     Greg




Sample Data:
A,X3,X2,X1,Y
1,A,1.6,A,43
1,A,1.6,B,56
2,A,1.6,A,43
2,A,1.6,B,50
3,A,1.6,A,46
3,A,1.6,B,64
4,A,2.0,A,42
4,A,2.0,B,42
5,A,2.0,A,36
5,A,2.0,B,44
6,A,2.0,A,31
6,A,2.0,B,44
7,A,2.0,A,39
7,A,2.0,B,46
8,A,2.4,A,47
8,A,2.4,B,52
9,A,2.4,A,50
9,A,2.4,B,52
10,A,2.4,A,43
10,A,2.4,B,56
11,B,1.6,A,35
11,B,1.6,B,48
12,B,1.6,A,32
12,B,1.6,B,45
13,B,1.6,A,37
13,B,1.6,B,50
14,B,2.0,A,40
14,B,2.0,B,55
15,B,2.0,A,33
15,B,2.0,B,46
16,B,2.0,A,38
16,B,2.0,B,53
17,B,2.0,A,42
17,B,2.0,B,55
18,B,2.4,A,39
18,B,2.4,B,56
19,B,2.4,A,43
19,B,2.4,B,56
20,B,2.4,A,45
20,B,2.4,B,62
21,C,1.6,A,33
21,C,1.6,B,41
22,C,1.6,A,33
22,C,1.6,B,44
23,C,1.6,A,31
23,C,1.6,B,46
24,C,2.0,A,38
24,C,2.0,B,48
25,C,2.0,A,34
25,C,2.0,B,45
26,C,2.0,A,36
26,C,2.0,B,51
27,C,2.0,A,40
27,C,2.0,B,57
28,C,2.4,A,39
28,C,2.4,B,53
29,C,2.4,A,34
29,C,2.4,B,51
30,C,2.4,A,36
30,C,2.4,B,52