Simultaneous Comparison of Adjusted Means in ANCOVA With DMAR

1. The Worked Example From Bryant & Bruvold (1980)

The test_market data ship with DMAR. A company compared six marketing strategies (“panels”) for a brand across four blocks of retail outlets, with the remaining category movement in each outlet as a random covariate.

data(test_market)
descriptives(test_market[c("brand_movement", "category_movement")])
#> $descriptives
#>            variable    type  n n_missing prop_missing      mean median
#> 1    brand_movement numeric 24         0            0  4.220833   4.25
#> 2 category_movement numeric 24         0            0 12.196667  12.66
#>          sd  min   max    q25    q75     skewness   kurtosis
#> 1 0.6999995 2.98  5.61 3.8750  4.575  0.020064949 -0.3790459
#> 2 2.8304196 7.88 16.51 9.5925 14.235 -0.005730893 -1.5356582
#> 
#> $correlations
#> NULL

fit <- lm(brand_movement ~ panel + block + category_movement, data = test_market)

s_ancova <- summary(fit)$sigma           # ANCOVA error SD = sqrt(MS error)
nu       <- fit$df.residual              # error degrees of freedom
c(slope = unname(coef(fit)["category_movement"]),
  error_MS = s_ancova^2, df = nu)
#>       slope    error_MS          df 
#>  0.40788407  0.01325859 14.00000000

xbar <- mean(test_market$category_movement)
adj_means <- vapply(levels(test_market$panel), function(p) {
  nd <- data.frame(panel = factor(p, levels = levels(test_market$panel)),
                   block = factor(1:4, levels = levels(test_market$block)),
                   category_movement = xbar)
  mean(predict(fit, nd))
}, numeric(1))
round(adj_means, 3)     # 3.595 3.619 4.102 4.515 4.618 4.876
#>     1     2     3     4     5     6 
#> 3.595 3.619 4.102 4.515 4.618 4.876

bp <- ci_c_ancova_bp(adj_means = adj_means, s_ancova = s_ancova,
                     n = 4, num_covariates = 1, df = nu)
bp

attr(bp, "critical_value")                      # q_.05;1,6,14 = 4.83
#> [1] 4.829856
unique(round((bp$upper_limit - bp$lower_limit) / 2, 3))   # 0.278
#> [1] 0.278

2. Simultaneous vs. Per-Comparison Intervals

DMAR already provides ci_c_ancova() for a single contrast of adjusted means. It uses a t critical value and includes the per-pair covariate term in the standard error, correct for one pre-planned comparison, but it does not protect the familywise error rate across many comparisons, and it treats the covariate adjustment as fixed. ci_c_ancova_bp() is its familywise, random-covariate counterpart.

# Per-comparison interval for panel 1 vs. panel 4 (needs the covariate means
# and the within-group SS of the covariate).
cov_means <- tapply(test_market$category_movement, test_market$panel, mean)
SSwx <- sum((test_market$category_movement -
             ave(test_market$category_movement, test_market$panel))^2)

per_comparison <- ci_c_ancova(
  adj_means = adj_means, s_ancova = s_ancova,
  c_weights = c(1, 0, 0, -1, 0, 0), n = 4,
  cov_means = cov_means, SSwithin_x = SSwx)
per_comparison

Confidence level: 95%

# The same contrast, but as one member of the simultaneous family:
ci_c_ancova_bp(adj_means = adj_means, s_ancova = s_ancova, n = 4,
               num_covariates = 1, df = nu,
               c_weights = c(1, 0, 0, -1, 0, 0))

Confidence level: 95%

The simultaneous interval is wider; that is the price of protecting all 15 comparisons at once rather than just this one.

3. A One-Way ANCOVA With ancova()

For the common one-way case (a single treatment factor with several levels and one or more covariates), DMAR’s ancova() gives a tidy fit (omnibus F, effect size CIs, adjusted means, and a homogeneity-of-regression check) and feeds directly into ci_c_ancova_bp(). Here is a four-group example.

set.seed(113)
k <- 4; n <- 25
group <- factor(rep(c("control", "low", "medium", "high"), each = n),
                levels = c("control", "low", "medium", "high"))
x <- rnorm(k * n, 50, 10)                          # random covariate
mu <- c(control = 0, low = 1.5, medium = 3, high = 3.2)   # adjusted effects
y <- 0.5 * (x - 50) + mu[as.integer(group)] + rnorm(k * n, 0, 5)
dat <- data.frame(group, x, y)

fit_tidy <- ancova(dat, outcome = "y", treatment = "group", covariates = "x")
fit_tidy

Confidence level: 95%

Pull the adjusted means and the ANCOVA error SD out of the tidy table (or any fitted model) and pass them to ci_c_ancova_bp(). With one covariate and N - k - p error degrees of freedom, the default df is correct, so it need not be supplied.

adj <- fit_tidy$value[grepl("^adjusted_mean", fit_tidy$term)]
names(adj) <- levels(group)

# ANCOVA error SD from the fitted model.
s_yx <- summary(lm(y ~ group + x, data = dat))$sigma

ci_c_ancova_bp(adj_means = adj, s_ancova = s_yx, n = n, num_covariates = 1)

Confidence level: 95%

The intervals that exclude zero identify the groups whose adjusted means differ, with familywise 95% protection and, through the Bryant-Paulson critical value, an accounting for the extra sampling variability the random covariate introduces.

4. Critical Values on Their Own

If you only need the critical value (for a table, a power calculation, or a hand computation), call qbryant_paulson() directly. It is the ANCOVA analogue of qtukey().

# 95% critical value for k = 5 groups, p = 2 covariates, nu = 40 error df:
qbryant_paulson(0.95, num_covariates = 2, num_groups = 5, df = 40)
#> [1] 4.145129

# How much larger than the (incorrect) Tukey value?
qbryant_paulson(0.95, 2, 5, 40) / qtukey(0.95, nmeans = 5, df = 40)
#> [1] 1.026245

References

Bryant, J. L., & Paulson, A. S. (1976). An extension of Tukey’s method of multiple comparisons to experimental designs with random concomitant variables. Biometrika, 63, 631–638.

Bryant, J. L., & Bruvold, N. T. (1980). Multiple comparison procedures in the analysis of covariance. Journal of the American Statistical Association, 75(372), 874–880.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge.