Doubly Robust Estimation with the R Package drgee

2.1 Theory

We assume that data consist of iid observations of (Y,A,L). The simple model [1] contains only main effects. We here consider a more general model which allows for interactions between covariates, and for interactions between the exposure and (some of) the covariates. Towards this end we let X(L) be a p×1-dimensional function of L, and V(L) be a q×1-dimensional function of L. We consider the model

where g is a link function, β is a p×1-dimensional parameter, AX(L) is the element-wise multiplication of A with X(L), and γ is a q-dimensional parameter. Typically, both V(L) and X(L) will contain the element “1”, so that there is an intercept in the model, and a main effect of the exposure on the outcome. Arguing as in Section 1 we can decompose the model in eq. [4] into the main model

and the outcome nuisance model

Our main interest lies in the main model and β is our target parameter. The outcome nuisance model is of secondary interest. However, in O-estimation both models are fitted simultaneously, in such a way that misspecification of the outcome nuisance model may cause bias in the estimate of β. In O-estimation we use the estimating function

When both models [5] and [6] are correct, we have that

It then follows from standard theory for estimating equations (see Appendix) that a consistent and asymptotically normal estimator βˆO can be obtained by solving the estimating equation

for β and γ, where we have used Eˆ(⋅) for sample mean. This estimating equation is identical to the ML score equation for a GLM with canonical link function. It is also identical to the estimating equation for a GEE with “independent working correlation matrix”. The standard error of βˆO can be obtained through the “sandwich formula” (see Appendix).

2.2 Example 1

To demonstrate how the drgee package can be used for O-estimation, we use the dataset SLID from the car package (Fox and Weisberg, 2011). This publically available dataset contains data from the 1994 wave of Canadian Survey of Labour and Income dynamics, with variables wages, sex, education, age and language for composite hourly wage (Canadian dollar), sex, years of education, age and native language respectively. The variables sex and language are factors with levels ("Female", "Male") and ("English", "French", "Other"), respectively. We refer to the car package for a more thorough documentation of the dataset.

Suppose that we wish to use these data study if there is a direct effect of sex on wages, not mediated through education level. Such a direct effect is clearly of substantive interest, since it would be an indication of sex discrimination. To eliminate the mediated effect we wish to control for education level. However, to avoid bias we must then additionally control for covariates that can be confounders for the mediator and the outcome (see Valeri and VanderWeele, 2013 and references therein). It is obvious that age can be such a mediator-outcome confounder, since age is likely to be associated with both education level and wages. Native language may also be a mediator-outcome confounder, by being associated with education level and wages through ethnicity and socio-economic status. We thus use a model for the mean of wages conditional on sex, education, age and language:

where the factors sex and language are recoded as dummy variables, with reference levels "Female" and "English", respectively. In this model, Y=wage, A=sexMale, L=(education,age,languageFrench,languageOther), X(L)=1 and V(L)=(1,L). The target parameter is β and the main model is

The nuisance parameter is γ=(γ0,γ1,γ2,γ3,γ4) and the outcome nuisance model is

To use O-estimation based on these models we type:

> library(car)
> fit <- drgee(oformula=wages~education + age + language,
+ exposure="sex", estimation.method="o", data=SLID)

By setting the estimation.method argument to "o" we tell the drgee function to use O-estimation. The oformula argument specifies the outcome nuisance model and the exposure argument specifies the exposure. There is no need to explicitly specify the outcome, since this is identified by the drgee function as the response variabel in oformula. To summarize the results we type:

> summary(fit)

which gives us the output

Call: drgee(exposure = "sex", oformula = wages ~ education + age + language, data
   = SLID, estimation.method = "o")
Outcome: wages
Exposure: sexMale
Covariates: education,age,languageFrench,languageOther
Main model: wages ~ sexMale
Outcome nuisance model: wages ~ education + age + languageFrench + languageOther
Outcome link function: identity
   Estimate Std. Error z value Pr(>|z|)
sexMale 3.4554 0.2091 16.53 <2e-16 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
(Note: The estimated parameters quantify the conditional exposure-outcome association, given the covariates included in the nuisance models)
3987  complete observations used

Only the result for the main parameters are shown in the output. Since drgee only uses complete observations, 3428 incomplete observations in the original dataset SLID were not used in the calculations. The estimate of β indicates that men earned 3.46 dollar more per hour than women in the target population, even when controlling for education level, age and native language. This observed sex difference is highly significant, with a p-value much smaller than the nominal 0.05-level.

