Anova 3-czynnikowa
Wprowadzenie
Trójczynnikowa ANOVA jest rozszerzeniem dwukierunkowej ANOVA do oceny, czy istnieje efekt interakcji pomiędzy trzema niezależnymi zmiennymi kategorycznymi (jakościowymi) na ciągłą zmienną wynikową (ilościową).
Wykorzystamy zbiór danych headache [z pakietu datarium], który zawiera ankietę dot. miary bólu w epizodzie migrenowego bólu głowy u 72 uczestników leczonych trzema różnymi metodami. Wśród uczestników jest 36 mężczyzn i 36 kobiet. Mężczyźni i kobiety zostali dalej podzieleni na grupy o niskim lub wysokim ryzyku migreny.
Chcemy zrozumieć, jak każda niezależna zmienna (rodzaj leczenia, ryzyko migreny i płeć) współdziałają w przewidywaniu wyniku bólu.
Statystyki opisowe
headache %>%
group_by(gender, risk, treatment) %>%
get_summary_stats(pain_score, type = "mean_sd")## # A tibble: 12 x 7
## gender risk treatment variable n mean sd
## <fct> <fct> <fct> <chr> <dbl> <dbl> <dbl>
## 1 male high X pain_score 6 92.7 5.12
## 2 male high Y pain_score 6 82.3 5.00
## 3 male high Z pain_score 6 79.7 4.05
## 4 male low X pain_score 6 76.1 3.86
## 5 male low Y pain_score 6 73.1 4.76
## 6 male low Z pain_score 6 74.5 4.89
## 7 female high X pain_score 6 78.9 5.32
## 8 female high Y pain_score 6 81.2 4.62
## 9 female high Z pain_score 6 81.0 3.98
## 10 female low X pain_score 6 74.2 3.69
## 11 female low Y pain_score 6 68.4 4.08
## 12 female low Z pain_score 6 69.8 2.72wykres1 <- ggboxplot(
headache, x = "treatment", y = "pain_score",
color = "risk" , palette = "jco", facet.by = "gender"
)
wykres1
Założenia
Zaczynamy od oceny normalności i jednorodność wariancji badanych prób.
Obserwacje odstające
Identyfikacja wartości odstających według grup:
headache %>%
group_by(gender, risk, treatment) %>%
identify_outliers(pain_score)## # A tibble: 4 x 7
## gender risk treatment id pain_score is.outlier is.extreme
## <fct> <fct> <fct> <int> <dbl> <lgl> <lgl>
## 1 female high X 57 68.4 TRUE TRUE
## 2 female high Y 62 73.1 TRUE FALSE
## 3 female high Z 67 75.0 TRUE FALSE
## 4 female high Z 71 87.1 TRUE FALSEMożemy zauważyć, że próba zawiera 4 obserwacje odstajace w tym jedną wartość skrajnie odstającą (id = 57, kobieta z wysokim ryzykiem migreny przyjmująca lek X).
Normalność
model <- lm(pain_score ~ gender*risk*treatment, data = headache)
# Create a QQ plot of residuals
ggqqplot(residuals(model))
# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))## # A tibble: 1 x 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 residuals(model) 0.982 0.398Rozkład jest normalny, gdyż punkty układają się wzdłuż linii odniesienia. Potwierdza to test Shapiro-Wilka. Wartość p nie jest istotna (p = 0,398), więc możemy założyć normalność.
headache %>%
group_by(gender, risk, treatment) %>%
shapiro_test(pain_score)## # A tibble: 12 x 6
## gender risk treatment variable statistic p
## <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 male high X pain_score 0.958 0.808
## 2 male high Y pain_score 0.902 0.384
## 3 male high Z pain_score 0.955 0.784
## 4 male low X pain_score 0.982 0.962
## 5 male low Y pain_score 0.920 0.507
## 6 male low Z pain_score 0.924 0.535
## 7 female high X pain_score 0.714 0.00869
## 8 female high Y pain_score 0.939 0.654
## 9 female high Z pain_score 0.971 0.901
## 10 female low X pain_score 0.933 0.600
## 11 female low Y pain_score 0.927 0.555
## 12 female low Z pain_score 0.958 0.801Wyniki miary bólu mają rozkład normalny(p>0.05) oprócz grupa kobiet z wysokim ryzykiem migreny przyjmującej lek X, p = 0,0086.
ggqqplot(headache, "pain_score", ggtheme = theme_bw()) +
facet_grid(gender + risk ~ treatment, labeller = "label_both"
)
Jednorodność wariancji
headache %>% levene_test(pain_score ~ gender*risk*treatment)## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 11 60 0.179 0.998Wynik testu Levene’a nie jest istotny(p>0.05). Zakładamy jednerodność wariancji w różnych grupach.
Anova
res.aov <- headache %>% anova_test(pain_score ~ gender*risk*treatment)## Coefficient covariances computed by hccm()res.aov## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 gender 1 60 16.196 0.000163000000000 * 0.213
## 2 risk 1 60 92.699 0.000000000000088 * 0.607
## 3 treatment 2 60 7.318 0.001000000000000 * 0.196
## 4 gender:risk 1 60 0.141 0.708000000000000 0.002
## 5 gender:treatment 2 60 3.338 0.042000000000000 * 0.100
## 6 risk:treatment 2 60 0.713 0.494000000000000 0.023
## 7 gender:risk:treatment 2 60 7.406 0.001000000000000 * 0.198Wystąpiła statystycznie istotna trójczynnikowa interakcja między płcią, ryzykiem i leczeniem, F(2, 60) = 7,41, p = 0,001.
Testy post-hoc
Jeśli istnieje znaczący trójstronny efekt interakcji, można go rozłożyć na:
- Prostą interakcję dwukierunkową: uruchomić interakcję dwukierunkową na każdym poziomie zmiennej trzeciej,
- Prosty efekt główny: uruchom model jednokierunkowy na każdym poziomie drugiej zmiennej,
- proste porównania parami: przeprowadź porównania parami lub inne porównania post-hoc, jeśli to konieczne.
Jeśli nie masz statystycznie istotnej interakcji trójdrożnej, musisz określić, czy masz statystycznie istotną interakcję dwukierunkową z danych wyjściowych ANOVA. Możesz śledzić znaczącą interakcję dwukierunkową za pomocą prostych analiz efektów głównych i porównań parami między grupami, jeśli to konieczne.
model <- lm(pain_score ~ gender*risk*treatment, data = headache)
headache %>%
group_by(gender) %>%
anova_test(pain_score ~ risk*treatment, error = model)## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()## # A tibble: 6 x 8
## gender Effect DFn DFd F p `p<.05` ges
## * <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 male risk 1 60 50.0 0.00000000187 "*" 0.455
## 2 male treatment 2 60 10.2 0.000157 "*" 0.253
## 3 male risk:treatment 2 60 5.25 0.008 "*" 0.149
## 4 female risk 1 60 42.8 0.000000015 "*" 0.416
## 5 female treatment 2 60 0.482 0.62 "" 0.016
## 6 female risk:treatment 2 60 2.87 0.065 "" 0.087W przypadku mężczyzn wynik ten sugeruje, że wpływ leczenia na wynik bólu zależy od ryzyka migreny. Różnice w poziomie bólu głowy są instotne dla 4 zmiennych: mężczyzna vs risk, mężczyzna vs treatment, mężczyna vs risk;treatment, konieta vs risk.