The Effect of Error-in-Confounders on the Estimation of the Causal Parameter When Using Marginal Structural Models and Inverse Probability-of-Treatment Weights: A Simulation Study

2.1 Counterfactuals and the marginal structural model

Counterfactuals, or potential outcomes, are a useful tool for understanding MSMs and often, they are presented in the context of time-varying covariates and treatments [10–12]. With our focus on MSMs for point-treatment studies, we have a single interval setting where data are collected at the initiation of the study , and at a successive point in time, . We let L denote the vector of baseline covariates that were measured at before follow-up or final assessment at . We define the treatment decisions A to be the treatment assignment occurring between times and ; throughout, we shall take . Finally, we define Y as the continuous outcome measured at .

We want a model describing the mean of Y for each treatment history, a, so we let denote the potential response that would be observed if the subject followed treatment a. We begin with the model

where parameterizes the model and the function g is a link function. In general, the past treatments can be incorporated into the systematic component using some function of the treatment history, , thus . The observed data for the ith subject are . We note that the point-treatment study is generalizable to k intervals [13].

Under the assumptions of consistency, exchangeability [14], positivity [15], and time ordering where treatment precedes the outcome [16], we can obtain an unbiased estimate of . We accomplish this by fitting a model to the observed data where, under the conditions of consistency and exchangeability, the marginal expectation of the counterfactual is equal to the conditional expectation of the response given the treatment, . Additionally, we assume that the n subjects are a random sample from some superpopulation of subjects, which is nearly infinite in size, about whom we wish to make inference [14].

2.2 The inverse probability-of-treatment weight

The assumption of no unmeasured confounders is an essential condition that permits causal inference from observational data [12]. Under this assumption, and in the absence of model misspecification, statistical exogeneity implies causal exogeneity [14]. A measure of statistical exogeneity for each subject can be used to weight each observation so that each subject represents that which would have been observed in the superpopulation with an ancillary treatment process [12].

The weight, called the stabilized weight, can be interpreted as the inverse of the incremental effect of the time-varying confounders on the current treatment beyond other treatment determinants including previous history and baseline covariates [1, 2]. When subject information is not observable, such as loss to follow-up, a censoring process runs concurrent with both the treatment and the covariate processes. Following Hernán et al. [17], we define the censoring process at time , as C. It occurs between times and and precedes treatment allocation; with defining data observation. The weight for point-treatment studies is specified as , where

(1)

Robins [12] shows that the weight corrects the regression estimator for the effect of confounding and censoring due to L. It is naturally expressed by considering the true joint density of the factual and counterfactual data, and it encodes information ensuring that the treatment history is causally exogenous, thus permitting the equivalence .

2.3 Classical unbiased nondifferential measurement error

We apply the classical, unbiased, additive, nondifferential measurement error model to confounders. For the ith subject and kth confounder, we have where is the mismeasured version of , , and the errors are independent, for , and ; . Under the assumption that or equivalently and Y are independent given L, the model is nondifferential [18, 19]. Furthermore, we assume that if and are mismeasured, then the associated errors are independent, .

For this situation, the relevant point-treatment directed acyclic graph (DAG) for classical measurement error is shown in Figure 1, where denotes exogenous measurement error noise. Here, L d-separates from , thus by the probabilistic implications of d-separation [20, Theorem 1.2.4], is independent of conditional on L in every distribution compatible with the specified DAG (Figure 1). Furthermore, for all compatible distributions, this conditional independence implies surrogacy and nondifferential measurement error [18, 21].

When , rather than L, is observed, then the inverse probability-of-treatment weight with a censoring process is

(2)

Using rather than L, the causal estimator may be biased. We are solving and not , where is a function of the treatment history and response parameterized by .

Figure 1

Nondifferential unbiased confounder measurement error for point-treatment directed acyclic graph, where denotes the measurement error noise such that

3.1 Primary study

For our investigation, we focus on the point-treatment case. Figure 2(a) and 2(b) are the DAGs corresponding to the data generating mechanism with and without censoring, respectively. The dashed line in Figure 1 indicates that the censored data indicator was not included as a covariate in the data generating model, but affects what is observed.

Table 1 gives the six scenarios we considered. Scenarios S1–S3 had no censoring, using the data generating mechanism in Figure 2(a). Scenarios S4–S6 had censored data (Figure 2b). Scenarios S2 and S5 had only one mismeasured confounder, whereas scenarios S3 and S6 had both confounders mismeasured (Table 1).

Figure 2

DAGs for the simulation study data generation mechanisms, where and denotes the measurement error noise: (a) without censoring and (b) with censoring. The dashed line indicates that the censored data indicator was not included as a covariate in the data generating model, but that it affects what is observed

Table 1

Simulation study scenarios

