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

#use these libraries 
library(tidyverse)  # ggplot(), etc.
## ── Attaching packages ────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.1.3     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## Warning: package 'tidyr' was built under R version 3.6.2
## ── Conflicts ───────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(lmPerm)     # permutation tests: aovp(), lmp(), etc.
library(DescTools)  # for PostHocTests()
## Warning: package 'DescTools' was built under R version 3.6.2
#library(grid)       # for grobbing - in function defintions

DATA The lifetime data is available in the Sleuth3::case0501 dataframe.

case0501 <- Sleuth3::case0501
names(case0501)
## [1] "Lifetime" "Diet"
ggplot(data=case0501, aes(x=Diet, y=Lifetime, fill=Diet)) +
  geom_boxplot(outlier.color="red", outlier.size=2) +
  labs(title="Lifespan Data") +
  theme_classic()

ANALYSIS Are there signficanct differences between caloric restrictions? Which restrictions have the best results?

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
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

The post hoc test shows that N/R50, NP, and R/R50 are almost equally effective.

CHECKING MODEL ASSUMPTIONS

#modeling 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()
}
#modeling homogenity 
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)

Both assumptions are violated.We need to verify the result with two non-parametric tests instead.

A Permuation 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)
## [1] "aovp" "aov"  "lmp"  "lm"
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

The comparisons are in the same order as those in the previous section.

A 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)
## [1] "htest"

Unfortunately, there are no built-in post hoc tests to follow up the Kruskal-Wallis test. So do pairwise tests using a Bonferroni adjustment: for an experiment-wise α=0.05 and (62)=15 possible comparisons, do the pairwise tests with α=0.05/15=0.0033.

So, for example, compare treatments N/R50 and R/R50.

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 0
detach(case0501)

FINDINGS 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