Covariate Adjustment in Randomized Controlled Trials

Aim and statistical setting

Baseline covariates can improve precision in a randomized controlled trial (RCT) when they predict the outcome. Adjustment is most useful when its target, estimator, and variance calculation are defined together. RCTCovAdj implements the continuous-outcome procedures studied in the accompanying paper and makes these choices explicit.

For each participant, let \(O=(X,A,Y)\) denote the observed data, where \(X\) contains covariates measured before treatment, \(A\in\{0,1\}\) is the randomized assignment, and \(Y\) is the observed outcome. The main analysis assumes independent sampling from a target population and Bernoulli assignment with known probability

\[ \Pr(A=1\mid X)=\pi, \qquad 0<\pi<1. \]

Here \(Y(a)\) is the potential outcome under assignment \(a\), assignment \(A\) is independent of \((X,Y(0),Y(1))\), and consistency gives \(Y=AY(1)+(1-A)Y(0)\). Both potential outcomes have finite second moments.

The treatment-specific conditional means are \(m_a(x)=E(Y\mid A=a,X=x)\). The target is the marginal mean difference

\[ \psi=E\{m_1(X)-m_0(X)\}=E\{Y(1)-Y(0)\}. \]

Randomization identifies this effect without a model for \(m_0\) or \(m_1\). Covariate adjustment changes how precisely the effect is estimated; it does not change the marginal target when both treatment-specific means are averaged over the same covariate distribution.

The main theory also assumes that the outcome and every baseline variable used for adjustment are observed for every participant. Missing outcomes require assumptions that identify the intended estimand and an analysis aligned with the outcome-observation process. Covariate-adaptive or fixed-count randomization, other dependent assignment schemes, and finite-population targets require design- and target-specific repeated-sampling arguments. Simply adding covariates to a regression does not supply these arguments.

Where the precision gain comes from

Under known Bernoulli randomization, the efficient influence function for the marginal mean difference is

\[ D^*(O)=m_1(X)-m_0(X)-\psi +\frac{A}{\pi}\{Y-m_1(X)\} -\frac{1-A}{1-\pi}\{Y-m_0(X)\}. \]

The baseline signal that can be removed from the unadjusted estimator is governed by

\[ q^*(x)=(1-\pi)m_1(x)+\pi m_0(x). \]

The conditional mean in the less frequently assigned arm receives greater weight. This opposite-arm weighting matters under unequal allocation. It also shows why strong prediction within each arm does not by itself guarantee a precision gain: the two arm-specific signals can cancel in \(q^*\).

For any fixed square-integrable working function \(q\), the augmented estimating function

\[ \varphi_q(O)=\frac{AY}{\pi}-\frac{(1-A)Y}{1-\pi} -\frac{A-\pi}{\pi(1-\pi)}q(X) \]

has mean \(\psi\). Its first-order variance exceeds the efficiency bound by

\[ \frac{E[\{q(X)-q^*(X)\}^2]}{\pi(1-\pi)}. \]

The ratio difference in means corresponds at first order to the constant \(q_U=E\{q^*(X)\}\), whereas \(q=0\) gives a fixed-denominator contrast. It follows that the gain from the ratio difference in means to the efficient estimator is \(\operatorname{Var}\{q^*(X)\}/[\pi(1-\pi)]\).

Thus approximation quality has a direct variance interpretation. A poor working adjustment can lose precision. With honest cross-fitting, a fixed number of folds whose proportions remain bounded away from zero, and the stated moment conditions, an estimator whose learned adjustment converges in population mean square to the deterministic function \(q^*\) attains the semiparametric bound.

Four analysis strategies

rct_adjust() reports three fixed-dimensional procedures. Supplying the known design probability through allocation keeps the analysis aligned with the randomization law.

Unadjusted difference in means (DIM). The estimate is the observed treatment-arm mean minus the observed control-arm mean. Ratio normalization uses the realized arm sizes, although its first-order analysis depends on the design probability.
Common-slope analysis of covariance (ANCOVA). A pooled regression uses a single covariate slope. Under balanced allocation it has the same first-order slope weighting as interacted linear adjustment. Under unequal allocation and different arm-specific slopes, it can be less precise than the unadjusted estimator.
Interacted standardization. Separate regressions are fit in the two arms, and their fitted contrast is averaged over all observed baseline covariates. This estimates the variance-minimizing adjustment within the chosen linear feature span. Its influence-function standard error includes variation from sampling the covariate distribution used for standardization.
Cross-fitted adjustment. rct_crossfit() trains predictions outside each evaluation fold and averages one held-out estimating contribution per participant. With known allocation, mean-square convergence of the weighted prediction to a deterministic limit suffices for the first-order expansion; no propensity-score model is needed.

