Three ANOVAs
Sumaiya Shefa
11/14/2020
library(tidyverse) # ggplot(), etc.
library(lmPerm) # permutation tests: aovp(), lmp(), etc.
library(DescTools) # for PostHocTests()
library(kableExtra)Introduction
Researchers have studied the effect of caloric restriction on lifespan.In one study, female mice were randomly assigned to the following treatment groups:
#tabular summary
case0501 <-
data.frame(
Groups = c("NP", "N/N85", "N/R50", "R/R50", "N/R50 lopro", "N/R40"),
Description = linebreak(
c(
"unlimited amounts of a nonpurified, standard diet for lab mice",
"normal feeding before and after weaning ration was controlled at 85 kcal / week after weaning(control group)",
"normal feeding before weaning, reduced-calorie diet of 50 kcal/week after weaning",
"reduced-calorie diet of 50 kcal/week before and after weaning",
"normal feeding before weaning, reduced-calorie diet of 50 cal/week after weaning, and dietary protein reduced with advancing age",
"normal feeding before weaning, severely reducedcalorie diet of 40 kcal/week after weaning"
)
)
)
kable(case0501) %>%
kable_styling(latex_options = "hold_position")| Groups | Description |
|---|---|
| NP | unlimited amounts of a nonpurified, standard diet for lab mice |
| N/N85 | normal feeding before and after weaning ration was controlled at 85 kcal / week after weaning(control group) |
| N/R50 | normal feeding before weaning, reduced-calorie diet of 50 kcal/week after weaning |
| R/R50 | reduced-calorie diet of 50 kcal/week before and after weaning |
| N/R50 lopro | normal feeding before weaning, reduced-calorie diet of 50 cal/week after weaning, and dietary protein reduced with advancing age |
| N/R40 | normal feeding before weaning, severely reducedcalorie diet of 40 kcal/week after weaning |
Data
The lifetime data used in this example is from “The Retardation of Aging in Mice by Dietary Restriction: Longevity, Cancer, immunity, and Lifetime Energy Intake.”
case0501 <- Sleuth3::case0501
names(case0501)## [1] "Lifetime" "Diet"#graphical summary
ggplot(data = case0501, aes(x = Diet, y = Lifetime, fill = Diet)) +
geom_boxplot(outlier.color = "red", outlier.size = 2) +
labs(title = "Lifetime Data") +
theme_classic()
Analysis of Variance
Do the treatment groups have significantly different mean lifetimes? Which treatment groups are different?
#ANOVA model
summary(ANOVA.model <- aov(Lifetime ~ Diet, data = case0501))## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 5 12734 2546.8 57.1 <2e-16 ***
## Residuals 343 15297 44.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Post hoc test
PostHocTest(ANOVA.model, method="scheffe", conf.level=NA, ordered=FALSE)##
## Posthoc multiple comparisons of means: Scheffe Test
##
## $Diet
## N/N85 N/R40 N/R50 NP R/R50
## N/R40 < 2e-16 - - - -
## N/R50 1.1e-11 0.3289 - - -
## NP 0.0063 < 2e-16 < 2e-16 - -
## R/R50 9.7e-12 0.6644 0.9986 < 2e-16 -
## lopro 1.6e-05 0.0022 0.4440 1.4e-15 0.2695The results of the tests shows that N/R50, NP, and R/R50 are almost equally effective.
Checking Model Assumptions
#Normality
GraphNormality <- function(model) {
residual.data <- data.frame(e=model$residuals)
H <- shapiro.test(residual.data$e)
ggplot(residual.data, aes(sample=e)) +
stat_qq() +
geom_abline(color="blue", intercept=mean(residual.data$e), slope=sd(residual.data$e)) +
labs(title="Normally Distributed Residuals?",
subtitle = paste("Shapiro-Wilks test p-value =", signif(H$p.value,5)) ) +
theme_classic()
}
GraphNormality(ANOVA.model)
#Homogeneity of Variance
GraphHomogeneity <- function(response, predictor, dataset=NULL) {
H <- bartlett.test(response~predictor)
ggplot(data=dataset, aes(x=predictor, y=response, fill=predictor)) +
geom_boxplot(outlier.color = "red", outlier.size=3) +
geom_jitter(width=0.2) +
labs( title="Homogenous Variances?",
subtitle=paste("Bartlett's test p-value =", signif(H$p.value,5)) ) +
theme_classic()
}
GraphHomogeneity(case0501$Lifetime, case0501$Diet, dataset=case0501)
The graphs representing the tests show that the assumptions of both normality and homogeneity are violated. To verify the results of these tests, two non-parametric tests will be performed.
Non-parametric Tests
summary(perm.model <- aovp(Lifetime ~ Diet, data=case0501))## [1] "Settings: unique SS "## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## Diet 5 12734 2546.8 5000 < 2.2e-16 ***
## Residuals 343 15297 44.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1class(perm.model)## [1] "aovp" "aov" "lmp" "lm"#Post hoc test
PostHocTest(perm.model, method="scheffe", conf.level=NA, ordered=FALSE)##
## Posthoc multiple comparisons of means: Scheffe Test
##
## $Diet
## N/N85 N/R40 N/R50 NP R/R50
## N/R40 < 2e-16 - - - -
## N/R50 1.1e-11 0.3289 - - -
## NP 0.0063 < 2e-16 < 2e-16 - -
## R/R50 9.7e-12 0.6644 0.9986 < 2e-16 -
## lopro 1.6e-05 0.0022 0.4440 1.4e-15 0.2695#Kruskal-Wallis Test
(kw.test <- kruskal.test(Lifetime ~ Diet, data=case0501))##
## Kruskal-Wallis rank sum test
##
## data: Lifetime by Diet
## Kruskal-Wallis chi-squared = 159.01, df = 5, p-value < 2.2e-16class(kw.test)## [1] "htest"There are no built-in post hoc tests to follow up the Kruskal-Wallis test. So a pairwise tests using Bonferroni adjustment is followed.Using \(\alpha=0.05\) and \({6\choose2}=15\) possible comparisons, do the pairwise tests with \(\alpha=\frac{0.05}{15}=0.0033\).
For example, compare treatment groups N/R40 and NP.
attach(case0501)
wilcox.test(Lifetime[Diet=="N/R40"], Lifetime[Diet=="NP"], alternative="two.sided", conf.level=0.9967)##
## Wilcoxon rank sum test with continuity correction
##
## data: Lifetime[Diet == "N/R40"] and Lifetime[Diet == "NP"]
## W = 2854, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0detach(case0501)Findings
The results of the non-parametric tests confirm the result of the standard ANOVA.
Reference
R. Weindruch, R. L. Walford, S Fligiel, and D. Guthrie, “The Retardation of Aging in Mice by Dietary Restriction: Longevity, Cancer, Immunity, and Lifetime Energy Intake,” Journal of Nutrition 116(4) (1986):641-54