Module 6: Panel Data

Fixed Effects, Random Effects, and the Power of Watching Over Time

Author

Austin Tucker

Published

July 23, 2026

1 Panel Data: Watching Over Time

Suppose a district health office runs an annual school health screening. Every year, the same children have a finger-prick blood draw and their hemoglobin (hgb, in g/dL) is recorded. Some children are enrolled in an iron supplementation program (treat), and enrollment can switch on or off from year to year. After five years, we have five observations per child.

Data in which we observe the same individuals at multiple points in time are called longitudinal or panel data. The term “individual” is generic — children here, but the units could be clinics, counties, or countries. Panel data contrast with three other data types:

  • Cross-sectional data: each individual observed exactly once at a single point in time.
  • Repeated cross-sections: fresh (generally random) samples drawn from the same population at different times — different individuals each round.
  • Time series data: a small number of units (often one) observed at high frequency over many periods.

A panel is balanced if every individual is observed in every period, and unbalanced otherwise. Our screening panel is balanced: 1,000 children, 5 visits each.

1.1 Wide vs. long format

The same panel can be stored two ways. In wide format, each child is one row and each period is a column. In long format, each child-period is its own row:

Table 1: Wide format — one row per child
id hgb_1 hgb_2 hgb_3 hgb_4 hgb_5
1 10.2 9.7 13.8 10.0 11.2
2 9.0 11.1 10.0 9.3 13.6
Table 2: Long format — one row per child-period
id period treat hgb
1 1 0 10.2
1 2 1 9.7
\(\vdots\) \(\vdots\) \(\vdots\) \(\vdots\)
2 1 0 9.0
2 2 0 11.1

Panel estimators (and essentially all modern software) want long format, which is how our simulated population pop is stored:

1.2 Within vs. between variation

When the same child is measured repeatedly, hemoglobin varies for two distinct reasons:

  • Between variation: some children are persistently sicker or healthier than others — child 412 sits low every single year, child 87 sits high.
  • Within variation: a given child’s hemoglobin moves around her own average from year to year — a bout of illness, a growth spurt, starting supplementation.

We can decompose the total variation in hgb into these two pieces:

Roughly half of the spread in hemoglobin comes from which child you are (between) and half from year-to-year movement within a child (within). Keep this decomposition in mind — every estimator in this module is defined by which slice of variation it uses.

1.3 Exercise 1: Spaghetti plot of hemoglobin trajectories

The single most useful first look at panel data is a spaghetti plot: one line per individual, traced over time. Plot the trajectories for the first 20 children. Fill in the filter cutoff and the aesthetic that tells ggplot which observations belong to the same line.

NoteHint 1

You want the first 20 children, so filter to id <= 20.

NoteHint 2

The group aesthetic tells ggplot which rows to connect into one line — here, all rows sharing the same id.

TipSolution
kids <- pop %>% filter(id <= 20) ggplot(kids, aes(x = period, y = hgb, group = id, color = factor(id))) + geom_line(alpha = 0.8, linewidth = 0.7) + geom_point(size = 1.6) + scale_color_manual(values = rep(c(terracotta, sage, slate_blue, dusty_gold), 5)) + labs( title = "Hemoglobin Trajectories, 20 Children Over 5 Years", subtitle = "Each line is one child. Lines rarely cross their neighbors far away.", x = "Screening year", y = "Hemoglobin (g/dL)" ) + theme_house() + theme(legend.position = "none")

kids <- pop %>% filter(id <= 20)

ggplot(kids, aes(x = period, y = hgb, group = id, color = factor(id))) +
  geom_line(alpha = 0.8, linewidth = 0.7) +
  geom_point(size = 1.6) +
  scale_color_manual(values = rep(c(terracotta, sage, slate_blue, dusty_gold), 5)) +
  labs(
    title = "Hemoglobin Trajectories, 20 Children Over 5 Years",
    subtitle = "Each line is one child. Lines rarely cross their neighbors far away.",
    x = "Screening year",
    y = "Hemoglobin (g/dL)"
  ) +
  theme_house() +
  theme(legend.position = "none")

The lines wiggle year to year, but a child who starts near the bottom tends to stay near the bottom. Which kind of variation does the vertical spread between lines represent? Which kind does the wiggle along each line represent?

Discussion: The vertical spread between lines is between variation — persistent differences in child-level health (in our simulation, \(\alpha_i\)). The wiggle along each line is within variation — year-to-year movement around a child’s own average. The persistence of the ordering (low starters stay low) is visual evidence that a stable child-specific component is at work, and that component is about to cause us serious trouble.

2 The Panel Population Model

We index children by \(i\) and screening years by \(t\). The population model linking hemoglobin to supplementation is:

\[ \text{hgb}_{it} = \beta_0 + \beta_1 \, \text{treat}_{it} + \alpha_i + \varepsilon_{it} \]

Two error components deserve a careful look:

  • \(\alpha_i\) is unobserved child-level heterogeneity: everything stable about child \(i\) that affects hemoglobin — genetics, chronic parasitic infection, household diet, poverty. It carries only an \(i\) subscript: it varies between children but is constant within a child over time. You can read \(\beta_0 + \alpha_i\) as child \(i\)’s personal intercept.
  • \(\varepsilon_{it}\) is the idiosyncratic error: year-to-year shocks (an infection this spring, a measurement blip) that vary both between and within children.

In our simulated screening program, the true supplementation effect is \(\beta_1 = 0.8\) g/dL. But there is a catch, and it is the heart of this module.

NoteThe selection story: why \(\alpha_i\) is correlated with treatment

Who gets enrolled in an iron supplementation program? Not a random sample of children — the sicker ones. The school nurse looks at pale, tired, frequently absent children and refers them. In our DGP, the probability of treatment is \(\text{plogis}(-1.0 \, \alpha_i + 0.4\,(t - 3))\): children with low health endowments (\(\alpha_i < 0\)) are more likely to be treated. Across children, the correlation between a child’s average treatment rate and her \(\alpha_i\) is about \(-0.76\).

This is confounding by indication, panel edition: treatment is targeted at exactly the children whose untreated outcomes would be worst. Formally, \(\text{Cov}(\text{treat}_{it}, \alpha_i) \neq 0\).

Because we built this world ourselves, pop carries a column alpha that no real dataset would ever contain. We will occasionally peek at it to check our reasoning — a privilege of simulation, not of practice.

One more feature of the DGP, planted deliberately: hemoglobin drifts up by 0.15 g/dL per year for everyone (better school meals, say), and the supplementation program expands over time. Park that thought — it will not matter until Section 6.

3 Why Pooled OLS Fails

The most naive approach ignores the panel structure entirely: stack all 5,000 child-years and run OLS as if they were 5,000 unrelated observations. This is pooled OLS:

\[ \text{hgb}_{it} = \beta_0 + \beta_1 \, \text{treat}_{it} + v_{it}, \qquad v_{it} = \alpha_i + \varepsilon_{it} \]

Since we do not model \(\alpha_i\), it gets folded into a composite error \(v_{it}\). For pooled OLS to be unbiased we need the zero conditional mean assumption on the composite error:

\[ E[v_{it} \mid \text{treat}_{it}] = 0 \quad \Longleftrightarrow \quad \text{Cov}(\text{treat}_{it}, \alpha_i) = 0 \;\text{ and }\; \text{Cov}(\text{treat}_{it}, \varepsilon_{it}) = 0 \]

The first condition is exactly what our selection story violates.

NoteThis is Module 5 wearing a new outfit

Pooled OLS here is the short regression from Module 5, and \(\alpha_i\) is the omitted variable. Apply the OVB formula: the omitted \(\alpha_i\) raises hemoglobin (\(\beta_{\alpha} = +1\) by construction) and is negatively related to treatment (\(\delta < 0\)), so the bias \(\beta_{\alpha} \cdot \delta\) is negative — large enough here to drag the estimate below zero. The only novelty is the fix: because \(\alpha_i\) is time-invariant, panel data let us eliminate it without ever observing it.

3.1 Exercise 2: Pooled OLS on the screening panel

Draw a random sample of 200 children (keeping all 5 years for each — we sample children, not child-years) and fit pooled OLS of hemoglobin on treatment. The true effect is \(+0.8\). Brace yourself.

NoteHint 1

The dependent variable is hgb; the independent variable is treat.

NoteHint 2

The formula should be: hgb ~ treat

TipSolution
set.seed(20231005) ids_samp <- sample(unique(pop$id), 200) samp <- pop %>% filter(id %in% ids_samp) mod_pooled <- lm(hgb ~ treat, data = samp) summary(mod_pooled)
set.seed(20231005)
ids_samp <- sample(unique(pop$id), 200)
samp <- pop %>% filter(id %in% ids_samp)

mod_pooled <- lm(hgb ~ treat, data = samp)
summary(mod_pooled)

The true effect of supplementation is \(+0.8\) g/dL. What did pooled OLS give you, and how would a naive analyst interpret it?

Discussion: The pooled estimate is about \(-0.52\) (SE \(\approx 0.10\)) — not just biased but wrong-signed, and “highly significant.” A naive analyst would conclude that iron supplementation lowers hemoglobin. What OLS is really reporting is that treated child-years belong disproportionately to persistently anemic children: the comparison of treated vs. untreated observations is mostly a comparison of sick kids vs. healthy kids. Statistical significance offers no protection whatsoever against bias — the estimator is precisely estimating the wrong thing, and no sample size will rescue it.

