Locally Efficient Estimation of Marginal Treatment Effects When Outcomes Are Correlated: Is the Prize Worth the Chase?

1 Introduction

Semiparametric estimators are appealing because of their robustness to distributional assumptions and model misspecification. In the analysis of randomized trials, semiparametric theory has been used to develop estimators of treatment effects that improve efficiency of inferences by incorporating baseline covariates, where “baseline” describes data measured prior to randomization. In this paper, we present a semiparametric locally efficient estimator to improve efficiency of inferences in randomized experiments with correlated outcomes when baseline covariates are available. We begin with a review of current estimators for multivariate outcomes and then introduce our semiparametric locally efficient estimator.

Correlated outcomes are often observed in medical research studies, such as those that randomize clusters of subjects or that randomize individual subjects but collect repeated measures of the response. We denote the outcome for the ith independent randomized unit, i=1,…,m, in such studies by the ni-dimensional response vector Yi=(Yi1,Yi2,…,Yini)T, which may represent longitudinal measurements from a single subject or a set of responses from subjects within a cluster defined by a family, hospital, or class. Considering the substantial costs incurred by such studies, it is of interest to maximize efficiency in the estimation of treatment effects using all available data.

In general, studies collect data on i.i.d. observations Oi=(Yi,Ai,Xi), where Ai denotes a scalar treatment assignment to 1 of K possible treatments, and Xi is a matrix of baseline covariates. Throughout we allow ni to be fixed or random and assume ignorability when ni is random. Longitudinal data also include a time variable ti=(ti1,ti2,…,tini)T denoting time points at which outcomes are measured. As in the case of unit size ni, we allow ti to be either fixed or random but ignorable. When repeated measures are taken on the same subject, baseline covariates are measured at tij=0; thus Xij=Xi for all j=1,2,…,ni, resulting in a single level of baseline covariate information. Clustered data, however, may include pre-treatment covariates at the level of the group or the individual, creating two layers of auxiliary data. In the longitudinal context, we refer to the vector Yi as the subject, or independent unit and Yij as observation- or measurement-level data. For clustered data, we refer to Yi as cluster-level and Yij as individual-level observations.

Semiparametric estimation often involves specifying a restricted mean model. When estimating marginal treatment effects, a model for the expected outcomes given treatment assignment is usually assumed. Consequently, only data on the treatment and outcome are used in estimation. For example, in longitudinal studies, the marginal effect of treatment over time may be measured by assuming the restricted mean model

where f1(tij) is a function of ti. The main effect βA, which measures imbalance in E(Yij|Ai,tij) at baseline, is expected to be zero when randomization successfully balances covariate profiles across treatment arms. The post-baseline effect of treatment is measured by βA,t. Parameters βt and βA,t may be vector-valued, as the function describing the effect of time on expected outcomes may be of some polynomial form. Similarly, for clustered data, the semiparametric model

may be assumed, with treatment effects determined by inference on β1.

Estimating equations are determined by geometric arguments that distinguish parameters of interest, such as the treatment-outcome association (β) in the context of randomized studies, from other components needed to fully specify the data-generating distribution, which are represented by η. For parameters of interest, we aim to derive regular asymptotically linear (RAL) estimators, where an asymptotically linear estimator βˆ is one for which there exists a function ψ(Oi) such that

Regularity conditions ensure that variance bounds are well-defined and exclude superefficient estimators that have undesirable properties under local alternatives [1]. The function ψ(Oi) is called the influence function of βˆ and determines its limiting distribution. As eq. (3) suggests, any RAL estimator may be obtained by solving an influence function equation. To derive the class of estimating functions under an assumed model M, one first defines the nuisance scores ∂log(Lη)/∂η for the data-generating distribution Lη; one then determines the subspace defined by the closed linear span of all scores of smooth parametric submodels Lη in model M. This nuisance tangent space is denoted by Λnuis [2]. The orthogonal complement of Λnuis, Λnuis⊥ defines the set {ψ(h):h}, indexed by h, which contains the set of influence functions of all RAL estimators [2, 3]. For correlated outcomes, the geometric arguments of semiparametric theory may be viewed as a generalization of the quasilikelihood approach of Liang and Zeger [4] in deriving generalized estimating equations (GEE). We denote as M1 the set of distributions of Wi=(Yi,Ai) with known treatment process satisfying eq. (1). Under model M1, Λnuis1⊥=ψ(Wi;h,β)=h(Ai,ti){Yi−g(Ai,ti;β)} defines the estimating equations

