testy statystyczne

Anova 1-czynnikowa

niejednorodność wariancji

Przykad 1 LAB

Zmienną badaną jest mpg - liczba mil, którą można przejechać na galonie paliwa. Chcielibyśmy sprawdzić, czy rozkład mpg różni się istotnie wg różnych charakterystyk samochodów.

H0: Liczba mpg nie różni się istotnie dla różnej liczby cylindrów.

Ha: Co najmniej jedna para średnich różni się istotnie.

ggboxplot

ggboxplot(mtcars, x = "cyl", y = "mpg", 
          color = "cyl", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("4", "6", "8"),
          ylab = "Mil na galon", xlab = "Liczba cylindrów")

1. Norlamność

mtcars %>%
  group_by(cyl) %>%
  shapiro_test(mpg)
## # A tibble: 3 x 4
##     cyl variable statistic     p
##   <dbl> <chr>        <dbl> <dbl>
## 1     4 mpg          0.912 0.261
## 2     6 mpg          0.899 0.325
## 3     8 mpg          0.932 0.323

Wniosek: Dla każdego cyl rozkład jest normalny.

2. jednorodność wariancji


leveneTest(mpg~cyl)


Wniosek: p<<alfa, co świadczy o niejednorodnej wariancji w badanych próbach.

3. ANOVA 1-czynnikowa dla niejednorodnych wariancji

oneway.test(mpg ~ cyl,
  data = mtcars,
  var.equal = FALSE # wariancje niejednorodne
)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  mpg and cyl
## F = 31.624, num df = 2.000, denom df = 18.032, p-value = 1.271e-06

Wniosek: Odrzucamy H0 i stwierdzamy, że co najmniej jedna para średnich mpg różni się istotnie wg liczby cyl.

4. Post-hoc*


model <- aov(mpg~cyl,data=mtcars)

post_hoc<-TukeyHSD(model,“cyl”,ordered=TRUE)

print(post_hoc)


Wniosek: Dla wszystkich par średnich poszczególne testy t odrzuciły H0 o różnicy mpg wg liczby cylindrów. p<0,05

ANOVA 2-czynnikowa

Próby niezależne

Przykład 1 LAB

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(gender). Druga dotyczy ewentualnych różnic dla osob z różnym wykształceniem(labels).

ggboxplot

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

1. Normalność

jobsatisfaction %>% 
  group_by(gender, education_level) %>%
  shapiro_test(score)
## # A tibble: 6 x 5
##   gender education_level variable statistic     p
##   <fct>  <fct>           <chr>        <dbl> <dbl>
## 1 male   school          score        0.980 0.966
## 2 male   college         score        0.958 0.779
## 3 male   university      score        0.916 0.323
## 4 female school          score        0.963 0.819
## 5 female college         score        0.963 0.819
## 6 female university      score        0.950 0.674

Wykresy normalności wg grup

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

Wniosek: Nie ma podstaw, aby odrzucić hipotezę o normalności rozkładów ocen z satysfakcji z pracy dla każdej badanej grupy. Wyniki mają rozkłady normalne (p>0.05).

2. Obserwacje odstające

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)

Wniosek: Brak obserwacji odstajacych

3. Jednorodne wariancje

jobsatisfaction %>%
  levene_test(score~gender*education_level)
## # A tibble: 1 x 4
##     df1   df2 statistic      p
##   <int> <int>     <dbl>  <dbl>
## 1     5    52      2.20 0.0686

Wniosek: p-value > 0.05, a więc stwierdzamy brak istotnej różnicy wariancji (nie możemy odrzucić hipotezy zerowej). Zatem wariancje są jednorodne.

4. ANOVA 2-czynnikowa

wyniki <- jobsatisfaction %>% anova_test(score~gender*education_level)
## Coefficient covariances computed by hccm()
wyniki
## ANOVA Table (type II tests)
## 
##                   Effect DFn DFd       F        p p<.05   ges
## 1                 gender   1  52   0.745 3.92e-01       0.014
## 2        education_level   2  52 187.892 1.60e-24     * 0.878
## 3 gender:education_level   2  52   7.338 2.00e-03     * 0.220

