Design-based RCT estimators and central limit theorems for baseline subgroup and related analyses

3.1 ATE estimators

Under the potential outcomes framework and SUTVA, the data generating process for the observed outcome measure, y i , is a result of the random assignment process:

This relation states that we can observe y i = Y i ( 1 ) for those in the treatment group and y i = Y i ( 0 ) for those in the control group, but not both.

Rearranging (2) generates the following nominal full sample regression model:

where T ∼ i = ( T i − p ) is the centered treatment indicator; τ = Y ̅ ( 1 ) − Y ̅ ( 0 ) is the full sample ATE estimand; Y ̅ ( t ) = 1 n ∑ i = 1 n Y i ( t ) is the mean potential outcome for t ∈ { 1,0 } ; α = p Y ̅ ( 1 ) + ( 1 − p ) Y ̅ ( 0 ) is the intercept (expected outcome); and the “error” term, u i , is

We center the treatment indicator in (3) to facilitate the theory without changing the estimator.

In contrast to usual formulations of the regression model, the residual, u i , is random solely due to T i [16,17,20]. This framework allows individual-level treatment effects, τ i = Y i ( 1 ) − Y i ( 0 ) , to vary across the sample, and is nonparametric because it makes no assumptions about the potential outcome distributions. The model does not satisfy key assumptions of the usual regression model: over the randomization distribution ( R ), u i is heteroscedastic, and E R ( u i ) , Cov R ( u i , u i ′ ) , and E R ( T ∼ i u i ) are nonzero if τ i varies across the sample.

The model in (3) also applies to each subgroup due to randomization. Thus, if we combine each subgroup model using the G ik indicators, we obtain the following pooled model:

where τ k = Y ̅ k ( 1 ) − Y ̅ k ( 0 ) is the subgroup ATE estimand; α k = p Y ̅ k ( 1 ) + ( 1 − p ) Y ̅ k ( 0 ) is the subgroup intercept; and ϵ i = ∑ k = 1 K G ik u ik is the error term with u ik = T i ( Y i ( 1 ) − Y ̅ k ( 1 ) ) + ( 1 − T i ) ( Y i ( 0 ) − Y ̅ k ( 0 ) ) . We include model terms for all subgroups and exclude the grand intercept.

Consider the ordinary least squares (OLS) differences-in-mean estimator for τ k from (4) using data on the full sample:

where n k 1 = ∑ i = 1 n T i G ik and n k 0 = ∑ i = 1 n ( 1 − T i ) G ik are subgroup sizes in the treatment and control groups with sample shares, π k 1 = n k 1 / n 1 and π k 0 = n k 0 / n 0 . We see that τ ˆ k is a ratio estimator because n k 1 and n k 0 are random variables (with hypergeometric distributions).

The finite population CLT in Theorem 4 in the study by Li and Ding [23] applies to τ ˆ k conditional on n k 1 and n k 0 , which are ancillary to (independent of) the potential outcomes. There is a long-standing debate on the merits of conditional inference in such settings [30]. In our DB RCT context, we view repeated sampling over the randomization distribution as applying to the full sample, which leads to random subgroup allocations and the need for unconditional inference to capture what “could” occur. In contrast, a conditional analysis measures what “did” occur and parallels a blocked subgroup design with fixed subgroup sizes. Accordingly, we focus on an unconditional CLT for τ ˆ k , but compare variance estimators using both approaches in our theory and simulations. The key difference is that an unconditional analysis leads to nonlinear ratio estimators with random numerators and denominators (subgroup sizes), which complicates the asymptotic analysis.

Our CLT is provided in Section 3.2 for a more general covariate-adjusted estimator from a working model that includes in (4) a 1 xV vector of fixed, baseline covariates other than the subgroup indicators, x i , with parameter vector, β .

where x ∼ i = ( x i − ∑ k = 1 K G ik x ̅ k ) are centered covariates; x ̅ k = 1 n k ∑ i = 1 n G ik x i are subgroup covariate means; and e i is the error term. While the covariates do not enter the true RCT model in (4) and the ATE estimands do not change, they will increase precision to the extent they are correlated with the potential outcomes. We do not need to assume that the true conditional distribution of y i given x i is linear in x i . We define β in Section 3.2.

