Wieloczynnikowa analiza wariancji
Podczas Laboratorium 7 zapoznali¶my siê z jednoczynnikow± analiz± wariancji. W praktyce jednak rzadko siê zdarza aby pojedyncza zmienna pozwala³a wyja¶niæ dane zjawisko. Przyk³adowo, je¶li chcieliby¶my dowiedzieæ siê jak zmaksymalizowaæ efekt dzia³ania jakiego¶ lekarstwa to nale¿a³oby wzi±æ pod uwagê takie czynniki jak dawka, sposób podania, a byæ mo¿e równie¿ wiek i p³eæ pacjenta. W wieloczynnikowej analizie wariancji mo¿emy rozdzieliæ ca³kowit± zmienno¶æ na wariancjê wewn±trzgrupow± (nazywan± wariancj± b³êdu lub b³êdem ¶redniokwadratowym), wariancjê miêdzygrupow± (nazywan± wariancj± efektu lub efektem ¶redniokwadratowym) oraz interakcjê miêdzyczynnikow±. Za³o¿enia wieloczynnikowej analizy wariancji s± takie same jak w przypadku jednoczynnikowej, tzn. obserwacje wewn±trz grup s± rozdystrybuowane normalnie a wariancje w ka¿dej grupie s± takie same.
Rozwa¿my przyk³ad dwuczynnikowej analizy wariancji.
# Przyk³ad 9.1
pigs <- ToothGrowth
pigs$dose <- as.factor(pigs$dose)
table(pigs$supp, pigs$dose)##
## 0.5 1 2
## OJ 10 10 10
## VC 10 10 10Nasz± odpowiedzi± jest d³ugo¶æ komórek odpowiedzialnych za wzrost zêbów u ¶winek morskich. Zmienne wyja¶niaj±ce to dawka witaminy C (0.5, 1.0 lub 2.0 mg/dzieñ) oraz sposób dostarczenia witaminy (sok pomarañczowy zakodowany jako OJ oraz kwas askorbinowy jako VC). Mamy zatem 6 grup, które w tym przypadku s± równoliczne (balanced design), co jest istotne dla naszych badañ. W przypadku nierównolicznych grup nale¿y u¿yæ specjalnej wersji ANOVA. Dokonajmy wizualizacji tych danych przy pomocy wykresów pude³kowych.
library(ggplot2)
theme_set(theme_bw())
g <- ggplot(pigs)
g + geom_boxplot(aes(x=dose, y=len, color=supp)) + labs(x="dose [mg/day]", y="length of cells", color="Suppl. type")
Wykonany przez nas test ANOVA bêdzie sprawdza³ nastêpuj±ce hipotezy statystyczne:
- Przyjmowana dawka witaminy C nie wp³ywa na d³ugo¶æ komórek (¶rednie warto¶ci d³ugo¶ci komórek w grupach podzielonych ze wzglêdu na dawkê nie ró¿ni± siê istotnie),
- Sposób suplementacji witaminy C nie wp³ywa na d³ugo¶æ komórek (¶rednie warto¶ci d³ugo¶ci komórek w grupach podzielonych ze wzglêdu na sposób suplementacji nie ró¿ni± siê istotnie),
- Nie ma interakcji pomiêdzy dawk± a sposobem suplementacji witaminy C (czynniki te wp³ywaj± niezale¿nie na d³ugo¶æ komórek, o ile w ogóle wp³ywaj±).
# ANOVA bez interakcji
pigs.aov <- aov(len ~ supp+dose, data=pigs)
summary(pigs.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.4 14.02 0.000429 ***
## dose 2 2426.4 1213.2 82.81 < 2e-16 ***
## Residuals 56 820.4 14.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVA z interakcjami, obie komendy poni¿ej s± równowa¿ne
pigs.aov <- aov(len ~ supp+dose+supp:dose, data=pigs)
pigs.aov <- aov(len ~ supp*dose, data=pigs)
summary(pigs.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.4 15.572 0.000231 ***
## dose 2 2426.4 1213.2 92.000 < 2e-16 ***
## supp:dose 2 108.3 54.2 4.107 0.021860 *
## Residuals 54 712.1 13.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretacja wyników:
- P-warto¶æ dla zmiennej
suppwynosi 0.000429 (F=15.572) co oznacza, ¿e typ suplementacji istotnie wp³ywa na d³ugo¶æ komórek odpowiedzialnych za wzrost zêbów, - P-warto¶æ dla zmiennej
dosejest mniejsza ni¿ 2e-16 (F=92.0) co oznacza, ¿e dawka witaminy wp³ywa nawet bardziej na d³ugo¶æ komórek ni¿ typ suplementacji, - P-warto¶æ dla interakcji pomiêdzy
suppidosewynosi 0.02186 co oznacza, ¿e relacja miêdzy dawk± a d³ugo¶ci± komórek zale¿y od metody suplementacji.
Wykonanie testu Tukey'a pozwoli stwierdziæ, które grupy maj± statystycznie ró¿n± ¶redni±. Wykonamy go jedynie dla zmiennej dose, poniewa¿ zmienna supp ma jedynie dwa poziomy i ju¿ stwierdzili¶my poprzednim testem, ¿e ¶rednie dla nich s± statystycznie ró¿ne.
TukeyHSD(pigs.aov, which = "dose")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = len ~ supp * dose, data = pigs)
##
## $dose
## diff lwr upr p adj
## 1-0.5 9.130 6.362488 11.897512 0.0e+00
## 2-0.5 15.495 12.727488 18.262512 0.0e+00
## 2-1 6.365 3.597488 9.132512 2.7e-06Je¶li w analizowanym zbiorze rozwa¿ane grupy nie s± równoliczne (unbalanced design) korzystamy z funkcji Anova() z pakietu car. Funkcja ta posiada dwie metody obliczania sumy kwadratów: Type-II oraz Type-III. Zaleca siê u¿ywaæ metody Type-III.
library(car)
Anova(pigs.aov, type="III")## Anova Table (Type III tests)
##
## Response: len
## Sum Sq Df F value Pr(>F)
## (Intercept) 1750.33 1 132.730 3.603e-16 ***
## supp 137.81 1 10.450 0.002092 **
## dose 885.26 2 33.565 3.363e-10 ***
## supp:dose 108.32 2 4.107 0.021860 *
## Residuals 712.11 54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wielowymiarowa analiza wariancji
Czasem mamy do czynienia z wiêcej ni¿ jedn± zmienn± odpowiedzi równocze¶nie. Je¶li podejrzewamy (na podstawie innych testów), ¿e zmienne te s± niezale¿ne od siebie to mo¿emy wykonaæ zwyk³± analizê wariancji dla sumy tych zmiennych wynikowych. Je¶li jednak te zmienne s± ze sob± powi±zane lub po prostu zmienna wynikowa ma charakter wielowymiarowy to nale¿y skorzystaæ z wielowymiarowej analizy wariancji.
# Przyk³ad 9.2
irysy <- iris
cor(irysy[,-5])## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
## Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
## Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
## Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000Sprawdzimy czy istnieje istotna ró¿nica w d³ugo¶ciach p³atków i dzia³ek pomiêdzy ró¿nymi gatunkami irysów.
irysy.maov <- manova(cbind(Sepal.Length, Petal.Length) ~ Species, data = irysy)
summary(irysy.maov)## Df Pillai approx F num Df den Df Pr(>F)
## Species 2 0.9885 71.829 4 294 < 2.2e-16 ***
## Residuals 147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Istotnie, tak sformu³owana zmienna dwuwymiarowa posiada ró¿ne warto¶ci ¶rednie w zale¿no¶ci od gatunku kwiatu.
Zadanie punktowane
Wczytaj zbiór danych Salaries z biblioteki carData i zapoznaj siê z nim. Potraktuj zmienn± salary jako odpowied¼ a zmienne rank, discipline, sex jako zmienne wyja¶niaj±ce. Przy pomocy funkcji table sprawd¼ czy grupy s± równoliczne. Wykonaj dwa wykresy skrzynkowe, jeden dla mê¿czyzn a drugi dla kobiet. Porównaj ze sob± wyniki trójczynnikowych analiz wariancji wykonanych bez interakcji oraz z interakcjami. Zinterpretuj wyniki. Zadanie wykonaj dwukrotnie, za pierwszym razem przy za³o¿eniu balanced design, a za drugim razem dla unbalanced design.
## , , = Female
##
##
## A B
## AsstProf 6 5
## AssocProf 4 6
## Prof 8 10
##
## , , = Male
##
##
## A B
## AsstProf 18 38
## AssocProf 22 32
## Prof 123 125
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 6.980e+09 6.980e+09 13.62 0.000256 ***
## rank 2 1.371e+11 6.855e+10 133.72 < 2e-16 ***
## discipline 1 1.828e+10 1.828e+10 35.67 5.25e-09 ***
## Residuals 392 2.009e+11 5.126e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 6.980e+09 6.980e+09 13.460 0.000278 ***
## rank 2 1.371e+11 6.855e+10 132.185 < 2e-16 ***
## discipline 1 1.828e+10 1.828e+10 35.257 6.46e-09 ***
## sex:rank 2 2.352e+08 1.176e+08 0.227 0.797203
## sex:discipline 1 4.558e+08 4.558e+08 0.879 0.349069
## rank:discipline 2 4.748e+08 2.374e+08 0.458 0.632997
## sex:rank:discipline 2 1.324e+08 6.620e+07 0.128 0.880195
## Residuals 385 1.996e+11 5.186e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Anova Table (Type III tests)
##
## Response: salary
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.1903e+11 1 232.196 < 2.2e-16 ***
## sex 6.9407e+08 1 1.354 0.2453
## rank 1.4658e+11 2 142.976 < 2.2e-16 ***
## discipline 1.8283e+10 1 35.666 5.254e-09 ***
## Residuals 2.0094e+11 392
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Anova Table (Type III tests)
##
## Response: salary
## Sum Sq Df F value Pr(>F)
## (Intercept) 3.1916e+10 1 61.5463 4.309e-14 ***
## sex 8.0354e+06 1 0.0155 0.901000
## rank 6.0927e+09 2 5.8746 0.003068 **
## discipline 3.4557e+08 1 0.6664 0.414816
## sex:rank 3.4341e+08 2 0.3311 0.718322
## sex:discipline 1.7226e+06 1 0.0033 0.954069
## rank:discipline 3.5697e+08 2 0.3442 0.709011
## sex:rank:discipline 1.3239e+08 2 0.1277 0.880195
## Residuals 1.9965e+11 385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1