Adjusting for time-varying confounders in survival analysis using structural nested cumulative survival time models.

Accounting for time-varying confounding when assessing the causal effects of time-varying exposures on survival time is challenging. Standard survival methods that incorporate time-varying confounders as covariates generally yield biased effect estimates. Estimators using weighting by inverse probability of exposure can be unstable when confounders are highly predictive of exposure or the exposure is continuous. Structural nested accelerated failure time models (AFTMs) require artificial recensoring, which can cause estimation difficulties. Here, we introduce the structural nested cumulative survival time model (SNCSTM). This model assumes that intervening to set exposure at time t to zero has an additive effect on the subsequent conditional hazard given exposure and confounder histories when all subsequent exposures have already been set to zero. We show how to fit it using standard software for generalized linear models and describe two more efficient, double robust, closed-form estimators. All three estimators avoid the artificial recensoring of AFTMs and the instability of estimators that use weighting by the inverse probability of exposure. We examine the performance of our estimators using a simulation study and illustrate their use on data from the UK Cystic Fibrosis Registry. The SNCSTM is compared with a recently proposed structural nested cumulative failure time model, and several advantages of the former are identified.

INTRODUCTION

Observational studies that attempt to assess the effect of a time‐varying exposure on a survival outcome typically suffer from time‐varying confounding bias. Such bias is the result of time‐varying factors that both influence exposure and are associated with survival, thereby distorting the association between the two. For example, studies of the effect of hospital‐acquired pneumonia on time to death (since hospital admission) in critically ill patients are confounded by disease severity, because disease severity influences susceptibility to pneumonia infection and is strongly associated with mortality (Bekaert et al., 2010). Time‐varying confounders (eg, disease severity) are often affected by earlier exposures (eg, pneumonia infection). This induces feedback relationships between exposures and confounders over time that cannot be untangled via traditional survival analysis regression methods that adjust for time‐varying covariates, such as a history of exposure and confounders, at each timepoint (Robins et al., 2000). The reason for this is twofold. First, such adjustment procedures eliminate the indirect effects of early exposures on survival that are mediated through those confounders. For example, it would be undesirable to eliminate the effects of hospital‐acquired pneumonia on survival that are mediated through disease severity, as scientific interest is primarily in the overall effect of infection. Second, such adjustment procedures are prone to collider‐stratification biases that can render exposure and outcome dependent even in the absence of an exposure effect. See Daniel et al. (2013) for a review of these difficulties.

Time‐varying confounding has received much attention in the causal inference literature. For survival time outcomes, the two main approaches are based on structural nested accelerated failure time models (AFTMs) (Robins and Tsiatis, 1991; Robins and Greenland, 1994) and marginal structural models (MSMs) (Robins et al., 2000). The latter approach is more popular, because of its greater simplicity and flexibility. In particular, accounting for noninformative censoring in MSMs does not, unlike in structural nested AFTMs, require an “artificial recensoring” procedure in which originally uncensored subjects may become censored. Avoiding this recensoring is advantageous, because recensoring causes information loss, which can result in poor estimators and difficulties solving the estimating equations (Joffe et al., 2012). However, fitting MSMs relies on inverse weighting by the probability of exposure, which has it own drawback: estimators prone to large finite‐sample bias and variance when confounders are highly predictive of the exposure, or when the exposure is continuous or discrete with many levels.

More recently, Young et al. (2010) and Picciotto et al. (2012) proposed a new class of discrete‐time structural nested cumulative failure time models, which parameterize the effect of the exposure at each time on the outcome at each later time in terms of the ratio of two (possibly) counterfactual cumulative failure risks at that later time under exposure regimes that differ only at time . Their procedure has the desirable properties of structural nested AFTMs—viz. by avoiding inverse probability weighting, it handles continuous exposures without estimators being subject to large bias and variance, and it allows modeling of effect modification by time‐varying covariates—while avoiding the need for artificial recensoring.

Here, we use developments by Martinussen et al. (2011) and Dukes et al. (2019) (hereafter DMTV). The former showed how to adjust for time‐varying confounding when the effects of exposure and confounders are parameterized on the additive hazard scale. They focused on the simple setting where interest is in estimating the direct effect of a binary baseline exposure on a survival outcome, that is, the effect not mediated by a given intermediate variable, and where there are no baseline confounders. DMTV proposed an additive hazards model for the effect of a baseline exposure on survival time conditional on baseline confounders and derived the efficient score when (as assumed by Martinussen et al. 2011) the confounders act additively on the hazard; this additivity assumption is not needed for consistency of their estimators.

