Executive summary: How to obtain with lme the exact F-tests for the design
(C/Sujet)*A*B, where Sujet is a random factor.
Rod,
Thank you very much for your e-mail. To continue on this subject,
using anova(...,,type="marginal") with your models will give exactly
the same F-values and P-values as aov, and therefore the exact F-tests
(see example below). Another way to obtain the exact tests for A, C
and A:C is to average out on levels of factor B, and to run lme
(equivalent to aov response, see example below). Averaging over A
gives us B, C and B:C. Now the challenge is to obtain the exact test
for A:B and A:B:C. Does anyone has an idea on how to "average over
interaction". If the levels of A and B are 2 and 2, we can average
a1b1,a2b2 on one side and a1b2,a2b1 on the other side, but how to
generalise? Anyone?
Olivier
> ## average out on levels of factor B
> orientaAC <- reducefacdes("Reponse", c( "C", "A", "Sujet"),orienta)
> ## lme and aov give the correct F-values, dfs and p-values
> orient.fitAC <- lme(Reponse ~ A*C, data=orienta, random=~1|Sujet,
> data=orientaAC, method="REML")
> anova(orient.fitAC,type="marginal")
numDF denDF F-value p-value
(Intercept) 1 18 63.20697 <.0001
A 1 18 14.13650 0.0014
C 1 18 0.10342 0.7515
A:C 1 18 0.12901 0.7236
> summary (aov(formula = Reponse ~ A*C + Error(Sujet/(A)), data = orientaAC))
Error: Sujet
Df Sum of Sq Mean Sq F Value Pr(F)
C 1 435.60 435.600 0.1034248 0.7514624
Residuals 18 75811.61 4211.756
Error: A %in% Sujet
Df Sum of Sq Mean Sq F Value Pr(F)
A 1 15163.24 15163.24 14.13652 0.0014341
A:C 1 138.38 138.38 0.12901 0.7236356
Residuals 18 19307.31 1072.63
> reducefacdes
function(resp, fac, data, fct = mean)
{
ans <- tapply(data[, resp], interaction(data[, fac, drop = F], drop =
T), fct)
facdes <- unique(data[, fac, drop = F])
ans <- ans[order(names(ans))]
facdes <- facdes[order(as.character(interaction(facdes[, fac, drop =
F]))), , drop = F]
uni1 <- names(ans)
uni2 <- as.character(interaction(facdes[, fac, drop = F]))
if((length(uni1) != length(uni2)) | (!all(uni1 == uni2)))
stop("Problem. Contact Olivier.Renaud@pse.unige.ch with
example")
df <- cbind.data.frame(ans, facdes, row.names = NULL)
names(df)[1] <- resp
return(df)
}
> summary (aov(formula = Reponse ~ A*B*C + Error(Sujet/(A*B*C)), data =
> orienta))
Error: Sujet
Df Sum of Sq Mean Sq F Value Pr(F)
C 1 871.2 871.200 0.1034248 0.7514624
Residuals 18 151623.2 8423.513
Error: A %in% Sujet
Df Sum of Sq Mean Sq F Value Pr(F)
A 1 30326.47 30326.47 14.13652 0.0014341
A:C 1 276.77 276.77 0.12901 0.7236356
Residuals 18 38614.62 2145.26
Error: B %in% Sujet
Df Sum of Sq Mean Sq F Value Pr(F)
B 1 10296.72 10296.72 19.04355 0.0003742
B:C 1 62.66 62.66 0.11588 0.7374823
Residuals 18 9732.48 540.69
Error: A:B %in% Sujet
Df Sum of Sq Mean Sq F Value Pr(F)
A:B 1 1033.922 1033.922 2.069967 0.1673868
A:B:C 1 76.050 76.050 0.152256 0.7009675
Residuals 18 8990.768 499.487
> orient.fit2AB <- lme(Reponse ~ A*B*C,data=orienta,
+ random=list(Sujet = ~A, B= ~A, A= ~1))
> anova(orient.fit2AB)
numDF denDF F-value p-value
(Intercept) 1 36 58.26193 <.0001
A 1 36 14.13652 0.0006
B 1 18 18.94454 0.0004
C 1 18 0.12719 0.7255
A:B 1 36 2.06997 0.1589
A:C 1 36 0.12901 0.7216
B:C 1 18 0.11803 0.7352
A:B:C 1 36 0.15226 0.6987
> anova(orient.fit2AB,type="marginal")
numDF denDF F-value p-value
(Intercept) 1 36 63.20685 <.0001
A 1 36 14.13652 0.0006
B 1 18 19.04355 0.0004
C 1 18 0.10342 0.7515
A:B 1 36 2.06997 0.1589
A:C 1 36 0.12901 0.7216
B:C 1 18 0.11588 0.7375
A:B:C 1 36 0.15226 0.6987
At 02:13 27/05/03, Rod Ball wrote:
Cf. http://www.biostat.wustl.edu/archives/html/s-news/2002-05/msg00146.html
On Monday 26 May 2003 07:21 pm, Olivier Renaud wrote:
> Hi,
> I had a few responses, but none of them gave the "Psycologist's / SPSS"
> model. After having sent another mail to the authors of lme and to the
> support of Insightful, I received the attached e-mail.
Your original post was for
>summary (aov(formula = Reponse ~ A*B*C + Error(Block/(A*B)), data = ori))
It is possible to get the `correct' degrees of freedom Olivier wants from LME,
for most terms at least. The key to this is to realise that LME uses the
lowest level of grouping in which a term is estimable. However with the model
form given by Bernhard the only level of grouping is 'Block' hence the
degrees of freedom are all 54. Hence LME is testing the terms within Blocks.
(`Sujet in the original post). The degrees of freedom LME gives are correct
for the test.
So we need a grouping structure appropriate for the test Olivier wants. LME
supports only nested grouping structures---it is not possible to get all the
information with one structure. But Sujet/A/B or Sujet/B/A are close. Use
the former for tests of A, the latter for tests of B. We fill in the missing
terms at each level of grouping giving the models:-
# Model with grouping Sujet/B/A : use for tests of B and B:C
orient.fit2AB <- lme(Reponse ~ A*B*C,data=orienta,
random=list(Sujet = ~A, B= ~1, A= ~1))
anova(orient.fit2AB)
# Model with grouping Sujet/A/B : use for tests of A
orient.fit2BA <- lme(Reponse ~ A*B*C,data=orienta,
random=list(Sujet = ~B, A= ~1, B= ~1))
anova(orient.fit2BA)
Both models have the same terms, and the same loglikelihood.
See results below. The answers vary very little between the 'correct' and
'incorrect' models, as far as the anova tests go.
Incidentally I found the following model gives the best fit and may be all
you need.
# best model
orienta.fit0.1a <- update(orienta.fit, random = ~1|Sujet/A)
I would focus on the parameter estimates and mean effects and standard errors
of differences rather than the hypothesis tests. My mean.effects.table.lme()
function calculates these. See results below.
Regards,
Rod Ball
--
Dr Roderick D. Ball, Statistician,
NZ Forest Research Institute, P.B. 3020, Rotorua, New Zealand.
email: rod.ball@forestresearch.co.nz
phone: (07)343-5899 Ext 5413 fax: (07)348-0952
#####################################################
Then we have the results viz:-
# Splus anova table
> summary (aov(formula = Reponse ~ A*B*C + Error(Sujet/(A*B)), data =
orienta))
Error: Sujet
Df Sum Sq Mean Sq F value Pr(>F)
C 1 871 871 0.1034 0.7515
Residuals 18 151623 8424
Error: Sujet:A
Df Sum Sq Mean Sq F value Pr(>F)
A 1 30326 30326 14.136 0.001434 **
A:C 1 277 277 0.129 0.723636
Residuals 18 38615 2145
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Error: Sujet:B
Df Sum Sq Mean Sq F value Pr(>F)
B 1 10296.7 10296.7 19.0436 0.0003742 ***
B:C 1 62.7 62.7 0.1159 0.7374823
Residuals 18 9732.5 540.7
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Error: Sujet:A:B
Df Sum Sq Mean Sq F value Pr(>F)
A:B 1 1033.9 1033.9 2.0700 0.1674
A:B:C 1 76.1 76.1 0.1523 0.7010
Residuals 18 8990.8 499.5
> # model with grouping Sujet/B/A : use for tests of B
> orient.fit2AB <- lme(Reponse ~ A*B*C,data=orienta,
+ random=list(Sujet = ~A, B= ~A, A= ~1))
> anova(orient.fit2AB)
numDF denDF F-value p-value
(Intercept) 1 36 58.26196 <.0001
A 1 36 14.13653 0.0006
B 1 18 18.94448 0.0004
C 1 18 0.12719 0.7255
A:B 1 36 2.06997 0.1589
A:C 1 36 0.12901 0.7216
B:C 1 18 0.11803 0.7352
A:B:C 1 36 0.15226 0.6987
> # Model with grouping Sujet/A/B : use for tests of A
> orient.fit2BA <- lme(Reponse ~ A*B*C,data=orienta,
+ random=list(Sujet = ~B, A= ~1, B= ~1))
> anova(orient.fit2BA)
numDF denDF F-value p-value
(Intercept) 1 36 55.79459 <.0001
A 1 18 14.13676 0.0014
B 1 36 19.04177 0.0001
C 1 18 0.13421 0.7184
A:B 1 36 2.07009 0.1589
A:C 1 18 0.12902 0.7236
B:C 1 36 0.11587 0.7355
A:B:C 1 36 0.15227 0.6987
# Summary of selected model and mean effects tables.
> summary(orienta.fit0.1a)
Linear mixed-effects model fit by REML
Data: orienta
AIC BIC logLik
770.6683 795.7116 -374.3341
Random effects:
Formula: ~1 | Sujet
(Intercept)
StdDev: 39.59699
Formula: ~1 | A %in% Sujet
(Intercept) Residual
StdDev: 28.51672 22.80576
Fixed effects: Reponse ~ A * B * C
Value Std.Error DF t-value p-value
(Intercept) 42.10 17.03299 36 2.4716741 0.0183
Aa2 47.90 16.32976 18 2.9332940 0.0089
Bb2 29.70 10.19905 36 2.9120370 0.0061
CH 10.14 24.08828 18 0.4209515 0.6788
Aa2:Bb2 -10.48 14.42363 36 -0.7265855 0.4722
Aa2:CH -3.54 23.09377 18 -0.1532881 0.8799
Bb2:CH 0.36 14.42363 36 0.0249590 0.9802
Aa2:Bb2:CH -7.80 20.39809 36 -0.3823887 0.7044
Correlation:
(Intr) Aa2 Bb2 CH A2:Bb2 Aa2:CH Bb2:CH
Aa2 -0.479
Bb2 -0.299 0.312
CH -0.707 0.339 0.212
Aa2:Bb2 0.212 -0.442 -0.707 -0.150
Aa2:CH 0.339 -0.707 -0.221 -0.479 0.312
Bb2:CH 0.212 -0.221 -0.707 -0.299 0.500 0.312
Aa2:Bb2:CH -0.150 0.312 0.500 0.212 -0.707 -0.442 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.81651686 -0.48369671 -0.04471886 0.50940184 2.27429483
Number of Observations: 80
Number of Groups:
Sujet A %in% Sujet
20 40
> intervals(orienta.fit0.1a)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 7.555495 42.10 76.64450
Aa2 13.592438 47.90 82.20756
Bb2 9.015375 29.70 50.38462
CH -40.467610 10.14 60.74761
Aa2:Bb2 -39.732477 -10.48 18.77248
Aa2:CH -52.058219 -3.54 44.97822
Bb2:CH -28.892477 0.36 29.61248
Aa2:Bb2:CH -49.169249 -7.80 33.56925
Random Effects:
Level: Sujet
lower est. upper
sd((Intercept)) 25.18679 39.59699 62.25174
Level: A
lower est. upper
sd((Intercept)) 18.41118 28.51672 44.16899
Within-group standard error:
lower est. upper
18.10223 22.80576 28.73141
> mean.effects.table.lme(orienta.fit0.1a,A,table.name="A", data=orienta)
Table of A
*** mean effects ***
a1 a2
62.11 101.0
rep 40.00 40.0
*** seds ***
min av max
10.3598 10.3598 10.3598
> mean.effects.table.lme(orienta.fit0.1a,B,table.name="B", data=orienta)
Table of B
*** mean effects ***
b1 b2
70.23 92.92
rep 40.00 40.00
*** seds ***
min av max
5.099523 5.099523 5.099523
> mean.effects.table.lme(orienta.fit0.1a,C,table.name="C", data=orienta)
Table of C
*** mean effects ***
F H
78.28 84.88
rep 40.00 40.00
*** seds ***
min av max
20.51609 20.51609 20.51609
###########################################################33
> From: Bernhard Brabec <bernhard.brabec@insightful.com>
>
> >Subject: RE: Probably badly called "repeated measures" with lme
> > {[INBOX#34498]} To: Olivier.Renaud@pse.unige.ch
> >Cc: Geoffrey Pegler <geoff.pegler@predict.ch>, Shelp
> > <Shelp@uk.insightful.com>
> >
> >Dear Mr. Renaud,
> >
> >I have received the following answer to your question from our technical
> >support in Seattle:
> >
> >---
> >I apologize for the delay in responding.
> >
> >The correct model formula is:
> >> lme(fixed = Reponse ~ A*B*C,
> >
> > random = list(Block = pdBlocked(list(pdIdent(~1),
> > pdIdent(~B-1),
> > pdIdent(~C-1)))),
> > data = ori)
> >
> >but unfortunately there is no way to get the repeated measures
> >denominator degrees of freedom that the user wants.
> >This is one of the limitations of the lme() function. I will contact
> >the author of the lme() function and let him know that your
> >user feels this is a worthy enhancement request. I will also
> >pass along the users original E-mail..
> >---
> >
> >I hope this helps.
> >
> >Best regards
> > Bernhard Brabec, Insightful Switzerland
>
|