# Simulation from Kang and Shafer (2007) and Imai and Ratkovic (2014)
# ATE estimation
# True ATE is 10
tau <- 10
set.seed(12345)
n <- 1000
X <- matrix(stats::rnorm(n * 4, mean = 0, sd = 1), nrow = n, ncol = 4)
prop <- 1 / (1 + exp(X[, 1] - 0.5 * X[, 2] +
0.25 * X[, 3] + 0.1 * X[, 4]))
treat <- rbinom(n, 1, prop)
y <- 210 +
27.4 * X[, 1] + 13.7 * X[, 2] + 13.7 * X[, 3] + 13.7 * X[, 4] +
tau * treat + stats::rnorm(n = n, mean = 0, sd = 1)
ybinom <- (y > 210) + 0
df0 <- data.frame(X, treat, y, ybinom)
colnames(df0) <- c("x1", "x2", "x3", "x4", "treat", "y", "ybinom")
# Variables for a misspecified model
Xmis <- data.frame(x1mis = exp(X[, 1] / 2),
x2mis = X[, 2] * (1 + exp(X[, 1]))^(-1) + 10,
x3mis = (X[, 1] * X[, 3] / 25 + 0.6)^3,
x4mis = (X[, 2] + X[, 4] + 20)^2)
# Data frame and formulas for propensity score estimation
df <- data.frame(df0, Xmis)
formula_ps_c <- stats::as.formula(treat ~ x1 + x2 + x3 + x4)
formula_ps_m <- stats::as.formula(treat ~ x1mis + x2mis +
x3mis + x4mis)
# Formula for a misspecified outcome model
formula_y <- stats::as.formula(y ~ x1mis + x2mis + x3mis + x4mis)
# Set plot margins
oldpar <- par(no.readonly = TRUE) # Just for adjusting plot margins
par(mar = c(5.1, 5.1, 4.1, 2.1)) # Just for adjusting plot margins
# Misspecified propensity score model
# Distribution balancing weighting without regularization
fitdbwm <- dbw(formula_y = formula_y, formula_ps = formula_ps_m,
estimand = "ATE", method = "dbw",
method_y = "wls", data = df, normalize = TRUE,
vcov = TRUE, lambda = 0, weights = NULL,
clevel = 0.95)
plot(fitdbwm, addcov = ~ x1 + x2 + x3 + x4, threshold = c(0.1, 0.2))
# Non-absolute values
plot(fitdbwm, addcov = ~ x1 + x2 + x3 + x4, absolute = FALSE,
threshold = c(0.1, 0.2))
# No sorting
plot(fitdbwm, addcov = ~ x1 + x2 + x3 + x4,
threshold = c(0.1, 0.2), sort = FALSE)
# Covariate balancing weighting function without regularization
fitcbwm <- dbw(formula_y = formula_y, formula_ps = formula_ps_m,
estimand = "ATE", method = "cb",
method_y = "wls", data = df, normalize = TRUE,
vcov = TRUE, lambda = 0, weights = NULL,
clevel = 0.95)
plot(fitcbwm, addcov = ~ x1 + x2 + x3 + x4, threshold = c(0.1, 0.2))
# Entropy balancing weighting function without regularization
fitebwm <- dbw(formula_y = formula_y, formula_ps = formula_ps_m,
estimand = "ATE", method = "eb",
method_y = "wls", data = df, normalize = TRUE,
vcov = TRUE, lambda = 0, weights = NULL,
clevel = 0.95)
plot(fitebwm, addcov = ~ x1 + x2 + x3 + x4, threshold = c(0.1, 0.2))
# Standard logistic regression
fitmlem <- dbw(formula_y = formula_y, formula_ps = formula_ps_m,
estimand = "ATE", method = "mle",
method_y = "wls", data = df, normalize = FALSE,
vcov = TRUE, lambda = 0, weights = NULL,
clevel = 0.95)
plot(fitmlem, addcov = ~ x1 + x2 + x3 + x4, threshold = c(0.1, 0.2))
par(oldpar)
# Display the covariate balance matrix
cb <- plot(fitdbwm, addcov = ~ x1 + x2 + x3 + x4, plot = FALSE)
cb
Run the code above in your browser using DataLab