Here, we propose a novel class of semiparametric structural nested cumulative survival time models (SNCSTMs), of which the models of Martinussen et al. (2011) and DMTV are special cases, and propose three estimators of its parameters. Our model allows baseline and time‐varying confounders, binary or continuous exposure, any number of exposure measurement times and the option of constraining exposure effects to be common at different times; it does not parameterize the effects of confounders on the baseline hazard. It also allows investigation of exposure effect modification by time‐varying factors. The SNCSTM is closely related to Picciotto et al.’s model, and our estimators share the forementioned desirable properties of the latter. The SNCSTM generalizes Picciotto et al.’s model to continuous time and parameterizes relative survival risks instead of failure risks. Our approach has several advantages over that of Picciotto et al. One of our estimators (method 1) can be calculated using GLM software. Our other two estimators (methods 2 and 3) are more efficient, double robust and available in closed form. All three estimators automatically handle random censoring. Also, because of being parameterized in continuous time, SNCSTMs can handle irregular measurement times and allow the interpretation of parameters in terms of hazards.

We define notation and state fundamental assumptions in Section 2. A simple version of our SNCSTM is introduced in Section 3. In Section 4, we propose three methods for estimating its parameters. The general SNCSTM is described in Section 5. In Section 6, we discuss random censoring. A simulation study is described in Section 7. Section 8 describes an analysis of data from the UK Cystic Fibrosis (CF) Registry, looking at the effect of treatment with DNase on survival in people with CF. We conclude with a discussion in Section 9.

NOTATION AND ASSUMPTIONS

Consider a study in which, for each of subjects, a time‐varying exposure and vector of possibly time‐varying confounders are measured at time zero and at up to follow‐up visits. Until Section 5 we assume the follow‐up times are regular, that is, the same for all individuals, and (for notational simplicity) are , and that all individuals are administratively censored at time . Until Section 6 we assume there is no censoring apart from this administrative censoring. If visits are regular but not at times , or if administrative censoring occurs at a time different from or not at all, this can easily be accommodated by rescaling the time variable within each interval between consecutive visits.

Let denote individual ’s failure time, and and denote, respectively, his exposure and vector of confounders measured at time (). Let be the at‐risk indicator. If individual fails before his th visit, and are defined as zero. Let , and . The causal ordering of the variables is , where means the minimum of and .

Define as individual ’s (possibly) counterfactual failure time that would have applied if his exposures up to the th visit had been as observed and his exposures from the th visit onward had been set to zero by an intervention. We make the consistency assumption that with probability one for individuals with . Note if and only if for all , that is, intervening on an exposure can only affect survival after the time of that exposure. It follows that events and are equivalent when . We assume () are i.i.d. and henceforth omit the subscript unless needed.

We make the following sequential no unmeasured confounders assumption (NUC): (Robins, 1986). That is, among individuals who are still alive (or event‐free) at time , the assigned exposure at time may depend on and , but given these, has no residual dependence on the remaining lifetime that would apply if all future exposures were set to zero.

STRUCTURAL NESTED CUMULATIVE SURVIVAL TIME MODEL

We first introduce a simple version of the SNCSTM that does not allow for exposure effect modification. The more general SNCSTM is described in Section 5.

For each , let be the model defined by the restriction , where is a vector of unknown parameters, and equals if , equals if , and equals if , and so forth. So, for any equals .

Equation (1) means that among the survivors in the population at the th visit time, in the stratum defined by any given the proportion who survive to a later time when exposures from visit onward (ie, ) have already been set to zero would be multiplied by if exposure were also set to zero. Hence, is the (controlled) direct effect of on the probability of survival to time given survival to visit , that is, the effect of not mediated through the later exposures . For example, if are all positive and , then intervening to set is beneficial, that is, exposure is harmful. Conversely, if are all negative, exposure is beneficial. This SNCSTM assumes the direct effect is the same for any history . In Section 5 we extend the SNCSTM to allow the effect to depend on the history.

By taking logs of each side of Equation (1) and differentiating with respect to , it can be shown that Model can also be written as for (with ), where is the counting process for . Equation (2) can be interpreted as follows. In a stratum defined by and , the hazard of failure at any time between visits and () when have already been set equal to zero would be reduced by if were also set to zero.

