--- 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.