for estimating the p-dimensional vector β. The index or weight h(Ai,ti) is a p×ni function of a random treatment variable Ai and time ti, and g(Ai,ti;β)={g(Ai,t0;β),g(Ai,t1;β),…,g(Ai,tni;β)}T. We use bold g(Ai,ti;β) to denote the vector-valued mean function and g(Ai,tij;β) to represent its scalar components.

A locally efficient estimator of a semiparametric model is defined as an estimator that achieves the semiparametric efficiency bound (minimum asymptotic variance among all RAL estimators) at a given submodel for the data-generating law, but remains consistent outside the data-generating submodel [2], provided that the marginal model is correct. More explicitly, semiparametric models parametrize specific components of a data-generating process and leave others unspecified. Estimation may require working models of unspecified components; a locally efficient estimator achieves the semiparametric efficiency bound when such working models are correctly specified, but also remains consistent when only the parametric component is correctly specified. A semiparametric locally efficient estimator is determined by finding the optimal estimating function, referred to as the efficient score, which for GEE requires finding the optimal h(⋅). When no baseline covariates are observed, Chamberlain [5] showed that the efficient score of β is obtained by setting h(Ai,ti)=DiTVi−1, where Vi is the ni×ni variance–covariance matrix of Yi, and Di=∂g(Ai,ti;β)∂β. The estimator remains consistent, however, when a working covariance other than the true covariance is substituted into the estimating equations, thereby demonstrating that consistency is achievable outside of the data-generating law.

For model M2, defined by observations Oi, marginal model (1), and known treatment process, the set of influence functions was derived by Robins et al. [6] and arises by augmenting the influence functions of β under model M1. Augmented estimators are constructed by subtracting the orthogonal projection of the standard estimating function onto the span of the scores of the treatment mechanism from the standard estimating function [6, 7]. For correlated outcomes, Λnuis2⊥={ψ(Oi,h,γ,β):h,γ}, and augmented GEE are

where for K-level treatment Ai, P(Ai=a)=πa. Fixing h(Ai,ti), the most efficient estimating function sets γk(Xi)=γkopt(Xi)=h(k,t)E(Y|Ai=k,Xi,t)−g(k,t;β) [6, 8–10]. The augmentation therefore involves estimation of the conditional mean outcome regression model E(Yi|Xi,Ai), which may be related the marginal parameter of interest in the binary treatment setting by β=EE(Yi|Xi,Ai=1)−E(Yi|Xi,Ai=0) under an identity link and β=logit[EE(Yi|Xi,Ai=1)]−logit[E{E(Yi|Xi,Ai=0)}] under the logit link. When baseline covariates are predictive of the outcome augmentation reduces variability in estimated treatment effects, irrespective of the outcome distribution. For the longitudinal marginal model (1), if outcomes Yij are restricted to post-baseline measurements, the baseline measurement Yi0 may be utilized as a baseline covariate and included in Xi. The βA,t term is then no longer required to assess a post-baseline effect of treatment and may be removed from the model, leaving βA to capture the marginal treatment effect. The interaction term βA,t may still be required for correct model specification even when the baseline outcome is included as a covariate if the treatment effect varies in time.

Semiparametric locally efficient estimators of parameters in restricted mean models of marginal treatment effects have been implemented for univariate data in the presence of baseline covariates by Robins [8], Bang and Robins [11], van der Laan and Rubin [12], Tsiatis et al. [13], Zhang et al. [10], and Moore and van der Laan [14, 15]. In these developments, the choice of h(⋅) has no impact on the resulting asymptotic variance and is therefore not considered for deriving efficient estimators. For a univariate outcome, the model gs(Ai;β)=E[Yi|Ai] defined by a unique parameter for each treatment level is saturated, and the choice of h(⋅) is inconsequential. When Yi is multivariate, gs(Ai;β)=E[Yij|Ai] is not saturated because a single parameter β is shared across components of the vector E[Yi|Ai]. As a result, Λnuis2⊥ provides a larger set of estimating functions than the orthogonal complement of the nuisance tangent space of the corresponding restricted mean model for a univariate outcome. Each element in Λnuis2⊥ is indexed by h(⋅); the choice of h(⋅) impacts efficiency, and the optimal h(⋅) must be found to achieve minimum variance. A similar discussion of h(⋅), model saturation, and estimating functions was included in Neugebauer and van der Laan [16].