We focus on the pooled covariate model in (6) because it is commonly used in practice. In Section 3.3, we discuss extensions to models that interact T ∼ i with x ∼ i and G ik .

Using OLS to estimate the working model in (6) yields the following covariate-adjusted estimator for τ k that is produced by standard OLS statistical packages:

where x ̅ k 1 and x ̅ k 0 are subgroup covariate means for treatments and controls, and β ˆ is the OLS estimator for β (see [23] for a parallel result for full sample analyses).

Our DB theory is conditional on randomizations that yield n k 1 > 0 and n k 0 > 0 so that τ ˆ k and its variance can be defined [25]. These restrictions yield subgroup allocation distributions that are truncated, but these effects disappear as n and n k increase (where we assume p = n 1 / n and π k = n k 1 / n 1 have finite limits). To see this, consider the general case where n k ≤ n 1 and n k ≤ n 0 so that either n k 1 or n k 0 can equal 0 , and define a k 1 = E ( n k 1 | 0 < n k 1 < n k ) as the expected value of the associated Truncated hypergeometric ( n , n 1 , n k ) distribution. We can express a k 1 in terms of the non-truncated expectation, n k p , using a k 1 = ( n k pf ) a k 1 n k pf , where ( n k pf ) is the mean of the Truncated binomial ( n k , p ) distribution with f = 1 − p ( n k − 1 ) 1 − ( 1 − p ) n k − p n k (derivation not shown). As n and n k increase, then a k 1 converges to n k p because a k 1 n k pf converges to 1 (as a hypergeometric mean converges to a binomial mean) and f also converges to 1 . A similar argument holds for n k 0 and when n k > n 1 or n k > n 0 (where there is no truncation if both hold). Relatedly, in finite samples, the restrictions are likely to have little effect on our results as they will hold with probability near 1 for typical subgroup analyses. For instance, even for a very small subgroup with n k = 12 , n = 40 , and p = 0.5 , the restrictions will hold with probability 0.9999. Thus, to simplify notation, we omit the conditioning on positive subgroup allocations to each research group and use the unconditional expectations, n k p and n k ( 1 − p ) , in the analysis.

3.2 Main CLT result

To consider the asymptotic properties of τ ˆ k x (which also apply to τ ˆ k without covariates), we consider a hypothetical increasing sequence of finite populations where n → ∞ . Parameters should be subscripted by n , but we omit this notation for simplicity. We assume that n 1 / n → p * as n → ∞ , so the numbers of treatments and controls both increase with n . In addition, we assume that n k / n → π k * for all k , where π k * > 0 and ∑ k = 1 K π k * = 1 . This implies that each subgroup also grows with n , where the number of subgroups, K , is assumed fixed.

Our CLT builds on Schochet et al. [14] who provided CLTs for RCT ratio estimators with general weights for clustered, blocked designs. We adapt these methods to our setting by treating the subgroup indicators, G ik , as “weights” when computing the subgroup sample means.

Before presenting our CLT, we need to define several terms. First, for t ∈ { 1,0 } , let ε ik ( t ) = ( Y i ( t ) − Y ̅ k ( t ) − x ∼ i β ) denote model residuals for subgroup k , where R ik ( t ) = G ik π k ε ik ( t ) are scaled residuals using the normalized weights, w i / w ̅ = G ik / π k , that sum to n . Second, let S R k 2 ( t ) = 1 n − 1 ∑ i = 1 n R ik 2 ( t ) denote the variance of R ik ( t ) , and let S R k 2 ( 1,0 ) = 1 n − 1 ∑ i = 1 n R ik ( 1 ) R ik ( 0 ) denote the treatment-control covariance. Third, we define D ̅ k as the mean treatment-control difference in the R ik ( t ) residuals, with associated variance:

where S 2 ( τ k ) = 1 n − 1 ∑ i = 1 n ( R ik ( 1 ) − R ik ( 0 ) ) 2 is the heterogeneity of treatment effects. Fourth, we define the variance of G ik as S 2 ( G k ) = 1 n − 1 ∑ i = 1 n 1 π k 2 ( G ik − π k ) 2 = n ( n − 1 ) ( 1 − π k ) π k . Fifth, we require the variances of each covariate, S x k , v 2 = 1 n − 1 ∑ i = 1 n G ik π k 2 ( [ x ∼ i ] v ) 2 for v ∈ { 1 , … , V } , and the full variance-covariance matrix for the covariates, S x , k 2 = 1 n ∑ i = 1 n G ik π k x ∼ i ′ x ∼ i . Finally, we need two outcome-covariate variance-covariance matrices: S x , Y , k 2 ( t ) = 1 n ∑ i = 1 n G ik π k x ∼ i ′ Y i ( t ) and S x Y , k 2 ( t ) = 1 n ∑ i = 1 n 1 π k ( G ik x ∼ i ′ Y i ( t ) − θ ̅ k ) 2 , where θ ̅ k = 1 n ∑ i = 1 n G ik x ∼ i ′ Y i ( t ) is the mean covariance.

We now present our CLT theorem, proved in Supplementary Materials S1.

(C3) Letting g k ( t ) = max 1 ≤ i ≤ n { R ik 2 ( t ) } , as n → ∞ ,

(C4) f 1 = n 1 / n and f 0 = n 0 / n have limiting values, p * and ( 1 − p * ) , for 0 < p * < 1 .

(C5) The subgroup shares, n k / n , converge to π k * for 0 < π k * < 1 and ∑ k = 1 K π k * = 1 .

(C6) As n → ∞ ,

(C7) Letting h v ( t ) = max 1 ≤ i ≤ n G ik π k ( [ x ∼ i ] v ) 2 for all v ∈ { 1 , … , V } , as n → ∞ ,

(C8) S R k 2 ( t ) , S R k 2 ( 1,0 ) , S x k , v 2 , S x , k 2 , S x , Y , k 2 ( t ) , and S x Y , k 2 ( t ) have finite limiting values.

Then, as n → ∞ , τ ˆ k x is a consistent estimator for τ k , and

where Var ( D ̅ k ) is defined as in (8).

Remark 1. The Var ( D ̅ k ) expression in (8) is difficult to interpret because S R k 2 ( t ) and S 2 ( τ k ) are scaled by π k 2 ( n − 1 ) to facilitate the theory. To address this, we apply the following relations in (8): n 1 π k 2 ( n − 1 ) = n k p ( n k − π k ) and n 0 π k 2 ( n − 1 ) = n k ( 1 − p ) ( n k − π k ) , which yields,

where Ω R k 2 ( t ) = 1 n k − 1 ∑ i = 1 n G ik ε ik 2 ( t ) ; Ω 2 ( τ k ) = 1 n k − 1 ∑ i = 1 n G ik ( ε ik ( 1 ) − ε ik ( 0 ) ) 2 ; and ϕ k = ( n k − 1 ) / ( n k − π k ) ≤ 1 is a correction term that reflects the single treatment indicator “shared” by each subgroup (and can be ignored as it converges to 1). The Ω R k 2 ( t ) and Ω 2 ( τ k ) terms are population variances for those in subgroup k , and n k p and n k ( 1 − p ) are expected subgroup sizes in the two research groups. This variance expression is more intuitive as it parallels the full sample asymptotic results in Li and Ding [23], the key difference being that (9) is based on expected subgroup sizes rather than actual ones. Note that for ϕ k = 1 , (9) is the same as for an RCT that stratifies on subgroup k to select fixed subgroup sample sizes, n k p and n k ( 1 − p ) .

Remark 2. The first two terms in (9) pertain to separate variances for the two research groups because we allow for heterogeneous treatment effects. The third term pertains to the treatment-control covariance, Ω R k 2 ( 1,0 ) , expressed in terms of the heterogeneity of treatment effects, Ω 2 ( τ k ) , which cannot be identified from the data but can be bounded [31].

