What is Panel Data?

Panel data importantly tracks the same groups over-time, which is what distinguishes it from cross-sectional data alone. It can take the form of repeated measurements on individuals (where each individual serves as it’s own group, \(i\)) through time, measures of groups of individuals (e.g. prevalence of a particular disease, see Figure 1) through time, etc. The data can be balanced or unbalanced similar to other data types wherein missingness can occur in one or both dimensions of the data, and a common sub-type of panel data is large numbers of \(N\) groups, but short timeseries \(T\) across all groups. For example, a large number of companies might be observed for performance metrics quarterly for only a few years.

flowchart TB
    Title["<b>PANEL DATA STRUCTURE</b><br/>Multiple entities observed over multiple time periods"]
    
    EntitiesTitle["<b>Cross-Sectional Dimension</b><br/>e.g. Prevalence of a disease across counties"]
    TimeTitle["<b>Time Dimension</b><br/>e.g. Through years"]
    
    subgraph Entities[" "]
        direction LR
        E1["<b>Prevalence in county 1</b>"]
        E2["<b>Prevalence in county 2</b>"]
        E3["<b>Prevalence in county 3</b>"]
        EDots["<b>⋮</b>"]
        EN["<b>Prevalence in county N</b>"]
    end
    
    subgraph Time[" "]
        direction LR
        T1["<b>2020</b>"]
        T2["<b>2021</b>"]
        T3["<b>2022</b>"]
        TDots["<b>⋮</b>"]
        TT["<b>2025</b>"]
    end
    
    DataTitle["<b>Data Matrix: yᵢₜ</b>"]
    
    subgraph Data[" "]
        direction LR
        R1["y₁₁, y₁₂, y₁₃, ..., y₁ₜ"]
        R2["y₂₁, y₂₂, y₂₃, ..., y₂ₜ"]
        R3["y₃₁, y₃₂, y₃₃, ..., y₃ₜ"]
        RDots["⋮   ⋮   ⋮   ⋱   ⋮"]
        RN["yₙ₁, yₙ₂, yₙ₃, ..., yₙₜ"]
    end
    
    Title --> EntitiesTitle
    Title --> TimeTitle
    EntitiesTitle --> Entities
    TimeTitle --> Time
    Entities --> DataTitle
    Time --> DataTitle
    DataTitle --> Data
    
    Summary["Each cell yᵢₜ represents the observation for county i at time t<br/>Total observations = N counties × T time periods"]
    
    Data --> Summary
    
    style Title fill:#2563eb,stroke:#1e40af,stroke-width:3px,color:#fff
    style EntitiesTitle fill:#60a5fa,stroke:#3b82f6,stroke-width:2px,color:#fff
    style TimeTitle fill:#60a5fa,stroke:#3b82f6,stroke-width:2px,color:#fff
    style DataTitle fill:#f59e0b,stroke:#d97706,stroke-width:2px,color:#fff
    style Entities fill:#dbeafe,stroke:#3b82f6,stroke-width:2px
    style Time fill:#dbeafe,stroke:#3b82f6,stroke-width:2px
    style Data fill:#fef3c7,stroke:#f59e0b,stroke-width:2px
    style Summary fill:#d1fae5,stroke:#10b981,stroke-width:2px

Regardless of data quality, variables in panel data can be related at more than one level in hierarchy, with different relationships at those different levels. Suffice it to say that an important component of panel data is that there are often effects that occur at level 1, which we can term “within-group” effects². These are different than the effects that might be observed at level 2 which might be “contextual” or which we can term “between-group” effects³. A visualization of this is provided in Figure 2.

Show the code

library(ggplot2)
library(dplyr)

set.seed(42)

n_groups <- 5
n_timepoints <- 20
time_points <- 1:n_timepoints

group_slopes <- c(0.8, 0.8, 0.8, 0.8, 0.8)

data_list <- list()

for (group_id in 1:n_groups) {
    group_mean <- 10 + (group_id * 3)
    group_slope <- group_slopes[group_id]

    for (t in time_points) {
        within_effect <- group_slope * t
        noise <- rnorm(1, mean = 0, sd = 0.8)
        outcome_value <- group_mean + within_effect + noise

        data_list[[length(data_list) + 1]] <- data.frame(
            Group = paste0("Group ", group_id),
            Time = t,
            Outcome = outcome_value
        )
    }
}

df <- do.call(rbind, data_list)

overall_means <- df %>%
    group_by(Time) %>%
    summarise(Mean_Outcome = mean(Outcome))