Note that Model treats —which, by NUC, equals —as a totally unspecified “baseline” hazard, rather than parameterizing its dependence on and . One advantage of this is that the danger of incompatibility between Models is avoided. To illustrate this danger, suppose it were assumed that for all . This, together with NUC, implies holds. However, it also implies a restriction on the association between and , a restriction which might conflict with that of . Such conflict would be the result of there being no coherent overall model.

ESTIMATION METHODS

In order to estimate , we introduce nuisance Models (). Model is a generalized linear model (GLM) for given and with , where is an unknown finite‐dimensional parameter and is a known vector function of whose first element equals 1, For example, . The dispersion parameter is assumed not to depend on or , and is the canonical link function. The methods described in Sections 4.1 to 4.3 consistently estimate when Models and () are correctly specified. Method 1 can be applied using standard GLM software. Methods 2 and 3 improve on method 1 by using more efficient estimators that are closely related to that described by DMTV in the setting of a single baseline exposure. Method 3 gives consistent estimation under slightly weaker conditions than method 2, but is more computationally intensive.

Method 1: Fitting the GLM implied by Models and

Model states that given and obeys a GLM. Bayes’ rule shows (see Web Appendix A) that Models and NUC together imply that for any given and obeys the same GLM but with the intercept shifted by a function of . Specifically, for , where and (). Our first estimation method for involves fitting this GLM to estimate and calculating . We now explain in more detail.

First we estimate () as follows. For , events and are equivalent, and so Equation (3) implies for any . Hence, a consistent estimate of can be obtained as follows. For each of a number (we used 10) of equally spaced values of between and (including and ), identify the set of individuals with and, for each of these individuals, create a copy (a “pseudo‐individual”) with the same value of and with new random variable equal to . Fit the GLM to the resulting set of (up to ) pseudo‐individuals. A consistent estimate of is then . When is unknown, it can be estimated by fitting Model to those of the original individuals with . In the simulation study of Section 7, we also tried using 50 values of to construct the set of pseudo‐individuals instead of 10, but found this made very little difference to the estimates.

Next we estimate (). When , Equation (3) is . If were known for all could be estimated just as was, but it is not. However, as shown in Web Appendix B, imply that for , where . That is, within the population stratum defined by any given value of and by (or equivalently ), the proportion of individuals with is equal to the proportion of individuals with after weighting each individual by . Remembering that the first element of equals one for all individuals, it follows that a consistent estimate of can be obtained by fitting the GLM to a set of pseudo‐individuals constructed as described above but using 10 values of between and (rather than and ) and using weights . The weights depend on , which is unknown, and so we replace it by its previously obtained estimate . A consistent estimate of is then .

In general, () is estimated by , where is the estimate of obtained by fitting the GLM to a set of pseudo‐individuals constructed using 10 equally spaced values of between and and using weights , with replaced by . For later reference, we denote the fitted value of thus obtained as . This is an estimate of . Note that () must be calculated before . If is unknown, it is estimated by fitting Model to the original individuals with .

Although this estimation procedure involves weights , these are different from the inverse probability of exposure weights used to fit MSMs and do not suffer the same instability that can plague the latter weights. In particular, if , that is, has no direct effect on survival, then . The variance of the weights can be reduced by using modified (or “stabilized”) weights in place of , where and (). This may improve efficiency, especially when is precisely predicted by . The ratio depends only on and , and as model (5) is conditional on these, remains consistent. Since () is unknown, a working model is specified for it and its parameters estimated from the set of individuals still at risk at time . Note that does not need to be correctly specified for to be consistent; indeed need not be compatible with .

Method 2: G‐estimation

The principle underlying the following estimator of is that after removing the effects of and later exposures from the increment in the counting process , NUC implies that the resulting “blipped down” increment at any time is independent of conditional on and and being still at risk.

First estimate () by solving unbiased estimating equation where . The expectation is unknown, so we replace it by , obtained exactly as in method 1. The next paragraph provides a rationale for Equation (6).

NUC implies that the counting process for is conditionally independent of given and . We do not observe but Equation (2) relates to , the counting process for , and we do observe when , because then it equals , the counting process for the observed failure time . In particular, Equation (2) implies that, for any and conditional on , the expected increment in during short time interval given can be unbiasedly estimated by the corresponding mean of the observed increments in minus among the survivors at time . Hence, the adjusted observed increment should be uncorrelated with given and .

