Average treatment effects on binary outcomes with stochastic covariates

Authors: Christoph Kiefer, Marcella L. Woud, Simon E. Blackwell, Axel Mayer

Abstract

When evaluating the effect of psychological treatments on a dichotomous outcome variable in a randomized controlled trial (RCT), covariate adjustment using logistic regression models is often applied. In the presence of covariates, average marginal effects (AMEs) are often preferred over odds ratios, as AMEs yield a clearer substantive and causal interpretation. However, standard error computation of AMEs neglects sampling‐based uncertainty (i.e., covariate values are assumed to be fixed over repeated sampling), which leads to underestimation of AME standard errors in other generalized linear models (e.g., Poisson regression). In this paper, we present and compare approaches allowing for stochastic (i.e., randomly sampled) covariates in models for binary outcomes. In a simulation study, we investigated the quality of the AME and stochastic‐covariate approaches focusing on statistical inference in finite samples. Our results indicate that the fixed‐covariate approach provides reliable results only if there is no heterogeneity in interindividual treatment effects (i.e., presence of treatment–covariate interactions), while the stochastic‐covariate approaches are preferable in all other simulated conditions. We provide an illustrative example from clinical psychology investigating the effect of a cognitive bias modification training on post‐traumatic stress disorder while accounting for patients' anxiety using an RCT.

INTRODUCTION

In (clinical) psychology and social sciences, the effects of an intervention, prevention, or treatment on a dichotomous outcome variable are often investigated using randomized controlled trials (RCTs). Examples include the effect of psychological interventions in prison on criminal recidivism (Beaudry et al., 2021), the effect of a drinking prevention strategy for college students on abstinence and occurrence of heavy drinking episodes (Larimer et al., 2007), the effect of interventions on smoking cessation (Dijkstra et al., 1998; Osch et al., 2008), and the effect of cognitive behavioral therapy on panic disorder with agoraphobia (Gloster et al., 2011).

Covariate adjustment is common in many of these RCTs for dichotomous outcomes, typically based on logistic regression models. The reasons for covariate adjustment are it can help to examine moderating effects of covariates or heterogeneity of treatment effects conditional on the covariate (or even on an individual level; see, for example, Wester et al., 2022) or it can be used to obtain effect estimates taking into account a predictive covariate, which can increase the power and thus reduce sample size requirements (Hernández et al., 2004; Moore & van der Laan, 2009). While covariate adjustment in nonlinear models (including logistic regression models) can result in higher efficiency, it is often cautioned that it also introduces bias into the effect estimates in small samples and should be avoided in situations where precision is most important (Imbens & Rubin, 2015; Robinson & Jewell, 1991). Recently, Negi and Wooldridge (2021) put forward that covariate adjustment in logistic regression models can indeed improve both efficiency and precision if the model is correctly specified. In addition, simulation studies have repeatedly shown consistency of treatment effects even from slightly misspecified logistic regression models (e.g., Negi & Wooldridge, 2021; Rosenblum & van der Laan, 2010). In this context, Negi and Wooldridge (2021) recommend always including treatment–covariate interactions unless one is sure that the model is correctly specified without them.

Logistic regression models can be considered as a part of the family of generalized linear models (McCullagh & Nelder, 1998; Nelder & Wedderburn, 1972) with a logistic link function and binomial distributed random component. They can be used to model the probability of the outcome occurring (usually coded as 1) conditional on a treatment, covariates, and possible covariate–treatment interactions. In logistic regression models, it is common to inspect treatment effects defined as risk ratios, odds ratios, or log odds ratios. Recent research suggests that these ratio effects are often misinterpreted in applied research (Niu, 2020) and it is often cautioned (Greenland et al., 1999; Hanmer & Ozan Kalkan, 2013; Mood, 2010) that these ratio effects might not even have a causal interpretation, especially when adjusting for covariates in an RCT.

Average marginal effects (AMEs) are suggested as an alternative to the odds ratio (Greenland et al., 1999; Hanmer & Ozan Kalkan, 2013; Mood, 2010; Norton & Dowd, 2018). A marginal effect is the difference between the conditional probability of the outcome given treatment and given control for a given observation. The mean of such marginal effects is called the AME and can easily be interpreted at the probability scale of the dichotomous outcome. For example, an AME of 5% would indicate that the probability of the outcome occurring increased by about 5% on average. Functions for estimating AMEs from logistic regression models are implemented in several statistical software packages, for example, the margins command in Stata (Williams, 2012) or the margins package in R (Leeper et al., 2021).

It is important to note that the estimation of AMEs does not account for sampling‐based uncertainty (Abadie et al., 2020). Sampling‐based uncertainty refers to sources of variability that are due to the sampling process. For example, if we want to gather a random sample of 100 participants for our RCT we could either predetermine a gender proportion and use stratified sampling (e.g., exactly 50 female participants) or we could randomly sample persons. In the latter case, the gender proportions are stochastic and may vary from sample to sample. While Wooldridge (2010) stated that accounting for sampling‐based uncertainty might be technically correct, but that ‘the adjustment may have a small effect’ (p. 919), Mayer and Thoemmes (2019) showed that properly accounting for this sampling‐based uncertainty in group sizes of categorical covariates can improve statistical inferences on average treatment effects. The improvement of statistical inferences can also be observed when treating continuous covariates as stochastic, that is, the observed values might change from sample to sample. If the sampling‐based uncertainty is neglected, standard errors tend to be underestimated, coverage of confidence intervals (CIs) can drop below the nominal level, and Type I errors can be inflated, especially if a strong treatment–covariate interaction (and, hence, heterogeneous treatment effects) is present. This has been shown in linear (Chen, 2006; Liu et al., 2017) and Poisson regression models (Kiefer & Mayer, 2019). However, a thorough investigation of this phenomenon for stochastic covariates in logistic regression models is lacking.