plot <- ggplot() +
    geom_smooth(
        data = df, aes(x = Time, y = Outcome, colour = Group, group = Group),
        linewidth = 1.2, alpha = 0.1, method = "lm"
    ) +
    geom_point(
        data = df, aes(x = Time, y = Outcome, colour = Group),
        size = 3, alpha = 0.8
    ) +
    geom_smooth(
        data = overall_means, aes(x = Time, y = Mean_Outcome),
        linetype = "dashed", linewidth = 1.5, colour = "black", method = "lm"
    ) +
    geom_point(
        data = overall_means, aes(x = Time, y = Mean_Outcome),
        size = 3, colour = "black", alpha = 0.7
    ) +
    scale_x_continuous(breaks = time_points) +
    scale_colour_brewer(palette = "Set2") +
    labs(
        title = "Within-Group vs. Between-Group Effects in Panel Data",
        x = "Time Point",
        y = "Outcome Variable",
        colour = "Group"
    ) +
    theme_bw()
plot

In this guide we’ll go over the formalism behind extracting these various kinds of differences from panel data, including an overview of different methods and approaches available.

Theoretical Intuition

As we’ve established by now, panel data (or whatever we wish to call it, “longitudinal” or “repeated measures”) consists of repeated observations on the same groups over time, denoted as \(y_{it}\) for group \(i = 1, ..., N\) at time \(t = 1, ..., T\). The attentive reader may notice that this leads to an obvious challenge in the form of non-independence⁴. This problematic correlation structure requires special consideration as it creates both inferential challenges (affecting the precision of our estimates) and identification challenges (affecting the consistency and interpretation of our estimates).

To get a better sense of what this means under the hood we can formalize this notion and look to what the formalism guides us to by way of solution. Let’s say that the observed outcome in a panel data setting can be decomposed as:

\[\begin{equation} y_{it} = \boldsymbol{\beta}x^T_{it} + \varepsilon_{it} \end{equation}\] where \(y_{it}\) is the observed value of some outcome variable in group \(i\) at time \(t\). We also assume that there is a set of \(K\) explanatory variables \(x_{it}\) which are observed as well, and therefore a set of \(K\) coefficients \(\mathbf{\beta}\), which are estimated. There is then an error term \(\varepsilon_{it}\), which is assumed to have mean 0, and contains multiple components that reflect the hierarchical nature of the data:

\[\varepsilon_{it} = \alpha_i + u_{it}\]

The component \(\alpha_i\) represents the group-specific effect, which is time-invariant and captures all systematic, but unobserved characteristics of group \(i\) that remain constant over the observation period. The variance of that effect, \(Var(\alpha_i) = \tau^2\) quantifies the amount of variation in the outcome that is attributable to differences between groups. Additionally, \(u_{it}\) is a random variable that is independent and identically distributed (i.i.d.) with mean 0 and some variance. It represents the idiosyncratic error that varies across both groups and time periods, capturing random variation specific to group \(i\) at time \(t\).

This will become important later, but typically since we wish to characterize the within-group and between-group variation, we can account for variation within groups by any set of time-dependent covariates \(x_{it}\) and our error term \(\varepsilon_{it}\). On the other hand, variation between groups is accounted for by the covariates but also the group-specific effects, \(\alpha_i\).

The challenge in using standard regression approaches arises due to the that fact that for any two observations from the same group \(i\) at different time periods \(t\) and \(s\), the error terms are correlated:

\[\text{Cov}(\varepsilon_{it}, \varepsilon_{is}) = \text{Var}(\alpha_i) \neq 0 \quad \text{for } t \neq s\]

This correlation exists because both observations share the same group-specific effect \(\alpha_i\). This is known as the intraclass correlation or within-cluster correlation, which can be expressed as:

\[\begin{equation} \rho = \frac{\text{Var}(\alpha_i)}{\text{Var}(\alpha_i) + \text{Var}(u_{it})} = \frac{\tau^2_{\alpha}}{\tau^2_{\alpha} + \sigma^2_u} \end{equation}\]

This quantity represents the proportion of total variance that is attributable to differences between groups rather than within-group variation over time. The presence of positive intraclass correlation has important consequences for inference. Standard ordinary least squares (OLS) methods assume independent observations and will understate the true standard errors when within-group correlation is present, leading to overconfident inference and inflated Type I error rates.

A second and often more serious problem is related to the heterogeneity bias⁵ of the system. This challenge arises when the group-specific effect \(\alpha_i\) is correlated with the observed predictors. If \(\text{Cov}(\mathbf{x}_{it}, \alpha_i) \neq 0\), then the group-specific effect represents an omitted variable that is correlated with the included predictors, causing omitted variable bias. In this case, the result of this troublesome correlation is that the OLS estimator will be inconsistent:

\[\mathbb{E}[\hat{\beta}_{OLS}] = \boldsymbol{\beta} + \text{bias}\]

The bias arises because standard OLS incorrectly attributes variation in \(y\) that is actually caused by \(\alpha_i\) to the observed predictors \(\mathbf{x}_{it}\), since these are correlated with \(\alpha_i\). This endogeneity problem is distinct from the correlation problem discussed above, because even if we could somehow compute the correct standard errors, our point estimates would still be biased and inconsistent in the presence of endogeneity.

Now that we have a sense of what makes panel data different and why we can’t just use standard OLS to answer our questions, what do we actually do?? In the rest of this Methods section, we outline some of the options that exist to effectively do regression on panel data.

Option 1: Fixed-effects models

Fixed-effects models have been the “gold-standard” to deal with panel data for some time in the Economics literature, but in the interest of not burrying the lead, we do not recommend the use of fixed-effects models (FE models) as the default for panel data. However, due to their commonality, we outline here their formulation and use, along with important caveats.

The use of a FE model stems from the challenge we identified earlier as heterogeneity bias. If a group-specific effect \(\alpha_i\) is correlated with any explanatory variables, we cannot treat it as a random variable in the classic sense (but see Section 2.3). Therefore, we simply do not attempt to measure the variation between groups (for which we would need to use a group-specific effect, \(\alpha_i\)). The fixed effects model, using the notation of Bell and Jones (2015) is then defined as

\[\begin{equation} y_{it} = \sum_{j=1}^N \beta_{0j} D_{ji} + \mathbf{\beta_1}\mathbf{x_{it}} + u_{it} \end{equation}\]

where \(D_i\) is a dummy variable that equals one for all observations from group \(i\) and zero otherwise. In this formulation, each group receives its own intercept \(\beta_{0i}\), which captures that group’s specific effect \(\alpha_i\).

In this formulation, we’ve explicitly dealt with our problem by controlling for all time-invariant group-specific characteristics. By including these group-specific intercepts, the model accounts for all time-invariant characteristics of each group, whether observed or unobserved, eliminating the problematic correlation \(\text{Cov}(\mathbf{x}_{it}, \alpha_i) \neq 0\).

To see how this works, recall that our concern was the correlation between the group-specific effect \(\alpha_i\) and our covariates. In the fixed effects specification above, \(\alpha_i\) is no longer part of the error term—it has been absorbed into the estimated parameters \(\beta_{0i}\). The residual \(u_{it}\) now contains only the idiosyncratic, time-varying component of the error, which by assumption is uncorrelated with the covariates: \(\text{Cov}(\mathbf{x}_{it}, u_{it}) = 0\). This eliminates heterogeneity bias because we have explicitly modeled the source of the endogeneity.

Option 2: Random-effects models

The rest of this methods section will focus on (our) preferred method which is the various form that a random effects model can take. Random effects models take a fundamentally different approach to panel data than fixed effects models. Rather than controlling out all group-specific variation, RE models explicitly partition the total variance into components at different levels of the hierarchy. However, the standard formulation of RE models is susceptible to heterogeneity bias, which has led many researchers to abandon them in favor of FE models. As we will show, this abandonment is premature—a correctly specified RE model can eliminate heterogeneity bias while retaining substantial advantages over FE approaches.

R Implementation Options and Setup

Before diving in, let’s briefly survey the available tools for panel data analysis in R. The R package ecosystem offers several excellent packages, each with particular strengths.

For fixed effects models, the most popular choice is the plm package (Croissant and Millo 2008), which was designed specifically for panel data and offers a wide range of estimation methods, diagnostic tests, and convenience functions. The package handles both balanced and unbalanced panels efficiently and provides straightforward syntax for panel data structures. An alternative is the fixest package, which is particularly fast for large datasets and offers excellent support for high-dimensional fixed effects. For those already familiar with mixed-effects models, the lme4 package (Bates et al. 2014) can also fit fixed effects models by including group indicators as factors in the model formula, though this approach is less efficient for large numbers of groups.

For random effects models, including the REWB specification we advocate, the lme4 package is the standard tool in R. It provides robust and efficient estimation of mixed-effects models with flexible random effects structures. The package’s lmer() function handles the random effects naturally and provides excellent tools for inference. The nlme package is an older alternative that offers additional features for modeling complex correlation structures, though lme4 is generally preferred for standard panel data applications. Some researchers also use the plm package for random effects estimation, particularly when they want to keep all panel data analysis within a single package framework.