The term linear refers to a feature span, rather than an assertion that the true conditional means are linear. Predetermined transformations such as squares or splines can be supplied as columns of the covariate matrix.

A small simulated RCT

The following example has a continuous outcome, balanced assignment, and treatment-specific slopes. The code is intentionally small enough to run while the vignette is built.

set.seed(20260904)
n <- 600L
x <- rnorm(n)
a <- rbinom(n, size = 1L, prob = 0.5)
y <- 0.5 * a + (2 - 3 * a) * x + rnorm(n)

trial <- data.frame(outcome = y, treatment = a, x = x)
head(trial)
#>      outcome treatment          x
#> 1 -3.7188494         0 -1.2556547
#> 2  0.9776823         1  0.9486140
#> 3 -0.1957867         0  0.3078025
#> 4  2.0409626         1 -1.1732217
#> 5  1.6512907         0  1.0691456
#> 6  0.2436397         0  0.2817702

The estimand is the marginal treatment mean difference, which equals \(0.5\) under this data-generating law. The adjusted analyses use the same covariate and the same marginal target.

fit <- rct_adjust(
  outcome = trial$outcome,
  treatment = trial$treatment,
  covariates = trial["x"],
  allocation = 0.5,
  methods = c("unadjusted", "ancova", "interacted")
)
fit
#> Covariate-adjusted marginal mean differences
#>                      method estimate std.error conf.low conf.high
#>                  Unadjusted   0.3703    0.1487   0.0788    0.6618
#>         Common-slope ANCOVA   0.4288    0.1475   0.1397    0.7180
#>  Interacted standardization   0.4615    0.1464   0.1745    0.7484
#>                                 inference   n n.control n.treated allocation
#>              empirical influence function 600       284       316        0.5
#>                  HC0 coefficient sandwich 600       284       316        0.5
#>  joint standardization influence function 600       284       316        0.5

The columns labeled estimate, std.error, and conf.low–conf.high give the point estimate, its standard error, and a normal-reference confidence interval. Sampling variation may make one method look better in a single dataset. Precision comparisons concern repeated-sampling variances, not the distance of one realized estimate from the known simulation truth.

The cross-fitted linear procedure uses held-out predictions. Fold construction is deterministic once the seed is fixed; each fitted regression excludes the observations on which its estimating contributions are evaluated.

fit_cf <- rct_crossfit(
  outcome = trial$outcome,
  treatment = trial$treatment,
  covariates = trial["x"],
  allocation = 0.5,
  learner = "linear",
  folds = 2L,
  seed = 2718L
)
fit_cf
#> Covariate-adjusted marginal mean differences
#>                          method estimate std.error conf.low conf.high
#>  Cross-fitted linear adjustment   0.4583    0.1461   0.1719    0.7446
#>                        inference   n n.control n.treated allocation
#>  cross-fitted influence function 600       284       316        0.5

The built-in quadratic learner adds the square of each model-matrix column; it does not add pairwise interactions or spline bases. Other nonlinear structures can be supplied as prespecified engineered columns or through a custom learner.

nonlinear_trial <- simulate_rct_case("A", n = 600L, seed = 314L)
rct_crossfit(
  outcome = nonlinear_trial$outcome,
  treatment = nonlinear_trial$treatment,
  covariates = nonlinear_trial["x"],
  allocation = 0.5,
  learner = "quadratic",
  folds = 2L,
  seed = 2718L
)
#> Covariate-adjusted marginal mean differences
#>                             method estimate std.error conf.low conf.high
#>  Cross-fitted quadratic adjustment   0.5216    0.0878   0.3495    0.6936
#>                        inference   n n.control n.treated allocation
#>  cross-fitted influence function 600       269       331        0.5

Standard errors after empirical standardization

An interacted regression coefficient and an empirically standardized contrast are different statistical objects. When arm-specific slopes differ, averaging the fitted contrast over the sampled covariates contributes first-order variation for the fixed superpopulation target. The default interacted result from rct_adjust() uses the joint influence-function variance that includes this component.

Setting diagnostic_hc2 = TRUE additionally reports the coefficient-only heteroskedasticity-consistent type 2 (HC2) calculation used as a diagnostic in the paper.

rct_adjust(
  outcome = trial$outcome,
  treatment = trial$treatment,
  covariates = trial["x"],
  allocation = 0.5,
  methods = "interacted",
  diagnostic_hc2 = TRUE
)
#> Covariate-adjusted marginal mean differences
#>                           method estimate std.error conf.low conf.high
#>       Interacted standardization   0.4615    0.1464   0.1745    0.7484
#>  Interacted coefficient-only HC2   0.4615    0.0804   0.3038    0.6191
#>                                 inference   n n.control n.treated allocation
#>  joint standardization influence function 600       284       316        0.5
#>           coefficient-only HC2 diagnostic 600       284       316        0.5

