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.
set.seed(123)
data("headache", package = "datarium")
headache %>% sample_n_by(gender, risk, treatment, size = 1)## # A tibble: 12 x 5
## id gender risk treatment pain_score
## <int> <fct> <fct> <fct> <dbl>
## 1 21 male high X 92.7
## 2 30 male high Y 80.1
## 3 33 male high Z 81.3
## 4 2 male low X 76.8
## 5 8 male low Y 80.7
## 6 18 male low Z 74.7
## 7 57 female high X 68.4
## 8 65 female high Y 82.1
## 9 70 female high Z 80.4
## 10 42 female low X 78.1
## 11 48 female low Y 64.1
## 12 49 female low Z 69.0Statystyki 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.72Założenia
rysuneczek <- ggboxplot(
headache, x = "treatment", y = "pain_score",
color = "risk", palette = "jco", facet.by = "gender"
)
rysuneczek
Obserwacje odstające
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 FALSEWystępuje tylko jedna wartość odstająca.
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.801Dla większości wariantów p-value jest większe od alfy, czyli jest spełniony warunek o normalności rozkładów.
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.998p-value jest większe od alfy, czyli jest zachowana jednorodność wariancji.
Anova
wynik <- headache %>% anova_test(pain_score ~ gender*risk*treatment)## Coefficient covariances computed by hccm()wynik## 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.198Istotne różnice nie występują dla par czynników płeć i ryzyko oraz ryzyko i lek.
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.087