Introduction

Researchers have studied the effect of caloric restriction on lifespan. In one study, female mice were randomly assigned to the following treatment groups:

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

This example uses the lifetime data 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

Do the treatment groups have significantly different mean lifetimes? Are there clear “best” and “worst” diets, or are there others that are “just as good” or “just as bad?”

# 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.2695

Looks like diets N/R50, and R/R50 are equally effective. Both N/R40 and Lopro are closely as effective. NP was the least 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()
}

#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()
}
GraphNormality(ANOVA.model)

GraphHomogeneity(case0501$Lifetime, case0501$Diet, dataset=case0501)

We will verify the result with two non-parametric tests.

Permutation test

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 ' ' 1
class(perm.model) #check if obj has more info
## [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

Rank-Sum test

#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-16
class(kw.test) #check if obj has more info
## [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 \(\binom{6}{2}\) = 15 possible comparisons, do the pairwise tests with \(\alpha\) = \(\frac{0.05}{15}\) = 0.0033.

For example, we will compare treatment groups N/R50 and NP.

attach(case0501)
wilcox.test(Lifetime[Diet=="N/R50"], Lifetime[Diet=="NP"], alternative="two.sided", conf.level=0.9967)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Lifetime[Diet == "N/R50"] and Lifetime[Diet == "NP"]
## W = 3211.5, p-value = 3.939e-15
## alternative hypothesis: true location shift is not equal to 0
detach(case0501)

Finding

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