Wniosek: Istotne różnice w ocenie satysfakcji z pracy występują dla różnego poziomu wykształcenia. + Istotne interakcje między płcią a wykształceniem.

5. POst-hoc

wynik2<- jobsatisfaction %>%
  group_by(gender) %>%
  tukey_hsd(score~education_level)
wynik2
## # A tibble: 6 x 10
##   gender term     group1 group2  null.value estimate conf.low conf.high    p.adj
## * <fct>  <chr>    <chr>  <chr>        <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
## 1 male   educati~ school college          0    0.797   0.341       1.25 5.59e- 4
## 2 male   educati~ school univer~          0    3.87    3.42        4.31 3.28e-14
## 3 male   educati~ colle~ univer~          0    3.07    2.62        3.51 3.99e-14
## 4 female educati~ school college          0    0.722  -0.0163      1.46 5.62e- 2
## 5 female educati~ school univer~          0    2.66    1.93        3.40 4.3 e- 9
## 6 female educati~ colle~ univer~          0    1.94    1.20        2.68 1.58e- 6
## # ... with 1 more variable: p.adj.signif <chr>

p.value parami: ns - nie istotne * 0,05 ** 0,01 *** 0,001

Wniosek: 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.

Próby zależne

Przykład 1 LAB

W zbiorze “depression” zawarto wyniki poziom stanu depresjnego u 24 osób w czterech punktach czasowych. połowa badanych otrzymała lek, a pozostali są poddani kontroli. Badanie ma na celu orzec, czy na wynik depresji wpływa czas i leczenie.

data("depression", package = "datarium")
head(depression, 3)
## # A tibble: 3 x 6
##   id    treatment    t0    t1    t2    t3
##   <fct> <fct>     <dbl> <dbl> <dbl> <dbl>
## 1 1     ctr         296   175   187   242
## 2 2     ctr         376   329   236   126
## 3 3     ctr         309   238   150   173

Przekształcenie ramki w potzrebną do wykonania obliczeń ramkę danych na długi format:

depression <- depression %>%
  gather(key = "time", value = "score", t0, t1, t2, t3) %>%
  convert_as_factor(id, time, treatment)
head(depression, 3)
## # A tibble: 3 x 4
##   id    treatment time  score
##   <fct> <fct>     <fct> <dbl>
## 1 1     ctr       t0      296
## 2 2     ctr       t0      376
## 3 3     ctr       t0      309

1. Normalność

depression%>%
  group_by(time,treatment)%>%
  shapiro_test(score)
## # A tibble: 8 x 5
##   treatment time  variable statistic      p
##   <fct>     <fct> <chr>        <dbl>  <dbl>
## 1 ctr       t0    score        0.945 0.572 
## 2 treated   t0    score        0.836 0.0248
## 3 ctr       t1    score        0.918 0.268 
## 4 treated   t1    score        0.953 0.680 
## 5 ctr       t2    score        0.951 0.652 
## 6 treated   t2    score        0.957 0.734 
## 7 ctr       t3    score        0.873 0.0706
## 8 treated   t3    score        0.932 0.398

Wykres normalności:

ggqqplot(depression, "score") + facet_grid(time~treatment)

Wniosek: Tylko dla jednej zmiennej (osoba leczona w czasie 0) hipoteza zerowa zostaje odrzucona. Pozostałe podlegają normalności.

2. Obserwacje odstające

depression%>%
  group_by(time,treatment)%>%
  identify_outliers(score)
## # A tibble: 4 x 6
##   treatment time  id    score is.outlier is.extreme
##   <fct>     <fct> <fct> <dbl> <lgl>      <lgl>     
## 1 ctr       t0    5       150 TRUE       FALSE     
## 2 ctr       t0    8       447 TRUE       FALSE     
## 3 treated   t0    18      138 TRUE       FALSE     
## 4 treated   t0    22      150 TRUE       FALSE

Wniosek: Dane depression posiadają 4 zmienne odstają, które nie odstają na szczęście ekstremalnie. Ten sam wniosek, można wydedukować dzięki zamieszczonemu wyżej wykresowi ramka-wąs

3. Anova dla zmniennych zależnych

anova <- anova_test(data=depression, formula = score~treatment*time,wid = id, within = time)
## Coefficient covariances computed by hccm()
get_anova_table(anova)
## ANOVA Table (type II tests)
## 
##           Effect DFn DFd      F        p p<.05   ges
## 1      treatment   1  88 43.423 3.13e-09     * 0.330
## 2           time   3  88 19.713 7.31e-10     * 0.402
## 3 treatment:time   3  88  3.755 1.40e-02     * 0.113

Wniosek: Zgodnie z testem ANOVA na wynik depresji istotnie wpływa czas oraz leczenie.

Testy post-hoc

Co sprawiło, że nastapiła widoczna różnica?

depression %>%
  group_by(treatment) %>%
  emmeans_test(score~time)
## # A tibble: 12 x 10
##    treatment term  .y.   group1 group2    df statistic            p        p.adj
##  * <chr>     <chr> <chr> <chr>  <chr>  <dbl>     <dbl>        <dbl>        <dbl>
##  1 ctr       time  score t0     t1        88     1.60  0.113        0.676       
##  2 ctr       time  score t0     t2        88     2.98  0.00375      0.0225      
##  3 ctr       time  score t0     t3        88     3.60  0.000520     0.00312     
##  4 ctr       time  score t1     t2        88     1.38  0.172        1           
##  5 ctr       time  score t1     t3        88     2.00  0.0485       0.291       
##  6 ctr       time  score t2     t3        88     0.625 0.534        1           
##  7 treated   time  score t0     t1        88     6.25  0.0000000142 0.0000000855
##  8 treated   time  score t0     t2        88     6.05  0.0000000338 0.000000203 
##  9 treated   time  score t0     t3        88     5.84  0.0000000854 0.000000513 
## 10 treated   time  score t1     t2        88    -0.196 0.845        1           
## 11 treated   time  score t1     t3        88    -0.409 0.683        1           
## 12 treated   time  score t2     t3        88    -0.213 0.832        1           
## # ... with 1 more variable: p.adj.signif <chr>

Wniosek: Dla pięciu z dwunastu powiązanych zmiennych różniece w wynikach poziomu stanu depresyjnego okazują sie istotne (p_value<0.5). Są to: 1. osoba pod obserwacjąT0T2 2. osoba pod obserwacjąT0T3 3. osoba leczonaT0T1 4. osoba leczonaT0T2 5. osoba leczonaTOT3

ANOVA 3-czynnikowa

ANOVA do oceny, czy istnieje efekt interakcji pomiędzy trzema niezależnymi zmiennymi kategorycznymi (jakościowymi) na ciągłą zmienną wynikową (ilościową).

Próby niezależne

Przykład 1 LAB

Chcemy zrozumieć, jak każda niezależna zmienna (rodzaj leczenia, ryzyko migreny i płeć) współdziałają w przewidywaniu wyniku bólu. Hipoteza zerową jest stwierdzenie, że rodzaj leczenia, ryzyko migreny i płeć nie wpływają na wyniki bólu.

ggboxplot

ggboxplot(depression, x = "time", y = "score", color = "treatment", palette = "jco")

Opisówka:

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

1. Normalność

headache %>% 
  group_by(treatment, risk, gender) %>%
  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 female high  X         pain_score     0.714 0.00869
