Comment spécifier des contrastes spécifiques pour des mesures répétées ANOVA en voiture?

12

J'essaie d'exécuter une mesure répétée Anova dans R suivie de quelques contrastes spécifiques sur cet ensemble de données. Je pense que la bonne approche serait d'utiliser à Anova()partir du package de voiture.

Permet d'illustrer ma question avec l'exemple tiré de l' ?Anovautilisation des OBrienKaiserdonnées (Remarque: j'ai omis le facteur de genre de l'exemple):
Nous avons une conception avec un facteur entre les sujets, le traitement (3 niveaux: contrôle, A, B) et 2 répétés -mesures (au sein des sujets) facteurs, phase (3 niveaux: pré-test, post-test, suivi) et heure (5 niveaux: 1 à 5).

La table ANOVA standard est donnée par (à la différence de l'exemple (Anova), je suis passé au type 3 Sums of Squares, c'est ce que mon domaine veut):

require(car)
phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
levels=c("pretest", "posttest", "followup"))
hour <- ordered(rep(1:5, 3))
idata <- data.frame(phase, hour)
mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser)
av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour, type = 3)
summary(av.ok, multivariate=FALSE)

Maintenant, imaginez que l'interaction d'ordre supérieur aurait été significative (ce qui n'est pas le cas) et nous aimerions l'explorer davantage avec les contrastes suivants:
y a-t-il une différence entre les heures 1 et 2 par rapport aux heures 3 (contraste 1) et entre les heures 1 et 2 par rapport aux heures 4 et 5 (contraste 2) dans les conditions de traitement (A&B ensemble)?
En d'autres termes, comment spécifier ces contrastes:

  1. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) contre ((treatment %in% c("A", "B")) & (hour %in% 3))
  2. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) contre ((treatment %in% c("A", "B")) & (hour %in% 4:5))

Mon idée serait d'exécuter une autre ANOVA en omettant la condition de traitement non nécessaire (contrôle):

mod2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser, subset = treatment != "control")
av2 <- Anova(mod2, idata=idata, idesign=~phase*hour, type = 3)
summary(av2, multivariate=FALSE)

Cependant, je ne sais toujours pas comment configurer la matrice de contraste intra-sujet appropriée en comparant les heures 1 et 2 à 3 et 1 et 2 à 4 et 5. Et je ne sais pas si l'omission du groupe de traitement non nécessaire est en effet une bonne idée car elle modifie le terme d'erreur global.

Avant de partir, Anova()je pensais aussi à partir lme. Cependant, il existe de petites différences dans les valeurs F et p entre l'ANOVA du manuel et ce qui est renvoyé en anove(lme) raison de possibles écarts négatifs dans l'ANOVA standard (qui ne sont pas autorisés danslme ). De manière similaire, quelqu'un m'a indiqué glsqui permet d'ajuster des mesures répétées ANOVA, cependant, il n'a aucun argument de contraste.

Pour clarifier: je veux un test F ou t (utilisant des sommes de carrés de type III) qui répond si les contrastes souhaités sont significatifs ou non.


Mise à jour:

J'ai déjà posé une question très similaire sur R-help, il n'y avait pas de réponse .

Une question similaire a été posée sur R-help il y a quelque temps. Cependant, les réponses n'ont pas non plus résolu le problème.


Mise à jour (2015):

Comme cette question génère encore une certaine activité, la spécification des thèses et, fondamentalement, tous les autres contrastes peuvent maintenant être effectués relativement facilement avec le afexpackage en combinaison avec le lsmeanspackage comme décrit dans la vignette afex .

Henrik
la source
1
Avez-vous déjà décidé de ne pas utiliser de tests t? Ce que je veux dire, c'est 1) jeter les données du groupe témoin, 2) ignorer les différents niveaux de treatment, 3) pour chaque personne moyenne sur les niveaux de prePostFup, 4) pour chaque personne moyenne sur les heures 1,2 (= groupe de données 1) ainsi que sur les heures 3,4 (= groupe de données 2), 5) exécuter le test t pour 2 groupes dépendants. Étant donné que Maxwell et Delaney (2004) ainsi que Kirk (1995) découragent de faire des contrastes avec un terme d'erreur groupé dans les plans internes, cela pourrait être une alternative simple.
caracal
Je voudrais faire des analyses de contraste et non des tests t groupés. La raison en est que les contrastes (malgré leurs problèmes) semblent être la procédure standard en psychologie et sont ce que veulent les lecteurs / réviseurs / superviseurs. De plus, ils sont relativement simples à faire dans SPSS. Cependant, malgré mes 2 ans en tant qu'utilisateur R actif jusqu'à présent, je n'ai pas pu y arriver avec R. Maintenant, je dois faire quelques contrastes et je ne veux pas revenir à SPSS juste pour cela. Quand R est le futur (ce que je pense), les contrastes doivent être possibles.
Henrik

Réponses:

6

Cette méthode est généralement considérée comme «à l'ancienne», donc bien que cela soit possible, la syntaxe est difficile et je soupçonne que moins de personnes savent comment manipuler les commandes anova pour obtenir ce que vous voulez. La méthode la plus courante consiste à utiliser glhtun modèle basé sur la vraisemblance à partir de nlmeou lme4. (Je suis certainement le bienvenu si je me trompe par d'autres réponses.)

Cela dit, si je devais le faire, je ne m'embêterais pas avec les commandes anova; Je voudrais simplement ajuster le modèle équivalent en utilisant lm, choisir le bon terme d'erreur pour ce contraste et calculer moi-même le test F (ou de manière équivalente, le test t car il n'y a que 1 df). Cela nécessite que tout soit équilibré et sphérique, mais si vous ne l'avez pas, vous devriez probablement utiliser un modèle basé sur la vraisemblance de toute façon. Vous pourriez peut-être corriger quelque peu la non-sphéricité en utilisant les corrections de Greenhouse-Geiser ou Huynh-Feldt qui (je crois) utilisent la même statistique F mais modifient le df du terme d'erreur.

Si vous voulez vraiment l'utiliser car, vous pourriez trouver les vignettes heplot utiles; ils décrivent comment les matrices du carpackage sont définies.

En utilisant la méthode de caracal (pour les contrastes 1 & 2 - 3 et 1 & 2 - 4 & 5), j'obtiens

      psiHat      tStat          F         pVal
1 -3.0208333 -7.2204644 52.1351067 2.202677e-09
2 -0.2083333 -0.6098777  0.3719508 5.445988e-01

Voici comment j'obtiendrais ces mêmes valeurs de p:

Remodelez les données en format long et exécutez lmpour obtenir tous les termes SS.

library(reshape2)
d <- OBrienKaiser
d$id <- factor(1:nrow(d))
dd <- melt(d, id.vars=c(18,1:2), measure.vars=3:17)
dd$hour <- factor(as.numeric(gsub("[a-z.]*","",dd$variable)))
dd$phase <- factor(gsub("[0-9.]*","", dd$variable), 
                   levels=c("pre","post","fup"))
m <- lm(value ~ treatment*hour*phase + treatment*hour*phase*id, data=dd)
anova(m)

Faites une matrice de contraste alternative pour le terme d'heure.

foo <- matrix(0, nrow=nrow(dd), ncol=4)
foo[dd$hour %in% c(1,2) ,1] <- 0.5
foo[dd$hour %in% c(3) ,1] <- -1
foo[dd$hour %in% c(1,2) ,2] <- 0.5
foo[dd$hour %in% c(4,5) ,2] <- -0.5
foo[dd$hour %in% 1 ,3] <- 1
foo[dd$hour %in% 2 ,3] <- 0
foo[dd$hour %in% 4 ,4] <- 1
foo[dd$hour %in% 5 ,4] <- 0

Vérifiez que mes contrastes donnent le même SS que les contrastes par défaut (et les mêmes que sur le modèle complet).

anova(lm(value ~ hour, data=dd))
anova(lm(value ~ foo, data=dd))

Obtenez les SS et df pour juste les deux contrastes que je veux.

anova(lm(value ~ foo[,1], data=dd))
anova(lm(value ~ foo[,2], data=dd))

Obtenez les valeurs p.

> F <- 73.003/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 2.201148e-09
> F <- .5208/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 0.5445999

Ajustez éventuellement la sphéricité.

pf(F, 1*.48867, 52*.48867, lower=FALSE)
pf(F, 1*.57413, 52*.57413, lower=FALSE)
Aaron a laissé Stack Overflow
la source
Ça marche aussi! Et merci pour le lien vers la heplotsvignette, c'est vraiment un bon résumé de ce qui se passe en termes de modèle linéaire général.
caracal
Merci beaucoup. J'accepterai cette réponse (au lieu de l'autre excellente réponse), car elle inclut une réflexion sur la correction de la sphéricité.
Henrik
Note aux futurs lecteurs: La correction de sphéricité est également applicable à l'autre solution.
Aaron a quitté Stack Overflow
6

Si vous souhaitez / devez utiliser des contrastes avec le terme d'erreur regroupé de l'ANOVA correspondante, vous pouvez procéder comme suit. Malheureusement, ce sera long, et je ne sais pas comment le faire plus facilement. Pourtant, je pense que les résultats sont corrects, car ils sont vérifiés par rapport à Maxwell & Delaney (voir ci-dessous).