4 Fixed Effects: Exploiting Only Within Variation

The fixed effects insight: \(\alpha_i\) is constant within a child, so any transformation that subtracts away each child’s own average also subtracts away \(\alpha_i\)without our ever observing it.

4.1 The within transformation

Average the population model over child \(i\)’s five years:

\[ \overline{\text{hgb}}_i = \beta_0 + \beta_1 \, \overline{\text{treat}}_i + \alpha_i + \bar{\varepsilon}_i \]

Subtract this from the original equation — demeaning, or the within transformation:

\[ (\text{hgb}_{it} - \overline{\text{hgb}}_i) = \beta_1 (\text{treat}_{it} - \overline{\text{treat}}_i) + (\underbrace{\alpha_i - \alpha_i}_{=\,0}) + (\varepsilon_{it} - \bar{\varepsilon}_i) \]

The nuisance term \(\alpha_i\) is annihilated because it equals its own within-child mean. What survives is a regression run entirely on within variation: each child’s deviations from her own averages. Hence the fixed effects estimator’s other name, the within estimator. Identification now asks a much gentler question: in the years a given child was supplemented, was her hemoglobin higher than her own norm? Persistent between-child differences — the whole source of our bias — never enter.

The assumptions are correspondingly weaker than pooled OLS’s. We need \(\text{Cov}(\text{treat}_{it}, \varepsilon_{is}) = 0\) (strict exogeneity of the idiosyncratic error), but \(\text{Cov}(\text{treat}_{it}, \alpha_i)\) can be anything at all — the sicker-kids-get-treated selection is fully permitted.

4.2 Exercise 3: Fixed effects by manual demeaning

Demean hemoglobin and treatment within each child using group_by() + mutate(), then run OLS on the demeaned variables.

NoteHint 1

Group by the child identifier id, so that mean() computes each child’s own 5-year average.

NoteHint 2
  • group_by(id)
  • hgb - mean(hgb)
  • treat - mean(treat)
TipSolution
set.seed(20231005) ids_samp <- sample(unique(pop$id), 200) samp <- pop %>% filter(id %in% ids_samp) samp_dm <- samp %>% group_by(id) %>% mutate( hgb_dm = hgb - mean(hgb), treat_dm = treat - mean(treat) ) %>% ungroup() mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm) summary(mod_within)
set.seed(20231005)
ids_samp <- sample(unique(pop$id), 200)
samp <- pop %>% filter(id %in% ids_samp)

samp_dm <- samp %>%
  group_by(id) %>%
  mutate(
    hgb_dm   = hgb - mean(hgb),
    treat_dm = treat - mean(treat)
  ) %>%
  ungroup()

mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm)
summary(mod_within)

Compare the coefficient on treat_dm to the pooled estimate from Exercise 2 and to the truth (\(+0.8\)). Also: why is the intercept essentially zero?

Discussion: The within estimate is about \(+0.84\) — the sign flips back and we land close to the true \(0.8\). Subtracting each child’s own mean removed \(\alpha_i\) and with it the confounding by indication. The intercept is numerically zero because demeaned variables have mean zero by construction. (Two footnotes: the small remaining gap above \(0.8\) is partly sampling noise and partly the secular time trend we planted — resolved in Section 6. And the lm() standard errors here are slightly too small, since plain OLS does not know we “used up” 200 child means in the demeaning; dedicated panel software applies the degrees-of-freedom correction and, ideally, clusters standard errors by child.)

4.3 Least squares dummy variables (LSDV)

An algebraically equivalent route: instead of demeaning, estimate every child’s intercept explicitly by adding a dummy variable for each child:

\[ \text{hgb}_{it} = \beta_1 \, \text{treat}_{it} + \sum_{j} \delta_j \, \mathbb{1}[i = j] + \varepsilon_{it} \]

Each child’s dummy absorbs all time-invariant variation specific to that child, so \(\beta_1\) is again identified purely from within variation. By the Frisch–Waugh–Lovell theorem, LSDV and the within estimator produce the identical slope on treat.

4.4 Exercise 4: Fixed effects via LSDV

Fit the same model by adding factor(id) to a plain lm(), then verify the slope matches Exercise 3 exactly.

NoteHint 1

Wrap the child identifier in factor() so lm() creates one dummy variable per child rather than treating id as a number.

NoteHint 2

The formula should be: hgb ~ treat + factor(id)