Remark 3. (C3) and (C7) are Lindeberg-type conditions from Li and Ding [23] that control the tails of the potential outcome and covariate distributions. (C6) yields a weak law of large numbers for the observed subgroup shares so that π k t / π k → p 1 (using Theorem B in Scott and Wu [24]). This condition is used to account for the randomness in n k t , because, for instance, it allows us to express the sample mean for the treatment group as np n k 1 ∑ i = 1 n G ik T i Y i ( 1 ) np , where the bracketed term converges to 1 by (C6) and the denominator in the summation term is fixed, so we can apply the CLT results in Li and Ding [23] and Slutsky’s theorem. While (C6) is implied by (C4) and (C5), it facilitates the addition of other weights (Section 3.3). (C8) specifies limiting values of the variances and variance-covariance matrices.

Remark 4. Theorem 1 is proved in two stages by expressing the ATE estimator in (7) as,

where τ ˆ k x β = ( y ̅ k 1 − y ̅ k 0 ) − ( x ̅ k 1 − x ̅ k 0 ) β and β = ∑ k = 1 K π k S x , k 2 − 1 [ ∑ k = 1 K p π k S x , Y , k 2 ( 1 ) + ∑ k = 1 K ( 1 − p ) π k S x , Y , k 2 ( 0 ) ] is assumed known. This β parameter is the (hypothetical) population OLS coefficient that would result from a regression of [ p Y i ( 1 ) + ( 1 − p ) Y i ( 0 ) ] on the covariates. In the first stage, we obtain a CLT for τ ˆ k x β . In the second stage, we prove that τ ˆ k x has the same asymptotic distribution as τ ˆ k x β by showing that ( x ̅ k 1 − x ̅ k 0 ) ( β ˆ − β ) = o p ( n − 1 / 2 ) , which holds under our conditions because x ̅ k 1 and x ̅ k 0 are both asymptotically normal and β ˆ − β → p 0 .

Remark 5. Under (C1)–(C6), Theorem 1 also applies to τ ˆ k for the model without covariates by setting β = 0 in (6) and in the residuals, ε ik ( t ) and R ik ( t ) , which enter the variances in (8) and (9). In Supplement S2.1, we prove that τ ˆ k is unbiased using the approach reported by Miratrix et al. [25] and Schochet [18], where we condition on n k 1 > 0 and n k 0 > 0 and then average over possible subgroup allocations ( A ) to the two research groups to show that E R ( τ ˆ k ) = E A E R ( τ ˆ k | n k 1 , n k 0 ) = τ k . Similarly, using the law of total variance, Supplement S2.1 shows that

where Ω R k 2 and Ω 2 are defined in (9) with β = 0 . Note that lim n → ∞ E A 1 n k 1 = 1 E A ( n k 1 ) = 1 n k p and lim n → ∞ E A 1 n k 0 = 1 n k ( 1 − p ) which aligns with (9) in large samples. In finite samples, however, the variance in (11) is at least as large as in (9) because E A 1 n k 1 ≥ 1 n k p and E A 1 n k 0 ≥ 1 n k ( 1 − p ) by Jensen’s inequality. Our simulations include both sets of sample sizes (Section 5).

The following corollary to Theorem 1, proved in Supplement S1, provides the joint asymptotic distribution of the subgroup estimators, ( τ ˆ 1 x , … , τ ˆ K x ) .

This corollary is important for real-world applications because it supports the use of standard F-tests (or chi-square tests) to test the null hypothesis of equal subgroup effects.

3.3 Extensions to related estimators

This section outlines extensions of our CLT result to post-stratification estimators, models that interact T ∼ i with x ∼ i and G ik , BTs, and the use of nonresponse weights to adjust for missing outcome data.

Post-stratification estimators. Miratrix et al. [25] considered variance estimation for a DB post-stratification ATE estimator that obtains overall effects for the model without covariates by averaging τ ˆ k across subgroups. With model covariates, we can express this estimator as, τ ˆ PS x = 1 n ∑ k = 1 K n k τ ˆ k x . Corollary 1 from above can then be applied to yield a new CLT for this estimator: τ ˆ PS x − ( Y ¯ ( 1 ) − Y ¯ ( 0 ) ) Var ( D ̅ PS ) → d N ( 0,1 ) , where Var ( D ̅ PS ) = 1 n 2 ∑ k = 1 K n k 2 Var ( D ̅ k ) .

