## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 6,
  fig.height = 4,
  dpi = 90,
  message = FALSE,
  warning = FALSE
)

## ----setup--------------------------------------------------------------------
library(rdhte)
library(rdrobust)
library(broom)
data(rdhte_dataset)
attach(rdhte_dataset)

## ----binary-------------------------------------------------------------------
summary(rd_left <- rdhte(y = y, x = x, covs.hte = factor(w_left),
                          cluster = cluster_var))

## ----binary-lincom------------------------------------------------------------
rdhte_lincom(model = rd_left,
             linfct = c("`factor(w_left)1` - `factor(w_left)0` = 0"))

## ----binary-bwjoint-----------------------------------------------------------
summary(rdhte(y = y, x = x, covs.hte = w_left, cluster = cluster_var,
              bw.joint = TRUE))

## ----binary-bare--------------------------------------------------------------
summary(rd_left2 <- rdhte(y = y, x = x, covs.hte = w_left,
                           cluster = cluster_var))
rdhte_lincom(model = rd_left2,
             linfct = c("`w_left1` - `w_left0` = 0"))

## ----categorical-unordered----------------------------------------------------
summary(rd_ideology <- rdhte(y = y, x = x,
                              covs.hte = factor(w_ideology),
                              cluster = cluster_var))

## ----categorical-lincom-------------------------------------------------------
rdhte_lincom(model = rd_ideology,
             linfct = c("`factor(w_ideology)2` = 0",
                        "`factor(w_ideology)3` = 0",
                        "`factor(w_ideology)4` = 0"))

## ----categorical-ordered------------------------------------------------------
summary(rdhte(y = y, x = x, covs.hte = factor(w_strength_qrt),
              cluster = cluster_var))

## ----categorical-as-continuous------------------------------------------------
summary(rdhte(y = y, x = x, covs.hte = w_strength_qrt,
              cluster = cluster_var))

## ----factor-by-factor---------------------------------------------------------
summary(rdhte(y = y, x = x,
              covs.hte = factor(w_left):factor(w_strong),
              cluster = cluster_var))

## ----factor-by-factor-alt-----------------------------------------------------
interactions <- 1*(w_left==0)*(w_strong==1) +
                2*(w_left==0)*(w_strong==2) +
                3*(w_left==1)*(w_strong==1) +
                4*(w_left==1)*(w_strong==2)
summary(rdhte(y = y, x = x, covs.hte = factor(interactions),
              cluster = cluster_var))

## ----continuous---------------------------------------------------------------
summary(rd_continuous <- rdhte(y = y, x = x, covs.hte = w_strength,
                                kernel = "uni",
                                cluster = cluster_var))

## ----continuous-vs-lm---------------------------------------------------------
trt <- (x > 0)
new.data <- data.frame(y, x, w_strength, trt)
using.lm <- coef(lm(y ~ trt * x * w_strength,
                    data = new.data[abs(new.data$x) < rd_continuous$h[1], ]))
rd_continuous$Estimate[1] - using.lm["trtTRUE"]
rd_continuous$Estimate[2] - using.lm["trtTRUE:w_strength"]

## ----binary-by-continuous-----------------------------------------------------
summary(rd_interaction <- rdhte(y = y, x = x,
                                 covs.hte = "w_left*w_strength",
                                 cluster = cluster_var))

## ----binary-by-continuous-factor----------------------------------------------
summary(rdhte(y = y, x = x, covs.hte = "factor(w_left)*w_strength",
              cluster = cluster_var))

## ----binary-by-continuous-joint-----------------------------------------------
rdhte_lincom(model = rd_interaction,
             linfct = c("`T` = 0",
                        "`T:w_left` = 0",
                        "`T:w_strength` = 0",
                        "`T:w_left:w_strength` = 0"))

## ----binary-by-continuous-fixedbw---------------------------------------------
summary(rd_interaction <- rdhte(y = y, x = x,
                                 covs.hte = "w_left*w_strength",
                                 h = 0.1, cluster = cluster_var))
summary(rdhte(y = y[w_left == 0], x = x[w_left == 0],
              covs.hte = w_strength[w_left == 0], h = 0.1,
              cluster = cluster_var[w_left == 0]))
summary(rdhte(y = y[w_left == 1], x = x[w_left == 1],
              covs.hte = w_strength[w_left == 1], h = 0.1,
              cluster = cluster_var[w_left == 1]))
rdhte_lincom(model = rd_interaction,
             linfct = c("`T` + `T:w_left` = 0",
                        "`T:w_strength` + `T:w_left:w_strength` = 0"))

## ----rdbwhte------------------------------------------------------------------
bw <- rdbwhte(y = y, x = x, covs.hte = factor(w_ideology),
              cluster = cluster_var)
bw$h

## ----rdbwhte-joint------------------------------------------------------------
bw_joint <- rdbwhte(y = y, x = x, covs.hte = factor(w_ideology),
                    cluster = cluster_var, bw.joint = TRUE)
bw_joint$h