In this paper, we present and examine different approaches for estimating average treatment effects accounting for stochastic covariates using logistic regression models. First, we provide definitions of the average (causal) treatment effect on a dichotomous outcome as a difference in conditional probabilities in an RCT. In addition, we show that the variability of conditional effects is limited in this scenario. Second, we explain how these conditional probabilities can be estimated with and without accounting for a covariate, that is, we briefly recapitulate the statistics of contingency tables and logistic regression models. Third, we introduce three estimation approaches for the average treatment effect based on the estimated conditional probabilities. These are a simple difference‐in‐means estimator, approaches relying on the sample average of conditional effects (including the AME, but also a variant treating covariates as stochastic), and we propose a new moment‐based approach which accounts for sampling‐based uncertainty by definition. Fourth, we present a simulation study comparing the proposed estimators with a focus on differences between approaches with stochastic and fixed covariate. Fifth, we provide an empirical illustration of the estimators using data from clinical psychology (Woud et al., 2021) investigating the effect of cognitive bias modification training on symptoms of post‐traumatic stress disorder (PTSD) while accounting for patients' anxiety. Finally, we discuss the findings, implications, and limitations of this study.

DEFINITION AND TERMINOLOGY OF TREATMENT EFFECTS

Throughout this paper, we focus on the average effect of a randomized treatment on a dichotomous outcome variable. Our notation refers to the stochastic theory of causal effects (Steyer et al., 2022), which is similar to the Neyman–Rubin causal model (Rubin, 2005), but with a stronger emphasis on definitions and notation based on probability theory.

We consider a scenario with a binary outcome variable Y, a binary and randomized treatment X (with levels X=0 for control and X=1 for treatment) and a single continuous covariate Z. When considering a randomized experiment, the causal average treatment effect (ATE) of the treatment (X = 1) compared to the control group (X = 0) can be computed in different ways, two of which (1)ATE=P(Y = 1|X = 1)−P(Y = 1|X = 0) (2)=E[P(Y = 1|X = 1,Z)−P(Y = 1|X = 0,Z)⏟=CE(Z)].Equation (1) shows a simple difference‐in‐means formulation, where we are not conditioning on any covariates. For example, an ATE of .05 could reflect a 5% difference between a probability of the outcome of 50% under control and 55% under treatment, but also between 80% and 85%. When accounting for a (continuous) covariate Z,1the ATE can be computed as expectation of the Z‐conditional effects, as is shown in Equation (2). That is, the effect function CE(Z) represents treatment effects given specific values of Z. We will show below that this technique is used in the marginal effects approach. Note that Equations (1) and (2) hold because we assume randomized assignment to the treatment. Without randomized assignment, the right‐hand sides of these equations can be causally biased and would represent a prima facie effect (for more details, see, for example, Steyer et al., 2000).

Note that in clinical research the ATE is referred to as absolute risk reduction and a related measure is sometimes used, namely, the number needed to treat (NNT; Hutton, 2000). The NNT reflects the number of persons needed to be treated – on average – to yield one person more to benefit from the treatment compared to the control group. It can be computed as the inverse of the ATE, that is, NNT=1/ATE. For example, if we have an ATE of .2 (i.e., the probability of success is 20% higher in treatment than in control), the corresponding NNT is 5, meaning that we need to treat five persons, to get – on average – one beneficial outcome more than in the control group.

Both ways of computing the ATE – as presented in Equations (1) and (2) – are equivalent, but the two equations point to different possible ways to statistically model the dependencies among the outcome Y, treatment X, and possibly covariate Z and to estimate the ATE. This might lead to estimators with systematically different results for statistical inference. For example, the simple difference‐in‐means from Equation (1) does not contain any information about the variability of the individual treatment effects. That is, if a treatment works differently for different people, we have to condition on covariates in order to explore this variability, as is done via the effect function CE(Z) in Equation (2). This additional information can improve precision and efficacy in estimating the ATE, but the actual improvement depends on a number of factors, such as correct specification of the model and balance of the treatment groups (Negi & Wooldridge, 2021).

One factor, which is also important for statistical inferences on ATEs, is heterogeneity of the conditional treatment effects. If treatment effects conditional on a covariate vary a lot, accounting for this covariate will increase the efficacy of treatment effect point estimation (Hernández et al., 2004; Negi & Wooldridge, 2021), but it has also been shown to have an effect of coverage rates of CIs and the empirical detection rate (i.e., Type I and II errors; for example, for Poisson regressions, see Kiefer & Mayer, 2019). Thus, variability of conditional effects is an important concept that we will deepen in the next section.

Variance of conditional treatment effects

Obviously, the ATE for binary outcomes is always bounded between −1 and 1 (i.e., ATE∈[−1,1]). Consequently, the variance of conditional treatment effects is bounded for dichotomous outcomes, too. That is, regardless of the predictive power of the covariate, the maximum variance of the conditional effects for dichotomous outcomes is limited due to the binary nature of the outcome. Generally, the variability of the conditional treatment effects is given by the variance of the conditional effect function, Var[CE(Z)]=Var[P(Y = 1|X = 1,Z)−P(Y = 1|X = 0,Z)],which is generally bounded between 0 and 1 (i.e., Var[CE(Z)]∈[0,1]). The supremum (i.e., the lowest upper bound) of the variance of conditional effects given a specific ATE is(3)sup{Var[CE(Z)]:ATE∈[−1,1]}=1−ATE2.This means that the greater the absolute value of the ATE is, the less variance of conditional effects is possible. Conversely, an ATE of zero will yield the highest possible variance of conditional effects. A derivation of the supremum is given in Appendix A. This formula can help to decide what amount of variance is to be considered ‘small’ or ‘large’ in a certain scenario. For example, an actual effect variance of .4 might be ‘large’ if the maximum variance is also .4, but ‘small’ if the maximum variance is 1.0. In the simulation study below, we show that the proportion of effect variance relative to the maximum possible variance is an important measure in deciding whether effect variability affects point and standard error estimation.