Interacted models. Theorem 1 can be extended to a model that replaces x ∼ i β in (6) with the interaction terms, ∑ k = 1 K G ik ( 1 − T i ) x ∼ i β k 0 and ∑ k = 1 K G ik T i x ∼ i β k 1 , which allows covariate effects to differ by subgroup and treatment status. The ATE estimator for this model is, τ ˆ k x GT = [ y ̅ k 1 − ( x ̅ k 1 − x ̅ k ) β ˆ k 1 ] − [ y ̅ k 0 − ( x ̅ k 0 − x ̅ k ) β ˆ k 0 ] . Theorem 1 can then be applied by redefining the residuals as, R ik ( t ) = G ik π k ( Y i ( t ) − Y ̅ k ( t ) − ( x i − x ̅ k ) β k t ) , where β k t = ( S x , k 2 ) − 1 S x , Y , k 2 ( t ) . The proof (not shown) follows using the same arguments as in Remark 4 by replacing (10) with τ ˆ k x GT = τ ˆ k x GT β − ∑ t = 0 1 ( x ̅ k t − x ̅ k ) ( β ˆ k t − β k t ) , and noting that under our regularity conditions, β ˆ k t − β k t → p 0 and ( x ̅ k t − x ̅ k ) ( β ˆ k t − β k t ) = o p ( n − 1 / 2 ) for t ∈ { 1,0 } , so that ∑ t = 0 1 ( x ̅ k t − x ̅ k ) ( β ˆ k t − β k t ) = o p ( n − 1 / 2 ) .

BTs. Our results also extend to BTs where each sample member is independently randomized to the treatment group with probability p , leading to random treatment-control sizes. This design pertains, e.g., to an RCT with rolling study intake. First, we can show that our result in Remark 5 on the unbiasedness of τ ˆ k for the model without covariates also applies to BTs. To see this, consider the full sample analysis. Then, Bernoulli sampling has the same properties as a two-stage design that first randomly selects n 1 from a truncated binomial distribution, and then selects a simple random sample of size n 1 to the treatment group [32]. Thus, this setting parallels the one in Remark 5 that calculates sample moments by first conditioning on subgroup sizes. The only difference is that E A is now taken over a truncated binomial rather than truncated hypergeometric distribution. The same argument applies to the subgroup analysis as for the full sample analysis.

Second, we can adapt the CLT in Theorem 1 to BTs by using expected rather than actual sample sizes in the theorem conditions (i.e., by replacing n 1 and n 0 with np and n ( 1 − p ) ). To see this, consider again the full sample analysis, omitting the k = 1 subscript for simplicity. Let p 1 be the observed treatment share, and express the observed mean outcomes as, y ̅ 1 = y ̅ BT 1 g 1 and y ̅ 0 = y ̅ BT 0 g 0 , where y ̅ BT 1 = 1 np ∑ i = 1 n T i y i and y ̅ BT 0 = 1 n ( 1 − p ) ∑ i = 1 n ( 1 − T i ) y i are divided by expected sample sizes rather than actual ones, and g 1 = p p 1 and g 0 = ( 1 − p ) ( 1 − p 1 ) . Note that g t → p 1 for t ∈ { 1,0 } , so y ̅ t and y ̅ BT t have the same asymptotic distributions by Slutsky’s theorem. Then, under our conditions, Theorem 4 in the study by Li and Ding [23] provides a CLT for ( y ̅ BT 1 − y ̅ BT 0 ) , and the same proof as in Supplement S1.2 for Theorem 1 extends this CLT to ( y ̅ 1 − y ̅ 0 ) . A similar approach yields a CLT for the covariate-adjusted ATE estimator using the variance in (9) by applying x ̅ t = x ̅ BT t g t and noting that the asymptotic properties of β ˆ do not change from Theorem 1. The same argument applies to the subgroup analysis.

