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

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

## -----------------------------------------------------------------------------
coef(lm(lwage ~ nvstat + age + educ, drlate_sim))["nvstat"]

## -----------------------------------------------------------------------------
with(drlate_sim,
     (mean(lwage[rsncode == 1]) - mean(lwage[rsncode == 0])) /
     (mean(nvstat[rsncode == 1]) - mean(nvstat[rsncode == 0])))

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

## -----------------------------------------------------------------------------
# Propensity score misspecified; regressions correct
coef(drlate(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ 1,
            data = drlate_sim))[1]

# Regressions misspecified (intercept only); propensity score correct
coef(drlate(lwage ~ 1, nvstat ~ 1, rsncode ~ age + educ,
            data = drlate_sim))[1]

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

## ----balance------------------------------------------------------------------
plot(fit, type = "balance")

## -----------------------------------------------------------------------------
balance(fit)

## ----balance-density, fig.height = 3------------------------------------------
plot(fit, type = "balance_density", var = "age")

## ----weights------------------------------------------------------------------
plot(fit, type = "weights")

## ----balance-test-------------------------------------------------------------
balance_test(fit)

## ----compliers----------------------------------------------------------------
complier_means(fit)

## -----------------------------------------------------------------------------
# Binary outcome: logit keeps fitted probabilities in [0, 1]
coef(drlate(hijob ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
            data = drlate_sim, omodel = "logit"))

# Positive outcome: Poisson (quasi-likelihood; no distributional claim)
coef(drlate(kwage ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
            data = drlate_sim, omodel = "poisson"))

## -----------------------------------------------------------------------------
coef(drlate(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
            data = drlate_sim, ivmodel = "ipt"))[1]

# probit propensity score with a weighting estimator (Section 6)
coef(drlate(lwage ~ 1, nvstat ~ 1, rsncode ~ age + educ,
            data = drlate_sim, method = "kappa10",
            ivmodel = "probit"))[1]

## -----------------------------------------------------------------------------
coef(drlate(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
            data = drlate_sim, estimand = "latt"))

## ----kappa-menu---------------------------------------------------------------
cmp_w <- drlate_compare(lwage ~ 1, nvstat ~ 1, rsncode ~ age + educ,
                        data = drlate_sim,
                        methods = c("ipw", "kappa", "kappa0", "kappa10"))
cmp_w

## ----kappa-shift--------------------------------------------------------------
d_shift <- transform(drlate_sim, lwage = lwage + 100)
rbind(
  kappa10 = c(original = coef(drlate(lwage ~ 1, nvstat ~ 1,
                                     rsncode ~ age + educ,
                                     data = drlate_sim,
                                     method = "kappa10"))[[1]],
              shifted = coef(drlate(lwage ~ 1, nvstat ~ 1,
                                    rsncode ~ age + educ,
                                    data = d_shift,
                                    method = "kappa10"))[[1]]),
  kappa   = c(original = coef(drlate(lwage ~ 1, nvstat ~ 1,
                                     rsncode ~ age + educ,
                                     data = drlate_sim,
                                     method = "kappa"))[[1]],
              shifted = coef(drlate(lwage ~ 1, nvstat ~ 1,
                                    rsncode ~ age + educ,
                                    data = d_shift,
                                    method = "kappa"))[[1]])
)

## -----------------------------------------------------------------------------
set.seed(42)
d_weak <- drlate_sim[1:300, ]
d_weak$zweak <- sample(d_weak$rsncode)
fit_weak <- drlate(lwage ~ age, nvstat ~ age, zweak ~ age, data = d_weak)
print(fit_weak)

## -----------------------------------------------------------------------------
confint(fit, method = "fieller")     # strong instrument: ~ Wald interval

## -----------------------------------------------------------------------------
fit_boot <- drlate(lwage ~ age + educ, nvstat ~ age + educ,
                   rsncode ~ age + educ, data = drlate_sim,
                   vcov = "bootstrap", boot_reps = 199, boot_seed = 42)
summary(fit_boot)

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

## -----------------------------------------------------------------------------
d_os <- drlate_sim
d_os$nvstat[d_os$rsncode == 0] <- 0L   # impose one-sided noncompliance

dr_hausman(lwage ~ age + educ, nvstat ~ age + educ, rsncode ~ age + educ,
           data = d_os)