Scenario	Censoring	mismeasured	mismeasured
S1	No	No	No
S2	No	Yes	No
S3	No	Yes	Yes
S4	Yes	No	No
S5	Yes	Yes	No
S6	Yes	Yes	Yes

The vector of confounders, , , is assumed bivariate Normal, , where and

where and , the correlation between and . We chose uncorrelated confounders to eliminate, or at least minimize, any induced augmentation [21, section 2.4]. We let A be binary, with denoting the event of interest (e.g. new treatment) and indicating missing data.

Table 2 identifies the observed data vectors for the six scenarios. Logistic regression was used to estimate the probability models given in Table 2. For example, we used for scenario S1 and for scenario S4, we had the additional model for the censoring mechanism. These two models have the basic structure except for the specification of mismeasured covariates (Table 2). Inverse probability-of-treatment weights were estimated by incorporating the probability models specified in Table 2 into eqs (1) and (2) for a point-treatment process, for example, scenario S4 used the stabilized weight of eq. (1). These weights were used to estimate the causal effect of A on

Table 2

Observed data and associated probability models for the six simulation study scenarios

Scenario	Observed data	Probability model
S1
S2
S3
S4
S5
S6

The response is continuous and specified as , and . Given the distribution of the confounders, the marginal model parameters are

where , thus and . We let L be mismeasured such that , where , and , where , such that and .

We implemented an extensive simulation study, using 500 randomly generated parameter value combinations within a defined space, to investigate the effect of measurement error on the causal parameter. Parameters were selected from the space defined by . We set and randomly selected the remaining parameters from the uniform distribution, , rounded to three decimal places. By design, the product space associated with the seven free parameters permitted independent sampling within each dimension. For each of the 500 locations in the parameter space, a simulation study was implemented using 5,000 simulated data sets with each data set having a sample size of 10,000.

We estimate the standard error of the bias using the empirical sampling distribution of biases, and the bootstrap. The estimated bias is where and . The 95% simulation based confidence intervals for and , across the six scenarios and four cases, were derived using the respective empirical distributions. The simulation study was conducted using R [22] with 4-GB RAM at 1,067 MHz DDR3 and a 2.8-GHz Intel Core 2 Duo processor.

3.2 Secondary study

We conducted additional simulations to investigate further the results of the primary study. We wanted to see if the findings of Gustafson [21] applied to the estimation of a MSM causal parameter. Specifically, we investigated the effects of confounder measurement error variability and confounder correlation. We chose the sign reversal scenario, as discussed in Section 4, where = (0, 1.037, 0.893, 0,1.86, –0.683, 0, 2.924, 2.661, –0.455). This parameterization indicates that and . We varied the parameters , , and , permitting and , the standard deviations associated with the measurement error models, to take values in the set and specified correlations between and , .

4.1 Attenuation and augmentation

Attenuation and augmentation are known effects of measurement error [18, 21]. As expected with attenuation, there is an increase in the magnitude of the downward bias as confounder measurement error is introduced (scenarios S2 and S3 and S5 and S6 of Figure 3). As suggested by Gustafson [21], we observed greater attenuation as the number of mismeasured confounders increased (Table 3). Although bias is noted for scenarios S2 and S5, Table 4 suggests that the induced bias may not be meaningful with respect to the variability of the distribution of . The opposite is true for scenarios S3 and S6 (Table 4). We found that the standard error of the bias of was relatively small, suggesting a meaningful bias as well as the potential for extreme values to unduly influence the construction of simulation based confidence intervals when confounder measurement error occurs (Figure 3).

As with attenuation, the effect of augmentation was anticipated. The augmentation in Figure 4, scenarios S2 and S5, is slightly less pronounced than that of attenuation with much less variability of the empirical distribution of (Table 4). The augmentation effect is similar with and without censoring (Table 4). In scenarios S3 and S6 (Figure 4 and Table 4), the addition of an extra mismeasured confounder had little effect on the estimate of and on the percentile interval. Table 4 shows that in all four scenarios, the intervals fail to capture the true value of . Tables 3 and 4 reveal that attenuation and augmentation of the causal parameter, , do not translate to the same measurement error effect on intercept, .

4.2 Null effect and sign reversal

Two unanticipated effects were observed: a null effect and sign reversal. A null effect occurs when the use of mismeasured confounders has little or no impact on the estimator (Figure 5). There are two distinguishing features: unbiased parameter estimation and small variation in the empirical distribution of . The case shown in Figure 5 has little variation in the estimates across the 5,000 simulations (Table 3), which, from our set of simulations, does not appear to be atypical of the underlying distribution of for this effect. We observe no added attenuation or augmentation with more than one confounder measured with error.