DMTV derived the semiparametric efficient estimating equation for under Model assuming known distribution of given and . This equation involves inverse weighting by the hazard function; such weighting also features in efficient estimating equations of other additive hazards models. In practice, accurate estimation of the hazard function is difficult and increases the computational complexity of the procedure, and so this weighting is commonly omitted by standard fitting procedures for additive hazards models. Results of DMTV imply (see Web Appendix C) that if this is done with the semiparametric efficient equation for under Model and if for all , the result is Equation (6).

To make Equation (6) invariant to additive transformations of , we replace by . Since is a constant given , this does not affect the unbiasedness of the estimating equations. Let denote the resulting estimator of .

Next estimate using estimating equation . The unknown and are replaced by and . The next paragraph provides a rationale for this estimating equation.

Again we exploit the conditional independence of and (NUC) and the relation between and but now over time interval . An added complication is that is not observed when . However, we know from Equation (2) that when the intensities of and differ by and (as noted in Section 4.1) there are times as many individuals with in the population as there are with . So, we can unbiasedly estimate the expected increment in over small interval as the weighted mean of the increments in minus with weights . This justifies the above estimating equation but with in place of . We use instead for the same reason that we replaced by in Equation (6).

In general, the consistent estimator of () is obtained by solving with replaced by and () replaced by ; this requires that () be estimated before . The estimator is available in closed form (see Web Appendix E for formulae when is the identity or logistic link function).

In Web Appendix F we prove is double robust in the following sense. Let denote the limit as of , and let Model () be defined by the restriction where are unknown parameters. is consistent if (a) , (b) either or , and (c) for each , either or all of are correctly specified. The term in Model arises because of the use of , rather than , in Equation (7) (see proof). Note that if or is a linear regression, so that is a linear function of , it can be omitted. As in method 1, efficiency may be gained by using stabilized weights in place of in Equation (7). Also, to make invariant to additive transformations of , the term can be replaced by .

Method 3: Improved G‐estimation

If we use a different estimator of for the and terms in Equation (7), then the estimator solving (7) remains consistent under a more general version of Model . In methods 1 and 2, is calculated by fitting a single GLM to a set of pseudo‐individuals, with time since th visit, , as a covariate. In method 3, we instead fit a separate GLM at each time since the th visit. That is, for any , we calculate by fitting the GLM to the set of individuals with , using weights . This set changes only at times at which an individual leaves the risk set, and so the GLM needs to be fitted only at these times. This is the approach taken by DMTV, who denoted the resulting estimator of as “” and, on the basis of results from a simulation study, recommended it over three alternatives. As in method 2, we can use stabilized weights and replace by . As shown in Web Appendix F, method 3 has the same double robustness property as method 2 but with the parameters in Model now allowed to be a function of .

Constraining exposure effects

In some applications, it may be desirable to impose the constraint that for all , that is, the effect of exposure measured at one visit on the hazard visits later is the same for all . This reduces the number of parameters and, as we see in Section 7, increases the precision of their estimates. In Web Appendix G we explain how estimation may be performed under this constraint. See Vansteelandt and Sjolander (2016) for how to impose other constraints.

THE GENERAL SNCSTM

In this section, we extend the SNCSTM to allow visit times to be irregular, that is, to vary from one individual to another, and effect modification, that is, the effect of exposure on survival to depend on the exposure and confounder histories.

Let denote the time of individual ’s th follow‐up visit (), and let () and . Until now, we have assumed . We assume visit times are planned or randomly chosen at baseline using only baseline confounder information , and we modify NUC to be . Also, let denote an administrative censoring time common to all individuals (until now, we assumed ). If there is no such time, let . To allow effect modification, we define , where is a known (possibly vector) function of (“int” stands for “interactions”), and let .

For each , let be the model defined by the restriction where equals if , equals if , and equals if , and so forth. If and , Equation (8) reduces to Equation (1). Model can also be written as for .

The modifications to methods 1 and 2 needed to fit the general SNCSTM are simple (see Web Appendix D). Modifying method 3 is simple when visit times are regular; it is possible for irregular visit times but is fiddly. In the simulation study reported in Section 7 we found little benefit from method 3 relative to method 2 when visit times were regular, and so did not implement it for irregular times.

