Fixed and Random Effects Models: A Worked Example

REWB Models with the Cornwell–Rupert PSID Wages Panel

panel data
code
analysis
Author

Robert W. Walker

Published

July 31, 2026

Dataset

The examples use the Cornwell–Rupert PSID Males Wages Panel (Vella & Verbeek 1998), available as Males in R’s plm package: 545 male workers observed annually from 1980 to 1987 (N = 4,360, balanced panel).

Variable Type Description
wage Time-varying Log hourly wage (outcome)
exper Time-varying Years of work experience
union_d Time-varying Union member (1 = yes)
married_d Time-varying Married (1 = yes)
health_d Time-varying Poor health (1 = yes)
school Time-invariant Years of schooling
black_d Time-invariant Black ethnicity (1 = yes)
hisp_d Time-invariant Hispanic ethnicity (1 = yes)

The central research question is: how does work experience relate to log wages? We distinguish the within-person effect (does your wage grow as you gain experience year by year?) from the between-person effect (do workers with more average experience earn more?).

Load and prepare

library(lme4)       # mixed-effects models
library(plm)        # panel FE / RE

data("Males", package = "plm")
d           <- Males
d$id        <- d$nr
d$union_d   <- as.integer(d$union   == "yes")
d$married_d <- as.integer(d$married == "yes")
d$health_d  <- as.integer(d$health  == "yes")
d$black_d   <- as.integer(d$ethn   == "black")
d$hisp_d    <- as.integer(d$ethn   == "hisp")

# Decompose experience into within and between components
d$exper_bar <- ave(d$exper, d$id)      # person-specific mean
d$exper_win <- d$exper - d$exper_bar   # deviation from own mean

pdat <- pdata.frame(d, index = c("id", "year"))
cat(sprintf("Individuals: %d  |  Years: 1980–1987  |  N: %d\n",
            length(unique(d$id)), nrow(d)))
Individuals: 545  |  Years: 1980–1987  |  N: 4360
/* males_wages.csv should be in your working directory.        */
/* Alternatively: webuse nlswork, clear                        */
/* (then rename: ln_wage→wage, ttl_exp→exper, grade→school)   */

import delimited using "males_wages.csv", clear
xtset id year

label variable wage      "Log wage"
label variable exper     "Work experience (years)"
label variable school    "Years of schooling (time-invariant)"

/* Decompose experience into within and between components      */
bysort id: egen exper_bar = mean(exper)
gen exper_win = exper - exper_bar

label variable exper_bar "Person mean of experience (between)"
label variable exper_win "exper minus person mean   (within)"

di "Individuals: 545  |  Years: 1980-1987  |  N: " _N

Confirm within-SD of time-invariant variables

vars <- c("school", "black_d", "hisp_d", "exper", "union_d", "married_d")
wsd <- sapply(vars, function(v)
  mean(tapply(d[[v]], d$id, sd, na.rm = TRUE), na.rm = TRUE))

data.frame(Variable = vars,
           `Within-person SD` = round(wsd, 4),
           Role = c("Time-invariant", "Time-invariant", "Time-invariant",
                    "Time-varying", "Time-varying", "Time-varying"),
           check.names = FALSE) |>
  kable(align = "lcc")
