Jednofaktorová ANOVA

Zadanie: Chcete porovnať priemerný úbytok hmotnosti pre tri rôzne diétne plány (diéta A, diéta B a diéta C). Zhromažďujete údaje od skupín účastníkov dodržiavajúcich každý diétny plán počas jedného mesiaca. Otestujte na hladine významnosti 0.01 pomocou testu implementovaného priamo v Rku. Rozhodnite sa na základe P hodnoty testu. Ak je to potrebné, tak urobte aj post testy. Testujeme nulovú hypotézu oproti alternatívnej:

𝐻0: 𝜇1=𝜇2=𝜇3 𝐻1: ∃ 𝑖, 𝑗: 𝜇𝑖 ≠ 𝜇𝑗

# Dáta o úbytku hmotnosti pre tri diéty
diet_A <- c(2.5, 3.0, 2.8, 3.2, 2.7)
diet_B <- c(3.5, 3.8, 3.2, 4.0, 3.6)
diet_C <- c(2.0, 2.5, 2.2, 2.3, 2.1)

# Vytvorenie dátového rámca
df_diet <- data.frame(
  weight_loss = c(diet_A, diet_B, diet_C),
  diet = factor(rep(c("A", "B", "C"), each = 5))
)

# Jednofaktorová ANOVA
df_diet$diet<-as.factor(df_diet$diet)
head(df_diet, 15)
##    weight_loss diet
## 1          2.5    A
## 2          3.0    A
## 3          2.8    A
## 4          3.2    A
## 5          2.7    A
## 6          3.5    B
## 7          3.8    B
## 8          3.2    B
## 9          4.0    B
## 10         3.6    B
## 11         2.0    C
## 12         2.5    C
## 13         2.2    C
## 14         2.3    C
## 15         2.1    C

Grafické zobrazenie:

library(ggplot2)
ggplot(df_diet, aes(x= diet, y= weight_loss))+
  geom_boxplot(fill=c("#4472c4", "#FF6666", "#00cc99"))+
  labs(x= "Diétny plán", y= "Úbytok hmotnosti", title = "Porovnanie úbytky hmotnosti v závislosti od diétneho plánu")+
  theme(plot.title = element_text(hjust=0.5))

Z grafu vidíme, že dieta C má nižši úbytok hmotnosti ako ostatné diety a dieta A ma najväčší úbytok hmotnosti. Použitím analýzy rozptylu zistíme, či je tento rozdiel štatisticky významný.

Normalita- použijeme Shapiro Wilkov test. 𝐻0: 𝑋𝑖 -ty výber je z normálneho rozdelenia, 𝑖=1,2,…,𝑘 , kde 𝑘 - počet úrovní faktora.

𝐻1: ∃ 𝑋𝑖 , ktorý nepochádza z normálneho rozdelenia.

Keďže pre všetky 3 úrovne faktora je p-hodnota väčšia ako 0.01, 𝐻0 nemôžeme zamietnuť.

tapply(df_diet$weight_loss, df_diet$diet, shapiro.test)
## $A
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.98998, p-value = 0.9796
## 
## 
## $B
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.99174, p-value = 0.9854
## 
## 
## $C
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.97872, p-value = 0.9276

Pre všetky 3 úrovne faktora má úbytok hmotnosti normálne rozdelenie. Prvý predpoklad použitia parametrickej ANOVY je splnený.

Rovnosť rozptylov. Použijeme Bartlettov test. 𝐻𝑂: 𝜎1=𝜎2=…=𝜎𝑘

𝐻1: ∃ 𝑖,𝑗: 𝜎𝑖≠𝜎𝑗

bartlett.test(df_diet$weight_loss, df_diet$diet)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  df_diet$weight_loss and df_diet$diet
## Bartlett's K-squared = 0.74088, df = 2, p-value = 0.6904

Na základe p-hodnoty 𝐻0 nemôžeme zmietnuť na hladine významnosti 𝛼=0.01. Teda aj druhý predpoklad je dodržaný, rozptyly nie sú štatisticky významne odlišné pre rôzne úrovne faktora.

Samotná ANOVA

