Mesures répétées avec R
C’est une de mes interrogations depuis que je pratique R : Comment prendre en compte les mesures répétées dans une analyse de variance et faire une analyse post-hoc satisfaisante? Je vais synthétiser toutes les solutions que j’ai pu collecter sur le web.
Le jeu de données a été produit par une étude sur la privation de sommeil (Gregory Belenky, Nancy J. Wesensten, David R. Thorne, Maria L. Thomas, Helen C. Sing, Daniel P. Redmond, Michael B. Russo and Thomas J. Balkin (2003) Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: a sleep dose-response study. Journal of Sleep Research 12, 1–12)
Y = temps de réaction moyen X = Le nombre de jour de privation Subject = …
library(lme4) # install.packages("lme4")## Loading required package: Matrix
## Loading required package: Rcpphelp(sleepstudy) # info sur les données
summary(sleepstudy)## Reaction Days Subject
## Min. :194 Min. :0.0 308 : 10
## 1st Qu.:255 1st Qu.:2.0 309 : 10
## Median :289 Median :4.5 310 : 10
## Mean :298 Mean :4.5 330 : 10
## 3rd Qu.:337 3rd Qu.:7.0 331 : 10
## Max. :466 Max. :9.0 332 : 10
## (Other):120sleepstudy$Days=as.factor(sleepstudy$Days)Graphique par sujet
library(ggplot2) # install.packages("ggplot2")
p<-ggplot(sleepstudy, aes(x=Days, y=Reaction,group=Subject))
p + geom_line(aes(colour=Subject)) + ylab("Average reaction time (ms)") + xlab("Days of Sleep deprivation")
Graphique moyen
d <- ggplot(sleepstudy, aes(x=Days, y=Reaction))
d + stat_summary(fun.data = "mean_cl_boot", colour = "red") + ylab("Average reaction time (ms)") + xlab("Days of Sleep deprivation")
Anova classique sans correction sans analyses post-hoc à l’exception du pairwise t.test car cela ne marche pas avec Error()
library(multcomp) # TukeyHSD## Loading required package: mvtnorm
## Loading required package: TH.dataaov.sleep=aov(Reaction~ Days + Error(Subject/Days), data= sleepstudy)
summary(aov.sleep)##
## Error: Subject
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 17 250618 14742
##
## Error: Subject:Days
## Df Sum Sq Mean Sq F value Pr(>F)
## Days 9 166235 18471 18.7 <2e-16 ***
## Residuals 153 151101 988
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#TukeyHSD(aov.sleep, "Days") Erreur dans UseMethod("TukeyHSD") : pas de méthode pour 'TukeyHSD' applicable pour un objet de classe "c('aovlist', 'listof')"
# summary(glht(aov.sleep,linfct=mcp(Days="Tukey")),adjusted(type="bonferroni"))
pairwise.t.test( x= sleepstudy$Reaction , g= sleepstudy$Days , p.adjust.method= "bonferroni") #post hoc qui marche##
## Pairwise comparisons using t tests with pooled SD
##
## data: sleepstudy$Reaction and sleepstudy$Days
##
## 0 1 2 3 4 5 6 7 8
## 1 1.00000 - - - - - - - -
## 2 1.00000 1.00000 - - - - - - -
## 3 1.00000 1.00000 1.00000 - - - - - -
## 4 1.00000 1.00000 1.00000 1.00000 - - - - -
## 5 0.07358 0.32729 0.38160 1.00000 1.00000 - - - -
## 6 0.03447 0.16688 0.19643 1.00000 1.00000 1.00000 - - -
## 7 0.00803 0.04506 0.05393 1.00000 1.00000 1.00000 1.00000 - -
## 8 8.5e-05 0.00069 0.00086 0.05123 0.15771 1.00000 1.00000 1.00000 -
## 9 1.3e-06 1.4e-05 1.8e-05 0.00203 0.00784 0.44066 0.81451 1.00000 1.00000
##
## P value adjustment method: bonferroniProblème = Les analyses post-hoc ne marchent pas pour les anovas avec une strate (within) “Error(Subject/Days)” à l’exception
1ère Solution
aov.sleep2=aov(Reaction ~ Days, data=sleepstudy)
summary(aov.sleep2)## Df Sum Sq Mean Sq F value Pr(>F)
## Days 9 166235 18471 7.82 1.3e-09 ***
## Residuals 170 401719 2363
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(aov.sleep2)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Reaction ~ Days, data = sleepstudy)
##
## $Days
## diff lwr upr p adj
## 1-0 7.8440 -44.10879 59.80 1.0000
## 2-0 8.7101 -43.24264 60.66 0.9999
## 3-0 26.3402 -25.61253 78.29 0.8339
## 4-0 31.9976 -19.95512 83.95 0.6179
## 5-0 51.8666 -0.08609 103.82 0.0508
## 6-0 55.5264 3.57371 107.48 0.0258
## 7-0 62.0988 10.14604 114.05 0.0067
## 8-0 79.9777 28.02496 131.93 0.0001
## 9-0 94.1994 42.24668 146.15 0.0000
## 2-1 0.8661 -51.08659 52.82 1.0000
## 3-1 18.4963 -33.45648 70.45 0.9796
## 4-1 24.1537 -27.79907 76.11 0.8942
## 5-1 44.0227 -7.93004 95.98 0.1754
## 6-1 47.6825 -4.27024 99.64 0.1021
## 7-1 54.2548 2.30209 106.21 0.0329
## 8-1 72.1338 20.18101 124.09 0.0006
## 9-1 86.3555 34.40273 138.31 0.0000
## 3-2 17.6301 -34.32262 69.58 0.9853
## 4-2 23.2875 -28.66521 75.24 0.9137
## 5-2 43.1566 -8.79618 95.11 0.1974
## 6-2 46.8164 -5.13638 98.77 0.1168
## 7-2 53.3887 1.43595 105.34 0.0386
## 8-2 71.2676 19.31487 123.22 0.0008
## 9-2 85.4893 33.53659 137.44 0.0000
## 4-3 5.6574 -46.29532 57.61 1.0000
## 5-3 25.5264 -26.42629 77.48 0.8582
## 6-3 29.1862 -22.76649 81.14 0.7336
## 7-3 35.7586 -16.19416 87.71 0.4563
## 8-3 53.6375 1.68476 105.59 0.0369
## 9-3 67.8592 15.90648 119.81 0.0018
## 5-4 19.8690 -32.08370 71.82 0.9672
## 6-4 23.5288 -28.42390 75.48 0.9085
## 7-4 30.1012 -21.85157 82.05 0.6973
## 8-4 47.9801 -3.97265 99.93 0.0974
## 9-4 62.2018 10.24906 114.15 0.0066
## 6-5 3.6598 -48.29294 55.61 1.0000
## 7-5 10.2321 -41.72061 62.18 0.9998
## 8-5 28.1111 -23.84169 80.06 0.7740
## 9-5 42.3328 -9.61997 94.29 0.2200
## 7-6 6.5723 -45.38041 58.53 1.0000
## 8-6 24.4513 -27.50149 76.40 0.8870
## 9-6 38.6730 -13.27977 90.63 0.3405
## 8-7 17.8789 -34.07381 69.83 0.9838
## 9-7 32.1006 -19.85210 84.05 0.6134
## 9-8 14.2217 -37.73102 66.17 0.9969summary(glht(aov.sleep2,linfct=mcp(Days="Tukey")),adjusted(type="bonferroni"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Reaction ~ Days, data = sleepstudy)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 - 0 == 0 7.844 16.204 0.48 1.00000
## 2 - 0 == 0 8.710 16.204 0.54 1.00000
## 3 - 0 == 0 26.340 16.204 1.63 1.00000
## 4 - 0 == 0 31.998 16.204 1.97 1.00000
## 5 - 0 == 0 51.867 16.204 3.20 0.07358 .
## 6 - 0 == 0 55.526 16.204 3.43 0.03447 *
## 7 - 0 == 0 62.099 16.204 3.83 0.00803 **
## 8 - 0 == 0 79.978 16.204 4.94 8.5e-05 ***
## 9 - 0 == 0 94.199 16.204 5.81 1.3e-06 ***
## 2 - 1 == 0 0.866 16.204 0.05 1.00000
## 3 - 1 == 0 18.496 16.204 1.14 1.00000
## 4 - 1 == 0 24.154 16.204 1.49 1.00000
## 5 - 1 == 0 44.023 16.204 2.72 0.32729
## 6 - 1 == 0 47.682 16.204 2.94 0.16688
## 7 - 1 == 0 54.255 16.204 3.35 0.04506 *
## 8 - 1 == 0 72.134 16.204 4.45 0.00069 ***
## 9 - 1 == 0 86.355 16.204 5.33 1.4e-05 ***
## 3 - 2 == 0 17.630 16.204 1.09 1.00000
## 4 - 2 == 0 23.288 16.204 1.44 1.00000
## 5 - 2 == 0 43.157 16.204 2.66 0.38160
## 6 - 2 == 0 46.816 16.204 2.89 0.19643
## 7 - 2 == 0 53.389 16.204 3.29 0.05393 .
## 8 - 2 == 0 71.268 16.204 4.40 0.00086 ***
## 9 - 2 == 0 85.489 16.204 5.28 1.8e-05 ***
## 4 - 3 == 0 5.657 16.204 0.35 1.00000
## 5 - 3 == 0 25.526 16.204 1.58 1.00000
## 6 - 3 == 0 29.186 16.204 1.80 1.00000
## 7 - 3 == 0 35.759 16.204 2.21 1.00000
## 8 - 3 == 0 53.637 16.204 3.31 0.05123 .
## 9 - 3 == 0 67.859 16.204 4.19 0.00203 **
## 5 - 4 == 0 19.869 16.204 1.23 1.00000
## 6 - 4 == 0 23.529 16.204 1.45 1.00000
## 7 - 4 == 0 30.101 16.204 1.86 1.00000
## 8 - 4 == 0 47.980 16.204 2.96 0.15771
## 9 - 4 == 0 62.202 16.204 3.84 0.00784 **
## 6 - 5 == 0 3.660 16.204 0.23 1.00000
## 7 - 5 == 0 10.232 16.204 0.63 1.00000
## 8 - 5 == 0 28.111 16.204 1.73 1.00000
## 9 - 5 == 0 42.333 16.204 2.61 0.44066
## 7 - 6 == 0 6.572 16.204 0.41 1.00000
## 8 - 6 == 0 24.451 16.204 1.51 1.00000
## 9 - 6 == 0 38.673 16.204 2.39 0.81451
## 8 - 7 == 0 17.879 16.204 1.10 1.00000
## 9 - 7 == 0 32.101 16.204 1.98 1.00000
## 9 - 8 == 0 14.222 16.204 0.88 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)Problème : le plan d’analyse ne prends pas en compte la strate Sujet ni aucune correction.
2ème Solution CAR
On change la structure de nos données vers une version en colonne : cad une modalité par colonne
library(car)
sleepstudy.wide=reshape(sleepstudy, timevar="Days",idvar="Subject",direction="wide")
isleep=data.frame(Days=levels(sleepstudy$Days))
ilm=lm(as.matrix(sleepstudy.wide[,2:11])~1)
ianova<-Anova(ilm,idata=isleep,idesign=~Days)## Note: model has only an intercept; equivalent type-III tests substituted.summary(ianova)##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## Reaction.0 1
## Reaction.1 1
## Reaction.2 1
## Reaction.3 1
## Reaction.4 1
## Reaction.5 1
## Reaction.6 1
## Reaction.7 1
## Reaction.8 1
## Reaction.9 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 160392530
##
## Sum of squares and products for error:
## (Intercept)
## (Intercept) 2506181
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.98 1088 1 17 <2e-16 ***
## Wilks 1 0.02 1088 1 17 <2e-16 ***
## Hotelling-Lawley 1 64.00 1088 1 17 <2e-16 ***
## Roy 1 64.00 1088 1 17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Days
##
## Response transformation matrix:
## Days1 Days2 Days3 Days4 Days5 Days6 Days7 Days8 Days9
## Reaction.0 1 0 0 0 0 0 0 0 0
## Reaction.1 0 1 0 0 0 0 0 0 0
## Reaction.2 0 0 1 0 0 0 0 0 0
## Reaction.3 0 0 0 1 0 0 0 0 0
## Reaction.4 0 0 0 0 1 0 0 0 0
## Reaction.5 0 0 0 0 0 1 0 0 0
## Reaction.6 0 0 0 0 0 0 1 0 0
## Reaction.7 0 0 0 0 0 0 0 1 0
## Reaction.8 0 0 0 0 0 0 0 0 1
## Reaction.9 -1 -1 -1 -1 -1 -1 -1 -1 -1
##
## Sum of squares and products for the hypothesis:
## Days1 Days2 Days3 Days4 Days5 Days6 Days7 Days8 Days9
## Days1 159724 146423 144955 115061 105469 71779 65573 54430 24114
## Days2 146423 134231 132884 105480 96686 65802 60113 49897 22106
## Days3 144955 132884 131552 104422 95717 65142 59510 49397 21884
## Days4 115061 105480 104422 82888 75977 51708 47238 39210 17371
## Days5 105469 96686 95717 75977 69643 47397 43300 35941 15923
## Days6 71779 65802 65142 51708 47397 32257 29468 24460 10837
## Days7 65573 60113 59510 47238 43300 29468 26921 22346 9900
## Days8 54430 49897 49397 39210 35941 24460 22346 18548 8217
## Days9 24114 22106 21884 17371 15923 10837 9900 8217 3641
##
## Sum of squares and products for error:
## Days1 Days2 Days3 Days4 Days5 Days6 Days7 Days8 Days9
## Days1 56091 50059 50770 42221 33178 18544 31964 33318 7810
## Days2 50059 53671 54176 46818 36552 21655 33684 31455 8574
## Days3 50770 54176 62638 54081 42205 25455 43356 39314 13911
## Days4 42221 46818 54081 51885 42260 26557 46309 33511 14521
## Days5 33178 36552 42205 42260 37716 24208 42749 29214 13610
## Days6 18544 21655 25455 26557 24208 23072 35114 19754 14207
## Days7 31964 33684 43356 46309 42749 35114 77925 43451 29871
## Days8 33318 31455 39314 33511 29214 19754 43451 43795 17334
## Days9 7810 8574 13911 14521 13610 14207 29871 17334 17045
##
## Multivariate Tests: Days
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.850 5.663 9 9 0.00825 **
## Wilks 1 0.150 5.663 9 9 0.00825 **
## Hotelling-Lawley 1 5.663 5.663 9 9 0.00825 **
## Roy 1 5.663 5.663 9 9 0.00825 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## SS num Df Error SS den Df F Pr(>F)
## (Intercept) 16039253 1 250618 17 1088.0 <2e-16 ***
## Days 166235 9 151101 153 18.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Days 0.000219 5.15e-08
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Days 0.369 5.5e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Days 0.4694 7.077e-11library(phia)
testInteractions(ilm, idata=isleep, idesign=~Days,adjustment="bonf")## Multivariate Test: Pillai test statistic
## P-value adjustment method: bonf
## Value Df test stat approx F num Df den Df Pr(>F)
## 0-1 -7.8 1 0.103 2.0 1 17 1.00000
## 0-2 -8.7 1 0.074 1.4 1 17 1.00000
## 0-3 -26.3 1 0.347 9.0 1 17 0.35979
## 0-4 -32.0 1 0.402 11.4 1 17 0.16075
## 0-5 -51.9 1 0.535 19.6 1 17 0.01675 *
## 0-6 -55.5 1 0.442 13.5 1 17 0.08559 .
## 0-7 -62.1 1 0.676 35.5 1 17 0.00070 ***
## 0-8 -80.0 1 0.667 34.0 1 17 0.00090 ***
## 0-9 -94.2 1 0.740 48.4 1 17 0.00010 ***
## 1-2 -0.9 1 0.002 0.0 1 17 1.00000
## 1-3 -18.5 1 0.341 8.8 1 17 0.39174
## 1-4 -24.2 1 0.365 9.8 1 17 0.27762
## 1-5 -44.0 1 0.511 17.7 1 17 0.02640 *
## 1-6 -47.7 1 0.389 10.8 1 17 0.19400
## 1-7 -54.3 1 0.605 26.1 1 17 0.00396 **
## 1-8 -72.1 1 0.636 29.7 1 17 0.00194 **
## 1-9 -86.4 1 0.714 42.5 1 17 0.00024 ***
## 2-3 -17.6 1 0.468 15.0 1 17 0.05572 .
## 2-4 -23.3 1 0.380 10.4 1 17 0.22326
## 2-5 -43.2 1 0.491 16.4 1 17 0.03769 *
## 2-6 -46.8 1 0.423 12.5 1 17 0.11594
## 2-7 -53.4 1 0.649 31.4 1 17 0.00143 **
## 2-8 -71.3 1 0.638 30.0 1 17 0.00185 **
## 2-9 -85.5 1 0.677 35.7 1 17 0.00068 ***
## 3-4 -5.7 1 0.102 1.9 1 17 1.00000
## 3-5 -25.5 1 0.349 9.1 1 17 0.34640
## 3-6 -29.2 1 0.292 7.0 1 17 0.76197
## 3-7 -35.8 1 0.445 13.7 1 17 0.08087 .
## 3-8 -53.6 1 0.565 22.1 1 17 0.00932 **
## 3-9 -67.9 1 0.615 27.2 1 17 0.00318 **
## 4-5 -19.9 1 0.365 9.8 1 17 0.27753
## 4-6 -23.5 1 0.248 5.6 1 17 1.00000
## 4-7 -30.1 1 0.414 12.0 1 17 0.13300
## 4-8 -48.0 1 0.601 25.6 1 17 0.00437 **
## 4-9 -62.2 1 0.649 31.4 1 17 0.00143 **
## 5-6 -3.7 1 0.008 0.1 1 17 1.00000
## 5-7 -10.2 1 0.064 1.2 1 17 1.00000
## 5-8 -28.1 1 0.549 20.7 1 17 0.01288 *
## 5-9 -42.3 1 0.583 23.8 1 17 0.00641 **
## 6-7 -6.6 1 0.022 0.4 1 17 1.00000
## 6-8 -24.5 1 0.234 5.2 1 17 1.00000
## 6-9 -38.7 1 0.257 5.9 1 17 1.00000
## 7-8 -17.9 1 0.180 3.7 1 17 1.00000
## 7-9 -32.1 1 0.298 7.2 1 17 0.70714
## 8-9 -14.2 1 0.176 3.6 1 17 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13ème Solution nlme
library(nlme)##
## Attaching package: 'nlme'
##
## L'objet suivant est masqué from 'package:lme4':
##
## lmListlme.sleep=lme(Reaction ~ Days, random = ~ 1 | Subject/Days, data=sleepstudy)
anova(lme.sleep) ## numDF denDF F-value p-value
## (Intercept) 1 153 1088.0 <.0001
## Days 9 153 18.7 <.0001summary(glht(lme.sleep, linfct=mcp(Days="Tukey")),adjusted(type="bonf")) ##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lme.formula(fixed = Reaction ~ Days, data = sleepstudy, random = ~1 |
## Subject/Days)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 7.844 10.475 0.75 1.00000
## 2 - 0 == 0 8.710 10.475 0.83 1.00000
## 3 - 0 == 0 26.340 10.475 2.51 0.53640
## 4 - 0 == 0 31.998 10.475 3.05 0.10142
## 5 - 0 == 0 51.867 10.475 4.95 3.3e-05 ***
## 6 - 0 == 0 55.526 10.475 5.30 5.2e-06 ***
## 7 - 0 == 0 62.099 10.475 5.93 1.4e-07 ***
## 8 - 0 == 0 79.978 10.475 7.63 1.0e-12 ***
## 9 - 0 == 0 94.199 10.475 8.99 < 2e-16 ***
## 2 - 1 == 0 0.866 10.475 0.08 1.00000
## 3 - 1 == 0 18.496 10.475 1.77 1.00000
## 4 - 1 == 0 24.154 10.475 2.31 0.95055
## 5 - 1 == 0 44.023 10.475 4.20 0.00119 **
## 6 - 1 == 0 47.682 10.475 4.55 0.00024 ***
## 7 - 1 == 0 54.255 10.475 5.18 1.0e-05 ***
## 8 - 1 == 0 72.134 10.475 6.89 2.6e-10 ***
## 9 - 1 == 0 86.355 10.475 8.24 1.0e-14 ***
## 3 - 2 == 0 17.630 10.475 1.68 1.00000
## 4 - 2 == 0 23.288 10.475 2.22 1.00000
## 5 - 2 == 0 43.157 10.475 4.12 0.00171 **
## 6 - 2 == 0 46.816 10.475 4.47 0.00035 ***
## 7 - 2 == 0 53.389 10.475 5.10 1.6e-05 ***
## 8 - 2 == 0 71.268 10.475 6.80 4.6e-10 ***
## 9 - 2 == 0 85.489 10.475 8.16 1.0e-14 ***
## 4 - 3 == 0 5.657 10.475 0.54 1.00000
## 5 - 3 == 0 25.526 10.475 2.44 0.66677
## 6 - 3 == 0 29.186 10.475 2.79 0.23999
## 7 - 3 == 0 35.759 10.475 3.41 0.02885 *
## 8 - 3 == 0 53.637 10.475 5.12 1.4e-05 ***
## 9 - 3 == 0 67.859 10.475 6.48 4.2e-09 ***
## 5 - 4 == 0 19.869 10.475 1.90 1.00000
## 6 - 4 == 0 23.529 10.475 2.25 1.00000
## 7 - 4 == 0 30.101 10.475 2.87 0.18266
## 8 - 4 == 0 47.980 10.475 4.58 0.00021 ***
## 9 - 4 == 0 62.202 10.475 5.94 1.3e-07 ***
## 6 - 5 == 0 3.660 10.475 0.35 1.00000
## 7 - 5 == 0 10.232 10.475 0.98 1.00000
## 8 - 5 == 0 28.111 10.475 2.68 0.32780
## 9 - 5 == 0 42.333 10.475 4.04 0.00239 **
## 7 - 6 == 0 6.572 10.475 0.63 1.00000
## 8 - 6 == 0 24.451 10.475 2.33 0.88138
## 9 - 6 == 0 38.673 10.475 3.69 0.01002 *
## 8 - 7 == 0 17.879 10.475 1.71 1.00000
## 9 - 7 == 0 32.101 10.475 3.06 0.09814 .
## 9 - 8 == 0 14.222 10.475 1.36 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)Nous avions bien une anova prenant la source d’erreur “Subject” imbriquée dans notre variable répétée et une analyse post-hoc. La méthode de calcul n’est pas identique : lm qui gère aov fait un traitement par moindre carré et lme calcule par la méthode de maximisation du log de la vraisemblance.
4ème Solution lmer
library(lme4)
str(sleepstudy)## 'data.frame': 180 obs. of 3 variables:
## $ Reaction: num 250 259 251 321 357 ...
## $ Days : Factor w/ 10 levels "0","1","2","3",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...library(pbkrtest) # pour calculer le F## Loading required package: MASS
## Loading required package: parallelsleepstudy$Days=as.factor(sleepstudy$Days)
lme4.sleep=lmer(Reaction ~ Days + (1|Subject) +(1|Days),data=sleepstudy)
lme4.sleep1=lmer(Reaction ~ Days + (1|Subject) ,data=sleepstudy)
Anova(lme4.sleep,test.statistic="F")## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: Reaction
## F Df Df.res Pr(>F)
## Days 1 9 53714 0.44summary(glht(lme4.sleep1, linfct=mcp(Days="Tukey")),adjusted("bonf")) ##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = Reaction ~ Days + (1 | Subject), data = sleepstudy)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 7.844 10.475 0.75 1.00000
## 2 - 0 == 0 8.710 10.475 0.83 1.00000
## 3 - 0 == 0 26.340 10.475 2.51 0.53640
## 4 - 0 == 0 31.998 10.475 3.05 0.10142
## 5 - 0 == 0 51.867 10.475 4.95 3.3e-05 ***
## 6 - 0 == 0 55.526 10.475 5.30 5.2e-06 ***
## 7 - 0 == 0 62.099 10.475 5.93 1.4e-07 ***
## 8 - 0 == 0 79.978 10.475 7.63 1.0e-12 ***
## 9 - 0 == 0 94.199 10.475 8.99 < 2e-16 ***
## 2 - 1 == 0 0.866 10.475 0.08 1.00000
## 3 - 1 == 0 18.496 10.475 1.77 1.00000
## 4 - 1 == 0 24.154 10.475 2.31 0.95055
## 5 - 1 == 0 44.023 10.475 4.20 0.00119 **
## 6 - 1 == 0 47.682 10.475 4.55 0.00024 ***
## 7 - 1 == 0 54.255 10.475 5.18 1.0e-05 ***
## 8 - 1 == 0 72.134 10.475 6.89 2.6e-10 ***
## 9 - 1 == 0 86.355 10.475 8.24 1.0e-14 ***
## 3 - 2 == 0 17.630 10.475 1.68 1.00000
## 4 - 2 == 0 23.288 10.475 2.22 1.00000
## 5 - 2 == 0 43.157 10.475 4.12 0.00171 **
## 6 - 2 == 0 46.816 10.475 4.47 0.00035 ***
## 7 - 2 == 0 53.389 10.475 5.10 1.6e-05 ***
## 8 - 2 == 0 71.268 10.475 6.80 4.6e-10 ***
## 9 - 2 == 0 85.489 10.475 8.16 1.0e-14 ***
## 4 - 3 == 0 5.657 10.475 0.54 1.00000
## 5 - 3 == 0 25.526 10.475 2.44 0.66677
## 6 - 3 == 0 29.186 10.475 2.79 0.23999
## 7 - 3 == 0 35.759 10.475 3.41 0.02885 *
## 8 - 3 == 0 53.637 10.475 5.12 1.4e-05 ***
## 9 - 3 == 0 67.859 10.475 6.48 4.2e-09 ***
## 5 - 4 == 0 19.869 10.475 1.90 1.00000
## 6 - 4 == 0 23.529 10.475 2.25 1.00000
## 7 - 4 == 0 30.101 10.475 2.87 0.18266
## 8 - 4 == 0 47.980 10.475 4.58 0.00021 ***
## 9 - 4 == 0 62.202 10.475 5.94 1.3e-07 ***
## 6 - 5 == 0 3.660 10.475 0.35 1.00000
## 7 - 5 == 0 10.232 10.475 0.98 1.00000
## 8 - 5 == 0 28.111 10.475 2.68 0.32780
## 9 - 5 == 0 42.333 10.475 4.04 0.00239 **
## 7 - 6 == 0 6.572 10.475 0.63 1.00000
## 8 - 6 == 0 24.451 10.475 2.33 0.88138
## 9 - 6 == 0 38.673 10.475 3.69 0.01002 *
## 8 - 7 == 0 17.879 10.475 1.71 1.00000
## 9 - 7 == 0 32.101 10.475 3.06 0.09814 .
## 9 - 8 == 0 14.222 10.475 1.36 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)J’ai conscienceusement évité de rentrer dans la jungle des modèles lineaires mixtes d’éludant la possibilité de comparer les différents modèles proposés comme proposé dans ?lmer sur ces mêmes données. D’ailleurs, je n’ai pas mis le modèle proposé dans l’exemple ?lmer car lmer(Reaction ~ Days + (Days|Subject)) produit une erreur lorque l’on transforme la variable Days en facteur
5ème Solution ezANOVA
sleepstudy$Days=as.factor(sleepstudy$Days)
library(ez)
ez.sleep=ezANOVA(data=sleepstudy,dv=Reaction,wid=Subject,within=Days)
print(ez.sleep)## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Days 9 153 18.7 8.995e-21 * 0.2927
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 Days 0.000219 5.155e-08 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 Days 0.369 5.463e-09 * 0.4694 7.077e-11 *pairwise.t.test( x= sleepstudy$Reaction , g= sleepstudy$Days , p.adjust.method= "bonferroni") #post hoc##
## Pairwise comparisons using t tests with pooled SD
##
## data: sleepstudy$Reaction and sleepstudy$Days
##
## 0 1 2 3 4 5 6 7 8
## 1 1.00000 - - - - - - - -
## 2 1.00000 1.00000 - - - - - - -
## 3 1.00000 1.00000 1.00000 - - - - - -
## 4 1.00000 1.00000 1.00000 1.00000 - - - - -
## 5 0.07358 0.32729 0.38160 1.00000 1.00000 - - - -
## 6 0.03447 0.16688 0.19643 1.00000 1.00000 1.00000 - - -
## 7 0.00803 0.04506 0.05393 1.00000 1.00000 1.00000 1.00000 - -
## 8 8.5e-05 0.00069 0.00086 0.05123 0.15771 1.00000 1.00000 1.00000 -
## 9 1.3e-06 1.4e-05 1.8e-05 0.00203 0.00784 0.44066 0.81451 1.00000 1.00000
##
## P value adjustment method: bonferroni