library(tidyverse)
library(broom)
library(huxtable)
library(estimatr)
library(fabricatr)
set.seed(1234)
dat <- fabricate(
N = 40,
Y = rpois(N, lambda = 4),
Z = rbinom(N, 1, prob = 0.4),
D = Z * rbinom(N, 1, prob = 0.8),
X = rnorm(N),
G = sample(letters[1:4], N, replace = TRUE)
)
# Manual way
model_first <- lm(D ~ Z + X, data = dat)
dat_updated <- augment_columns(model_first, dat) %>%
rename(D_hat = .fitted)
model_second <- lm(Y ~ D_hat + X, data = dat_updated)
# Automatic way
model_automatic <- iv_robust(Y ~ D + X | Z + X, data = dat)
huxreg(list("First stage" = model_first,
"Second stage" = model_second,
"Second stage (auto)" = model_automatic))──────────────────────────────────────────────────────────────────
First stage Second stage Second stage (auto)
────────────────────────────────────────────────────
(Intercept) 0.000 3.750 *** 3.750 ***
(0.030) (0.366) (0.385)
Z 0.917 ***
(0.054)
X -0.001 0.644 * 0.644 *
(0.026) (0.309) (0.281)
D_hat -0.482
(0.717)
D -0.482
(0.597)
────────────────────────────────────────────────────
N 40 40 40
R2 0.885 0.113 0.114
logLik 18.761 -80.932
AIC -29.522 169.864
──────────────────────────────────────────────────────────────────
*** p < 0.001; ** p < 0.01; * p < 0.05.
