What is already known?

• • Externally controlled single-arm trials (SATs) are important for assessing treatment efficacy when randomized controlled trials are not feasible.

• • Unanchored indirect treatment comparisons (ITCs), often using the entropy balancing or matching-adjusted indirect comparison (MAIC) approach, can provide treatment effect estimates for disconnected networks in health technology assessment.

What is new?

• • A unified causal inference framework for both externally controlled SATs and unanchored ITCs is presented.

• • A doubly robust (DR) estimator that augments the MAIC approach by combining outcome model predictions with entropy balancing weights is proposed.

• • The proposed estimator provides increased protection against model misspecification, as demonstrated in a simulation study with binary outcomes.

• • The proposed estimator is extended to the setting with unavailable individual participant data for the external control and illustrated through an applied example.

Potential impact for RSM readers

• • DR augmented weighting estimators, particularly using entropy balancing weights, offer increased robustness to model misspecification than conventional non-augmented approaches for confounding adjustment.

• • Entropy balancing-based methods have desirable properties compared to standard “modeling” approaches to weighting but should be augmented to improve protection against bias and guarantee double robustness.

3.1 Data and assumptions

As per Section 2, let $T=t$ denote a time-fixed binary treatment, with $t\in \{0, 1\}$ , such that $T=1$ represents the active intervention and $T=0$ the control. Let $S=s$ denote the data source, with ${s \in \{0, 1\}}$ , such that $S=1$ represents the SAT and $S=0$ the external data source. In addition, let $\mathbf {X}$ denote vector-valued pre-treatment baseline covariates, for example, clinical or demographic characteristics, measured across the SAT and the external data source. Let Y denote the clinical outcome of interest. We assume that only distributional differences in $\mathbf {X}$ are preventing exchangeability between the SAT and external subjects, and that covariates and outcomes are defined and measured similarly across data sources.

The observed IPD consist of $(S_i, \mathbf {X}_i, T_i, Y_i)$ , $i=1, \dots , n_1+n_0$ , realizations of $(S, \mathbf {X}, T, Y)$ denoting the data source, baseline covariates, treatment assignment, and observed outcome for subject i. Here, the SAT and external data source have been stacked, with $n_1$ and $n_0$ as the sample sizes of the SAT and the external data source, respectively. It is assumed that all individuals in the SAT are under $T=1$ and all individuals in the external data source are under $T=0$ , such that the control group is fully external. To be clear, we have $S_i=T_i$ for all $i=1, \dots , n$ , where $n = n_1+n_0$ ; and we also have that $S_i=1$ and $T_i=1$ for all $i=1,\dots ,n_1$ , and $S_i=0$ and $T_i=0$ for all $i=n_1+1, \dots , n$ . We shall assume that there is no missingness or measurement error.

The observed outcome for subject i is $Y_i = Y_i^1T_i + Y_i^0(1-T_i)$ , where $Y_i^t$ is the potential outcome had subject i been assigned treatment $t \in \{0, 1\}$ , with $Y_i=Y_i^1$ if $i=1,\dots , n_1$ and $Y_i=Y_i^0$ if $i=n_1+1, \dots , n$ . Namely, the observed outcome for an individual in the SAT equals their potential outcome under the active intervention, and the observed outcome for an individual in the external data source equals their potential outcome under the control. Implicit in the notation is the stable unit treatment value assumption (SUTVA): that there is no interference between subjects and there is treatment version irrelevance, that is, one well-defined version of the active intervention and the control across all subjects and data sources.Reference Zhou, Zhu, Drake and Pang ⁵⁴ Also implicit is that there is no direct effect of trial participation.Reference Zhou, Zhu, Drake and Pang ⁵⁴ Namely, that trial participation—in the SAT or a historical trial, for that matter—does not affect the outcome except through treatment assignment itself, that is, there are no Hawthorne effects.Reference Braunholtz, Edwards and Lilford ^{55 ,} Reference Dahabreh, Robins, Haneuse and Hernán ⁵⁶

To estimate the ATC, we must construct estimators for $\mu ^1_0 = \text {E}(Y^1 \mid S=0)$ and $\mu ^0_0= \text {E}(Y^0 \mid S=0)$ . Outcomes for the subjects from the external data source have been generated under the control and we assume that there is no informative missingness or measurement error. Hence, unbiased estimation of $\mu ^0_0$ should be trivial using the sample mean, such that ${\widehat{\mu} }_{0}^{0}= \frac {1}{n_0}\sum _{i=n_1+1}^nY_i$ . Conversely, while the active intervention has been investigated in the SAT, its outcomes in the external control (sub) population are unobserved. Our challenge is therefore to produce a reliable estimate $\widehat {\mu }^1_0$ of the mean absolute outcome $\mu ^1_0$ under the active intervention in the external control (sub) population, based on the observed data.

Two causal identification conditions, together known as strong ignorability, are required to construct a valid estimator of $\mu _0^1$ . These ensure that the SAT and external control outcomes are comparable given adjustment for baseline covariates. The first assumption is conditional data source ignorability; formally, $Y_i^1 \perp \!\!\!\!\perp S_i \mid \mathbf {X}_i$ for all $i=1,\dots , n$ . Namely, conditional on baseline covariates, the potential outcome under the active intervention is independent of the data source. This is akin to the conditional constancy or exchangeability of absolute outcomes invoked for unanchored ITCs, used to transport mean absolute outcomes under $T=1$ from $S=1$ to $S=0$ .Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{25 ,} Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ²⁸ Conditional ignorability is a strong assumption, resting on the SAT and the external control capturing all variables that are prognostic of outcome under the active intervention.

The second assumption is positivity or overlap. That is, the support of the baseline covariates in the external control is contained within that of the SAT. Mathematically, the probability of SAT participation, conditional on the covariates necessary to ensure ignorability, should be bounded away from zero and one: $0 < \text {Pr}(S=1 \mid \mathbf {X}=\mathbf {x}) < 1$ for all $\mathbf {x}$ with positive density in the external control, that is, for all $\mathbf {x}$ such that $f(\mathbf {x} \mid S=0)>0$ . Hence, it is possible to have SAT subjects in all regions of the covariate space in $S=0$ .Reference Zhou, Zhu, Drake and Pang ^{54 ,} Reference Dahabreh, Robertson, Tchetgen, Stuart and Hernán ^{57 ,} Reference Buchanan, Hudgens and Cole ⁵⁸ Positivity violations can be deterministic or random. The former arise structurally, due to non-overlapping SAT and external control eligibility criteria. The latter arise empirically due to chance, particularly with small sample sizes.Reference Westreich and Cole ⁵⁹ To enforce positivity, analysts may subset the SAT based on the selection criteria of the external control.Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{25 ,} Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ²⁸ However, this further reduces the sample size of the SAT. Positivity is typically assessed by comparing the empirical distributions of the covariates in the SAT and the external control.Reference Glimm and Yau ⁶⁰ While outcome modeling-based approaches such as G-computation can overcome failures of positivity, they do so by potentially problematic and difficult-to-diagnose model-based extrapolation. Even minor model misspecification over the observed covariate space in the SAT may lead to poor extrapolation in unobserved regions of the covariate space.Reference Vo ⁵⁰

Analogously, targeting the ATT would require constructing estimators for $\mu ^1_1 = \text {E}(Y^1 \mid S=1)$ and $\mu ^0_1= \text {E}(Y^0 \mid S=1)$ . Here, the challenge is the estimation of $\mu ^0_1$ because outcomes under the control have not been generated in the SAT. The conditional ignorability assumption would formally be $Y_i^0 \perp \!\!\!\!\perp S_i \mid \mathbf {X}_i$ for all $i=1,\dots , n$ , and would rest on the SAT and the external control measuring all variables that are prognostic of outcome under the control. The positivity assumption would be $0 < \text {Pr}(S=0 \mid \mathbf {X}=\mathbf {x}) < 1$ for all $\mathbf {x}$ with positive density in the SAT, $f(\mathbf {x} \mid S=1)>0$ , such that the support of the baseline covariates in the SAT is contained within that of the external control and it is possible to have external control subjects in all regions of the SAT covariate distribution. In this setting, analysts may apply the SAT selection criteria to the external control to guarantee that there is sufficient overlap. While the methods in the next sections target the ATC, we present any required modifications to target the ATT in the Supplementary Material.

3.2 Inverse odds weighting

We first present a covariate adjustment method that models the data source assignment mechanism, conditional on baseline covariates, to estimate weights.Reference Ren, Cislo, Cappelleri, Hlavacek and DiBonaventura ^{34 ,} Reference Loiseau, Trichelair and He ^{35 ,} Reference Josey, Berkowitz, Ghosh and Raghavan ^{41 ,} Reference Dahabreh, Robertson, Steingrimsson, Stuart and Hernan ^{61 ,} Reference Colnet, Mayer and Chen ⁶² Where the target estimand is the ATC, SAT subjects are weighted by their inverse conditional odds of SAT participation—their conditional odds of external control participation—to transport the SAT outcomes to the external control (sub) population. Such “inverse odds” weights (IOW) are defined as

(1) $$ \begin{align} w_i = \frac{(1-e_i)S_i} {e_i} + (1-S_i), \end{align} $$

for subject $i=1,\dots , n$ , where the propensity score $e_i=e(\mathbf {X}_i)=\text {Pr}(S_i=1 \mid \mathbf {X}_i)$ denotes the conditional probability of SAT participation given covariates $\mathbf {X}_i$ for subject i. In Equation (1), note that the SAT subjects $(S_i=1)$ are weighted as $w_i=(1-e_i)/e_i$ , whereas the external control subjects $(S_i=0)$ are unweighted, that is, assigned a weight of $w_i=1$ .

In practice, the true propensity scores are unknown. Almost invariably, there are multiple baseline covariates and at least one of these is continuous, such that a data source assignment model is required to estimate the propensity scores. The model is often a logistic regression:

(2) $$ \begin{align} \operatorname{\mathrm{logit}}(e_i) = \alpha_0 + \mathbf{c}(\mathbf{X}_i)^\top \boldsymbol{\alpha}, \end{align} $$

