[R] lme problem
Henric Nilsson
henric.nilsson at statisticon.se
Wed Apr 13 18:16:52 CEST 2005
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.694
Coefficients:
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 ` ' 1
Residual 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 mzarkov at eunet.yu
> ******************************************************
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> 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 at 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
More information about the R-help
mailing list