This is an update of Paired t-test as a special case of linear model and hierarchical model

Figure 2A of the paper Meta-omics analysis of elite athletes identifies a performance-enhancing microbe that functions via lactate metabolism uses a paired t-test to compare endurance performance in mice treated with a control microbe (Lactobacillus bulgaricus) and a test microbe (Veillonella atypica) in a cross-over design (so each mouse was treated with both bacteria). The data are in the “Supplementary Tables” Excel file, which includes the raw data for the main paper.

Endurance performance of each mouse was measured three times under each treatment level. The authors used two different analyses of the data: 1) a paired t-test between the treatments using the maximum performance over the three trials, 2) a linear mixed model with mouse ID as a random intercept (and sequence and day as covariates). A problem with the estimation of a treatment effect using the maximum response (the authors used the maximum of three endurance trials) is that the variance of the estimate of a maximum is bigger than the variance of the estimate of a mean. Consequently, the variance of the effect will be bigger, with more inflated large effects (Type M error) and more effects with the wrong sign (Type S error).

Here I repeat the analysis from the earlier post (which focussed on the paired t-test as a special case of a linear model) but add the analysis of the means and a simulation to show the effect on the variance of the estimated treatment effect.



bookdown_it <- TRUE
  data_path <- "../data"
  out_path <- "../output"
  data_path <- "../content/data"
  out_path <- "../content/output"


