Instrumental Variable Estimation with the R Package ivtools

4.1 Two-stage estimation

As the name suggests, two-stage estimation requires fitting two regression models. In the first stage, a regression model is fitted for the exposure, using L and Z as regressors. This model is used to create a prediction Xˆ=Eˆ(X|Z,L) for each subject. In the second stage, another regression model is fitted for the outcome, using L and mT(L)Xˆ as regressors. The choice of regression model in the second stage mimics the causal target model. That is, for causal models (1), (2) and (4) we use a GLM with link function η, a Cox proportional hazards model, and an additive hazards model, respectively, in the second stage. Two-stage estimation of an unrestricted function B(t) in model (3) is currently not available. The estimated coefficient vector for mT(L)Xˆ in the second stage regression model is the two-stage estimate of ψ in the corresponding causal model. The standard error of the two-stage estimate may be obtained by stacking all estimating equations involved in the estimation process, and applying the ‘sandwich formula’ to the whole equation system (Stefanski and Boos 2002).

The rationale behind two-stage estimation is most easily understood in linear models, where the outcome can be thought of as a sum of two terms; the exposure X and an ‘error term’ ε. By confounding, the error term is associated with the exposure. Intuitively then, the first stage can be thought of as extracting the part of the variation in X that is independent of ε, since, under the IV assumptions, the IV is (conditionally) independent of the confounders (given L). In the second stage, we then use this ‘exogenous’ part of the variation in X to estimate the causal effect ψ.

An advantage of two-stage estimation is that it is computationally simple. It gives unbiased estimates in linear models, that is, in model (1) with identity link, in which case it reduces to standard TSLS, and in model (4) if the first stage model is a linear regression (Tchetgen Tchetgen et al. 2015). However, several authors have noted that two-stage estimation generally gives biased estimates in nonlinear models, such as model (1) with logit link, and in model (2). This bias essentially arises from the non-collapsibility of non-linear effect measures, such as the odds ratio and the hazard ratio (Greenland et al. 1999). Under some circumstances, typically when the outcome is rare, the non-collapsibility may not be severe, in which case the bias may be small. The bias can sometimes be reduced by adding the ‘control function’ R=X−Eˆ(X|Z,L) to the set of regressors in the second stage model (Vansteelandt et al. 2011; Tchetgen Tchetgen et al. 2015).

4.2 G-estimation

G-estimation is a fairly general estimation technique, which includes, for instance, semiparametrically efficient and doubly robust estimators. We here describe a special case of G-estimation, which is implemented in the ivtools package. We refer to Robins (1989, 1994); Vansteelandt and Goetghebeur (2003); Martinussen et al. (2017, 2018) for more general expositions on G-estimation in IV analyses.

As for two-stage estimation, the rationale behind G-estimation is most easily understood in linear models. Under the causal model (1) with identity link, we may define the residual hi(ψ)=Yi−Xiψ. When the causal effect Xiψ has been removed from the factual outcome, this residual can be thought of as a prediction of the potential outcome Y0i. Under the IV assumptions, this potential outcome is (conditionally) independent of the IV (given L). Thus, under the IV assumptions we may estimate ψ as the value which makes this independence happen in the sample.

Formally, for model (1) a G-estimator of ψ can be defined as the solution to the estimating equation

where

and d(L, Z) is an arbitrary function with the same dimension as ψ, such that E{d(L,Z)|L}=0, see Robins (1989, 1994) and Vansteelandt and Goetghebeur (2003).

[Correction added on July 20, 2019 after online publication: In the above equation “logit” was added before E^(Y|Li,Zi,Xi).]

By similar reasoning, a G-estimator of ψ in model (2) can be defined as the solution to the estimating equation

where t is a fixed time-point, and

see Martinussen et al. (2018). The G-estimator in (6) depends on the choice of t. If both the causal model (2) and the model for p(T>t|L,Z,X) are correct, then the choice of t may affect the precision, but not the unbiasedness, of the G-estimator. Martinussen et al. (2018) argued that t may be chosen as the time-point that minimizes the estimated variance of the G-estimate.