where $\operatorname {\mathrm {logit}}(e_i)=\ln \left ((e_i)/(1-e_i) \right )$ , $\alpha _0=\ln \left (\text {Pr}(S_i=1|\mathbf {c}(\mathbf {X}_i)=0)/(\text {Pr}(S_i=0|\mathbf {c}(\mathbf {X}_i)=0) \right )$ is an intercept term, $\boldsymbol {\alpha }$ is a vector of regression parameters, and $\mathbf {c}(\mathbf {X}_i)= [c_1(\mathbf {X}_{i}), c_2(\mathbf {X}_{i}), \dots , c_k(\mathbf {X}_i)]^\top $ is a vector of covariate “balance functions” for subject $i=1,\dots n$ . This is the set of functions containing the distributional features to be balanced between the SAT and the external control,Reference Josey, Berkowitz, Ghosh and Raghavan ^{41 ,} Reference Josey, Yang, Ghosh and Raghavan ⁶³ potentially including sensible transformations of the covariates, for example, polynomials and interaction terms.

The logistic regression is typically fitted to the concatenated IPD using maximum-likelihood estimation, with the regression coefficient point estimates denoted by $\hat {\alpha_0 }$ and $\boldsymbol{ \hat {\alpha }}$ , and model-based propensity scores for subject $i=1, \dots , n_1$ , predicted by $\hat {e_i}=\operatorname {\mathrm {logit}}^{-1}\left ( \hat {\alpha_0 } + \mathbf {c}(\mathbf {X}_i)^\top \boldsymbol {\hat {\alpha }} \right )=\operatorname {\mathrm {expit}} \left ( \hat {\alpha_0 } + \mathbf {c}(\mathbf {X}_i)^\top \boldsymbol {\hat {\alpha }} \right )$ , where $\operatorname {\mathrm {expit}}(\cdot )=\exp (\cdot )/\left (1 + \exp (\cdot ) \right )$ . Weight estimates $\hat {w_i}$ for $i=1,\dots , n_1$ are derived by plugging the corresponding propensity score predictions into Equation (1). With correct specification of the model in Equation (2), such that the log-odds of SAT participation are linear across the balance functions of the covariates, $\hat {e_i}$ and $\hat {w_i}$ consistently estimate the true conditional probability and inverse odds of SAT participation, respectively.

The ATC is estimated by contrasting the weighted average of observed outcomes under the active intervention with the unweighted average of observed outcomes for the external control. As per Section 2, mean absolute outcomes are converted to the additive scale imposed by link function $g(\cdot )$ prior to taking the difference between treatments on such scale, leading to the IOW estimator for the ATC:

(3) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \frac{1}{n_0}\sum_{i=1}^{n_1} \hat{w_i} Y_i \right ) }_{\widehat{\mu}_0^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) }_{\widehat{\mu}_0^0}, \end{align} $$

The mean absolute outcome estimate for the active intervention can be bounded within its feasible range, for example, between 0 and 1 for probabilities, by normalizing or “stabilizing” the weights so that they sum to one.Reference Dahabreh, Robertson, Tchetgen, Stuart and Hernán ^{57 ,} Reference Dahabreh, Robertson, Steingrimsson, Stuart and Hernan ⁶¹ This results in the alternative normalized IOW estimator for the ATC:

(4) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \frac{ \sum_{i=1}^{n_1} \hat{w_i} Y_i}{\sum_{i=1}^{n_1} \hat{w_i}} \right ) }_{\widehat{\mu}_0^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) }_{\widehat{\mu}_0^0}, \end{align} $$

which should provide improved finite sample properties and more stable and precise estimation.Reference Busso, DiNardo and McCrary ^{64 ,} Reference Hernán, Brumback and Robins ⁶⁵ Drawing an analogy with survey sampling, the estimator in Equation (3) is a Horvitz–Thompson-type estimator and that in Equation (4) is a Hajek-type estimator, with the latter typically considered to improve the performance of the former in the literature.Reference Särndal, Swensson and Wretman ⁶⁶

In expectation, if the model in Equation (2) is correctly specified, the estimated weights ( $\hat {w_{i}}$ , for i in $1,\ldots ,n_{1}$ ) will balance the covariate distribution of the SAT with respect to that of the external control, enabling consistent estimation of mean absolute outcome $\mu ^1_0$ and the ATC. A mathematical derivation showing the consistency of the IOW estimators is provided in the Appendix.

Unfortunately, this “modeling approach” to weighting,Reference Chattopadhyay, Hase and Zubizarreta ^{67 ,} Reference Filla, Schwender and Kuss ⁶⁸ where propensity scores are explicitly modeled as a function of baseline covariates by a logistic regression, then estimated by maximizing the fit of such regression, has certain limitations:Reference Chattopadhyay, Hase and Zubizarreta ^{67 ,} Reference Filla, Schwender and Kuss ⁶⁸

• • The resulting weights do not produce adequate covariate balance if the propensity score model is misspecified, and even a correctly specified model does not guarantee balance in finite samples.

• • Propensity score predictions that are close to zero produce extreme and highly variable weights, which lead to unstable and imprecise ATC estimation, particularly where overlap is poor or the sample size of the SAT is small.

• • There is limited applicability when covariate IPD for the external control are unavailable and only marginal summary moments from published tables of baseline characteristics are available.

3.3 Entropy balancing (matching-adjusted indirect comparison)

The limitations of the weighting methods in Section 3.2 motivate alternative “balancing” or “calibration” approaches to weighting. These estimate weights under the condition that covariates are balanced, viewing balance as an optimization problem, without explicitly modeling the propensity score. Generally, balancing approaches to weighting are: (1) less susceptible to bias by directly enforcing covariate balance; (2) produce more stable weights, which translate into larger ESSs and more precise treatment effect estimation; and (3) are applicable where only aggregate-level marginal covariate moments are available for the external control.Reference Chattopadhyay, Hase and Zubizarreta ^{67 –} Reference Wang ⁷⁰

Our focus here is on an entropy balancing approachReference Hainmueller ³⁸ called MAIC,Reference Signorovitch, Wu and Yu ^{39 –} Reference Josey, Berkowitz, Ghosh and Raghavan ⁴¹ but see Chattopadhyay et al. and Filla et al. for details about similar balancing techniques.Reference Chattopadhyay, Hase and Zubizarreta ^{67 ,} Reference Filla, Schwender and Kuss ⁶⁸ MAIC is the most popular balancing method in the context of externally controlled SATs. It has many features that are considered desirable: “linear double robustness,” minimally dispersed weights, and the estimation of odds weights that are guaranteed to be positive, resulting in increased interpretability and sample boundedness, that is, interpolating the observed data as opposed to extrapolating beyond its support.Reference Filla, Schwender and Kuss ^{68 ,} Reference Tan ⁷¹ We review the main steps of MAIC, building on prior literature.Reference Signorovitch, Wu and Yu ^{39 ,} Reference Wang ^{70 ,} Reference Phillippo, Dias, Ades and Welton ^{72 –} Reference Jiang, Cappelleri, Gamalo, Chen, Thomas and Chu ⁷⁵

While MAIC does not explicitly model the propensity score as a function of baseline covariates, it implicitly assumes the following logistic model for data source assignment:

(5) $$ \begin{align} \ln(v_i) \propto \ln \Bigg ( \frac{(1-e_i)}{e_i} \Bigg ) = \gamma_0 + \mathbf{c}(\mathbf{X}_i)^\top \boldsymbol{\gamma}, \end{align} $$

where $v_i$ is a weight proportional to the inverse conditional odds of SAT participation for subject $i=1,\dots , n$ , $\gamma _0$ is an intercept term parameter, and $\boldsymbol {\gamma }$ is a vector of model parameters. In Equation (5), it is the log-odds of external control participation, $\operatorname {\mathrm {logit}}(1-e_i)$ , that are linear across the covariate balance functions. Because $\operatorname {\mathrm {logit}}(1-e_i)=-\operatorname {\mathrm {logit}}(e_i)$ , this implies that the log-odds of SAT participation, $\operatorname {\mathrm {logit}}(e_i)$ , also vary linearly with $\mathbf {c}(\mathbf {X}_i)$ as per Equation (2).

Signorovitch et al. Reference Signorovitch, Wu and Yu ³⁹ propose using the “method of moments” to estimate the model in Equation (5), such that:

(6) $$ \begin{align}\frac{ \sum_{i=1}^{n_1} v_i \mathbf{c}(\mathbf{X}_i) } { \sum_{i=1}^{n_1} v_i } &= \frac{1}{n_0} \sum_{i=n_1+1}^{n} \mathbf{c}(\mathbf{X}_i) \end{align} $$

(7)

(8) $$ \begin{align} \frac{ \sum_{i=1}^{n_1} \exp \left (\mathbf{c}(\mathbf{X}_i)^\top \boldsymbol{\gamma} \right ) \mathbf{c}(\mathbf{X}_i) } { \sum_{i=1}^{n_1} \exp \left (\mathbf{c}(\mathbf{X}_i)^\top \boldsymbol{\gamma} \right ) } &= \boldsymbol{\hat{\theta}}, \end{align} $$

where $ \boldsymbol {\hat {\theta }} = ( \boldsymbol{ \hat {\theta }}_1, \boldsymbol {\hat {\theta }}_2,\dots , \boldsymbol{ \hat{\theta }}_k)^\top$ is a vector of covariate balance function moments $j=1,\dots , k$ for the external control sample, with $\boldsymbol {\hat {\theta }}_j = \frac {1}{n_0}\sum _{i=n_1+1}^n c_j(\mathbf {X}_i)$ assumed to be a consistent estimator for  ${\boldsymbol {\theta }_j}$ . Equation (6) is a constraint enforcing that the covariate distributional features of the weighted SAT subjects are exactly balanced with respect to those of the unweighted external control subjects. Equation (7) follows from introducing the assumed model in Equation (5) into the balancing constraint, and Equation (8) results from the exponentiated intercept terms canceling out.

