Absolutely! You can totally replicate the fixed effects (FE) and random effects (RE) models from Wooldridge's Example 14.4 using the lme4 package. Let's break this down step-by-step to align with what you've already done in plm.

First, Setup: Load Packages and Data #

Assuming you have the wooldridge package (to access the wagepan dataset) installed, start by loading the necessary tools:

# Load required packages
library(lme4)
library(wooldridge)
library(dplyr) # For demeaning data (optional but helpful)

# Load the example dataset
data(wagepan)

1. Random Effects (RE) Model #

In plm, your RE model uses the model = "random" argument. With lme4, we model the individual-specific intercept as a random effect (this aligns with the GLS-based RE estimation from plm).

Basic RE Model (REML, lme4's Default) #

# Random Effects model with random individual intercepts
re_model <- lmer(
  lwage ~ educ + black + hisp + exper + expersq + married + union + (1 | nr),
  data = wagepan
)

# View results
summary(re_model)

Match plm's ML-Based RE Estimation #

By default, plm uses maximum likelihood (ML) for RE models, while lme4 uses restricted maximum likelihood (REML). To align results exactly, add REML = FALSE:

re_model_ml <- lmer(
  lwage ~ educ + black + hisp + exper + expersq + married + union + (1 | nr),
  data = wagepan,
  REML = FALSE
)

summary(re_model_ml)

You’ll see the coefficients match what you got from plm's RE model perfectly.

2. Fixed Effects (FE) Model #

In plm, the FE model (model = "within") removes time-invariant variables (like educ, black, hisp) by demeaning the data. With lme4, we can replicate this in two straightforward ways:

Option 1: Use Individual Fixed Effects as Dummy Variables #

Treat the individual ID (nr) as a fixed factor in a linear model. This directly estimates individual-specific intercepts, just like the within estimator:

# Fixed Effects model with individual dummy variables
fe_model <- lm(
  lwage ~ exper + expersq + married + union + factor(nr),
  data = wagepan
)

summary(fe_model)

Note: The output will show coefficients for each factor(nr) level (individual intercepts), but the key coefficients for time-varying variables (exper, married, etc.) will match your plm FE results.

Option 2: Demean Data First (Mirror plm's Within Estimator) #

If you prefer to replicate the demeaning process that plm uses, use dplyr to demean time-varying variables, then fit an OLS model:

# Demean time-varying variables by individual
wagepan_demeaned <- wagepan %>%
  group_by(nr) %>%
  mutate(
    lwage = lwage - mean(lwage, na.rm = TRUE),
    exper = exper - mean(exper, na.rm = TRUE),
    expersq = expersq - mean(expersq, na.rm = TRUE),
    married = married - mean(married, na.rm = TRUE),
    union = union - mean(union, na.rm = TRUE)
  ) %>%
  ungroup()

# Fit FE model on demeaned data (no intercept needed)
fe_model_demeaned <- lm(
  lwage ~ exper + expersq + married + union - 1,
  data = wagepan_demeaned
)

summary(fe_model_demeaned)

This will produce identical coefficient estimates to plm's within model, since we’re using the same demeaning logic.

Key Notes #

• lme4 is primarily designed for mixed-effects models, so it shines with RE and more complex hierarchical models. For simple FE models, lm (with dummies) or plm might feel more streamlined—but both approaches work!
• Always double-check that your time-invariant variables are excluded from the FE model (they get absorbed by individual fixed effects, just like in plm).

内容的提问来源于stack exchange，提问作者Eric Fail