TipSolution
set.seed(20231005) ids_samp <- sample(unique(pop$id), 200) samp <- pop %>% filter(id %in% ids_samp) # LSDV: one dummy per child mod_lsdv <- lm(hgb ~ treat + factor(id), data = samp) # Within estimator from Exercise 3, for comparison samp_dm <- samp %>% group_by(id) %>% mutate(hgb_dm = hgb - mean(hgb), treat_dm = treat - mean(treat)) %>% ungroup() mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm) cat("LSDV coefficient on treat: ", round(coef(mod_lsdv)["treat"], 4), "\n") cat("Within coefficient on treat:", round(coef(mod_within)["treat_dm"], 4), "\n") cat("Number of coefficients in LSDV model:", length(coef(mod_lsdv)), "\n")
set.seed(20231005)
ids_samp <- sample(unique(pop$id), 200)
samp <- pop %>% filter(id %in% ids_samp)

# LSDV: one dummy per child
mod_lsdv <- lm(hgb ~ treat + factor(id), data = samp)

# Within estimator from Exercise 3, for comparison
samp_dm <- samp %>%
  group_by(id) %>%
  mutate(hgb_dm = hgb - mean(hgb), treat_dm = treat - mean(treat)) %>%
  ungroup()
mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm)

cat("LSDV coefficient on treat:  ", round(coef(mod_lsdv)["treat"], 4), "\n")
cat("Within coefficient on treat:", round(coef(mod_within)["treat_dm"], 4), "\n")
cat("Number of coefficients in LSDV model:", length(coef(mod_lsdv)), "\n")

The two slopes agree to every decimal place, yet the LSDV model estimated 201 coefficients while the within regression estimated 2. If you compared their \(R^2\) values, which would be larger, and would that make LSDV a “better” model?

Discussion: The slopes are identical by the Frisch–Waugh–Lovell theorem — demeaning and dummying-out are the same projection. The LSDV \(R^2\) would be far larger, but only because 200 child dummies mechanically soak up all between-child variation; it says nothing about how well we explain within-child changes, which is the comparison we actually care about. This is also a computational warning: with thousands of units, factor(id) builds an enormous design matrix. Demeaning scales gracefully; production packages (see below) use it under the hood.

4.5 What fixed effects cost you

The within transformation is a scorched-earth tactic: it kills \(\alpha_i\), but it kills every other time-invariant variable along with it. Our data contain rural, an observed time-invariant characteristic. Watch what happens:

Demeaned rural has zero variance, and in LSDV it is perfectly collinear with the child dummies, so lm() returns NA. The general lesson:

ImportantThe fixed effects bargain

FE buys robustness to any correlation between treatment and time-invariant confounders. The price:

  • Time-invariant covariates are swept away. You cannot estimate the effect of rural, sex, or baseline genotype in an FE model — their coefficients are absorbed into the child intercepts.
  • Only within variation identifies \(\beta_1\). A child treated in all 5 years, or none, contributes nothing to the estimate. If treatment rarely changes within units, FE estimates are imprecise.
  • Only time-invariant confounding is removed. A time-varying omitted variable (a diarrheal outbreak that both lowers hemoglobin and triggers enrollment mid-panel) still biases FE.

5 Random Effects: The Efficiency Bargain

Suppose — hypothetically — that supplementation had been assigned by lottery, so \(\text{Cov}(\text{treat}_{it}, \alpha_i) = 0\). Then pooled OLS would be unbiased, and throwing away all between variation (as FE does) would be wasteful. But pooled OLS would still have a problem: serial correlation. Because \(\alpha_i\) sits in the composite error \(v_{it} = \alpha_i + \varepsilon_{it}\) of every one of child \(i\)’s five observations, her errors are correlated over time:

\[ \text{Corr}(v_{it}, v_{is}) = \frac{\sigma_\alpha^2}{\sigma_\alpha^2 + \sigma_\varepsilon^2} \quad \text{for } t \neq s \]

The random effects (RE) estimator models this correlation structure and reweights accordingly. Mechanically, RE applies quasi-demeaning: subtract only a fraction \(\lambda\) of each child’s mean,

\[ (\text{hgb}_{it} - \lambda\, \overline{\text{hgb}}_i) = \beta_0 (1 - \lambda) + \beta_1 (\text{treat}_{it} - \lambda\, \overline{\text{treat}}_i) + (v_{it} - \lambda \bar{v}_i) \]

with

\[ \lambda = 1 - \sqrt{\frac{\sigma_\varepsilon^2}{\sigma_\varepsilon^2 + T \sigma_\alpha^2}} \]

where \(T\) is the number of periods per child. RE is a compromise dial:

  • \(\lambda = 0\): no demeaning — RE collapses to pooled OLS.
  • \(\lambda = 1\): full demeaning — RE collapses to fixed effects.
  • In between, RE blends within and between variation, buying efficiency (smaller standard errors) by using both.