ESTIMATION OF CONDITIONAL PROBABILITIES

In the previous section we provided nonparametric effect definitions based on the conditional probability P(Y = 1|X) or P(Y = 1|X,Z), respectively. In an applied scenario, we need a statistical model for estimating these conditional probabilities, for example, based on a contingency table or using a logistic regression model, respectively.

Contingency table

Given a sample of N independent and identically distributed (i.i.d.) observations of a binary outcome variable Y and a binary treatment variable X, we can compute the absolute frequencies Nyx of observations of each of the four possible combinations of X = x and Y = y:N00:=∑i=1N(1−Yi)·(1−Xi),N10:=∑i=1NYi·(1−Xi),N01:=∑i=1N(1−Yi)·Xi,N11:=∑i=1NYi·Xi.These frequencies are typically illustrated with a contingency table as in Table 1.

The marginal frequencies Ny· for Y = y and N·x for X = x are defined as column or row sums, respectively. Based on these absolute cell and marginal frequencies, we can estimate the conditional probabilities used in Equation (1) as relative (4)P^(Y=1|X=0)=N10N·0,P^(Y=1|X=1)=N11N·1.This procedure is equivalent to computing the group‐specific means of Y for the treatment and control group.

Before we move on to estimation and statistical inference for the ATE from these quantities, we consider how to model the conditional probabilities given in Equation (2) involving a continuous covariate.

Logistic regression model

In addition to Y and X, we now also consider an i.i.d. sampled continuous covariate Z. The conditional probability of Yi=1 of observation i, given the treatment and the covariate, is often modelled with a logistic regression model with treatment variable, covariate, and a treatment–covariate interaction as πi:=P(Yi=1|Xi,Zi)=expγ00+γ10·Xi+γ01·Zi+γ11·Xi·Zi1+expγ00+γ10·Xi+γ01·Zi+γ11·Xi·Ziwith the vector of parameters γ=(γ00,γ10,γ01,γ11).

Commonly, the regression coefficients γ are estimated using a maximum likelihood approach based on the generalized linear model framework (McCullagh & Nelder, 1998; Nelder & Wedderburn, 1972).2The log‐likelihood function for the logistic regression model islogL(γ)=∑i=1NYi·log(πi)+(1−Yi)·log(1−πi)which has to be solved iteratively. For a semiparametric estimation alternative, see, for example, Basu and Rathouz (2005).

The standard errors can be obtained via the Hessian matrix HLL of the log‐likelihood function evaluated at the maximum likelihood estimate, that is,HLL(γ^)=−∑i=1Nπ^i·(1−π^i)·xixiTwhere xi=(1,Xi,Zi,Xi·Zi)T (Miller, 2021), and the covariance matrix of the estimator is given byVar(γ^)=1N−1N·HLL(γ^)−1where the square root of the diagonal elements of Var(γ^) gives the respective standard errors. Negi and Wooldridge (2021) recommend using a robust estimator for the covariance matrix if it is used for treatment effect estimation – for example, a sandwich estimator as described by Stefanski and Boos (2002) which is implemented in many statistical software packages (e.g., the sandwich package in R by Zeileis, 2006).

Notice that the direct interpretation of the regression coefficients γ can be challenging in logistic regressions in presence of covariates and treatment–covariate interactions (Mood, 2010). Thus, it is often recommended to subsequently estimate the ATE (e.g., Hanmer & Ozan Kalkan, 2013) In the following section, we show how the ATE can be estimated from both the contingency table and the logistic regression model.

ESTIMATES AND STANDARD ERRORS OF TREATMENT EFFECTS

In this section, we present the estimation of and statistical inference on the ATE based on the estimated conditional probabilities from the previous section.

Simple difference‐in‐means

The simple difference‐in‐means estimator is fairly simple to obtain by plugging the observed frequencies from Equation (4) into the effect formula from Equation (1):ATE^SDM=P^(Y=1|X=1)−P^(Y=1|X=0)=N11N·1−N10N·0. ATE^SDM is a simple and consistent estimator of the ATE. The corresponding standard error isSESDM=N11·N01N·13+N10·N00N·03(see, for example, Fleiss et al., 2003, p. 60) and can, for example, be used to construct a 95% CI around ATE^SDM. In addition, one can test the hypothesis that both conditional probabilities are equal (which corresponds to an ATE of zero) using a two‐proportions test as described by Fleiss et al. (2003, p. 54). For an overview and comparison of further approaches to testing this hypothesis, see Newcombe (1998).

However, effect estimation is more challenging for the ATE when a continuous covariate within a logistic regression model is involved. In the following, we present two general approaches to estimating the ATE as the expectation of the conditional effect function CE(Z), namely, by means of computing a sample average or by integration over the covariate's domain.

Sample average over conditional effects