folder <- "Data from Meta-omics analysis of elite athletes identifies a performance-enhancing microbe that functions via lactate metabolism"
file_i <- "41591_2019_485_MOESM2_ESM.xlsx"
file_path <- paste(data_path, folder, file_i, sep="/")
sheet_i <- "ST3 Veillonella run times"
range_vec <- c("a5:d21", "f5:i21", "a24:d40", "f24:i40")
fig2a.orig <- data.table(NULL)
poop <- c("Lb", "Va")
for(week_i in 1:4){
  range_i <- range_vec[week_i]
  if(week_i %in% c(1, 3)){
    treatment <- rep(c(poop[1], poop[2]), each=8)
    treatment <- rep(c(poop[2], poop[1]), each=8)
  fig2a.orig <- rbind(fig2a.orig, data.table(week=week_i, treatment=treatment, data.table(
    read_excel(file_path, sheet=sheet_i, range=range_i))))

# clean column names
setnames(fig2a.orig, old="...1", new="ID")
setnames(fig2a.orig, old=colnames(fig2a.orig), new=clean_label(colnames(fig2a.orig)))

# add treatment sequence
fig2a.orig[, sequence:=rep(rep(rep(c("LLLVVV", "VVVLLL"), each=8), 2), 2)]

# wide to long
fig2a <- melt(fig2a.orig, id.vars=c("week", "ID", "treatment", "sequence"),

Because of the cross-over design, the ID is not in order within each level of treatment. Make a wide data table with new treatment columns matched by ID.

# get max for each week x ID x treatment combo
fig2a.max <- fig2a[, .(time_max=max(time)), by=.(week, ID, treatment)]

# match the two treatments applied to each ID
va <- fig2a.max[treatment=="Va", ]
lb <- fig2a.max[treatment=="Lb", ]
fig2a_wide <- merge(va, lb, by="ID")


paired t-test

# paired t
res.t <- t.test(fig2a_wide[, time_max.x], fig2a_wide[, time_max.y], paired=TRUE)
##  Paired t-test
## data:  fig2a_wide[, time_max.x] and fig2a_wide[, time_max.y]
## t = 2.4117, df = 31, p-value = 0.02199
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   20.47859 244.89641
## sample estimates:
## mean of the differences 
##                132.6875

linear model

# as a linear model
y <- fig2a_wide[, time_max.x] - fig2a_wide[, time_max.y] # Va - Lb
res.lm <- lm(y ~ 1)
##             Estimate Std. Error  t value   Pr(>|t|)
## (Intercept) 132.6875   55.01749 2.411733 0.02198912

hierarchical (linear mixed) model with random intercept

# as multi-level model with random intercept
res.lmer <- lmer(time_max ~ treatment + (1|ID), data=fig2a.max)
##              Estimate Std. Error       df   t value     Pr(>|t|)
## (Intercept) 1023.7187   43.37421 59.71689 23.602014 5.703668e-32
## treatmentVa  132.6875   55.01752 30.99995  2.411732 2.198920e-02

Re-analysis – What about the mean and not max response?

Paired t-test of means

# t test
# get mean for each week x ID x treatment combo
fig2a.mean <- fig2a[, .(time_mean=mean(time)), by=.(week, ID, treatment)]
va <- fig2a.mean[treatment=="Va", ]
lb <- fig2a.mean[treatment=="Lb", ]
fig2a_wide <- merge(va, lb, by="ID")
res.t <- t.test(fig2a_wide[, time_mean.x], fig2a_wide[, time_mean.y], paired=TRUE)
##  Paired t-test
## data:  fig2a_wide[, time_mean.x] and fig2a_wide[, time_mean.y]
## t = 1.7521, df = 31, p-value = 0.08964
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -12.38639 163.40722
## sample estimates:
## mean of the differences 
##                75.51042

The p-value of the t-test using max performance was 0.02 and that using mean performance is 0.09. These aren’t too different (do not be distracted by the fact that they fall on opposite sides of 0.05) but neither inspires much confidence that the effect will reproduce.

Hierarchical model

# hierarchical model
res2.lmer <- lmer(time ~ treatment + (1|ID), data=fig2a)
##              Estimate Std. Error        df   t value     Pr(>|t|)
## (Intercept) 848.25000   32.09252  65.82282 26.431394 9.129879e-37
## treatmentVa  75.51042   36.83330 159.00005  2.050059 4.200108e-02

The p-value is 0.04, which differs from that for the t-test because this model is fit to all the data and not just the mean values.


A simulation of Type I, II, M and S error that arises when estimating the treatment effect from differences in the maximum response.

do_sim <- FALSE
  n <- 32
  p <- 3 # number of trials
  mu <- coef(summary(res2.lmer))[1, "Estimate"]
  beta <- coef(summary(res2.lmer))[2, "Estimate"]
  sigma <-[2, "sdcor"]
  sigma_id <-[1, "sdcor"]
  icc <- sigma_id^2/(sigma_id^2 + sigma^2)
  icov <- matrix(icc*sigma_id^2, nrow=3, ncol=3) # covariance
  diag(icov) <- sigma^2
  niter <- 5000
  b.max <- numeric(niter)
  b.mean <- numeric(niter)
  b.lmer <- numeric(niter)
  p.max <- numeric(niter)
  p.mean <- numeric(niter)
  p.lmer <- numeric(niter)
  mean.max <- numeric(niter)
  mean.mean <- numeric(niter)
  sim_table <- data.table(NULL)
  mu_i <- rnorm(n, mean=mu, sd=sigma_id) # true performance in treatment A
  for(b1 in c(0, beta)){
    for(iter in 1:niter){
      # performance_A <- matrix(rnorm(n*p, mean=mu, sd=sqrt(sigma^2 + sigma_id^2)), nrow=n)
      # performance_B <- matrix(rnorm(n*p, mean=mu, sd=sqrt(sigma^2 + sigma_id^2)), nrow=n) + b1
      # performance_A <- matrix(rnorm(n*p, mean=mu_i, sd=sigma), nrow=n)
      # performance_B <- matrix(rnorm(n*p, mean=mu_i, sd=sigma), nrow=n) + b1
      performance_A <- rmvnorm(n, mean=rep(0, p), sigma=icov) + matrix(mu_i, nrow=n, ncol=p)
      performance_B <- rmvnorm(n, mean=rep(0, p), sigma=icov) + matrix(mu_i, nrow=n, ncol=p) + b1
      # t test of max
      A <- apply(performance_A, 1, max)
      B <- apply(performance_B, 1, max)
      mean.max[iter] <- mean(A)
      t_res <- t.test(B, A, paired=TRUE)
      b.max[iter] <- t_res$estimate
      p.max[iter] <- t_res$p.value
      # t test of mean
      A <- apply(performance_A, 1, mean)
      B <- apply(performance_B, 1, mean)
      mean.mean[iter] <- mean(A)
      t_res <- t.test(B, A, paired=TRUE)
      b.mean[iter] <- t_res$estimate
      p.mean[iter] <- t_res$p.value
      # # hierarchical model
      performance <- rbind(data.table(treatment="A", ID=1:n, performance_A),
                           data.table(treatment="B", ID=1:n, performance_B))
      performance <- melt(performance, id.vars=c("ID", "treatment"),
      sim.lmer <- lmer(time ~ treatment + (1|ID), data=performance)
      b.lmer[iter] <- coef(summary(sim.lmer))[2, "Estimate"]
      p.lmer[iter] <- coef(summary(sim.lmer))[2, "Pr(>|t|)"]
    sim_table <- rbind(sim_table, data.table(effect = b1,
                                             b_max = b.max,
                                             b_mean = b.mean,
                                             b_lmer = b.lmer,
                                             p_max = p.max,
                                             p_mean = p.mean,
                                             p_lmer = p.lmer,
                                             mean_max = mean.max,
                                             mean_mean = mean.mean))
  write.table(sim_table, paste(out_path,"june-25-2019.txt",sep="/"), 
              row.names = FALSE, quote=FALSE, sep="\t")
}else{ # read file from simulation
  sim_table <- fread(paste(out_path,"june-25-2019.txt",sep="/"))

The SD among times using max performance is 34.1 and using mean performance is 26.5.

The increased variance of estimating the treatment using the maximum performance decreases power. The hierarchical model, which uses all the data, has the most power (at some cost to Type I?).

niter <- nrow(sim_table)/2
error_table <- data.table(Method=c("max", "mean", "lmm"),
                     "Type 1" = c(sum(sim_table[effect==0, p_max < 0.05])/niter,
                     sum(sim_table[effect==0, p_mean < 0.05])/niter,
                     sum(sim_table[effect==0, p_lmer < 0.05])/niter),
                     "Power" = c(sum(sim_table[effect!=0, p_max < 0.05])/niter,
                     sum(sim_table[effect!=0, p_mean < 0.05])/niter,
                     sum(sim_table[effect!=0, p_lmer < 0.05])/niter))
knitr::kable(error_table, digits=c(NA, 3, 3))
Method Type 1 Power
max 0.050 0.323
mean 0.048 0.489
lmm 0.057 0.530

Type S and M error can be seen by looking at the distribution of effects (more inflated effects, and more sign effects)

sim_table_long <- melt(sim_table, id.vars=c("effect"), 
                       measure.vars=c("b_max", "b_mean", "b_lmer"),
gg <- ggplot(data=sim_table_long, aes(x=method, y=estimate)) +
  geom_boxplot() +
  facet_grid(.~effect, labeller = "label_both") +

Run through a statistical significance filter, the mean effects, of those that are significant, estimated using the max response, the mean response, and a hierarchical model, are:

p_table_long <- melt(sim_table, id.vars=c("effect"), 
                       measure.vars=c("p_max", "p_mean", "p_lmer"),
knitr::kable(sim_table_long[effect!=0 & estimate > 0 & p_table_long[, p] < 0.05, 
               .(effect=mean(estimate)), by=method], digits=c(NA, 1))
method effect
b_max 126.8
b_mean 105.4
b_lmer 103.7

which is inflated for all three (the true effect is 75.5) but most inflated using maximum performance.