---
title: "Panel data: random intercepts, fixed effects, and the Mundlak device"
author: "Benjamin E. Bagozzi"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Panel data: random intercepts, fixed effects, and the Mundlak device}
  %\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)
```

Most applications of the inflated ordered models are panels -- country-years,
agreement-years, campaign-years -- and every estimator in `iop` offers the same
four ways of handling unit heterogeneity:

| Tool                                 | What it does                                                      | Use when                                               |
|--------------------------------------|-------------------------------------------------------------------|--------------------------------------------------------|
| `cluster = "unit"`                   | cluster-robust standard errors, point estimates unchanged         | always sensible as a baseline                          |
| `re = "unit"`                        | a unit random intercept integrated out by adaptive quadrature     | the unit effect is uncorrelated with the covariates    |
| `mundlak()` + `re =` or `cluster =`  | unit means of the covariates (correlated random effects)          | the unit effect may be correlated with the covariates; short panels |
| `fe = "unit"`                        | a dummy per unit (with `fe_correction = "jackknife"`)             | long panels                                            |

The package's Monte Carlo, summarized at the end, is the reason for the "use
when" column.

## Cluster-robust standard errors

`cluster =` names a column and switches `se` to `"cluster"`; the point
estimates are those of the pooled model.

```{r cluster}
library(iop)
data(bp)
f <- violence ~ loggdppc + parliament + disaster | loggdppc + parliament + disaster
m_pool <- iop(f, data = bp, inflate = "bottom")
m_cl   <- iop(f, data = bp, inflate = "bottom", cluster = "country")
round(cbind(estimate = coef(m_pool), se_iid = sqrt(diag(vcov(m_pool))),
            se_cluster = sqrt(diag(vcov(m_cl)))), 3)
```

A nonparametric bootstrap is the other route: `se = "bootstrap"` refits the
model `nboot` times on resamples of the clusters (when `cluster =` is given),
of the random-intercept units, or of the rows, keeps the replicate estimates
in `$boot`, and `confint(type = "percentile")` then gives percentile
intervals. Dale and Sirchenko (2021) find it better calibrated than the
asymptotic standard errors for the error correlation of the correlated models
in small samples. It costs `nboot` refits (`cores =` parallelizes them), so it
is not run here:

```{r boot-se, eval = FALSE}
m_boot <- iop(f, data = bp, inflate = "bottom", se = "bootstrap", cluster = "country",
              nboot = 200, cores = 4)
confint(m_boot, type = "percentile")
```

## Random intercepts

`re = "unit"` adds a normal random intercept to the outcome equation and
integrates it out by adaptive Gauss--Hermite quadrature (`nAGQ` nodes per
unit; the default 15 is accurate for the panels in this package, and fewer
nodes can leave the fit short of a clean optimum, which the fit reports).

```{r re}
m_re <- oprobit(violence ~ loggdppc + parliament + disaster, data = bp, re = "country")
m_re
```

The same fit from `ordinal::clmm(violence ~ loggdppc + parliament + disaster +
(1 | country), link = "probit", nAGQ = 15)` has log-likelihood -1077.816 and
random-intercept SD 1.306, which is what `iop` reports; the package's tests
check this agreement on simulated panels. `ranef()` returns the empirical-Bayes
unit effects:

```{r ranef}
head(ranef(m_re), 4)
```

Reported probabilities and first differences are *marginal* over the random
intercept (population-averaged); `type = "prob_conditional"` gives the
probabilities for a unit with a zero intercept:

```{r re-predict}
head(cbind(marginal = predict(m_re)[, "civil war"],
           conditional = predict(m_re, type = "prob_conditional")[, "civil war"]), 3)
```

Random intercepts work identically for `iop()` and `iol()`, where
`re_inflation = TRUE` also gives the inflation equation an independent unit
intercept (two-dimensional quadrature). A fit whose random-intercept SD goes
to zero is flagged as a boundary case; see `vignette("model")`.

## Unit fixed effects

`fe = "unit"` adds a dummy per unit to the outcome equation (and, with
`fe_inflation = TRUE`, to the inflation equation). Two things happen
automatically:

* units whose response never leaves an extreme category have no finite fixed
  effect and are dropped with a message;
* when the median unit has fewer than 10 observations the fit warns about the
  incidental-parameters bias of maximum-likelihood dummies (Greene 2004).

```{r fe, message = TRUE}
m_fe <- oprobit(violence ~ loggdppc + disaster, data = bp, fe = "country")
m_fe
```

A covariate that does not vary within units is collinear with the dummies and
is refused with a message naming it -- `parliament`, which changes in only one
country of this panel, is the example here:

```{r fe-tiv, eval = FALSE}
oprobit(violence ~ loggdppc + parliament + disaster, data = bp, fe = "country")
#> Error: Covariate(s) parliament do not vary within the units of 'country' on the
#> estimation rows and are collinear with the unit fixed effects. Drop them from the
#> outcome equation, or keep them with re = "country" or mundlak() instead of fe =.
```

### The split-panel jackknife

`fe_correction = "jackknife"` applies the split-panel jackknife of Dhaene and
Jochmans (2015): the model is refit on the first and second half of every
unit's observations (ordered by `time`), and the common parameters are
bias-corrected as `2 * full - (half1 + half2) / 2`. On a simulated short panel
(T = 8, so the fit issues its short-panel warning) with a covariate correlated
with the unit effect, the correction removes most of the incidental-parameters
bias:

```{r jackknife}
set.seed(42)
G <- 100; Tn <- 8
alpha <- rnorm(G); unit <- rep(1:G, each = Tn)
x1 <- rnorm(G * Tn) + 0.5 * alpha[unit]; x2 <- rnorm(G * Tn)
ystar <- 0.8 * x1 - 0.5 * x2 + alpha[unit] + rnorm(G * Tn)
d <- data.frame(y = findInterval(ystar, c(-1, 0.3, 1.2)), x1, x2, unit, t = rep(1:Tn, G))
m_jk <- oprobit(y ~ x1 + x2, data = d, fe = "unit", fe_correction = "jackknife", time = "t")
rbind(uncorrected = m_jk$coefficients_uncorrected[c("x1", "x2")],
      jackknife   = coef(m_jk)[c("x1", "x2")],
      truth       = c(0.8, -0.5))
