Jak wykonać ANCOVA w R.

Chcę przeprowadzić analizę ANCOVA danych dotyczących gęstości epifitów roślinnych. Na początku chciałbym wiedzieć, czy istnieje jakakolwiek różnica w gęstości rośliny między dwoma zboczami, jednym N i jednym S, ale mam inne dane, takie jak wysokość, otwartość czaszy i wysokość rośliny żywicielskiej. Wiem, że moją zmienną towarzyszącą musiałyby być dwa zbocza (N i S). Zbudowałem ten model, który działa w R i chociaż nie mam pojęcia, czy działa dobrze. Chciałbym również wiedzieć, jaka jest różnica, jeśli użyję symbolu +lub *.

model1 <- aov(density~slope+altitude+canopy+height)
summary(model1)
model1

Podstawowym narzędziem do tego jest lm; zauważ, że aovjest to opakowanie dla lm.

W szczególności, jeśli masz jakąś zmienną grupującą (czynnik), $g$ oraz ciągłą zmienną towarzyszącą $x$ , model y ~ x + gpasowałby do modelu ANCOVA z głównymi efektami, natomiasty ~ x * g pasowałby do modelu, który obejmuje interakcję z zmienną towarzyszącą. aovprzyjmie te same formuły.

Zwróć szczególną uwagę Notena pomoc waov .

Jeśli chodzi o +vs *, russellpierce prawie to obejmuje, ale polecam przyjrzeć się, ?lma ?formulazwłaszcza sekcji 11.1 instrukcji Wprowadzenie do R, która pochodzi z R (lub możesz znaleźć go online, jeśli nie zorientowałeś się, jak to zrobić znaleźć go na komputerze; najłatwiej polega na znalezieniu rozwijanego menu „Pomoc” w R lub RStudio).

Polecam pobieranie i czytanie Discovering Statistics using R by Field. Ma fajną sekcję o ANCOVA.

Aby uruchomić ANCOVA w R, załaduj następujące pakiety:

car
compute.es
effects
ggplot2
multcomp
pastecs
WRS

Jeśli używasz lmlub aov(używam aov), upewnij się, że ustawiłeś kontrasty za pomocą funkcji „kontrasty” przed wykonaniem którejkolwiek aovlublm . R domyślnie używa kontrastów nieortogonalnych, które mogą zepsuć wszystko w ANCOVA. Jeśli chcesz ustawić kontrasty ortogonalne, użyj:

contrasts(dataname$factorvariable)=contr.poly(# of levels, i.e. 3) 

następnie uruchom swój model jako

model.1=aov(dv~covariate+factorvariable, data=dataname)

Aby wyświetlić model, użyj:

Anova(model.1, type="III") 

Upewnij się, że używasz tutaj kapitału „A” Anova, a nieanova . To da wyniki przy użyciu SS typu III.

summary.lm(model.1)poda kolejne podsumowanie i obejmie R-sq. wynik.

posth=glht(model.1, linfct=mcp(factorvariable="Tukey"))  ##gives the post-hoc Tukey analysis
summary(posth) ##shows the output in a nice format.

Jeśli chcesz przetestować homogeniczność nachyleń regresji, możesz również dołączyć termin interakcji dla IV i kowariancji. To byłoby:

model=aov(dv~covariate+IV+covariate:IV, data=dataname)

Jeśli termin interakcji jest znaczący, to nie masz jednorodności.

Oto uzupełniająca dokumentacja http://goo.gl/yxUZ1R procedury sugerowanej przez @Butorovich. Ponadto, moje spostrzeżenie jest takie, że gdy zmienna towarzysząca jest binarna, użycie podsumowania (lm.object) dałoby takie samo oszacowanie IV jak wygenerowane przez Anova (lm.object, type = "III").

Używamy analizy regresji do tworzenia modeli opisujących wpływ zmienności zmiennych predykcyjnych na zmienną odpowiedzi. Czasami, jeśli mamy zmienną kategorialną o wartościach takich jak Tak / Nie lub Mężczyzna / Kobieta itp., Prosta analiza regresji daje wiele wyników dla każdej wartości zmiennej jakościowej. W takim scenariuszu możemy zbadać wpływ zmiennej jakościowej, używając jej wraz ze zmienną predykcyjną i porównując linie regresji dla każdego poziomu zmiennej jakościowej. Taka analiza jest określana jako Analiza kowariancji zwana również ANCOVA.

Przykład
RozważR wbudowany zestaw danych mtcars. Obserwujemy w nim, że pole amreprezentuje typ transmisji (automatyczny lub ręczny). Jest to zmienna kategorialna o wartościach 0 i 1. Wartość mil na galon ( mpg) samochodu może również zależeć od niej oprócz wartości mocy konia ( hp). Badamy wpływ wartości amna regresję między mpgi hp. Odbywa się to za pomocą aov()funkcji, po której następujeanova() funkcja porównywania wielu regresji.

Dane wejściowe
Tworzenie ramkę danych zawierającą pola mpg, hpa amze zbioru danych mtcars. Przyjmujemy tutaj mpgjako zmienną odpowiedzi, zmienną hppredykcyjną i amzmienną kategorialną.

input <- mtcars[,c("am","mpg","hp")]
head(input)

Wykonanie powyższego kodu powoduje następujący wynik:

                  am  mpg  hp
Mazda RX4          1 21.0 110
Mazda RX4 Wag      1 21.0 110
Datsun 710         1 22.8  93
Hornet 4 Drive     0 21.4 110
Hornet Sportabout  0 18.7 175
Valiant            0 18.1 105

Analiza ANCOVA
Tworzymy model regresji przyjmujący hpjako zmienną predykcyjną i mpgjako zmienną odpowiedzi uwzględniającą interakcję między ami hp.

Model z interakcją między zmienną kategorialną a zmienną predykcyjną

Utwórz model regresji 1

result1 <- aov(mpg~hp*am,data=mtcars)
summary(result1)

Wykonanie powyższego kodu powoduje następujący wynik:

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4  77.391 1.50e-09 ***
am           1  202.2   202.2  23.072 4.75e-05 ***
hp:am        1    0.0     0.0   0.001    0.981    
Residuals   28  245.4     8.8                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ten wynik pokazuje, że zarówno moc konia, jak i rodzaj transmisji mają znaczący wpływ na mile na galon, ponieważ wartość p w obu przypadkach jest mniejsza niż 0,05. Ale interakcja między tymi dwiema zmiennymi nie jest znacząca, ponieważ wartość p jest większa niż 0,05.

Model bez interakcji między zmienną kategorialną a zmienną predykcyjną

Utwórz model regresji 2

result2 <- aov(mpg~hp+am,data=mtcars)
summary(result2)

Wykonanie powyższego kodu powoduje następujący wynik:

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4   80.15 7.63e-10 ***
am           1  202.2   202.2   23.89 3.46e-05 ***
Residuals   29  245.4     8.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ten wynik pokazuje, że zarówno moc konia, jak i rodzaj transmisji mają znaczący wpływ na mile na galon, ponieważ wartość p w obu przypadkach jest mniejsza niż 0,05.

Porównywanie dwóch modeli
Teraz możemy porównać dwa modele, aby stwierdzić, czy interakcja zmiennych jest naprawdę nieistotna. Do tego używamy anova()funkcji.

 anova(result1,result2)

 Model 1: mpg ~ hp * am
 Model 2: mpg ~ hp + am
   Res.Df    RSS Df  Sum of Sq     F Pr(>F)
 1     28 245.43                           
 2     29 245.44 -1 -0.0052515 6e-04 0.9806

Ponieważ wartość p jest większa niż 0,05, dochodzimy do wniosku, że interakcja między mocą konia a rodzajem przenoszenia nie jest znacząca. Tak więc przebieg na galon będzie zależał w podobny sposób od mocy samochodu w trybie automatycznym i manualnym.