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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Hallo Peter,
Thank you for the advice, now I have to update my table:
SAS R
random statement random subj(program); random = ~ 1 |
Subj
-2*loglik 1420.8 1420.820
random effects
variance(Intercept) 9.6033 9.603331
variance(residual) 1.1969 1.196873
the first 3 fixed effects
intercept 83.0952 83.09524
ProgramCont -3.4952 -3.49524
ProgramRI -1.9702 -1.97024
... ... ...
Everything looks nice. Perhaps Douglas could update the help file in
SASmixed,
where I copied the misleading code!
Kind regards,
Dominik
> -----Original Message-----
> From: Peter Dalgaard BSA [mailto:p.dalgaard at biostat.ku.dk]
> Sent: mercredi, 9. octobre 2002 12:37
> To: Grathwohl,Dominik,LAUSANNE,NRC/NT
> Subject: Re: [R] proc mixed vs. lme
>
>
> "Grathwohl,Dominik,LAUSANNE,NRC/NT"
> <dominik.grathwohl at rdls.nestle.com> writes:
>
> > 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 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 )
>
> It helps considerably if you fit the same model! Try:
>
> fm1Weight <- lme( strength ~ factor(Program) * factor(Time),
> data = Weights, random = ~ 1 | Subj)
>
> (or change the variables from numeric to factor inside Weights).
>
> --
> O__ ---- Peter Dalgaard Blegdamsvej 3
> c/ /'_ --- Dept. of Biostatistics 2200 Cph. N
> (*) \(*) -- University of Copenhagen Denmark Ph:
> (+45) 35327918
> ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX:
> (+45) 35327907
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
"Grathwohl,Dominik,LAUSANNE,NRC/NT" <dominik.grathwohl at rdls.nestle.com> writes:> 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?Different contrasts. Try setting options(contrasts = c(contr.SAS, contr.poly)) and doing the analysis in R again. Note that all of these examples are available in the SASmixed package for R. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
"Grathwohl,Dominik,LAUSANNE,NRC/NT" <dominik.grathwohl at rdls.nestle.com> writes:> Hallo Peter, > > Thank you for the advice, now I have to update my table: > > SAS R > random statement random subj(program); random = ~ 1 | > Subj > -2*loglik 1420.8 1420.820 > random effects > variance(Intercept) 9.6033 9.603331 > variance(residual) 1.1969 1.196873 > the first 3 fixed effects > intercept 83.0952 83.09524 > ProgramCont -3.4952 -3.49524 > ProgramRI -1.9702 -1.97024 > ... ... ... > > Everything looks nice. Perhaps Douglas could update the help file in > SASmixed, > where I copied the misleading code!Umm, my version of that help file has \examples{ options( contrasts = c(unordered = "contr.SAS", ordered = "contr.poly")) data(Weights) fm1Weight <- lme( strength ~ Program * Time, data = Weights, random = ~ 1 | Subj) summary( fm1Weight ) # compare with output 3.1, p. 91 VarCorr( fm1Weight ) anova( fm1Weight ) fm2Weight <- update( fm1Weight, random = ~ Time | Subj ) anova( fm1Weight, fm2Weight ) summary( fm2Weight ) VarCorr( fm2Weight ) intervals( fm2Weight ) } Notice the first line. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._