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

Greg Tarpinian sasprog474474 at yahoo.com
Tue Jan 31 02:30:15 CET 2006

```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,

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,

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

```