an1<-aov(df_diet$weight_loss ~ df_diet$diet)
summary(an1)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## df_diet$diet  2  4.921  2.4607   36.55 7.87e-06 ***
## Residuals    12  0.808  0.0673                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(an1, alpha=0.01)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = df_diet$weight_loss ~ df_diet$diet)
## 
## $`df_diet$diet`
##      diff        lwr        upr     p adj
## B-A  0.78  0.3421668  1.2178332 0.0012611
## C-A -0.62 -1.0578332 -0.1821668 0.0068724
## C-B -1.40 -1.8378332 -0.9621668 0.0000054
library(agricolae)
## Warning: пакет 'agricolae' был собран под R версии 4.4.2
an11<-aov(weight_loss ~ diet, data=df_diet)
result<- scheffe.test(an11,"diet", alpha = 0.01)
result
## $statistics
##      MSerror Df        F     Mean       CV  Scheffe CriticalDifference
##   0.06733333 12 6.926608 2.893333 8.968433 3.721991            0.61083
## 
## $parameters
##      test name.t ntr alpha
##   Scheffe   diet   3  0.01
## 
## $means
##   weight_loss       std r       se Min Max Q25 Q50 Q75
## A        2.84 0.2701851 5 0.116046 2.5 3.2 2.7 2.8 3.0
## B        3.62 0.3033150 5 0.116046 3.2 4.0 3.5 3.6 3.8
## C        2.22 0.1923538 5 0.116046 2.0 2.5 2.1 2.2 2.3
## 
## $comparison
## NULL
## 
## $groups
##   weight_loss groups
## B        3.62      a
## A        2.84      b
## C        2.22      c
## 
## attr(,"class")
## [1] "group"

Na základe oboch post testov je štatisticky významný rozdiel medzi diétnymi plánmi A, B a C. Diétny plán B je účinnejší ako diétny plán A a C a diétny plán A je účinnejší ako diétny plán C.

Jednofaktorová ANOVA s F testom

Zadanie: Chcete otestovať, či tri rôzne hnojivá (hnojivo A, B a C) vedú k rôznym výnosom plodín. Zhromažďujete údaje o výnosoch z polí, na ktorých sa aplikovalo každé hnojivo, a zaujíma vás, či existuje významný rozdiel medzi priemernými výnosmi. Otestujte na hladine významnosti 0.1 pomocou testu implementovaného priamo v Rku. Rozhodnite sa na základe F štatistiky. Ak je to potrebné, tak urobte aj post testy.

fertilizer_A <- c(10.2, 10.5, 9.8, 10.3, 10.4)  
fertilizer_B <- c(11.1, 11.5, 11.3, 10.8, 11.0)  
fertilizer_C <- c(9.5, 9.8, 9.7, 9.6, 9.9) 

df_diet2<-data.frame(yield= c(fertilizer_A, fertilizer_B, fertilizer_C),
                fertilizer = rep(c("A", "B", "C"), c(length(fertilizer_A), length(fertilizer_B), length(fertilizer_C)))
                )
df_diet2$fertilizer<-as.factor(df_diet2$fertilizer)
head(df_diet2, 15)
##    yield fertilizer
## 1   10.2          A
## 2   10.5          A
## 3    9.8          A
## 4   10.3          A
## 5   10.4          A
## 6   11.1          B
## 7   11.5          B
## 8   11.3          B
## 9   10.8          B
## 10  11.0          B
## 11   9.5          C
## 12   9.8          C
## 13   9.7          C
## 14   9.6          C
## 15   9.9          C

Grafické zobrazenie:

ggplot(df_diet2, aes(x= fertilizer, y= yield))+
  geom_boxplot(fill=c("#4472c4", "#FF6666", "#00cc99"))+
  labs(x= "Hnojivo", y= "Výnos plodín", title = "Porovnanie výnosov plodín v závislosti od hnojív")+
  theme(plot.title = element_text(hjust=0.5))

Z grafu vidíme, že hnojivo C má nižši výnos plodin ako ostatné hnojíva a hnojivo A ma najväčší výnos plodin. Použitím analýzy rozptylu zistíme, či je tento rozdiel štatisticky významný.

Predpoklady

normalita- použijeme Shapiro Wilkov test. 𝐻0: 𝑋𝑖 -ty výber je z normálneho rozdelenia, 𝑖=1,2,…,𝑘 , kde 𝑘 - počet úrovní faktora.

𝐻1: ∃ 𝑋𝑖 , ktorý nepochádza z normálneho rozdelenia.