The catch is the fine print of the bargain: RE leaves \(\alpha_i\) in the error term, so it is consistent only if \(\text{Cov}(\text{treat}_{it}, \alpha_i) = 0\) — the very assumption our sicker-kids-get-treated world violates. In our DGP, RE inherits a diluted version of pooled OLS’s bias.

The variances \(\sigma_\varepsilon^2\) and \(\sigma_\alpha^2\) are unknown, so in practice \(\lambda\) is estimated from the data — making RE a feasible generalized least squares (FGLS) estimator. We estimate \(\sigma_\varepsilon^2\) from the within (FE) residuals and \(\sigma_\alpha^2\) from the between regression (child means on child means) residuals, using \(\sigma_\alpha^2 = \sigma_B^2 - \sigma_\varepsilon^2 / T\).

NotePedagogical FGLS vs. production tools

The RE estimator we build below is a deliberately simplified, hand-rolled FGLS — the same Swamy–Arora logic real packages use, laid bare so you can see every moving part. In real analyses, use plm::plm(..., model = "random") or the fixest package (and lme4::lmer for the mixed-model formulation); they handle unbalanced panels, degrees-of-freedom corrections, and proper (clustered) standard errors. The same advice applies to fixed effects: plm(..., model = "within") or fixest::feols(y ~ x | id) are the production tools — we demeaned by hand because hand-rolling teaches more, and because they run anywhere (including this browser).

5.1 Exercise 5: Random effects by hand (FGLS)

Build the RE estimator in three steps: estimate the variance components, compute \(\lambda\), quasi-demean, and fit OLS. Fill in the \(\lambda\) formula and the quasi-demeaning.

NoteHint 1

In the \(\lambda\) formula, the denominator inside the square root is \(\sigma_\varepsilon^2 + T \sigma_\alpha^2\): use T_periods and sigma_a2.

NoteHint 2
  • lambda <- 1 - sqrt(sigma_e2 / (sigma_e2 + T_periods * sigma_a2))
  • Quasi-demeaning subtracts lambda * mean(hgb) and lambda * mean(treat).
TipSolution
set.seed(20231005) ids_samp <- sample(unique(pop$id), 200) samp <- pop %>% filter(id %in% ids_samp) T_periods <- 5 n_kids <- length(unique(samp$id)) # Step 1a: sigma_e^2 from the within (FE) residuals samp_dm <- samp %>% group_by(id) %>% mutate(hgb_dm = hgb - mean(hgb), treat_dm = treat - mean(treat)) %>% ungroup() mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm) sigma_e2 <- sum(residuals(mod_within)^2) / (nrow(samp) - n_kids - 1) # Step 1b: sigma_alpha^2 from the between regression (child means) between <- samp %>% group_by(id) %>% summarize(hgb_bar = mean(hgb), treat_bar = mean(treat), .groups = "drop") mod_between <- lm(hgb_bar ~ treat_bar, data = between) sigma_B2 <- sum(residuals(mod_between)^2) / (n_kids - 2) sigma_a2 <- max(sigma_B2 - sigma_e2 / T_periods, 0) # Step 2: lambda = 1 - sqrt( sigma_e^2 / (sigma_e^2 + T * sigma_alpha^2) ) lambda <- 1 - sqrt(sigma_e2 / (sigma_e2 + T_periods * sigma_a2)) # Step 3: quasi-demean by lambda and fit OLS samp_q <- samp %>% group_by(id) %>% mutate( hgb_q = hgb - lambda * mean(hgb), treat_q = treat - lambda * mean(treat) ) %>% ungroup() mod_re <- lm(hgb_q ~ treat_q, data = samp_q) cat("Estimated sigma_e^2: ", round(sigma_e2, 3), "\n") cat("Estimated sigma_alpha^2:", round(sigma_a2, 3), "\n") cat("Estimated lambda: ", round(lambda, 3), "\n") cat("RE estimate of beta_1: ", round(coef(mod_re)["treat_q"], 3), "\n")
set.seed(20231005)
ids_samp <- sample(unique(pop$id), 200)
samp <- pop %>% filter(id %in% ids_samp)

T_periods <- 5
n_kids    <- length(unique(samp$id))

# Step 1a: sigma_e^2 from the within (FE) residuals
samp_dm <- samp %>%
  group_by(id) %>%
  mutate(hgb_dm = hgb - mean(hgb), treat_dm = treat - mean(treat)) %>%
  ungroup()
mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm)
sigma_e2 <- sum(residuals(mod_within)^2) / (nrow(samp) - n_kids - 1)

# Step 1b: sigma_alpha^2 from the between regression (child means)
between <- samp %>%
  group_by(id) %>%
  summarize(hgb_bar = mean(hgb), treat_bar = mean(treat), .groups = "drop")
