Anova tiga arah
Yudhoro Kunto Wibisono
2024-10-14
#Descriptive Outlier
library(readr)
data <- read.csv("D:/download/data_fix.csv")
data## gender Diet weight weight6weeks
## 1 M B a 60.0
## 2 M B c 103.0
## 3 F A a 54.2
## 4 F A a 54.0
## 5 F A a 63.3
## 6 F A a 61.1
## 7 F A a 62.2
## 8 F A a 64.0
## 9 F A a 65.0
## 10 F A a 60.5
## 11 F A a 68.1
## 12 F A a 66.9
## 13 F A b 70.5
## 14 F A a 69.0
## 15 F A a 68.4
## 16 F A b 81.1
## 17 F B a 60.1
## 18 F B a 56.0
## 19 F B a 57.3
## 20 F B a 56.7
## 21 F B a 55.0
## 22 F B a 62.4
## 23 F B a 60.3
## 24 F B a 59.4
## 25 F B a 62.0
## 26 F B a 64.0
## 27 F B a 63.8
## 28 F B a 63.3
## 29 F B b 72.7
## 30 F B b 77.5
## 31 F C a 53.0
## 32 F C a 56.4
## 33 F C a 60.6
## 34 F C a 58.2
## 35 F C a 58.2
## 36 F C a 61.6
## 37 F C a 60.2
## 38 F C a 61.8
## 39 F C a 63.0
## 40 F C a 62.7
## 41 F C b 71.1
## 42 F C a 64.4
## 43 F C a 68.9
## 44 F C a 68.7
## 45 F C b 71.0
## 46 M A b 71.6
## 47 M A b 70.9
## 48 M A a 69.5
## 49 M A b 73.9
## 50 M A b 71.0
## 51 M A b 77.6
## 52 M A b 79.1
## 53 M A b 81.5
## 54 M A b 81.9
## 55 M A b 84.5
## 56 M B a 66.8
## 57 M B b 72.6
## 58 M B a 69.2
## 59 M B b 72.5
## 60 M B b 72.7
## 61 M B b 76.3
## 62 M B b 73.6
## 63 M B b 72.9
## 64 M B b 71.1
## 65 M B b 81.4
## 66 M B b 75.7
## 67 M C a 68.5
## 68 M C b 72.1
## 69 M C b 72.5
## 70 M C b 77.5
## 71 M C b 75.2
## 72 M C a 69.4
## 73 M C b 74.5
## 74 M C b 80.2
## 75 M C b 79.9
## 76 M C b 79.7
## 77 M C b 77.8
## 78 M C b 81.9str(data)## 'data.frame': 78 obs. of 4 variables:
## $ gender : chr "M" "M" "F" "F" ...
## $ Diet : chr "B" "B" "A" "A" ...
## $ weight : chr "a" "c" "a" "a" ...
## $ weight6weeks: num 60 103 54.2 54 63.3 61.1 62.2 64 65 60.5 ...library(ggpubr)## Loading required package: ggplot2ggboxplot(data, x = "Diet", y = "weight6weeks", color = "weight", facet.by = "gender")
ms <- data %>%
group_by(gender, Diet, weight) %>%
get_summary_stats(weight6weeks, type = "mean_sd")
library(kableExtra)
library(formattable)
ms$sd <- color_bar("lightgreen")(ms$sd)
ms$mean <- color_tile("yellow", "red")(ms$mean)
ms$gender <- ifelse(
ms$gender == "male",
cell_spec(ms$gender, color = "steelblue", bold = T),
cell_spec(ms$gender, color = "hotpink", bold = T)
)
ms %>%
kbl(escape = F) %>%
kable_material_dark()| gender | Diet | weight | variable | n | mean | sd |
|---|---|---|---|---|---|---|
| F | A | a | weight6weeks | 12 | 63.058 | 5.049 |
| F | A | b | weight6weeks | 2 | 75.800 | 7.495 |
| F | B | a | weight6weeks | 12 | 60.025 | 3.173 |
| F | B | b | weight6weeks | 2 | 75.100 | 3.394 |
| F | C | a | weight6weeks | 13 | 61.362 | 4.482 |
| F | C | b | weight6weeks | 2 | 71.050 | 0.071 |
| M | A | a | weight6weeks | 1 | 69.500 | NA |
| M | A | b | weight6weeks | 9 | 76.889 | 5.210 |
| M | B | a | weight6weeks | 3 | 65.333 | 4.772 |
| M | B | b | weight6weeks | 9 | 74.311 | 3.117 |
| M | B | c | weight6weeks | 1 | 103.000 | NA |
| M | C | a | weight6weeks | 2 | 68.950 | 0.636 |
| M | C | b | weight6weeks | 10 | 77.130 | 3.406 |
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterlibrary(kableExtra)
tab <- data %>%
group_by(gender, Diet, weight) %>%
identify_outliers(weight6weeks)
tab$is.extreme <- cell_spec(tab$is.extreme, color = ifelse(tab$is.extreme == TRUE, "red", "white"))
tab %>%
kbl(caption = "Outliers identification", escape = F) %>%
kable_material_dark("striped")| gender | Diet | weight | weight6weeks | is.outlier | is.extreme |
|---|---|---|---|---|---|
| M | B | b | 81.4 | TRUE | FALSE |
#Normality
data## gender Diet weight weight6weeks
## 1 M B a 60.0
## 2 M B c 103.0
## 3 F A a 54.2
## 4 F A a 54.0
## 5 F A a 63.3
## 6 F A a 61.1
## 7 F A a 62.2
## 8 F A a 64.0
## 9 F A a 65.0
## 10 F A a 60.5
## 11 F A a 68.1
## 12 F A a 66.9
## 13 F A b 70.5
## 14 F A a 69.0
## 15 F A a 68.4
## 16 F A b 81.1
## 17 F B a 60.1
## 18 F B a 56.0
## 19 F B a 57.3
## 20 F B a 56.7
## 21 F B a 55.0
## 22 F B a 62.4
## 23 F B a 60.3
## 24 F B a 59.4
## 25 F B a 62.0
## 26 F B a 64.0
## 27 F B a 63.8
## 28 F B a 63.3
## 29 F B b 72.7
## 30 F B b 77.5
## 31 F C a 53.0
## 32 F C a 56.4
## 33 F C a 60.6
## 34 F C a 58.2
## 35 F C a 58.2
## 36 F C a 61.6
## 37 F C a 60.2
## 38 F C a 61.8
## 39 F C a 63.0
## 40 F C a 62.7
## 41 F C b 71.1
## 42 F C a 64.4
## 43 F C a 68.9
## 44 F C a 68.7
## 45 F C b 71.0
## 46 M A b 71.6
## 47 M A b 70.9
## 48 M A a 69.5
## 49 M A b 73.9
## 50 M A b 71.0
## 51 M A b 77.6
## 52 M A b 79.1
## 53 M A b 81.5
## 54 M A b 81.9
## 55 M A b 84.5
## 56 M B a 66.8
## 57 M B b 72.6
## 58 M B a 69.2
## 59 M B b 72.5
## 60 M B b 72.7
## 61 M B b 76.3
## 62 M B b 73.6
## 63 M B b 72.9
## 64 M B b 71.1
## 65 M B b 81.4
## 66 M B b 75.7
## 67 M C a 68.5
## 68 M C b 72.1
## 69 M C b 72.5
## 70 M C b 77.5
## 71 M C b 75.2
## 72 M C a 69.4
## 73 M C b 74.5
## 74 M C b 80.2
## 75 M C b 79.9
## 76 M C b 79.7
## 77 M C b 77.8
## 78 M C b 81.9library(ggpubr)
model <- lm(weight6weeks ~ gender * weight, data = data)
ggqqplot(residuals(model))
library(stats)
shapiro.test(residuals(model))##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.9817, p-value = 0.325data## gender Diet weight weight6weeks
## 1 M B a 60.0
## 2 M B c 103.0
## 3 F A a 54.2
## 4 F A a 54.0
## 5 F A a 63.3
## 6 F A a 61.1
## 7 F A a 62.2
## 8 F A a 64.0
## 9 F A a 65.0
## 10 F A a 60.5
## 11 F A a 68.1
## 12 F A a 66.9
## 13 F A b 70.5
## 14 F A a 69.0
## 15 F A a 68.4
## 16 F A b 81.1
## 17 F B a 60.1
## 18 F B a 56.0
## 19 F B a 57.3
## 20 F B a 56.7
## 21 F B a 55.0
## 22 F B a 62.4
## 23 F B a 60.3
## 24 F B a 59.4
## 25 F B a 62.0
## 26 F B a 64.0
## 27 F B a 63.8
## 28 F B a 63.3
## 29 F B b 72.7
## 30 F B b 77.5
## 31 F C a 53.0
## 32 F C a 56.4
## 33 F C a 60.6
## 34 F C a 58.2
## 35 F C a 58.2
## 36 F C a 61.6
## 37 F C a 60.2
## 38 F C a 61.8
## 39 F C a 63.0
## 40 F C a 62.7
## 41 F C b 71.1
## 42 F C a 64.4
## 43 F C a 68.9
## 44 F C a 68.7
## 45 F C b 71.0
## 46 M A b 71.6
## 47 M A b 70.9
## 48 M A a 69.5
## 49 M A b 73.9
## 50 M A b 71.0
## 51 M A b 77.6
## 52 M A b 79.1
## 53 M A b 81.5
## 54 M A b 81.9
## 55 M A b 84.5
## 56 M B a 66.8
## 57 M B b 72.6
## 58 M B a 69.2
## 59 M B b 72.5
## 60 M B b 72.7
## 61 M B b 76.3
## 62 M B b 73.6
## 63 M B b 72.9
## 64 M B b 71.1
## 65 M B b 81.4
## 66 M B b 75.7
## 67 M C a 68.5
## 68 M C b 72.1
## 69 M C b 72.5
## 70 M C b 77.5
## 71 M C b 75.2
## 72 M C a 69.4
## 73 M C b 74.5
## 74 M C b 80.2
## 75 M C b 79.9
## 76 M C b 79.7
## 77 M C b 77.8
## 78 M C b 81.9#Homogenitas
data %>%
levene_test(weight6weeks~gender*Diet*weight) %>%
kbl(caption = "Levene's test") %>%
kable_material_dark()| df1 | df2 | statistic | p |
|---|---|---|---|
| 12 | 65 | 1.33497 | 0.2211714 |
##Fit linear model and plot
plot(lm(weight6weeks~gender*Diet*weight, data=data),1)
##Anova
library(formattable)
library(kableExtra)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(rstatix)
library(report)res.aov <- data %>% anova_test(weight6weeks ~ gender*Diet*weight)
res.aov %>% kbl(caption = "ANOVA table") %>%
kable_material_dark() %>%
row_spec(c(7), color ="red", bold = T)| Effect | DFn | DFd | F | p | p<.05 | ges |
|---|---|---|---|---|---|---|
| gender | 1 | 65 | 10.239 | 0.002 |
|
0.136 |
| Diet | 2 | 65 | 2.776 | 0.070 | 0.079 | |
| weight | 2 | 65 | 56.401 | 0.000 |
|
0.634 |
| gender:Diet | 2 | 65 | 1.037 | 0.360 | 0.031 | |
| gender:weight | 1 | 65 | 2.531 | 0.116 | 0.037 | |
| Diet:weight | 2 | 65 | 0.468 | 0.628 | 0.014 | |
| gender:Diet:weight | 2 | 65 | 0.296 | 0.745 | 0.009 |
results <- aov(weight6weeks ~ gender*Diet*weight, data = data)
anova(results) %>% kbl %>% kable_material_dark()| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 2854.20155 | 2854.201553 | 162.0065705 | 0.0000000 |
| Diet | 2 | 67.01552 | 33.507760 | 1.9019250 | 0.1574962 |
| weight | 2 | 1974.09789 | 987.048944 | 56.0256210 | 0.0000000 |
| gender:Diet | 2 | 23.79940 | 11.899701 | 0.6754358 | 0.5124701 |
| gender:weight | 1 | 41.60905 | 41.609054 | 2.3617604 | 0.1291957 |
| Diet:weight | 2 | 16.50181 | 8.250904 | 0.4683273 | 0.6281445 |
| gender:Diet:weight | 2 | 10.41805 | 5.209024 | 0.2956680 | 0.7450297 |
| Residuals | 65 | 1145.15788 | 17.617814 | NA | NA |
report(results)## The ANOVA (formula: weight6weeks ~ gender * Diet * weight) suggests that:
##
## - The main effect of gender is statistically significant and large (F(1, 65) =
## 162.01, p < .001; Eta2 (partial) = 0.71, 95% CI [0.62, 1.00])
## - The main effect of Diet is statistically not significant and small (F(2, 65)
## = 1.90, p = 0.157; Eta2 (partial) = 0.06, 95% CI [0.00, 1.00])
## - The main effect of weight is statistically significant and large (F(2, 65) =
## 56.03, p < .001; Eta2 (partial) = 0.63, 95% CI [0.51, 1.00])
## - The interaction between gender and Diet is statistically not significant and
## small (F(2, 65) = 0.68, p = 0.512; Eta2 (partial) = 0.02, 95% CI [0.00, 1.00])
## - The interaction between gender and weight is statistically not significant
## and small (F(1, 65) = 2.36, p = 0.129; Eta2 (partial) = 0.04, 95% CI [0.00,
## 1.00])
## - The interaction between Diet and weight is statistically not significant and
## small (F(2, 65) = 0.47, p = 0.628; Eta2 (partial) = 0.01, 95% CI [0.00, 1.00])
## - The interaction between gender, Diet and weight is statistically not
## significant and very small (F(2, 65) = 0.30, p = 0.745; Eta2 (partial) =
## 9.02e-03, 95% CI [0.00, 1.00])
##
## Effect sizes were labelled following Field's (2013) recommendations.#Post hoc test
model <- lm(weight6weeks ~ gender*Diet*weight, data = data)
data %>%
group_by(gender) %>%
anova_test(weight6weeks ~ Diet*weight, error = model) %>%
kbl(caption = "Two-way interactions") %>%
kable_material_dark()| gender | Effect | DFn | DFd | F | p | p<.05 | ges |
|---|---|---|---|---|---|---|---|
| F | Diet | 2 | 65 | 1.607 | 0.208 | 0.047 | |
| F | weight | 1 | 65 | 45.717 | 0.000 |
|
0.413 |
| F | Diet:weight | 2 | 65 | 0.714 | 0.493 | 0.022 | |
| M | Diet | 2 | 28 | 2.046 | 0.148 | 0.127 | |
| M | weight | 2 | 28 | 38.469 | 0.000 |
|
0.733 |
| M | Diet:weight | 2 | 28 | 0.055 | 0.947 | 0.004 |
weight.effect <- data %>%
group_by(gender, Diet) %>%
anova_test(weight6weeks ~ weight, error = model)
weight.effect %>% kbl(caption = "weight effect") %>%
kable_material_dark() %>%
row_spec(c(1,2), color ="red", bold = T)| gender | Diet | Effect | DFn | DFd | F | p | p<.05 | ges |
|---|---|---|---|---|---|---|---|---|
| F | A | weight | 1 | 65 | 15.797 | 1.79e-04 |
|
0.196 |
| F | B | weight | 1 | 65 | 22.113 | 1.38e-05 |
|
0.254 |
| F | C | weight | 1 | 65 | 9.235 | 3.00e-03 |
|
0.124 |
| M | A | weight | 1 | 65 | 2.789 | 1.00e-01 | 0.041 | |
| M | B | weight | 2 | 65 | 30.214 | 0.00e+00 |
|
0.482 |
| M | C | weight | 1 | 65 | 6.330 | 1.40e-02 |
|
0.089 |
library(ggplot2)
library(kableExtra)
library(dplyr)
library(rstatix)
library(report)
pwc <- data %>%
group_by(gender, weight) %>%
emmeans_test(weight6weeks ~ Diet, p.adjust.method = "bonferroni") %>%
select(-df, -statistic, -p) # Remove details
# Show comparison results for male at high weight
pwc %>% filter(gender == "male", weight == "96-110") %>%
kbl(caption = "Pairwise comparisons") %>%
kable_material_dark()| gender | weight | term | .y. | group1 | group2 | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|