---
title: "Correlated random effects in linear and nonlinear models"
author: "Xiang Ao"
date: "2025-11-11"
---
```{r}
#| label: setup
#| include: false
library(dplyr)
library(fixest)
library(bacondecomp)
library(lme4)
library(marginaleffects)
```
## The CRE idea
Suppose a panel-data model has an unobserved individual effect $v_i$:
$$y_{it} = \alpha + X_{it}\beta + v_i + \epsilon_{it}$$
The usual random-effects model assumes that $v_i$ is independent of the covariates, so $E[v_i \mid X_i] = 0$. This is often hard to defend. Fixed effects allows $v_i$ to be correlated with $X_{it}$ by removing $v_i$ through demeaning or first-differencing, but then we cannot estimate coefficients on time-invariant covariates. The fixed-effects approach is also less convenient in nonlinear models.
Correlated random effects (CRE) provides another way to allow this correlation (Mundlak, 1978; Chamberlain, 1982; Wooldridge, 2010). We write the projection of $v_i$ on the history of $X_i$ as
$$v_i = \psi + \bar{X}_i \xi + a_i, \qquad E[a_i \mid X_i] = 0$$
where $\bar{X}_i = \frac{1}{T}\sum_{t=1}^T X_{it}$ is the unit mean of the time-varying covariates. Substituting it into the outcome model gives
$$y_{it} = (\alpha + \psi) + X_{it}\beta + \bar{X}_i \xi + a_i + \epsilon_{it}$$
So the practical step is simple: add the unit means $\bar{X}_i$ to the regression.
## Linear models
In a balanced linear panel, pooled OLS or random-effects GLS with $\bar{X}_i$ gives the same coefficient on $X_{it}$ as the within estimator.
### An example using the `castle` data
The `castle` data form a balanced panel of 50 states observed for 11 years ($N = 550$). I use log homicide as the outcome and poverty and log police employment as the covariates.
```{r}
#| label: fix1
#| message: false
#| warning: false
data("castle")
# 1. Two-way Fixed Effects
fe_model <- feols(l_homicide ~ poverty + l_police | state + year, data = castle)
# 2. Standard Random Effects (assumes Cov(X, v) = 0)
re_model <- lmer(l_homicide ~ poverty + l_police + factor(year) + (1 | state), data = castle)
# 3. CRE Pooled OLS (Mundlak device)
castle2 <- castle |>
group_by(state) |>
mutate(poverty_mean = mean(poverty, na.rm = TRUE),
l_police_mean = mean(l_police, na.rm = TRUE))
cre_pols <- feols(l_homicide ~ poverty + poverty_mean + l_police + l_police_mean + factor(year),
data = castle2, cluster = ~state)
# 4. CRE Random Effects
cre_re <- lmer(l_homicide ~ poverty + poverty_mean + l_police + l_police_mean + factor(year) + (1 | state),
data = castle2)
# Compare estimates on poverty and l_police
res_table <- data.frame(
Estimator = c("Two-way Fixed Effects", "Standard Random Effects", "CRE Pooled OLS", "CRE Random Effects"),
beta_poverty = c(coef(fe_model)["poverty"], fixef(re_model)["poverty"],
coef(cre_pols)["poverty"], fixef(cre_re)["poverty"]),
beta_police = c(coef(fe_model)["l_police"], fixef(re_model)["l_police"],
coef(cre_pols)["l_police"], fixef(cre_re)["l_police"])
)
res_table[, 2:3] <- round(res_table[, 2:3], 6)
knitr::kable(res_table, caption = "Comparison of linear panel estimators.")
```
### Comparing the estimates
The two-way fixed-effects model, CRE pooled OLS, and CRE random-effects model give the same coefficients: $\beta_{\text{poverty}} = -0.027068$ and $\beta_{\text{police}} = 0.066033$. Standard random effects, which omits the unit means, gives $\beta_{\text{poverty}} = -0.000799$. The difference is large because the standard model rules out correlation with the state effect. For pooled OLS, we should cluster the standard errors by state because observations within a state are correlated.
## Binary outcomes
For a binary outcome $y_{it} \in \{0, 1\}$, the latent model is:
$$P(y_{it} = 1 \mid X_{it}, v_i) = \Phi(X_{it}\beta + v_i)$$
The exact equivalence between fixed effects and CRE no longer holds in nonlinear models. Conditional logit requires conditional independence of $y_{it}$ over time given $X_i$ and $v_i$; serial correlation in the idiosyncratic errors violates this assumption. It also conditions $v_i$ out of the likelihood, so we cannot recover partial effects that depend on $v_i$.
An alternative is pooled probit with the unit means included. With standard errors clustered by unit, this approach allows arbitrary serial correlation (Lin and Wooldridge, 2019), and average marginal effects are easy to calculate.
### An example with a binary outcome
For an illustration, I define `high_homicide` to equal one when the homicide rate is above the sample mean:
```{r}
#| label: fix2
#| message: false
#| warning: false
castle2 <- castle2 |>
mutate(high_homicide = ifelse(l_homicide > mean(l_homicide, na.rm = TRUE), 1, 0))
# 1. CRE Pooled Probit
cre_probit <- glm(high_homicide ~ poverty + poverty_mean + l_police + l_police_mean + factor(year),
data = castle2, family = binomial(link = "probit"))
# 2. Linear Probability Model with CRE
cre_lpm <- lm(high_homicide ~ poverty + poverty_mean + l_police + l_police_mean + factor(year),
data = castle2)
# 3. CRE Random Effects Probit
cre_re_probit <- glmer(high_homicide ~ poverty + poverty_mean + l_police + l_police_mean + factor(year) + (1 | state),
data = castle2, family = binomial(link = "probit"))
# Compute Average Marginal Effects (AMEs)
ame_probit <- avg_slopes(cre_probit, variables = c("poverty", "l_police"))
ame_lpm <- avg_slopes(cre_lpm, variables = c("poverty", "l_police"))
ame_re <- avg_slopes(cre_re_probit, variables = c("poverty", "l_police"))
ame_comp <- data.frame(
Model = c("CRE Pooled Probit", "CRE Linear Probability Model", "CRE Random Effects Probit"),
AME_poverty = c(ame_probit$estimate[ame_probit$term == "poverty"],
ame_lpm$estimate[ame_lpm$term == "poverty"],
ame_re$estimate[ame_re$term == "poverty"]),
SE_poverty = c(ame_probit$std.error[ame_probit$term == "poverty"],
ame_lpm$std.error[ame_lpm$term == "poverty"],
ame_re$std.error[ame_re$term == "poverty"])
)
ame_comp[, 2:3] <- round(ame_comp[, 2:3], 5)
knitr::kable(ame_comp, caption = "Average Marginal Effects across binary CRE models.")
```
All three specifications give an average marginal effect of poverty close to $-0.064$, with a standard error of about $0.034$. In this example, CRE pooled probit gives us the partial effect without imposing the serial-independence assumption of the nonlinear random-effects model.
## What we learn
For a balanced linear panel, adding the unit averages $\bar{X}_i$ to pooled OLS or random effects recovers the fixed-effects coefficients. In binary and count models, CRE is useful because we can calculate average marginal effects and avoid estimating a fixed effect for every unit. Pooled estimation with unit-clustered standard errors is especially attractive when we do not want to assume serial independence.