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

12

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 .

Henrik
źródło
1
Czy już zdecydowałeś się nie używać testów t? Mam na myśli to: 1) wyrzucenie danych z grupy kontrolnej, 2) zignorowanie różnych poziomów treatment, 3) dla każdej osoby średnio ponad poziomy prePostFup, 4) dla każdej osoby ze średniej godziny 1,2 (= grupa danych 1) jak również w ciągu godzin 3,4 (= grupa danych 2), 5) przeprowadzić test t dla 2 grup zależnych. Ponieważ Maxwell i Delaney (2004), a także Kirk (1995) zniechęcają do robienia kontrastów ze zbiorczym terminem błędu w projektach, może to być prosta alternatywa.
caracal
Chciałbym wykonać analizy kontrastu, a nie połączone testy t. Powodem jest to, że kontrasty (pomimo ich problemów) wydają się być standardową procedurą w psychologii i są tym, czego chcą czytelnicy / recenzenci / superwizorzy. Co więcej, są stosunkowo łatwe do zrobienia w SPSS. Jednak pomimo moich 2 lat jako aktywnego użytkownika R do tej pory nie byłem w stanie tego osiągnąć za pomocą R. Teraz muszę zrobić kilka kontrastów i nie chcę wracać do SPSS tylko po to. Kiedy R jest przyszłością (którą myślę, że jest), kontrasty muszą być możliwe.
Henrik

Odpowiedzi:

6

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)
Aaron opuścił Stack Overflow
źródło
To też działa! I dzięki za link do heplotswiniety, to naprawdę miłe podsumowanie tego, co się dzieje pod względem ogólnego modelu liniowego.
caracal
Wielkie dzięki. Przyjmuję tę odpowiedź (zamiast innej świetnej odpowiedzi), ponieważ zawiera ona pewne przemyślenia na temat korekcji sferyczności.
Henrik
Uwaga dla przyszłych czytelników: Korekta sferyczności ma również zastosowanie do innych rozwiązań.
Aaron opuścił Stack Overflow
6

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=ψ^-0||do||M.S.mido||do||ψ^=k=1qdokM..kM.S.mihour: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)

α metodach korekcji , np. Bonferroni.

Anova()carϵ^

karakal
źródło
Niezła odpowiedź. To mniej więcej to, co zrobiłbym, gdybym miał cierpliwość, aby to wszystko rozwiązać.
Aaron opuścił stos przepełnienia
Dziękuję za szczegółową odpowiedź. Chociaż w praktyce wydaje się to trochę nieprzydatne.
Henrik