Milos Zarkovic said the following on 2005-04-12 16:40:
I have recently started using R. For the start I have tried to repeat examples from Milliken & Johnson "Analysis of Messy Data - Analysis of Covariance", but I can not replicate it in R. The example is chocolate chip experiment. Response variable vas time to dissolve chocolate chip in seconds (time), covariate was time to dissolve butterscotch chip (bstime), and type was a type of chocolate chip. Problem is that I obtain different degrees of freedom compared to one in the book. Could it be sum of squares problem (type III vs. type I)? Milliken & Johnson use SAS for calculations and this is program the used:
proc mixed data=mmacov method=type3; class type; model time=type bstime*type/solution noint.
The PROC MIXED code above doesn't correspond to the R code below.
My R code is:
LME.1=lme(time~bstime:type+type-1,data=CCE,random=~1|type)
In your `lme' call, a random effect for each of the levels of `type' has been added to the model.
Since the analysis performed by PROC MIXED doesn't have any random effects it can be reproduced in R using the `lm' function. The results below match those of Milliken & Johnson p. 49 (using PROC MIXED) and the results on p. 43 (using PROC GLM).
> fit <- lm(time ~ bstime:type + type - 1, data = CCE) > summary(fit)
Call: lm(formula = time ~ bstime:type + type - 1, data = CCE)
Residuals:
Min 1Q Median 3Q Max
-16.982 -3.196 -0.250 1.400 21.694Coefficients:
Estimate Std. Error t value Pr(>|t|)
typeBlue M&M 17.9744 16.1923 1.110 0.27845
typeButton 21.5719 10.7832 2.001 0.05738 .
typeChoc Chip 16.9167 15.1673 1.115 0.27622
typeRed M&M 26.5760 13.1722 2.018 0.05545 .
typeSmall M&M 22.1977 29.0849 0.763 0.45310
typeSnow Cap 8.7000 9.4131 0.924 0.36495
bstime:typeBlue M&M 1.0641 0.6187 1.720 0.09887 .
bstime:typeButton 1.3352 0.3743 3.567 0.00164 **
bstime:typeChoc Chip 1.1667 0.7302 1.598 0.12373
bstime:typeRed M&M 0.5300 0.5564 0.953 0.35075
bstime:typeSmall M&M 0.1919 0.9881 0.194 0.84775
bstime:typeSnow Cap 0.9000 0.3999 2.250 0.03428 *
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1Residual standard error: 8.196 on 23 degrees of freedom Multiple R-Squared: 0.9774, Adjusted R-squared: 0.9656 F-statistic: 82.8 on 12 and 23 DF, p-value: 5.616e-16
and summary is:
Value Std.Error DF t-value p-value
typeBlue M&M 18.0 18.5 0 0.97 NaN
typeButton 21.6 14.1 0 1.53 NaN
typeChoc Chip 16.9 17.7 0 0.96 NaN
typeRed M&M 26.6 16.0 0 1.66 NaN
typeSmall M&M 22.2 30.5 0 0.73 NaN
typeSnow Cap 8.7 13.1 0 0.67 NaN
bstime:typeBlue M&M 1.1 0.6 24 1.72 0.098
bstime:typeButton 1.3 0.4 24 3.57 0.002
bstime:typeChoc Chip 1.2 0.7 24 1.60 0.123
bstime:typeRed M&M 0.5 0.6 24 0.95 0.350
bstime:typeSmall M&M 0.2 1.0 24 0.19 0.848
bstime:typeSnow Cap 0.9 0.4 24 2.25 0.034
However in Milliken & Johnson all df are 23. Values (estimates) are almost identical, but there are some small differences in SE and t.
Using
anova(LME.1)
I obtain
numDF denDF F-value p-value type 6 0 18.19 NaN bstime:type 6 24 4.04 0.0061
but in the book it is:
numDF denDF F-value p-value type 6 23 2.00 0.1075 bstime:type 6 23 4.04 0.0066
The tests reported by Milliken & Johnson are based on so called "Type III" sums of squares. If you want to reproduce these, try the `Anova' function in John Fox's indispensable `car' package.
> library(car)
> options(contrasts = c("contr.sum", "contr.poly"))
> Anova(fit, type = "III")
Anova Table (Type III tests)Response: time
Sum Sq Df F value Pr(>F)
type 805.13 6 1.9976 0.107510
bstime:type 1628.79 6 4.0412 0.006557 **
Residuals 1545.01 23
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
HTH, Henric
Data are at the end of the letter.
I am not sure what I did wrong.
Sincerely,
Milos Zarkovic
****************************************************** Milos Zarkovic MD, Ph.D. Associate Professor of Internal Medicine Institute of Endocrinology Dr Subotica 13 11000 Beograd Serbia
Tel +381-63-202-925 Fax +381-11-685-357
Email [EMAIL PROTECTED] ******************************************************
type,person,bstime,time Button,1,27,53 Choc Chip,2,17,36 Blue M&M,3,28,60 Blue M&M,4,30,45 Red M&M,5,20,30 Choc Chip,6,29,51 Small M&M,7,30,25 Button,8,16,47 Small M&M,9,32,25 Blue M&M,10,19,38 Blue M&M,11,33,48 Button,12,19,39 Snow Cap,13,15,20 Blue M&M,14,19,34 Choc Chip,15,20,40 Blue M&M,16,24,42 Snow Cap,17,21,29 Button,18,35,90 Red M&M,19,35,45 Small M&M,20,30,33 Button,21,34,65 Button,22,40,58 Small M&M,23,22,26 Snow Cap,24,16,23 Button,25,28,72 Blue M&M,26,25,48 Choc Chip,27,14,34 Button,28,23,45 Snow Cap,28,40,44 Blue M&M,30,28,48 Snow Cap,31,19,26 Snow Cap,32,21,29 Small M&M,33,32,30 Red M&M,34,16,32 Red M&M,35,19,47
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