2.3 Example 2

The linear and interaction-free model in Example 1 is simple, but may not be entirely realistic. Wage distributions are often right skewed, and therefore a linear model may not fit data very well. Furthermore, inference for direct effects may be misleading if significant exposure-mediator interactions are omitted from the model (see Valeri and VanderWeele, 2013 and references therein).

Therefore, suppose that we want to use a log-linear model instead, including an interaction between sex and education. We then replace the main model [7] with

in which X(L)= (1,education), and the outcome nuisance model [8] with

To use O-estimation based on these models we type:

> fit <- drgee(oformula=wages~education+age+language,
+ exposure="sex", iaformula=~education, olink="log",
+ estimation.method="o", data=SLID)

The olink argument specifies the link function g. It defaults to "identity" which gives a linear model. By setting the olink argument equal to "log" we tell the drgee function to fit a log-linear model instead. The iaformula argument specifies X(L), i.e. the set of covariates that are assumed to interact with the exposure. Summarizing the results gives:

Call: drgee(exposure = "sex", oformula = wages ~ education + age + language,
  iaformula = ~education, olink = "log", data = SLID, estimation.method = "o")
Outcome: wages
Exposure: sexMale
Covariates: education,age,languageFrench,languageOther
Main model: wages ~ sexMale + sexMale:education
Outcome nuisance model: wages ~ education + age + languageFrench + languageOther
Outcome link function: log
              Estimate Std. Error z value Pr(>|z|)
sexMale           0.581088  0.062848  9.246 < 2e-16 ***
sexMale:education     –0.026175  0.004511  –5.802 6.54e-09 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
(Note: The estimated parameters quantify the conditional exposure-outcome association, given the covariates included in the nuisance models)
3987  complete observations used

We observe that both the main effect of sex and the sex-education interaction are highly significant.

3.1 Theory: E-estimation when g is the identity or log link function

When g is the identity or log link function, we use a model for the mean of the exposure conditional on the covariates:

In this model, h the identity, log or logit link function, α an r×1-dimensional nuisance parameter and Z(L) an r×1-dimensional function of L. We will refer to eq. [11] as the exposure nuisance model. To find an E-estimator of β, we can proceed as follows. Define

When the main model [5] is correct with true parameter β, S(Y,A,L;β) predicts E(Y|A=0,L). If g is the identity link function, then

where the main model [5] was used in the last equality. If g is the log link function, then

where the main model [5] was used in the second equality. Since A−h−1{αTZ(L)} is a constant conditional on A and L and since X(L), E(Y|A=0,L) and h−1{αTZ(L)} are constants conditional on L, we have (by the law of total expectation) that

The last expression equals 0 when the exposure nuisance model [11] is correct. This motivates the estimating function

When both the main model [5] and the exposure nuisance model [11] are correct, we have that

It then follows from standard theory for estimating equations (see Appendix) that a consistent and asymptotically normal estimator βˆE can be obtained by solving the estimating equation

for β and α. This approach is often referred to as “G-estimation” in the causal inference literature, and was initially proposed by Robins et al. (1992). The standard error of βˆE can be obtained through the “sandwich formula” (see Appendix).

3.2 Example 3

Continuing Example 2 (Section 2.3), suppose that we want to combine the main model [9] with the exposure nuisance model

To use E-estimation based on these models we would type:

> fit <- drgee(outcome="wages",
+ eformula=sex~education + age + language,
+ iaformula=~education, olink="log", elink="logit",
+ estimation.method="e", data=SLID)

By setting the estimation.method argument to "e" we tell the drgee function to use E-estimation. We then use an eformula argument to specify the exposure nuisance model. Now there is no need to explicitly specify the exposure, since this is identified by the drgee function as the response variable in eformula. The elink argument specifies the link function h for the exposure nuisance model. Since no outcome nuisance model is used in E-estimation, the oformula argument may be omitted. However, we then need to specify the outcome through the outcome argument. Summarizing the results gives:

Call: drgee(outcome = "wages", eformula = sex ~ education + age
  + language, iaformula = ~education, olink = "log", elink =
  "logit", data = SLID, estimation.method = "e")
Outcome: wages
Exposure: sexMale
Covariates: education,age,languageFrench,languageOther
Main model: wages ~ sexMale + sexMale:education
Outcome link function: log
Exposure nuisance model: sexMale ~ education + age + languageFrench + languageOther
Exposure link function: logit
              Estimate Std. Error z value Pr(>|z|)