One possible way to translate the expectation over conditional effects from Equation (2) into a statistical procedure is by taking the sample average of the conditional effects (SACE) as an estimator for the ATE:(5) ATE^SACE is a consistent estimator of the ATE. The derivation of the estimator for this case, but also for cases with multiple covariates, can be found, for example, in Wooldridge (2010, Ch. 21).

For ATE^SACE, the covariate Z is still treated as an i.i.d. sampled random variable, which we call a stochastic covariate. The standard error for ATE^SACE is given by(6)SESACE=∇ATE^SACEγ^·Var(γ^)·∇ATE^SACEγ^T+1NVar^[CE(Zi)],where ∇ denotes the gradient of a function. This standard error formula has been derived and proven, for example, by Bartlett (2018), Basu and Rathouz (2005), Terza (2016), and Wooldridge (2010).

However, there exists a simplified version of ATE^SACE which is predominantly used in the econometrics and biostatistics literature and is more commonly implemented in statistical software – e.g., in the margins package in R (Leeper et al., 2021) or the margins command in Stata (Williams, 2012) – than ATE^SACE. In this approach, the AME is used as estimator for the ATE. The AME estimator is based on the estimated regression coefficients γ^ and observed values zi of Z in a specific ATE^AMEγ^=1N·∑i=1Nexpγ^00+γ^10+γ^01·zi+γ^11·zi1+expγ^00+γ^10+γ^01·zi+γ^11·zi−expγ^00+γ^01·zi1+expγ^00+γ^01·zi.This estimator looks very similar to the one given in Equation (5) and, in fact, the point estimates from both estimators are identical. Thus, the AME is also a consistent estimator of the ATE (Greene, 2012). However, there are two important differences between ATE^SACE and ATE^AME. First, the lower‐case zi emphasizes that only a certain set of observed covariate values is considered and not a set of random variables. The covariate is treated as fixed by design. Therefore, the AME is not meant to generalize beyond the sample at hand. Second, it is possible to derive the standard error for ATE^AME by simply using the delta method (for an introduction, see Raykov & Marcoulides, 2004). However, this approach neglects sampling‐based uncertainty, because the observed zi are not considered to be stochastic (i.e., as randomly sampled):SEAME=∇ATE^AMEγ^·Var(γ^)·∇ATE^AMEγ^T.The terms under the root are identical to the first part in the standard error formula for the SACE given in Equation (6), but the second part is omitted here. Thus, one can immediately see that these standard errors provide different results if there is considerable variance in the conditional effects (i.e., if the treatment works differently for different persons). In these cases, neglecting sampling‐based uncertainty will yield underestimated standard errors for the AME and, in turn, results in flawed inferences, for instance, poor coverage rates and inflation of Type I errors (as previously shown for Poisson regressions by Kiefer & Mayer, 2019). Conversely, if conditional effects are homogeneous (i.e., there is no effect variance), both standard error formulas should provide similar estimates even if the covariate was randomly sampled.

Integral over conditional effects

An alternative to the above‐presented sample averages over the conditional effects, is to estimate the ATE via an integral using moments of the covariate Z. In linear regression, this procedure simplifies to computing the conditional effect at the expectation of Z (Liu et al., 2017; Mayer et al., 2016). In this case, it suffices to estimate the mean μZ of Z and to evaluate the respective conditional effect function at this value. This procedure has the additional benefit that sampling‐based uncertainty in the covariate can be accounted for by including the standard error of μ^Z in the computation.

In nonlinear regression models, such moment‐based approaches usually require a distributional assumption for the covariate and potentially more than one moment (e.g., expectation and variance; for a moment‐based approach for Poisson regression models, see Kiefer & Mayer, 2019, 2021b). In addition, there is not always a convenient analytical solution as in the linear and Poisson regression case. To our knowledge, no moment‐based approach for logistic regressions has been proposed. Such approaches have been found to outperform sample average‐based estimators under specific conditions (e.g., Kiefer & Mayer, 2019), because the distributional assumption provides additional information unless it is violated. Thus, we derive such a moment‐based estimator for logistic regression models in the following.

For the logistic regression model, the ATE from Equation (2) based on the i.i.d. sampled variables can be rewritten to emphasize the meaning of the unconditional expectation as an (7)where pZ(z) denotes the density or probability density function of the covariate Z. Note that this approach treats the covariate Z as stochastic by construction, as we are not looking at observed conditional effects, but consider all conditional effects on the domain of the covariate weighted by its density function. The moment‐based estimator for the ATE can then be written as(8)ATE^MOM(γ^,θ^)=∫zCE(z,γ^)·pZ(z,θ^)dz (9)≈∑i=1Mwi∗·CE(zi∗,γ^),where Equation (8) shows estimation of the ATE with the integral over the product of the conditional effect function and the density function of Z using the parameter estimates γ^ and θ^. This is a consistent estimator of the ATE (for a proof, see Appendix B). However, the computation of this estimator usually requires numerical integration as for most densities there is no closed‐form solution. For example, if Z is normally distributed with estimated parameters θ^=(μ^,σ^) the density is well known, but no analytical solution exists to compute the integral.

In Equation (9) a numerical approximation of ATE^MOM is shown using a finite sum over M integration points zi∗ and weights wi∗. There exist various techniques to compute the integration points and weights, with several variants of Gaussian quadrature (e.g., Gauss–Kronrod, Gauss–Hermite) being readily implemented in most statistical software packages, for instance, the integrate function in R. However, the approximation of this estimator shown in Equation (9) can introduce bias if the numerical integration technique is inadequate. In Appendix C, we provide examples for unidimensional integration techniques of the proposed moment‐based approaches and also discuss settings with multiple covariates, where multidimensional integration is required.

