[R] proc mixed vs. lme

Grathwohl,Dominik,LAUSANNE,NRC/NT dominik.grathwohl at rdls.nestle.com
Wed Oct 9 10:28:44 CEST 2002


Dear All,

Comparing linear mixed effect models in SAS and R, I found the following
discrepancy:

                            SAS                           R
random statement            random subj(program);         random = ~ 1 |
Subj
-2*loglik                   1420.8                        1439.363
random effects
variance(Intercept)         9.6033                        9.604662
variance(residual)          1.1969                        1.187553
the first 3 fixed effects
intercept                   83.0952                       81.10544
ProgramCont                -3.4952                       -1.11526
ProgramRI                  -1.9702                       -1.04517
...                        ...                            ...

Can somebody explain me this different results?

thanks in advance,

Dominik



The two examples can be found in the internet:
http://ftp.sas.com/samples and http://cran.r-project.org

#My work around:

> R.version
         _              
platform i386-pc-mingw32
arch     i386           
os       mingw32        
system   i386, mingw32  
status                  
major    1              
minor    6.0            
year     2002           
month    10             
day      01             
language R

#The SAS code from "http://ftp.sas.com/samples/A55235":

data weights;
   input subj program$ s1 s2 s3 s4 s5 s6 s7;
   datalines;
 1      CONT      85    85    86    85    87    86    87
 2      CONT      80    79    79    78    78    79    78
 3      CONT      78    77    77    77    76    76    77
 4      CONT      84    84    85    84    83    84    85
 5      CONT      80    81    80    80    79    79    80
 6      CONT      76    78    77    78    78    77    74
 7      CONT      79    79    80    79    80    79    81
 8      CONT      76    76    76    75    75    74    74
 9      CONT      77    78    78    80    80    81    80
10      CONT      79    79    79    79    77    78    79
11      CONT      81    81    80    80    80    81    82
12      CONT      77    76    77    78    77    77    77
13      CONT      82    83    83    83    84    83    83
14      CONT      84    84    83    82    81    79    78
15      CONT      79    81    81    82    82    82    80
16      CONT      79    79    78    77    77    78    78
17      CONT      83    82    83    85    84    83    82
18      CONT      78    78    79    79    78    77    77
19      CONT      80    80    79    79    80    80    80
20      CONT      78    79    80    81    80    79    80
 1      RI        79    79    79    80    80    78    80
 2      RI        83    83    85    85    86    87    87
 3      RI        81    83    82    82    83    83    82
 4      RI        81    81    81    82    82    83    81
 5      RI        80    81    82    82    82    84    86
 6      RI        76    76    76    76    76    76    75
 7      RI        81    84    83    83    85    85    85
 8      RI        77    78    79    79    81    82    81
 9      RI        84    85    87    89    88    85    86
10      RI        74    75    78    78    79    78    78
11      RI        76    77    77    77    77    76    76
12      RI        84    84    86    85    86    86    86
13      RI        79    80    79    80    80    82    82
14      RI        78    78    77    76    75    75    76
15      RI        78    80    77    77    75    75    75
16      RI        84    85    85    85    85    83    82
 1      WI        84    85    84    83    83    83    84
 2      WI        74    75    75    76    75    76    76
 3      WI        83    84    82    81    83    83    82
 4      WI        86    87    87    87    87    87    86
 5      WI        82    83    84    85    84    85    86
 6      WI        79    80    79    79    80    79    80
 7      WI        79    79    79    81    81    83    83
 8      WI        87    89    91    90    91    92    92
 9      WI        81    81    81    82    82    83    83
10      WI        82    82    82    84    86    85    87
11      WI        79    79    80    81    81    81    81
12      WI        79    80    81    82    83    82    82
13      WI        83    84    84    84    84    83    83
14      WI        81    81    82    84    83    82    85
15      WI        78    78    79    79    78    79    79
16      WI        83    82    82    84    84    83    84
17      WI        80    79    79    81    80    80    80
18      WI        80    82    82    82    81    81    81
19      WI        85    86    87    86    86    86    86
20      WI        77    78    80    81    82    82    82
21      WI        80    81    80    81    81    82    83
;
               
/*---Data Set 3.2(b)---*/

data weight2; 
   set weights;
   time=1; strength=s1; output;
   time=2; strength=s2; output;
   time=3; strength=s3; output;
   time=4; strength=s4; output;
   time=5; strength=s5; output;
   time=6; strength=s6; output;
   time=7; strength=s7; output;
   keep subj program time strength;
run;
    
proc sort data=weight2; 
   by program time;
run;


/*---Data Set 3.2(c)---*/

proc means data=weight2 noprint; 
   by program time; 
   var strength;
   output out=avg mean=strength;
run;
               
               
/*---produces Output 3.1 on pages 91-92---*/

proc mixed data=weight2;
   class program subj time;
   model strength = program time program*time /s;
   random subj(program);
run;