CENSORING

We now allow for censoring before the administrative censoring time. Let and denote individual ’s censoring and failure times, respectively. Redefine and as and ; is unchanged except that has this new meaning. With these changes, methods 1 to 3 remain valid, provided two further conditions hold (Vansteelandt and Sjolander, 2016). First, the censoring hazard does not depend on the exact failure time or future exposures or confounders. That is, the counting process, , for the censoring time satisfies , where and are the exposure and confounder histories up to time and is some function only of . The second condition, which can be weakened by using censoring weights (see Web Appendix H), is that , so censoring depends only on baseline confounders.

SIMULATION STUDY

We used a simulation study to investigate the bias and efficiency of the methods. There were visits and two time‐dependent confounders (ie, dim). These and the exposure were generated as: and , where . The hazard of failure during the interval between the th and th visits was . For this data‐generating mechanism, () is correctly specified with no effect modification (ie, ) and the true exposure effects are and .

We considered three scenarios: two with regular and one with irregular visit times. For regular visits, . For irregular visits, inter‐visit times were independently uniformly distributed on . There was administrative censoring at time 4. In one of the regular visit scenarios, there was no random censoring. In the other, and in the irregular visit scenario, there was an exponentially distributed random censoring time with mean 5. For the regular visit scenario without random censoring, the expected percentages of individuals observed to fail between visits 0 and 1, 1 and 2, 2 and 3, and between visit 3 and time 4 were 20%, 14%, 11%, and 9%, respectively. For the regular and irregular visit scenarios with random censoring, these percentages were 18%, 10%, 6%, and 4%, and the corresponding expected percentages of individuals censored were 16%, 11%, 8%, and 5%. For each scenario, we generated 1000 data sets, each with individuals. Estimation was done with and without the constraint, which is true here, that .

Tables 1 and 2 show for the regular visit scenario without and with random censoring, respectively, the mean estimates and standard errors (SEs) for methods 1 to 3. Results for the irregular visit scenario are in Web Appendix L. We see that all the estimators are approximately unbiased, though there is some bias for , and , parameters for which there is relatively little information in the data. Comparing SEs, we see that methods 2 and 3 give very similar results and that these methods are more efficient than method 1. This difference in efficiency is much greater when there is random censoring (it is even greater when visit times are irregular—see Web Appendix L). This may be because method 1, unlike 2 and 3, does not distinguish between failure and censoring (or occurrence of next visit). Although methods 2 and 3 use fitted values from the same GLM that is used in method 1, the estimating equations for methods 2 and 3 involve increments , which equal one only when a failure occurs. For methods 1 and 2, coverage of 95% bootstrap confidence intervals (using 1000 bootstraps) was close to 95% (see Table 3). Coverage was not evaluated for method 3, as it is computationally intensive to bootstrap this method for 1000 simulated data sets. Imposing the constraint that reduced SEs, as expected.

Table 1

Means () and SEs () of parameter estimates when , visits are regular and the only censoring is administrative

Mtd	Con	ψ_0(0)	ψ_0(1)	ψ_0(2)	ψ_0(3)	ψ_1(1)	ψ_1(2)	ψ_1(3)	ψ_2(2)	ψ_2(3)	ψ_3(3)
True		0.400	0.100	0.040	0.020	0.400	0.100	0.040	0.400	0.100	0.400
Means
1	No	0.393	0.098	0.031	0.025	0.391	0.096	0.034	0.403	0.098	0.383
2	No	0.396	0.100	0.032	0.024	0.394	0.097	0.033	0.408	0.100	0.392
3	No	0.395	0.100	0.031	0.023	0.392	0.096	0.033	0.406	0.099	0.388
P	No	0.394	0.107	0.030	0.021	0.394	0.094	0.049	0.408	0.102	0.387
1	Yes	0.386	0.096	0.032	0.024	0.386	0.096	0.032	0.386	0.096	0.386
2	Yes	0.397	0.099	0.032	0.023	0.397	0.099	0.032	0.397	0.099	0.397
3	Yes	0.395	0.098	0.032	0.023	0.395	0.098	0.032	0.395	0.098	0.395
P	Yes	0.394	0.104	0.030	0.029	0.394	0.104	0.030	0.394	0.104	0.394
SEs
1	No	0.177	0.187	0.199	0.218	0.243	0.254	0.260	0.251	0.273	0.272
2	No	0.169	0.180	0.191	0.204	0.237	0.246	0.253	0.240	0.262	0.267
3	No	0.169	0.179	0.190	0.204	0.236	0.245	0.252	0.239	0.260	0.265
P	No	0.196	0.290	0.349	0.397	0.265	0.376	0.452	0.270	0.384	0.300
1	Yes	0.113	0.131	0.158	0.217	0.113	0.131	0.158	0.113	0.131	0.113
2	Yes	0.109	0.129	0.151	0.203	0.109	0.129	0.151	0.109	0.129	0.109
3	Yes	0.109	0.128	0.150	0.203	0.109	0.128	0.150	0.109	0.128	0.109
P	Yes	0.126	0.206	0.306	0.494	0.126	0.206	0.306	0.126	0.206	0.126

