Mastering Nonparametric Repeated Measures: A Comprehensive Guide to the Friedman Test in R

In the evolving landscape of data science and biostatistics, researchers frequently encounter datasets that defy the strict assumptions of traditional parametric testing. When analyzing quantitative variables measured across multiple time points or related conditions, the repeated measures ANOVA is often the default choice. However, when data fail to meet the stringent criteria of normality and sphericity—or when dealing with ordinal scales—analysts must pivot toward robust alternatives.

Enter the Friedman test: the nonparametric counterpart to the repeated measures ANOVA. Originally introduced by economist Milton Friedman in his landmark 1937 paper, "The use of ranks to avoid the assumption of normality implicit in the analysis of variance," this statistical procedure relies on ranks rather than raw values. By evaluating measurements within individual subjects, it accounts for within-subject dependency without demanding normally distributed residuals.

This article provides an in-depth, end-to-end guide on executing, interpreting, and visualizing the Friedman test in R, complete with post-hoc analyses and real-world clinical context.

Main Facts: Understanding the Friedman Test Framework

To contextualize the Friedman test within broader statistical methodologies, it is helpful to understand how it maps against other common omnibus tests used for comparing three or more groups:

While the Kruskal-Wallis test pools and ranks all observations together across independent samples, the Friedman test isolates its ranking within each subject. This intra-subject ranking mechanism effectively neutralizes baseline differences between participants, mirroring how a repeated measures ANOVA removes inter-subject variability from its error term.

Core Hypotheses

• Null Hypothesis ($H_0$): The distributions of the quantitative or ordinal variable are identical across all $k$ related conditions or time points.
• Alternative Hypothesis ($H_1$): At least one of the conditions tends to yield different values compared to the others.

Crucially, just like an ANOVA or a Kruskal-Wallis test, a significant Friedman test is an omnibus test. It signals that a difference exists among the groups, but identifying precisely which groups differ requires subsequent post-hoc testing.

Chronology and Evolution of Nonparametric Repeated Measures

The historical trajectory of statistical testing shifted dramatically in the early 20th century as researchers sought alternatives to strict Gaussian (normal distribution) assumptions.

1. 1920s–1930s: The development of Analysis of Variance (ANOVA) by Ronald Fisher revolutionized experimental design, providing robust methods to test group differences. However, these parametric tests assumed normal distributions and homoscedasticity (equal variances).
2. 1937: Economist Milton Friedman published his seminal work introducing rank-based methods for analysis of variance, allowing researchers to circumvent the normality assumption entirely in related-samples designs.
3. Late 20th Century: The popularization of computational languages like S and subsequently R transformed theoretical statistics into accessible, applied tools. Packages such as rstatix, PMCMRplus, and ggstatsplot integrated these classical rank-based procedures into modern tidyverse workflows.
4. Modern Era: Statisticians continue to refine post-hoc methodologies (such as Benavoli et al., 2016, critiquing mean-rank post-hoc tests), reinforcing the importance of pairing omnibus evaluations with reliable pairwise comparisons like Wilcoxon signed-rank tests or Conover tests.

Supporting Data: A Clinical Trial Simulation in R

To demonstrate the practical application of the Friedman test, let us examine a simulated clinical pilot study involving 16 chronic pain patients. Pain scores are recorded on a continuous scale from 0 (no pain) to 100 (unbearable pain) across three distinct milestones: before, during, and one month after a novel treatment.

Unlike standard textbook examples that assume clean normal distributions, this dataset simulates real-world clinical anomalies, including right-skewed flare-ups governed by an exponential distribution (rexp()).

# Set seed for reproducibility
set.seed(42)

# Number of patients in the pilot study
n <- 16

# Baseline patient-specific effects
patient_effect <- rnorm(n, mean = 0, sd = 6)

# Pain scores with right-skewed fluctuations (flare-ups)
before <- 60 + patient_effect + rexp(n, rate = 1 / 6)
during <- 48 + patient_effect + rexp(n, rate = 1 / 6)
after  <- 46 + patient_effect + rexp(n, rate = 1 / 6)

# Constructing the dataset in long format
dat <- data.frame(
  patient = factor(rep(1:n, times = 3)),
  time = factor(rep(c("before", "during", "after"), each = n),
    levels = c("before", "during", "after")
  ),
  pain = round(c(before, during, after), 1)
)

head(dat)

Descriptive Statistics

Because rank-based tests prioritize medians and interquartile ranges (IQRs) over means and standard deviations, we summarize our data accordingly:

library(dplyr)