The efficient score in model M2 does not generally have the same optimal index h(Ai,ti) as the efficient score in model M1. When incorporating auxiliary covariates in the estimation of marginal treatment effects via augmented GEE, the choice h(Ai,ti)=DiTVi−1, while resulting in a consistent estimator, is no longer optimal in model M2. The efficient score is determined by optimizing over all p×ni index functions h(Ai,ti) [6, 7, 9]. Robins [7] established general theory for deriving the efficient score of treatment effects in marginal structural models (MSMs) of time-dependent exposures, including the case of multivariate outcomes. Application of the Robins [7] theory to establish locally efficient estimators in specific settings, such as for randomized trials with correlated outcomes, requires further derivation. Additionally, the locally efficient estimators of Robins [7] were neither implemented nor evaluated for practical use. Models (1) and (2) may be viewed as examples of MSMs for a point exposure; the Robins [7] theory therefore equally applies. Although the efficient score may be obtained theoretically, it is often computationally intensive to calculate. Consequently, inefficient estimators are typically used. The suboptimal estimator based on augmenting GEE with h(Ai,ti)=DiTVi−1 was shown to improve efficiency by Zhang et al. [10] within the context of the linear mixed model and Stephens et al. [17] for general continuous and binary outcomes. In subsequent text, we refer to unaugmented GEE (4) under model M1 with the index function h(Ai,ti)=DiTVi−1 as standard GEE, and the suboptimal estimator obtained by augmenting standard GEE is referred to as simple augmented GEE. Here we show how to further improve on simple augmented GEE by deriving the corresponding semiparametric locally efficient estimator for model M2. We then evaluate the feasibility of achieving such improvement in practice.

The following section presents the efficient score and derives a locally efficient estimator of marginal treatment effects in randomized trials with correlated outcomes when auxiliary data are available as in model M2. We also discuss an implementation procedure detailing how to appropriately estimate each component of the efficient score. In Sections 3 and 4, we compare the derived semiparametric locally efficient estimator to standard and simple augmented GEE through a simulation study and application to the AIDS Clinical Trial Group Study 398, a randomized longitudinal HIV intervention trial.

2 Methods

3 Simulation study

Semiparametric locally efficient GEE for model M2 were compared to standard and simple augmented GEE through a simulation study. Simulations were completed for clustered data with continuous and binary outcomes and longitudinal data with continuous outcomes. Results are based on 1,000 Monte Carlo datasets.

3.1 Continuous outcomes

3.1.1 Clustered data

Data for m=500 clusters were generated, with ni = 2, 4, 6, 8, 10, 12 with equal probability for the first set of simulations and ni = 10, 20, 30, 40, 50 in the second set. Auxiliary covariates Xij1, Xij2, and Xij3 were each generated from a multivariate normal distribution with Var(Xij1)=2, Var(Xij2)=6, and Var(Xij3)=5. Correlation was induced among individual-level covariates within the same cluster by setting Cov(Xij1,Xij′1)=ςX1, Cov(Xij2,Xij′2)=ςX2, and Cov(Xij3,Xij′3)=1. Covariance terms ςX1 and ςX2 were varied from 0.5 to 2 and 1.5 to 6, respectively, to evaluate the effect of auxiliary covariate correlation on the performance of locally efficient augmented GEE. At ςX1=0.5 and ςX2=1.5, covariates were weakly correlated among individuals in the same cluster, while at ςX1=5 and ςX2,4=6, covariates were perfectly correlated, thereby becoming cluster-level. The exact values considered for ςX1 and ςX2 were (0.5, 1, 1.5, 2) and (1.5, 3, 4.5, 6), for simulation sets 1–4 at each set of cluster sizes. Within the jth individual in the ith cluster, auxiliary covariates were independent. The treatment variable Ai was drawn from the Bernoulli distribution with p = 1/2. Clustered responses were generated from the following model, with individual-level error terms εij∼N(0,40) and cluster-level effects bi∼N(0,σb2):Yij|Ai,Xij,bi=1.0+1.1Xij12+0.9Xij2+0.5Ai+bi+εij. The proportion of variability in Yij explained by auxiliary covariates Xij was held fixed at roughly 25%. Simulations were completed with σb2=0 and σb2=6, representing the case in which covariates account for all between-cluster heterogeneity (V(Y|X,A) independent) and the alternative of some intracluster correlation caused by an unmeasured variable (V(Y|X,A) exchangeable), respectively.