mod_between <- lm(hgb_bar ~ treat_bar, data = between)
sigma_B2  <- sum(residuals(mod_between)^2) / (n_kids - 2)
sigma_a2  <- max(sigma_B2 - sigma_e2 / T_periods, 0)

# Step 2: lambda = 1 - sqrt( sigma_e^2 / (sigma_e^2 + T * sigma_alpha^2) )
lambda <- 1 - sqrt(sigma_e2 / (sigma_e2 + T_periods * sigma_a2))

# Step 3: quasi-demean by lambda and fit OLS
samp_q <- samp %>%
  group_by(id) %>%
  mutate(
    hgb_q   = hgb - lambda * mean(hgb),
    treat_q = treat - lambda * mean(treat)
  ) %>%
  ungroup()
mod_re <- lm(hgb_q ~ treat_q, data = samp_q)

cat("Estimated sigma_e^2:    ", round(sigma_e2, 3), "\n")
cat("Estimated sigma_alpha^2:", round(sigma_a2, 3), "\n")
cat("Estimated lambda:       ", round(lambda, 3), "\n")
cat("RE estimate of beta_1:  ", round(coef(mod_re)["treat_q"], 3), "\n")

Where does the RE estimate land relative to pooled OLS (\(\approx -0.52\)) and fixed effects (\(\approx +0.84\))? Why there?

Discussion: You should find \(\hat{\lambda} \approx 0.50\) and an RE estimate around \(+0.36\) — squarely between pooled OLS and FE, exactly as the quasi-demeaning algebra predicts (a \(\lambda\)-weighted compromise between no demeaning and full demeaning). And that is the indictment: because the RE assumption \(\text{Cov}(\text{treat}_{it}, \alpha_i) = 0\) fails in this world, RE retains a diluted share of pooled OLS’s bias. Its smaller standard errors are cold comfort — it is more precisely estimating the wrong number. (Reassuringly, this hand-rolled FGLS matches plm::plm(model = "random") on the same sample to three decimals.)

5.2 Choosing between pooled OLS, RE, and FE

Table 3: Decision framework for panel estimators
Estimator Variation used Key assumption on \(\alpha_i\) Handles serial correlation? In our DGP (\(\beta_1 = 0.8\))
Pooled OLS Within + between, unweighted \(\text{Cov}(x_{it}, \alpha_i) = 0\) No \(\approx -0.52\) (badly biased)
Random effects Within + between, \(\lambda\)-weighted \(\text{Cov}(x_{it}, \alpha_i) = 0\) Yes (modeled) \(\approx +0.36\) (biased)
Fixed effects Within only None — \(\alpha_i\) eliminated Partially (cluster SEs still advised) \(\approx +0.84\) (consistent)

The practical logic, informally known as Hausman logic: if RE’s assumption holds, FE and RE both estimate the same \(\beta_1\) consistently, so their estimates should be close (differing only by noise), and you would prefer RE for efficiency. If the estimates are far apart, the RE assumption is the prime suspect — distrust RE and report FE. The formal version is the Hausman test (plm::phtest(mod_fe, mod_re)), which turns the gap into a test statistic; here we will just eyeball the gap against the estimates’ precision.

5.3 Exercise 6: The FE–RE gap (informal Hausman check)

Compute all three estimates on the same sample and the FE\(-\)RE gap. Fill in the gap calculation.

NoteHint 1

The gap is the fixed effects estimate minus the random effects estimate.

NoteHint 2

gap <- b_fe - b_re

TipSolution
set.seed(20231005) ids_samp <- sample(unique(pop$id), 200) samp <- pop %>% filter(id %in% ids_samp) # Pooled OLS b_pooled <- coef(lm(hgb ~ treat, data = samp))["treat"] # Fixed effects (within) samp_dm <- samp %>% group_by(id) %>% mutate(hgb_dm = hgb - mean(hgb), treat_dm = treat - mean(treat)) %>% ungroup() mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm) b_fe <- coef(mod_within)["treat_dm"] # Random effects (FGLS, exactly as in Exercise 5) T_periods <- 5 n_kids <- length(unique(samp$id)) sigma_e2 <- sum(residuals(mod_within)^2) / (nrow(samp) - n_kids - 1) between <- samp %>% group_by(id) %>% summarize(hgb_bar = mean(hgb), treat_bar = mean(treat), .groups = "drop") sigma_B2 <- sum(residuals(lm(hgb_bar ~ treat_bar, data = between))^2) / (n_kids - 2) sigma_a2 <- max(sigma_B2 - sigma_e2 / T_periods, 0) lambda <- 1 - sqrt(sigma_e2 / (sigma_e2 + T_periods * sigma_a2)) samp_q <- samp %>% group_by(id) %>% mutate(hgb_q = hgb - lambda * mean(hgb), treat_q = treat - lambda * mean(treat)) %>% ungroup() b_re <- coef(lm(hgb_q ~ treat_q, data = samp_q))["treat_q"] # The informal Hausman check: how far apart are FE and RE? gap <- b_fe - b_re cat("Pooled OLS estimate: ", round(b_pooled, 3), "\n") cat("Random effects (FGLS):", round(b_re, 3), "\n") cat("Fixed effects (within):", round(b_fe, 3), "\n") cat("True beta_1: 0.800\n\n") cat("FE - RE gap: ", round(gap, 3), "\n") if (abs(gap) > 0.2) { cat("Large gap -> the RE assumption Cov(treat, alpha_i) = 0 looks untenable.\n") cat("Distrust RE; report fixed effects.\n") } else { cat("Small gap -> FE and RE agree; RE's efficiency makes it attractive.\n") }
set.seed(20231005)
ids_samp <- sample(unique(pop$id), 200)
samp <- pop %>% filter(id %in% ids_samp)