Let’s load the packages we’ll use for our examples:

# Load required packages
library(plm) # For fixed effects models
library(lme4) # For random effects REWB models
library(dplyr) # For data manipulation
library(ggplot2) # For visualization
library(gridExtra) # For arranging plots
library(magrittr) # for the pipes

# Set seed for reproducibility
set.seed(42)

Simulated School Data

To most clearly demonstrate what each modeling approach can and cannot tell us, we use a single dataset throughout this section. We simulate data for forty schools observed over twelve years, where schools differ in both their funding levels and in whether they are public or private institutions.

In our simulation, we create several important features. Changes in funding within schools over time have a modest effect on test scores (the within effect, which we set to 0.6 points per thousand dollars). However, schools with higher average funding across all years show substantially higher test scores (the between effect, which we set to 1.8 points per thousand dollars). This difference might arise because historically well-funded schools have accumulated advantages beyond just current funding—better facilities, experienced teachers, engaged communities. Additionally, private schools score higher on average than public schools (we set this effect to 10 points), even after accounting for funding differences. This represents unmeasured differences in resources, selectivity, or educational approaches between school types.

panel_data <- readr::read_csv(
    paste0(
        here::here("./research-guides/panel-regression/assets"),
        "/data/school_panel.csv"
    ),
    show_col_types = FALSE
) |>
    dplyr::mutate(
        school_id = factor(school_id),
        school_type = factor(school_type)
    )

panel_pdata <- plm::pdata.frame(panel_data, index = c("school_id", "year"))

Let’s visualize our data to understand its structure. We’ll create two plots: one showing the overall relationship between funding and test scores, and another highlighting the differences between public and private schools:

library(gridExtra)

# Calculate school-level summaries for colouring
school_summaries <- panel_data %>%
    dplyr::group_by(school_id) %>%
    dplyr::summarise(mean_funding = mean(funding))

panel_data <- panel_data %>%
    dplyr::left_join(school_summaries, by = "school_id")

# Create two plots
p1 <- ggplot(panel_data, aes(x = funding, y = test_score, group = school_id)) +
    geom_line(aes(colour = mean_funding), alpha = 0.4) +
    geom_point(aes(colour = mean_funding), alpha = 0.3, size = 1) +
    scale_colour_gradient(
        low = "blue", high = "red",
        name = "School's\nMean\nFunding"
    ) +
    labs(
        title = "Test Scores vs. Funding",
        subtitle = "Each line represents one school over time",
        x = "Per-Student Funding ($1000s)",
        y = "Test Score"
    ) +
    theme_bw()

p2 <- ggplot(panel_data, aes(x = funding, y = test_score, group = school_id)) +
    geom_line(aes(colour = school_type), alpha = 0.4) +
    geom_point(aes(colour = school_type), alpha = 0.3, size = 1) +
    scale_colour_manual(
        values = c("Public" = "steelblue", "Private" = "darkred"),
        name = "School\nType"
    ) +
    labs(
        title = "Differences by School Type",
        subtitle = "Private schools tend to have higher funding and scores",
        x = "Per-Student Funding ($1000s)",
        y = "Test Score"
    ) +
    theme_bw()

grid.arrange(p1, p2, ncol = 2)

OK so what are we actually looking at here?? This plot shows us the multi-level structure of our data. Within individual schools (along each line), test scores increase modestly with funding. Comparing across schools, those with higher average funding show substantially higher average scores. Additionally, private schools generally occupy the upper range of both funding and achievement. One can easily imagine how these simulated data could represent data found in the real world. Our modeling goal is to disentangle these relationships and understand what each level tells us about the role of funding and school type in student achievement.

Fixed Effects Model

We begin with the fixed effects approach, which we would like to highlight is useful for it’s ability to control for all time-invariant confounding. Let’s fit the model!

# Fit fixed effects model
fe_model <- plm(
    test_score ~ funding,
    data = panel_pdata,
    model = "within"
)

# Display results
summary(fe_model)

Oneway (individual) effect Within Model

Call:
plm(formula = test_score ~ funding, data = panel_pdata, model = "within")

Balanced Panel: n = 40, T = 12, N = 480

Residuals:
    Min.  1st Qu.   Median  3rd Qu.     Max. 
-20.3846  -4.1995   0.0621   4.5808  18.6010 

Coefficients:
        Estimate Std. Error t-value  Pr(>|t|)    
funding 0.545973   0.058296  9.3655 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    27751
Residual Sum of Squares: 23130
R-Squared:      0.16653
Adj. R-Squared: 0.090586
F-statistic: 87.7129 on 1 and 439 DF, p-value: < 2.22e-16