sexMale           0.370139  0.064773  5.714  1.1e-08 ***
sexMale:education      –0.010613  0.004738  –2.240  0.0251 *
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
(Note: The estimated parameters quantify the conditional exposure-outcome association, given the covariates included in the nuisance models)
3987  complete observations used

We observe that the obtained E-estimates are quite different from the O-estimates (see Section 2.3). This indicates that at least one of the nuisance models [10] and [13] is misspecified.

3.3 Theory: E-estimation when g is the logit link function

When g is the logit link function and main model [5] holds, E{S(Y,A,L;β)|A,L} is not generally equal to E(Y|A=0,L). This implies that G-estimation (see Section 3.1) does not work. However, when both A and Y are binary we can instead use the symmetry of the odds ratio to perform E-estimation. In this case, the main model [5] is equivalent with the model

Combining this with the exposure nuisance model

where δ is an r×1-dimensional parameter and W(L) is an r×1-dimensional function of L, we arrive at the model

which has the same structural form as the model [4] used in O-estimation when g is the logit link function. Therefore, we can use the estimating function ψO to obtain perform E-estimation if we let A and Y “switch place” in the estimating function. It then follows from standard theory for estimating equations (see Appendix) that a consistent and asymptotically normal estimator βˆE can be obtained by solving the estimating equation

for β and δ. The standard error of βˆE can be obtained through the “sandwich formula” (see Appendix).

3.4 Example 4

To perform E-estimation with logit outcome link function, we need both outcome and exposure to be binary. Suppose that we recode wages in the SLID dataset as a binary variable:

> SLID$highWage<-ifelse(SLID$wages< = 14,0,1)

We can then use the logistic main model

and the logistic exposure nuisance model

To use E-estimation based on these models we type:

> fit <-
drgee(outcome = "highWage", eformula = sex~education + age + language,
+ iaformula = ~education, olink = "logit", elink = "logit",
+ estimation.method = "e", data = SLID)
> summary(fit)
Call: drgee(outcome = "highWage", eformula = sex ~ education + age +
  language, iaformula = ~education, olink = "logit", elink = "logit",
  data = SLID, estimation.method = "e")
Outcome: highWage
Exposure: sexMale
Covariates: education,age,languageFrench,languageOther
Main model: highWage ~ sexMale + sexMale:education
Outcome link function: logit
Exposure nuisance model: sexMale ~ education + age + languageFrench + languageOther
Exposure link function: logit
              Estimate Std. Error z value Pr(>|z|)
sexMale           1.93520  0.30980  6.247  4.2e-10 ***
sexMale:education      –0.06361  0.02287  –2.782  0.00541 **
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
(Note: The estimated parameters quantify the conditional exposure-outcome association, given the covariates included in the nuisance models)
3987  complete observations used

Again, both main effect and interaction are highly significant.

4.1 Theory: DR-estimation when g is the identity or log link function

Let γTV(L) be a model for g{E(Y|A=0,L)} and αTZ(L) a model for h{E(A|L)} as in eqs [6] and [11], respectively. Also, let S(Y,A,L;β) be as in eq. [12]. When the main model [5] is correct and when g is the identity or log link function, the results in Section 3 imply that

We then have that

The last expression equals 0 when at least one of the nuisance models [6] and [11] is correct. This motivates the estimating function

When the main model [5] is correct and either the outcome nuisance model [6] or the exposure nuisance model [14] is correct, we have that

It then follows from standard theory for estimating equations (see Appendix) that a consistent and asymptotically normal estimator βˆDR can be obtained by solving the equation

for (βT,β†T,γT,αT)T. The standard error of βˆDR can be obtained through the “sandwich formula” (see Appendix).

4.2 Example 5

Continuing Example 2 (Section 2.3) and Example 3 (Section 3.2), suppose that we want to combine the main model [9] with the outcome nuisance model [10] and the exposure nuisance model [13]. To use DR-estimation based on these models we type:

> fit <- drgee(oformula=wages~education+age+language,
+ eformula=sex~education+age+language,
+ iaformula=~education, olink="log", elink="logit",
+ estimation.method="dr", data=SLID)

By setting the estimation.method argument to "dr" we tell the drgee function to use DR-estimation. The outcome and exposure arguments may be omitted, since the drgee function identifies the outcome and exposure as the response variables in oformula and eformula, respectively. Summarizing the results gives:

