Podział anova i porównań (ortogonalny pojedynczy df) w r

Question

Chcę zrobić kontrast ortogonalny pojedynczego df w anova (model stały lub mieszany). Oto właśnie przykład:

require(nlme) 
data (Alfalfa) 
    Variety: a factor with levels Cossack, Ladak, and Ranger 
    Date : a factor with levels None S1 S20 O7 
    Block: a factor with levels 1 2 3 4 5 6 
    Yield : a numeric vector 

Dane te są opisane w Snedecora i Cochrana (1980) jako przykład z split-plot. Struktura leczenia użyta w eksperymencie była pełną silnią 3 \ times4, z trzema odmianami lucerny i czterema datami trzeciego cięcia w 1943 roku. Jednostki eksperymentalne zaaranżowano na sześć bloków, z których każdy podzielono na cztery działki. Odmiany lucerny (Cossac, Ladak i Ranger) przydzielono losowo do bloków i daty trzeciego cięcia (brak, S1-1 września, S20-wrzesień 20, i O7-7 października) losowo przydzielono do działki. Wszystkie cztery daty zostały użyte w każdym bloku.

model<-with (Alfalfa, aov(Yield~Variety*Date +Error(Block/Date/Variety))) 

    > summary(model) 

Error: Block 
      Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 5 4.15 0.83 

Error: Block:Date 
      Df Sum Sq Mean Sq F value Pr(>F) 
Date  3 1.9625 0.6542 17.84 3.29e-05 *** 
Residuals 15 0.5501 0.0367 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Block:Date:Variety 
      Df Sum Sq Mean Sq F value Pr(>F) 
Variety  2 0.1780 0.08901 1.719 0.192 
Variety:Date 6 0.2106 0.03509 0.678 0.668 
Residuals 40 2.0708 0.05177 

Chcę wykonać pewne porównanie (ortogonalnych kontrastów w grupie), na przykład do tej pory dwa kontrasty:

(a) S1 vs others (S20 O7) 
    (b) S20 vs 07, 

dla czynnika różnorodność dwóch kontrastów:

(c) Cossack vs others (Ladak and Ranger) 
    (d) Ladak vs Ranger 

Zatem Wyjście anova wyglądałoby następująco:

Error: Block 
       Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 5 4.15 0.83 

Error: Block:Date 
      Df Sum Sq Mean Sq F value Pr(>F) 
Date  3 1.9625 0.6542 17.84 3.29e-05 *** 
     (a) S1 vs others ?  ? 
     (b) S20 vs 07  ?  ? 
Residuals 15 0.5501 0.0367 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Block:Date:Variety 
      Df Sum Sq Mean Sq F value Pr(>F) 
Variety  2 0.1780 0.08901 1.719 0.192 
    (c) Cossack vs others ?  ? ? 
    (d) Ladak vs Ranger ?  ?  ? 
Variety:Date 6 0.2106 0.03509 0.678 0.668 
Residuals 40 2.0708 0.05177 

Jak mogę to wykonać? ....................

Answer 1

Po pierwsze, dlaczego warto używać ANOVA? Możesz użyć lme z pakietu nlme, a oprócz testów hipotez aov otrzymasz również możliwe do zinterpretowania oszacowania rozmiarów efektów i kierunków efektów. W każdym razie przychodzą na myśl dwa podejścia:

• Ręczne określ kontrasty zmiennych, zgodnie z wyjaśnieniami here.
• Zainstaluj pakiet multcomp i użyj glht.

glht jest trochę opinii na temat modeli, które są wielowymiarowe w swoich predyktorach. Krótko mówiąc, gdyby utworzyć przekątną macierz cm0 o tych samych wymiarach i wymiarach co model vcov w swoim modelu (załóżmy, że jest to pasujący lme o nazwie model0), wówczas summary(glht(model0,linfct=cm0)) powinien dać takie same oszacowania, SE i test statystyki jako summary(model0)$tTable (ale niepoprawne wartości p-). Teraz, jeśli masz kłopoty z liniowymi kombinacjami wierszy od cm0 i tworz nowe macierze z taką samą liczbą kolumn, co cm0, ale te liniowe kombinacje jako wiersze, w końcu odkryjesz wzór, tworząc matrycę, która da ci przechwycenie oszacuj dla każdej komórki (sprawdź ją pod predict(model0,level=0)). Teraz inna macierz z różnicami między różnymi rzędami tej macierzy da ci odpowiednie różnice między grupami. To samo podejście, ale z wartościami liczbowymi ustawionymi na 1 zamiast 0, można wykorzystać do uzyskania oszacowań nachylenia dla każdej komórki. Następnie różnice między tymi oszacowaniami nachylenia można wykorzystać do uzyskania różnic nachylenia między grupami.

trzy rzeczy, o których należy pamiętać:

• Jak powiedziałem wartości p będą złe dla modeli innych niż lm, (ewentualnie nie próbowałem) aov i niektórych modeli przetrwania. Wynika to z tego, że glht zakłada domyślnie dystrybucję z zamiast t (z wyjątkiem lm). Aby uzyskać poprawne wartości p, wziąć statystyka badania glht oblicza i ręcznie zrobić 2*pt(abs(STAT),df=DF,lower=F) uzyskać dwustronny wartość p gdzie STAT jest statystyka badania zwrócony przez glht i DF jest df z odpowiedniego typu domyślnego kontrastu w summary(model0)$tTable.
• Twoje kontrasty prawdopodobnie nie testują już niezależnych hipotez i konieczna jest wielokrotna korekta testowania, o ile jeszcze nie była. Uruchom wartości p przez p.adjust.
• To jest moja własna destylacja z handwevingu od profesorów i kolegów, oraz dużo czytania Crossvalidated i Stackoverflow na pokrewne tematy. Mogę się mylić na wiele sposobów, a jeśli tak, mam nadzieję, że ktoś bardziej kompetentny poprawi nas obu.