Keďže pre všetky 3 úrovne faktora je p-hodnota väčšia ako 0.1, 𝐻0 nemôžeme zamietnuť.

tapply(df_diet2$yield, df_diet2$fertilizer, shapiro.test)
## $A
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.90309, p-value = 0.4272
## 
## 
## $B
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.98998, p-value = 0.9796
## 
## 
## $C
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.98676, p-value = 0.9672

Pre všetky 3 úrovne faktora má úbytok hmotnosti normálne rozdelenie. Prvý predpoklad použitia parametrickej ANOVY je splnený.

Rovnosť rozptylov. Použijeme Bartlettov test. 𝐻𝑂: 𝜎1=𝜎2=…=𝜎𝑘

𝐻1: ∃ 𝑖,𝑗: 𝜎𝑖≠𝜎𝑗

bartlett.test(df_diet2$yield, df_diet2$fertilizer)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  df_diet2$yield and df_diet2$fertilizer
## Bartlett's K-squared = 1.1857, df = 2, p-value = 0.5528

Na základe p-hodnoty 𝐻0 nemôžeme zmietnuť na hladine významnosti 𝛼=0.01. Teda aj druhý predpoklad je dodržaný, rozptyly nie sú štatisticky významne odlišné pre rôzne úrovne faktora.

Samotná ANOVA

an2<-aov(df_diet2$yield ~ df_diet2$fertilizer)
summary(an2)
##                     Df Sum Sq Mean Sq F value   Pr(>F)    
## df_diet2$fertilizer  2  5.292   2.646   46.42 2.25e-06 ***
## Residuals           12  0.684   0.057                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Podľa výsledku testu vidíme, že F-hodnota (46.42) je rádovo väčšia ako jednota. Teda hypotézu o nulovosti faktora môžeme zamietnuť. Môžeme teda povedať, že rozdiely vo výnosoch plodín v závislosti od hnojív sú štatisticky významné na hladine významnosti 0.1.

Na určenie, ktorý plan diety sa štatisticky významne odlišujé, budeme testovať kontrasty (alebo inak post hoc test, viacnásobné porovnávanie priemerov). Použijeme Tukeyov a Scheffeho test.

TukeyHSD(an2, alpha=0.1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = df_diet2$yield ~ df_diet2$fertilizer)
## 
## $`df_diet2$fertilizer`
##      diff        lwr        upr     p adj
## B-A  0.90  0.4971614  1.3028386 0.0001806
## C-A -0.54 -0.9428386 -0.1371614 0.0098537
## C-B -1.44 -1.8428386 -1.0371614 0.0000017
an21<-aov(yield ~ fertilizer, data=df_diet2)
result<- scheffe.test(an21,"fertilizer", alpha = 0.1)
result
## $statistics
##   MSerror Df        F  Mean       CV  Scheffe CriticalDifference
##     0.057 12 2.806796 10.36 2.304505 2.369302          0.3577567
## 
## $parameters
##      test     name.t ntr alpha
##   Scheffe fertilizer   3   0.1
## 
## $means
##   yield       std r        se  Min  Max  Q25  Q50  Q75
## A 10.24 0.2701851 5 0.1067708  9.8 10.5 10.2 10.3 10.4
## B 11.14 0.2701851 5 0.1067708 10.8 11.5 11.0 11.1 11.3
## C  9.70 0.1581139 5 0.1067708  9.5  9.9  9.6  9.7  9.8
## 
## $comparison
## NULL
## 
## $groups
##   yield groups
## B 11.14      a
## A 10.24      b
## C  9.70      c
## 
## attr(,"class")
## [1] "group"

Na základe oboch post testov je štatisticky významný rozdiel medzi hnojivami A, B a C. Hnojivo B je účinnejšie ako hnojivo A a C a hnojivo A je účinnejšie ako hnojivo C.

Dvojfaktorová ANOVA

Zadanie: Predpokladajme, že študujete účinky typu cvičenia a typu stravy na chudnutie. Pre každý faktor máte dve úrovne: Typ cvičenia: Kardio a sila Typ stravy: Nízkosacharidový a nízkotučný Máte podozrenie, že medzi typom cvičenia a typom stravy môže existovať interakcia na chudnutie, t. j. účinnosť cvičenia sa môže líšiť v závislosti od typu dodržiavanej diéty.