# Pooled OLS
b_pooled <- coef(lm(hgb ~ treat, data = samp))["treat"]

# Fixed effects (within)
samp_dm <- samp %>%
  group_by(id) %>%
  mutate(hgb_dm = hgb - mean(hgb), treat_dm = treat - mean(treat)) %>%
  ungroup()
mod_within <- lm(hgb_dm ~ treat_dm, data = samp_dm)
b_fe <- coef(mod_within)["treat_dm"]

# Random effects (FGLS, exactly as in Exercise 5)
T_periods <- 5
n_kids    <- length(unique(samp$id))
sigma_e2  <- sum(residuals(mod_within)^2) / (nrow(samp) - n_kids - 1)
between   <- samp %>%
  group_by(id) %>%
  summarize(hgb_bar = mean(hgb), treat_bar = mean(treat), .groups = "drop")
sigma_B2  <- sum(residuals(lm(hgb_bar ~ treat_bar, data = between))^2) / (n_kids - 2)
sigma_a2  <- max(sigma_B2 - sigma_e2 / T_periods, 0)
lambda    <- 1 - sqrt(sigma_e2 / (sigma_e2 + T_periods * sigma_a2))
samp_q    <- samp %>%
  group_by(id) %>%
  mutate(hgb_q = hgb - lambda * mean(hgb),
         treat_q = treat - lambda * mean(treat)) %>%
  ungroup()
b_re <- coef(lm(hgb_q ~ treat_q, data = samp_q))["treat_q"]

# The informal Hausman check: how far apart are FE and RE?
gap <- b_fe - b_re

cat("Pooled OLS estimate: ", round(b_pooled, 3), "\n")
cat("Random effects (FGLS):", round(b_re, 3), "\n")
cat("Fixed effects (within):", round(b_fe, 3), "\n")
cat("True beta_1:            0.800\n\n")
cat("FE - RE gap:          ", round(gap, 3), "\n")

if (abs(gap) > 0.2) {
  cat("Large gap -> the RE assumption Cov(treat, alpha_i) = 0 looks untenable.\n")
  cat("Distrust RE; report fixed effects.\n")
} else {
  cat("Small gap -> FE and RE agree; RE's efficiency makes it attractive.\n")
}

The gap is about \(0.49\) — enormous relative to standard errors of roughly \(0.08\). If instead supplementation had been assigned by lottery, what would you expect the gap to look like, and which estimator would you then prefer?

Discussion: Under lottery assignment, \(\text{Cov}(\text{treat}_{it}, \alpha_i) = 0\), so FE and RE would both be consistent for \(0.8\) and the gap would be small — attributable to sampling noise. You would then prefer RE: it uses both within and between variation, delivering tighter standard errors, and it can still estimate time-invariant covariates like rural that FE sweeps away. The formal Hausman test makes this comparison rigorous, but the reasoning is exactly what you just did by eye: a big FE–RE gap is the data telling you that the “random” in random effects is wishful thinking.

6 Two-Way Fixed Effects: A Preview

Time to cash in the thought we parked in Section 2. Our DGP includes a secular trend: everyone’s hemoglobin drifts up by \(0.15\) g/dL per year, and the supplementation program expands over time. So within a child, later years are both more-treated and mechanically higher-hemoglobin — child fixed effects do nothing about this, because the trend varies over time, not between children. That is why our FE estimates have hovered slightly above \(0.8\): one-way FE hands the trend’s credit to treatment.

The symmetric fix: add time fixed effects alongside the child fixed effects,

\[ \text{hgb}_{it} = \beta_1 \, \text{treat}_{it} + \alpha_i + \gamma_t + \varepsilon_{it} \]