HC2 changes residual weights to account for leverage. It does not restore the variation omitted by treating the empirical standardization distribution as fixed. Consequently, coefficient-only HC2 inference is not the package default for the superpopulation target.

Application 1: postoperative throat-pain scores

The paper applies these methods to 30-minute postoperative throat-pain scores from a licorice-gargle RCT. RCTCovAdj distributes neither the trial records nor results derived from them. The source portal states that permission from the data contributor or corresponding author is required before the records are used in a new publication. Installing the package grants no right to analyze, publish, or redistribute the records or derived output.

After obtaining the required permission and an authorized copy of the reviewed RDA file, a user can verify and analyze it locally. read_licorice_data() checks both the expected file fingerprint and the required structure before returning the analysis data. This code is not run during package checks because the input is not part of the package.

licorice <- read_licorice_data(file.choose())
licorice_fit <- analyze_licorice(licorice)
licorice_fit$estimates
plot_licorice_results(licorice_fit)

The optional medicaldata package supplies the source-format data through a separate distribution. Its availability does not establish permission for a proposed use. A user whose permission covers that copy may prepare it as follows:

licorice <- prepare_licorice_data()
licorice_fit <- analyze_licorice(licorice)

Calling prepare_licorice_data() does not itself grant permission to analyze, publish, or redistribute the participant records or derived output. In either route, analyze_licorice() returns the estimates, outcome and baseline summaries, bounded-score completion analysis, numerical diagnostics, and working-model coefficients needed to reproduce the paper’s application.

Application 2: power and sample size

For a regular, consistently studentized analysis with influence-function variance factor \(V\), the magnitude of the normal-approximation noncentrality parameter at total sample size \(n\) and planning effect \(\Delta\) is \(\sqrt{n}|\Delta|/\sqrt{V}\). Holding the effect, allocation, type I error, test alternative, and target power fixed therefore gives

\[ \frac{n_j}{n_k}=\frac{V_j}{V_k} \]

before integer rounding and attrition inflation.

The package evaluates one- and two-sided power under the normal approximation and numerically inverts it. The following values reproduce the balanced nonlinear planning law (Case A) used in the paper.

variance <- c(
  unadjusted = 12.5,
  interacted_linear = 8.5,
  efficient_quadratic = 4.0
)

rct_sample_size(
  effect = 0.5,
  variance = variance,
  power = 0.80,
  alpha = 0.05,
  alternative = "two.sided"
)
#> [1] 393 267 126

The corresponding rounded totals are 393, 267, and 126. The variance factor of pooled ANCOVA can exceed the unadjusted factor under unequal allocation. The second planning law (Case C) makes this visible.

design <- paper_power_design()
design$sample_sizes
#>   case         analysis      curve variance_factor relative_variance
#> 1    A Unadjusted (DIM) unadjusted           12.50         1.0000000
#> 2    A      ANHECOVA-IF interacted            8.50         0.6800000
#> 3    A     CF-quadratic  efficient            4.00         0.3200000
#> 4    C Unadjusted (DIM) unadjusted           27.50         1.0000000
#> 5    C      ANHECOVA-IF interacted           15.25         0.5545455
#> 6    C    Pooled ANCOVA     pooled           35.50         1.2909091
#>   required_n_80 achieved_power_80 required_n_90 achieved_power_90 power_n_300
#> 1           393         0.8005559           526         0.9003400   0.6877704
#> 2           267         0.8002038           358         0.9005937   0.8438674
#> 3           126         0.8013024           169         0.9014810   0.9911099
#> 4           864         0.8002839          1156         0.9000452   0.3789962
#> 5           479         0.8001797           641         0.9000210   0.6016954
#> 6          1115         0.8001624          1493         0.9001803   0.3065879
#>   power_n_500 power_n_800
#> 1   0.8853791   0.9793266
#> 2   0.9695941   0.9980784
#> 3   0.9998584   0.9999998
#> 4   0.5683196   0.7693904
#> 5   0.8167451   0.9516901
#> 6   0.4667904   0.6604236
plot_power_design(design)

Two-panel normal-approximation power curves under balanced nonlinear Case A and unequal-allocation Case C, comparing unadjusted, interacted, and efficient quadratic or pooled ANCOVA analyses; symbols mark 80 percent power.

Normal-approximation power under planning Cases A and C. Symbols mark the smallest integer total sample sizes attaining 80% power for each analysis.

These calculations are design consequences of specified planning laws. The semiparametric bound is an optimistic planning variance unless the proposed estimator can learn the relevant prognostic signal at the intended sample size. Use credible external information or prespecified sensitivity scenarios for the variance factors, and assess the complete analysis under the actual assignment scheme.