dat %>%
  group_by(time) %>%
  summarise(
    n = n(),
    median = median(pain),
    IQR = IQR(pain),
    mean = round(mean(pain), 2)
  ) %>%
  as.data.frame()

• Before treatment: Median = 68.50, IQR = 15.33, Mean = 71.14
• During treatment: Median = 57.45, IQR = 12.03, Mean = 59.45
• After treatment: Median = 52.10, IQR = 7.48, Mean = 55.54

The raw descriptive metrics indicate a downward trend in pain scores following the intervention. However, evaluating whether this trend holds true for the broader population requires rigorous diagnostic and statistical testing.

Evaluating Assumptions and Executing the Test in R

Checking Normality of Residuals

Although the Friedman test is nonparametric and does not require normality, comparing it against a linear model helps justify why a nonparametric approach is necessary.

library(car)

# Linear model incorporating patient and time effects
res_lm <- lm(pain ~ time + patient, data = dat)

# Shapiro-Wilk normality test on residuals
shapiro.test(residuals(res_lm))

The Shapiro-Wilk test yields a $p$-value well below 0.05 ($p = 0.0004984$), and the Q-plot deviates significantly from the theoretical confidence bands. Given our small sample size ($n = 16$), the Central Limit Theorem cannot salvage these non-normal residuals, confirming that the Friedman test is the correct methodological choice.

Performing the Friedman Test

We can execute the Friedman test in base R using the formula interface (dependent ~ condition | subject) or via the rstatix package for tidy-compatible data frames:

# Base R implementation
friedman.test(pain ~ time | patient, data = dat)

# Tidyverse implementation via rstatix
library(rstatix)

dat %>%
  friedman_test(pain ~ time | patient)

# Calculating effect size (Kendall’s W)
dat %>%
  friedman_effsize(pain ~ time | patient)

Output Interpretation:

• Test Statistic ($Q$ / Friedman chi-squared): $17.375$
• Degrees of Freedom ($textdf$): $2$
• $p$-value: $0.0001689$
• Kendall’s $W$ (Effect Size): $0.543$ (interpreted as a large effect size, where values $ge 0.5$ indicate strong concordance across conditions).

Because the $p$-value is substantially lower than the standard significance threshold ($alpha = 0.05$), we reject the null hypothesis and conclude that pain scores vary significantly across the three measurement periods ($chi^2(2) = 17.38, p < 0.001, W = 0.54$).

Post-Hoc Testing: Uncovering Pairwise Differences

To determine which specific time points differ from one another, we employ adjusted pairwise post-hoc comparisons.

1. Pairwise Wilcoxon Signed-Rank Tests

dat %>%
  pairwise_wilcox_test(
    pain ~ time,
    paired = TRUE,
    p.adjust.method = "holm"
  )

• Before vs. During: Adjusted $p = 0.0184$ (*)
• Before vs. After: Adjusted $p = 0.0101$ (*)
• During vs. After: Adjusted $p = 0.1750$ ($textns$)

2. Conover Test

Using the PMCMRplus package, we can execute the Conover test—which estimates error terms via ranked $t$-distributions—yielding consistent conclusions:

library(PMCMRplus)

frdAllPairsConoverTest(
  y = dat$pain,
  groups = dat$time,
  blocks = dat$patient,
  p.adjust.method = "holm"
)

Implications and Visualizing Results

The clinical implications of these findings are profound. The novel treatment is associated with an immediate, statistically significant reduction in median chronic pain scores (dropping from $68.50$ down to $57.45$). Furthermore, this therapeutic benefit remains stable one month post-treatment (median $52.10$), with no statistically significant rebound occurring between the active treatment phase and the one-month follow-up.

To communicate these findings effectively, researchers can leverage the ggstatsplot package to combine descriptive boxplots, individual subject trajectories, and exact test statistics into a single cohesive visualization:

library(ggstatsplot)

ggwithinstats(
  data = dat,
  x = time,
  y = pain,
  subject.id = patient,
  type = "nonparametric",
  pairwise.display = "all",
  bf.message = FALSE
)

Summary

When analyzing repeated measures data characterized by non-normal distributions, ordinal scaling, or severe outliers (such as clinical flare-ups), the Friedman test serves as an indispensable tool in the data scientist’s toolkit. By evaluating intra-subject ranks rather than raw parametric values, it preserves statistical power while sidestepping the rigid assumptions of sphericity and normality. Paired with robust post-hoc procedures like the Conover test or pairwise Wilcoxon signed-rank tests, it provides a comprehensive framework for longitudinal data analysis in R.