Note: “Mtd” is method (“P” is Picciotto et al.’s method—see Section 9) and “Con” is whether constraint is imposed.

Table 2

Means () and SEs () of parameter estimates when , visits are regular and censoring is random

Mtd	Con	ψ_0(0)	ψ_0(1)	ψ_0(2)	ψ_0(3)	ψ_1(1)	ψ_1(2)	ψ_1(3)	ψ_2(2)	ψ_2(3)	ψ_3(3)
True		0.400	0.100	0.040	0.020	0.400	0.100	0.040	0.400	0.100	0.400
Means
1	No	0.394	0.108	0.021	0.054	0.396	0.105	0.055	0.403	0.111	0.383
1cw	No	0.396	0.102	0.020	0.054	0.393	0.096	0.054	0.408	0.097	0.383
2	No	0.396	0.104	0.036	0.033	0.399	0.096	0.038	0.411	0.098	0.393
3	No	0.396	0.103	0.036	0.033	0.396	0.095	0.038	0.407	0.096	0.385
P	No	0.397	0.117	0.024	0.050	0.399	0.095	0.078	0.405	0.117	0.390
1	Yes	0.391	0.106	0.031	0.053	0.391	0.106	0.031	0.391	0.106	0.391
1cw	Yes	0.392	0.099	0.031	0.054	0.392	0.099	0.031	0.392	0.099	0.392
2	Yes	0.398	0.099	0.037	0.032	0.398	0.099	0.037	0.398	0.099	0.398
3	Yes	0.396	0.099	0.037	0.032	0.396	0.099	0.037	0.396	0.099	0.396
P	Yes	0.395	0.108	0.035	0.051	0.395	0.108	0.035	0.395	0.108	0.395
SEs
1	No	0.265	0.313	0.372	0.467	0.400	0.483	0.569	0.462	0.563	0.577
1cw	No	0.201	0.234	0.373	0.469	0.298	0.346	0.572	0.348	0.424	0.406
2	No	0.180	0.211	0.252	0.304	0.276	0.313	0.380	0.317	0.385	0.373
3	No	0.180	0.211	0.251	0.303	0.275	0.310	0.375	0.314	0.380	0.367
P	No	0.219	0.389	0.571	0.728	0.334	0.557	0.855	0.385	0.652	0.457
1	Yes	0.186	0.241	0.311	0.463	0.186	0.241	0.311	0.186	0.241	0.186
1cw	Yes	0.140	0.179	0.313	0.465	0.140	0.179	0.313	0.140	0.179	0.140
2	Yes	0.130	0.162	0.211	0.303	0.130	0.162	0.211	0.130	0.162	0.130
3	Yes	0.130	0.161	0.210	0.301	0.130	0.161	0.210	0.130	0.161	0.130
P	Yes	0.157	0.282	0.475	0.802	0.157	0.282	0.475	0.157	0.282	0.157

Note: “Mtd” is method (“1cw” is method 1 with censoring weights; ‘P’ is Picciotto et al.’s method—see Section 9) and ‘Con’ is whether constraint is imposed.

Table 3

Coverage (%) of 95% bootstrap confidence intervals for methods 1, 2, and 1cw (ie, method 1 with censoring weights) when n = 1000, visits are regular, either there is only administrative censoring or there is random censoring, and the constraint is not imposed