For each dataset, the marginal effect of treatment was estimated by fitting model (2) through standard, simple augmented, and locally efficient GEE for M2. The impact of misspecification on the locally efficient estimator and its efficiency relative to simple augmented and standard GEE was evaluated by fitting various models to estimate E(Y|X,A). The correct model for E(Y|X,A), denoted by “C” in tables and figures, was E(Yij|Xij,Ai)=ξ0+ξ1Xij12+ξ2Xij2+ξ3Ai, and two additional models: Model 1, “M1” = E(Yij|Xij,Ai)=ξ0+ξ1Xij1+ξ2Xij2+ξ3Ai, and Model 2, “M2” = E(Yij|Xij,Ai)=ξ0+ξ1Xij12+ξ2Xij2+ξ3Xij3+ξ4Ai. “Model 1” evaluated the impact of misspecifying the functional form of Xij1, while “Model 2” examined the effect of adding noise to the outcome regression model. All working covariance matrices were fit under exchangeable structure.

Efficiency comparisons relative to standard GEE are summarized in Figure 1(a) and 1(b), while the Monte Carlo Relative Efficiency (MCRE) of the locally efficient estimator for M2 to simple augmented GEE is given in Table 1. Small cluster figures are included in the supplementary material. Across all levels of correlation, augmented estimators resulted in increased efficiency compared to the unaugmented estimator (MCRE 1.25–11.6). For low correlation among Xij simple augmented and locally efficient augmented estimators performed similarly. Simple augmented GEE and locally efficient GEE for M2 also resulted in similar efficiency when the conditional mean model did not include the data-generating quadratic term Xij12, or the true conditional variance was exchangeable (MCRE locally efficient to simple augmented GEE 0.99–1.01). When correlation was increased among Xij within a cluster, the assumed conditional mean model included all important covariates in the correct functional form, and baseline covariates accounted for all within-subject correlation, locally efficient GEE for M2 gained in efficiency over the simple augmented GEE (MCRE locally efficient to simple augmented GEE 1.04–1.22). Increased covariance among auxiliary covariates also resulted in greater efficiency gains for any augmented GEE relative to the standard estimator. Trends were more pronounced for large average cluster size (average ni = 30 vs average ni = 7).

Figure 1

MCRE of locally efficient and simple augmented GEE relative to standard (unaugmented) GEE. Continuous clustered outcomes. Estimators corresponding to each curve are denoted by “estimator-outcome regression” using the abbreviations: Loc Eff – locally efficient, Simp – simple augmented, Std – standard, C – correct, M1–Model 1, M2–Model 2. All estimators use exchangeable working covariance for V(Y|A) and V{E(Y|X,A)}. The order of curves in the legend follows the order of curves in the figure, with sets of superimposed curves denoted by “()”, “[]”, or “{}”

Table 1

MCRE of locally efficient augmented GEE to suboptimal augmented GEE: continuous clustered outcomes

Cluster size = 2, 4, 6, 8, 10, 12
	Correlation among Xij
WMCov/OR	0.25	0.50	0.75	1.00
Exch/C	1.0115	1.0450	1.0907	1.1464
	1.0036	0.9991	1.0010	1.0085
Exch/M1	1.0062	1.0089	1.0064	1.0038
	1.0006	1.0008	1.0018	1.0019
Exch/M2	1.0114	1.0448	1.0905	1.1462
	1.0036	0.9990	1.0009	1.0083
Cluster size = 10, 20, 30, 40, 50
	Correlation amongXij
Cov/OR	0.25	0.50	0.75	1.00
ExchC	1.0356	1.1096	1.1563	1.2259
	1.0005	0.9999	1.0002	1.0011
ExchM1	1.0126	1.0081	1.0050	1.0032
	1.0000	1.0000	1.0001	1.0003
ExchM2	1.0352	1.1090	1.1556	1.2247
	1.0006	0.9998	1.0001	1.0009

^[1]

3.1.2 Longitudinal responses

For each Monte Carlo dataset, m = 500 longitudinal response vectors Yi were generated from the model Yij=1.5+1.1Xi12+0.9Xi2+1.0tij+1.0Ai+εij, where εij∼N(0,σε2), and Cov(εij,εij′) had an AR1 structure with correlation parameter α=0.1,0.3, or 0.5 for different sets of simulations. Covariates Xi1 and Xi2 were normally distributed with mean 0 and variance σX12 and σX22, respectively. Variance parameters σε2, σX12, and σX22 were varied so that baseline covariates accounted for 10–60% of the variability in Y|A in increments of 10%. Subjects were randomly assigned to treatment (Ai = 1) with probability 1/2. For each subject ti=(ti1=1,ti2=2,…,tini=ni), where ni varied from 1 to 8, as might be the case in a longitudinal study with staggered entry.

