You've run your model. The output looks clean. The coefficients make sense. Now you just need that one number — the p-value — to tell you whether what you're seeing is real or just noise.
Here's the thing: R gives you p-values in about a dozen different ways, depending on what you're testing and which package you're using. So non-standard models. You'll hit edge cases. Most tutorials show you one method and call it a day. Messy data. But in practice? Situations where the default output doesn't have what you need.
This guide covers the main ways to extract p-values in R — and the traps that catch people off guard.
What Is a P-Value (Quick Refresher)
Skip this if you know it cold. But a surprising number of people use p-values daily without being able to explain them in plain English.
A p-value answers: If the null hypothesis were true, how likely would I be to see data this extreme (or more extreme) just by chance?*
That's it. Because of that, low p-value = "this would be really unlikely if nothing was going on. " High p-value = "totally plausible that this is just random variation.
The conventional threshold is 0.A p-value of 0.051 isn't meaningfully different from 0.049. 05. But that's arbitrary. Anyone who treats 0.05 as a magic line is doing science a disservice.
The Basics: Built-In Model Summaries
Most of the time, you don't need to calculate* a p-value. You just need to find* it in the output.
Linear Models (lm)
model <- lm(mpg ~ wt + hp, data = mtcars)
summary(model)
Look at the Pr(>|t|) column. Those are your p-values for each coefficient. The intercept gets one too — though honestly, nobody cares about the intercept p-value unless you're testing whether the baseline is different from zero.
Generalized Linear Models (glm)
Same deal:
model <- glm(am ~ mpg + wt, data = mtcars, family = binomial)
summary(model)
Now you're looking at Pr(>|z|) — z-statistics instead of t, because GLMs use Wald tests by default. No workaround needed.
ANOVA Tables (aov)
model <- aov(mpg ~ factor(cyl), data = mtcars)
summary(model)
This gives you an F-test p-value for the overall factor. So different question: "Does cylinder count as a whole* explain variation in mpg? " rather than "Is this specific coefficient non-zero?
Extracting P-Values Programmatically
Interactive work? Because of that, summary() is fine. But if you're writing a function, building a table, or running 500 models in a loop, you need to pull p-values out as objects.
The broom Package (Best All-Around)
library(broom)
tidy(model)
Returns a tibble with term, estimate, std.On the flip side, clean. In practice, value. On top of that, error, statistic, p. Consistent across model types. Works with lm, glm, lmer, glmer, coxph, survreg, and dozens more.
If you only want the p-values:
tidy(model)$p.value
Base R Extraction (No Dependencies)
Sometimes you can't add packages. Here's how to dig them out manually:
# For lm/glm
coef(summary(model))[, "Pr(>|t|)"] # lm
coef(summary(model))[, "Pr(>|z|)"] # glm
# For aov
summary(model)[[1]][["Pr(>F)"]][1]
The aov one is ugly. summary.aov returns a list of data frames. The first element is the ANOVA table. The p-value is in the Pr(>F) column, first row. You'll forget this syntax. Bookmark it.
Hypothesis Tests Outside Regression
Not everything is a regression coefficient. Here are the common ones:
t-Tests
t.test(mpg ~ am, data = mtcars)$p.value
Two-sample, paired, one-sample — all return a p.value element in the list.
Chi-Square Tests
chisq.test(table(mtcars$cyl, mtcars$am))$p.value
Wilcoxon / Mann-Whitney
wilcox.test(mpg ~ am, data = mtcars)$p.value
Correlation Tests
cor.test(mtcars$mpg, mtcars$wt)$p.value
Shapiro-Wilk (Normality)
shapiro.test(mtcars$mpg)$p.value
Mixed Models and Complex Objects
This is where people get stuck. lme4 models (lmer, glmer) don't have p-values in summary() by design. The authors argue — correctly — that denominator degrees of freedom are ambiguous in mixed models.
Option 1: lmerTest (Adds p-values to summary)
library(lmerTest)
model <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
summary(model)
Now you get Pr(>|t|) with Satterthwaite or Kenward-Roger df approximation. Easiest path if you just want the table.
If you found this helpful, you might also enjoy what is the smell of rain called or what is gummy candy made of.
Option 2: broom.mixed
library(broom.mixed)
tidy(model, effects = "fixed")
Returns p-values using the same approximations. Cleaner for programmatic extraction.
Option 3: Likelihood Ratio Tests (More Principled)
model_full <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
model_null <- lmer(Reaction ~ 1 + (Days | Subject), data = sleepstudy)
anova(model_null, model_full)
Compares nested models. Which means gives a chi-square p-value. This is often better* than coefficient-wise p-values for mixed models — but it answers a different question.
Multiple Testing: The Problem Nobody Talks About Enough
Run 20 tests at α = 0.05. Expect one false positive by definition*.
R makes correction easy. The hard part is remembering to do it.
Bonferroni (Conservative)
p.adjust(p_values, method = "bonferroni")
Benjamini-Hochberg (FDR Control, Standard for Genomics/High-Dim)
p.adjust(p_values, method = "BH")
Holm (Step-Down Bonferroni, Uniformly More Powerful)
p.adjust(p_values, method = "holm")
Use BH for exploratory work with many tests. Bonferroni when you really* can't afford a false positive. Holm is a solid default — it's never worse than Bonferroni and often better.
Permutation and Bootstrap P-Values
When assumptions fail (non-normality, small samples, weird statistics), asymptotic p-values lie. Resampling saves you.
Permutation Test (Two Groups)
obs_diff <- mean(mtcars$mpg[mtcars$am == 1]) - mean(mtcars$mpg[mtcars$am == 0])
perm_diffs <- replicate(10000, {
shuffled <- sample(mtcars$am)
mean(mtcars$mpg[shuffled == 1]) - mean(mtcars$mpg[shuffled == 0])
})
p_val <- mean(abs(perm_diffs) >= abs(obs_diff
```r
p_val <- mean(abs(perm_diffs) >= abs(obs_diff))
This permutation p‑value estimates the probability of observing a difference as extreme as the one in the data under the null hypothesis that group labels are exchangeable. Because of that, with 10 000 replications the Monte‑Carlo error is roughly ±0. 01, which is usually sufficient for exploratory work; increase the number of replicates for stricter precision.
Bootstrap Confidence Intervals and p‑Values
When the sampling distribution of a statistic is unknown or skewed, the bootstrap provides an empirical alternative. A common approach is to compute a two‑sided p‑value from the proportion of bootstrap replicates that fall on the opposite side of zero relative to the observed estimate.
set.seed(123)
boot_est <- replicate(5000, {
idx <- sample(seq_len(nrow(mtcars)), replace = TRUE)
diff <- mean(mtcars$mpg[idx][mtcars$am[idx] == 1]) -
mean(mtcars$mpg[idx][mtcars$am[idx] == 0])
diff
})
# Two‑sided bootstrap p‑value
p_boot <- 2 * min(
mean(boot_est <= 0),
mean(boot_est >= 0)
)
p_boot
The boot package streamlines this workflow and offers bias‑corrected and accelerated (BCa) intervals:
library(boot)
boot_fun <- function(data, indices) {
d <- data[indices, ]
mean(d$mpg[d$am == 1]) - mean(d$mpg[d$am == 0])
}
boot_out <- boot(mtcars, boot_fun, R = 5000)
boot.ci(boot_out, type = c("perc", "bca"))
The resulting confidence interval can be inverted to obtain a p‑value: if the interval does not contain zero, the corresponding two‑sided p‑value is < 0.05.
Choosing Between Asymptotic, Permutation, and Bootstrap p‑Values
| Situation | Recommended approach |
|---|---|
| Large sample, standard test assumptions satisfied | Classical parametric p‑values (t, χ², F) |
| Small sample or obvious non‑normality | Permutation test (exact under exchangeability) |
| Complex statistic, no closed‑form null | Bootstrap (percentile or BCa) |
| Mixed‑effects models with ambiguous df | Likelihood‑ratio test or lmerTest/broom.mixed approximations |
| Many simultaneous tests | Adjust raw p‑values with BH (FDR) or Holm |
Practical Tips
- Set a seed before any resampling to ensure reproducibility.
- Check convergence of mixed models; warnings often indicate that the approximation used for p‑values may be unreliable.
- Report both the raw statistic and the method used to obtain the p‑value (e.g., “Welch’s t‑test, p = 0.032” or “Permutation test (10 000 reps), p = 0.018”).
- Document the number of resamples; 5 000–10 000 is a good default, increase if the p‑value is near your decision threshold.
- Visualise the resampling distribution (histogram of
perm_diffsorboot_est) alongside the observed value to convey the intuition behind the p‑value.
Conclusion
Extracting p‑values in R is straightforward for many classic tests, but modern analyses frequently involve models or data structures where the usual asymptotic approximations break down or are ambiguous. Practically speaking, mixed, and the built‑in resampling tools, we can obtain reliable p‑values through Satterthwaite/Kenward‑Roger degrees of freedom, likelihood‑ratio comparisons, permutation, or bootstrap methods. At the end of the day, the choice of method should be guided by the study design, sample size, and the specific hypothesis being tested—never by convenience alone. Because of that, coupled with thoughtful multiple‑testing correction, these strategies make sure the inferential statements we make are both statistically sound and transparent to readers. By leveraging packages such as lmerTest, broom.When assumptions are doubtful, resampling provides a dependable, intuitive fallback that aligns the p‑value directly with the empirical data at hand.
