## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 4.2, dpi = 110,
                      out.width = "100%")

## -----------------------------------------------------------------------------
library(drlate)
data(drlate_sim)

fit <- drlate(lwage ~ age + educ,      # outcome model
              nvstat ~ age + educ,     # treatment model
              rsncode ~ age + educ,    # instrument propensity score model
              data = drlate_sim)
summary(fit)

## -----------------------------------------------------------------------------
coef(fit)
confint(fit)

## -----------------------------------------------------------------------------
# AIPW with unnormalized moments
drlate(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
       data = drlate_sim, method = "aipw", normalized = FALSE)

# IPW: no covariates in the outcome/treatment equations
drlate(lwage ~ 1, nvstat ~ 1, rsncode ~ age + educ,
       data = drlate_sim, method = "ipw")

# Regression adjustment: no instrument covariates
drlate(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ 1,
       data = drlate_sim, method = "ra")

## -----------------------------------------------------------------------------
# Normalized Abadie kappa (kappalate tau_a,10); reports the LATE only,
# since the estimator is a difference of two ratios
drlate(lwage ~ 1, nvstat ~ 1, rsncode ~ age + educ,
       data = drlate_sim, method = "kappa10")

# Unnormalized Abadie kappa (tau_a); Fieller sets available
fit_k <- drlate(lwage ~ 1, nvstat ~ 1, rsncode ~ age + educ,
                data = drlate_sim, method = "kappa")
confint(fit_k, method = "fieller")

## -----------------------------------------------------------------------------
# LATT with an inverse-probability-tilted instrument propensity score
drlate(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
       data = drlate_sim, estimand = "latt", ivmodel = "ipt")

# Poisson outcome model for the positive wage level
drlate(kwage ~ age + educ, nvstat ~ age + educ, rsncode ~ 1,
       data = drlate_sim, method = "ra", omodel = "poisson")

## -----------------------------------------------------------------------------
drlate(lwage ~ age, nvstat ~ age, rsncode ~ age, data = drlate_sim,
       cluster = drlate_sim$educ)

## ----diag-balance-------------------------------------------------------------
fit <- drlate(lwage ~ age + educ, nvstat ~ age + educ,
              rsncode ~ age + educ, data = drlate_sim)
plot(fit, type = "balance")
balance(fit)

## ----diag-overlap-------------------------------------------------------------
plot(fit, type = "overlap")

## ----diag-postest-------------------------------------------------------------
complier_means(fit)
balance_test(fit)

## -----------------------------------------------------------------------------
# Weak-instrument-robust Fieller confidence set (may be unbounded when
# the first stage is weak -- that is the honest answer)
confint(fit, method = "fieller")

# Nonparametric bootstrap (percentile CIs; clusters resampled whole
# when `cluster` is supplied)
fit_b <- drlate(lwage ~ age + educ, nvstat ~ age + educ,
                rsncode ~ age + educ, data = drlate_sim,
                vcov = "bootstrap", boot_reps = 199, boot_seed = 1)
confint(fit_b)

## -----------------------------------------------------------------------------
d_os <- drlate_sim
d_os$nvstat[d_os$rsncode == 0] <- 0L
dr_hausman(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
           data = d_os)

## ----compare-plot-------------------------------------------------------------
cmp <- drlate_compare(lwage ~ age + educ, nvstat ~ age + educ,
                      rsncode ~ age + educ, data = drlate_sim)
cmp
plot(cmp)

## ----eval = FALSE-------------------------------------------------------------
# sipp <- haven::read_dta("https://people.brandeis.edu/~tslocz/sipp.dta")
# sipp <- subset(as.data.frame(sipp),
#                !is.na(kwage) & !is.na(educ) & rsncode != 999)
# sipp$lwage <- log(sipp$kwage)
# 
# # Stata: drlate (lwage age_5) (nvstat age_5) (rsncode age_5)
# drlate(lwage ~ age_5, nvstat ~ age_5, rsncode ~ age_5, data = sipp)
# 
# # Stata: drlate (lwage age_5) (nvstat age_5) (rsncode age_5, ipt), latt
# drlate(lwage ~ age_5, nvstat ~ age_5, rsncode ~ age_5, data = sipp,
#        ivmodel = "ipt", estimand = "latt")
# 
# # Stata: kappalate lwage (nvstat = rsncode) age_5, zmodel(logit) which(all)
# drlate(lwage ~ 1, nvstat ~ 1, rsncode ~ age_5, data = sipp,
#        method = "kappa")     # tau_a; likewise "kappa0", "kappa10",
#                              # and method = "ipw" for tau_u / tau_a,1
# 
# # Stata: kappalate lwage (nvstat = rsncode) age_5, zmodel(probit)
# drlate(lwage ~ 1, nvstat ~ 1, rsncode ~ age_5, data = sipp,
#        method = "kappa", ivmodel = "probit")