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

## ----fit----------------------------------------------------------------------
library(iop)
data(bp)
m <- iop(violence ~ loggdppc + parliament + disaster + major_oil + major_primary |
           loggdppc + parliament + disaster + major_oil + major_primary,
         data = bp, inflate = "bottom")

## ----predict-types------------------------------------------------------------
head(predict(m), 3)
head(predict(m, type = "class"), 3)
head(predict(m, type = "prob_outcome"), 3)
summary(predict(m, type = "regime"))

## ----posterior----------------------------------------------------------------
post <- predict(m, type = "posterior")
summary(post[bp$violence == "none"])

## ----zeros--------------------------------------------------------------------
z <- predict(m, type = "zeros")
head(cbind(z, total = rowSums(z), P_none = predict(m)[, "none"]), 3)
colMeans(z)

## ----zeros-fd-----------------------------------------------------------------
first_difference(m, "loggdppc", from = 7, to = 9, decompose = TRUE)

## ----newdata------------------------------------------------------------------
nd <- data.frame(loggdppc = c(6, 8, 10), parliament = 0, disaster = 0,
                 major_oil = 0, major_primary = 0)
p <- predict(m, newdata = nd, se.fit = TRUE)
round(p$fit, 3)
round(p$se.fit, 3)
predict(m, newdata = nd, type = "inflated", se.fit = TRUE)

## ----fd-----------------------------------------------------------------------
first_difference(m, "loggdppc", from = 7, to = 9)

## ----fd-stage-----------------------------------------------------------------
first_difference(m, "loggdppc", from = 7, to = 9, stage = "outcome")
first_difference(m, "loggdppc", from = 7, to = 9, stage = "inflation")

## ----fd-avg-------------------------------------------------------------------
first_difference(m, "loggdppc", from = 7, to = 9, average = TRUE)
first_difference(m, "major_oil", from = 0, to = 1, ci = "sim", R = 500)

## ----fd-plot------------------------------------------------------------------
fd <- first_difference(m, "disaster", from = 0, to = 3,
                       newdata = data.frame(loggdppc = 8, parliament = 1, disaster = 0,
                                            major_oil = 0, major_primary = 0))
plot(fd, main = "Three disasters vs none, parliamentary democracy at log GDP 8")

## ----ame, fig.height = 5------------------------------------------------------
a <- ame(m, vars = c("loggdppc", "disaster", "major_oil"))
a
plot(a)

## ----ame-stage----------------------------------------------------------------
ame(m, vars = "loggdppc", stage = "inflation")

## ----dpqr---------------------------------------------------------------------
diord(0:2, eta = 0.3, tau = c(-0.5, 0.8), a = 0.4, k = 0)        # P(y = j) at one profile
piord("repression", object = m, newdata = bp[1:3, ])              # P(y <= repression)
qiord(0.5, object = m, newdata = bp[1:3, ])                       # median category
table(riord(nrow(bp), object = m))                                # one draw per observation

## ----broom, eval = has("broom")-----------------------------------------------
broom::tidy(m, conf.int = TRUE)[1:4, ]
broom::glance(m)

## ----texreg, eval = has("texreg")---------------------------------------------
m_op <- oprobit(violence ~ loggdppc + parliament + disaster + major_oil + major_primary, data = bp)
texreg::screenreg(list(m_op, m), custom.model.names = c("Ordered probit", "ZiOP"),
                  include.cutpoints = FALSE, digits = 3)

## ----modelsummary, eval = has("modelsummary") && has("broom"), results = "asis"----
modelsummary::modelsummary(list("Ordered probit" = m_op, "ZiOP" = m), output = "markdown",
                           stars = TRUE, gof_map = c("nobs", "logLik", "AIC", "BIC"))

## ----classification-----------------------------------------------------------
m_op <- oprobit(violence ~ loggdppc + parliament + disaster + major_oil + major_primary, data = bp)
classification(m)
c(`ordered probit` = classification(m_op)$brier, ZiOP = classification(m)$brier)

## ----dharma, eval = has("DHARMa")---------------------------------------------
sims <- simulate(m, nsim = 250)
res <- DHARMa::createDHARMa(simulatedResponse = as.matrix(sims), observedResponse = m$y,
                            fittedPredictedResponse = as.numeric(fitted(m) %*% (0:2)),
                            integerResponse = TRUE)
plot(res)