Replacing $\boldsymbol {\gamma }$ with estimate $\boldsymbol {\hat {\gamma }}$ in Equation (8) and centering the SAT covariate balance functions on their external control means, one obtains

(9) $$ \begin{align} \frac{ \sum_{i=1}^{n_1} \exp \left (\mathbf{c^{*}}(\mathbf{X}_i)^\top \boldsymbol{\hat{\gamma}} \right ) \mathbf{c^{*}}(\mathbf{X}_i) } { \sum_{i=1}^{n_1} \exp \left (\mathbf{c^{*}}(\mathbf{X}_i)^\top \boldsymbol{\hat{\gamma}} \right ) } = \mathbf{0}, \end{align} $$

where $\mathbf {0}$ is a vector of zeros and $\mathbf {c^{*}}(\mathbf {X}_i)= \mathbf {c(\mathbf {X}_i)} - \boldsymbol {\hat {\theta }}$ is a vector of centered covariate balance functions for subject $i=1,\dots ,n_1$ in the SAT. Then, because the denominator is positive, Equation (9) is equal to $\sum _{i=1}^{n_1} \exp \left (\mathbf {c^{*}}(\mathbf {X}_i)^\top \boldsymbol {\hat {\gamma }} \right ) \mathbf {c^{*}}(\mathbf {X}_i) = \mathbf {0}$ . Solving for $\boldsymbol {\hat {\gamma }}$ is equivalent to minimizing the objective function:

(10) $$ \begin{align} Q(\boldsymbol{\hat{\gamma}}) = \sum_{i=1}^{n_1} \exp \left ( \mathbf{c^{*}}(\mathbf{X}_i)^\top \boldsymbol{\hat{\gamma}} \right ),\end{align} $$

as the derivative of $Q(\boldsymbol {\hat {\gamma }})$ with respect to $\boldsymbol {\hat {\gamma }}$ is $\sum _{i=1}^{n_1} \exp \left (\mathbf {c^{*}}(\mathbf {X}_i)^\top \boldsymbol {\hat {\gamma }} \right ) \mathbf {c^{*}}(\mathbf {X}_i)$ . The objective function in Equation (10) is strictly convex and can be minimized using standard Newton-type convex optimization algorithms,Reference Hainmueller ³⁸ yielding a unique finite solution corresponding to the global minimum of $Q(\boldsymbol {\hat {\gamma }})$ . We have $v_i \propto \exp \left ( \mathbf {c} (\mathbf {X}_i)^\top \boldsymbol {\gamma } \right ) \propto \exp \left ( \left ( \mathbf {c}(\mathbf {X}_i)- \boldsymbol {\theta } \right ) ^\top \boldsymbol {\gamma } \right ) = \exp \left ( \mathbf {c^{*}}(\mathbf {X}_i)^\top \boldsymbol {\gamma } \right ) $ . Subject to the normalization constraint $\sum _{i=1}^{n_1} \hat {v}_i=1$ , weights for each individual $i=1,\dots ,n_1$ in the SAT are estimated as

(11) $$ \begin{align} \hat{v_i} = \frac{ \exp \left ( \mathbf{c^{*}}(\mathbf{X}_i)^\top \boldsymbol{ \hat{\gamma}} \right ) } { \sum_{i=1}^{n_1} \exp \left ( \mathbf{c^{*}}(\mathbf{X}_i)^\top \boldsymbol{\hat{\gamma}} \right ) }.\end{align} $$

Note that this definition of the entropy balancing weights coincides with the definitions presented by Jiang et al.Reference Jiang, Cappelleri, Gamalo, Chen, Thomas and Chu ⁷⁵ and by Jackson et al.Reference Jackson, Rhodes and Ouwens ⁷⁴

Similar to Equation (4), the MAIC estimator for the ATC contrasts absolute outcomes on the additive scale:

(12) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \sum_{i=1}^{n_1} \hat{v_i} Y_i \right ) } _{\widehat{\mu}_0^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) } _{\widehat{\mu}_0^0} ,\end{align} $$

where the weights have already been normalized to sum to one. Alternatively, fitting a weighted univariable regression of outcome on treatment to the concatenated IPD (with weights $\hat {v_i}$ for ${i=1,\dots ,n_1}$ and unit weights for $i=n_1+1,\dots , n$ ) has been proposed, with the treatment coefficient of the fitted model yielding an ATC estimate.Reference Ishak, Proskorovsky and Benedict ^{49 ,} Reference Park, Campbell and Towle ^{76 –} Reference Chandler and Proskorovsky ⁷⁹ We discourage this approach. The weights estimated by Equation (11) are relative; their arbitrary rescaling by a constant of proportionality, for example, a normalization constant, will also balance the specified covariate functions and does not affect $\widehat {\mu }_0^1$ and $\widehat {\text {ATC}}$ in Equation (12).Reference Phillippo, Dias, Ades and Welton ^{72 ,} Reference Jackson, Rhodes and Ouwens ⁷⁴ Nevertheless, rescaling the weights while retaining weights of one for the external control subjects results in different fitted models.Reference Jackson, Rhodes and Ouwens ⁷⁴

The method-of-moments MAIC estimator gives the Lagrangian dual solution to an entropy balancing primal problem: minimizing the negative entropy of the weights.Reference Josey, Berkowitz, Ghosh and Raghavan ^{41 ,} Reference Cheng, Tchetgen and Signorovitch ^{43 ,} Reference Phillippo, Dias, Ades and Welton ⁷² Namely, minimizing the objective function in Equation (10) is equivalent to minimizing the negative entropy $\sum _{i=1}^{n_1} v_i \ln (v_i)$ with Lagrange multipliers.Reference Cheng, Tchetgen and Signorovitch ^{43 ,} Reference Phillippo, Dias, Ades and Welton ⁷² The dual optimization problem is easier to solve than the primalReference Hainmueller ³⁸ and has been formulated in different ways.Reference Josey, Berkowitz, Ghosh and Raghavan ^{41 ,} Reference Phillippo, Dias, Ades and Welton ^{72 ,} Reference Jiang, Cappelleri, Gamalo, Chen, Thomas and Chu ⁷⁵ These may perform differently computationally—for example, if minimization is performed on the log scaleReference Phillippo, Dias, Ades and Welton ⁷² —but result in the same unique dual solution, up to numerical error, due to strict convexity. Because the negative entropy measures the distance of the weights from a uniform distribution, its minimization should produce less disperse weights than the modeling approach in Section 3.2.

MAIC (entropy balancing) has a “linear double robustness” property as it is consistent under two distinct underlying data-generating modelsReference Josey, Berkowitz, Ghosh and Raghavan ^{41 ,} Reference Cheng, Tchetgen and Signorovitch ^{43 ,} Reference Zhao and Percival ⁴⁴ : (1) the true log-odds of the propensity score are linear across the specified balance functions or (2) the true potential outcome under the active intervention is linear across the specified balance functions. Respectively, either $\operatorname {\mathrm {logit}}(e_i) = \alpha _0 + \mathbf {c}(\mathbf {X}_i)^\top \boldsymbol {\alpha }$ for some parameters $\alpha _0$ and $\boldsymbol {\alpha }$ ; or, $\text {E}(Y^1_i \mid \mathbf {X}_i) = \beta _0 + \mathbf {c}(\mathbf {X}_i)^\top \boldsymbol {\beta }$ for some parameters $\beta _0$ and $\boldsymbol {\beta }$ .

For continuous covariates, it is common practice to only balance first-order sample moments (means) by setting $\mathbf {c}(\mathbf {X}_i)=c_1(\mathbf {X}_i)=\mathbf {X}_i$ and $\boldsymbol {\hat {\theta }} = \boldsymbol {\hat {\theta }}_1 = \frac {1}{n_0}\sum _{i=n_1+1}^n \mathbf {X}_i$ , in which case consistency is guaranteed as long as either the true propensity score model is $\operatorname {\mathrm {logit}}(e_i) = \alpha _0 + \mathbf {X}_i^\top \boldsymbol {\alpha }$ ; or, the potential outcome for the active intervention varies linearly with the covariates $\mathbf {X}_i$ that are mean-balanced: $\text {E}(Y^1_i \mid \mathbf {X}_i) = \beta _0 + \mathbf {X}_i^\top \boldsymbol {\beta }$ . An alternative strategy is to also balance second-order sample moments (variances) by setting $\mathbf {c}(\mathbf {X}_i)=[c_1(\mathbf {X}_{i}), c_2(\mathbf {X}_{i})]^\top =[\mathbf {X}_{i},\mathbf {X}^2_{i}]^\top $ , $\boldsymbol {\hat {\theta }} = (\boldsymbol {\hat {\theta }}_1, \boldsymbol {\hat {\theta }}_2)^\top $ , $\hat {\boldsymbol {\theta }}_1 = \frac {1}{n_0}\sum _{i=n_1+1}^n \mathbf {X}_i$ and $\boldsymbol {\hat {\theta }}_2 = \frac {1}{n_0}\sum _{i=n_1+1}^n \mathbf {X}_i^2$ , in which case consistency is achieved if either the true propensity score model is $\text {logit}(e_i) = \alpha _0 + \mathbf {X}_i^\top \boldsymbol {\alpha _{1}} + (\mathbf {X}_i^{2})^\top \boldsymbol {\alpha _{2}}$ ; or, the potential outcome for the active intervention varies quadratically with $\mathbf {X}_i$ , such that $\text {E}(Y^1_i \mid \mathbf {X}_i) = \beta _0 + \mathbf {X}_i^\top \boldsymbol {\beta }_1 + (\mathbf {X}_i^2)^\top \boldsymbol {\beta }_2$ . Balancing the means of the covariates and the squared covariates enforces that variances are balanced because $\text {Var}(X) = \text {E}(X^2) - \text {E}(X)^2$ .

One could conjecture more flexible implicit data-generating models by considering other non-linear transformations of $\mathbf {X}_i$ , for example, higher-order polynomial terms and flexible basis functions such as splines, and balancing on the means of the transformed covariates. Moreover, one could go beyond the marginal moments of individual covariates and also balance joint covariate distributions; for instance, by including quantiles of interactions for pairs or triples of covariates.Reference Zubizarreta ⁸⁰ However, pursuing these more ambitious balancing strategies is often infeasible:

• • As the number of balancing conditions increases, it is more likely that $\boldsymbol {\theta }$ falls outside the convex hull of $\mathbf {c}(\mathbf {X}_i)$ for $i=1,\dots ,n_1$ .Reference Glimm and Yau ⁶⁰ This implies that a feasible weighting solution to the convex optimization problem does not exist: there is no set of positive weights that can enforce balance in the required distributional features and the numerical optimization algorithm will fail to converge.Reference Chattopadhyay, Hase and Zubizarreta ⁶⁷

• • Increasing the number of balancing conditions leads to further reductions in ESS and precision, which are particularly problematic with low sample sizes and poor covariate overlap.

• • Where covariate IPD for the external control are unavailable, only first- and second-order marginal moments from published tables of baseline characteristics are often available for balancing. Higher-order moments and the means of transformed covariates are rarely reported.

3.4 G-computation

An alternative to weighting involves postulating a model $m \left (\mathbf {X}_i; \boldsymbol {\beta } \right )$ for the potential outcome expectation under the active intervention, conditional on covariates $\mathbf {X}_i$ :

(13) $$ \begin{align} q\left( \text{E}(Y_i^1 \mid \mathbf{X}_i; \boldsymbol{\beta}) \right)=m \left(\mathbf{X}_i; \boldsymbol{\beta} \right ), \end{align} $$

where $q(\cdot )$ denotes an appropriate link function and $\boldsymbol {\beta }$ is a vector of model parameters encoding the covariate–outcome relationships. For instance, a logistic regression could be used for binary outcomes, such that the model is binomial, the link function is the logit, and the potential outcome expectation is on the probability scale. We have assumed that the model for the conditional outcome expectation is parametric, but this need not necessarily be the case.Reference Kern, Stuart, Hill and Green ⁸¹

The model in Equation (13) is fitted to the SAT participants $i=1,\dots , n_1$ using maximum-likelihood estimation, such that the fitted model $m \left (\mathbf {X}_i; \boldsymbol {\hat {\beta }} \right )$ with parameter estimates $\boldsymbol {\hat {\beta }}$ is an estimator of the (transformed) conditional expectation $q\left (\text {E}(Y_i^1 \mid \mathbf {X}_i, S_i=1)\right )$ . Based on the model, potential outcomes under the active intervention are predicted for each subject in the external control:

(14) $$ \begin{align} \widehat{Y}_i^1 = q^{-1} \left (m({\mathbf{X}_i};\boldsymbol{\hat{\beta}})\right). \end{align} $$

In this case, $i=n_1+1,\dots ,n$ , and the predicted outcomes are counterfactual because the subjects in the external control have not received the active intervention.

By averaging (“marginalizing”) the potential outcome predictions generated by Equation (14) over the empirical covariate distribution of the external control, one obtains the G-computation estimator,Reference Loiseau, Trichelair and He ^{35 ,} Reference Wang, Nianogo and Arah ^{45 –} Reference Remiro-Azócar, Heath and Baio ⁴⁸ also known as “regression standardization” or “plug-in G-formula,” for the mean potential outcome $\mu _0^1$ had subjects in the external control received the active intervention:

(15) $$ \begin{align} \widehat{\mu}_0^1= \frac{1} {n_0} \sum_{i=n_1+1}^{n} \widehat{Y}_i^1. \end{align} $$

Consequently, mean absolute outcomes are converted to the additive scale imposed by link function $g(\cdot )$ , prior to taking the difference between treatments on such scale to derive the G-computation estimator for the ATC:

(16) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \frac{1} {n_0} \sum_{i=n_1+1}^{n} \widehat{Y}_i^1 \right ) }_{\widehat{\mu}_0^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) }_{\widehat{\mu}_0^0}. \end{align} $$

Note that the link function $g(\cdot )$ used for summarizing the treatment effect does not necessarily need to match the link function $q(\cdot )$ used for modeling.Reference Remiro-Azócar ⁸²

The readership may be more familiar with the well-known G-computation estimator for the ATE, which contrasts averages of predicted potential outcomes between both treatment groups.Reference Ren, Cislo, Cappelleri, Hlavacek and DiBonaventura ³⁴ Our G-computation estimator in Equation (16) contrasts the average of predicted counterfactual outcomes under the active intervention with the average of observed factual outcomes for the external control group. To estimate the ATC, only potential outcomes for the SAT subjects need to be predicted because all potential outcomes under the control have already been observed for the subjects in the external control.Reference Wang, Nianogo and Arah ⁴⁵

The G-computation estimators in Equations (15) and (16) rely on the outcome model in Equation (13) being correctly specified, in order to be consistent for the mean absolute outcome ${\mu }_0^1$ and the ATC, respectively. Assuming correct model specification for all estimators, G-computation is more precise and efficient than weighting, particularly when poor overlap leads to large reductions in ESS.Reference Loiseau, Trichelair and He ^{35 ,} Reference Remiro-Azócar, Heath and Baio ⁴⁸ However, the increase in precision is achieved by implicit extrapolation into non-overlapping regions of the covariate space, hiding underlying failures of the positivity assumption.Reference Vo ⁵⁰ Model misspecification bias is almost impossible to diagnose in extrapolated regions, and there is typically no inflation of the variance to reflect the extrapolation uncertainty.

3.5 Doubly robust augmented weighting estimators

The limitations of the methods in Sections 3.2–3.4 motivate the explicit augmentation of the weighting estimators,Reference Funk, Westreich, Wiesen, Stürmer, Brookhart and Davidian ⁸³ allowing for a less restrictive outcome model that permits non-linear link functions and/or outcomes that depend on non-linear covariate transformations. Our proposed augmented weighting estimators will not use the outcome model to extrapolate, but to gain bias-robustness—and, potentially, some precisionReference Lunceford and Davidian ³⁶ —with respect to their weighting counterparts. To construct such augmented estimators, we proceed as follows.

Firstly, based on Funk et al. and Shinozaki and Matsuyama,Reference Funk, Westreich, Wiesen, Stürmer, Brookhart and Davidian ^{83 ,} Reference Shinozaki and Matsuyama ⁸⁴ we propose combining the modeling approach to inverse odds weighting, described in Section 3.2, with G-computation. Specifically, suppose we have fitted the outcome model in Equation (13) to the SAT. We now let $\hat {Y_i}^1= q^{-1} \left ( m (\mathbf {X}_i; \boldsymbol {\hat {\beta }} ) \right )$ be a prediction of potential outcome $Y_i^1$ for the active intervention based on the fitted regression, not only for the external control subjects but for all subjects $i=1,\dots ,n$ in the SAT and the external control. The propensity score weights $\hat {w_i}$ derived for $i=1,\dots , n_1$ are used to add an error-correcting term to the G-computation estimator. The resulting DR augmented IOW estimator for the mean absolute outcome $\mu _0^1$ is

(17) $$ \begin{align} \widehat{\mu}_0^1 =& \frac{1} {n_0} \left ( \sum_{i=1}^{n_1} \hat{w_i} \left ( Y_i - \hat{Y_i}^1 \right ) \right ) + \frac{1}{n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1 \nonumber \\ =& \frac{1} {n_0} \left ( \sum_{i=1}^{n_1} \hat{w_i} \epsilon_i^1 + \sum_{i=n_1+1}^{n} \hat{Y_i}^1 \right ), \end{align} $$

where $\epsilon _i^1=Y_i -\hat {Y_i}^1$ is a residual term for subject $i=1,\dots ,n_1$ in the SAT. Note that this estimator exactly corresponds to the DR estimator proposed by Shinozaki and Matsuyama,Reference Shinozaki and Matsuyama ⁸⁴ except that their estimator is for $\text {E}(Y^{0}| S = 1)$ while the estimator in Equation (17) is for $\text {E}(Y^{1}| S = 0)$ . The estimand ultimately targeted by Shinozaki and Matsuyama is the ATT but ours is the ATC.

In Equation (17), the G-computation estimator has been augmented with a weighted average of the residuals for the SAT subjects. We shall explain later in this section why such a weighted average is an error-correcting term for the potential bias of the G-computation estimator.Reference Benkeser, Carone, Laan and Gilbert ⁸⁵ The corresponding DR augmented IOW estimator for the ATC is

(18) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \frac{1} {n_0} \left ( \sum_{i=1}^{n_1} \hat{w_i} \epsilon_i^1 + \sum_{i=n_1+1}^{n} \hat{Y_i}^1 \right ) \right ) } _{\widehat{\mu}_0^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) }_{\widehat{\mu}_0^0},\end{align} $$

on the additive scale imposed by link function $g(\cdot )$ . As for the non-augmented IOW estimators in Section 3.2, we can normalize or stabilize the weights so that they sum to one, thereby ensuring bounded estimates and improving the finite sample properties of the estimator.Reference Dahabreh, Robertson, Steingrimsson, Stuart and Hernan ⁶¹ Using the normalized weights, we obtain the DR augmented normalized IOW estimators:

