Anova wieloczynnikowa

Anova 2-czynnikowa

W Anovie 2-czynnikowej badamy 2 hipotezy. Pierwsza dotyczy ewentualnej różnicy między wynikami (score) ankiety nt. satysfakcji z pracy osób o różnej płci. Druga dotyczy ewentualnych różnic dla osob z różnym wykształceniem.

Badać możemy również trzecią - dodatkową hipotezę zerową: interakcje między czynnikami (płcią i wykształceniem).

Statystyki opisowe

Zobaczmy jak wygląda rozkład wyników - satysfakcji z pracy pracowników o różnej płci i wykształceniu.

ggplot(jobsatisfaction, aes(x=score)) +
  geom_histogram(bins=10) +
  facet_grid(gender~education_level)

pudelkowy<- ggboxplot(jobsatisfaction, x="gender",
          y="score", color="education_level")
pudelkowy

Zerknijmy, jak wyglądają poszczególne statystyki opisowe wg grup płci i wykształcenia dla ocen z satysfakcji z pracy:

Założenia

Sprawdźmy, czy mamy prawo używać Anovy parametrycznej.

jobsatisfaction %>% 
  group_by(gender, education_level) %>%
  identify_outliers(score)

## [1] gender          education_level id              score          
## [5] is.outlier      is.extreme     
## <0 wierszy> (lub 'row.names' o zerowej długości)

jobsatisfaction %>% 
  group_by(gender, education_level) %>%
  shapiro_test(score)

## # A tibble: 6 × 5
##   gender    education_level variable statistic     p
##   <fct>     <fct>           <chr>        <dbl> <dbl>
## 1 Mężczyzna Średnie         score        0.980 0.966
## 2 Mężczyzna Licencjat       score        0.958 0.779
## 3 Mężczyzna Mgr i wyższe    score        0.916 0.323
## 4 Kobieta   Średnie         score        0.963 0.819
## 5 Kobieta   Licencjat       score        0.963 0.819
## 6 Kobieta   Mgr i wyższe    score        0.950 0.674

ggqqplot(jobsatisfaction, "score") +
  facet_grid(gender~education_level)

jobsatisfaction %>%
  levene_test(score~gender*education_level)

## # A tibble: 1 × 4
##     df1   df2 statistic      p
##   <int> <int>     <dbl>  <dbl>
## 1     5    52      2.20 0.0686

Anova

wyniki <- jobsatisfaction %>% anova_test(score~gender*education_level)
wyniki

## ANOVA Table (type II tests)
## 
##                   Effect DFn DFd       F                           p p<.05
## 1                 gender   1  52   0.745 0.3920000000000000150990331      
## 2        education_level   2  52 187.892 0.0000000000000000000000016     *
## 3 gender:education_level   2  52   7.338 0.0020000000000000000416334     *
##     ges
## 1 0.014
## 2 0.878
## 3 0.220

W powyższej tabeli ANOVY widzimy, że istotne różnice w ocenie satysfakcji z pracy występują dla różnego poziomu wykształcenia. Dodatkowo zauważamy, że istnieją istotne interakcje między płcią a wykształceniem.

ges - oznacza uogólniony (generalized) efekt standaryzowany - czyli ilość wariancji, zróżnicowania objaśnianego przez czynnik.

Testy post-hoc

W przypadku odrzucenia hipotez zerowych w Anovie zwykle przeprowadza się tzw. testy post-hoc aby sprawdzić, co sprawiło, iż pojawia się istotna różnica między średnimi - czyli aby zobaczyć, które dokładnie z nich różnią się istotnie.

Sprawdźmy, jak dokładnie wyniki ocen satysfakcji z pracy różnią się wg poziomu wykształcenia i płci:

wynik2<- jobsatisfaction %>%
  group_by(gender) %>%
  tukey_hsd(score~education_level)
wynik2

## # A tibble: 6 × 10
##   gender    term  group1 group2 null.…¹ estim…² conf.…³ conf.…⁴    p.adj p.adj…⁵
## * <fct>     <chr> <chr>  <chr>    <dbl>   <dbl>   <dbl>   <dbl>    <dbl> <chr>  
## 1 Mężczyzna educ… Średn… Licen…       0   0.797  0.341     1.25 5.59e- 4 ***    
## 2 Mężczyzna educ… Średn… Mgr i…       0   3.87   3.42      4.31 3.28e-14 ****   
## 3 Mężczyzna educ… Licen… Mgr i…       0   3.07   2.62      3.51 3.99e-14 ****   
## 4 Kobieta   educ… Średn… Licen…       0   0.722 -0.0163    1.46 5.62e- 2 ns     
## 5 Kobieta   educ… Średn… Mgr i…       0   2.66   1.93      3.40 4.3 e- 9 ****   
## 6 Kobieta   educ… Licen… Mgr i…       0   1.94   1.20      2.68 1.58e- 6 ****   
## # … with abbreviated variable names ¹​null.value, ²​estimate, ³​conf.low,
## #   ⁴​conf.high, ⁵​p.adj.signif

Jak widzimy powyżej, dla prawie wszystkich par średnich (dla kobiet i mężczyzn vs. wykształcenie) różnice w ocenach satysfakcji z pracy okazały się istotne (p<0.05). Dla jednej pary różnice nie są istotne (ns - non significant) - dla kobiet z wykształceniem college vs. school.

Podsumowanie

Podsumujmy testy Anova za pomocą wykresu ujmującego razem wszystkie poziomy badanych czynników z etykietami i podpisami.

wynik2 <- wynik2 %>% add_xy_position(x="gender")
pudelkowy +
  stat_pvalue_manual(wynik2) +
  labs(
    subtitle=get_test_label(wyniki,detailed=TRUE),
    caption=get_pwc_label(wynik2)
  )

Jak widać, interakcje między czynnikami są istotne (F=7.34, p=0.002).

attach(jobsatisfaction)
interaction.plot(education_level,gender,score,type="b",col=c(2:3),leg.bty="o ",leg.bg="beige",lwd=2,pch=c(18,24) ,xlab="Wykształcenie",ylab="Satysfakcja z pracy",main="Wykres interakcji")

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

ggplot(headache, aes(x=pain_score)) +
  geom_histogram(bins=10) +
  facet_grid(gender~treatment)

pudelkowy<- ggboxplot(headache, x="gender",
          y="pain_score", color="treatment")
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

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

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)
    )