Mtd	ψ_0(0)	ψ_0(1)	ψ_0(2)	ψ_0(3)	ψ_1(1)	ψ_1(2)	ψ_1(3)	ψ_2(2)	ψ_2(3)	ψ_3(3)
No censoring
1	96.0	96.0	95.5	94.7	94.4	95.5	96.6	95.4	95.7	94.5
2	96.5	96.4	95.4	95.7	94.9	95.6	96.5	96.0	95.8	94.7
Random censoring
1	95.0	95.6	96.4	94.8	95.3	95.5	95.9	95.6	96.0	95.4
1cw	96.5	96.8	96.6	95.2	95.9	97.9	95.9	97.1	97.8	97.7
2	95.7	95.7	95.9	96.1	94.9	95.9	96.7	95.9	96.6	96.1

In this simulation study, censoring times are independent of exposures and confounders, and so censoring weights (Section 6) are not required for consistent estimation of the ’s. However, applying method 1 with censoring weights improved its efficiency (see method 1cw in Tables 1 and 2), probably because chance associations between exposures and censoring events are reduced in the weighted sample. Coverage of bootstrap confidence intervals (Table 3) was close to 95% for most parameters, but there was overcoverage for some parameters. Using censoring weights had no effect on the efficiency of method 2.

Web Appendix L shows results for or for a shorter follow‐up time with times between visits divided by four and administrative censoring at time 1 (and so fewer failures). These are qualitatively similar to the results in Tables 1 and 2, but with the relative inefficiency of method 1 being even more marked in the scenarios with shorter follow‐up time. Web Appendix L also describes a simulation study that demonstrates the double robustness of methods 2 and 3.

ANALYSIS OF CYSTIC FIBROSIS REGISTRY DATA

The UK CF Registry records health data on nearly all people with CF in the UK at designated approximately annual visits (Taylor‐Robinson et al., 2018). To illustrate the use of the SNCSTM, we used data on 2386 individuals observed during 2008 to 2016 to investigate the causal effect of the drug DNase on survival. DNase has been found to have a beneficial effect on lung function, including using Registry data (Newsome et al., 2019), but its effect on survival has not been studied. Baseline visit was defined as an individual's first visit during 2008 to 2015, and there were up to follow‐up visits. The (irregular) visit times were defined as years after baseline visit; the median time between visits was 1.00 years (interquartile range 0.93–1.07). Individuals were defined as “treated” if they had used DNase since the previous visit and “untreated” otherwise. Individuals treated at a visit prior to their baseline visit were excluded, as were visits prior to age 18. Administrative censoring was applied at the end of 2016 and nonadministrative censoring when an individual had a transplant or had not been seen for 18 months. The percentage of treated patients increased from 14% at the baseline visit to 52% at visit 8, and most patients who began using DNase continued to use it. There were 137 deaths during follow‐up and 653 nonadministrative censorings (including 36 transplants). Of those who died, 74 (63) were treated (untreated) at the time of death. Total follow‐up was 12 380 person‐years (py), and death rates while treated and untreated were, respectively, 0.019 (74/3930) and 0.0075 (63/8450) . The ratio of the probabilities of surviving for one year while treated and untreated is thus . However, this may be due to confounding: sicker patients being more likely to receive treatment.

We estimated the effect on survival of delaying initiation of treatment by one year. To do this, we (re)defined as for those treated at visit , and for those untreated. Now represents the multiplicative causal effect of intervening to start treatment at visit rather than at visit on the probability of surviving for at least one year after visit , among patients who survive to, and are untreated at, visit and conditional on confounder history . More generally, is the effect on the probability of surviving at least years after visit if visits are exactly annual. We imposed the constraint . (Potential) confounders at visit were baseline variables sex, age, and genotype class (low, high, and not assigned), and time‐varying variables %, body mass index, days of IV antibiotic use, and binary indicators for four infections (Pseudomonas aeruginosa, Staphylococcus aureus, Burkholderia cepacia complex, and Aspergillus), CF‐related diabetes, smoking, and use of other mucoactive treatments and oxygen therapy. The same variables (and treatment) were included in models for the inverse probability of censoring weights.

Figure 1A shows estimates of from method 2. These suggest that starting treatment now rather than waiting may cause a small decrease in the probability of survival, at least for the first 5 years: , 0.996, 0.997, 0.994, and 0.988 for , respectively. However, the confidence intervals (obtained by bootstrapping) include 1, that is, no treatment effect. This lack of a significant treatment effect may be because we have focused on a subset of the population (adults not previously treated with DNase) and/or because there are unmeasured confounders. As expected, method 1 was very inefficient in this situation of irregular visits and substantial censoring. The confidence intervals were between 4 and 9 times wider than those from method 2.