##  3 male   low   X         pain_score     0.982 0.962  
##  4 female low   X         pain_score     0.933 0.600  
##  5 male   high  Y         pain_score     0.902 0.384  
##  6 female high  Y         pain_score     0.939 0.654  
##  7 male   low   Y         pain_score     0.920 0.507  
##  8 female low   Y         pain_score     0.927 0.555  
##  9 male   high  Z         pain_score     0.955 0.784  
## 10 female high  Z         pain_score     0.971 0.901  
## 11 male   low   Z         pain_score     0.924 0.535  
## 12 female low   Z         pain_score     0.958 0.801

Wniosek: Hipoteza zerowa odrzucona jest jedynie dla jednej zmiennej - kobiety o wysokim ryzyku migreny oraz leczeniem lekiem X. Wszytskie pozostałe przypadki posiadają p>0.05 zatem nie ma dla nich podstaw do odrzucenia hipotez zerowych.

Wykres normalności wg podgrup

ggqqplot(headache, x = ("pain_score"), color = "blue", title = "Wykres norlamności według podgrup" )+
  facet_grid(gender~risk~treatment)

2. Obserwacje odstające

headache %>% 
  group_by(treatment, risk, gender) %>%
  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

Wniosek: Badana próba zawiera 4 obserwacje odstajace, w tym jedna to obserwacja ekstremalnie odstająca.

3. Jednrorodnośc 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

Wniosek: p-value > 0.05, Zatem wariancje są jednorodne.

4. ANOVA 3-czynnikowa

headache %>% 
  anova_test(pain_score~gender*risk*treatment)
## Coefficient covariances computed by hccm()
## ANOVA Table (type II tests)
## 
##                  Effect DFn DFd      F        p p<.05   ges
## 1                gender   1  60 16.196 1.63e-04     * 0.213
## 2                  risk   1  60 92.699 8.80e-14     * 0.607
## 3             treatment   2  60  7.318 1.00e-03     * 0.196
## 4           gender:risk   1  60  0.141 7.08e-01       0.002
## 5      gender:treatment   2  60  3.338 4.20e-02     * 0.100
## 6        risk:treatment   2  60  0.713 4.94e-01       0.023
## 7 gender:risk:treatment   2  60  7.406 1.00e-03     * 0.198

Wniosek:

Analizując tabele ANOVY widzimy, że istotne różnice w poziomie bólu głowy występują dla wszytskich trzech czynników: płci, ryzyka oraz podawanego leku. Dodatkowo zauważamy, że istnieją istotne interakcje między płcią a otrzymywanym lekiem oraz trójczynnikowa interakcja płci, leku i ryzyka.

5. Post-hoc lm ()

lm(pain_score~gender*risk*treatment,data=headache)
## 
## Call:
## lm(formula = pain_score ~ gender * risk * treatment, data = headache)
## 
## Coefficients:
##                     (Intercept)                     genderfemale  
##                          92.739                          -13.874  
##                         risklow                       treatmentY  
##                         -16.687                          -10.397  
##                      treatmentZ             genderfemale:risklow  
##                         -13.058                           11.978  
##         genderfemale:treatmentY          genderfemale:treatmentZ  
##                          12.708                           15.228  
##              risklow:treatmentY               risklow:treatmentZ  
##                           7.484                           11.462  
## genderfemale:risklow:treatmentY  genderfemale:risklow:treatmentZ  
##                         -15.589                          -18.009
headache %>%
  group_by(gender) %>%
  anova_test(pain_score~risk*treatment)
## 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    30 45.0   0.000000195 "*"     0.6  
## 2 male   treatment          2    30  9.15  0.000788    "*"     0.379
## 3 male   risk:treatment     2    30  4.72  0.016       "*"     0.24 
## 4 female risk               1    30 48.2   0.000000104 "*"     0.616
## 5 female treatment          2    30  0.542 0.587       ""      0.035
## 6 female risk:treatment     2    30  3.23  0.054       ""      0.177

Wniosek: 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.

Próby zależne

Przykład 2

Czy osiąnięcia uczestnika egzaminu istotnie różnią się w zależności od czasu, płci oraz poziomu stresu?

