knitr::opts_chunk$set(echo = TRUE)
causalsandwich
for estimating causal effects with empirical sandwich variance estimatorCurrently implemented as of version 0.0.0.9000: - Estimators of the average treament effect - IPTW - G-formula - Doubly Robust estimator
causalsandwich
? Easy implementation with valid inference!geex
?geex is a package designed for easy implementation of estimating equations. The causalsandwich
package is powered by geex
The package can be installed from Github:
devtools::install_github("BarkleyBG/causalsandwich")
Whenver possible, the causal effect estimate is listed last.
First, some data:
library(causalsandwich) n <- 200 data <- data.frame( Covar1 = rnorm(n), Covar2 = rnorm(n) ) trtprobs <- plogis(0.2 + 0.5*data$Covar1 + 0.1*data$Covar2 + 0.2*data$Covar1 * data$Covar2) data$BinaryTrt <- rbinom(n, 1, trtprobs) outprobs <- plogis(0.5 + 0.2*data$Covar1 - 0.2*data$Covar2 + 0.2*data$Covar1 * data$Covar2 -0.2*data$BinaryTrt) data$BinaryOutcome <- rbinom(n, 1, outprobs)
Until rigorous tests and defensive programming are added, users are recommended to use the following coding schemes:
Estimating the Average Treatment Effect only takes specifying a few arguments. The below code will estimate propensity scores with the logistic GLM model and estimate the ATE with the Horwitz-Thompson IPTW estimator.
The formula specified is a mixture of two formulas. That is, BinaryTrt ~ Covar1*Covar2
will be passed into the logistic model to fit a model for treatment with main effects and interaction terms with the two covariates. The first section of the multi-part formula indicates that the column named BinaryOutcome
in data
corresponds to the outcome variable of interest.
formula <- BinaryOutcome | BinaryTrt ~ Covar1 * Covar2 IPTW <- estimateIPTW( data = data, formula = formula )
Output includes point estimates and the sandwich variance estimator:
(ests <- IPTW@estimates) (vcov <- IPTW@vcov)
The first four parameters correspond to the treatment model parameters. It is the fifth parameter that corresponds to the ATE. Thus, a point estimate and 95\% confidence interval for this parameter are:
param_num <- 5 ests[param_num] ests[param_num] + stats::qnorm(c(0.025,0.975))*vcov[param_num,param_num]
Since the treatment model is correctly specified, the estimator is consistent (as is the variance estimator).
outcome_regression_formula
to the formula
argument.treatment_var_name = "BinaryTrt"
.outcome_regression_formula <- BinaryOutcome ~ BinaryTrt + Covar1 * Covar2 GF <- estimateGF( data = data, treatment_var_name = "BinaryTrt", formula = outcome_regression_formula, model_method = "logistic" ) (ests <- GF@estimates) (vcov <- GF@vcov)
An estimate and 95\% confidence interval is then
param_num <- 6 ests[param_num] ests[param_num] + stats::qnorm(c(0.025,0.975))*vcov[param_num,param_num]
deriv_control = geex::setup_deriv_control(method="simple")
. This results in quicker (but probably less accurate) derivatives in the variance computations. This is passed directly to the geex
guts.outcome_regression_formula <- BinaryOutcome ~ BinaryTrt + Covar1 * Covar2 treatment_model_formula <- BinaryTrt ~ Covar1 * Covar2 DRIPTW <- estimateDRIPTW( data = data, outcome_formula = outcome_regression_formula, treatment_formula = treatment_model_formula, outcome_model_method = "logistic", treatment_model_method = "logistic", deriv_control = geex::setup_deriv_control(method="simple") ) (ests <- DRIPTW@estimates) (vcov <- DRIPTW@vcov)
An estimate and 95\% confidence interval is then
param_num <- 10 ests[param_num] ests[param_num] + stats::qnorm(c(0.025,0.975))*vcov[param_num,param_num]
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.