where \(\gamma_t\) absorbs anything that hits all children in year \(t\) (school meal reforms, a new lab analyzer, a flu season). This is the two-way fixed effects (TWFE) estimator: \(\alpha_i\) removes stable child-level confounding, \(\gamma_t\) removes common shocks, and \(\beta_1\) is identified from variation in treatment within children, relative to what happened to everyone else that year.

6.1 Exercise 7: Adding time fixed effects

Extend the LSDV model with period dummies and compare the estimate to one-way FE.

NoteHint 1

Time fixed effects are dummies for each period, just as child fixed effects are dummies for each id.

NoteHint 2

The formula should be: hgb ~ treat + factor(id) + factor(period)

TipSolution
set.seed(20231005) ids_samp <- sample(unique(pop$id), 200) samp <- pop %>% filter(id %in% ids_samp) mod_oneway <- lm(hgb ~ treat + factor(id), data = samp) mod_twfe <- lm(hgb ~ treat + factor(id) + factor(period), data = samp) cat("One-way FE (child only): ", round(coef(mod_oneway)["treat"], 3), "\n") cat("Two-way FE (child + year):", round(coef(mod_twfe)["treat"], 3), "\n") cat("True beta_1: 0.800\n\n") cat("Estimated year effects (relative to year 1):\n") print(round(coef(mod_twfe)[grep("period", names(coef(mod_twfe)))], 3))
set.seed(20231005)
ids_samp <- sample(unique(pop$id), 200)
samp <- pop %>% filter(id %in% ids_samp)

mod_oneway <- lm(hgb ~ treat + factor(id), data = samp)
mod_twfe   <- lm(hgb ~ treat + factor(id) + factor(period), data = samp)

cat("One-way FE (child only):  ", round(coef(mod_oneway)["treat"], 3), "\n")
cat("Two-way FE (child + year):", round(coef(mod_twfe)["treat"], 3), "\n")
cat("True beta_1:               0.800\n\n")
cat("Estimated year effects (relative to year 1):\n")
print(round(coef(mod_twfe)[grep("period", names(coef(mod_twfe)))], 3))

Which direction did the treatment estimate move when you added time effects, and why? What do the estimated year effects trace out?

Discussion: The estimate falls from about \(0.84\) to about \(0.68\) — adding \(\gamma_t\) strips out the upward secular trend that one-way FE had been crediting to the expanding program (the remaining wobble around \(0.8\) is sampling noise). The year effects climb across periods, tracing the planted \(0.15\)/year drift. TWFE is the workhorse specification behind difference-in-differences: comparing changes within treated units against contemporaneous changes in untreated units is exactly what child-plus-year fixed effects operationalize — a bridge we will cross properly in the causal inference modules.

7 Summary

Table 4: Panel data toolkit at a glance
Concept What it does Watch out for
Pooled OLS Ignores panel structure entirely Biased if \(\text{Cov}(x_{it}, \alpha_i) \neq 0\); serially correlated errors understate SEs
Within/FE Demeans by unit; eliminates \(\alpha_i\) Sweeps away time-invariant covariates; needs within variation in \(x\); time-varying confounding survives
LSDV One dummy per unit; identical slope to within Inflated \(R^2\); computationally heavy with many units
Random effects Quasi-demeans by \(\hat{\lambda}\); efficient Consistent only if \(\text{Cov}(x_{it}, \alpha_i) = 0\)
Hausman logic Big FE–RE gap discredits RE Small gap is necessary, not proof, of RE validity
Two-way FE Adds \(\gamma_t\) for common time shocks Foundation of DiD; complications arise with staggered treatment timing
TipWhat comes next?

Fixed effects gave us our first true identification strategy: a way to get consistent estimates despite unobserved confounding, by being deliberate about which variation we use. The causal inference sequence sharpens this into a research design toolkit — instrumental variables (Module 9), regression discontinuity (Module 10), and difference-in-differences (Module 11), where the two-way fixed effects machinery you just built becomes the estimating equation. The recurring question stays the same: what comparison, exactly, identifies the effect?

References

  • Wooldridge, Jeffrey M. 2020. Introductory Econometrics: A Modern Approach. 7th ed. Cengage Learning. (Chapters 13–14.)
  • Angrist, Joshua D., and Jorn-Steffen Pischke. 2009. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton University Press. (Chapter 5.)
  • Cameron, A. Colin, and Douglas L. Miller. 2015. “A Practitioner’s Guide to Cluster-Robust Inference.” Journal of Human Resources 50 (2): 317–372.
  • Croissant, Yves, and Giovanni Millo. 2008. “Panel Data Econometrics in R: The plm Package.” Journal of Statistical Software 27 (2).
  • Cunningham, Scott. 2021. Causal Inference: The Mixtape. Yale University Press.