(19) $$ \begin{align} \widehat{\mu}_0^1 = \frac{ \sum_{i=1}^{n_1} \hat{w_i} \epsilon_i^1 } { \sum_{i=1}^{n_1} \hat{w_i} } + \frac{1}{n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1, \end{align} $$

(20) $$ \begin{align}\widehat{\text{ATC}} = g \underbrace { \left ( \frac{ \sum_{i=1}^{n_1} \hat{w_i} \epsilon_i^1 } { \sum_{i=1}^{n_1} \hat{w_i} } + \frac{1}{n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1 \right ) } _{\widehat{\mu}_0^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) }_{\widehat{\mu}_0^0}. \end{align} $$

Equations (19) and (20) should provide more stable and precise estimation than Equations (17) and (18), respectively.Reference Dahabreh, Robertson, Steingrimsson, Stuart and Hernan ⁶¹

Our main contribution is combining the entropy balancing-based MAIC approach, described in Section 3.3, with the G-computation estimator. Again, based on the outcome model fitted to the SAT, let $\hat {Y_i}^1=q^{-1} \left ( m (\mathbf {X}_i; \boldsymbol {\hat {\beta }} ) \right )$ be a prediction of the potential outcome $Y_i^1$ under the active intervention for all subjects $i=1,\dots ,n$ in the SAT and the external control. We now employ the MAIC weights $\hat {v_i}$ for $i=1,\dots ,n_1$ derived in Equation (11), instead of the “modeling” IOW weights, to estimate the error-correcting term used to augment the G-computation estimator. The resulting DR augmented MAIC estimator of the mean absolute outcome $\mu _0^1$ is

(21) $$ \begin{align} \widehat{\mu}_0^1 =& \sum_{i=1}^{n_1} \hat{v _i } \left ( Y_i - \hat{Y_i}^1 \right ) + \frac{1}{n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1 \nonumber \\ =& \sum_{i=1}^{n_1} \hat{v _i} \epsilon_i^1 + \frac{1}{n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1,\end{align} $$

where the G-computation estimator has been augmented with a weighted average of the residuals ${\epsilon _i^1=Y_i -\hat {Y_i}^1}$ for $i=1,\dots ,n_1$ , but this time the weighted average has been computed using the MAIC weights.

The corresponding DR augmented MAIC estimator for the ATC is

(22) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \sum_{i=1}^{n_1} \hat{v_i} \epsilon_i^1 + \frac{1}{n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1 \right ) } _{\widehat{\mu}_0^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) }_{\widehat{\mu}_0^0}, \end{align} $$

on the additive scale imposed by link function $g(\cdot )$ . We conjecture that the DR augmented MAIC estimators in Equations (21) and (22) will perform better statistically than the augmented estimators based on the modeling approach to weighting, which could exhibit erratic performance with highly variable weights, particularly if these are combined with a misspecified outcome model.Reference Robins, Sued, Lei-Gomez and Rotnitzky ^{86 ,} Reference Kang and Schafer ⁸⁷ We expect the DR augmented MAIC estimators to inherit the more attractive properties of the entropy balancing weights: (1) lower susceptibility to bias by directly enforcing covariate balance and (2) greater stability, translating into larger ESSs after weighting and enhanced precision in estimation.

The augmented weighting estimators in Equations (17), (19), and (21) are DR for the mean potential outcome $\mu ^1_0$ . That is, they estimate $\mu ^1_0$ consistently as long as either the propensity score model for data source assignment or the outcome model is correctly specified, but not necessarily both. Note that, contrary to the approach that we shall describe in Section 3.6, this double robustness does not depend on using a canonical link function for the outcome model. In the Appendix, we provide an intuitive heuristic to demonstrate double robustness, which illustrates why the weighted average of residuals is an error-correcting term for the potential bias of G-computation. The heuristic is summarized as follows. Consider that all the augmented weighting estimators have the general form:

(23) $$ \begin{align} \widehat{\mu}_0^1 = \sum_{i=1}^{n_1} \hat{u_{i}} \left ( Y_i - \hat{Y_i }^1\right ) + \frac{1}{n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1, \end{align} $$

for a generic weight estimate $\hat {u_i}$ , where $\hat {u_i}=\hat {w_i}/n_0$ for the DR augmented IOW estimator in Equation (17), $\hat {u}_i=\hat {w}_i/ \sum _{i=1}^{n_1} \hat {w}_i$ for the DR augmented normalized IOW estimator in Equation (19), and $\hat {u_i}=\hat {v_i}$ for the DR augmented MAIC estimator in Equation (21).

If the outcome model is correctly specified, the expectation of the first summation in Equation (23) converges to zero as $n_1 \rightarrow \infty $ because $\hat {Y_i}^1 \rightarrow Y_{i}$ and the terms inside the summation cancel out, irrespective of any postulated propensity score model for data source assignment. The second summation is equivalent to the G-computation estimator and is consistent for $\mu _0^1$ because the outcome model is correct. Consequently, $\widehat {\mu }^1_0 \rightarrow \mu _0^1$ and $\widehat {\text {ATC}} \rightarrow \text {ATC}$ , under the assumption that $\widehat {\mu }_0^0=\frac {1}{n_0}\sum _{i=n_1+1}^n Y_i$ is consistent for $\mu _0^0$ . Conversely, if the propensity score model is correctly specified but the outcome model is incorrect, the first summation consistently cancels out the bias produced by the G-computation estimator in the second summation and the remainder term is exactly equal to the non-augmented weighting estimator, which converges to $\mu _0^1$ as $n_1 \rightarrow \infty $ because the propensity score model is correct. Similarly, $\widehat {\text {ATC}} \rightarrow \text {ATC}$ (assuming $\widehat {\mu }_0^0 \rightarrow \mu _0^0$ ). Hence, the augmented weighting estimators in Equations (18), (20), and (22) are DR for the ATC.

While all the augmented weighting estimators described in this section are DR, the DR augmented MAIC estimators defined in Equations (21) and (22) are arguably more robust to model misspecification bias because they are consistent under a greater number of distinct underlying data-generating mechanisms. Namely, DR augmented MAIC is consistent as long as either: (1) the log-odds of the propensity score are linear across the covariate balance functions; (2) the potential outcome under the active intervention is linear across the covariate balance functions; or (3) the explicit augmentation model for the potential outcome under the active intervention is correctly specified. The estimation of the weights is consistent as long as either the first or the second condition holds. Conversely, the augmented estimators based on the modeling approach to weighting (defined in Equations (17)–(20)) are consistent as long as either the first or the third condition holds, with the first condition being necessary for consistent estimation of the weights.

3.6 Weighted G-computation

In Section 3.5, we proposed augmented estimators that combine the predictions of an unweighted outcome model with weights in a weighted average. Nevertheless, there are other ways of constructing augmented weighting estimators. One approach popularized by Park et al.Reference Park, Campbell and Towle ⁷⁶ that has been claimed to be DR consists of G-computation based on the predictions of a weighted outcome model.Reference Dahabreh, Robertson, Steingrimsson, Stuart and Hernan ^{61 ,} Reference Park, Campbell and Towle ⁷⁶ Where the target estimand is the ATC, this involves: (1) estimating weights using the methods described in Sections 3.2 and 3.3; (2) fitting a weighted model for the conditional outcome expectation to the SAT participants; and (3) marginalizing the outcome predictions of the weighted regression over the external control covariate distribution. The resulting estimator for the mean absolute outcome $\mu _0^1$ would be

(24) $$ \begin{align} \widehat{\mu}_0^1= \frac{1} {n_0} \sum_{i=n_1+1}^{n} \hat{Y_i}^1 = \frac{1} {n_0} \sum_{i=n_1+1}^{n} q^{-1} \left (m({\mathbf{X}_i};\boldsymbol{\hat{\beta_v}})\right), \end{align} $$

where $m({\mathbf {X}_i};\hat {\boldsymbol {\beta_v }})$ indexes the fitted weighted regression with vector $\hat {\boldsymbol {\beta }}_v$ of parameter estimates. The ATC would be estimated by substituting Equation (24) into Equation (16).

Such estimators are only DR where the outcome model is a generalized linear model (GLM) with a canonical link function $q(\cdot )$ .Reference Dahabreh, Robertson, Steingrimsson, Stuart and Hernan ^{61 ,} Reference Robins, Sued, Lei-Gomez and Rotnitzky ^{86 ,} Reference Gabriel, Sachs and Martinussen ^{88 ,} Reference Wooldridge ⁸⁹ The estimator in Equation (24) and the corresponding ATC estimator would not be DR where the GLM link function is non-canonical,Reference Gabriel, Sachs and Martinussen ⁸⁸ or where the outcome model is a Cox proportional hazards model or a parametric survival model in the time-to-event setting.Reference Gabriel, Sachs and Waernbaum ⁹⁰ Nevertheless, results by Gabriel et al. suggest asymptotic equivalence and similar finite sample performance to the augmented weighting estimators in Section 3.5 for GLMs with canonical link functions fitted via maximum-likelihood,Reference Gabriel, Sachs and Martinussen ⁸⁸ provided that the same weights are used and correct model specification. We note that the target of the investigations by Gabriel et al. is the ATE and the modeling approach to weighting.Reference Gabriel, Sachs and Martinussen ⁸⁸

3.7 Variance estimation

To estimate the variance and construct confidence intervals (CIs) for $\widehat {\mu }_{0}^1$ and $\widehat {\text {ATC}}$ , it is possible to use empirical sandwich-type (“robust”) variance estimators to account for the correlation induced by weighting.Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{25 ,} Reference Signorovitch, Wu and Yu ³⁹ In the specific context of non-randomized comparisons, such as the externally controlled SATs and unanchored ITCs explored in this article, these estimators have exhibited either under-precision or over-precision for the ATT (or the ATC)Reference Reifeis and Hudgens ^{91 ,} Reference Kostouraki, Hajage and Rachet ⁹² and under-precision for the ATE.Reference Austin ^{37 ,} Reference Hernán, Brumback and Robins ⁶⁵ This is because most implementations ignore the estimation of the propensity score model or the weights, assuming the weights to be fixed quantities.Reference Reifeis and Hudgens ^{91 ,} Reference Kostouraki, Hajage and Rachet ⁹²

Analytic expressions that incorporate weight estimation could be derived,Reference Cheng, Ayyagari and Signorovitch ^{40 ,} Reference Reifeis and Hudgens ^{91 ,} Reference Kostouraki, Hajage and Rachet ⁹² but we propose a practical alternative based on the ordinary non-parametric bootstrap,Reference Austin ³⁷ explicitly accounting for uncertainty in the weight estimation. This involves resampling with replacement the concatenated IPD consisting of the SAT and external control data. In each bootstrap iteration, the weight estimation and/or outcome modeling procedures are performed, and $\mu _{0}^1$ , $\mu _{0}^0$ and $\text {ATC}$ are re-estimated. Standard errors for $\widehat {\mu }_{0}^1$ , $\widehat {\mu }_{0}^0$ , and $\widehat {\text {ATC}}$ are given by the standard deviations across the bootstrap resamples. Subsequently, Wald-type CIs can be constructed. Alternatively, one can directly calculate CIs from the percentiles of the bootstrap resamples, for example, 2.5% and 97.5% for the 95% CI.

3.8 External controls with unavailable individual participant data

In the context of unanchored ITCs, the external control is often a historical trial for which IPD are unavailable, due to privacy or confidentiality reasons. In this case, only published AD are available for the external control.Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{25 ,} Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{28 ,} Reference Lambert, Lengliné, Porcher, Thiébaut, Zohar and Chevret ³³ Such data consist of marginal summary moments $\boldsymbol {\hat{ \theta }}$ from reported tables of baseline characteristics, typically only including means and standard deviations (for continuous covariates), and an estimate $\widehat {\mu }_0^0$ of the mean absolute outcome under the control in the external data source.Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{25 ,} Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{28 ,} Reference Signorovitch, Wu and Yu ^{39 ,} Reference Wang ^{70 ,} Reference Jackson, Rhodes and Ouwens ⁷⁴ An important shortfall of this scenario is the need to assume that $\hat {\boldsymbol {\theta_j }} = \boldsymbol {\theta }_j$ with zero variability for the covariate balance function moments $j=1,\dots ,k$ , that is, that the external control covariate distributional data are fixed.Reference Josey, Berkowitz, Ghosh and Raghavan ⁴¹ While this may be reasonable with large sample sizes for the external control, it can otherwise result in overly precise inferences and inflated Type I error rates.Reference Josey, Berkowitz, Ghosh and Raghavan ⁴¹

In this setting, for all methods except (non-augmented) MAIC, one must simulate M individual-level covariate profiles from the assumed covariate distribution of the external control based on published summary statistics.Reference Remiro-Azócar, Heath and Baio ^{48 ,} Reference Remiro-Azócar ^{82 ,} Reference Ren, Ren, Welton and Strong ⁹³ The number M of hypothetical subject profiles should be relatively large, for example, $M=1,000$ , to minimize sampling variability and random seed sensitivity, and does not necessarily need to match the original sample size $n_0$ of the external control.Reference Remiro-Azócar, Heath and Baio ^{48 ,} Reference Remiro-Azócar ^{82 ,} Reference Ren, Ren, Welton and Strong ⁹³ Necessary information to infer the joint covariate distribution of the external control, for example, distributional forms and correlation structures, is rarely published. Hence, this must be borrowed from other data sources or selected based on theoretical properties, following recommendations in the literature.Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{25 ,} Reference Phillippo, Ades, Dias, Palmer, Abrams and Welton ^{28 ,} Reference Remiro-Azócar, Heath and Baio ^{48 ,} Reference Remiro-Azócar ^{82 ,} Reference Ren, Ren, Welton and Strong ⁹³

The notation and procedures for Sections 3.2 and 3.4–3.6 change as follows. The observed IPD for the SAT are stacked with the simulated subject-level covariate data for the external control. The concatenated dataset is now $(S_i, \mathbf {X}_i, T_i, Y_i)$ for $i=1,\dots , n_1, n_1+1, \dots , n_1+M$ . For the SAT subjects $i=1,\dots , n_1$ , we have $S_i=1$ and $T_i=1$ , with $\mathbf {X}_i$ and $Y_i$ corresponding to the actual covariate and outcome values observed in the trial. For the hypothetical external controls $i=n_1+1,\dots ,n_1+M$ , we have $S_i=0$ and $T_i=0$ , the values of $\mathbf {X}_i$ are simulated, and $Y_i$ are unavailable but not required for the analysis. This is because the target estimand is the ATC and the outcomes for the external group under the control have already been observed factually, with the mean estimate $\widehat {\mu }_0^0$ available from published results.

The general form of the IOW estimators for the ATC, described in Section 3.2, is now:

$$ \begin{align*} \widehat{\text{ATC}} = g \underbrace { \left ( \frac{1}{K} \sum_{i=1}^{n_1} \hat{w_i} Y_i \right ) } _{\widehat{\mu}_0^1} - g \left ( \widehat{\mu}_0^0 \right ), \end{align*} $$

where K is a constant. There is only a change in notation here given that individual-level outcomes under the control are now unavailable for the external data source, and cannot be included in the concatenated dataset.

For the G-computation estimator for the ATC, outlined in Section 3.4, we now have

(25) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \frac{1} {M} \sum_{i=n_1+1}^{n_1+M} \hat{Y_i}^1 \right ) }_{\widehat{\mu}_0^1} - g \left ( \widehat{\mu}_0^0 \right ), \end{align} $$