```

The reported covariance after the correction is the full-sample one: Dhaene
and Jochmans (2015) show that the split-panel jackknife removes the leading
bias without changing the first-order asymptotic variance, so it is the
asymptotically valid covariance for the corrected estimator (`summary()`
says so); in finite samples the corrected estimator is somewhat noisier, so
treat its intervals as approximate. A computational note: the dummies enter
the parameter vector one per unit, and the exact-Hessian steps cost roughly
the square of the parameter count (about 0.5 s at 20 units, 3 s at 100, 18 s
at 300 for an ordered probit on 4,000 rows, more for an inflated model's
multistart), so for panels with many hundreds of units prefer `re =` or
`mundlak()`. The correction is only as good as the two half-panel fits, which
are stored in `m_jk$jackknife`. When a covariate is identified mainly by a within-unit
trend -- as log GDP per capita is in `bp`, where 99 percent of its
within-country variation is linear in time -- the half-panels barely identify
it, their estimates are far from the full-sample one, and the "correction"
inherits that noise. Inspect the half-panel coefficients before reporting a
jackknifed estimate:

```{r jackknife-halves}
m_jk$jackknife$half_coefficients[, c("x1", "x2")]
```

## The Mundlak device

`mundlak()` augments the data with the unit means of the time-varying numeric
covariates of both equations and returns the augmented formula and data. Any
estimator fit on the result implements the correlated-random-effects
specification: the coefficients on the original covariates are the within-unit
effects, and the coefficients on the unit means capture -- and test -- the
correlation between the covariates and the unit effect. It pairs naturally with
cluster-robust standard errors or a random intercept, and it keeps
time-invariant covariates estimable.

```{r mundlak}
md <- mundlak(violence ~ loggdppc + parliament + disaster | loggdppc + parliament + disaster,
              data = bp, unit = "country")
md$formula
md$added
m_md <- oprobit(violence ~ loggdppc + parliament + disaster + loggdppc_mean + disaster_mean,
                data = md$data, cluster = "country")
round(summary(m_md)$coefficients[, 1:2], 3)
```

(The unit mean of `parliament` is not added to the outcome equation above
because it is time-invariant in all but one country, so its mean would be
nearly collinear with the variable itself.)

## The Monte Carlo: which device for which panel

`system.file("mc", package = "iop")` holds a Monte Carlo (`fe_bias.R`) that
compares the four devices on an ordered probit with G = 100 units, a N(0, 1)
unit effect, and a covariate correlated with it (`x1 = N(0,1) + 0.5 * alpha`),
for T = 4, 8, 16, 32 and 100 replications, together with an audit script
(`fe_bias_verify.R`) that checks the pooled bias against its closed form and
each estimator against `ordinal::clm()` / `clmm()` on identical draws. The
percent bias of the `x1` coefficient (true value 0.8):

```{r mc, echo = FALSE}
mc <- read.csv(system.file("mc", "fe_bias_results.csv", package = "iop"))
mc$estimator <- factor(mc$estimator, levels = c("pooled", "fe", "re", "mundlak"),
                       labels = c("pooled", "unit dummies (fe)", "random intercept (re)", "Mundlak"))
tab <- reshape(mc[, c("T", "estimator", "bias_pct")], idvar = "estimator", timevar = "T", direction = "wide")
names(tab) <- c("estimator", paste0("T = ", sort(unique(mc$T))))
tab[-1] <- lapply(tab[-1], function(v) sprintf("%+.1f", v))
knitr::kable(tab, row.names = FALSE, align = c("l", rep("r", 4)),
             caption = "Percent bias of the x1 coefficient, 100 replications per cell")
```

The pooled and random-intercept estimators are biased by the omitted
correlated effect (the pooled bias has a closed form under this design,
+11.8 percent); the dummies carry the incidental-parameters bias that fades
with T; the Mundlak device is essentially unbiased at every T. Hence the
recommendation: `mundlak()` (with `cluster =` or `re =`) for short panels,
`fe =` -- with the jackknife if the half-panels identify the parameters -- for
long ones.

## References

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

Dhaene, G. and Jochmans, K. (2015). Split-panel jackknife estimation of
fixed-effect models. *Review of Economic Studies*, 82, 991-1030.

Greene, W. (2004). The behaviour of the maximum likelihood estimator of
limited dependent variable models in the presence of fixed effects.
*Econometrics Journal*, 7, 98-119.

Mundlak, Y. (1978). On the pooling of time series and cross section data.
*Econometrica*, 46, 69-85.

Wooldridge, J.M. (2010). *Econometric Analysis of Cross Section and Panel
Data*, 2nd ed. MIT Press.