library(datarium)
data("performance")

Konwertacja:

performance <- performance %>%
  gather(key = "time", value = "score", t1, t2) %>%
  convert_as_factor(id, gender, time)

1. Normalność

performance %>%
  group_by(gender,stress,time) %>%
  shapiro_test(score)
## # A tibble: 12 x 6
##    gender stress   time  variable statistic      p
##    <fct>  <fct>    <fct> <chr>        <dbl>  <dbl>
##  1 male   low      t1    score        0.942 0.574 
##  2 male   low      t2    score        0.966 0.849 
##  3 male   moderate t1    score        0.848 0.0547
##  4 male   moderate t2    score        0.958 0.761 
##  5 male   high     t1    score        0.915 0.319 
##  6 male   high     t2    score        0.925 0.403 
##  7 female low      t1    score        0.898 0.207 
##  8 female low      t2    score        0.886 0.154 
##  9 female moderate t1    score        0.946 0.626 
## 10 female moderate t2    score        0.865 0.0880
## 11 female high     t1    score        0.989 0.996 
## 12 female high     t2    score        0.930 0.452

Wniosek: Dla każdej zmiennej p.value > 0.05 więc nie ma podstaw do odrzucenia hipotezy zerowej.

2. Wartości odstające

performance %>%
  group_by(gender, stress, time)%>%
  identify_outliers(score)
## # A tibble: 1 x 7
##   gender stress time  id    score is.outlier is.extreme
##   <fct>  <fct>  <fct> <fct> <dbl> <lgl>      <lgl>     
## 1 female low    t2    36     6.15 TRUE       FALSE

Tylko jedna zmienna jest wartością odstającą - kobieta o niskim poziomie stresu i czasie t2

3. ANOVA dla prob zależnych

wyniki<-anova_test(data=performance, formula = score~gender*time*stress,wid = id, within = time)
## Coefficient covariances computed by hccm()
get_anova_table(wyniki)
## ANOVA Table (type II tests)
## 
##               Effect DFn DFd      F        p p<.05      ges
## 1             gender   1 108  2.487 1.18e-01       0.023000
## 2               time   1 108  0.061 8.06e-01       0.000564
## 3             stress   2 108 21.878 1.05e-08     * 0.288000
## 4        gender:time   1 108  4.571 3.50e-02     * 0.041000
## 5      gender:stress   2 108  1.607 2.05e-01       0.029000
## 6        time:stress   2 108  1.759 1.77e-01       0.032000
## 7 gender:time:stress   2 108  5.896 4.00e-03     * 0.098000

Tabela anovy przedstawia widoczne różnice w wynikach na sprawdzianie, zależne od poziomu stresu. Istotna jest też interakcja płci i czasu. Ponadto osiąnięcia uczestnika egzaminu istotnie różnią się w zależności od czasu, płci oraz poziomu stresu.

post-hoc

performance %>%
  group_by(gender, stress) %>%
  emmeans_test(score~time)
## # A tibble: 6 x 11
##   gender stress   term  .y.   group1 group2    df statistic       p   p.adj
## * <chr>  <chr>    <chr> <chr> <chr>  <chr>  <dbl>     <dbl>   <dbl>   <dbl>
## 1 female high     time  score t1     t2       108     2.13  0.0351  0.0351 
## 2 female low      time  score t1     t2       108    -1.93  0.0561  0.0561 
## 3 female moderate time  score t1     t2       108     2.72  0.00767 0.00767
## 4 male   high     time  score t1     t2       108    -1.95  0.0537  0.0537 
## 5 male   low      time  score t1     t2       108     0.266 0.791   0.791  
## 6 male   moderate time  score t1     t2       108    -0.631 0.529   0.529  
## # ... with 1 more variable: p.adj.signif <chr>

Jak widzimy powyżej, dla TYLKO DWÓCH par średnich (dla kobiet vs. poziom stresu) 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).