Figure 1

Estimates of the ratio of the survival probabilities when treatment is initiated immediately compared to initiation being delayed by one year. A: from the model with no interaction. B, C and D: from the model with interaction between treatment and

For illustration, we also fitted an SNCSTM with an interaction between treatment and %. Figure 1B to 1D shows the estimated ratios of survival probabilities for three values of %: 40, 75, and 100 (the 10th, 50th, and 90th centiles of the distribution at baseline). Figure 1D suggests the ratio may actually be greater than 1 for % , that is, starting treatment now may be better than waiting for patients with high %. However, the interaction terms are not significant.

DISCUSSION

One advantage of SNCSTMs is that, in contrast to MSMs, they can cope well with situations where the inverse probabilities of exposure are highly variable. Indeed, they can even be used when the so‐called experimental treatment assignment assumption is violated, that is, when some individuals are, on the basis of their time‐varying covariate information, excluded from receiving particular exposure levels. For these individuals, , meaning they do not contribute to the estimating functions of methods 1 to 3.

Another advantage of SNCSTMs is that they can be used to investigate time‐varying modification of exposure effects on survival time. Although it is, in principle, possible to do this using structural nested AFTMs, estimation difficulties caused by artificial recensoring mean that such models are usually kept simple and interactions are not explored.

The SNCSTM can also be used to estimate the counterfactual exposure‐free survivor function, that is, , as . This is because Equations (4) and (8) imply . If there is censoring before time should be inversely weighted by an estimate of .

A limitation is that, like other additive hazards models, the SNCSTM does not constrain hazards to be nonnegative, and so does not exclude survival probabilities greater than one. Similarly, Picciotto et al.'s (2012) structural nested cumulative failure time model does not exclude failure probabilities greater than one.

Method 1 appears to be less efficient than methods 2 and 3 but has the appeal that it can be applied using standard GLM software. In our simulation study, the efficiency loss was fairly small when the only censoring was administrative and visit times were regular. This method became much less competitive, however, when there was random censoring, and even more so when visit times were irregular. By not distinguishing between failure and censoring, method 1 may also be more sensitive than methods 2 and 3 to violation of the assumption that . Of the three, method 3 gives consistent estimation under the weakest assumptions. However, it needs more computation than methods 1 and 2, especially when visit times are irregular and the exposure is binary. In our simulation study, methods 2 and 3 performed similarly, and so the theoretical advantage of method 3 may not be worth the extra computation. An R function for implementing our methods, with examples, is described in Web Appendix I.

DMTV discuss the close connection between their model for a point exposure (which is equivalent to the SNCSTM with ) and Picciotto et al.'s (2012) cumulative failure time model. Although the latter is a discrete‐time model for the probability of failure, it is easy to finely discretize time so as to approximate continuous time and (as Picciotto et al. note) to reformulate it as a model for the probability of survival. As DMTV explain, a drawback of Picciotto et al.'s method is the difficulty of deriving the efficient estimating equation. This difficulty arises because their class of estimating functions uses correlated survival indicators. By instead using independent increments of a counting process, DMTV were able to derive the efficient estimating function. Methods 2 and 3 are extensions to time‐varying exposures of DMTV's recommended method, and are, therefore, expected also to be more efficient than Picciotto et al.'s method. In Web Appendix J we elaborate on DMTV's discussion of Picciotto et al.'s model and reformulate it as a model for the probability of survival. Tables 1 and 2 show mean estimates and SEs for the resulting Picciotto et al. estimator (described in Web Appendix J and denoted “Method P” in tables). The SEs are larger than those of methods 2 and 3, suggesting methods 2 and 3 are indeed more efficient. Methods 2 and 3 also have the advantages of using closed‐form estimators, handling random censoring automatically (because estimating functions are framed in terms of increments, which are observable up to the time of censoring), and being double robust. Picciotto et al. use an iterative Nelder‐Mead algorithm, employ inverse probability of censoring weighting to handle random censoring, even when this censoring is completely at random, and their estimator is not double robust.

In Web Appendix K we outline how the SNCSTM can handle competing risks.

Web Appendices referenced in Sections 4–7 and 9 are available with this paper at the Biometrics website on Wiley Online Library.