Anova wieloczynnikowa

Challenge

Wykorzystamy zestaw danych dotyczących bólu głowy [pakiet datarium], który zawiera miary wyniku bólu epizodu migrenowego bólu głowy u 72 uczestników leczonych trzema różnymi metodami.

Wśród uczestników znalazło się 36 mężczyzn i 36 kobiet. Mężczyźni i kobiety byli dalej podzieleni w podgrupy, zależnie od tego czy byli na niskim lub wysokim ryzyku migreny.

Chcemy zrozumieć, jak każda zmienna niezależna (rodzaj zabiegów, ryzyko migreny i płeć) współdziałają w przewidywaniu wyniku bólu.

Załaduj dane i skontroluj jeden losowy rząd według kombinacji grup:

set.seed(123)
data("headache", package = "datarium")
headache %>% sample_n_by(gender, risk, treatment, size = 1)

## # A tibble: 12 × 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.0

pudelkowy <- ggboxplot(
  headache, x = "treatment", y = "pain_score", 
  color = "risk", palette = "jco", facet.by = "gender"
  )
pudelkowy

ggqqplot(headache, "pain_score", ggtheme = theme_bw()) +
  facet_grid(gender + risk ~ treatment, labeller = "label_both")

headache %>% levene_test(pain_score ~ gender*risk*treatment)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1    11    60     0.179 0.998

wyniki <- headache %>% anova_test(pain_score~gender*risk*treatment)
wyniki

## 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.198

Testy post-hoc

library(emmeans)

## Warning: pakiet 'emmeans' został zbudowany w wersji R 4.2.2

parami <- headache %>%
  group_by(gender, risk) %>%
  emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni") %>%
  select(-df, -statistic, -p) 

parami %>% filter(gender == "male", risk == "high")

## # A tibble: 3 × 8
##   gender risk  term      .y.        group1 group2      p.adj p.adj.signif
##   <chr>  <chr> <chr>     <chr>      <chr>  <chr>       <dbl> <chr>       
## 1 male   high  treatment pain_score X      Y      0.000386   ***         
## 2 male   high  treatment pain_score X      Z      0.00000942 ****        
## 3 male   high  treatment pain_score Y      Z      0.897      ns

Przeprowadzono trójczynnikową ANOVA w celu określenia wpływu płci, ryzyka i leczenia na pain_score epizodu migrenowego bólu głowy.

Istniała statystycznie istotna trójstronna interakcja między płcią, ryzykiem i leczeniem, F(2, 60) = 7,41, p = 0,001.

Istniała statystycznie istotna prosta dwukierunkowa interakcja między ryzykiem a leczeniem dla mężczyzn, F(2, 60) = 5,2, p = 0,008, ale nie dla kobiet, F(2, 60) = 2,8, p = 0,065.

Istniał statystycznie istotny prosty efekt główny leczenia dla mężczyzn z wysokim ryzykiem migreny, F(2, 60) = 14,8, p < 0,0001), ale nie dla mężczyzn z niskim ryzykiem migreny, F(2, 60) = 0,66, p = 0,521.

Wystąpiła statystycznie istotna średnia różnica między leczeniem X a leczeniem Y. Jednak różnica między leczeniem Y a leczeniem Z, nie była statystycznie istotna.

parami <- parami %>% add_xy_position(x = "treatment")
filtrowane_dane <- parami %>% filter(gender == "male", risk == "high")
pudelkowy +
  stat_pvalue_manual(
    filtrowane_dane, color = "risk", linetype = "risk", hide.ns = TRUE,
    tip.length = 0, step.increase = 0.1, step.group.by = "gender"
  ) +
  labs(
    subtitle = get_test_label(wyniki, detailed = TRUE),
    caption = get_pwc_label(parami)
    )