Data nonresponse weights. Theorem 1 also applies, with additional assumptions, to a “subgroup” analysis that adjusts for missing outcome data using respondents only with nonresponse weights, w i r . To show this, let R i ( T i ) denote an indicator of potential data response in the treatment or control condition, where R i ( T i ) = 1 for a respondent and 0 for a nonrespondent. Further, let r i = R i ( T i ) denote the observed response indicator. If baseline covariates are available for the full sample, a common approach is to set w i rt = 1 / e t ( x i ) as the weight for a respondent in research group t ∈ { 1,0 } , where e t ( x i ) = Pr ( r i = 1 | x i , T i = t ) is the propensity score [33].

We invoke two missing data assumptions for each subgroup: (i) data are missing at random for each research group conditional on covariates [34]: Y i ( 1 ) , Y i ( 0 ) ⫫ r i | T i = t , G ik = 1 , x i ; and (ii) data response and treatment status are independent conditional on the covariates: r i ⫫ T i | G ik = 1 , x i . The first ignorability (selection-on-observables) condition – which is commonly invoked for observational studies using inverse probability weighting methods [35,36] – ensures that weighted ATE estimators using the respondent sample will consistently estimate τ k . The second condition implies that the response, r i = R i ( T i ) = R i , will be identical in the treatment and control conditions, so that e t ( x i ) = e ( x i ) and w i rt = w i r are independent of t . Stated differently, this condition implies that respondents and their weights are randomly allocated to the two research groups. Thus, the summed weights, ∑ i : T i = t w i r , which enter the denominators of the weighted differences-in-means estimators become random, which parallels the subgroup analysis from above with random subgroup sizes in the two research groups.

Accordingly, we can apply Theorem 1, assuming known weights, where the respondent sample and weighted least squares are used to obtain the ATE estimator, τ ˆ k r x , using (6). We assume e ( x i ) is known and converges to e * ( x i ) as n → ∞ , where 0 < e * ( x i ) ≤ 1 for all x i in its finite population support so that (C6) holds. For the proof, we replace the weights, w i = G ik , in the theorem with w i = G ik r i w i r to define the variables and regularity conditions. The resulting variance for the CLT has the same form as (9) but is based on expected subgroup respondent sizes (Supplement S2.2). Developing a finite population CLT that allows for estimated nonresponse weights rather than known weights is a topic for future research.

3.4 Variance estimation

To obtain consistent variance estimators for (9) (and model variants), we can either use expected subgroup sizes, n k p and n k ( 1 − p ) , or actual ones, n k 1 and n k 0 , as for the conditional analysis. Our simulations find very similar results using either approach (Section 5). This occurs because the difference between the hypergeometric random variable, π k t , and its expected value, π k , decreases exponentially with n t [37,38]. For instance, Figure 1a shows that for modest n t , there is a high probability that | π k t − π k | / π k ≤ c for small c (defined as 10 or 20%).

Figure 1

(a) Probabilities for the differences, ( π k t – π k ) , relative to π k , and (b) SE ratios using actual to expected subgroup sizes, as a function of δ 1 / 5 n k for δ 1 = ( n k 1 – n k p ) . Note: See text for definitions and formulas. (b) assumes n = 100 , φ = 1.1 , and ϑ = 0.05 . SE = Standard error.

Further, Figure 1b shows that for small c , the ratios of SEs using actual to expected subgroup sizes in (9) are close to 1, leading to similar confidence interval coverage. For example, for n = 100 , p = 0.5 , and π k = 0.5 , the ratios range only from 1 to 1.026 as c ranges from 0 to 0.2 (assuming Ω R k 2 ( 1 ) = φ Ω R k 2 ( 0 ) and Ω 2 ( τ k ) = ϑ Ω R k 2 ( 0 ) with plausible values, φ = 1.1 and ϑ = 0.05 ). In expectation, the SE ratios are greater than 1 for all values of φ and ϑ (Remark 5 above), but the differences are small for typical subgroup sizes used in practice.