where the potential outcome predictions $\hat {Y_i}^1$ under the active intervention are generated for each hypothetical external control subject $i=n_1+1,\dots ,n_1+M$ , and averaged over the simulated covariate profiles.

The general form of the DR augmented weighting estimators for the ATC, proposed in Section 3.5, is now:

$$ \begin{align*} \widehat{ \text{ ATC}} = g \underbrace{ \left ( \frac{1}{K} \sum_{i=1}^{n_1} \hat{\omega_i} \epsilon_i^1 + \frac{1}{M} \sum_{i=n_1+1}^{n_1+M} \hat{Y_i}^1 \right ) }_{\widehat{\mu}_0^1} - g \left ( \widehat{\mu}_0^0 \right ), \end{align*} $$

where K is a constant, $\epsilon _i^1=Y_i -\hat {Y_i}^1$ and $\hat {\omega }_i$ are a residual term and a weight estimate, respectively, for $i=1,\dots ,n_1$ , and the potential outcome predictions $\hat {Y_i}^1$ under the active intervention are generated for all SAT subjects and hypothetical external controls $i=1,\dots ,n_1,n_1+1,\dots , n_1+M$ .

For the weighted G-computation estimator in Section 3.6, the outcome predictions would be averaged over the simulated covariates for the external control. The resulting estimator for the mean absolute outcome $\mu _0^1$ is $\widehat {\mu }_0^1= \frac {1} {M} \sum _{i=n_1+1}^{n_1+M} \hat {Y_i}^1 = \frac {1} {M} \sum _{i=n_1+1}^{n_1+M} q^{-1} \left (m({\mathbf {X}_i};\hat {\boldsymbol {\beta }}_v)\right )$ , which would be input into Equation (25) for estimation of the ATC.

The unavailability of IPD for the external control entails some changes to the non-parametric bootstrap procedure described in Section 3.7 for variance estimation. In this case, only the SAT data, $(S_i, \mathbf {X}_i, T_i, Y_i)$ for $i=1,\dots , n_1$ , are resampled to re-estimate $g(\mu _0^1)$ in each bootstrap iteration, with the standard error, $\text {SE} \left (g(\widehat {\mu }^1_0)\right )$ , computed as the standard deviation over the bootstrap resamples. Then, the decomposition

(26) $$ \begin{align} \text{SE} \left(\widehat{\text{ATC}}\right ) = \sqrt{\big(\text{SE} \left ( g(\widehat{\mu}^1_0) \right )\big)^2 + \big(\text{SE} \left ( g(\widehat{\mu}^0_0) \right )\big)^2} \end{align} $$

is used to estimate the standard error of the ATC, where $\text {SE} \left (g(\widehat {\mu }^0_0) \right )$ is derived from published aggregate-level results.Reference Remiro-Azócar, Heath and Baio ^{48 ,} Reference Chandler and Proskorovsky ^{79 ,} Reference Ren, Ren, Welton and Strong ⁹³ A limitation of the above formula is that it assumes that the mean absolute outcomes are statistically independent. Moreover, while computing $\text {SE} \left (g(\widehat {\mu }^0_0) \right )$ is a trivial exercise for continuous and binary outcomes (e.g., there is a closed-form formula for the standard error of the log-odds using the Delta method), it can be challenging for other outcomes such as those in the time-to-event setting.Reference Chandler and Proskorovsky ^{79 ,} Reference Ren, Ren, Welton and Strong ⁹³ Once $\text {SE} \left (\widehat {\text {ATC}}\right )$ is computed, Wald-type CIs can be readily constructed.

4.1 Aims

We conducted a simulation study to evaluate the performance of various estimators under different conditions. The simulation study design was planned following the structured “ADEMP” approach outlined by Morris et al.,Reference Morris, White and Crowther ⁹⁴ to ensure reproducibility and meaningful conclusions. Specifically, we clearly defined research aims, data-generating mechanisms under controlled scenarios and estimands, and assessed the performance of several estimators using relevant performance measures: bias, empirical standard error (ESE), and coverage. All simulations and analyses were performed using R statistical software version 4.3.1. ⁹⁵ The files and code required to run the simulations are publicly available on Github at https://github.com/harlanhappydog/DRAWE-.

4.2 Data-generating mechanisms

We simulated data inspired by the data-generating mechanisms in a simulation study by Kang and Schafer.Reference Kang and Schafer ⁸⁷ Some modifications were required since Kang and Schafer considered continuous-valued outcomes,Reference Kang and Schafer ⁸⁷ while we consider binary outcomes. The simulated data consist of $(\mathbf {X}_i, \mathbf {Z}_i, T_i, S_i, Y_i)$ for $i=1,\ldots ,n$ , with the control group fully external such that $S_i=T_i$ , with $n_1= \sum _{i=1}^{n}{S_{i}}$ , and $n_0 =\sum _{i=1}^{n}{(1-S_{i})}$ , as detailed in Section 3.1. While $\mathbf {X}_i$ is observed, $\mathbf {Z}_i$ is unobserved. To generate the data, $\mathbf {X}_i$ is distributed as $\mathrm {Normal}(0,I_4)$ , for i in $1,\ldots ,n$ , and $\mathbf {Z}_i$ is then obtained by applying the following transformations:

$$ \begin{align*} Z_{i1} &= \text{scale}(\exp(X_{i1}/2)), \\ Z_{i2} &= \text{scale}(X_{i2}^{2}),\\ Z_{i3} &= \text{scale}((X_{i1}X_{i3} + 0.6)^3), \\ Z_{i4} &= \text{scale}((X_{i2} + X_{i4} + 20)^2), \end{align*} $$