Note that the integration points zi∗ themselves are fixed values, derived by the integration technique used. Thus, ATE^MOM is computed as a function of estimated parameters only – similarly as in the AME estimator. Consequently, the standard error for the moment‐based ATE estimate can be derived using the Delta method and is(10)SEMOM=∇ATE^MOM(γ^,θ^)·Var(γ^,θ^)·∇ATE^MOM(γ^,θ^)T,but the sampling‐based uncertainty in Z is naturally included via the parameters θ^ and by integrating over all possible conditional effects given the distribution of Z (for a proof, see Appendix B).

SIMULATION STUDY

The aforementioned ATE estimators and their respective standard errors can all be shown to be consistent, that is, they are unbiased with sample sizes increasing to infinity. In psychological RCTs, sample sizes tend to be rather small. Thus, we conducted a simulation study to investigate the finite‐sample properties of the different ATE estimators and corresponding standard errors under various possible scenarios. We were especially interested in the statistical inference for ATE estimators when incorporating or neglecting sampling‐based uncertainty. Thus, a central design factor of our simulation study is the variance of conditional treatment effects as this term reflects the key difference between the standard error formulas of the AME and SACE. While effect heterogeneity might stem from the influence of multiple covariates in applied settings, we used only a single covariate accounting for the variance of conditional treatment effects for simplicity. As the focus of this study is on consequences of different variances of the conditional treatment effects, we would argue that it is (at least from a technical perspective) of minor importance whether the effect heterogeneity is generated from one (comprehensive) covariate or from multiple covariates.

As effect variance is actually limited by the ATE, as was shown in Equation (3), we simulated different proportions of variance in relation to the maximum possible variance. For example, for an ATE of 0, the maximum possible effect variance would be equal to 1. So, we simulated effects with 80% (i.e., .8) variance, 50% variance and so on. Thus, we do not investigate absolute variance of treatment effects, but in relation to variance possible given an ATE.

From the existing literature we derived two additional design factors. First, we vary the total sample size from a very small sample of N=25 to a very large sample (N=5000). While such large samples are uncommon for psychological RCTs, we try to investigate at what point consistency of the estimator starts to show. Second, we differentiate between balanced designs (i.e., 50% probability of being treated) and unbalanced designs (i.e., 33% probability of being treated). Negi and Wooldridge (2021) found an effect of different degrees of balance on the root mean squared error of effect estimates, that is, the ATE is estimated more inefficiently in unbalanced designs. Especially for small samples, we would expect unbalanced designs to have a negative effect on all estimators, as the precision and efficacy within the treatment group should be reduced. This should also affect statistical inferences. An overview of our simulation study design is given in Table 2. The design results in a total of 686 conditions and each condition was replicated R=5000 times.

We refrained from investigating these factors for different distributions of the covariate, as this would be expected to cause biased point estimates under extreme conditions, but no additional information regarding the statistical inferences would be generated (see, for example, Kiefer & Mayer, 2019, examining different distributions for the moment‐based approach for Poisson regression models). In addition, the distributional assumption of the moment‐based approach always coincides with the true distribution of the covariate. Thus, the findings from this simulation study reflect the statistical inferences when this assumption is correct. In applied settings, it should be expected that both point estimates and standard errors can be biased if the distributional assumption is violated.

For point estimation, we computed ATE^SDM, ATE^AME, and ATE^MOM. The estimator ATE^SACE is not shown separately, as it is computationally identical to ATE^AME. We computed the absolute and relative bias of the estimators for each condition. For standard error estimation, we computed SE^SDM, SE^AME, SE^SACE, and SE^MOM, and computed the coverage of the 95% CIs as well as the empirical detection rate (EDR; i.e., Type I error for ATEs equal to zero and power otherwise) for each condition.

