We evaluated the finite sample behavior of our estimation procedure when a key confounder was missing among those who die. We considered a case with two confounding covariates, X1 (continuous; fully observed) and X2 (binary; missing on those who die). For our simulation, we generated data (under Assumptions 1-6) for an individual according to the following scheme: (1) generate X1 as Normal with mean 0 and variance 1; (2) generate X2 as Bernoulli with P[X2 = 1|X1] = expit(0.5 + 0.25X1); (3) generate Z as Bernoulli with P[Z = 1|X1,X2] = expit(-2 + 2.5X1 + 3.5X2); (4) generate D(0) as Bernoulli with P[D(0) = 1|X1,X2] = expit(-0.5 + 2.5X1 + 3X2); (5) if D(0) = 0 set D(1) = 0; otherwise, generate D(1) as Bernoulli with P[D(1) = 1|D(0) = 1, X1, X2] = expit(-2.25 + X1 + 1.75X2); (6) if D(0) = 0, generate Y (0) as Normal with E[Y (0)|D(0) = 0, X1,X2] = 40 + 5X1 + 40X2 and variance 1; (7) if D(1) = 0 generate Y (1) as Normal with E[Y(1)|D(1) = 0,X1,X2] = 33 + X1 + 20X2, and variance 1; (8) set D = D(Z) and Y = Y (Z) if D = 0.
Using the analysis model logitP[D = 1|Z,X1,X2] = β0 + β1Z + β2X1 + β3X1Z + αX2 + τX2Z, the true values of α and τ are approximately 3.0 and -1.0. Since these values are unknown when X2 is missing on those who die, we conducted a sensitivity analysis, where we considered combinations of these parameters in a 4 × 4 square, centered at the truth. We used a linear (analysis) model for the conditional (on Z, X1, and X2) mean of Y among observed survivors, where the functional form of the linear predictor was the same as that used in death model above. For fixed α and τ, we estimated SACE and derived confidence intervals using the large sample theory approximation for the variance as described in Section 6. We evaluated the proportion of confidence intervals that contained the true value of SACE = -10.4.
In our data generation scheme, higher values of Y (0) and Y (1) were considered worse outcomes and higher values of X1 and X2 were indicative of poorer health and worse outcomes. Those with poorer health were more likely to be assigned to the treatment group. The probability of being assigned to treatment was approximately 53%. The parameters resulted in approximately 53% of the sample being assigned to Z = 1, and 62% having X2 = 1. In , we present treatment-specific distributional summaries of key variables.
Distributional summaries of simulated variables, stratified by treatment assignment
The unadjusted naive effect of the intervention indicates that the intervention is harmful (E[Y|D = 0, Z = 1] - E[Y|D = 0, Z = 0] = 2.8). Adjusting for X1 and X2 using an analysis model for E[Y|D = 0,Z,X1,X2] which is linear in Z, X1, and X2 and includes no interactions between Z and (X1, X2), yields an overly optimistic treatment effect of -16.7. The effect among the survivors under each treatment (E[Y (1)|D(1) = 0] - E[Y (0)|D(0) = 0]) is -5.99, a more modest effect. This example demonstrates the importance of estimating SACE even when observed death rates between the treatment arms are relatively similar (38% vs. 43%).
For each of three sample sizes, 500, 1000, and 2000, we generated 1000 simulated datasets. displays the results. Coverage of the 95% confidence intervals was good when α and τ are fixed at their true values: 94.9%, 93.4%, and 93.7% for the three respective sample sizes. The coverage worsened when α and τ were fixed at values different than the truth. Coverage was particularly poor when α was assumed to be smaller than the truth. In terms of bias, the trends were very similar: there was no evidence of bias when the sensitivity analysis parameters were correctly specified and absolute bias increases substantially when α was assumed to be smaller than the truth.
Figure 2 Results of simulations with the 95% coverage indicated using contour lines and the absolute bias of the point estimates represented by shading. The figures are centered at the true values of the sensitivity parameters, denoted α* and τ (more ...)
The top left corner of the figures depicts the region in which α and τ are incorrectly assumed to equal zero. It is in this region where the bias and empirical coverage are worst. Specifically, the results are biased towards a more beneficial treatment benefit. The reason that the treatment estimates are more beneficial as α decreases is that such an assumption reduces the probability that those with high values of X2 will die. Since higher values of X2 are associated with a greater treatment effect, the point estimates of the intervention become stronger as more individuals with higher values of X2 are assumed to survive.
Our analysis model for the conditional probability of death, logitP[D = 1|Z,X1,X2], might be incorrectly specified in this simulation. The true values of α and τ are those that minimize the Kullback-Leibler distance between the true and (possibly) misspecified model. They are found by simulating a massive dataset and estimating the value of these parameters via maximum likelihood. The coverage probabilities and bias at these values of α and τ suggest that our analysis model works well. Further, when we performed simulations using a more flexible analysis model that entered X1 into the model using restricted cubic splines, our results did not substantively change. This suggests that misspecification of the analysis model did not affect our simulation inferences.
Overall, the simulations suggest the utility of the sensitivity analysis approach. There is little bias and good coverage when the sensitivity parameters are correctly specified. Since the true value of the sensitivity parameters are not identifiable without additional assumptions and/or data, it is important to present inferences over a range of these parameters.