ANOVA Test - Exercise
Elvriana Elvani(20184920012)
November 20, 2020
1 Exercise 1 : One-way ANOVA
In this exercise you have to consider data set named PlantGrowth from R environment. It contains the weight of plants obtained under a control and two different treatment conditions. Please consider the following steps:
- Data Preparation (One-way ANOVA)
- One-way ANOVA Test
- Check Assumptions (One-way ANOVA)
- Pairwise-comparison (One-way ANOVA)
- Tukey HSD (One-way ANOVA)
1.1 Data Preparation
We’ll use the PlantGrowth dataset.
## weight group
## 1 4.17 ctrl
## 2 5.58 ctrl
## 3 5.18 ctrl
## 4 6.11 ctrl
## 5 4.50 ctrl
## 6 4.61 ctrl
## 7 5.17 ctrl
## 8 4.53 ctrl
## 9 5.33 ctrl
## 10 5.14 ctrl## [1] "ctrl" "trt1" "trt2"Before we begin to build the model, it’s important to consider the summary statistics and also identify if there is some outliers. We ca apply box plot methods or implement by using the R function identify_outliers().
##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterlibrary(ggplot2) # for visualization
# Get the summary statistics
PlantGrowth %>%
group_by(group) %>%
get_summary_stats(weight, type = "mean_sd") ## # A tibble: 3 x 5
## group variable n mean sd
## <fct> <chr> <dbl> <dbl> <dbl>
## 1 ctrl weight 10 5.03 0.583
## 2 trt1 weight 10 4.66 0.794
## 3 trt2 weight 10 5.53 0.443We can visualize the data by plotting using boxplots per level of group.
# boxplot of `weight` by `group`
ggplot(PlantGrowth, aes(x = group, y = weight, fill = group)) +
geom_boxplot() +
geom_jitter(shape = 15,
color = "orange",
position = position_jitter(0.21)) +
theme_minimal()
From the boxplot we can assume that each group have a difference of means.
We can identify outliers by identify_outliers() in the R function.
## # A tibble: 2 x 4
## group weight is.outlier is.extreme
## <fct> <dbl> <lgl> <lgl>
## 1 trt1 5.87 TRUE FALSE
## 2 trt1 6.03 TRUE FALSEFrom the table, we know that there were no extreme outliers.
1.2 One-way ANOVA Test
The aov() function is used to obtain the relevant sums of squares. Using the summary() function on the output from aov() creates the desired ANOVA table.
## Coefficient covariances computed by hccm()## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 group 2 27 4.846 0.016 * 0.264## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 3.766 1.8832 4.846 0.0159 *
## Residuals 27 10.492 0.3886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In the table above, the column ges corresponds to the generalized eta squared (effect size). It measures the proportion of the variability in the outcome variable (here plant weight) that can be explained in terms of the predictor (here, group). An effect size of 26% change can be accounted for each group. However, we notice that the p-value of this test is incredibly low, so using any reasonable significance level we would reject the null hypothesis. Thus we believe the group had an effect on plant weight.
Then, we create a dataframe with a row for each group. By predicting on this dataframe, we obtain the sample means of each group.
groups = data.frame(group = unique(PlantGrowth$group))
data.frame(groups, weight = predict(weight_aov, groups))## group weight
## 1 ctrl 5.032
## 2 trt1 4.661
## 3 trt2 5.5261.3 Check Assumptions
The ANOVA test assumes that, the data are normally distributed and the variance across groups are homogeneous. We can check that with some diagnostic plots.
## # A tibble: 3 x 4
## group variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 ctrl weight 0.957 0.747
## 2 trt1 weight 0.930 0.452
## 3 trt2 weight 0.941 0.564Points 4, 15, 17 are detected as outliers, which can severely affect normality and homogeneity of variance. It can be useful to remove outliers to meet the test assumptions. I recommend Levene’s test, which is less sensitive to departures from normal distribution.
## Loading required package: carData## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.1192 0.3412
## 27From the output above, we can see that the p-value is \(>0.05\), which is not significant. This means that, there is not significant difference between variances across groups. Therefore, we can assume the homogeneity of variances in the different treatment groups.
##
## One-way analysis of means (not assuming equal variances)
##
## data: weight and group
## F = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.017391.4 Pairwise-comparison
## # A tibble: 3 x 9
## term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 group ctrl trt1 0 -0.371 -1.06 0.320 0.391 ns
## 2 group ctrl trt2 0 0.494 -0.197 1.19 0.198 ns
## 3 group trt1 trt2 0 0.865 0.174 1.56 0.012 *The output contains the following columns:
estimate: estimate of the difference between means of the two groupsconf.low, conf.high: the lower and the upper end point of the confidence interval at 95% (default)p.adj: p-value after adjustment for the multiple comparisons.
1.5 Tukey HSD
Tukey’s Honest Significance difference can be applied directly to an object which was created using aov(). It will adjust the p-values of the pairwise comparisons of the means to control the FWER, in this case, for 0.05. Notice it also gives confidence intervals for the difference of the means.
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = weight ~ group, data = PlantGrowth)
##
## $group
## diff lwr upr p adj
## trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
## trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960
## trt2-trt1 0.865 0.1737839 1.5562161 0.0120064Based on these results, we see the difference between trt2 and trt1 are significant but no other group differences were statistically significant.
Also, nicely, we can easily produce a plot of these confidence intervals.
2 Exercise 2 : Two-way ANOVA
Please apply Two-way ANOVA procedures to the rats data from the faraway package. There are two factors here: poison and treat. You can use the levels() function to extract the levels of a factor variable. Give your opinions accordingly!
2.1 Data Preparation
We’ll use the rats dataset from faraway package.
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car##
## Attaching package: 'faraway'## The following objects are masked from 'package:car':
##
## logit, vif## time poison treat
## 1 0.31 I A
## 2 0.82 I B
## 3 0.43 I C
## 4 0.45 I D
## 5 0.45 I A
## 6 1.10 I B
## 7 0.45 I C
## 8 0.71 I D
## 9 0.46 I A
## 10 0.88 I B## [1] "I" "II" "III"## [1] "A" "B" "C" "D"Before we begin to build the model, it’s important to consider the summary statistics and also identify if there is some outliers. We can apply box plot methods or implement by using the R function identify_outliers().
library(tidyr) # to create tidy data
library(rstatix) # use for basic Anova/statistical tests
library(ggplot2) # for visualization
# Get the summary statistics
rats %>%
group_by(poison, treat) %>%
get_summary_stats(time, type = "mean_sd") ## # A tibble: 12 x 6
## poison treat variable n mean sd
## <fct> <fct> <chr> <dbl> <dbl> <dbl>
## 1 I A time 4 0.412 0.069
## 2 I B time 4 0.88 0.161
## 3 I C time 4 0.568 0.157
## 4 I D time 4 0.61 0.113
## 5 II A time 4 0.32 0.075
## 6 II B time 4 0.815 0.336
## 7 II C time 4 0.375 0.057
## 8 II D time 4 0.667 0.271
## 9 III A time 4 0.21 0.022
## 10 III B time 4 0.335 0.047
## 11 III C time 4 0.235 0.013
## 12 III D time 4 0.325 0.026We can visualize the data by plotting using boxplots per level of treat.
# boxplot of `time` by `treat and poison`
ggplot(rats, aes(x = treat, y = time, fill = poison)) +
geom_boxplot() +
geom_jitter(shape = 15,
color = "orange",
position = position_jitter(0.21)) +
theme_minimal()
From the boxplot we can assume that each group have a difference of means.
We can visualize too by interaction plot.
The interaction plot suggests that there might be an interaction between treat and poison as the effect of poison is less pronounced in A.
We can identify outliers by identify_outliers() in the R function.
## # A tibble: 1 x 5
## poison treat time is.outlier is.extreme
## <fct> <fct> <dbl> <lgl> <lgl>
## 1 I A 0.31 TRUE FALSEFrom the table, we know that there were no extreme outliers.
2.2 Two-way ANOVA Test
rats_int = aov(time ~ poison * treat, data = rats) # interaction model
rats_add = aov(time ~ poison + treat, data = rats) # additive model
rats_poison= aov(time ~ poison, data = rats) # single factor model
rats_treat = aov(time ~ treat, data = rats) # single factor model
rats_null = aov(time ~ 1, data = rats) # null model
summary(rats_int)## Df Sum Sq Mean Sq F value Pr(>F)
## poison 2 1.0330 0.5165 23.222 3.33e-07 ***
## treat 3 0.9212 0.3071 13.806 3.78e-06 ***
## poison:treat 6 0.2501 0.0417 1.874 0.112
## Residuals 36 0.8007 0.0222
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Df Sum Sq Mean Sq F value Pr(>F)
## poison 2 1.0330 0.5165 20.64 5.7e-07 ***
## treat 3 0.9212 0.3071 12.27 6.7e-06 ***
## Residuals 42 1.0509 0.0250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Df Sum Sq Mean Sq F value Pr(>F)
## poison 2 1.033 0.5165 11.79 7.66e-05 ***
## Residuals 45 1.972 0.0438
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Df Sum Sq Mean Sq F value Pr(>F)
## treat 3 0.9212 0.30707 6.484 0.000992 ***
## Residuals 44 2.0839 0.04736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 47 3.005 0.06394Using a significance level of \(\alpha=0.05\), we see that either interaction and additive models are significantly different (we reject the Null hypothesis). Noted that the one-way ANOVA model for variable gender is significantly each group has not a difference of means.
## Coefficient covariances computed by hccm()## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 poison 2 36 23.222 3.33e-07 * 0.563
## 2 treat 3 36 13.806 3.78e-06 * 0.535
## 3 poison:treat 6 36 1.874 1.12e-01 0.238There is a statistically significant interaction between poison and treat for rats time, F(6, 36) = 1.874, p = 1.12e-01.
2.3 Check Assumptions
The ANOVA test assumes that, the data are normally distributed and the variance across groups are homogeneous. We can check that with some diagnostic plots.
## # A tibble: 12 x 5
## poison treat variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 I A time 0.782 0.0741
## 2 I B time 0.947 0.700
## 3 I C time 0.895 0.405
## 4 I D time 0.899 0.427
## 5 II A time 0.971 0.848
## 6 II B time 0.948 0.701
## 7 II C time 0.983 0.921
## 8 II D time 0.981 0.907
## 9 III A time 0.927 0.577
## 10 III B time 0.831 0.171
## 11 III C time 0.993 0.972
## 12 III D time 0.946 0.689Points 24, 26, 30 are detected as outliers, which can severely affect normality and homogeneity of variance. It can be useful to remove outliers to meet the test assumptions. I recommend Levene’s test, which is less sensitive to departures from normal distribution. The score were normally distributed (p > 0.05) for each cell, as assessed by Shapiro-Wilk’s test of normality.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 4.1323 0.0005833 ***
## 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 11 36 4.13 0.000583##
## One-way analysis of means (not assuming equal variances)
##
## data: time and poison * treat
## F = 14.848, num df = 11.000, denom df = 13.801, p-value = 8.107e-06From the output above, we can see that the p-value is \(<0.05\), which is significant. This means that, there is significant difference between variances across groups. Therefore, we can assume the homogeneity of variances in the different treatment groups.
2.4 Pairwise-comparison
library(emmeans)
pwc2 <- rats %>%
group_by(poison) %>%
emmeans_test(time ~ treat, p.adjust.method = "bonferroni")
head(pwc2, 10) # 10 of 18 rows## # A tibble: 10 x 10
## poison term .y. group1 group2 df statistic p p.adj p.adj.signif
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 I treat time A B 36 -4.43 8.37e-5 5.02e-4 ***
## 2 I treat time A C 36 -1.47 1.50e-1 9.02e-1 ns
## 3 I treat time A D 36 -1.87 6.92e-2 4.15e-1 ns
## 4 I treat time B C 36 2.96 5.37e-3 3.22e-2 *
## 5 I treat time B D 36 2.56 1.48e-2 8.88e-2 ns
## 6 I treat time C D 36 -0.403 6.89e-1 1.00e+0 ns
## 7 II treat time A B 36 -4.69 3.81e-5 2.29e-4 ***
## 8 II treat time A C 36 -0.522 6.05e-1 1.00e+0 ns
## 9 II treat time A D 36 -3.30 2.22e-3 1.33e-2 *
## 10 II treat time B C 36 4.17 1.82e-4 1.09e-3 **All pairwise comparisons were analyzed between the different treat groups organized by poison. There was a significant difference of rats time between all groups for both I, II and III (p < 0.05).
2.5 Tukey HSD
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = time ~ poison * treat, data = rats)
##
## $poison
## diff lwr upr p adj
## II-I -0.073125 -0.2020091 0.05575913 0.3583151
## III-I -0.341250 -0.4701341 -0.21236587 0.0000005
## III-II -0.268125 -0.3970091 -0.13924087 0.0000339
##
## $treat
## diff lwr upr p adj
## B-A 0.36250000 0.19852116 0.52647884 0.0000047
## C-A 0.07833333 -0.08564550 0.24231217 0.5772283
## D-A 0.22000000 0.05602116 0.38397884 0.0048556
## C-B -0.28416667 -0.44814550 -0.12018783 0.0002333
## D-B -0.14250000 -0.30647884 0.02147884 0.1077087
## D-C 0.14166667 -0.02231217 0.30564550 0.1107678
##
## $`poison:treat`
## diff lwr upr p adj
## II:A-I:A -0.0925 -0.460578867 0.27557887 0.9989739
## III:A-I:A -0.2025 -0.570578867 0.16557887 0.7395476
## I:B-I:A 0.4675 0.099421133 0.83557887 0.0041387
## II:B-I:A 0.4025 0.034421133 0.77057887 0.0220694
## III:B-I:A -0.0775 -0.445578867 0.29057887 0.9998038
## I:C-I:A 0.1550 -0.213078867 0.52307887 0.9393871
## II:C-I:A -0.0375 -0.405578867 0.33057887 0.9999999
## III:C-I:A -0.1775 -0.545578867 0.19057887 0.8641326
## I:D-I:A 0.1975 -0.170578867 0.56557887 0.7673377
## II:D-I:A 0.2550 -0.113078867 0.62307887 0.4203275
## III:D-I:A -0.0875 -0.455578867 0.28057887 0.9993833
## III:A-II:A -0.1100 -0.478078867 0.25807887 0.9953340
## I:B-II:A 0.5600 0.191921133 0.92807887 0.0003212
## II:B-II:A 0.4950 0.126921133 0.86307887 0.0019661
## III:B-II:A 0.0150 -0.353078867 0.38307887 1.0000000
## I:C-II:A 0.2475 -0.120578867 0.61557887 0.4645761
## II:C-II:A 0.0550 -0.313078867 0.42307887 0.9999936
## III:C-II:A -0.0850 -0.453078867 0.28307887 0.9995291
## I:D-II:A 0.2900 -0.078078867 0.65807887 0.2439183
## II:D-II:A 0.3475 -0.020578867 0.71557887 0.0790981
## III:D-II:A 0.0050 -0.363078867 0.37307887 1.0000000
## I:B-III:A 0.6700 0.301921133 1.03807887 0.0000139
## II:B-III:A 0.6050 0.236921133 0.97307887 0.0000893
## III:B-III:A 0.1250 -0.243078867 0.49307887 0.9868993
## I:C-III:A 0.3575 -0.010578867 0.72557887 0.0634939
## II:C-III:A 0.1650 -0.203078867 0.53307887 0.9105785
## III:C-III:A 0.0250 -0.343078867 0.39307887 1.0000000
## I:D-III:A 0.4000 0.031921133 0.76807887 0.0234641
## II:D-III:A 0.4575 0.089421133 0.82557887 0.0054007
## III:D-III:A 0.1150 -0.253078867 0.48307887 0.9932581
## II:B-I:B -0.0650 -0.433078867 0.30307887 0.9999650
## III:B-I:B -0.5450 -0.913078867 -0.17692113 0.0004903
## I:C-I:B -0.3125 -0.680578867 0.05557887 0.1618705
## II:C-I:B -0.5050 -0.873078867 -0.13692113 0.0014940
## III:C-I:B -0.6450 -1.013078867 -0.27692113 0.0000285
## I:D-I:B -0.2700 -0.638078867 0.09807887 0.3379251
## II:D-I:B -0.2125 -0.580578867 0.15557887 0.6808722
## III:D-I:B -0.5550 -0.923078867 -0.18692113 0.0003699
## III:B-II:B -0.4800 -0.848078867 -0.11192113 0.0029569
## I:C-II:B -0.2475 -0.615578867 0.12057887 0.4645761
## II:C-II:B -0.4400 -0.808078867 -0.07192113 0.0085466
## III:C-II:B -0.5800 -0.948078867 -0.21192113 0.0001822
## I:D-II:B -0.2050 -0.573078867 0.16307887 0.7252287
## II:D-II:B -0.1475 -0.515578867 0.22057887 0.9563176
## III:D-II:B -0.4900 -0.858078867 -0.12192113 0.0022537
## I:C-III:B 0.2325 -0.135578867 0.60057887 0.5569107
## II:C-III:B 0.0400 -0.328078867 0.40807887 0.9999998
## III:C-III:B -0.1000 -0.468078867 0.26807887 0.9979417
## I:D-III:B 0.2750 -0.093078867 0.64307887 0.3126107
## II:D-III:B 0.3325 -0.035578867 0.70057887 0.1086625
## III:D-III:B -0.0100 -0.378078867 0.35807887 1.0000000
## II:C-I:C -0.1925 -0.560578867 0.17557887 0.7938454
## III:C-I:C -0.3325 -0.700578867 0.03557887 0.1086625
## I:D-I:C 0.0425 -0.325578867 0.41057887 0.9999996
## II:D-I:C 0.1000 -0.268078867 0.46807887 0.9979417
## III:D-I:C -0.2425 -0.610578867 0.12557887 0.4949176
## III:C-II:C -0.1400 -0.508078867 0.22807887 0.9695827
## I:D-II:C 0.2350 -0.133078867 0.60307887 0.5413031
## II:D-II:C 0.2925 -0.075578867 0.66057887 0.2335610
## III:D-II:C -0.0500 -0.418078867 0.31807887 0.9999976
## I:D-III:C 0.3750 0.006921133 0.74307887 0.0426208
## II:D-III:C 0.4325 0.064421133 0.80057887 0.0103742
## III:D-III:C 0.0900 -0.278078867 0.45807887 0.9992007
## II:D-I:D 0.0575 -0.310578867 0.42557887 0.9999899
## III:D-I:D -0.2850 -0.653078867 0.08307887 0.2655759
## III:D-II:D -0.3425 -0.710578867 0.02557887 0.08807633 Exercise 3 : Three-way ANOVA
In this exercise you will use the headache dataset [in datarium package], which contains the measures of migraine headache episode pain score in 72 participants treated with three different treatments. The participants include 36 males and 36 females. Males and females were further subdivided into whether they were at low or high risk of migraine. You have to understand how each independent variable (type of treatments, risk of migraine and gender) interact to predict the pain score.
3.1 Data Preparation
We’ll use the headache dataset from datarium package.
## # A tibble: 10 x 5
## id gender risk treatment pain_score
## <int> <fct> <fct> <fct> <dbl>
## 1 1 male low X 79.3
## 2 2 male low X 76.8
## 3 3 male low X 70.8
## 4 4 male low X 81.2
## 5 5 male low X 75.1
## 6 6 male low X 73.1
## 7 7 male low Y 68.2
## 8 8 male low Y 80.7
## 9 9 male low Y 75.3
## 10 10 male low Y 73.4## [1] "I" "II" "III"## [1] "A" "B" "C" "D"Before we begin to build the model, it’s important to consider the summary statistics and also identify if there is some outliers.
# Get the summary statistics
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.72We can visualize the data by plotting using boxplots.
library(ggpubr)
ggboxplot(headache, x = "treatment", y = "pain_score", color = "risk", palette="jco", facet.by="gender") 
We can identify outliers by identify_outliers() in the R function.
# Identify outliers
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 FALSEFrom the table, we know that there was one extreme outliers (id = 57).
3.2 Check Assumptions
We’ll check the assumption by analizing the model residuals. This is use QQ Plot and Shapiro-Wilk test.
model <- lm(pain_score ~ gender * risk * treatment, data = headache)
ggqqplot(residuals(model)) # Create a QQ plot of residuals
## # A tibble: 1 x 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 residuals(model) 0.982 0.398In the QQ Plot and Shapiro-Wilk test, the p-value isn’t significant (0.398). So, we can assume that is normality.
We can check the assumption by groups too.
## # 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.801From that, we know that the pain scores were normally distributed (p > 0.05), except for one group (female at high risk of migraine taking drug X, p = 0.0086), as assessed by Shapiro-Wilk’s test of normality.
Now, we create QQ plot for each cell of design :
ggqqplot(headache, "pain_score", ggtheme = theme_bw()) +
facet_grid(gender + risk ~ treatment, labeller = "label_both")
From that, all the points fall approximately along the reference line, except for one group (female at high risk of migraine taking drug X).
Then, we can check the Homogneity of variance assumption using the Levene’s test.
## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 11 60 0.179 0.9983.3 Three-way ANOVA Test
## 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.198There is a statistically significant three-way interaction between gender, risk and treatment, F(2, 60) = 7.406, p = 1.00e-03.
3.4 Pairwise-comparison
library(emmeans)
pwc3 <- headache %>%
group_by(gender, risk) %>%
emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni") %>%
select(-df, -statistic, -p) # Remove details
# Show comparison results for male at high risk
pwc3 %>%
filter(gender == "male", risk == "high")## # A tibble: 3 x 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# Estimated marginal means with 95% confidence interval
get_emmeans(pwc3) %>%
filter(gender == "male", risk == "high")## # A tibble: 3 x 9
## gender risk treatment emmean se df conf.low conf.high method
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 male high X 92.7 1.80 60 89.1 96.3 Emmeans test
## 2 male high Y 82.3 1.80 60 78.7 85.9 Emmeans test
## 3 male high Z 79.7 1.80 60 76.1 83.3 Emmeans testFrom that, we know for male at high risk, there is a statistically significant mean difference between treatment X and Y of 10.4 (p.adj < 0.001), and between treatment X and Z of 13.1 (p.adj < 0.0001). However, the difference between treatment Y and Z (2.66) is not statistically significant, p.adj = 0.897.