The ivtools package allows for two choices of d(L, Z):

and

The function in (7) has the advantage of requiring less modeling, i. e. it only requires a model for E(Z|L) whereas the function in (8) requires models for both E(X|L,Z) and E(X|L). However, the function in (7) requires Z to be a scalar variable, whereas function in (8) allows for Z to a vector of arbitrary length. In the real data analysis (Section 5) we only present G-estimation with (7), and in Appendix C we repeat the analysis with (8).

The estimating functions in (5) and (6) require additional regression modeling. Specifically, they require estimates of E(Z|L) (for d(L, Z) in (7)), E(X|L,Z) and E(X|L) (for d(L, Z) in (8)), E(Y|L,Z,X) (for (5)) with logit link) and p(T>t|L,Z,X) (for (6)), which are obtained by regression models. To avoid bias due to model misspecification, these models may be saturated, if possible.

The standard error of the G-estimator estimate may be obtained by stacking all estimating equations involved in the estimation process, and applying the ‘sandwich formula’ to the whole equation system (Stefanski and Boos 2002).

Martinussen et al. (2017) derived G-estimators of B(t) and ψ in models (3) and (4), respectively. These estimators follow similar intuition as the G-estimators of ψ in models (1) and (2). However, they are more complex, as they are defined as recursive estimators of counting process integrals, rather than solutions to estimating equations. We refer to Martinussen et al. (2017) for details on these G-estimators.

We end this section with a remark on terminology. We have here defined ‘G-estimator’ as the class of estimators which make the potential outcome Y0 (or T0) independent of Z in the sample, conditionally on L. Under this definition, all estimators in this section are indeed ‘G-estimators’. However, there appears to be no uniform agreement on the definition, and some authors seem to reserve the term for estimators that solely rely on a model for the IV. Under this more restrictive definition, the estimators of ψ in model (1) with logit link and model (2) are not G-estimators, since they additionally require models for E(Y|L,Z,X) and p(T>t|L,Z,X), respectively.

5.1 The Vitamin D data

All functions that we use for the real data analysis are included in the ivtools and survival packages; we load these as usual by typing

We use a data set borrowed from Martinussen et al. (2018), on a cohort study of vitamin D and mortality. To estimate the causal effect of vitamin D on mortality, these authors carried out an IV analysis, using mutations in the filaggrin gene as an instrument. These mutations have been shown to be associated with higher vitamin D status (Skaaby et al. 2013). For ethical reasons, we were not able to make the original data publicly available. Thus, we have mutilated the original data slightly, by adding random noise to each variable. These mutilated data are included in the ivtools package, as the VitD data frame. We load these data by typing

The VitD data frame contains 2571 subjects and 5 variables: age (at baseline), filaggrin (binary indicator of whether filaggrin mutations are present) vitd (vitamin D level at baseline, measured as serum 25-OH-D (nmol/L)), time (follow-up time), and death (indicator of whether the subjects died during follow-up). The median follow-up time is 16.1 years, and 23.5 % of all subjects died during follow-up.

In Sections 5.2 and 5.3 we analyze death during follow-up as a binary point outcome. This analysis is mainly for illustration; since the exact time of death is known it is more natural to view the outcome as a time-to-event. We thus proceed in Sections 5.4 and 5.5 with more appropriate time-to-event analyses. We control for age at baseline in all analyses. We follow Martinussen et al. (2018), and define the vitamin D exposure as a continuous variable centered around 20 and normed with 20:

The analyses below require fitting several regression models. We use the naming convention fitA.B, where A is the name of the left-hand side variable in the model, and B are the names of the right-hand side variables in the model. For instance, fitY.LX is a model with Y on the LHS and (L, X) on the RHS, whereas fitY.LZX is a model with Y on the LHS and (L, Z, X) on the RHS.

5.2 Two-stage estimation with a point outcome

We start by using standard logistic regression to estimate the association between vitamin D and mortality, while controlling for baseline age. We fit the model, store the results in an object called fitY.LX, and summarize by typing