## ----covs-eff-----------------------------------------------------------------
m_no_eff <- rdhte(y = y, x = x, covs.hte = factor(w_ideology),
                   cluster = cluster_var)
m_eff    <- rdhte(y = y, x = x, covs.hte = factor(w_ideology),
                   covs.eff = w_strength, cluster = cluster_var)
data.frame(cell       = m_no_eff$W.lev,
           se_no_eff  = round(m_no_eff$se.rb, 4),
           se_eff     = round(m_eff$se.rb,    4),
           pct_change = round(100 * (m_eff$se.rb - m_no_eff$se.rb)
                              / m_no_eff$se.rb, 1))

## ----plot-basic, fig.alt="rdhte forest plot, ideology buckets"----------------
plot(rd_ideology)

## ----plot-sort, fig.alt="sorted forest plot"----------------------------------
plot(rd_ideology, sort = TRUE)

## ----plot-custom, fig.alt="custom-themed forest plot"-------------------------
p <- plot(rd_ideology, sort = TRUE,
          title = "Heterogeneity by ideology bucket",
          ylab  = "Sharp RD ITT")
if (requireNamespace("ggplot2", quietly = TRUE)) {
  p + ggplot2::theme_minimal(base_size = 11) +
      ggplot2::coord_flip()
}

## ----table-tidy---------------------------------------------------------------
tab_A <- tidy(rd_ideology)
tab_A

## ----table-A-kable------------------------------------------------------------
knitr::kable(
  tab_A[, c("term", "estimate", "std.error", "conf.low", "conf.high",
            "n.eff.left", "n.eff.right", "h.left", "h.right")],
  digits = c(NA, 3, 3, 3, 3, 0, 0, 3, 3),
  col.names = c("Cell", "Estimate", "SE", "CI lo", "CI hi",
                "Nh-", "Nh+", "h-", "h+"),
  caption = "Per-cell RD effects by ideology bucket"
)

## ----table-glance-------------------------------------------------------------
glance(rd_ideology)

## ----table-spec-compare-------------------------------------------------------
specs <- list(HC0 = list(vce = "hc0"),
              HC2 = list(vce = "hc2"),
              HC3 = list(vce = "hc3"),
              CR1 = list(vce = "cr1", cluster = cluster_var))

fit_one <- function(args) {
  args <- c(list(y = y, x = x, covs.hte = factor(w_ideology)), args)
  do.call(rdhte, args)
}

cells  <- sort(unique(w_ideology))
mat_pt <- sapply(specs, function(s) tidy(fit_one(s))$estimate)
mat_se <- sapply(specs, function(s) tidy(fit_one(s))$std.error)
rownames(mat_pt) <- rownames(mat_se) <- as.character(cells)

mat_pt
mat_se

## ----table-spec-stacked-------------------------------------------------------
n_cells <- nrow(mat_pt)
out     <- data.frame(matrix(NA_character_, nrow = 2 * n_cells,
                             ncol = ncol(mat_pt)))
for (i in seq_len(n_cells)) {
  out[2 * i - 1, ] <- sprintf("%.3f", mat_pt[i, ])
  out[2 * i,     ] <- sprintf("(%.3f)", mat_se[i, ])
}
out <- cbind(Cell = rep(rownames(mat_pt), each = 2), out)
out$Cell[seq(2, nrow(out), by = 2)] <- ""
colnames(out)[-1] <- colnames(mat_pt)
knitr::kable(out, caption = "Per-cell estimates by variance option (SE in parentheses)")

## ----table-latex, eval = requireNamespace("xtable", quietly = TRUE)-----------
tex <- xtable::xtable(
  tab_A[, c("term", "estimate", "std.error", "conf.low", "conf.high")],
  digits  = c(0, 0, 3, 3, 3, 3),
  caption = "Per-cell RD effects by ideology bucket",
  label   = "tab:rdhte-ideology"
)
print(tex, include.rownames = FALSE, booktabs = TRUE,
      caption.placement = "top")
# writeLines(capture.output(print(tex, include.rownames = FALSE,
#                                booktabs = TRUE,
#                                caption.placement = "top")),
#            con = "rdhte_table_A.tex")

## ----rdrobust-defaults--------------------------------------------------------
summary(rdhte(y = y, x = x))
summary(rdrobust(y = y, x = x))

## ----rdrobust-match-ate-------------------------------------------------------
summary(rdhte(y = y, x = x, h = 0.1, vce = "hc3"))
summary(rdrobust(y = y, x = x, h = 0.1, rho = 1, vce = "hc3"))

## ----rdrobust-match-cells-----------------------------------------------------
summary(rdhte(y = y, x = x, covs.hte = w_left,
              h = c(0.078, 0.116), vce = "hc3"))
summary(rdrobust(y = y[w_left == 1], x = x[w_left == 1],
                 h = 0.116, rho = 1, vce = "hc3"))

## ----cleanup, include = FALSE-------------------------------------------------
detach(rdhte_dataset)