---
title: "Predicted probabilities, first differences, and marginal effects"
author: "Benjamin E. Bagozzi"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Predicted probabilities, first differences, and marginal effects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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)
```

An inflated ordered model has more moving parts than a plain ordered model:
category probabilities, regime probabilities, and the posterior probability
that an observation in the inflated category is an inflated case, each of
which can change when a covariate moves in one equation, the other, or both.
This vignette walks through the tools that turn a fit into those quantities
and into tables and figures. The running example is the zero-inflated ordered
probit of political violence from `vignette("iop")`:

```{r 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")
```

Everything below works the same way for `oprobit()`, `ologit()`, and `iol()`
fits; the regime-related quantities are simply absent (or equal to one) for the
uninflated models.

## Predicted probabilities

`predict()` has one `type` per quantity:

| `type`               | Returns                                                                 |
|----------------------|-------------------------------------------------------------------------|
| `"prob"` (default)   | the n x J matrix of category probabilities P(y = j)                     |
| `"class"`            | the modal category                                                      |
| `"prob_outcome"`     | the ordered-stage probabilities P(y = j \| ordered regime)              |
| `"regime"`           | P(ordered regime)                                                       |
| `"inflated"`         | P(inflated regime) = 1 - P(ordered regime)                              |
| `"posterior"`        | P(inflated regime \| y = inflated category), 0 for other categories     |
| `"prob_conditional"` | for random-intercept fits, P(y = j) at a zero random intercept          |
| `"link"`, `"link_inflation"` | the two linear predictors                                       |

```{r predict-types}
head(predict(m), 3)
head(predict(m, type = "class"), 3)
head(predict(m, type = "prob_outcome"), 3)
summary(predict(m, type = "regime"))
```

The posterior probability separates, among the observed zeros, those the model
attributes to the inflated regime (structurally peaceful country-years) from
those it attributes to the ordered regime (at-risk country-years that stayed
calm):

```{r posterior}
post <- predict(m, type = "posterior")
summary(post[bp$violence == "none"])
```

### The two types of zeros

Harris and Zhao (2007) decompose the probability of the inflated category into
its two sources: `type = "zeros"` returns, for every observation, the mass
that arrives through the inflation process (\(1 - P(s = 1)\), the
structurally peaceful) and the mass that arrives through the ordered stage
(\(P(s = 1, y^* = k)\), at-risk country-years that stayed calm); the two
columns sum to \(P(y = k)\):

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

`first_difference()` and `ame()` report effects on both components with
`decompose = TRUE`, which is how Harris and Zhao present the marginal effect
on the zero probability:

```{r zeros-fd}
first_difference(m, "loggdppc", from = 7, to = 9, decompose = TRUE)
```

Two further summaries are `type = "mean"` (the expected category index
\(\sum_j j\,P(y = j)\) on the 0, ..., J-1 scale) and `type = "cumulative"`
(\(P(y \le j)\)).

### Covariate profiles and standard errors

`newdata` takes any data frame with the model's covariates, and
`se.fit = TRUE` adds delta-method standard errors for every probability type:

```{r 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)
```

## First differences

`first_difference()` moves one covariate from one value to another, holding
the other covariates at a profile -- by default the weighted means of numeric
covariates and the modal levels of factors -- and reports the change in every
category probability and, for inflated models, in the probability of the
ordered regime, with delta-method intervals:

```{r fd}
first_difference(m, "loggdppc", from = 7, to = 9)
```

For inflated models the covariate can be moved in the outcome equation only,
in the inflation equation only, or in both (the total effect, the default):

```{r fd-stage}
first_difference(m, "loggdppc", from = 7, to = 9, stage = "outcome")
first_difference(m, "loggdppc", from = 7, to = 9, stage = "inflation")
```

Averaging over the estimation data instead of fixing a profile
(`average = TRUE`) lets every observation keep its own other covariates (the
average honours the fit's `weights`, as does `ame()`), and
`ci = "sim"` replaces the delta method by simulation from the asymptotic
distribution of the parameters (Krinsky and Robb 1986; King, Tomz and
Wittenberg 2000):

```{r 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)
```

A one-row `newdata` sets the profile explicitly, and the result has a `plot()`
method:

```{r 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")
```

## Average marginal effects

`ame()` reports, for every covariate, the average (over the estimation data)
effect on every category probability and on the regime probability: a
derivative for continuous covariates, the discrete change 0 to 1 for binary
ones, and each level against the base level for factors. A covariate is moved
in every equation in which it appears.

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

Like `first_difference()`, `ame()` can move a covariate in one equation only
(`stage = "outcome"` or `"inflation"`); the two stage-specific effects add up
to the total effect to first order, and a covariate absent from the selected
equation has a zero effect there:

```{r ame-stage}
ame(m, vars = "loggdppc", stage = "inflation")
```

## The response distribution: `diord`, `piord`, `qiord`, `riord`

The probability mass, cumulative probability, quantile, and random-draw
functions of the (inflated) ordered response are available in the usual
`d/p/q/r` form, either from explicit parameters -- the outcome linear
predictor `eta`, the cutpoints `tau`, and for inflated models the split
predictor `a`, the inflated category `k`, and `rho` -- or from a fitted
model at the estimation data or at `newdata`:

```{r 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
```

`riop()` remains the full data-generating simulator (it also draws the
covariates); `simulate()` draws replicate responses for a fit.

## Regression tables

`tidy()`, `glance()`, and `augment()` methods are registered with
[broom](https://broom.tidymodels.org) when it is installed, which is also what
`modelsummary` uses; `texreg` is served by an `extract()` method. Inflation
coefficients keep an `infl_` prefix in `tidy()` and are labelled
"Inflation: ..." in `texreg`, so the two equations stay distinct.

```{r broom, eval = has("broom")}
broom::tidy(m, conf.int = TRUE)[1:4, ]
broom::glance(m)
```

```{r 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)
```

```{r 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 and accuracy scores

`classification()` summarizes how well the predicted probabilities reproduce
the observed categories: the classification table of observed against modal
predicted categories, the share correctly classified, the Brier score and the
ranked probability score (proper scoring rules, smaller is better), the mean
log score, and the per-category precision, recall, and adjusted
noise-to-signal ratio -- the suite Dale and Sirchenko (2021) use to compare
ordered and inflated fits. It works on the estimation data or on new data:

```{r 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)
```

## Diagnostics

Ordinal outcomes have no canonical residual. `residuals()` offers the
observed-minus-expected category index (`"response"`) and its Pearson version,
but the recommended check is simulation: `simulate()` draws replicate response
vectors from the fitted probabilities, which feed
[DHARMa](https://CRAN.R-project.org/package=DHARMa)'s simulated-residual
diagnostics:

```{r 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)
```

For random-intercept fits, `simulate()` draws fresh unit intercepts in every
replicate, so the diagnostics are marginal over the random effects; see
`vignette("panels")`.

## References

Dale, D. and Sirchenko, A. (2021). Estimation of nested and zero-inflated
ordered probit models. *Stata Journal*, 21, 3-38.

Harris, M.N. and Zhao, X. (2007). A zero-inflated ordered probit model, with an
application to modelling tobacco consumption. *Journal of Econometrics*, 141,
1073-1099.

King, G., Tomz, M. and Wittenberg, J. (2000). Making the most of statistical
analyses: improving interpretation and presentation. *American Journal of
Political Science*, 44, 347-361.

Krinsky, I. and Robb, A.L. (1986). On approximating the statistical properties
of elasticities. *Review of Economics and Statistics*, 68, 715-719.
