Jak określić konkretne kontrasty dla ANOVA z powtarzanymi pomiarami za pomocą samochodu?

Próbuję uruchomić powtarzane miary Anova w R, a następnie pewne specyficzne kontrasty dla tego zestawu danych. Myślę, że właściwym podejściem byłoby skorzystanie Anova()z pakietu samochodowego.

Zilustrujmy moje pytanie na przykładzie zaczerpniętym z ?Anovaużycia OBrienKaiserdanych (uwaga: odrzuciłem czynnik płci z przykładu):
Mamy wzór z jednym czynnikiem między podmiotami, leczeniem (3 poziomy: kontrola, A, B) i 2 powtórzeniami -pomiary (w obrębie badanych) czynniki, faza (3 poziomy: przedtestowy, posttestowy, kontrolny) i godzina (5 poziomów: 1 do 5).

Standardową tabelę ANOVA podaje (w przeciwieństwie do przykładu (Anova) przełączyłem się na Sumę Kwadratów Typu 3, czyli tego chce moje pole):

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)

Teraz wyobraź sobie, że interakcja najwyższego rzędu byłaby znacząca (co nie jest prawdą) i chcielibyśmy zbadać ją dalej z następującymi kontrastami:
Czy istnieje różnica między godzinami 1 i 2 a godzinami 3 (kontrast 1) oraz między godzinami 1 i 2 a godziny 4 i 5 (kontrast 2) w warunkach leczenia (A&B razem)?
Innymi słowy, jak określić te kontrasty:

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

Moim pomysłem byłoby uruchomienie kolejnej ANOVA, w zależności od niepotrzebnego warunku leczenia (kontrola):

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)

Jednak nadal nie mam pojęcia, jak ustawić odpowiednią matrycę kontrastu w obrębie podmiotu, porównując godziny 1 i 2 z 3 i 1 i 2 z 4 i 5. I nie jestem pewien, czy pominięcie niepotrzebnej grupy terapeutycznej jest rzeczywiście dobrym pomysłem, ponieważ zmienia ogólny termin błędu.

Przed wyjazdem Anova()myślałem również o tym lme. Istnieją jednak niewielkie różnice w wartościach F i p między ANOVA podręcznika a tym, co jest zwracane z anove(lme) powodu możliwych ujemnych odchyleń w standardowej ANOVA (które nie są dozwolone wlme ). W związku z tym ktoś wskazał mi, glsktóry pozwala na dopasowanie ANOVA z powtarzanymi pomiarami, jednak nie ma argumentu kontrastowego.

Wyjaśnienie: Chcę testu F lub t (wykorzystującego sumy kwadratów typu III), który odpowiada, czy pożądane kontrasty są znaczące, czy nie.

Aktualizacja:

Zadałem już bardzo podobne pytanie na temat pomocy R, nie było odpowiedzi .

Podobne pytania zadawano R-help jakiś czas temu. Jednak odpowiedzi również nie rozwiązały problemu.

Aktualizacja (2015):

Ponieważ to pytanie wciąż generuje pewną aktywność, określenie tez i zasadniczo wszystkich innych kontrastów można teraz zrobić stosunkowo łatwo z afexopakowaniem w połączeniu z lsmeansopakowaniem, jak opisano w afekcie winietowej .

Ta metoda jest ogólnie uważana za „staromodną”, więc chociaż może być to możliwe, składnia jest trudna i podejrzewam, że mniej osób wie, jak manipulować poleceniami anova, aby uzyskać to, czego chcesz. Bardziej powszechną metodą jest używanie glhtz modelem opartym na prawdopodobieństwie z nlmelub lme4. (Z pewnością jestem mile widziany, gdy udowodnili mi, że popełniają błędy inne odpowiedzi).

To powiedziawszy, gdybym musiał to zrobić, nie zawracałbym sobie głowy poleceniami anova; Po prostu pasowałbym do modelu równoważnego, wybrałem lmodpowiedni termin błędu dla tego kontrastu i sam obliczę test F (lub równoważnie test t, ponieważ jest tylko 1 df). Wymaga to zrównoważenia wszystkiego i posiadania sferyczności, ale jeśli go nie masz, prawdopodobnie i tak powinieneś użyć modelu opartego na prawdopodobieństwie. Być może będziesz w stanie nieco poprawić niesferyczność za pomocą poprawek Greenhouse-Geisera lub Huynha-Feldta, które (jak sądzę) używają tej samej statystyki F, ale modyfikują df terminu błędu.

Jeśli naprawdę chcesz użyć car, pomocne mogą być winiety heplot ; opisują, jak macierze wcar definiowania pakiecie.

Używam metody karakala (dla kontrastów 1 i 2 - 3 oraz 1 i 2 - 4 i 5)

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

W ten sposób uzyskałem te same wartości p:

Przekształć dane w długi format i uruchom, lmaby uzyskać wszystkie warunki 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)

Stwórz alternatywną matrycę kontrastu dla okresu godzinowego.

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

Sprawdź, czy moje kontrasty dają takie same SS jak kontrasty domyślne (i takie same jak w przypadku pełnego modelu).

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

Zdobądź SS i df tylko dla dwóch kontrastów, których chcę.

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

Uzyskaj wartości 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

Opcjonalnie dostosuj sferyczność.

pf(F, 1*.48867, 52*.48867, lower=FALSE)
pf(F, 1*.57413, 52*.57413, lower=FALSE)

Jeśli chcesz / musisz użyć kontrastów z połączonym terminem błędu z odpowiedniej ANOVA, możesz wykonać następujące czynności. Niestety, to potrwa długo i nie wiem, jak to zrobić wygodniej. Mimo to uważam, że wyniki są prawidłowe, ponieważ są one weryfikowane w stosunku do Maxwella i Delaneya (patrz poniżej).

Chcesz porównać grupy pierwszego czynnika wewnątrz hourw układzie SPF-p.qr (notacja z Kirk (1995): Projekt podziału-wykresu -czynnikowy 1 między współczynnikiem treatmentz grupami p, pierwszy wewnątrz czynnika hourz grupami q, drugi wewnątrz czynnika prePostFupz grupy r). Poniższe zakłada treatmentgrupy o identycznej wielkości i sferyczność.

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

Najpierw zauważ, że główny efekt hourjest taki sam po uśrednieniu prePostFup, a zatem przejście do prostszej konstrukcji SPF-pq, która zawiera tylko treatmenti hourjako IV.

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))

Zauważmy teraz, że w ANOVA SPF-pq efekt hourjest testowany pod kątem interakcji id:hour, tj. Ta interakcja zapewnia termin błędu dla testu. Teraz kontrasty dla hourgrup mogą być testowane tak jak w jednokierunkowej ANOVA między podmiotami, po prostu zastępując termin błędu i odpowiadające mu stopnie swobody. Najłatwiejszym sposobem uzyskania SS i df tej interakcji jest dopasowanie do modelu 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

Ale obliczmy też tutaj wszystko ręcznie.

# 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 = \frac{\hat{\psi} - 0}{||c|| \sqrt{MS_{E}}}$$c$$||c||$$\hat{\psi} = \sum\limits_{k=1}^{q} c_{k} M_{.k}$$MS_{E}$hour: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)

$\alpha$ metodach korekcji , np. Bonferroni.

Anova()car$\hat{\epsilon}$