where $\text {scale}(\dot )$ indicates normalization such that $Z_{i1}, Z_{i2}, Z_{i3}$ , and $Z_{i4}$ each have mean of $0$ and standard deviation of $1$ , that is, $Z = \text {scale}(f(X)) = ( f(X) - \text {mean}(f(X)))/{\text {sd}(f(X))}$ . Note that these transformations are similar to the ones detailed by Kang and SchaferReference Kang and Schafer ⁸⁷ but not identical, with changes made to highlight the consequences of model misspecification.

We consider four different scenarios. For each, we generated 10,000 simulated datasets of size $n=200$ and $n=1,000$ . Note that, in all four scenarios, the distribution of S is approximately balanced such that $n_{1} \approx n_{0}$ . The four scenarios are defined as:

• • KS1: $Y_i$ is generated from a Bernoulli distribution with

  $\text {Pr}(Y_i=1 \mid \mathbf {X}_i, T_i) = \mathrm {expit}({X}_{1i} - 1.50{X}_{2i} + 0.5{X}_{3i} - 0.5{X}_{4i} + 1.50T_{i} - 0.50T_{i}{X}_{1i})$ ,

  where $T_i=S_i$ , and $S_{i}$ is generated from a Bernoulli distribution with

  $\text {Pr}(S_i=1 \mid \mathbf {X}_i) = \mathrm {expit}(-X_{i1} + 0.5 X_{i2} - 0.25 X_{i3} -0.5 X_{i4})$ .

  The distribution of the covariates is such that overlap between the two groups is relatively high, with overlap proportions of 0.68, 0.85, 0.92, and 0.85 for $X_{1}$ , $X_{2}$ , $X_{3}$ , and $X_{4}$ , respectively (see Figure 3 in the Supplementary Material).

• • KS2: $Y_i$ is generated from a Bernoulli distribution with

  $\text {Pr}(Y_i=1 \mid \mathbf {Z}_i, T_i) = \mathrm {expit}({Z}_{1i} - 1.50{Z}_{2i} + 0.5{Z}_{3i} - 0.5{Z}_{4i} + 1.50T_{i} - 0.50T_{i}{Z}_{1i})$ ,

  where $T_i=S_i$ , and $S_{i}$ is generated from a Bernoulli distribution with

  $\text {Pr}(S_i=1 \mid \mathbf {X}_i) = \mathrm {expit}(-X_{i1} + 0.5 X_{i2} - 0.25 X_{i3} -0.5 X_{i4})$ .

  The relevant covariate adjustment approaches would fit an outcome model to the observed $\mathbf {X}_i$ , as the $\mathbf {Z}_i$ used for the true outcome-generating process are unobserved. The distribution of the covariates is such that overlap between the two groups is relatively high, with overlap proportions of 0.68, 0.85, 0.92, and 0.85 for $X_{1}$ , $X_{2}$ , $X_{3}$ , and $X_{4}$ , respectively (see Figure 3 in the Supplementary Material).

• • KS3: $Y_i$ is generated from a Bernoulli distribution with

  $\text {Pr}(Y_i=1 \mid \mathbf {X}_i, T_i) = \mathrm {expit}({X}_{1i} - 1.50{X}_{2i} + 0.5{X}_{3i} - 0.5{X}_{4i} + 1.50T_{i} - 0.50T_{i}{X}_{1i})$ ,

  where $T_i=S_i$ , and $S_{i}$ is generated from a Bernoulli distribution with

  $\text {Pr}(S_i=1 \mid \mathbf {Z}_i) = \mathrm {expit}(-Z_{i1} + 0.5 Z_{i2} - 0.25 Z_{i3} - 0.5 Z_{i4})$ .

  The relevant covariate adjustment approaches would balance or fit a propensity score model to the observed $\mathbf {X}_i$ , as the $\mathbf {Z}_i$ used for the true data source assignment process are unobserved. The distribution of the covariates is such that overlap between the two groups is relatively high, with overlap proportions of 0.71, 0.84, 0.99, and 0.89 for $X_{1}$ , $X_{2}$ , $X_{3}$ , and $X_{4}$ , respectively (see Figure 4 in the Supplementary Material).

• • KS4: $Y_i$ is generated from a Bernoulli distribution with

  $\text {Pr}(Y_i=1 \mid \mathbf {Z}_i, T_i) = \mathrm {expit}({Z}_{1i} - 1.50{Z}_{2i} + 0.5{Z}_{3i} - 0.5{Z}_{4i} + 1.50T_{i} - 0.50T_{i}{Z}_{1i})$ ,

  where $T_i=S_i$ , and $S_{i}$ is generated from a Bernoulli distribution with

  $\text {Pr}(S_i=1 \mid \mathbf {Z}_i) = \mathrm {expit}(-Z_{i1} + 0.5 Z_{i2} - 0.25 Z_{i3} - 0.5 Z_{i4})$ .

  The distribution of the covariates is such that overlap between the two groups is relatively high, with overlap proportions of 0.71, 0.84, 0.99, and 0.89 for $X_{1}$ , $X_{2}$ , $X_{3}$ , and $X_{4}$ , respectively (see Figure 4 in the Supplementary Material).

We assume that there is unlimited access to subject-level data for the SAT and the external control, such that $(\mathbf {X}_i, T_i, S_i, Y_i)$ are observed for all $i=1,\ldots ,n$ .

4.3 Estimands

The estimand of interest is the ATC, as defined in Section 2. We adopt the logit link function ${g(p) =\ln (p/(1-p))}$ for marginal outcome probability p, such that the ATC is on the marginal log-odds ratio scale. The values of the ATC estimands were calculated numerically, by simulating 10 million binary outcomes using the true data-generating mechanisms outlined in Section 4.2. Data-generating mechanisms KS1, KS2, KS3, and KS4 correspond to true ATCs of 1.116, 1.215, 1.068, and 1.181, respectively.

4.4 Methods

We compared 16 estimators:

1. 1. The naïve estimator, which does not perform covariate adjustment:

   (27) $$ \begin{align} \widehat{\text{ATC}} = g \underbrace { \left ( \frac{1}{n_1}\sum_{i=1}^{n_1} Y_i \right ) }_{\widehat{\mu}_1^1} - g \underbrace { \left ( \frac{1}{n_0}\sum_{i=n_1+1}^n Y_i \right ) }_{\widehat{\mu}_0^0}.\end{align} $$

2. 2. The IOW estimator as per Equation (3) (Section 3.2), with weights derived using the “modeling” approach with a logistic regression data source assignment model as per Equation (2).

3. 3. The normalized IOW estimator with normalized weights as per Equation (4) (Section 3.2) derived using the “modeling” approach with a logistic regression data source assignment model as per Equation (2).

4. 4. The MAIC (entropy balancing) estimator, as per Equation (12) (Section 3.3).

5. 5. The G-computation estimator, as per Equation (16) (Section 3.4) with the outcome model defined with a logistic link function.

6. 6. The DR augmented IOW estimator, as per Equation (18) (see Section 3.5) with the outcome model defined with a logistic link function.

7. 7. The DR augmented normalized IOW estimator, as per Equation (20) (see Section 3.5) with the outcome model defined with a logistic link function.

8. 8. The DR augmented MAIC (augmented entropy balancing) estimator, as per Equation (22) (our main contribution, see Section 3.5) with the outcome model defined with a logistic link function.

9. 9. The weighted G-computation (normalized IOW weights) estimator described in Section 3.6 with the outcome model defined with a logistic link function, and using the (normalized) “modeling” IOW weights.

10. 10. The weighted G-computation (MAIC weights) estimator described in Section 3.6 with the outcome model defined with a logistic link function, but using the MAIC (entropy balancing) weights instead of the “modeling” IOW weights.

11. 11–16. Estimators 11–16 are the same as Estimators 5–10 but using a non-canonical Cauchit link function, $q(x)= \text {tan}(\pi (x-0.5))$ . See Morgan and SmithReference Morgan and Smith ⁹⁶ for an example of model fitting with the Cauchit.

To be clear, for Estimators 1–10, both the outcome model and the propensity score model for data source assignment will be correctly specified in KS1. For Estimators 11–16, the propensity score model for data source assignment will be correctly specified in KS1, but not the outcome model. In KS2, the outcome model will be incorrectly specified for all estimators fitting an outcome model (Estimators 5–16). Finally, in KS3, the propensity score model will be incorrectly specified for all estimators considering a propensity score model, and in KS4, both the outcome model and the propensity score model will be incorrectly specified.

For all estimators, we computed 95% CIs using the non-parametric bootstrap approach described in Section 3.7. Specifically, we used $B=200$ bootstrap resamples of the concatenated SAT and external control IPD to approximate the standard error of $\widehat {\text {ATC}}$ and subsequently constructed Wald-type CIs.

4.5 Performance measures

To assess the performance of the estimators in our simulation study, we computed several key metrics: bias, ESE, and 95% CI coverage. Bias was calculated as the difference between the average of the point estimates across simulations and the true estimand value, providing a measure of systematic error. The ESE was computed as the standard deviation of the point estimates across simulations, reflecting the precision or variability of the different estimators. The 95% CI coverage was determined as the proportion of simulated datasets in which the constructed 95% CI contained the true estimand value, evaluating the quality of interval estimation. In addition, we estimated Monte Carlo standard errors (MCSEs) using the formulas provided by Morris et al. to quantify the uncertainty in the performance measures due to using a finite number of simulations.Reference Morris, White and Crowther ⁹⁴

4.6 Results

The complete results for the simulation study are displayed in Tables 1–4. Note that, for the $n=200$ simulations, the MCSEs are less than 0.018 for bias, less than 0.013 for ESE, and less than 0.006 for 95% CI coverage; for the $n=1,000$ simulations, the MCSEs are less than 0.011 for bias, less than 0.008 for ESE, and less than 0.005 for 95% CI coverage.

Table 1

Results from Scenario KS1, where both the logit-link outcome model and the propensity score model are correctly specified

Table 1 Long description

The table presents results for Scenario K S 1 where logit-link and propensity score models are correctly specified. It is divided into two primary sections based on sample size.