We observe that higher vitamin D levels are associated with lower mortality, with a decrease of 0.17 units in the log odds of death during follow-up, for each 1 unit increase in standardized vitamin D. This association is statistically significant, with a p-value equal to 9.39×10−5.

Since we only controlled for age in the logistic regression, we may worry that the observed association is partly or fully due to unmeasured confounding. We thus proceed with an IV analysis, fitting the causal model (1) with η being the logit link and m(L) = 1. For this purpose we use the ivglm function in the ivtools package, which allows for both two-stage estimation and G-estimation.

To carry out two-stage estimation we need one model for vitamin D, regressed on age and filaggrin (the first stage model), and one model for death, regressed on age and vitamin D (the second stage model). For the latter, we use the logistic regression model fitted above. For the former, we use an ordinary linear regression. We fit this model, store the results in an object called fitX.LZ, and summarize by typing

We observe that presence of mutations in the filaggrin gene are associated with higher vitamin D levels, as expected. We now use two-stage estimation to fit the causal model (1), by typing

The estmethod argument specifies the desired estimation method; "ts" or "g" for two-stage estimation and G-estimation, respectively. For two-stage estimation, three additional arguments must be specified; fitX.LZ, fitY.LX and data, which specify the exposure model, the outcome model and the data frame, respectively. The exposure and outcome models may be arbitrary GLMs, as fitted by the glm function. The ivglm function uses the exposure model to create predictions, and re-fits the outcome model with the true exposure levels replaced by these predictions. The link function η and the interaction function m(L) in the target causal model (1) are implicitly defined by the outcome model fitY.LX.

We summarize as usual, by typing

We observe that the estimated effect of vitamin D is quite substantial, with a decrease of 1.39 units in the log odds of death during follow-up, for each 1 unit increase in standardized vitamin D. However, this effect is not statistically significant.

The summary method for "ivglm" objects produces the same statistics as the summary method for "glm" objects. It would often be desirable to compute confidence intervals as well. This can be done with the confint method:

We remind the reader that the two-stage estimate of ψ in model (1) is generally biased, unless η is the identity link. The bias can often be reduced by adding the ‘control function’ R=X−Eˆ(X|Z,L) to the set of regressors in the second stage model (Vansteelandt et al. 2011). To use this control function we set the optional argument ctrl equal to TRUE:

We observe that the control function approach gives a slightly larger estimated effect of vitamin D.

5.3 G-estimation with a point outcome

G-estimation with d(L,Z)=Z−E(Z|L) requires one model for filaggrin, regressed on age, and one model for death, regressed on age, filaggrin and vitamin D. We use logistic regression models for these, fitted as

We now use G-estimation to the fit the causal model (1) with a logit link, by typing

The X argument specifies the names of the exposure; this must be specified when estmethod="g". The fitZ.L and fitY.LZX arguments specify the fitted models for the IV and the outcome, respectively. Finally, the link argument specifies the desired link function η in model (1); either "identity", "log" or "logit". For the identity and log links, G-estimation does not require a model for the outcome, in which case the argument fitY.LZX is not used.

We summarize and compute confidence intervals, by typing

We observe that the estimated effect of vitamin D is similar for G-estimation as for two-stage estimation, with a decrease of 1.42 units in the log odds of death during follow-up, for each 1 unit increase in standardized vitamin D. However, as for two-stage estimation the effect is not statistically significant.

The two stage estimate and the G-estimate are very similar (1.39 vs 1.42), but the standard errors are strikingly different (0.93 vs 1.84). We investigate the possible explanations for this discrepancy in Appendix D, where we do a bootstrap analysis.

Several authors have noted that the estimating eq. (5) does not always have a solution (Vansteelandt et al. 2011; Burgess et al. 2014). The ivglm function prints a warning if no solution is found, but it may also be desirable to manually inspect the estimating function H(ψ) in (5), to make sure that it is equal to 0 at the reported solution. For this purpose we may use the estfun function, which computes H(ψ) for a range of values for ψ around the G-estimate. We type

which gives the plot in Figure 2, where the vertical dashed line indicates the G-estimate. We observe that H(ψ) is indeed equal to 0 at the reported G-estimate −1.423.