weight_loss <- c( 
# Cardio + Low-Carb 
rnorm(10, mean = 5, sd = 1),   
# Cardio + Low-Fat 
rnorm(10, mean = 5, sd = 1), 
# Strength + Low-Carb 
rnorm(10, mean = 5, sd = 1), 
# Strength + Low-Fat 
rnorm(10, mean = 7, sd = 1)   
) 
exercise <- factor(rep(c("Cardio", "Strength"), each = 20)) 
diet <- factor(rep(c("Low-Carb", "Low-Fat"), times = 10, each = 10)) 

df_diet3 <- data.frame(
  weight_loss,
  exercise,
  diet
)
df_diet3$exercise<-as.factor(df_diet3$exercise)
df_diet3$diet<-as.factor(df_diet3$diet)

Urobme si prvotnu grafickú analýzu:

ggplot(df_diet3, aes(x= diet, y= weight_loss))+
  geom_boxplot(aes(fill=diet))+
  scale_fill_brewer(palette = "Greys", name="")+
  labs(x= "Diéta", y= "Сhudnutie [kg]", title = "Porovnanie chudnutia vzhľadom na diéty")+
  theme(plot.title = element_text(hjust=0.5))

ggplot(df_diet3, aes(x= exercise, y= weight_loss))+
  geom_boxplot(aes(fill=exercise))+
  scale_fill_brewer(palette = "Greys", name="")+
  labs(x= "Сvičenie", y= "Сhudnutie [kg]", title = "Porovnanie chudnutia vzhľadom na cvičenia")+
  theme(plot.title = element_text(hjust=0.5))

ggplot(df_diet3, aes(x = diet, y= weight_loss, fill = exercise))+ 
  geom_boxplot()+
  scale_fill_brewer(palette = "Greys", name="Сvičenie")+
  labs(x= "Diéta", y= "Сhudnutie [kg]", title = "Porovnanie možných interakcií")+
  theme(plot.title = element_text(hjust=0.5))

Z boxplotov sa zdá že kombinácia nízkotučnej diety a silových cvičení dáva lepšie chudnutie.

an3<-aov(weight_loss ~ diet + exercise + diet * exercise, df_diet3)
summary(an3)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## diet            1 103.56  103.56  110.31  < 2e-16 ***
## exercise        1  18.28   18.28   19.47 1.69e-05 ***
## diet:exercise   1  61.77   61.77   65.80 5.36e-14 ***
## Residuals     196 184.00    0.94                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Čím je F vyššie, tým silnejšie tento faktor ovplyvňuje konečný výsledok, preto najsilnejšie pôsobí na chudnutie kombinácia stravy a cvičenia.

Na určenie, ktorý plan diety sa štatisticky významne odlišujé, budeme testovať kontrasty (alebo inak post hoc test, viacnásobné porovnávanie priemerov). Použijeme Tukeyov test.

TukeyHSD(an3)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = weight_loss ~ diet + exercise + diet * exercise, data = df_diet3)
## 
## $diet
##                      diff     lwr      upr p adj
## Low-Fat-Low-Carb 1.439157 1.16893 1.709385     0
## 
## $exercise
##                      diff      lwr       upr    p adj
## Strength-Cardio 0.6045807 0.334353 0.8748084 1.69e-05
## 
## $`diet:exercise`
##                                          diff        lwr          upr     p adj
## Low-Fat:Cardio-Low-Carb:Cardio      0.3276983 -0.1744235  0.829820053 0.3311285
## Low-Carb:Strength-Low-Carb:Cardio  -0.5068783 -1.0090000 -0.004756484 0.0469135
## Low-Fat:Strength-Low-Carb:Cardio    2.0437380  1.5416162  2.545859732 0.0000000
## Low-Carb:Strength-Low-Fat:Cardio   -0.8345765 -1.3366983 -0.332454757 0.0001528
## Low-Fat:Strength-Low-Fat:Cardio     1.7160397  1.2139179  2.218161459 0.0000000
## Low-Fat:Strength-Low-Carb:Strength  2.5506162  2.0484944  3.052737996 0.0000000

Vidíme, že najlepšie výsledky prináša nízkotučná dieta v kombinácii so silovými cvičeniami. (má najmenšie p-hodnoty a kladny rozdiel)