Section 1: n equals 200, A T C equals 1.116.

* The naive estimator shows the highest bias at 0.623 and lowest coverage at 0.539.

* I O W, Normalized I O W, and M A I C show bias ranging from 0.027 to 0.055 with coverage near 0.94.

* Methods 5 through 10 (Logit link outcome model) show low bias (0.019 to 0.031) and stable coverage (0.942 to 0.956).

* Methods 11 through 16 (Cauchit link outcome model) show similar performance, though Weighted G-computation with normalized I O W weights has a higher E S E of 0.914 and wider C I of 2.696.

Section 2: n equals 1,000, A T C equals 1.116.

* The naive estimator remains highly biased at 0.604 with very poor coverage at 0.008.

* All other methods (2 through 16) show significantly reduced bias, mostly below 0.010, and coverage levels closely approaching the nominal 0.95.

* G-computation and D R augmented methods show the lowest E S E values, ranging from 0.150 to 0.180.

* Under the Cauchit link, G-computation and Weighted G-computation methods show slight negative bias (minus 0.020 to minus 0.003).

Table 2

Results from Scenario KS2, where the outcome model is incorrectly specified and the propensity score model is correctly specified

Table 2 Long description

The table presents results for Scenario K S 2 where the outcome model is incorrectly specified and the propensity score model is correctly specified. Both sample size blocks share an A T C value of 1.215.

Block 1: n = 200

* 1. The naive estimator: Bias 0.230, E S E 0.324, 95% C I coverage 0.909, Average width 1.279.

* 2. I O W: Bias 0.023, E S E 0.642, 95% C I coverage 0.929, Average width 2.393.

* 3. Normalized I O W: Bias 0.060, E S E 0.511, 95% C I coverage 0.940, Average width 1.952.

* 4. M A I C: Bias 0.061, E S E 0.517, 95% C I coverage 0.965, Average width 2.778.

* Logit link methods (5-10): Bias ranges from 0.036 to 0.088; C I coverage ranges from 0.937 to 0.953.

* Cauchit link methods (11-16): Bias ranges from 0.039 to 0.225; C I coverage ranges from 0.950 to 0.972.

Block 2: n = 1,000

* 1. The naive estimator: Bias 0.226, E S E 0.141, 95% C I coverage 0.658, Average width 0.557.

* 2. I O W: Bias 0.020, E S E 0.283, 95% C I coverage 0.950, Average width 1.085.

* 3. Normalized I O W: Bias 0.016, E S E 0.220, 95% C I coverage 0.939, Average width 0.834.

* 4. M A I C: Bias 0.015, E S E 0.205, 95% C I coverage 0.936, Average width 0.779.

* Logit link methods (5-10): Bias ranges from 0.014 to 0.075; C I coverage ranges from 0.935 to 0.942.

* Cauchit link methods (11-16): Bias ranges from 0.012 to 0.228; C I coverage ranges from 0.817 to 0.942.

Table 3

Results from Scenario KS3, where the propensity score model is incorrectly specified. The logit-link outcome model is correctly specified; however, when the Cauchit link is used for the outcome model, both models are incorrectly specified

Table 3 Long description

The table presents results for Scenario K S 3. It is divided into two primary sections based on sample size.

Section 1: n equals 200, A T C equals 1.068.

* Methods 1 to 4 (Naive, I O W, Normalized I O W, M A I C) show bias ranging from minus 0.032 to 0.147.

* With logit link (Methods 5 to 10): G-computation has the lowest bias at 0.015. D R augmented methods and Weighted G-computation show bias between 0.018 and 0.028.

* With Cauchit link (Methods 11 to 16): Bias ranges from minus 0.031 to 0.018. Weighted G-computation with normalized I O W weights shows a high E S E of 0.701 and C I width of 2.176.

Section 2: n equals 1,000, A T C equals 1.068.

* Methods 1 to 4: Bias ranges from minus 0.044 to 0.118. M A I C has the lowest 95% C I coverage at 0.894.

* With logit link (Methods 5 to 10): All methods show very low bias (0.005 to 0.007) and consistent C I coverage around 0.94.

* With Cauchit link (Methods 11 to 16): Bias is generally negative, ranging from minus 0.053 to minus 0.010. G-computation (Method 11) shows the highest negative bias at minus 0.053.

Table 4

Results from Scenario KS4, where both the outcome model and the propensity score model are incorrectly specified

Table 4 Long description

The table compares 16 methods across two sample sizes, n equals 200 and n equals 1,000, where A T C equals 1.181. Columns include Bias, E S E, 95% C I coverage, and Average 95% C I width.

For n equals 200:

* The naive estimator: Bias 0.521, E S E 0.334, Coverage 0.685.

* I O W: Bias 0.807, E S E 0.778, Coverage 0.912.

* Normalized I O W: Bias 0.582, E S E 0.470, Coverage 0.760.

* M A I C: Bias 0.613, E S E 0.463, Coverage 0.790.

* Under Logit link (Methods 5-10): D R augmented M A I C shows the lowest bias at 0.517.

* Under Cauchit link (Methods 11-16): D R augmented M A I C shows the lowest bias at 0.400.

For n equals 1,000:

* The naive estimator: Bias 0.495, E S E 0.146, Coverage 0.069.

* I O W: Bias 0.778, E S E 0.350, Coverage 0.340.

* Normalized I O W: Bias 0.514, E S E 0.197, Coverage 0.256.

* M A I C: Bias 0.546, E S E 0.190, Coverage 0.174.

* Under Logit link (Methods 5-10): D R augmented M A I C shows the lowest bias at 0.482.

* Under Cauchit link (Methods 11-16): D R augmented M A I C shows the lowest bias at 0.350 and highest coverage at 0.403.

Under $n=1,000$ , all methods perform as expected in terms of bias. The naïve estimator appears biased in all four scenarios (but in KS3, the bias is very small), the DR estimators appear unbiased when either model is correct (Scenarios KS1–KS3), and the singly robust estimators avoid bias only when the corresponding propensity score or outcome model is correct. Unlike the three DR methods, the weighted G-computation estimators show bias in Scenario KS2 when the Cauchit link function is used in the outcome model (but not when the canonical logistic link is used). This suggests that the weighted G-computation estimators may be DR for the ATC when the canonical logit link is used, but not otherwise. The trends observed for the bias are similar under $n=200$ with the caveat that the (augmented and non-augmented) weighting estimators that appeared unbiased under $n=1,000$ exhibit some small-sample bias in the corresponding scenarios. This is particularly notable in Scenario KS2 and, to a lesser extent, in Scenario KS1, and is probably due to small ESSs after weighting.

Under correct specification of the outcome model (KS1 and KS3), G-computation is the most precise covariate-adjusted estimator, but the augmented estimators are almost as precise (e.g., compare the G-computation estimator which obtains ESE = 0.150 to the DR augmented MAIC estimator which obtains ESE = 0.170 for KS1 with $n=1,000$ ; see Table 1). Moreover, the augmented estimators generally produce precision gains versus their respective non-augmented weighting counterparts. When both the outcome model and the propensity score model are correctly specified (KS1), all augmented estimators have increased precision compared to the non-augmented weighting estimators based on modeling weights, but not necessarily against MAIC (any increase in precision for $n=1,000$ is modest). When only the propensity score model is correctly specified (KS2), outcome model misspecification does not induce any meaningful loss in precision for the augmented estimators compared to their non-augmented weighting counterparts.

There have been some concerns in the literature about DR augmented estimators amplifying bias and variance when misspecified weights are combined with a misspecified outcome model.Reference Kang and Schafer ⁸⁷ Such amplification is not apparent in our simulation study. In KS4, our proposed DR augmented MAIC estimator is the least biased of all estimators and the most precise of the augmented and non-augmented weighting estimators. Under $n=1,000$ , this advantage is modest when the logit link is used for the outcome model (bias of 0.482 vs. 0.512 for G-computation), but more pronounced when the Cauchit link is used (bias of 0.350 vs. 0.524 for G-computation).

In Section 3.3, we hypothesized that entropy balancing weights, like those employed by MAIC, can lead to more stable and precise ATC estimation than inverse odds modeling weights. This appears to be confirmed for the non-augmented estimators in our simulation study; MAIC exhibits greater precision than the approaches using (normalized or non-normalized) IOW modeling weights in all simulation scenarios. Additionally, the precision gains have been inherited by the augmented approaches. For the methods highlighted in Section 3.5, estimators using MAIC weights display enhanced precision compared to those using IOW modeling weights in all simulation scenarios, while producing similar levels of bias, even lower bias under dual model misspecification.

Assuming unbiasedness, interval estimation is appropriate if the coverage is approximately equal to 0.95; poor coverage can arise due to bias or to inadequate variance/interval estimation. Coverage is generally close to 0.95 for all covariate adjustment methods in the simulation scenarios under which they are unbiased, which suggests that our proposed non-parametric bootstrap approach for variance estimation is adequate. In the cases in which the DR methods are unbiased (i.e., KS1–KS3), coverage rates are between 0.935 and 0.947 across $n=1,000$ scenarios. Interestingly, coverage rates seem to increase for these methods under $n=200$ despite the small-sample bias, lying between 0.942 and 0.965. Note that, due to computational limitations, the non-parametric bootstrap approach in the simulation study was conducted with only $B=200$ resamples. This may have impacted the observed coverage rates and we suspect that coverage might be more appropriate when using a larger number of resamples.

While some covariate adjustment methods display bias-induced undercoverage in the scenarios under which they are biased (e.g., MAIC in KS3 under $n=1,000$ or all estimators in KS4), they may also display adequate coverage (e.g., normalized IOW with $n=200$ ) because of excessively large standard errors, probably due to low ESSs after weighting. As observed for KS4, bias-induced undercoverage tends to worsen with higher sample sizes, as interval estimates around the wrong target value become narrower. The naïve estimator displays discernible undercoverage in KS1, KS2, and KS4 (particularly under $n=1,000$ ), not only due to bias but also due to overprecise standard errors that do not account for covariate differences.