Figure 2:

Estimating function for the Vitamin D data, using G-estimation under the causal logistic model (1).

So far, we have considered causal models without interactions, i. e. models with m(L) = 1. For two-stage estimation, interactions are defined implicitly by the fitted outcome model. For G-estimation, interactions have to be defined explicitly, by using an optional argument formula in the ivglm function. This argument specifies m(L) as the right-hand side of a formula on the usual R-format. For instance, to allow for an interaction between vitamin D and age, we type

Here we observe that, for this more complex model, no solution to the estimating equations was found. However, the performance of numerical equation solvers may often be sensitive to specific control parameter settings, such as the choice of starting values and maximum number of iterations. For G-estimation, ivglm and ivcoxph internally use the nleqslv function from the nleqslv package to solve the estimating equations. This function has a wide range of control parameters (see the help page), which can be specified in the call to ivglm and ivcoxph, and are then passed on to nleqslv. For instance, the starting values for the equation solver in nleqslv are specified by an argument x, which defaults to 0 for all elements of ψ . Specifying starting values −7 and 0.1 for the main effect and the interaction term resolves the convergence problem in this particular case:

Figure 3 shows H(ψ) for the main effect and the interaction term, as computed by estfun. When the causal model has more than one parameter, estfun computes H(ψ) over a range of values for each parameter separately, at the G-estimate for the remaining parameters. We observe that H(ψ) is indeed equal to 0 at the reported G-estimates −9.20 and 0.13 for the main effect and the interaction, respectively.

Figure 3:

Estimating function for the Vitamin D data, using G-estimation under the causal logistic model (1), with an interaction term between vitamin D and age.

5.4 Two-stage estimation with a time-to-event outcome

We start by using standard Cox proportional hazards regression to estimate the association between vitamin D and mortality. We fit the model, store the results in an object called fitT.LX, and summarize by typing

We observe that higher vitamin D levels are associated with lower mortality, with a decrease of 0.14 units in the log hazard of death during follow-up, for each 1 unit increase in standardized vitamin D. This association is statistically significant, with a p-value equal to 3.69×10−5.

As before, we may worry that the observed association is partly or fully due to unmeasured confounding. We thus proceed with an IV analysis, fitting the causal model (2) with m(L) = 1. For this purpose we use the ivcoxph function in the ivtools package, which allows for both two-stage estimation and G-estimation.

To carry out two-stage estimation we need one model for vitamin D, regressed on age and filaggrin (the first stage model), and one model for death, regressed on age and vitamin D (the second stage model). For the latter, we use the Cox proportional hazards model fitted above. For the former, we use the same model fitX.LZ as in Section 5.2. We now use two-stage estimation to the fit the causal model (2), by typing

For two-stage estimation, the ivcoxph function has almost identical syntax to the ivglm function. An exception is that the argument specifying the outcome model is called fitT.LX instead of fitY.LX, to emphasize that the outcome is a time-to-event. The exposure model may be an arbitrary GLM, as fitted by the glm function, and the outcome model must be a Cox proportional hazards model, as fitted by the coxph function.

We summarize and compute confidence intervals by typing

We observe that the estimated effect of vitamin D is quite substantial, with a decrease of 1.07 units in the log hazard of death during follow-up, for each 1 unit increase in standardized vitamin D. However, this effect is not statistically significant. We remind the reader that the two-stage estimate of ψ in model (2) is generally biased. The bias can often be reduced by adding the ‘control function’ R=X−Eˆ(X|Z,L) to the set of regressors in the second stage model (Tchetgen Tchetgen et al. 2015). Using the control function approach gives a slightly larger estimate:

The Cox proportional hazards model is by far the most commonly used model for time-to-event outcomes in epidemiologic research. However, the ivtools package also allows for IV analyses with additive hazards models. To fit an additive hazards model we use the ah function in the ivtools package. This function is essentially a wrapper around the ahaz function in the ahaz package, which has a less standard R input/output interface. We fit the model and summarize by typing