#The SAS results:

                                      Covariance Parameter
                                            Estimates

                                    Cov Parm          Estimate

                                    subj(program)       9.6033
                                    Residual            1.1969

                                      The Mixed Procedure

                                         Fit Statistics

                              -2 Res Log Likelihood          1420.8
                              AIC (smaller is better)        1424.8
                              AICC (smaller is better)       1424.9
                              BIC (smaller is better)        1428.9


                                    Solution for Fixed Effects

                                                     Standard
      Effect          program    time    Estimate       Error      DF    t
Value    Pr > |t|

      Intercept                           83.0952      0.7171      54
115.87      <.0001
      program         CONT                -3.4952      1.0268      54
-3.40      0.0013
      program         RI                  -1.9702      1.0906      54
-1.81      0.0764
      program         WI                        0           .       .
.         .
      time                       1        -2.0476      0.3376     324
-6.06      <.0001
      time                       2        -1.4286      0.3376     324
-4.23      <.0001
      time                       3        -1.1905      0.3376     324
-3.53      0.0005
      time                       4        -0.5714      0.3376     324
-1.69      0.0915
      time                       5        -0.4762      0.3376     324
-1.41      0.1594
      time                       6        -0.3810      0.3376     324
-1.13      0.2600
      time                       7              0           .       .
.         .
      program*time    CONT       1         2.1976      0.4834     324
4.55      <.0001
      program*time    CONT       2         1.7786      0.4834     324
3.68      0.0003
      program*time    CONT       3         1.5905      0.4834     324
3.29      0.0011
      program*time    CONT       4         1.0214      0.4834     324
2.11      0.0354
      program*time    CONT       5         0.6762      0.4834     324
1.40      0.1628
      program*time    CONT       6         0.3810      0.4834     324
0.79      0.4312
      program*time    CONT       7              0           .       .
.         .
      program*time    RI         1         0.6101      0.5134     324
1.19      0.2356
      program*time    RI         2         0.8661      0.5134     324
1.69      0.0926
      program*time    RI         3         0.8780      0.5134     324
1.71      0.0882
      program*time    RI         4         0.4464      0.5134     324
0.87      0.3852
      program*time    RI         5         0.6012      0.5134     324
1.17      0.2425
      program*time    RI         6         0.3810      0.5134     324
0.74      0.4586
      program*time    RI         7              0           .       .
.         .
      program*time    WI         1              0           .       .
.         .
      program*time    WI         2              0           .       .
.         .
      program*time    WI         3              0           .       .
.         .
      program*time    WI         4              0           .       .
.         .
      program*time    WI         5              0           .       .
.         .
      program*time    WI         6              0           .       .
.         .
      program*time    WI         7              0           .       .
.         .



#The R code:

library(nlme)
library(SASmixed)
options( contrasts = c(unordered = "contr.SAS", ordered = contr.poly")) 
data(Weights) 
fm1Weight <- lme( strength ~ Program * Time, data = Weights, random = ~ 1 |
Subj) 
summary( fm1Weight )
VarCorr( fm1Weight ) 



#The R results:

> summary( fm1Weight )               
Linear mixed-effects model fit by REML
 Data: Weights 
       AIC      BIC    logLik
  1455.363 1487.153 -719.6815

Random effects:
 Formula: ~1 | Subj
        (Intercept) Residual
StdDev:    3.099139 1.089749

Fixed effects: strength ~ Program * Time 
                    Value Std.Error  DF   t-value p-value
(Intercept)      81.10544 0.7001315 339 115.84315  <.0001
ProgramCONT      -1.11526 1.0024358  54  -1.11255  0.2708
ProgramRI        -1.04517 1.0646834  54  -0.98168  0.3306
Time              0.15986 0.0224702 339   7.11447  <.0001
ProgramCONT:Time -0.18397 0.0321725 339  -5.71827  <.0001
ProgramRI:Time   -0.05495 0.0341703 339  -1.60822  0.1087
 Correlation: 
                 (Intr) PrCONT PrgrRI Time   PCONT:
ProgramCONT      -0.698                            
ProgramRI        -0.658  0.459                     
Time             -0.225  0.157  0.148              
ProgramCONT:Time  0.157 -0.225 -0.103 -0.698       
ProgramRI:Time    0.148 -0.103 -0.225 -0.658  0.459

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.11669139 -0.60963281  0.03963523  0.63074983  3.33520406 

Number of Observations: 399
Number of Groups: 57 

> VarCorr( fm1Weight )
Subj = pdLogChol(1) 
            Variance StdDev  
(Intercept) 9.604662 3.099139
Residual    1.187553 1.089749



Dominik Grathwohl 
Biostatistician 
Nestlé Research Center 
PO Box 44, CH-1000 Lausanne 26 
Phone: + 41 21 785 8034 
Fax: + 41 21 785 8556 
e-mail: dominik.grathwohl at rdls.nestle.com 

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list