## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(seine)
data(elec_1968)

print(elec_1968)

## -----------------------------------------------------------------------------
elec_1968_turn = ei_proportions(elec_1968, turnout = pres_total,
                                .total = vap, clamp = 0.01)

subset(elec_1968_turn, select = c(fips, state, county, turnout, vap, .other))

## ----echo=FALSE---------------------------------------------------------------
set.seed(35005)
subs = subset(elec_1968, select=c(state, county, pop_city:pop_rural, farm:educ_coll,
                                  inc_00_03k:inc_25_99k, vap_white:vap_other, pres_total))
km = kmeans(model.matrix(~ 0 + . - county - vap_white - vap_black -
                             vap_other - pres_total, subs), centers = 400)
cl = which.max(tapply(subs$vap_black, km$cluster, mad)) # largest Black sd
d_z = subset(subs, km$cluster == cl, c(state, county, pop_urban, farm, educ_elem:inc_25_99k))
d_xn = subset(subs, km$cluster == cl, c(state, county, vap_white:pres_total))

knitr::kable(d_z, digits=4)

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(d_xn, digits=4)

## -----------------------------------------------------------------------------
spec = ei_spec(
    elec_1968,
    predictors = vap_white:vap_other,
    outcome = pres_dem_hum:pres_abs,
    total = pres_total,
    covariates = c(state, pop_city:pop_rural, farm:educ_coll, inc_00_03k:inc_25_99k),
    preproc = function(x) {
        x = model.matrix(~ 0 + ., x) # convert factors to dummies
        bases::b_bart(x, trees = 200)
    }
)

print(spec)

## -----------------------------------------------------------------------------
m = ei_ridge(spec)

print(m)

## -----------------------------------------------------------------------------
m_form = ei_ridge(
    cbind(pres_dem_hum, pres_rep_nix, pres_ind_wal, pres_abs) ~
        vap_white + vap_black + vap_other |
        state * (pop_urban + pop_rural + farm + educ_hsch + educ_coll +
                     inc_03_08k + inc_08_25k + inc_25_99k),
    data = elec_1968, total = pres_total
)

print(m_form)

## -----------------------------------------------------------------------------
summary(m)

## -----------------------------------------------------------------------------
rr = ei_riesz(spec, penalty = m$penalty)

## -----------------------------------------------------------------------------
rr_form = ei_riesz(
    ~ vap_white + vap_black + vap_other |
        state * (pop_urban + pop_rural + farm + educ_hsch + educ_coll +
                     inc_03_08k + inc_08_25k + inc_25_99k),
    data = elec_1968, total = pres_total, penalty = m_form$penalty
)

## -----------------------------------------------------------------------------
summary(rr)

## -----------------------------------------------------------------------------
est = ei_est(m, rr, spec, conf_level = 0.95)
print(est)

## -----------------------------------------------------------------------------
est_form = ei_est(m_form, rr_form, elec_1968)

## -----------------------------------------------------------------------------
est_c = ei_est(m, rr, spec, contrast = list(predictor = c(1, -1, 0)), conf_level = 0.95)
print(est_c)

## -----------------------------------------------------------------------------
as.matrix(est)

as.matrix(est, which = "conf.low")

## -----------------------------------------------------------------------------
as.matrix(ei_est(m, rr, spec, subset = pop_city >= 0.9))
as.matrix(ei_est(m, rr, spec, subset = state == "Mississippi"))

## -----------------------------------------------------------------------------
# Not recommended
est_m = ei_est(regr = m, data = spec)
est_rr = ei_est(riesz = rr, data = spec)

sd(est_m$estimate - est_rr$estimate) # estimates (here) are close
sd(est_m$std.error - est_rr$std.error) # standard errors are very different