Again we observe a protective, statistically significant effect of high vitamin D levels. However, this may be partly or fully due to unmeasured confounding. To carry out two-stage estimation with the fitted additive hazards model, together with the previously fitted model fitX for the exposure, we use the ivah function:

We remind the reader that, since we have used a linear regression model at the first stage, the two-stage estimate of ψ is consistent when using an additive hazards model at the second stage, as opposed to when using a Cox proportional hazards model (Tchetgen Tchetgen et al. 2015).

5.5 G-estimation with a time-to-event outcome

G-estimation with d(L,Z)=Z−E(Z|L) requires one model for filaggrin, regressed on age, and one model for death, regressed on age, filaggrin and vitamin D. For the former, we use the same logistic regression fitZ.L as in Section 5.3. For the latter, we use a Cox proportional hazards model, fitted as

We now use G-estimation to the fit the causal model (2), by typing

For G-estimation, the ivcoxph function has similar syntax as the ivglm function. The two exceptions are that the outcome is called T instead of Y, and the the outcome model is called fitT instead of fitY. We note that T is the name of the follow-up time variable, not the event indicator. We also note that the estimating eq. (6) always requires an estimate of p(T>t|L,Z,X) obtained from an outcome model, and thus fitT.LZX must be specified. This model can be either a Cox proportional hazards model or a set of nonparametric Kaplan-Meier curves, as fitted by the coxph and survfit functions, respectively.

We summarize and compute confidence intervals by typing

We observe that the estimated effect of vitamin D is smaller for G-estimation than for two-stage estimation, with a decrease of 0.86 units in the log hazard of death during follow-up, for each 1 unit increase in standardized vitamin D. This effect is also statistically significant. We inspect the estimating function by typing

which gives the plot in Figure 4, where the vertical dashed line indicates the G-estimate. We observe that H(ψ) is indeed equal to 0 at the reported G-estimate −0.86.

Figure 4:

Estimating function for the Vitamin D data, using G-estimation under the causal Cox proportional hazards model (2).

The summary output tells us that the estimation equation is solved at t = 16.46. By default, ivcoxph solves the equation at the value of t that minimizes the estimated variance of the G-estimate, as suggested by Martinussen et al. (2018). This behavior can be over-ridden by specifying a specific time-point, through an optional argument t. Figure 5 shows the G-estimate (solid line) as a function of t, together with 95 % confidence limits (dashed lines). We observe that, for this particular case, the estimate is not very sensitive to the choice of t; it increases from −1.10 at t = 10 to −0.84 at t = 20. We further observe that the estimate does not change when t > 17; this is because there are no deaths in the data set beyond this point. The estimating equation in (6) only depends on t through pˆ(T>t|Zi,Xi,Li), which is flat beyond the last observed event when obtained from a Cox proportional hazards model.

Figure 5:

G-estimate (solid line) as a function of t, together with 95 % confidence limits (dashed lines).

To carry out G-estimation under the causal additive hazards model (3) we type

The syntax for ivah is identical to the syntax for ivcoxph, except for the additional arguments T and event, which specify the name of the time-to-event variable and event indicator, respectively. To display the G-estimate of B(t) we type

which gives the plot in Figure 6. The solid curve is the estimate of B(t), the dashed curves are point-wise lower and upper 95 % confidence limits. The solid straight line corresponds to a constant exposure effect, dB(t) = ψ. From Figure 6, the assumption of a constant exposure effect seems fairly reasonable. Summarizing the model fit gives:

Figure 6:

Estimate of B(t) for the Vitamin D data, under the causal additive hazards model (3).

This model summary is somewhat different than the other summaries we have presented above, which reflects that the causal model (3) is also somewhat different than the causal models (1) and (2). On the top, two p-values are displayed. The first p-value tests the null hypothesis of no exposure effect, B(t) = 0, and the second p-value tests the null hypothesis of a constant exposure effect, dB(t) = ψ. Both p-values are high, and thus we observe no strong evidence against either of the null hypotheses. On the bottom, the usual summary statistics are given for ψ, assuming a constant effect.