---
title: "Fixed or Random Effect, or Both?"
date: "2019-05-23"
---
## Panel data
With panel data — repeated observations over time or observations clustered at a higher level — the standard choice is between fixed and random effects. Consider
$$ y_{it} = \beta_0 + \beta_1 x_{it} + c_i + \epsilon_{it} $$ {#eq-fixed-random-effect-1}
for individuals $i = 1, \ldots, n$ measured at $t = 1, \ldots, T$, where $c_i$ is an unobserved time-invariant individual effect. The two approaches differ in how they handle $c_i$.
Fixed effects eliminate $c_i$ by demeaning or by including individual dummies (the two are equivalent in a linear model). In nonlinear models only Poisson avoids the incidental-parameter problem; the standard solution for other models is conditional likelihood (conditional logit, for example), which absorbs the fixed effects in the likelihood. When no conditional likelihood exists, the incidental-parameter bias shrinks with panel depth but does not vanish.
Random effects treat $c_i$ as part of the error term. This is more efficient but requires $x_{it}$ to be uncorrelated with $c_i$ — an assumption economists are generally skeptical of.
## Time-invariant variables
Sometimes we want the effect of a time-invariant variable, so the model becomes
$$ y_{it} = \beta_0 + \beta_1 x_{it} + c_i + \gamma z_i+ \epsilon_{it} $$ {#eq-fixed-random-effect-2}
Fixed effects cannot identify $\gamma$ because $z_i$ is perfectly collinear with $c_i$. Random effects can still estimate both, treating $z_i$ as another covariate.
## Between-within model
The textbook advice is a Hausman test: random effects are more efficient when the uncorrelatedness assumption holds; if not, fixed effects are still consistent. The between-within model (BW) sidesteps the choice. [Neuhaus and Kalbfleisch (1998)](https://www.ncbi.nlm.nih.gov/pubmed/9629647) introduced the BW estimator,
$$ y_{it} = \beta_0 + \beta_1 (x_{it} - \bar x_i) + \beta_2 \bar x_i + c_i + \gamma z_i+ \epsilon_{it} $$ {#eq-fixed-random-effect-3}
Here $\beta_1$ is the within effect — identical to the fixed-effect coefficient. $\beta_2$ is the between effect (mean of $x$ on mean of $y$). $\gamma$ is the effect of the time-invariant variable.
The other specification of BW estimator is
$$ y_{it} = \beta_0 + \beta_1 x_{it} + \beta_2 \bar x_i + c_i + \gamma z_i+ \epsilon_{it} $$ {#eq-fixed-random-effect-4}
This is the same model reparametrized. $\beta_1$ is unchanged, and $\beta_2$ is now the between-minus-within difference (substitute $\beta_1(x_{it}-\bar x_i) + \beta_2 \bar x_i = \beta_1 x_{it} + (\beta_2-\beta_1)\bar x_i$ to see this). This is the "contextual model"; $\beta_2$ is the contextual effect and acts as an embedded Hausman test.
The BW model has two advantages. It recovers the fixed-effect coefficient while estimating time-invariant covariates. And it extends naturally to cross-level interactions, random slopes, and other multilevel structures. Implementation is straightforward: run a random-effects model on either equation above.
### Is this just Mundlak?
Yes. @eq-fixed-random-effect-4 *is* the Mundlak device: Mundlak's regression puts
$x_{it}$ and $\bar x_i$ on the right-hand side of a random-effects model, which is
exactly what that equation does. So the between-within model and the Mundlak
device are one model in two parameterisations, fitted by the same random-effects
command. It is worth doing the algebra, because the identity is easy to state
loosely and easy to get wrong on the inference side.
Write $d_{it} = x_{it} - \bar x_i$ for the within deviation, and label the two
specifications by their coefficients: $b_1, b_2$ for the hybrid form
@eq-fixed-random-effect-3 and $p_1, p_2$ for the contextual form
@eq-fixed-random-effect-4. Expanding the hybrid regressor block,
$$
b_1 (x_{it} - \bar x_i) + b_2 \bar x_i = b_1 x_{it} + (b_2 - b_1)\bar x_i ,
$$ {#eq-fixed-random-effect-5}
so $p_1 = b_1$ and $p_2 = b_2 - b_1$. Substitution shows the map, but the reason
behind it is worth stating more strongly: the two regressor matrices are related
by an invertible linear transformation. With $X = [\,x_{it}, \bar x_i\,]$ and
$W = [\,d_{it}, \bar x_i\,]$,
$$
W = XA, \qquad
A = \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}, \qquad
A^{-1} = \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix},
$$ {#eq-fixed-random-effect-6}
and $\det A = 1$, so the two designs span the same column space. Everything that
depends only on that space is numerically identical: fitted values, residuals,
the likelihood, the estimated variance components, and the GLS quasi-demeaning
weight $\theta$. Only the coordinates change, by
$b = A^{-1}p$, and the covariance matrix by $\text{Var}(b) = A^{-1}\text{Var}(p)A^{-1\prime}$.
Componentwise that is
$$
\text{Var}(p_1) = \text{Var}(b_1), \qquad
\text{Var}(p_2) = \text{Var}(b_1) + \text{Var}(b_2) + 2\,\text{Cov}(b_1,b_2),
$$ {#eq-fixed-random-effect-7}
which already says something useful: the standard error printed next to the
between effect is not the standard error you want for the contextual effect.
#### Why the within coefficient is the fixed-effect estimator
This part needs a proof rather than a substitution. The relevant fact is that
$\sum_t d_{it} = 0$ within every unit by construction, so $d$ is orthogonal to
any group-invariant column: for $v_i$ constant within $i$,
$$
\sum_{i,t} d_{it} v_i = \sum_i v_i \sum_t d_{it} = 0 ,
$$ {#eq-fixed-random-effect-8}
which covers the intercept, $\bar x_i$ and $z_i$ alike.
Random effects is OLS after quasi-demeaning, $\tilde w_{it} = w_{it} - \theta \bar w_i$
with $\theta = 1 - \sqrt{\sigma^2_\epsilon / (\sigma^2_\epsilon + T\sigma^2_c)}$.
Apply it to the hybrid regressors. Because $\bar d_i = 0$ we get
$\tilde d_{it} = d_{it}$: the within-deviation regressor passes through the GLS
transform untouched, whatever $\theta$ is. The remaining regressors are
group-invariant and stay so after scaling by $1-\theta$. Orthogonality therefore
survives the transform, and by Frisch-Waugh-Lovell
$$
b_1 = \frac{\sum d_{it}\tilde y_{it}}{\sum d_{it}^2}
= \frac{\sum d_{it} y_{it}}{\sum d_{it}^2}
= \frac{\sum (x_{it}-\bar x_i)(y_{it}-\bar y_i)}{\sum (x_{it}-\bar x_i)^2},
$$ {#eq-fixed-random-effect-9}
the second and third equalities both by @eq-fixed-random-effect-8. That last
expression is the within estimator. Note that $\theta$ has dropped out, so the
identity holds for *any* $\theta$ — for random effects, for pooled OLS
($\theta = 0$) and for fixed effects ($\theta = 1$). That is why running a
random-effects model on either equation returns the fixed-effect coefficient. It
also gives $\text{Cov}(b_1,b_2) = 0$ exactly, which combined with
@eq-fixed-random-effect-7 yields
$\text{Var}(p_2) = \text{Var}(b_1) + \text{Var}(b_2)$ and
$\text{Cov}(p_1,p_2) = -\text{Var}(b_1)$.
All of this is checkable on the chapter's own data. We take a single regressor,
weeks worked, fit both parameterisations by `plm`, and compare against the
within estimator and against the variance map:
```{r}
#| label: bw-mundlak-identity
#| echo: true
#| message: false
#| warning: false
library(plm); library(panelr); data("WageData")
d <- WageData
d$xbar <- ave(d$wks, d$id) # unit mean
d$xdev <- d$wks - d$xbar # within deviation
fe <- plm(lwage ~ wks, data = d, index = c("id", "t"), model = "within")
hyb <- plm(lwage ~ xdev + xbar, data = d, index = c("id", "t"), model = "random")
mun <- plm(lwage ~ wks + xbar, data = d, index = c("id", "t"), model = "random")
b <- coef(hyb); p <- coef(mun); Vh <- vcov(hyb); Vm <- vcov(mun)
data.frame(
check = c("p1 - b1", "p2 - (b2 - b1)", "b1 - FE within",
"sum(xdev * xbar)", "Cov(b1, b2)",
"Var(p2) - [Var(b1) + Var(b2)]", "Cov(p1, p2) + Var(b1)",
"theta(hybrid) - theta(Mundlak)"),
value = c(p["wks"] - b["xdev"],
p["xbar"] - (b["xbar"] - b["xdev"]),
b["xdev"] - coef(fe)["wks"],
sum(d$xdev * d$xbar),
Vh["xdev", "xbar"],
Vm["xbar", "xbar"] - (Vh["xdev", "xdev"] + Vh["xbar", "xbar"]),
Vm["wks", "xbar"] + Vh["xdev", "xdev"],
hyb$ercomp$theta - mun$ercomp$theta),
row.names = NULL)
```
Every entry is zero to machine precision. The two parameterisations are the same
fit, the within coefficient is the fixed-effect estimate, the hybrid regressors
are orthogonal, and the variance map holds as derived.
#### So why keep both
The fit is one object, but the parameterisation settles three practical things.
*Which quantity is printed.* $b_2$ is the between effect and $p_2$ is between
minus within. They are different parameters with non-interchangeable standard
errors, and whichever form you fit, the other costs a `lincom`.
*Which hypothesis is one line.* $p_2 = 0$ says between equals within, which is
exactly the random-effects assumption that the regressors are uncorrelated with
$c_i$. The Mundlak form prints that as a single $t$; in hybrid form it is the
contrast $b_1 = b_2$. The two tests agree numerically, but only one is printed.
This is the Mundlak form of the Hausman test, asymptotically equivalent to the
classical version and with the practical advantage of working directly with
cluster-robust standard errors. The [next chapter](mundlak-device.qmd) develops
it and its two-way extension for difference-in-differences.
*What the second coefficient is for.* In multilevel work $p_2$ is not a
diagnostic but the estimand: holding a pupil's own status fixed, what does the
school mean do? The hybrid's $b_2$ answers a different question, the total
between-unit association. Two literatures ask two questions of one regression,
and each named the parameterisation that printed its answer. That is the honest
reason both names survive.
There is also a reason to prefer the hybrid form in extensions. Its two
regressors are exactly orthogonal, where $x_{it}$ and $\bar x_i$ are correlated.
That does not change the fit, but if you want a random slope on the within
component you need the deviation as its own column; a random slope on $x_{it}$
would mix within and between variation.
#### Where the equivalence stops
The reparameterisation is exact for anything with a linear index, logit and
Poisson mixed models included, since $XA$ spans the same space and the index is
unchanged. The two names stay interchangeable there.
The fixed-effect identity does not carry over. The proof above used two
linear-model facts: that random effects is OLS on quasi-demeaned data, and that
$d$ is orthogonal to group constants. A logit or Poisson mixed model has no such
transform, the unit effects do not sweep out, and the hybrid within coefficient
is not the conditional-logit estimate but an approximation to it. Keep the two
claims separate: the parameterisations are always the same model, while "within
coefficient equals fixed effects" is a linear-model result. The last section of
this chapter returns to that point.
## BW model in R
R's [`panelr`](https://panelr.jacob-long.com/articles/wbm.html) package implements BW models directly. The data are `WageData`, a balanced panel of 595 men observed over 7 years. We model log wage on weeks worked, union membership, marital status and occupation, all time-varying, plus two time-invariant covariates — black and female — that a fixed-effect model could not touch.
```{r}
#| echo: true
#| message: false
library(panelr)
data("WageData")
wages <- panel_data(WageData, id = id, wave = t)
model1 <- wbm(lwage ~ wks + union + ms + occ | blk + fem, data = wages)
summary(model1)
```
The within estimates are 0.00 for weeks, 0.06 for union, $-0.08$ for married
and $-0.08$ for occupation. The between estimates for the same four variables
are 0.01, 0.15, 0.17 and $-0.41$ — two to five times larger, and for `ms` the
opposite sign. That gap is the reason the fixed-versus-random choice matters
here: a random-effects model without the decomposition would blend the two and
report neither. The time-invariant covariates come out at $-0.15$ for black and
$-0.32$ for female, quantities a fixed-effect model cannot produce at all.
Compare with `lfe`.
```{r}
#| echo: true
#| message: false
library(lfe)
model2 <- felm(lwage ~ wks + union + ms + occ | id, data = wages)
summary(model2)
```
`felm` gives 0.001083, 0.064320, $-0.082905$ and $-0.077507$ — the same four
numbers `panelr` reported as within effects, identical to machine precision
(`WageData` is a balanced panel: 7 waves, 595 individuals). That is the BW
model's selling point: it reproduces the fixed-effect estimates exactly while
also estimating the time-invariant covariates `blk` and `fem`, which `felm`
cannot. `lfe` supports clustered standard errors, which `panelr` does not.
```{r}
#| echo: true
#| message: false
model3 <- felm(lwage ~ wks + union + ms + occ | id | 0 | id, data = wages)
summary(model3)
```
Clustering on `id` leaves every coefficient untouched and changes the standard
errors substantially: `union` goes from 0.0254 to 0.0409, `ms` from 0.0322 to
0.0474, `wks` from 0.00102 to 0.00133. Two of the four results do not survive.
`union` falls from $p = 0.011$ to $p = 0.117$ and `ms` from $0.010$ to $0.081$;
only `occ` stays significant at 5%. With 7 observations per man the residuals
are serially correlated, and the unclustered standard errors were treating 4,165
observations as though they were independent draws.
## BW model in Stata
Stata has no dedicated BW command, but `xtreg` suffices. We switch to `nlswork`,
an unbalanced panel of 4,710 women over 28,510 person-years, and regress log
wage on age.
```{r}
#| echo: false
#| message: false
library(Statamarkdown)
```
```{stata}
*| cache: true
webuse nlswork
xtset idcode
xtreg ln_w age, fe cluster(idcode)
```
The fixed-effect estimate is $.0181349$ with a clustered standard error of
$.0006099$: within a woman, each additional year of age raises log wage by about
1.8%. Keep that number in view — the BW model has to reproduce it.
Generate the group mean and run the BW estimation.
```{stata}
*| cache: true
webuse nlswork
xtset idcode
bysort idcode: center age, prefix(d) mean(m)
xtreg ln_w dage mage i.race, re cluster(idcode)
```
In this BW model, we use the centered (within-group-demeaned) `dage` together with the group mean `mage`, which gives the standard within-between (Mundlak) decomposition: the coefficient on `dage`, .0181349, is exactly the fixed-effect ("within") coefficient on age reported above — the same seven digits, with a standard error of .00061 against .0006099 — and the coefficient on `mage`, .022558, is directly the between effect. Women who are older on average earn about 2.3% more per year of average age, against a within-woman return of 1.8%; the between effect is larger because it also carries cohort and selection differences that the within estimate sweeps out. And we have the effect of time-invariant covariate race estimated: Black women earn 11.9% less ($-.1190246$, SE .0127), while the "Other" category is not distinguishable from white ($.0975$, $p = 0.114$). The advantage of using xtreg is that we have clustered standard errors implemented.
Note: if we instead ran `xtreg ln_w age mage i.race, re cluster(idcode)` — using the raw (uncentered) `age` together with `mage` — we would get the *contextual model* instead: the coefficient on `age` stays .0181, but the coefficient on `mage` becomes .0044, the "contextual effect" (the *additional* between-effect on top of the within effect). The two models are algebraically related: contextual-model mage = Mundlak-model mage − Mundlak-model dage (here, .0226 − .0181 ≈ .0044).
## BW model in non-linear models
[Allison](https://statisticalhorizons.com/between-within-contextual-effects) applies the BW model to binary outcomes. The bias relative to conditional logit is an open question, but if the linear probability model is a reasonable approximation the BW approach should perform similarly for binary panel data.