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

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

Poziom bólu głowy osób o różnej płci, ryzyku migreny i rodzaju leczenia

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

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

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.72

Założenia

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

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       FALSE

Występują cztery obserwacje odstające,a jedna z nich to obserwacja odstająca ekstremalnie.

Normalność

Sprawdzenie normalności rozkładu zmiennych.

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.801

Tylko dla jednej grupy występuje odrzucenie H0 o normalności rozkładu(female,high,X). Dla reszty jest spełniona normalność (p>0.05).

ggqqplot(headache, "pain_score") +
  facet_grid(gender~risk~treatment)

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.998

Występuje jednorodnosć warianacji w różnych grupach na podstawie testu Lelevena.(p > 0,05) ### Anova

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

## Coefficient covariances computed by hccm()

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

Istnieje istotna statystycznie trzy czynnikowa interakcja między płcią, ryzykiem i leczeniem(p = 0,001<0,05).

Testy post-hoc

wyniki2<-lm(pain_score~gender*risk*treatment,data=headache)
headache %>%
  group_by(gender) %>%
  anova_test(pain_score~risk*treatment,error=wyniki2)

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

Jak widzimy powyżej, dla czterech par średnich różnice w ocenach poziomu bólu głowy okazały się istotne (p<0.05). Dla dwóch pary różnice nie są istotne (female-treatment, female-risk:treatment).

wyniki3  <- lm(pain_score ~ gender*risk*treatment, data = headache)
 emmip(wyniki3,gender~risk~treatment,cov.reduce=range,xlab="treatment",ylab="pain score")