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 ?Anova
użycia
OBrienKaiser
danych (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:
((treatment %in% c("A", "B")) & (hour %in% 1:2))
przeciw((treatment %in% c("A", "B")) & (hour %in% 3))
((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, gls
któ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 afex
opakowaniem w połączeniu z lsmeans
opakowaniem, jak opisano w afekcie winietowej .
treatment
, 3) dla każdej osoby średnio ponad poziomyprePostFup
, 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.Odpowiedzi:
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
glht
z modelem opartym na prawdopodobieństwie znlme
lublme4
. (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
lm
odpowiedni 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)
W ten sposób uzyskałem te same wartości p:
Przekształć dane w długi format i uruchom,
lm
aby uzyskać wszystkie warunki SS.Stwórz alternatywną matrycę kontrastu dla okresu godzinowego.
Sprawdź, czy moje kontrasty dają takie same SS jak kontrasty domyślne (i takie same jak w przypadku pełnego modelu).
Zdobądź SS i df tylko dla dwóch kontrastów, których chcę.
Uzyskaj wartości p.
Opcjonalnie dostosuj sferyczność.
źródło
heplots
winiety, to naprawdę miłe podsumowanie tego, co się dzieje pod względem ogólnego modelu liniowego.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
hour
w układzie SPF-p.qr (notacja z Kirk (1995): Projekt podziału-wykresu -czynnikowy 1 między współczynnikiemtreatment
z grupami p, pierwszy wewnątrz czynnikahour
z grupami q, drugi wewnątrz czynnikaprePostFup
z grupy r). Poniższe zakładatreatment
grupy o identycznej wielkości i sferyczność.Najpierw zauważ, że główny efekt
hour
jest taki sam po uśrednieniuprePostFup
, a zatem przejście do prostszej konstrukcji SPF-pq, która zawiera tylkotreatment
ihour
jako IV.Zauważmy teraz, że w ANOVA SPF-pq efekt
hour
jest testowany pod kątem interakcjiid:hour
, tj. Ta interakcja zapewnia termin błędu dla testu. Teraz kontrasty dlahour
grup 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 modelulm()
.Ale obliczmy też tutaj wszystko ręcznie.
hour:id
Anova()
car
źródło