# Extract the coefficient and standard error
fe_coef <- coef(fe_model)["funding"]
fe_se <- sqrt(diag(vcov(fe_model)))["funding"]

# Construct 95% confidence interval
ci_lower <- fe_coef - 1.96 * fe_se
ci_upper <- fe_coef + 1.96 * fe_se

cat(sprintf(
    "Within-school effect of funding: %.3f (95%% CI: [%.3f, %.3f])\n",
    fe_coef, ci_lower, ci_upper
))

Within-school effect of funding: 0.546 (95% CI: [0.432, 0.660])

Visualizing FE Estimates

Let’s further understand what the fixed effects model has done with a couple plots.

# Add fitted values to data
panel_data$fitted_fe <- predict(fe_model)

ggplot(panel_data, aes(x = funding, y = test_score)) +
    geom_point(aes(colour = school_type), alpha = 0.3, size = 1.5) +
    geom_line(aes(y = fitted_fe, group = school_id, colour = school_type),
        linewidth = 0.6, alpha = 0.6
    ) +
    scale_colour_manual(
        values = c("Public" = "steelblue", "Private" = "darkred"),
        name = "School Type"
    ) +
    labs(
        title = "Fixed Effects Model: School-Specific Predictions",
        subtitle = "Each line shows one school's predicted trajectory (all with same slope)",
        x = "Per-Student Funding ($1000s)",
        y = "Test Score"
    ) +
    theme_bw()

What we can see here perhaps most obviously are the parallel nature of the fixed effects model’s predictions. Every school gets its own baseline (intercept), but they all share the same slope—the within effect of funding.

REWB Model

Now we turn to the Random Effects Within-Between model, which can in theory address all three limitations of the fixed effects approach while not succumbing to heterogeneity bias. Let’s spend some time building this model up.

# Create within and between components of funding
library(magrittr)
library(dplyr)

panel_data_wb <- panel_data %>%
    dplyr::group_by(school_id) %>%
    dplyr::mutate(
        funding_mean = mean(funding), # Between component
        funding_centered = funding - funding_mean # Within component
    ) %>%
    dplyr::ungroup()
head(panel_data_wb)

# A tibble: 6 × 9
  school_id  year school_type funding test_score mean_funding fitted_fe
  <fct>     <dbl> <fct>         <dbl>      <dbl>        <dbl> <pseries>
1 1             1 Public         45.3       576.         47.6 574.8135 
2 1             2 Public         51.1       581.         47.6 577.9822 
3 1             3 Public         48.1       571.         47.6 576.3325 
4 1             4 Public         47.5       583.         47.6 576.0022 
5 1             5 Public         43.8       570.         47.6 574.0233 
6 1             6 Public         49.3       579.         47.6 576.9932 
# ℹ 2 more variables: funding_mean <dbl>, funding_centered <dbl>

# Fit basic REWB model (without school type)
rewb_model_basic <- lmer(
    test_score ~ funding_centered + funding_mean + (1 | school_id),
    data = panel_data_wb,
    REML = TRUE
)

# Display results
summary(rewb_model_basic)

Linear mixed model fit by REML ['lmerMod']
Formula: test_score ~ funding_centered + funding_mean + (1 | school_id)
   Data: panel_data_wb

REML criterion at convergence: 3386.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.74480 -0.59808 -0.00035  0.60619  2.65033 

Random effects:
 Groups    Name        Variance Std.Dev.
 school_id (Intercept) 86.30    9.290   
 Residual              52.69    7.259   
Number of obs: 480, groups:  school_id, 40

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      499.6550     6.3670  78.476
funding_centered   0.5460     0.0583   9.366
funding_mean       1.9033     0.1204  15.811

Correlation of Fixed Effects:
            (Intr) fndng_c
fndng_cntrd  0.000        
funding_men -0.972  0.000 

test_score ~ funding_centered + funding_mean + (1 | school_id)

# Extract coefficients
within_coef <- fixef(rewb_model_basic)["funding_centered"]
between_coef <- fixef(rewb_model_basic)["funding_mean"]

# Extract standard errors
se <- sqrt(diag(vcov(rewb_model_basic)))
within_se <- se["funding_centered"]
between_se <- se["funding_mean"]

# Create a results data frame
rewb_results <- data.frame(
    Effect = c("Within-school", "Between-school"),
    Estimate = c(within_coef, between_coef),
    SE = c(within_se, between_se),
    Interpretation = c(
        "When a school's funding increases by $1000, test scores increase by this amount",
        "Schools with $1000 higher average funding score this much higher on average"
    )
)