The commented R code of the data‐generating process and the simulation study can be found on OSF (https://osf.io/gu37s). The simulation study was carried out using R (R Core Team, 2021) with the SimDesign package (Chalmers & Adkins, 2020). Note that the SimDesign package generates “working” replications of each condition, that is, simulated data sets leading to non‐convergent models or related issues are automatically excluded from the analyses.

Results

An overview of the median results over all approaches and criteria is provided in Table 3.

Bias and relative bias of ATE estimates

First, we looked at the bias B of each of the point estimators. Note that ATE^AME and ATE^SACE provide identical point estimates, so we do not report them separately here. In general, we can see that bias was largest in conditions with very small sample size (i.e., N=25), large absolute ATE (i.e., |ATE|≥.3), and high relative effect variance. However, the bias quickly vanishes with increasing sample size. This finding is in line with derivations of Negi and Wooldridge (2021), that is, the estimators are consistent given the model is correctly specified (as was the case in our simulation). Figure 1 provides an overview of the results.

The median bias of ATE^SDM was |BM|<.0001 with 95% of the values lying in the interval [−.008,.008]. Thus, ATE^SDM was unbiased in most scenarios. The most extreme values of bias for ATE^SDM, Bmin=−.021 and Bmax=.020, were found with very small sample size (N=25). The median bias of ATE^AME (and ATE^SACE) was also |BM|<.0001 with 95% of the values lying in the interval [−.004,.005]. Thus, ATE^AME and ATE^SACE were unbiased in most scenarios and slightly more precise than ATE^SDM. The most extreme values of bias for ATE^AME and ATE^SACE were Bmin=−.012 and Bmax=.011, again slightly better than the results of the SDM estimator and only observed for very small sample sizes. Finally, the median bias of ATE^MOM was also |BM|<.0001 with 95% of the values lying in the interval [−.004,.004]. Thus, ATE^MOM was also unbiased in most scenarios with precision comparable to ATE^AME. The most extreme values of bias for ATE^MOM were Bmin=−.010 and Bmax=.008, thus showing the smallest range of bias among all conditions.

We also examined the relative bias RB, that is, the bias in relation to the true ATE. Note that a positive relative bias means that the absolute value of the ATE is overestimated, regardless of the direction of the effect. In general, we found that relative bias was within a range of ±5% for at least 95% of all conditions, and within a range of ±10% for all conditions and estimators. As for bias, the relative bias diminishes quickly with increasing sample size, meaning that the highest and lowest values of relative bias were found for very small sample sizes (N=25). In contrast to the bias, the extreme values for relative bias were found for small absolute ATE (i.e., .1). Figure 2 provides an overview of the results.

The median relative bias of ATE^SDM was RBM=−.04% with 95% of the values lying in the interval [−4.13%,1.77%]. The most extreme values of relative bias for ATE^SDM were RBmin=−9.48% and RBmax=4.94%, respectively slightly over and under the traditional cutoff of ±5%, and these were only found with very small sample size (N=25). The median relative bias of ATE^AME and ATE^SACE was RBM=−.03% with 95% of the values lying in the interval [−3.23%,1.64%]. Thus, ATE^AME and ATE^SACE tended to be slightly more precise than ATE^SDM. The most extreme values of bias for ATE^AME and ATE^SACE were RBmin=−9.45% and RBmax=5.39%, and were similar to the results from the SDM estimator. Finally, the median relative bias of ATE^MOM was RBM=−.02% with 95% of the values lying in the interval [−2.55%,2.04%]. The most extreme values of bias for ATE^MOM were RBmin=−8.44% and RBmax=5.20%. Overall, these results are comparable to the results of the SDM and AME/SACE estimators.

Overall, these findings suggest that bias and relative bias of all three estimators is acceptable across all conditions. The minimum values of relative bias were observed in the same condition for all three estimators. Given the low value of ATE=.1 in this condition, it seems plausible that this extreme value is due to chance, especially as such large relative bias was not systematically observed under similar conditions.

Coverage of 95% confidence intervals

In a next step, we computed the (symmetric) 95% CIs based on each ATE estimator and the respective standard errors and examined the coverage rates (i.e., how often the CIs actually include the true ATE at a nominal level of 95%). In general, we found that coverage rates C were acceptable for most conditions in the SDM, SACE, and MOM estimators, except for the combination of very small sample size (N=25) and an unbalanced design. In contrast, the AME estimator had too low coverage rates under all conditions. Figure 3 provides an overview of the results.

The median coverage of the SDM estimator was CM=.943 with 95% of the values lying in the interval [.898,.954], which is acceptably close to the nominal level of 95%. The most extreme values of coverage for the SDM were Cmin=.880 and Cmax=.959. The lowest value was slightly under a traditional cutoff of 90%, but this value was only found with very small sample size (N=25) and an unbalanced design. The median coverage of the AME estimator was CM=.821 with 95% of the values lying in the interval [.617,.921] which is substantially below the nominal level of 95%. The most extreme values of coverage for the AME were Cmin=.583 and Cmax=.928. Thus, the AME did not meet the nominal level of 95% in any scenario, but produced CIs that were too narrow under all conditions. The median coverage of the SACE was CM=.946 with 95% of the values lying in the interval [.902,.961], which is acceptably close to the nominal level of 95%. The most extreme values of coverage for the SACE were Cmin=.866 and Cmax=.968. Again, the lowest coverage rate was found in conditions with unbalanced design and a very small sample size (N=25). Finally, the median coverage of the MOM estimator was CM=.947 with 95% of the values lying in the interval [.914,.976]. The most extreme values of coverage for the MOM estimator were Cmin=.884 and Cmax=.989. As for the SDM and SACE estimators, the lowest coverage rate was found in conditions with unbalanced design and a very small sample size (N=25). The highest coverage rate was found in conditions with balanced design, a very small sample size, and large relative effect variance.

Empirical detection rate

Finally, we looked at the empirical detection rate, which is the Type I error rate T at a nominal level of 5% for an ATE of zero and the power otherwise. In general, we found that Type I error rates were acceptable under most conditions for the SDM, SACE, and MOM estimators, with an exception of the MOM estimator in very small sample sizes (N=25) when large relative effect variance was present. The AME estimator showed an inflated Type I error rate under all conditions, with values up to almost 40%. Figure 4 provides an overview of the results.

The median Type I error rate of the SDM was TM=.052 with 95% of the values lying in the interval [.045,.052]. The most extreme values of Type I error rate for the SDM, Tmin=.042 and Tmax=.060. These results reflect an acceptable error Type I error rate under all simulated conditions. The median Type I error rate of the AME was TM=.166 with 95% of the values lying in the interval [.081,.378]. The most extreme values of Type I error rate for the AME estimator were Tmin=.075 and Tmax=.392. Thus, the AME did not meet the nominal level of 5% Type I error rate in any scenario, but produced inflated Type I error rates under all conditions. The median Type I error rate of the SACE was TM=.052 with 95% of the values lying in the interval [.036,.070]. The most extreme values of Type I error rate for the SACE were Tmin=.030 and Tmax=.085. These results reflect an acceptable error Type I error rate under all simulated conditions. Finally, the median Type I error rate of the MOM estimator was TM=.052 with 95% of the values lying in the interval [.023,.064]. The most extreme values of Type I error rate for the MOM estimator were Tmin=.010 and Tmax=.070. While these Type I error rates are acceptable in most scenarios, the MOM showed a tendency for deflated Type I error rates in scenarios with very small sample sizes (N=75), balanced design, and large relative effect variance.

Our simulation study yielded varying results regarding the power of the different approaches. Broadly speaking, the SDM estimator showed substantially higher power than the SACE and MOM in scenarios with small sample sizes (N≤100), large ATE (|ATE|=.5), and 50% or more relative effect variance. For smaller ATEs, large relative effect variance, and sample sizes ≥250, the MOM showed higher power than the SDM and SACE. However, in most scenarios the power of SACE and MOM was similar to that for the SDM, meaning the increase in power due to accounting for a covariate was moderate at best. It is noteworthy that power was generally higher in conditions with balanced designs compared to their unbalanced counterparts. That is, the precision gained through a larger control group did not counterbalance the loss of precision for estimating the parameters in a smaller treatment group. We do not consider the AME here, as it did not meet the nominal level of 5% Type I error rate under any conditions and, therefore, the empirical detection rate cannot be interpreted as power.

EMPIRICAL ILLUSTRATION

In this section, we give an empirical illustration of how the above‐mentioned ATE estimators can be applied to real data from psychotherapy research using data from Woud et al. (2021), who examined the effects of cognitive bias modification (CBM) training on PTSD. They used an RCT and compared CBM to sham training (i.e., a control group). A total of N=80 participants were randomized, with N=65 providing the outcome data used in our analyses.

In our analysis, we investigated the ATE of the CBM training (X = 1) on a categorical PTSD outcome (estimated PTSD diagnosis) at the six‐week follow‐up assessment. While the ATE is usually expected to be larger at the immediate post‐training assessment, the effect variability is expected to be larger with later assessments due to differences between patients in the extent to which they retain and implement the learning from the intervention. As diagnostic interviews were only administered pre‐treatment, we operationalized the outcome variable Y by means of reaching a certain cut‐off score on the German version of the PTSD checklist for DSM‐5 (PCL‐5; Krüger‐Gottschalk et al., 2017). We applied a cut‐off score of ≥33 as recommended by the test developers (see Krüger‐Gottschalk et al., 2017). That is, Y = 1 indicates a high probability for a PTSD diagnosis and Y = 0 indicates the opposite.

We used the patients' score on the Beck Anxiety Inventory (BAI; Beck et al., 2012) measured at pre‐training as our covariate. While this covariate was chosen largely for demonstration purposes, pre‐training anxiety levels could plausibly be expected both to prognostic of treatment outcomes and also interact with treatment condition. High scores on the BAI reflect general (i.e., not PTSD‐specific) anxiety severity, which could interfere with patients' ability to engage with and benefit from treatment as usual, leading to worse outcomes; patients with higher levels of anxiety at baseline may therefore particularly benefit from an adjunctive treatment (such as CBM).

We estimated the ATE using all four aforementioned SDM, SACE, AME, and MOM. For the covariate‐based estimators, we estimated a logistic regression model with robust standard errors, to account for potential misspecification as recommended by Negi and Wooldridge (2021). The commented R code for our analyses can be found on OSF (https://osf.io/gu37s).

Results

First, we estimated ATE^SDM based on the contingency table given in Table 4. The estimate for ATE^SDM was −.164 (Δ=−.376), which means that the probability for a PTSD diagnosis six weeks after treatment was about 16.4% lower in the CBM training group compared to the control group. That is, ATE^SDM reflects the difference between the probability of a PTSD diagnosis under control (75.8%) and under training (59.4%). The estimated ATE also reflects a number needed to treat (NNT) of about 6.1, meaning that 6.1 more patients would have to be treated (than non‐treated), to achieve one diagnosed person less – on average. A summary of the results is given in Table 5.

Second, we estimated a logistic regression model with treatment variable (γ^10=4.034, p=.074), covariate (γ^01=.155, p=.023), and treatment–covariate interaction (γ^11=−.164, p=.026). Based on the logistic regression model, we estimated both ATE^AME (−.144, Δ=−.330; same values for ATE^SACE), and ATE^MOM (−.139, Δ=−.318). Both covariate‐based ATE estimates were slightly smaller than the SDM estimate, which is also reflected in the corresponding effect sizes.

Third, we investigated the standard errors, 95% CIs, and p‐values of the respective ATE estimates. While all four p‐values would result in keeping the null hypothesis of ATE of zero, the approaches differ with regard to their statistical confidence of doing so. For AME and SACE, the standard errors are lower (SE^AME=.107; SE^SACE=.112) than for the SDM (SE^SDM=.114) and, therefore, the range of the corresponding CIs is smaller, reflecting a higher confidence in the statistical inferences. However, based on our simulation study we would conclude that the AME is overconfident in this case by simply neglecting sampling‐based uncertainty in the covariate. In contrast, the SACE takes this uncertainty into account and still improves beyond the SDM, but in a more moderate way.

In this application, the moment‐based approach results in both the highest standard error (SE^MOM=.123) and the widest 95% CI of all four approaches. Thus, we additionally used a Shapiro–Wilk test to examine whether the normality assumption for the covariate was reasonable and got a non‐significant result (W=.984, p=.574).

In sum, we compared the simple estimate of the ATE (i.e., not including any covariates, just comparing group means, called SDM) to several different ways (i.e., AME, SACE, MOM) of calculating the ATE based on inclusion of a covariate that is prognostic of treatment outcomes and interacts with the treatment condition. The comparison illustrates two important aspects. First, with regard to point estimates, all four estimators yield very similar results, which is to be expected in an RCT. Second, adjusting for a covariate and ignoring sampling‐based uncertainty (i.e., the AME) seemingly provides favourable standard errors, p‐values, and CIs, but these are due to an underestimation of uncertainty and would lead to inflated Type I error rates and undercoverage of CIs. Accounting for sampling‐based uncertainty in the covariate leads to statistical inferences closer to the simple ATE estimate (without covariates). While the overall implications and inferences would be similar for all four estimators in this application, the example also illustrates that substantial differences are possible even in a very simple scenario like this.

DISCUSSION

In this paper, we provided a thorough investigation of the adjustment for stochastic covariates in estimation of average treatment effects in logistic regression models in psychological RCTs. We presented four adjusted and unadjusted estimators of the ATE: a simple difference‐in‐means estimator, a sample average of conditional treatment effects treating the covariate as stochastic (i.e., SACE) and one treating it as fixed (i.e., AME), and a newly developed moment‐based approach, which also treats the covariate as stochastic. We discussed the conditions under which these estimators differ with regard to their standard errors and examined the finite‐sample properties of all four estimators in a simulation study. Finally, we offered an empirical example from clinical psychology illustrating the different results provided by the four estimators.

The most important finding from our simulation study is that statistical inferences from the AME estimator (i.e., treating the covariate as fixed) were problematic under almost all conditions. This is in contrast to the statement of Wooldridge (2010) that omitting the part accounting for sampling‐based uncertainty in the standard error formula might only have a small effect. We found that for the AME actual coverage rates did not meet the nominal level of 95%, and the actual Type I error rates did not meet the nominal level of 5% in any condition. These rates were close to their nominal levels when the conditional treatment effects were close to homogeneous. However, homogeneous treatment effects are only possible in logistic regression models if the covariate is not predictive of the outcome at all, and in this case covariate adjustment would only introduce noise into the analysis (Negi & Wooldridge, 2021).

Another important finding is that covariate adjustment did not automatically lead to increased power of an estimator, even though the logistic regression model was correctly specified in our simulation. The increased empirical detection rates of the AME estimator must not be misinterpreted as power – otherwise, one would wrongly conclude that covariate adjustment tremendously improves the power. This is generally not the case. We found that covariate adjustment (with SACE and MOM) was only beneficial in larger samples, when estimating a small to moderate ATE with large variance of conditional effects. In scenarios with larger ATEs, the SDM showed higher power than the covariate‐adjusted approaches. One possible reason for this finding is that the maximum likelihood estimation of the regression parameters might come with much estimation uncertainty in the presence of strong interactions in small samples, which in turn reduces the power of the effect estimate. However, in larger samples the power of the SDM is already quite high, so there is not much room for covariate adjustment to improve beyond. This finding is noteworthy especially because the SDM only relies on the assumption of i.i.d. observations, whereas the other estimators require additional assumptions, such as a correct model specification.

These findings emphasize the importance of accounting for sampling‐based uncertainty when estimating (average) treatment effects based on regression models. They are in line with findings from previous work, for example, on accounting for stochastic group sizes in an ANOVA‐like framework (Mayer & Thoemmes, 2019), stochastic covariates in linear regression models (Chen, 2006; Liu et al., 2017), and stochastic covariates in Poisson regression models (Kiefer & Mayer, 2019). We contribute to this line of research by illustrating the severity of neglecting sampling‐based uncertainty for logistic regression models.

In addition, we proposed a new moment‐based approach which also accounts for stochastic covariates. In a previous study, Kiefer and Mayer (2019) showed that a moment‐based approach can outperform the SACE in Poisson regression models with regard to bias of point estimates and accuracy of coverage rates in the presence of strong interaction effects. However, we could not find these advantages in logistic regression models. While the moment‐based approach yield slightly better coverage rates in very small samples than the SACE, both estimators yield similar results with regard to bias. One possible reason for these findings is that we only examined a normally distributed covariate, while Kiefer and Mayer (2019) found a Poisson‐distributed covariate as leading to huge performance differences between the two estimators. Nevertheless, the moment‐based approach can be a useful in scenarios where the SACE is not directly applicable, for example, with latent covariates (Kiefer & Mayer, 2021a).

We have shown that statistical inferences on the ATE can be difficult, when treatment effects are heterogeneous. In the same vein, it is noteworthy that the ATE itself might be substantively less interesting in cases with heterogeneous treatment effects, as the actual (conditional) treatment effects deviate from the ATE. In our illustrative example, CBM training was found to be more effective for persons with higher anxiety at baseline and less effective for persons with lower anxiety. This information can be useful for tailoring treatments to persons for whom it is beneficial. For another example, see the study on differential effects of a classroom intervention Flunger et al. (2019). Based on covariate selection, there also exist approaches to (approximately) estimate individual treatment effects (Mayer et al., 2020; Wester et al., 2022).

In sum, we showed that accounting for sampling‐based uncertainty is important when doing covariate adjustment for a binary outcome in a psychological RCT. Neglecting this uncertainty can lead to severely inflated Type I error rates and overestimation of power. Instead of the AME we recommend using the SACE, which provides identical point estimates. If treatment effects are heterogeneous, the statistical inferences provided by the SACE are more accurate than those provided by the AME, and if treatment effects are homogeneous, both approaches would provide identical results. The moment‐based approach can be a viable alternative to the SACE, but might be inaccurate in the presence of strong interactions in small samples.

AUTHOR CONTRIBUTIONS

Christoph Kiefer: conceptualization; methodology; software; writing – original draft. Marcella L. Woud: resources; writing – review and editing. Simon E. Blackwell: resources; writing – review and editing. Axel Mayer: conceptualization; writing – review and editing; supervision.