The second is sign reversal (Figure 6). Although sign reversal (trend reversal) is discussed in the context of misclassification [3, 23, 24], it is not well documented for classical measurement error models and perhaps never observed due to the strong analytic and empirical results that support attenuation and augmentation [18, 19, 21, 25]. In scenarios S2 and S5 of Figure 6, we observe that is very close , but the sign is reversed, being positive rather than negative (Table 4). Scenarios S3 and S6 exhibit an increased effect of attenuation with the use of a second mismeasured confounder, suggesting that the anticipated effects of attenuation and augmentation are functioning in an intuitive manner despite the reversed sign. There is a tighter distribution of the empirical distribution of for scenarios S2, S3, S5, and S6, when compared to scenarios S1 and S3 (Table 3).

We used the 500 studies to estimate the occurrence of the four measurement error effects and defined the effects to be mutually exclusive. Null events are locations where the confidence interval covered the truth. The sign reversal events are defined as all those locations where the confidence interval did not capture the truth and there was a sign reversal. The augmentation and attenuation events are mutually exclusive by definition, but here they are defined conditional on being neither a null nor a sign reversal event.

We observe that augmentation events are the most frequent across all the scenarios and are followed closely by attenuation events (Table 5). When going from a model with only one mismeasured confounder to a model with two mismeasured confounders, there is a reduction in the proportion of null events and an increase in the other three event types. Attenuation events have the largest increase.

4.3 Secondary study

Using the sign reversal example, scenario S6, we investigate the observed effects of confounder measurement error on the bias of as a function of , , and , (Figure 7). To orient ourselves within the contour plots, scenario S6 corresponds to Figure 7 where (top left), and ; scenario S6 specified a location between the 0.6 and 0.8 contours. This location corresponds to region where we should expect to see sign reversal combined with attenuation. Furthermore, we make the general observation that as increases, the distance widens between the contour lines.

Figure 7

Contour graph of the bias as a function of and when for the sign reversal example

We are able to demarcate four effect areas with three boundary lines. The first boundary occurs along the contour line where the bias is 0. To the left of this line is the area for which combinations of and , given the specified regression coefficients, will result in augmentation. To the right is the area that yields attenuation. This area, in our example, is bounded on the right by the contour line where the bias is 0.45. To the right of this boundary line is the area where we would see attenuated sign reversal, that is the magnitude of the estimated regression parameter is less than that of the true parameter and the sign is opposite, and . This region is bounded on the right by the contour defined at 0.9. Here, we transition from attenuated sign reversal to augmented sign reversal, and , with equality holding on the boundary.

In our secondary study, changes in appear to be more influential than changes in , across all levels of , for determining the type of effect seen. This may be related to the chosen parameterization of the censoring model. For the sign reversal example, we had , thus changes in the imprecision of measurement for the more important variable have more of an impact on the type of effect observed. This is analogous to the effect of the measurement error for multiple mismeasured predictors using linear regression when the more important variable is measured with more error than others in the model [21, section 2.4].

The secondary study afforded the opportunity to explore the stabilized weight profile of each of the four types of measurement error effects within the sign reversal scenario (Figure 6). With this approach, we are able to fix and investigate how the stabilized weights change across the four types of measurement error. We obtained the 50th, 90th, 95th, and 99th percentile distribution of the weights using the 5,000 simulations; these were summarized using the median. The results for the 90th and 95th percentiles do not deviate much from those of the 50th percentile, so we only focus on the 50th and 99th.

The profile of the 50th percentile of the stabilized weights, , shows an ordering for the weight profile as a function of (Figure 8). All four effects exhibit a positive trend as increases. We observe an ordering of the minimum median stabilized weight: null, augmentation, attenuation, and sign reversal. This ordering is maintained for the maximum median stabilized weight. Despite these trends, the weights are similar in magnitude across all four effects.

Figure 8

Scatterplot of the 50th percentile of the stabilized weights, , for the four measurement error effects as a function of the confounder correlation, , where M denotes the median and denotes the 50th percentile

The differences between the four measurement error effects become more pronounced when we consider the 99th percentile, the tail of the stabilized weight distributions (Figure 9). All four effects show that the stabilized weights trends are decreasing as increases, which is opposite to the trend observed with the 50th percentile (Figure 8). We observe the reverse ordering of sign reversal to null effect with sign reversal having the smallest median at each correlation. Furthermore, we observe a larger difference in the magnitude of the median of the 99th percentile at each level of . The unexpected result was that the sign reversal had the smallest median weights (Figure 9), which suggests that the profile of the tail of the weight distribution may be important in determining the observed measurement error effect. These empirical results suggest that the effect of confounder measurement error is a function of the model parameters, covariance structure of the confounders, and the noise of the measurement error models, which conceptually agrees with Gustafson [21, Result 2.3, eq. (2.8)] and with Buonaccorsi [18, eq. (6.38)].

Figure 9

Scatterplot of the 99th percentile of the stabilized weights, , for the four measurement error effects as a function of the confounder correlation, , where M denotes the median and denotes the 99th percentile