Skip to content

lme problem

2 messages · Milos Zarkovic, Henric Nilsson

Milos Zarkovic
# 
Sorry for the long letter!

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.

My R code is:

LME.1=lme(time~bstime:type+type-1,data=CCE,random=~1|type)

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

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
1 day later
Henric Nilsson
# 
Milos Zarkovic said the following on 2005-04-12 16:40:
The PROC MIXED code above doesn't correspond to the R code below.
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
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