Panel data: random intercepts, fixed effects, and the Mundlak device

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:

The package’s Monte Carlo, summarized at the end, is the reason for the “use when” column.

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

Cluster-robust standard errors

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

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)
#>                      estimate se_iid se_cluster
#> loggdppc                0.005  0.056      0.170
#> parliament             -0.085  0.172      0.392
#> disaster                0.269  0.033      0.059
#> none|repression         0.519  0.407      1.231
#> repression|civil war    1.411  0.416      1.259
#> infl_(Intercept)       19.411  3.515      4.819
#> infl_loggdppc          -2.123  0.376      0.492
#> infl_parliament        -0.427  0.374      0.751
#> infl_disaster          -0.120  0.141      0.249

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:

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

m_re <- oprobit(violence ~ loggdppc + parliament + disaster, data = bp, re = "country")
m_re
#> Ordered probit, random intercept by country 
#> Response levels (in order): none < repression < civil war 
#> Call:  oprobit(formula = violence ~ loggdppc + parliament + disaster, 
#>     data = bp, re = "country")
#> 
#> Outcome coefficients:
#>   loggdppc parliament   disaster 
#>    -0.4874    -1.0941     0.1132 
#> 
#> Cutpoints:
#>      none|repression repression|civil war 
#>              -3.0180              -1.6385 
#> 
#> Random-intercept SD: 1.3052  over 113 units
#> 
#> logLik: -1077.82   N: 1984

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:

head(ranef(m_re), 4)
#>          unit          u      u_sd
#> 1 Afghanistan  0.4995294 0.4238394
#> 2    Albania  -1.2891857 0.6875121
#> 3    Algeria   0.4122190 0.2640184
#> 4      Angola  1.0791696 0.3109542

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:

head(cbind(marginal = predict(m_re)[, "civil war"],
           conditional = predict(m_re, type = "prob_conditional")[, "civil war"]), 3)
#>        marginal conditional
#> [1,] 0.03423599 0.001369266
#> [2,] 0.05203033 0.003762048
#> [3,] 0.03249352 0.001205750

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

m_fe <- oprobit(violence ~ loggdppc + disaster, data = bp, fe = "country")
#> 37 of 113 units have no within-unit variation that identifies their fixed effect (response constant at an extreme category); their 579 observations are dropped.
m_fe
#> Ordered probit 
#> Response levels (in order): none < repression < civil war 
#> Call:  oprobit(formula = violence ~ loggdppc + disaster, data = bp, 
#>     fe = "country")
#> 
#> Outcome coefficients:
#> loggdppc disaster 
#>  -0.4947   0.0920 
#> (75 unit fixed-effect dummies by country not shown; see coef())
#> 
#> Cutpoints:
#>      none|repression repression|civil war 
#>              -3.6579              -2.2147 
#> 
#> logLik: -905.57   N: 1405

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:

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:

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")
#> Warning in .iord_fit(formula, data, link, inflate, correlated, parallel, : Unit
#> fixed effects with few observations per unit (median 8): maximum-likelihood
#> dummies carry incidental-parameters bias in nonlinear models (Greene 2004);
#> consider re = "unit" or mundlak() with cluster = "unit".
rbind(uncorrected = m_jk$coefficients_uncorrected[c("x1", "x2")],
      jackknife   = coef(m_jk)[c("x1", "x2")],
      truth       = c(0.8, -0.5))
#>                    x1         x2
#> uncorrected 1.0148638 -0.6726296
#> jackknife   0.8556086 -0.5460368
#> truth       0.8000000 -0.5000000

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:

m_jk$jackknife$half_coefficients[, c("x1", "x2")]
#>             x1         x2
#> half1 1.324288 -0.8560121
#> half2 1.023950 -0.7424329

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.

md <- mundlak(violence ~ loggdppc + parliament + disaster | loggdppc + parliament + disaster,
              data = bp, unit = "country")
md$formula
#> violence ~ loggdppc + parliament + disaster + loggdppc_mean + 
#>     parliament_mean + disaster_mean | loggdppc + parliament + 
#>     disaster + loggdppc_mean + parliament_mean + disaster_mean
md$added
#> [1] "loggdppc_mean"   "parliament_mean" "disaster_mean"
m_md <- oprobit(violence ~ loggdppc + parliament + disaster + loggdppc_mean + disaster_mean,
                data = md$data, cluster = "country")
round(summary(m_md)$coefficients[, 1:2], 3)
#>                      Estimate Std. Error
#> loggdppc               -0.263      0.326
#> parliament             -0.610      0.345
#> disaster                0.052      0.028
#> loggdppc_mean           0.057      0.338
#> disaster_mean           0.323      0.115
#> none|repression        -0.954      0.776
#> repression|civil war   -0.099      0.794

(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):

Percent bias of the x1 coefficient, 100 replications per cell
estimator	T = 4	T = 8	T = 16	T = 32
pooled	+11.9	+11.6	+11.7	+11.8
unit dummies (fe)	+26.4	+10.7	+5.3	+2.5
random intercept (re)	+20.6	+10.8	+5.8	+2.9
Mundlak	+0.3	-0.1	+0.2	-0.0

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.