# Print formatted table
knitr::kable(rewb_results,
    digits = 3,
    col.names = c("Effect Type", "Estimate", "Std. Error", "Interpretation")
)

# Test whether within and between effects differ
diff <- between_coef - within_coef
diff_se <- sqrt(within_se^2 + between_se^2)
z_stat <- diff / diff_se
p_value <- 2 *
    (1 - pnorm(abs(z_stat)))
cat(p_value)

0

Visualizing Within and Between Effects

Let’s visualize what the REWB model has estimated. We’ll show both the within-school relationships (individual school trajectories) and the between-school relationship (how average funding relates to average test scores):

# Generate fitted values for within-school relationships
panel_data$fitted_rewb_basic <- fitted(rewb_model_basic)

# For between-school relationship, create predictions at school means
school_means <- panel_data_wb %>%
    group_by(school_id) %>%
    summarise(
        funding_mean = first(funding_mean),
        test_score_mean = mean(test_score),
        school_type = first(school_type)
    )

# Generate population-level predictions for between effect
funding_range <- seq(min(school_means$funding_mean),
    max(school_means$funding_mean),
    length.out = 100
)
between_pred <- predict(
    rewb_model_basic,
    newdata = data.frame(
        funding_centered = 0,
        funding_mean = funding_range,
        school_id = NA
    ),
    re.form = NA
)

between_line <- data.frame(
    funding_mean = funding_range,
    test_score = between_pred
)

# Create visualization
ggplot() +
    # Within-school fitted lines
    geom_line(
        data = panel_data,
        aes(x = funding, y = fitted_rewb_basic, group = school_id),
        alpha = 0.3, colour = "gray50", linewidth = 0.5
    ) +
    # Between-school relationship (thick red line)
    geom_line(
        data = between_line,
        aes(x = funding_mean, y = test_score),
        colour = "red", linewidth = 2
    ) +
    # School means as points
    geom_point(
        data = school_means,
        aes(x = funding_mean, y = test_score_mean),
        size = 3, alpha = 0.6
    ) +
    labs(
        title = "REWB Model: Within and Between Effects",
        subtitle = "Gray lines: within-school effects | Red line: between-school effect",
        x = "Funding ($1000s per student)",
        y = "Test Score"
    ) +
    theme_bw()

This visualization shows (quite elegantly) what we’ve gained by using REWB instead of FE. The gray lines show the within effect that fixed effects also estimates—modest changes within schools as funding varies. But now we also see the red line showing the much stronger between-school relationship: schools with higher average funding score substantially higher overall.

# Fit full REWB model including school type
rewb_model_full <- lmer(
    test_score ~ funding_centered + funding_mean + school_type + (1 | school_id),
    data = panel_data_wb,
    REML = TRUE
)

# Display results
summary(rewb_model_full)

Linear mixed model fit by REML ['lmerMod']
Formula: test_score ~ funding_centered + funding_mean + school_type +  
    (1 | school_id)
   Data: panel_data_wb

REML criterion at convergence: 3372.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.70280 -0.59850  0.00283  0.60940  2.66352 

Random effects:
 Groups    Name        Variance Std.Dev.
 school_id (Intercept) 66.01    8.125   
 Residual              52.69    7.259   
Number of obs: 480, groups:  school_id, 40

Fixed effects:
                  Estimate Std. Error t value
(Intercept)       513.5881     6.9073  74.355
funding_centered    0.5460     0.0583   9.366
funding_mean        1.7353     0.1167  14.877
school_typePublic -10.1020     2.9219  -3.457

Correlation of Fixed Effects:
            (Intr) fndng_c fndng_m
fndng_cntrd  0.000                
funding_men -0.960  0.000         
schl_typPbl -0.583  0.000   0.416 

test_score ~ funding_centered + funding_mean + school_type + (1 | school_id)

# Create a summary table of all effects
coef_summary <- data.frame(
    Effect = c("Within (Funding)", "Between (Funding)", "School Type (Private vs Public)"),
    Estimate = c(
        fixef(rewb_model_full)["funding_centered"],
        fixef(rewb_model_full)["funding_mean"],
        fixef(rewb_model_full)["school_typePublic"]
    ),
    `Std. Error` = sqrt(diag(vcov(rewb_model_full)))[c("funding_centered", "funding_mean", "school_typePrivate")],
    `True Value` = c(0.600, 1.800, 10.00),
    check.names = FALSE
)

knitr::kable(coef_summary,
    digits = 3,
    caption = "Full REWB Model: All Fixed Effects"
)