> summary(fit)
Call: drgee(oformula = wages ~ education + age + language,
  eformula = sex ~ education + age + language,
  iaformula = ~education, olink = "log",
  elink = "logit", data = SLID, estimation.method = "dr")
Outcome: wages
Exposure: sexMale
Covariates: education,age,languageFrench,languageOther
Main model: wages ~ sexMale + sexMale:education
Outcome nuisance model: wages ~ education + age + languageFrench + languageOther
Outcome link function: log
Exposure nuisance model: sexMale ~ education + age + languageFrench + languageOther
Exposure link function: logit
              Estimate Std. Error z value Pr(>|z|)
sexMale           0.57752  0.06333  9.119  < 2e-16 ***
sexMale:education      –0.02591  0.00455  –5.696  1.23e-08 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
(Note: The estimated parameters quantify the conditional exposure-outcome association, given the covariates included in the nuisance models)
3987  complete observations used

We observe that the DR estimates are very similar to the O-estimates in Example 2 (Section 2.3), but quite different from the E-estimates in Example 3 (Section 3.2). This indicates that the exposure nuisance model may not be well specified. However, it does not prove that the outcome nuisance model is correct; we demonstrate in Section 6.1 that all methods may agree well even though both nuisance models are misspecified.

4.3 Theory: DR-estimation when g is the logit link function

When g is the logit link function, DR-estimation can be performed by the drgee package if the exposure A is binary and h is the logit link function. In this case the main model [5] is equivalent to the odds ratio model

By combining this model with the outcome nuisance model

and the exposure nuisance model [14] we can construct the estimating function

where

when the main model [18] is correct and either the outcome nuisance model [19] or the exposure nuisance model [14] is correct. It then follows from standard theory for estimating equations (see Appendix) that a consistent and asymptotically normal estimator βˆDR can be obtained by solving the equation

for (βT,β†T,β‡T,γT,δT)T. Tchetgen Tchetgen et al. (2010) showed that this estimator attains the semiparametric efficiency bound when all three models are correct, i.e. under the intersection of models [14], [18] and [19]. We note that this intersection model is not fully parametric, since it does not parametrize the distribution of L.

4.4 Example 6

Continuing Example 4 (Section 3.4), suppose that we want to combine main model [15], the exposure nuisance model [16] and the logistic outcome nuisance model

To use DR-estimation based on these models we type:

> fit <- drgee(oformula=highWage~education+age+language,
+ eformula=sex~education+age+language, iaformula=~education,
+ olink="logit", elink="logit", estimation.method="dr", data=SLID)
> summary(fit)
Call: drgee(oformula = highWage ~ education + age + language,
  eformula = sex ~ education + age + language,
  iaformula = ~education, olink = “logit", elink = “logit",
  data = SLID, estimation.method = “dr")
Outcome: highWage
Exposure: sexMale
Covariates: education,age,languageFrench,languageOther
Main model: highWage ~ sexMale + sexMale:education
Outcome nuisance model: highWage ~ education + age + languageFrench + languageOther
Outcome link function: logit
Exposure nuisance model: sexMale ~ education + age + languageFrench +
                languageOther
Exposure link function: logit
               Estimate Std. Error z value Pr(>|z|)
sexMale            2.9050  0.4015  7.236  4.62e-13 ***
sexMale:education       –0.1341  0.0295  –4.547  5.44e-06 ***
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
(Note: The estimated parameters quantify the conditional exposure-outcome association, given the covariates included in the nuisance models)
3987  complete observations used

We observe that the DR estimates are quite different from the E-estimates in Example 4 (Section 3.4). This indicates that the exposure nuisance model [16] may not be well specified.

5.1 Theory

When data are clustered, the parameter estimates obtained in Sections 2–4 are still consistent. However, their standard errors must be corrected for within-cluster correlations. Suppose that data contain m independent clusters with ni observations Oi1,...,Oini in cluster i for i=1,...,m. Corrected standard errors can then be obtained by replacing the generic estimating function ψ(Oij,θ) in the sandwich formula with ∑j=1niψ(Oij,θ) for i=1,...,ni. This modification corrects the standard errors, but does not affect the parameter estimates.

5.2 Example 7

To demonstrate estimation with clustered data, we use the dataset ohio in geepack package (Højsgaard et al., 2006). This dataset contains data to study the health effects of air pollution on children. The children were examined annually at ages 7–10. Four variables are included in the dataset: resp, smoke, age and id. The variables resp (binary), smoke (binary) and age (continuous) indicate wheezing status, maternal smoking and age, respectively, at each examination. id is the child’s identification number. Suppose that we are interested in the association between maternal smoking and wheezing status conditional on age. Since all subjects are from the same city, we assume that the effect of air pollution is similar for all subjects in the sample. Assuming the main model

the outcome nuisance model

and the exposure nuisance model

we can perform DR-estimation with cluster-corrected standard errors by setting the argument clusterid to "id":

> library(geepack)
> data(ohio)
> fit <- drgee(oformula=resp~age,eformula=smoke~age,
+ olink="logit",elink="logit", estimation.method="dr",
+ data=ohio, clusterid="id")
> summary(fit)
Call: drgee(oformula = resp ~ age, eformula = smoke ~ age,
  olink = "logit", elink ="logit", data = ohio,
  estimation.method = "dr", clusterid="id")
Outcome: resp
Exposure: smoke
Covariates: age
Main model: resp ~ smoke
Outcome nuisance model: resp ~ age
Outcome link function: logit
Exposure nuisance model: smoke ~ age
Exposure link function: logit
     Estimate Std. Error z value Pr(>|z|)
smoke   0.2721  0.1781  1.528  0.127
(Note: The estimated parameters quantify the conditional exposure-outcome association, given the covariates included in the nuisance models)
 2148  complete observations used
Cluster-robust Std. errors
using 537 clusters defined by levels of id

6.1 Simulation study 1

In the first simulation we generated data from the following model:

Under this model, we generated 1,000 samples with 500 observations each. Each sample was analyzed with O-estimation, E-estimation and DR-estimation. For each estimation method we carried out four analyses. In the first analysis, β was estimated with both nuisance models correctly specified. In the second analysis we used the misspecified outcome nuisance model

In the third analysis we used the misspecified exposure nuisance model

In the fourth analysis we used both misspecified nuisance models [20] and [21]. For each estimation method, we calculated the mean (over the 1,000 samples) estimate of β, the mean standard error, and the empirical coverage probability of a 95% Wald confidence interval for β. For comparison we also calculated the empirical standard error, i.e. the standard deviation of the estimate over the 1,000 samples. The results are shown in Table 1.

In the first analysis, all three mean estimates are close to the true value 1.5. This demonstrates that all three methods give unbiased estimates of β when all modelling assumptions are correct. The standard error is larger for the E-estimator than for the O-estimator. This is not a coincidence; Robins et al. (1992) showed that the O-estimator is always at least as efficient as the E-estimator. Note that the standard error of the DR estimator is close to the standard error of the O-estimator. In the second analysis, the assumed outcome nuisance model is misspecified, leading to bias in the O-estimator. The E-estimator does not use the outcome nuisance model and is therefore unaffected by this misspecification. Despite using a misspecified model for the outcome nuisance, the DR estimator is unbiased. Furthermore, its standard error is nearly identical to the standard error of the E-estimator. In the third analysis, the assumed exposure nuisance model is misspecified, leading to bias in the E-estimator. The DR estimator is unbiased for β and its standard error is nearly identical to the standard error of the O-estimator. In the fourth analysis, both nuisance models are misspecified, leading to bias for all three methods. In this analysis, both the mean estimates and standard errors are very similar for all three methods. This demonstrates that agreement between the three methods does not imply that both nuisance models are correctly specified. We finally observe that the mean standard errors agree well with the empirical standard errors for all methods and analyzes, and that the 95% confidence interval has close to 95% coverage probability as long as the corresponding estimator is unbiased.

6.2 Simulation study 2

In the second simulation we generated data from the following model:

Under this model, we generated 1,000 samples with 500 observations each. Each sample was analyzed with O-estimation, E-estimation and DR-estimation. For each estimation method we carried out four analyses. In the first analysis, both outcome and exposure nuisance models were correct. In the second analysis, we used the misspecified outcome nuisance model

In the third analysis, we used the misspecified exposure nuisance model

In the fourth analysis we used both misspecified nuisance models [22] and [23]. We calculated the same summary statistics as in the first simulation. The results are shown in Table 2.

The results in Table 2 again demonstrate that the DR estimator is indeed doubly robust, i.e. that it is unbiased as long as at least one of the nuisance models is correctly specified. For the main effect β1, the standard error of the DR estimator is very similar to the standard errors of the O- and E-estimators. However, for the interaction β2, the standard errors of the DR estimator is generally larger than the standard errors of the O- and E-estimators.