Vous voulez comparer des groupes de votre premier facteur intra hourdans un plan SPF-p.qr (notation de Kirk (1995): Plan factoriel divisé en deux 1 entre facteur treatmentavec p groupes, d'abord dans facteur houravec q groupes, second dans facteur prePostFupavec groupes r). Ce qui suit suppose des treatmentgroupes de taille identique et une sphéricité.

Nj    <- 10                                             # number of subjects per group
P     <- 3                                              # number of treatment groups
Q     <- 5                                              # number of hour groups
R     <- 3                                              # number of PrePostFup groups
id    <- factor(rep(1:(P*Nj), times=Q*R))                                  # subject
treat <- factor(rep(LETTERS[1:P], times=Q*R*Nj), labels=c("CG", "A", "B")) # treatment
hour  <- factor(rep(rep(1:Q, each=P*Nj), times=R))                         # hour
ppf   <- factor(rep(1:R, each=P*Q*Nj), labels=c("pre", "post", "fup"))     # prePostFup
DV    <- round(rnorm(Nj*P*Q*R, 15, 2), 2)               # some data with no effects
dfPQR <- data.frame(id, treat, hour, ppf, DV)           # data frame long format

summary(aov(DV ~ treat*hour*ppf + Error(id/(hour*ppf)), data=dfPQR)) # SPF-p.qr ANOVA

Notons d' abord que l'effet principal hourest le même après une moyenne de plus prePostFup, passant ainsi à la conception SPF-pq plus simple qui ne contient que treatmentet hourque IVs.

dfPQ <- aggregate(DV ~ id + treat + hour, FUN=mean, data=dfPQR)  # average over ppf
# SPF-p.q ANOVA, note effect for hour is the same as before
summary(aov(DV ~ treat*hour + Error(id/hour), data=dfPQ))

Notez maintenant que dans l'ANOVA SPF-pq, l'effet de hourest testé par rapport à l'interaction id:hour, c'est-à-dire que cette interaction fournit le terme d'erreur pour le test. Maintenant, les contrastes pour les hourgroupes peuvent être testés comme dans une ANOVA inter-sujets à sens unique en remplaçant simplement le terme d'erreur et les degrés de liberté correspondants. Le moyen le plus simple d'obtenir le SS et le df de cette interaction est d'adapter le modèle avec lm().

(anRes <- anova(lm(DV ~ treat*hour*id, data=dfPQ)))
SSE    <- anRes["hour:id", "Sum Sq"]     # SS interaction hour:id -> will be error SS
dfSSE  <- anRes["hour:id", "Df"]         # corresponding df

Mais calculons également tout manuellement ici.

# substitute DV with its difference to cell / person / treatment group means
Mjk   <- ave(dfPQ$DV,           dfPQ$treat, dfPQ$hour, FUN=mean)  # cell means
Mi    <- ave(dfPQ$DV, dfPQ$id,                         FUN=mean)  # person means
Mj    <- ave(dfPQ$DV,           dfPQ$treat,            FUN=mean)  # treatment means
dfPQ$IDxIV <- dfPQ$DV - Mi - Mjk + Mj                             # interaction hour:id
(SSE  <- sum(dfPQ$IDxIV^2))               # SS interaction hour:id -> will be error SS
dfSSE <- (Nj*P - P) * (Q-1)               # corresponding df
(MSE  <- SSE / dfSSE)                     # mean square

t=ψ^0||c||MSEc||c||ψ^=k=1qckM.kMSEhour:id

Mj     <- tapply(dfPQ$DV, dfPQ$hour, FUN=mean)  # group means for hour
Nj     <- table(dfPQ$hour)                      # cell sizes for hour (here the same)
cntr   <- rbind(c(1, 1, -2,  0, 0),
                c(1, 1, -1, -1, 0))             # matrix of contrast vectors
psiHat <- cntr   %*% Mj                         # estimates psi-hat
lenSq  <- cntr^2 %*% (1/Nj)                     # squared lengths of contrast vectors
tStat  <- psiHat / sqrt(lenSq*MSE)              # t-statistics
pVal   <- 2*(1-pt(abs(tStat), dfSSE))           # p-values
data.frame(psiHat, tStat, pVal)

α

Anova()carϵ^

caracal
la source
Bonne réponse. C'est plus ou moins ce que j'aurais fait si j'avais eu la patience de tout régler.
Aaron a quitté Stack Overflow
Merci pour votre réponse détaillée. Bien que cela semble un peu maladroit dans la pratique.
Henrik