Standard GEE, simple augmented GEE, and locally efficient GEE for M2 were applied to each Monte Carlo dataset to estimate marginal treatment effects. All GEE were fit based on the marginal mean model E(Yij|Ai)=β0+β1Ai+β2tij with inferences on the treatment effect completed through β1. Standard and simple augmented GEE were applied to each Monte Carlo dataset with AR1, exchangeable, and true working covariance structures, with the true structure under the marginal model being a summation of AR1 and exchangeable matrices as described in Section 2. Locally efficient GEE for M2 were fit under the true covariance structure and a misspecified marginal AR1 working covariance. Baseline covariates were incorporated fitting several outcome regression models. We use “C” to denote the correct model E(Yij|Xij,Ai,tij)=ξ0+ξ1Ai+ξ2tij+ξ3Xi12+ξ4Xi2, which corresponds to the true data-generating mechanism; “M1” indicates the model E(Yij|Xi,Ai,tij)=ξ0+ξ1Ai+ξ2tij+ξ3Xi1+ξ4Xi2, omitting the exponent on Xi1; and “M2” is the model that includes a noisy covariate Xi3, such that E(Yij|Xi,Ai,tij)=ξ0+ξ1Ai+ξ2tij+ξ3Xi12+ξ4Xi2+ξ5Xi3.

Efficiency comparisons are summarized in Figure 2 and Table 2. Additional figures are given in the supplementary material. For well-specified variance components and conditional mean models, the locally efficient GEE for M2 was more efficient than the simple augmented GEE, with the difference in efficiency increasing with the percent variability explained by Xi (MCRE of locally efficient to simple augmented GEE 1.0–1.27). Similarly, all augmented estimators were more efficient than standard GEE, with efficiency gains from augmenting increasing with correlation in Y and X (MCRE of augmented GEE to standard GEE 1.36–5.28). For poorly specified conditional mean models, locally efficient GEE for M2 and simple augmented GEE were nearly equally efficient (MCRE of locally efficient to simple augmented GEE 0.97–1.0), but when the marginal variance was also misspecified locally efficient GEE were less efficient than simple augmented GEE (MCRE 0.88–0.99). This demonstrates that the locally efficient GEE for M2 is a bit more sensitive to working marginal covariance misspecification than simple augmented GEE. Among the simple augmented estimators, the estimator with the incorrect marginal AR1 working covariance resulted in the β1 estimate with the lowest variability. This illustrates an important distinction between locally efficient and suboptimal estimating functions. Considering estimators using a suboptimal index, misspecified models for parameters in the index may result in more efficient inferences than correctly specified models. For the locally efficient estimator, semiparametric asymptotic efficiency is achieved only in the absence of model misspecification of all parameters in hopt(⋅).

Figure 2

MCRE of locally efficient and simple augmented GEE relative to standard (unaugmented) GEE. Continuous longitudinal outcomes. Estimators corresponding to each curve are denoted by “estimator (marginal working covariance) outcome regression” using the abbreviations: Loc Eff – locally efficient, Simp – simple augmented, Std – standard, AR1–autoregressive(1) V(Y|A), true-exchangeable/AR1 for V{E(Y|X,A)} and V(Y|X,A); C – correct, M1–Model 1;α = 0.3. The order of curves in the legend follows the order of curves in the figure, with the set of superimposed curves denoted by “[]”

Table 2

MCRE of locally efficient augmented GEE to suboptimal augmented GEE: continuous longitudinal outcomes

	Correlation between Y and X
WMCov/OR	10	20	30	40	50	60
True/C	1.0281	1.0700	1.1168	1.1662	1.2175	1.2702
	1.0166	1.0425	1.0728	1.1055	1.1398	1.1752
	1.0090	1.0234	1.0409	1.0603	1.0811	1.1028
True/M1	0.9995	0.9929	0.9851	0.9783	0.9735	0.9717
	1.0006	0.9974	0.9930	0.9887	0.9854	0.9837
	1.0009	0.9999	0.9982	0.9961	0.9943	0.9931
True/M2	1.0284	1.0703	1.1171	1.1664	1.2176	1.2701
	1.0168	1.0428	1.0731	1.1058	1.1401	1.1754
	1.0092	1.0237	1.0412	1.0606	1.0814	1.1031
AR1/M1	0.9916	0.9645	0.9300	0.8902	0.8832	0.8887
	0.9972	0.9858	0.9707	0.9567	0.9481	0.9481
	0.9996	0.9958	0.9903	0.9849	0.9811	0.9802

^[2]