# Compare funding effects across models
funding_comparison <- data.frame(
    Effect = c("Within Effect", "Between Effect"),
    `Basic Model` = c(
        fixef(rewb_model_basic)["funding_centered"],
        fixef(rewb_model_basic)["funding_mean"]
    ),
    `With School Type` = c(
        fixef(rewb_model_full)["funding_centered"],
        fixef(rewb_model_full)["funding_mean"]
    ),
    Change = c(
        fixef(rewb_model_full)["funding_centered"] - fixef(rewb_model_basic)["funding_centered"],
        fixef(rewb_model_full)["funding_mean"] - fixef(rewb_model_basic)["funding_mean"]
    ),
    check.names = FALSE
)

knitr::kable(funding_comparison,
    digits = 3,
    caption = "Stability of Funding Effects After Adding School Type"
)

Visualizing the Complete Model

Finally, let’s create a visualization that shows all three effects estimated by our full REWB model:

# Generate fitted values
panel_data_wb$fitted_rewb_full <- fitted(rewb_model_full)

# Create separate between lines for public and private schools
between_lines <- expand.grid(
    funding_mean = seq(min(school_means$funding_mean),
        max(school_means$funding_mean),
        length.out = 100
    ),
    school_type = factor(c("Public", "Private"))
)

between_lines$test_score <- predict(
    rewb_model_full,
    newdata = data.frame(
        funding_centered = 0,
        funding_mean = between_lines$funding_mean,
        school_type = between_lines$school_type,
        school_id = NA
    ),
    re.form = NA
)

# Create visualization
ggplot() +
    geom_point(
        data = panel_data_wb,
        aes(x = funding, y = test_score, colour = school_type),
        alpha = 0.3
    ) +
    # Within-school fitted lines (coloured by type)
    geom_line(
        data = panel_data_wb,
        aes(x = funding, y = fitted_rewb_full, group = school_id, colour = school_type),
        alpha = 0.3, linewidth = 0.5
    ) +
    # Between-school relationships (thick lines by type)
    geom_line(
        data = between_lines,
        aes(x = funding_mean, y = test_score, colour = school_type),
        linewidth = 2
    ) +
    scale_colour_manual(
        values = c("Public" = "steelblue", "Private" = "darkred"),
        name = "School Type"
    ) +
    labs(
        title = "Complete REWB Model: Within, Between, and School Type Effects",
        subtitle = "Thin lines: within-school effects | Thick lines: between-school effects by type",
        x = "Funding ($1000s per student)",
        y = "Test Score"
    ) +
    theme_bw() +
    theme(legend.position = "right")

This final visualization captures the full richness of the model as far as we’ve specified it here. Here we can see that the within-group effects (the thin lines for each school are in fact different), and that the between-group effects (thick lines for each school type) show how average funding relates to average achievement, nested by school type. That’s a lot of information. We’ve successfully moved through most of the components of a panel regression now, utilizing all of the relevant parts of the data on hand. However, this is certainly far from all-encompassing, and there are more details and FAQs found below.

How many time periods do I need?

The minimum is two time periods—otherwise you don’t have panel data! However, more periods are strongly recommended, particularly for within-group effects. With only two periods, within-group estimates are based on a single difference for each unit, making them sensitive to measurement error and transitory shocks. With 4-5 or more periods, you have multiple within-unit observations that provide much more stable estimates. However, this is still the rough “minimum”.

What if my panel is unbalanced?

These models handle unbalanced panels naturally—units can have different numbers of observations without requiring any special treatment. The estimation procedures automatically account for this. However, you should investigate why your panel is unbalanced. If observations are missing at random (for example, some surveys were randomly not conducted in certain years), this typically poses no problems. If missingness is related to the outcome or predictors (for example, low-performing schools systematically drop out of the study), you may have selection bias that neither FE nor REWB can address without additional modeling assumptions.

Can I use these methods with binary or count outcomes?

Yes, but with some important caveats. For fixed effects with non-continuous outcomes, you need specialized estimators: conditional logit for binary outcomes, conditional Poisson for counts. These are available in most statistical packages but have some limitations compared to the linear case. For REWB models, you use generalized linear mixed models (GLMMs): logistic regression with random effects for binary outcomes, Poisson or negative binomial with random effects for counts. As discussed in the Bell and Jones (2015) paper we referenced earlier, the REWB decomposition in GLMMs doesn’t perfectly eliminate all between-unit confounding (unlike in the linear case), but simulations suggest the bias is typically modest. The advantage of REWB—allowing time-invariant predictors and variance components—often outweighs this small potential for bias.