Variable Within-person SD Role
school school 0.0000 Time-invariant
black_d black_d 0.0000 Time-invariant
hisp_d hisp_d 0.0000 Time-invariant
exper exper 2.4495 Time-varying
union_d union_d 0.1954 Time-varying
married_d married_d 0.2576 Time-varying 
foreach v in school black_d hisp_d exper union_d married_d {
    bysort id: egen `v'_sd = sd(`v')
    sum `v'_sd, meanonly
    di "`v' within-SD: " r(mean)
    drop `v'_sd
}

Model Estimation

Model 1 — Pooled OLS

The coefficient on exper blends the within and between effects. Because clustering is ignored, standard errors are underestimated.

m_ols <- lm(wage ~ exper + school + union_d + married_d +
              health_d + black_d + hisp_d, data = d)
summary(m_ols)$coefficients[, 1:2] |>
  round(4) |> kable()
Estimate Std. Error
(Intercept) 0.0255 0.0632
exper 0.0501 0.0029
school 0.1035 0.0046
union_d 0.1831 0.0171
married_d 0.1125 0.0157
health_d -0.0515 0.0566
black_d -0.1428 0.0236
hisp_d 0.0124 0.0208 
regress wage exper school union_d married_d health_d black_d hisp_d

Model 2 — Fixed Effects

Within-group demeaning eliminates all between-person variation. Time-invariant variables (school, black_d, hisp_d) are absorbed by the person dummies and cannot be estimated.

m_fe <- plm(wage ~ exper + union_d + married_d + health_d,
            data   = pdat,
            model  = "within",
            effect = "individual")
summary(m_fe)$coefficients[, 1:2] |> round(4) |> kable()
Estimate Std. Error
exper 0.0599 0.0026
union_d 0.0836 0.0194
married_d 0.0608 0.0183
health_d -0.0183 0.0475 
/* school, black_d, hisp_d cannot be included -- time-invariant */
xtreg wage exper union_d married_d health_d, fe

Model 3 — Standard Random Effects

Assumes the within and between effects of experience are equal (\(\beta_W = \beta_B\)). The Hausman-style test below shows this assumption is rejected.

m_re <- lmer(wage ~ exper + school + union_d + married_d +
               health_d + black_d + hisp_d + (1 | id),
             data = d, REML = FALSE)
round(summary(m_re)$coefficients[, 1:2], 4) |> kable()
Estimate Std. Error
(Intercept) -0.0466 0.1114
exper 0.0580 0.0025
school 0.1081 0.0089
union_d 0.1093 0.0179
married_d 0.0753 0.0168
health_d -0.0234 0.0466
black_d -0.1412 0.0481
hisp_d 0.0159 0.0430 
# plm_re <- plm(wage ~ exper + school + union_d + married_d + health_d + black_d + hisp_d, data = pdat, model="random")
# summary(plm_re) nearly identical
xtreg wage exper school union_d married_d health_d black_d hisp_d, re

Hausman-Style Test

The test checks whether the within and between effects of experience differ — i.e., whether the contextual effect \(\beta_{2C} = 0\). A significant result means the two effects are distinct and standard RE produces an uninterpretable blend.

m_rewb_h <- lmer(wage ~ exper_win + exper_bar + school +
                   union_d + married_d + health_d +
                   black_d + hisp_d + (1 | id),
                 data = d, REML = FALSE)

lrt <- anova(m_re, m_rewb_h)
data.frame(
  Test       = "LRT: Standard RE vs REWB",
  `Chi-sq`   = round(lrt$Chisq[2], 3),
  df         = 1,
  `p-value`  = round(lrt$`Pr(>Chisq)`[2], 4),
  Conclusion = ifelse(lrt$`Pr(>Chisq)`[2] < 0.05,
                      "Reject H₀: within ≠ between → use REWB",
                      "Cannot reject H₀: standard RE acceptable"),
  check.names = FALSE
) |> kable()
Test Chi-sq df p-value Conclusion
LRT: Standard RE vs REWB 5.122 1 0.0236 Reject H₀: within ≠ between → use REWB 
quietly xtreg wage exper union_d married_d health_d, fe
estimates store fe_stored
quietly xtreg wage exper school union_d married_d health_d black_d hisp_d, re
hausman fe_stored ., sigmamore
/* Significant result => within ≠ between => use REWB           */

Model 4 — Within-Between RE (REWB)

The recommended model. exper_win captures the within-person effect of gaining experience year by year; exper_bar captures the between-person effect of having more average experience. Time-invariant variables are estimable — a key advantage over FE.

m_rewb <- lmer(wage ~ exper_win + exper_bar + school +
                 union_d + married_d + health_d +
                 black_d + hisp_d + (1 | id),
               data = d, REML = FALSE)
round(summary(m_rewb)$coefficients[, 1:2], 4) |> kable()
Estimate Std. Error
(Intercept) 0.2725 0.1791
exper_win 0.0591 0.0025
exper_bar 0.0332 0.0112
school 0.0946 0.0107
union_d 0.1104 0.0179
married_d 0.0752 0.0167
health_d -0.0226 0.0466
black_d -0.1340 0.0480
hisp_d 0.0173 0.0428 
/* exper_win and exper_bar computed earlier via:                */
/*   bysort id: egen exper_bar = mean(exper)                    */
/*   gen exper_win = exper - exper_bar                          */
mixed wage exper_win exper_bar school ///
      union_d married_d health_d black_d hisp_d || id:, reml

Model 5 — Mundlak Parameterisation

Algebraically equivalent to REWB. Uses raw exper rather than the demeaned form, and adds exper_bar as a control. The coefficient on exper is the within effect \(\hat{\beta}_W\); the coefficient on exper_bar is the contextual effect \(\hat{\beta}_{2C}\).

The identity \(\hat{\beta}_W + \hat{\beta}_{2C} = \hat{\beta}_{2B}\) (REWB between effect) must hold exactly.

m_mndl <- lmer(wage ~ exper + exper_bar + school +
                 union_d + married_d + health_d +
                 black_d + hisp_d + (1 | id),
               data = d, REML = FALSE)
coefs <- round(summary(m_mndl)$coefficients[, 1:2], 4)
kable(coefs)
Estimate Std. Error
(Intercept) 0.2725 0.1791
exper 0.0591 0.0025
exper_bar -0.0259 0.0114
school 0.0946 0.0107
union_d 0.1104 0.0179
married_d 0.0752 0.0167
health_d -0.0226 0.0466
black_d -0.1340 0.0480
hisp_d 0.0173 0.0428 
bW  <- fixef(m_mndl)["exper"]
b2C <- fixef(m_mndl)["exper_bar"]
b2B <- fixef(m_rewb)["exper_bar"]

data.frame(
  Quantity          = c("β_W  (Mundlak, raw exper)",
                        "β_2C (Mundlak, exper_bar)",
                        "Sum  (β_W + β_2C)",
                        "β_2B (REWB,    exper_bar)"),
  Estimate = round(c(bW, b2C, bW + b2C, b2B), 5),
  check.names = FALSE
) |> kable(caption = "Mundlak identity: β_W + β_2C = β_2B")
Mundlak identity: β_W + β_2C = β_2B
Quantity Estimate
β_W (Mundlak, raw exper) 0.05905
β_2C (Mundlak, exper_bar) -0.02587
Sum (β_W + β_2C) 0.03318
β_2B (REWB, exper_bar) 0.03318 
xtreg wage exper exper_bar school ///
      union_d married_d health_d black_d hisp_d, re

/* Verify identity manually:                                    */
scalar bW  = _b[exper]
scalar b2C = _b[exper_bar]
di "beta_W + beta_2C = " bW + b2C
/* Compare with beta_2B from the REWB mixed model above        */

Model 6 — REWB with Random Slopes

Allows the within-experience effect to vary across individuals. Bell et al. (2019, Section 4.1, Figure 1) show that omitting random slopes when they exist in the true DGP produces anti-conservative standard errors in both FE and RE models.

m_rs <- lmer(wage ~ exper_win + exper_bar + school +
               union_d + married_d + health_d +
               black_d + hisp_d + (1 + exper_win | id),
             data    = d,
             REML    = FALSE,
             control = lmerControl(optimizer = "bobyqa"))

# Fixed effects
round(summary(m_rs)$coefficients[, 1:2], 4) |> kable(caption = "Fixed effects")
Fixed effects
Estimate Std. Error
(Intercept) 0.2481 0.1789
exper_win 0.0590 0.0034
exper_bar 0.0359 0.0112
school 0.0952 0.0107
union_d 0.1101 0.0179
married_d 0.0758 0.0174
health_d -0.0507 0.0451
black_d -0.1306 0.0479
hisp_d 0.0163 0.0428 
# Random effects variance components
as.data.frame(VarCorr(m_rs))[, c("grp","var1","var2","vcov","sdcor")] |>
  transform(vcov  = round(vcov,  5),
            sdcor = round(sdcor, 5)) |>
  kable(caption = "Random effects variance components")
Random effects variance components
grp var1 var2 vcov sdcor
id (Intercept) NA 0.10935 0.33069
id exper_win NA 0.00318 0.05635
id (Intercept) exper_win 0.00136 0.07295
Residual NA NA 0.10593 0.32547 
mixed wage exper_win exper_bar school ///
      union_d married_d health_d black_d hisp_d ///
      || id: exper_win, reml covariance(unstructured)

Results Summary

make_row <- function(model, w, b, sch, blk) {
  data.frame(Model = model,
             `exper within` = round(w,   4),
             `exper between` = if (is.na(b)) NA else round(b, 4),
             school         = if (is.na(sch)) NA else round(sch, 4),
             black_d        = if (is.na(blk)) NA else round(blk, 4),
             check.names    = FALSE)
}

summary_tbl <- rbind(
  make_row("OLS",
           coef(m_ols)["exper"],       NA,
           coef(m_ols)["school"],
           coef(m_ols)["black_d"]),
  make_row("Fixed Effects",
           as.numeric(coef(m_fe)["exper"]), NA, NA, NA),
  make_row("Standard RE",
           fixef(m_re)["exper"],        NA,
           fixef(m_re)["school"],
           fixef(m_re)["black_d"]),
  make_row("REWB",
           fixef(m_rewb)["exper_win"],
           fixef(m_rewb)["exper_bar"],
           fixef(m_rewb)["school"],
           fixef(m_rewb)["black_d"]),
  make_row("Mundlak",
           fixef(m_mndl)["exper"],
           fixef(m_mndl)["exper_bar"],
           fixef(m_mndl)["school"],
           fixef(m_mndl)["black_d"]),
  make_row("REWB + Random Slopes",
           fixef(m_rs)["exper_win"],
           fixef(m_rs)["exper_bar"],
           fixef(m_rs)["school"],
           fixef(m_rs)["black_d"])
)

summary_tbl |>
  kable(
    caption = "Key coefficients across models. FE cannot estimate time-invariant variables.",
    col.names = c("Model", "exper (within/blend)",
                  "exper_bar (between/contextual)",
                  "school", "black_d"),
    align = "lcccc"
  ) |>
  kable_styling(full_width = FALSE) |>
  row_spec(4, bold = TRUE, background = "#f0f8ff") |>
  add_footnote(c(
    "OLS and Std RE: exper coefficient blends within and between effects.",
    "FE: within effect only; time-invariant variables absorbed by person dummies.",
    "REWB (highlighted): separates both; within ≈ FE; between is new information.",
    "Mundlak: equivalent to REWB; exper_bar = contextual (not between) effect.",
    "For Mundlak: exper_bar column shows contextual effect β_2C, not β_2B."
  ), notation = "alphabet")
Key coefficients across models. FE cannot estimate time-invariant variables.
Model exper (within/blend) exper_bar (between/contextual) school black_d
exper OLS 0.0501 NA 0.1035 -0.1428
1 Fixed Effects 0.0599 NA NA NA
exper1 Standard RE 0.0580 NA 0.1081 -0.1412
exper_win REWB 0.0591 0.0332 0.0946 -0.1340
exper2 Mundlak 0.0591 -0.0259 0.0946 -0.1340
exper_win1 REWB + Random Slopes 0.0590 0.0359 0.0952 -0.1306
a OLS and Std RE: exper coefficient blends within and between effects.
b FE: within effect only; time-invariant variables absorbed by person dummies.
c REWB (highlighted): separates both; within ≈ FE; between is new information.
d Mundlak: equivalent to REWB; exper_bar = contextual (not between) effect.
e For Mundlak: exper_bar column shows contextual effect β_2C, not β_2B.

Key findings:

  • Within effect (REWB exper_win ≈ FE): each additional year of own experience raises wages by approximately 5.9%. The small difference from FE (0.0599 vs 0.0591) reflects GLS vs OLS estimation, not a modelling difference.
  • Between effect (REWB exper_bar): workers with one more year of average experience earn 3.3% more. This is materially smaller than the within effect — between-person experience differences are partly confounded with unobserved stable traits (ability, career selection) that also affect wages.
  • OLS blends both to ~5.0% — an uninterpretable weighted average.
  • Education (time-invariant): FE cannot estimate this; REWB recovers ~9.5% per year of schooling.
  • Black ethnicity (time-invariant): FE cannot estimate this; REWB estimates a wage penalty of approximately 13%.

A Cautionary Note: Nonstationarity

WarningAre wages and experience stationary?

Standard FE, RE, and REWB models assume that the variables in the model are stationary (i.e., their statistical properties do not systematically change over time). If wages or key time-varying predictors follow a unit root (stochastic trend), estimates from level-based models can be spurious. This concern is worth investigating before interpreting the results above.

Experience as a deterministic time trend

The most immediate issue is not stochastic at all. In the Cornwell–Rupert data, exper grows by exactly one year per person per year with no exceptions:

d_sorted <- d[order(d$id, d$year), ]
d_sorted$d_exper <- ave(d_sorted$exper, d_sorted$id,
                        FUN = function(x) c(NA, diff(x)))
cat("Year-to-year change in exper (Δexper):\n")
Year-to-year change in exper (Δexper):
table(d_sorted$d_exper, useNA = "always")

   1 <NA> 
3815  545

Because \(\Delta \text{exper}_{it} = 1\) for every observation, exper is a deterministic linear time trend indexed to when each person entered the labour market. It is perfectly collinear with year within each person.

ImportantWhat the within effect of experience actually identifies

In a within-person (FE or REWB) regression of log wages on experience:

\[y_{it} = \beta_0 + \beta_W \cdot \underbrace{(\text{exper}_{it} - \bar{\text{exper}}_i)}_{\text{exper\_win}} + \ldots\]

because \(\text{exper}_{it} = \text{exper}_{i,t-1} + 1\) always, the demeaned predictor exper_win is proportional to \((t - \bar{t})\) — a person-centred time index. \(\hat{\beta}_W\) therefore estimates the average annual wage growth rate, not the causal return to one additional year of a worker’s accumulated human capital.

This matters for interpretation. The within estimate (≈ 5.9% per year) conflates genuine returns to experience with any systematic change in real wages over the 1980–1987 period (e.g., business cycle, labour-market shifts). Including year dummies would be a natural control, but since exper increments by exactly 1 each year, year dummies and exper_win are perfectly collinear within persons and cannot both be included.

The between effect of exper_bar is also partly a life-cycle / cohort effect: workers who are older on average in 1980–1987 have more average experience, but may also differ in cohort-specific ways.

Panel unit root tests on wages

Even setting aside the experience issue, log wages themselves may contain a stochastic trend. We apply two panel unit root tests with opposite null hypotheses:

# Im-Pesaran-Shin (2003): H0 = unit root in all panels
ips <- purtest(pdat[, "wage"], test = "ips", exo = "intercept", lags = 1)

# Hadri (2000) KPSS-type: H0 = stationarity in all panels
hadri <- purtest(pdat[, "wage"], test = "hadri", exo = "intercept")

data.frame(
  Test       = c("Im–Pesaran–Shin (IPS)",
                 "Hadri KPSS"),
  `H₀`       = c("All panels have a unit root",
                  "All panels are stationary"),
  Statistic  = c(round(ips$statistic$statistic, 3),
                 round(hadri$statistic$statistic, 3)),
  `p-value`  = c(format.pval(ips$statistic$p.value,   digits = 3),
                 format.pval(hadri$statistic$p.value,  digits = 3)),
  Decision   = c("Reject H₀ → stationary",
                 "Reject H₀ → unit root"),
  check.names = FALSE
) |> kable(caption = "Conflicting panel unit root test results for log wages")
Conflicting panel unit root test results for log wages
Test H₀ Statistic p-value Decision
Wtbar Im–Pesaran–Shin (IPS) All panels have a unit root -25.378 <2e-16 Reject H₀ → stationary
z Hadri KPSS All panels are stationary 36.010 <2e-16 Reject H₀ → unit root

The tests give contradictory conclusions — a common finding in short panels:

  • IPS rejects the unit root null, suggesting stationarity.
  • Hadri KPSS rejects the stationarity null, suggesting a unit root.

This contradiction is unsurprising given the very short time dimension (T = 8). Both tests have well-known size and power problems in small T settings (Harris & Tzavalis 1999; Hlouskova & Wagner 2006), and with only 8 observations per individual, individual-level ADF tests are uninformative.

# Check whether first-differenced wages are more clearly stationary
ips_fd <- purtest(diff(pdat[, "wage"]), test = "ips",
                  exo = "intercept", lags = 1)
cat(sprintf(
  "IPS on Δwage: statistic = %.3f, p = %.2e\n → %s\n",
  ips_fd$statistic$statistic,
  ips_fd$statistic$p.value,
  if (ips_fd$statistic$p.value < 0.05) "Rejects unit root (Δwage stationary)"
  else "Cannot reject unit root"))
IPS on Δwage: statistic = -63.897, p = 0.00e+00
 → Rejects unit root (Δwage stationary)

First-differenced wages strongly reject the unit root null, consistent with wages being integrated of order one, I(1) — i.e., the level of wages is non-stationary but the change in wages is stationary.

Implications and practical remedies

NoteWhen does nonstationarity matter most?

The concern is most acute for cross-sectional (between) estimates. If wages are I(1) and individual means \(\bar{y}_i\) vary systematically (as they do when T is small and wages have trends), the between-person regression can be spurious. The within estimator is more robust because the within transformation removes a large portion of trending variation, though it does not guarantee consistency under all unit root specifications.

Practical options, in order of preference:

  1. Include year dummies (two-way FE/REWB). Control for any common time trend in wages. However, as noted above, year dummies and exper_win are perfectly collinear within persons in this dataset — so one must be omitted. A reasonable choice is to replace exper with a quadratic in age or a time counter and include year fixed effects.
# Two-way FE: person + year dummies (year_f as factor)
d$year_f <- factor(d$year)

# REWB with year dummies but WITHOUT exper_win (collinear with year_f within person)
m_twoway <- lmer(wage ~ year_f + exper_bar + school +
                   union_d + married_d + health_d +
                   black_d + hisp_d + (1 | id),
                 data = d, REML = FALSE)

# Show year fixed effects and between effect of experience
fe_coefs <- fixef(m_twoway)
cat("Year effects (relative to 1980):\n")
Year effects (relative to 1980):
round(fe_coefs[grepl("year_f", names(fe_coefs))], 4)
year_f1981 year_f1982 year_f1983 year_f1984 year_f1985 year_f1986 year_f1987 
    0.1121     0.1649     0.2071     0.2738     0.3227     0.3820     0.4402 
cat(sprintf("\nBetween effect of experience (exper_bar): %.4f\n",
            fe_coefs["exper_bar"]))

Between effect of experience (exper_bar): 0.0333
/* Two-way FE: absorb person and year effects                   */
/* exper_win omitted -- collinear with year dummies within id   */
quietly tabulate year, gen(yr)   /* year dummies               */

mixed wage yr2-yr8 exper_bar school ///
      union_d married_d health_d black_d hisp_d || id:, reml
/* Or equivalently: xtreg wage yr2-yr8 exper_bar ..., re       */
  1. First-difference the outcome. Regress \(\Delta \text{wage}_{it}\) on time-varying changes. Since \(\Delta \text{exper} = 1\) always, experience is absorbed into the intercept (average wage growth). This is principled but loses the ability to estimate level effects.
# First-difference model
d_fd <- d[order(d$id, d$year), ]
d_fd$d_wage    <- ave(d_fd$wage,    d_fd$id, FUN = function(x) c(NA, diff(x)))
d_fd$d_union   <- ave(d_fd$union_d, d_fd$id, FUN = function(x) c(NA, diff(x)))
d_fd$d_married <- ave(d_fd$married_d, d_fd$id, FUN = function(x) c(NA, diff(x)))
d_fd$d_health  <- ave(d_fd$health_d,  d_fd$id, FUN = function(x) c(NA, diff(x)))
d_fd <- d_fd[!is.na(d_fd$d_wage), ]

# OLS on first differences (exper absorbed into intercept)
m_fd <- lm(d_wage ~ d_union + d_married + d_health, data = d_fd)
cat("First-difference model (Δwage ~ Δcontrols):\n")
First-difference model (Δwage ~ Δcontrols):
round(summary(m_fd)$coefficients[, 1:2], 4) |> kable()
Estimate Std. Error
(Intercept) 0.0648 0.0073
d_union 0.0418 0.0197
d_married 0.0426 0.0229
d_health -0.0507 0.0429 
cat("\nIntercept = average annual wage growth rate: ",
    round(coef(m_fd)["(Intercept)"], 4), "\n")

Intercept = average annual wage growth rate:  0.0648 
cat("(experience absorbed into intercept since Δexper = 1 always)\n")
(experience absorbed into intercept since Δexper = 1 always)
/* First-difference estimator                                   */
/* exper drops out since Δexper = 1 (absorbed by _cons)        */
xtset id year
gen d_wage    = D.wage
gen d_union   = D.union_d
gen d_married = D.married_d
gen d_health  = D.health_d

regress d_wage d_union d_married d_health
/* Intercept = average annual log-wage growth rate              */

  1. Test for cointegration. If wages and predictors are all I(1) but cointegrated, the level-based REWB model remains valid as an error-correction representation. Pedroni (1999) panel cointegration tests are available in R via the cointReg package and in Stata via xtpedroni.

  2. Interpret cautiously. With T = 8, the evidence on stationarity is genuinely ambiguous. The most defensible position is to report the REWB estimates in levels alongside the two-way FE and first-difference specifications as robustness checks, and note that the within-person estimate of experience captures average annual wage growth rather than the structural return to human capital accumulation.

Summary of the nonstationarity issue

Concern Mechanism Practical effect
exper is a deterministic time trend \(\Delta\text{exper} = 1\) always Within effect ≡ average annual wage growth; collinear with year dummies
Wages may be I(1) Stochastic trend in log wages Level-based regressions potentially spurious
Short T (= 8) Low power of unit root tests Cannot reliably distinguish I(0) from I(1)
Conflicting test results IPS rejects unit root; Hadri rejects stationarity Model specification uncertainty

Bottom line: the REWB results are informative and internally consistent, but the within-person estimate of the experience effect should be interpreted as an average annual wage growth rate rather than a pure return to experience. Year fixed effects (in place of exper_win) are a straightforward robustness check when the goal is to separate individual wage growth from common macroeconomic trends.

References

Bell, A., Fairbrother, M. & Jones, K. (2019). Fixed and random effects models: making an informed choice. Quality & Quantity, 53, 1051–1074. https://doi.org/10.1007/s11135-018-0802-x

Bell, A. & Jones, K. (2015). Explaining fixed effects: random effects modelling of time-series cross-sectional and panel data. Political Science Research and Methods, 3(1), 133–153. Replication files: https://doi.org/10.7910/DVN/23415

Cornwell, C. & Rupert, P. (1988). Efficient estimation with panel data: an empirical comparison of instrumental variables estimators. Journal of Applied Econometrics, 3, 149–155.

Hadri, K. (2000). Testing for stationarity in heterogeneous panel data. Econometrics Journal, 3(2), 148–161.

Im, K. S., Pesaran, M. H. & Shin, Y. (2003). Testing for unit roots in heterogeneous panels. Journal of Econometrics, 115(1), 53–74.

Jordan, S. & Philips, A. Q. (2023). Improving the interpretation of random effects regression results. Political Studies Review, 21(1), 210–220. Code: https://github.com/andyphilips/qdmean

Mundlak, Y. (1978). Pooling of time-series and cross-section data. Econometrica, 46(1), 69–85.

Vella, F. & Verbeek, M. (1998). Whose wages do unions raise? A dynamic model of unionism and wage rate determination for young men. Journal of Applied Econometrics, 13(2), 163–183.