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

## ----cluster------------------------------------------------------------------
library(iop)
data(bp)
f <- violence ~ loggdppc + parliament + disaster | loggdppc + parliament + disaster
m_pool <- iop(f, data = bp, inflate = "bottom")
m_cl   <- iop(f, data = bp, inflate = "bottom", cluster = "country")
round(cbind(estimate = coef(m_pool), se_iid = sqrt(diag(vcov(m_pool))),
            se_cluster = sqrt(diag(vcov(m_cl)))), 3)

## ----boot-se, eval = FALSE----------------------------------------------------
# m_boot <- iop(f, data = bp, inflate = "bottom", se = "bootstrap", cluster = "country",
#               nboot = 200, cores = 4)
# confint(m_boot, type = "percentile")

## ----re-----------------------------------------------------------------------
m_re <- oprobit(violence ~ loggdppc + parliament + disaster, data = bp, re = "country")
m_re

## ----ranef--------------------------------------------------------------------
head(ranef(m_re), 4)

## ----re-predict---------------------------------------------------------------
head(cbind(marginal = predict(m_re)[, "civil war"],
           conditional = predict(m_re, type = "prob_conditional")[, "civil war"]), 3)

## ----fe, message = TRUE-------------------------------------------------------
m_fe <- oprobit(violence ~ loggdppc + disaster, data = bp, fe = "country")
m_fe

## ----fe-tiv, eval = FALSE-----------------------------------------------------
# oprobit(violence ~ loggdppc + parliament + disaster, data = bp, fe = "country")
# #> Error: Covariate(s) parliament do not vary within the units of 'country' on the
# #> estimation rows and are collinear with the unit fixed effects. Drop them from the
# #> outcome equation, or keep them with re = "country" or mundlak() instead of fe =.

## ----jackknife----------------------------------------------------------------
set.seed(42)
G <- 100; Tn <- 8
alpha <- rnorm(G); unit <- rep(1:G, each = Tn)
x1 <- rnorm(G * Tn) + 0.5 * alpha[unit]; x2 <- rnorm(G * Tn)
ystar <- 0.8 * x1 - 0.5 * x2 + alpha[unit] + rnorm(G * Tn)
d <- data.frame(y = findInterval(ystar, c(-1, 0.3, 1.2)), x1, x2, unit, t = rep(1:Tn, G))
m_jk <- oprobit(y ~ x1 + x2, data = d, fe = "unit", fe_correction = "jackknife", time = "t")
rbind(uncorrected = m_jk$coefficients_uncorrected[c("x1", "x2")],
      jackknife   = coef(m_jk)[c("x1", "x2")],
      truth       = c(0.8, -0.5))

## ----jackknife-halves---------------------------------------------------------
m_jk$jackknife$half_coefficients[, c("x1", "x2")]

## ----mundlak------------------------------------------------------------------
md <- mundlak(violence ~ loggdppc + parliament + disaster | loggdppc + parliament + disaster,
              data = bp, unit = "country")
md$formula
md$added
m_md <- oprobit(violence ~ loggdppc + parliament + disaster + loggdppc_mean + disaster_mean,
                data = md$data, cluster = "country")
round(summary(m_md)$coefficients[, 1:2], 3)

## ----mc, echo = FALSE---------------------------------------------------------
mc <- read.csv(system.file("mc", "fe_bias_results.csv", package = "iop"))
mc$estimator <- factor(mc$estimator, levels = c("pooled", "fe", "re", "mundlak"),
                       labels = c("pooled", "unit dummies (fe)", "random intercept (re)", "Mundlak"))
tab <- reshape(mc[, c("T", "estimator", "bias_pct")], idvar = "estimator", timevar = "T", direction = "wide")
names(tab) <- c("estimator", paste0("T = ", sort(unique(mc$T))))
tab[-1] <- lapply(tab[-1], function(v) sprintf("%+.1f", v))
knitr::kable(tab, row.names = FALSE, align = c("l", rep("r", 4)),
             caption = "Percent bias of the x1 coefficient, 100 replications per cell")