What about three-level data structures?

Not a big deal! For example, with students nested in classrooms nested in schools, you might have time-varying predictors at the student level, which can be decomposed into within-student, within-classroom (between-student), and between-classroom components. The lme4 package in R handles such structures with formulas we’ve not covered here, but there is ample documentation for, like (1 | school_id/classroom_id) for nested random effects.

What if my treatment variable barely varies within units?

This is a serious challenge for these models, and ultimately is a question for a power analysis. If you’re trying to measure a very small effect, you will need sufficiently large amounts of data to find it.

Can’t I just use pooled OLS and ignore the panel structure?

You could, but you almost certainly shouldn’t. Pooled OLS treats all observations as independent, which means:

1. Biased coefficients: If there are time-invariant confounders correlated with your predictors (heterogeneity bias), your estimates will be biased. This is the fundamental problem that both FE and REWB address.

2. Wrong standard errors: Even if coefficients happen to be unbiased, your standard errors will be too small because OLS doesn’t account for within-unit correlation. Your confidence intervals and p-values will be misleading.

3. You’re possibly posing the wrong answers to your questions: You can’t distinguish within from between effects, which are often substantively different and interesting.

The minimal responsible approach is to at least use cluster-robust standard errors with pooled OLS, but even this doesn’t solve the bias problem. Unless you have strong theoretical reasons to believe there is no heterogeneity bias and no interesting difference between within and between effects, use one of these approaches.

What about random slopes? When should I allow effects to vary across units?

This is an important extension we mentioned earlier but didn’t implement in our main examples. Random slopes allow the effect of a predictor to vary across higher-level units. In our school funding example, this would mean allowing the within-school effect of funding to differ from school to school—perhaps funding is more effective in some schools than others.

When to consider random slopes:

1. You have theoretical reasons to expect heterogeneous effects (e.g., policies might work differently in different contexts)
2. Your outcome shows different slopes across groups in exploratory plots
3. You want to test whether a group-level variable moderates a lower-level effect

You could implement this easily with:

In lme4, change (1 | school_id) to (1 + funding_centered | school_id). This allows both intercepts and the effect of centered funding to vary across schools. Be aware that this substantially increases model complexity—you’re now estimating the variance of the slopes and their covariance with the intercepts.

Important consideration:

It is definitely true that failing to include random slopes when they exist can lead to anti-conservative standard errors (confidence intervals that are too narrow). If you have a key predictor where heterogeneous effects are plausible, it’s worth testing whether allowing random slopes substantially changes your inferences. However, random slopes require sufficient data—you need enough within-unit observations to reliably estimate unit-specific slopes. With only 2-3 time points per unit, random slopes models often fail to converge or produce unreliable estimates.

Should I test whether within and between effects are different?

A common question is whether one should always include both within and between components in the model. The answer depends on your research question and the data. If \(\beta_W = \beta_B\), then the simpler standard RE model (which constrains them to be equal) may be adequate. However, we recommend starting with the REWB specification as the default, since you can then just formally test the hypothesis.

Simulation code

The code below defines the exact function used to generate school_panel.csv. Copy it into an R script, then call

panel_data <- simulate_school_data()

to reproduce the dataset.

simulate_school_data <- function(seed = 123, n_schools = 40, n_years = 12) {
    set.seed(seed)

    total_obs <- n_schools * n_years

    school_id <- rep(1:n_schools, each = n_years)
    year <- rep(1:n_years, times = n_schools)

    school_types <- sample(c("Public", "Private"),
        size = n_schools,
        replace = TRUE,
        prob = c(0.6, 0.4)
    )
    school_type <- rep(school_types, each = n_years)

    school_random_effects <- rnorm(n_schools, mean = 0, sd = 8)

    school_mean_funding <- rnorm(n_schools,
        mean = ifelse(school_types == "Private", 55, 45),
        sd = 12
    )
    funding_means <- rep(school_mean_funding, each = n_years)

    within_school_variation <- rnorm(total_obs, mean = 0, sd = 6)
    funding <- funding_means + within_school_variation

    within_effect <- 0.6
    between_effect <- 1.8
    school_type_effect <- 10

    test_scores <- 500 +
        rep(school_random_effects, each = n_years) +
        within_effect * (funding - funding_means) +
        between_effect * funding_means +
        school_type_effect * (school_type == "Private") +
        rnorm(total_obs, mean = 0, sd = 7)

    data.frame(
        school_id = as.factor(school_id),
        year = year,
        school_type = as.factor(school_type),
        funding = funding,
        test_score = test_scores
    )
}