## ----setup, include=FALSE-----------------------------------------------------
options(marginaleffects_print_ncols = 10)
options(marginaleffects_print_type = FALSE)
options(width = 80)
knitr::opts_chunk$set(
  fig.width = 5,
  fig.height = 5 * 0.618,
  fig.align = "center",
  cache = FALSE
)


## ----include = FALSE----------------------------------------------------------
library("brms")
library("mgcv")
library("ggdist")
library("ggplot2")
library("marginaleffects")
library("modelsummary")
dat <- read.csv("https://marginaleffects.com/data/impartiality.csv")
head(dat)
# load locally from the supplemental information package
# dat <- read.csv("impartiality.csv")

## ----echo=FALSE, warnings=FALSE-----------------------------------------------
m <- glm(
  impartial ~ equal * democracy + continent, 
  data = dat, family = binomial
)

modelsummary(m, title = "Logistic regression model with impartiality as 
             the outcome.")

## ----calc-background-stuff, include=FALSE, cache = FALSE----------------------
library("janitor")
library("dplyr")
library("knitr")
library("broom")

details <- list(
  n_countries = nrow(dat)
)

fmt_pvalue_p <- scales::label_pvalue(prefix = c("p < ", "p = ", "p > "))
fmt_pvalue <- scales::label_pvalue(prefix = c("< ", "", "> "))


## ----calc-nice-predictions, echo=FALSE, cache=FALSE, warning=FALSE------------
avg_values <- datagrid(model = m) %>%
  mutate(equal_nice = round(equal, 1))

all_predictions <- predictions(m) %>%
  mutate(estimate_nice = round(estimate, 3)) %>%
  mutate(country_name = dat$country_name)

avg_preds_mean <- predictions(m, newdata = "mean") %>%
  mutate(estimate_nice = round(estimate, 3))

avg_preds <- avg_predictions(m) %>%
  mutate(estimate_nice = round(estimate, 3))

preds_regime_eqdr_estimands <- predictions(m,
  newdata = datagrid(democracy = unique, equal = c(30, 90))
)

preds_regime_eqdr <- preds_regime_eqdr_estimands %>%
  mutate(estimate_nice = round(estimate, 3)) %>%
  split(~ democracy + equal)

preds_by_regime <- predictions(m, by = "democracy") %>%
  mutate(estimate_nice = round(estimate, 3)) %>%
  split(~democracy)

preds_hypo <- avg_predictions(m,
  by = c("democracy"), vcov = ~continent, hypothesis = "revpairwise"
) %>%
  mutate(
    estimate_pp = round(estimate * 100, 1),
    p_nice = fmt_pvalue(p.value),
    p_nice_eq = fmt_pvalue_p(p.value)
  )

preds_hypo_regime <- bind_rows(
  autocracy = hypotheses(preds_regime_eqdr_estimands, "b3 = b4"),
  democracy = hypotheses(preds_regime_eqdr_estimands, "b1 = b2"),
  .id = "democracy"
) %>%
  mutate(
    estimate_pp = round(estimate * 100, 1),
    p_nice = fmt_pvalue(p.value),
    p_nice_eq = fmt_pvalue_p(p.value)
  ) %>%
  split(~democracy)


## ----pred-step2-predictions-all, warning=FALSE--------------------------------
p <- predictions(m)
p


## -----------------------------------------------------------------------------
library("dplyr")
p$estimate[1:4]
p[2, "p.value"]
p |> filter(estimate == min(estimate))


## -----------------------------------------------------------------------------
options(marginaleffects_print_column_names = FALSE)


## ----warning=FALSE------------------------------------------------------------
predictions(m, newdata = head(dat, 2), type = "link")


## ----pred-step2-predictions-mean, warning=FALSE-------------------------------
predictions(m, newdata = "mean") |> data.frame()


## ----pred-step2-grid-example, warning=FALSE-----------------------------------
datagrid(model = m, democracy = unique, equal = c(30, 90))


## ----pred-step2-grid-predictions, warning=FALSE-------------------------------
predictions(m,
  newdata = datagrid(democracy = unique, equal = c(30, 90))
)


## ----include=FALSE------------------------------------------------------------
p <- avg_predictions(m, type = "response")


## ----pred-step3-avg-predictions-----------------------------------------------
avg_predictions(m, type = "response")


## -----------------------------------------------------------------------------
mean(predict(m, newdata = dat, type = "response"))


## ----pred-step3-avg-by-democracy----------------------------------------------
predictions(m, by = "democracy")


## ----plotpredictions1, fig.cap="Average predicted probability of impartiality by regime type and continent."----
plot_predictions(m, by = c("democracy", "continent"))


## ----plotpredictions2, fig.height = 4, fig.cap="Predicted probability of impartiality by levels of equality and democracy.", out.width="70%"----
library("ggplot2")
theme_set(theme_bw(base_family = "serif"))
okabeito <- c("#E69F00", "#56B4E9")

plot_predictions(m, condition = c("equal", "democracy")) +
  labs(color = NULL, fill = NULL,
       x = "Equality", y = "Predicted probability of impartiality") +
  scale_colour_manual(values = okabeito) +
  scale_fill_manual(values = okabeito) +
  theme(legend.position = "bottom")


## ----pred-step4-robust-ses----------------------------------------------------
avg_predictions(m,
  by = "democracy",
  vcov = "HC3",
  conf_level = .99
)


## ----pred-step4-cluster-robust-ses--------------------------------------------
avg_predictions(m,
  by = "continent",
  vcov = ~continent,
  conf_level = .95
)


## ----pred-step4-simulation-ses, messages = FALSE, warnings = FALSE------------
set.seed(1024)
avg_predictions(m, by = "democracy") |>
  inferences(method = "simulation")


## -----------------------------------------------------------------------------
coef(m)[c("continentAsia", "continentAmericas")]

hypotheses(m, hypothesis = "continentAsia = continentAmericas")


## -----------------------------------------------------------------------------
avg_predictions(m,
  by = "democracy",
  type = "response")


## -----------------------------------------------------------------------------
avg_predictions(m,
  by = "democracy",
  type = "response",
  hypothesis = "revpairwise")


## -----------------------------------------------------------------------------
predictions(m,
  by = "democracy",
  type = "response",
  hypothesis = "b2 = b1 * 2")


## ----include=FALSE------------------------------------------------------------
tmp <- avg_predictions(m, by = "democracy", type = "response") 
res <- sprintf("%.4f", tmp$estimate[2] - 2 * tmp$estimate[1])
tmp <- transform(tmp, estimate = sprintf("%.4f", estimate))


## -----------------------------------------------------------------------------
predictions(m,
  by = "democracy",
  type = "response",
  hypothesis = "b2 = b1 * 2",
  equivalence = c(-.2, .2))


## ----calc-nice-comparisons, include=FALSE-------------------------------------
cmp_average <- avg_comparisons(m,
  variables = list(democracy = c("Democracy", "Autocracy"))
) %>%
  mutate(contrast = janitor::make_clean_names(contrast)) %>%
  mutate(estimate_pp = round(estimate * 100, 1))

cmp_regime <- comparisons(m,
  variables = list(equal = c(30, 90)),
  newdata = datagrid(democracy = unique)
) %>%
  mutate(estimate_pp = round(estimate * 100, 1)) %>%
  split(~democracy)

cmp_regime_hypo <- comparisons(m,
  variables = list(equal = c(30, 90)),
  newdata = datagrid(democracy = unique),
  hypothesis = "b2 = b1"
) %>%
  mutate(
    estimate_pp = round(estimate * 100, 1),
    p_nice = fmt_pvalue(p.value),
    p_nice_eq = fmt_pvalue_p(p.value)
  )


## ----comparisons-regime-original-data-----------------------------------------
comparisons(m, variables = "democracy")


## ----comparisons-regime-avg---------------------------------------------------
avg_comparisons(m)


## -----------------------------------------------------------------------------
dat_lo <- transform(dat, democracy = "Autocracy")
dat_hi <- transform(dat, democracy = "Democracy")
pred_lo <- predict(m, newdata = dat_lo, type = "response")
pred_hi <- predict(m, newdata = dat_hi, type = "response")
mean(pred_hi - pred_lo)


## ----eval = FALSE-------------------------------------------------------------
## avg_comparisons(m, variables = list("equal" = 4))
## avg_comparisons(m, variables = list("equal" = "sd"))
## avg_comparisons(m, variables = list("equal" = "iqr"))
## avg_comparisons(m, variables = list("equal" = c(30, 90)))


## -----------------------------------------------------------------------------
avg_comparisons(m, variables = "democracy", comparison = "ratio")


## -----------------------------------------------------------------------------
avg_comparisons(m,
  comparison = "lnor",
  transform = exp)


## -----------------------------------------------------------------------------
avg_comparisons(m,
  variables = "equal",
  comparison = \(hi, lo) mean(hi) / mean(lo))


## ----include=FALSE------------------------------------------------------------
cmp <- comparisons(m,
  by = "democracy",
  variables = list(equal = c(30, 90)))


## -----------------------------------------------------------------------------
cmp <- comparisons(m,
  by = "democracy",
  variables = list(equal = c(30, 90)))
cmp


## ----fig.cap = "Effect of a change from 30 to 90 in resource equality on the predicted probability of having impartial public institutions.\\label{fig:comparisons-democracy}"----
plot_comparisons(m,
  by = "democracy",
  variables = list(equal = c(30, 90))) +
  labs(x = NULL, y = "Risk Difference (Impartiality)")


## -----------------------------------------------------------------------------
cmp <- comparisons(m,
  by = "democracy",
  variables = list(equal = c(30, 90)),
  hypothesis = "pairwise")
cmp


## ----tangents, echo=FALSE, fig.cap = "Tangents to the prediction function at 25 and 50.", fig.pos="h"----
p = predictions(m, datagrid(equal = c(25, 50)))
s = slopes(m, datagrid(equal = c(25, 50)), variables = "equal")
tan1 <- data.frame(x = c(25, 50), y = c(0.726, 0.726 + 0.01667 * (50 - 25)))
tan2 <- data.frame(x = c(25, 75), y = c(0.956, 0.956 + 0.00355 * (75 - 50)))
plot_predictions(m, condition = "equal", vcov = FALSE) +
  geom_abline(slope = s$estimate[1], intercept = p$estimate[1] - 25 * s$estimate[1], color = "red", linetype = 2) +
  geom_abline(slope = s$estimate[2], intercept = p$estimate[2] - 50 * s$estimate[2], color = "red", linetype = 2) +
  theme_bw()


## -----------------------------------------------------------------------------
slopes(m, newdata = datagrid(equal = c(25, 50)), variables = "equal")


## -----------------------------------------------------------------------------
avg_slopes(m, variables = "equal")


## -----------------------------------------------------------------------------
slopes(m, variables = "equal", newdata = "mean")


## -----------------------------------------------------------------------------
slopes(m, variables = "equal", newdata = "median", by = "democracy")


## -----------------------------------------------------------------------------
avg_slopes(m, variables = "equal", slope = "eyex")

sessionInfo()