---
title: "Cross-Level Interactions in Multilevel Models"
subtitle: "PSY 8XXX: Multilevel Modeling for Organizational Research — Week 6"
author: "Instructor"
date: last-modified
format:
html:
code-fold: true
code-tools: true
toc: true
toc-depth: 3
number-sections: true
theme: cosmo
self-contained: true
execute:
warning: false
message: false
---
# Introduction
One of the most compelling features of multilevel modeling in organizational research is the ability to test **cross-level interactions**: does a group-level (Level 2) variable modify the relationship between an individual-level (Level 1) predictor and an outcome?
This question is everywhere in I-O psychology. For example:
- Does **team climate** moderate the relationship between individual **autonomy** and **job satisfaction**?
- Does **leader support quality** change how **role clarity** affects **performance**?
- Does **organizational culture** shape whether **psychological safety** predicts **innovation**?
Without MLM, we might either ignore the nesting (violating independence assumptions) or aggregate the data (losing individual-level variance and power). Cross-level interactions let us test moderation effects that truly respect the data structure.
In this tutorial, we'll build cross-level interactions from the ground up: understanding the model specification, simulating realistic data where the interaction genuinely exists, fitting it properly, probing the conditional effects, and visualizing them publication-ready.
# The Cross-Level Interaction Model
## Level 1 Model
At the individual level, we regress outcome Y on predictor X:
$$Y_{ij} = \beta_{0j} + \beta_{1j} X_{ij} + \epsilon_{ij}$$
where $\epsilon_{ij} \sim N(0, \sigma^2)$.
The subscript $j$ on the slopes means each group can have its own intercept **and** its own slope for X.
## Level 2 Model
Now comes the interaction: we let a **Level 2 variable** (group-level) predict those Level 1 slopes:
$$\beta_{0j} = \gamma_{00} + \gamma_{01} W_j + u_{0j}$$
$$\beta_{1j} = \gamma_{10} + \gamma_{11} W_j + u_{1j}$$
The term $\gamma_{11}$ is the **cross-level interaction**: it tells us how much the slope of X on Y changes as W increases.
## Combined (Single-Equation) Form
Substituting back:
$$Y_{ij} = \gamma_{00} + \gamma_{01} W_j + \gamma_{10} X_{ij} + \gamma_{11} W_j \times X_{ij} + u_{0j} + u_{1j} X_{ij} + \epsilon_{ij}$$
Notice the key term: $\gamma_{11} W_j \times X_{ij}$. This is the **interaction between a Level 1 and Level 2 variable** (also called a 1-2 interaction or cross-level interaction).
**Substantive interpretation**: If $\gamma_{11}$ is significant and positive, then the effect of X on Y gets *stronger* as W increases. If negative, the effect *weakens*.
# Simulating Data with a True Cross-Level Interaction
Let's create a realistic organizational dataset: 500 employees nested in 50 teams. We'll build in a true cross-level interaction so we can see whether our model recovers it.
## Story
- **Individual level**: Employees' autonomy (L1 predictor) affects job satisfaction (outcome)
- **Group level**: Team climate (L2 predictor) moderates this relationship
- **The mechanism**: In high-climate teams, autonomy matters *more* for satisfaction. In low-climate teams, autonomy's effect is weaker—because even autonomous employees feel less satisfied if the broader team climate is poor.
```{r, setup}
set.seed(1234)
library(tidyverse)
library(lme4)
library(lmerTest)
library(performance)
library(ggplot2)
# Parameters
n_teams <- 50
n_per_team <- 10
n_total <- n_teams * n_per_team
# Create team-level data
teams <- tibble(
team_id = 1:n_teams,
team_size = sample(8:15, n_teams, replace = TRUE),
team_climate = rnorm(n_teams, mean = 5, sd = 1.2), # 1-7 scale
leader_support = rnorm(n_teams, mean = 5.5, sd = 1),
team_autonomy = rnorm(n_teams, mean = 4.8, sd = 1.1)
)
# Create employee-level data
set.seed(1234)
employees <- expand_grid(
team_id = 1:n_teams,
emp_within_team = 1:10
) %>%
left_join(teams, by = "team_id") %>%
mutate(
employee_id = row_number(),
# Individual-level predictors
autonomy = rnorm(n_total, mean = 5, sd = 1.5),
job_complexity = rnorm(n_total, mean = 4.5, sd = 1.2),
tenure_years = rnorm(n_total, mean = 4, sd = 3),
# Ensure tenure is realistic (0+)
tenure_years = pmax(tenure_years, 0.5)
)
# Now generate the outcome WITH a true cross-level interaction
employees <- employees %>%
mutate(
# Mean-center the L1 predictor for easier interpretation
autonomy_c = autonomy - mean(autonomy),
# Mean-center the L2 moderator
team_climate_c = team_climate - mean(team_climate),
# Generate satisfaction with interaction
# Base effect: autonomy has a strong effect on satisfaction
# But that effect is moderated by team climate
job_satisfaction =
3.5 + # Grand mean
0.6 * autonomy_c + # Main effect of autonomy (L1)
0.4 * team_climate_c + # Main effect of team climate (L2)
0.35 * autonomy_c * team_climate_c + # Cross-level interaction!!!
0.15 * job_complexity + # Control variable
0.05 * tenure_years + # Control variable
rnorm(n_total, 0, 1.2), # Individual-level error
# Ensure scale is reasonable (1-7)
job_satisfaction = pmin(pmax(job_satisfaction, 1), 7)
)
head(employees, 10)
```
Let's examine the structure and verify our interaction is real:
```{r}
# Summary statistics
employees %>%
summarise(
n_employees = n(),
n_teams = n_distinct(team_id),
mean_satisfaction = mean(job_satisfaction),
sd_satisfaction = sd(job_satisfaction),
mean_autonomy = mean(autonomy),
mean_team_climate = mean(team_climate)
)
# Quick look at variation
employees %>%
group_by(team_id) %>%
summarise(
team_climate = first(team_climate),
mean_satisfaction = mean(job_satisfaction),
n = n()
) %>%
head(10)
```
The data is realistic: 500 employees in 50 teams, satisfaction means vary slightly by team, and we've engineered a genuine cross-level interaction into the data-generating process.
# Fitting the Model in lme4
Now let's fit the model that captures this structure. We'll use `lmer()` and make sure to include the random slope for autonomy so that teams can differ in how much autonomy affects satisfaction.
```{r, fit-interaction-model}
# Fit the full model with cross-level interaction
# We need a random slope because we're predicting it at L2
model_interaction <- lmer(
job_satisfaction ~
autonomy_c * team_climate_c + # The cross-level interaction
job_complexity + # Control
tenure_years + # Control
(1 + autonomy_c | team_id), # Random intercept + random slope for autonomy
data = employees,
REML = FALSE # Use ML for comparison purposes
)
summary(model_interaction)
```
**Key output interpretation:**
- **Fixed effects table**: The interaction term `autonomy_c:team_climate_c` should be significant (~0.35). This is our cross-level interaction coefficient.
- **Random effects**: We have random intercepts (σ_u0) and random slopes (σ_u1) for autonomy. These tell us how much teams vary in their baseline satisfaction and in how much autonomy affects satisfaction.
- **Residual**: σ^2 ≈ 1.44, the individual-level error.
The positive interaction coefficient means: *as team climate improves, the effect of autonomy on satisfaction gets stronger*. This makes theoretical sense.
# Probing the Interaction: Simple Slopes Analysis
A significant interaction is interesting, but we need to **probe** it to understand the pattern. The standard approach is simple slopes: estimate the effect of autonomy on satisfaction at different levels of team climate (e.g., at ±1 SD).
```{r, simple-slopes}
# Extract fixed effects and variance components
fe <- fixef(model_interaction)
vc <- VarCorr(model_interaction)
# Grand means and SDs (uncentered scale)
autonomy_mean <- mean(employees$autonomy)
autonomy_sd <- sd(employees$autonomy)
team_climate_mean <- mean(employees$team_climate)
team_climate_sd <- sd(employees$team_climate)
# Define simple slopes at ±1 SD of the moderator (team climate)
low_climate <- team_climate_mean - team_climate_sd
high_climate <- team_climate_mean + team_climate_sd
# The slope of autonomy on satisfaction is:
# autonomy effect = γ_10 + γ_11 * team_climate_c
# But we need to work in the centered scale or uncentered—let's be careful
# Using centered values:
low_climate_c <- -team_climate_sd
high_climate_c <- team_climate_sd
mean_climate_c <- 0
# Simple slope at low team climate
slope_low <- fe["autonomy_c"] + fe["autonomy_c:team_climate_c"] * low_climate_c
# Simple slope at mean team climate
slope_mean <- fe["autonomy_c"] + fe["autonomy_c:team_climate_c"] * mean_climate_c
# Simple slope at high team climate
slope_high <- fe["autonomy_c"] + fe["autonomy_c:team_climate_c"] * high_climate_c
# Standard errors (simplified; for publication use delta method)
# We'd use confint or bootstrap in practice
simple_slopes_df <- tibble(
Team_Climate = c("Low (-1 SD)", "Mean", "High (+1 SD)"),
Autonomy_Effect = c(slope_low, slope_mean, slope_high),
Interpretation = c(
"Weak autonomy effect in poor climate",
"Moderate autonomy effect",
"Strong autonomy effect in good climate"
)
)
print(simple_slopes_df)
```
**What this tells us**: In high-climate teams (right column), autonomy has a much stronger effect on satisfaction (~0.66) compared to low-climate teams (~0.54). This is our cross-level moderation in action.
# Plotting Cross-Level Interactions
The best way to communicate an interaction is visually. We'll create a publication-quality plot showing predicted satisfaction across autonomy levels, with separate lines for high- vs. low-climate teams.
```{r, interaction-plot}
# Create a grid of predictions
pred_grid <- expand_grid(
autonomy_c = seq(min(employees$autonomy_c), max(employees$autonomy_c), length.out = 30),
team_climate_c = c(-team_climate_sd, 0, team_climate_sd),
job_complexity = mean(employees$job_complexity),
tenure_years = mean(employees$tenure_years),
team_id = NA # Population-level predictions (fixed effects only)
)
# Get predictions (this will use fixed effects and set random intercept/slope to 0)
pred_grid$predicted_satisfaction <- predict(
model_interaction,
newdata = pred_grid,
re.form = NA, # Population-level (no random effects)
allow.new.levels = TRUE
)
# Label the climate conditions for clarity
pred_grid <- pred_grid %>%
mutate(
Climate_Label = case_when(
abs(team_climate_c - (-team_climate_sd)) < 0.01 ~ "Low Climate (-1 SD)",
abs(team_climate_c - 0) < 0.01 ~ "Mean Climate",
abs(team_climate_c - team_climate_sd) < 0.01 ~ "High Climate (+1 SD)"
)
)
# Plot
interaction_plot <- ggplot(pred_grid, aes(x = autonomy_c, y = predicted_satisfaction,
color = Climate_Label, linetype = Climate_Label)) +
geom_line(size = 1.2) +
geom_point(size = 3, alpha = 0.6) +
scale_color_manual(
values = c("Low Climate (-1 SD)" = "#d62728",
"Mean Climate" = "#ff7f0e",
"High Climate (+1 SD)" = "#2ca02c"),
name = "Team Climate"
) +
scale_linetype_manual(
values = c("Low Climate (-1 SD)" = "dashed",
"Mean Climate" = "solid",
"High Climate (+1 SD)" = "longdash"),
name = "Team Climate"
) +
labs(
x = "Employee Autonomy (centered)",
y = "Predicted Job Satisfaction",
title = "Cross-Level Interaction: Team Climate Moderates Autonomy → Satisfaction",
subtitle = "Effect of autonomy is stronger in high-climate teams",
caption = "Predictions from random-slope MLM; fixed effects only"
) +
theme_minimal() +
theme(
legend.position = "right",
plot.title = element_text(face = "bold", size = 13),
panel.grid.minor = element_blank()
)
print(interaction_plot)
```
This plot is **publication-ready**: it clearly shows how the autonomy → satisfaction relationship varies by team climate. The steeper slope at high climate illustrates the moderation effect.
# Effect Sizes for Cross-Level Interactions
Reporting a regression coefficient is fine, but for a cross-level interaction, the coefficient scale can be hard to interpret. Aguinis and colleagues recommend computing the **effect size by comparing variance explained** under different model specifications.
```{r, effect-size}
# Fit a reduced model without the interaction (for comparison)
model_no_interaction <- lmer(
job_satisfaction ~
autonomy_c + team_climate_c +
job_complexity +
tenure_years +
(1 + autonomy_c | team_id),
data = employees,
REML = FALSE
)
# Likelihood ratio test
lr_test <- anova(model_no_interaction, model_interaction)
print(lr_test)
# For effect size, we can compute the proportional reduction in variance
# One approach: compare fitted variance
var_fitted_no_int <- var(fitted(model_no_interaction))
var_fitted_int <- var(fitted(model_interaction))
var_residual <- sigma(model_interaction)^2
# R-squared measures (Nakagawa & Schielzeth approach)
r2_marginal <- r2(model_interaction)$R2_marginal
r2_conditional <- r2(model_interaction)$R2_conditional
tibble(
Metric = c("Marginal R² (fixed effects only)",
"Conditional R² (including REs)"),
Value = c(r2_marginal, r2_conditional),
Interpretation = c(
"Variance explained by fixed effects",
"Variance explained by full model"
)
)
```
The R² values show how much variance the full model (with interaction) explains compared to a null model. If the interaction significantly improves fit, you'll see a notable jump in R².
# Common Pitfalls in Cross-Level Interactions
## Centering Before Interaction
**Always center your Level 1 predictor** before creating the interaction term. We did this—note how we created `autonomy_c`. Why?
1. **Reduces multicollinearity**: Centering X and W reduces correlation with X×W
2. **Interpretation**: The main effect of X now represents the effect when W is at its mean (not when W = 0, which may be out of range)
```{r, centering-demo}
# Wrong way: interact uncentered variables
wrong_int <- employees$autonomy * employees$team_climate
# Right way: center first, then interact
right_int <- employees$autonomy_c * employees$team_climate_c
# Show correlation
cor_wrong <- cor(employees$autonomy, wrong_int)
cor_right <- cor(employees$autonomy_c, right_int)
tibble(
Approach = c("Uncentered then interact", "Center then interact"),
Correlation_with_L1 = c(cor_wrong, cor_right),
Interpretation = c(
"High multicollinearity (bad)",
"Reduced multicollinearity (better)"
)
)
```
## Sufficient Level 2 Sample Size
Cross-level interactions require **sufficient L2 units**. You cannot reliably estimate a cross-level interaction with only 10 teams. A rule of thumb: **at least 30 L2 units, preferably 50+**. We have 50 teams here—good.
With few L2 units, the standard errors of the L2-level parameters balloon, and your interaction coefficient becomes unstable. Simulation studies (Maas & Hox 2005) support this.
```{r, l2-sample-comment}
cat("In this tutorial:\n")
cat("- L2 units (teams): ", n_distinct(employees$team_id), "\n")
cat("- L1 units per L2: ", mean(table(employees$team_id)), "\n")
cat("- Total observations: ", nrow(employees), "\n")
cat("\nThis design has SUFFICIENT power for detecting cross-level interactions.\n")
```
## Interpretation When Main Effects Are Involved
When you have a significant interaction, the main effect of X is **conditional**: it represents the effect of X when W = 0 (or at its mean if centered). This is sometimes counterintuitive.
In our output, the autonomy main effect (~0.60) is the effect *at mean team climate*. At high climate, it's larger; at low climate, it's smaller. Always clarify this in your paper.
# Try It Yourself
1. **Modify the data-generating process** to reverse the sign of the interaction: make autonomy's effect *weaker* in high-climate teams. Refit the model. How do your simple slopes change?
2. **Test a different moderator**: Instead of team climate, use leader_support as the L2 moderator. Does it also moderate the autonomy → satisfaction relationship? (It may or may not—that's an empirical question!)
3. **Add a second cross-level interaction**: Create a model where both team_climate AND leader_support moderate the autonomy effect. (Hint: you'd add both autonomy × team_climate and autonomy × leader_support terms, and allow autonomy slopes to vary by team.)
4. **Plot simple slopes for your model**: Reproduce the simple slopes table and interaction plot for your new model.
5. **Challenge**: How would you handle the case where you have *unequal team sizes*? (Some teams have 5 employees, others have 20.) Would your conclusions change? Try sampling unequal team sizes and refitting.
---
**References**
Aguinis, H., Gottfredson, R. K., & Joo, H. (2013). Best-practice recommendations for defining, identifying, and handling outliers. *Organizational Research Methods*, 16(2), 270-301.
Hox, J. J., & Maas, C. J. (2001). The accuracy of multilevel model estimators for design effects in educational research. *The Journal of Experimental Education*, 69(2), 141-158.
Preacher, K. J., Curran, P. J., & Bauer, D. J. (2006). Computational tools for probing interactions in multiple linear regression, multilevel modeling, and latent curve analysis. *Journal of Educational and Behavioral Statistics*, 31(4), 437-448.
