Ignoring the Covariate

    t.test(y ~ arm, data = dat)

## 
##  Welch Two Sample t-test
## 
## data:  y by arm
## t = 0.29696, df = 7.9985, p-value = 0.7741
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -68.57353  88.84470
## sample estimates:
## mean in group 1 mean in group 2 
##        171.4535        161.3179

No significance, because that gargantuan covariate effect which adds to the observed variability of the outcom is overshadowing the treatment difference. The linear model ignoring the covariate confirms this:

    mod <- lm(y ~ arm, data = dat)
    emmeans(mod, "arm") %>%
      contrast("pairwise")

##  contrast    estimate   SE df t.ratio p.value
##  arm1 - arm2     10.1 34.1  8   0.297  0.7741

Including the Covariate

    mod_cov <- lm(y ~ arm + cov, data = dat)
    emmeans(mod_cov, "arm") %>%
      contrast("pairwise")

##  contrast    estimate   SE df t.ratio p.value
##  arm1 - arm2    -9.54 2.19  7  -4.356  0.0033
## 
## Results are averaged over the levels of: cov

Magically, your effect is back in the limelight. By adjusting for the covariate, you reduce the unexplained variation in your model and more precisely isolate the treatment effect. Applause.

Simulation

Below is a very nice illustration of the phenomenon.

    set.seed(999)
    n_sim <- 200
    p_cov_combined <- data.frame()

    # pb <- txtProgressBar(min = 0, max = n_sim, style = 3)

    for (i in 1:n_sim) {
      dat <- sim_cov_study(n_g = 10, eff = 10)
      
      # Model w/o the covariate
      mod_nocov <- lm(y ~ arm, data = dat)
      res_nocov <- emmeans(mod_nocov, "arm") %>%
        contrast("pairwise") %>%
        as.data.frame()
      res_nocov$model <- "No covariate"
      res_nocov$rep <- i

      # Model w/ the covariate
      mod_withcov <- lm(y ~ arm + cov, data = dat)
      res_withcov <- emmeans(mod_withcov, "arm") %>%
        contrast("pairwise") %>%
        as.data.frame()
      res_withcov$model <- "With covariate"
      res_withcov$rep <- i

      p_cov_combined <- rbind(p_cov_combined, res_nocov, res_withcov)
      # setTxtProgressBar(pb, i)
    }
    # close(pb)

    p_cov_combined <- p_cov_combined %>%
      mutate(
        ci_low = estimate - 1.96 * SE,
        ci_upp = estimate + 1.96 * SE
      )

    p_cov_combined <- p_cov_combined %>%
      group_by(model) %>%
      arrange(estimate, .by_group = TRUE) %>%
      mutate(id = row_number()) %>%
      ungroup()

      ggplot(p_cov_combined, aes(x = id, y = estimate, ymin = ci_low, ymax = ci_upp)) +
      geom_pointrange() +
      facet_wrap(~ model, scales = "free_x") +
      geom_hline(yintercept = -10, linetype = "dashed", color = "salmon4") +
      theme_minimal() +
      labs(
        x = "Simulation Index",
        y = "Estimated Contrast (arm2 - arm1)",
        title = "Covariate? No vs Yes",
        caption = "Dashed line indicates the true effect of 10 (shown as -10)."
      )

Notice how “No covariate” bounces around widely (with “clusters” denoting the per arm distribution of the covariate), whereas “With covariate” aligns with the known effect with a confidence resembling that of the scenario described before.

Comparing Models

1. Naive: just look at change ~ arm.

2. Baseline correction: look at y_after ~ arm + y_before, which can reduce variance when y_before is highly predictive.

The second setup results in smaller standard errors.

    mod_simp <- lm(change ~ arm, data = dat)
    mod_simp %>% effects::predictorEffects(partial.residuals = TRUE) %>% plot(main="No Baseline in Model")

    mod <- lm(y_after ~ arm + y_before, data = dat)
    mod %>% effects::predictorEffects(partial.residuals = TRUE) %>% plot(main="With Baseline in Model")

    # Quick look at the estimated effect
    emmeans(mod_simp, ~ arm) %>%
      contrast("pairwise")

##  contrast    estimate  SE df t.ratio p.value
##  arm1 - arm2    -9.12 1.5 18  -6.090  <.0001

    emmeans(mod, ~ arm) %>%
      contrast("pairwise")

##  contrast    estimate   SE df t.ratio p.value
##  arm1 - arm2    -8.51 1.34 17  -6.363  <.0001

Simulation

    set.seed(999)
    n_sim <- 200
    res_baseline <- data.frame()

    # pb <- txtProgressBar(min = 0, max = n_sim, style = 3)

    for (i in 1:n_sim) {
      dat <- sim_baseline_study(n_g = 5, eff = 10, eff_corr = 0.3, sd = 3)
      
      # 1) no baseline correction
      mod_nobc <- lm(change ~ arm, data = dat)
      tmp_nobc <- emmeans(mod_nobc, ~ arm) %>%
        contrast("pairwise") %>%
        as.data.frame()
      tmp_nobc$model <- "Measuring raw change"
      tmp_nobc$rep   <- i

      # 2) with baseline correction
      mod_bc <- lm(y_after ~ arm + y_before, data = dat)
      tmp_bc <- emmeans(mod_bc, ~ arm) %>%
        contrast("pairwise") %>%
        as.data.frame()
      tmp_bc$model <- "With baseline correction"
      tmp_bc$rep   <- i

      res_baseline <- rbind(res_baseline, tmp_nobc, tmp_bc)
      # setTxtProgressBar(pb, i)
    }
    # close(pb)

    res_baseline <- res_baseline %>%
      mutate(ci_low = estimate - 1.96 * SE,
             ci_upp = estimate + 1.96 * SE)

    res_baseline <- res_baseline %>%
      group_by(model) %>%
      arrange(estimate, .by_group = TRUE) %>%
      mutate(id = row_number()) %>%
      ungroup()

      ggplot(res_baseline, aes(x = id, y = estimate, ymin = ci_low, ymax = ci_upp)) +
      geom_pointrange() +
      facet_wrap(~ model, scales = "free_x") +
      geom_hline(yintercept = -10, linetype = "dashed", color = "salmon4") +
      theme_minimal() +
      labs(
        x = "Simulation Index (sorted by estimate)",
        y = "Estimated Contrast (arm2 - arm1)",
        caption = "Dashed line shows the true effect of -10"
      )

Yup, baseline correction results in more consistent estimations overall. Using baseline as a covariate typically shrinks standard errors (too).