The parametric bootstrap and a boundary problem

resampling

GLM

ecology tutorial

Simulate from a fitted model to build a null distribution the asymptotic chi-squared gets wrong, and test Poisson against negative binomial overdispersion.

Author

Tidy Ecology

Published

2026-06-07

Modified

2026-09-30

Corrected 30 September 2026: the post said the parametric bootstrap hits the nominal level by construction and gives the right p-value; it now explains that the match on the same simulations is automatic and that the real level is only approximately nominal, in the text, the figure caption and the closing section. It also no longer reads the raw variance-to-mean ratio as a hint of overdispersion, since the gradient alone inflates that ratio.

The previous post resampled the data. The parametric bootstrap resamples from a model instead: fit the model you believe holds under the null hypothesis, simulate fresh datasets from it, and recompute your statistic on each. The spread of those simulated statistics is the null distribution you need. This is the right tool when the data-resampling bootstrap does not apply, and above all when a textbook asymptotic distribution is quietly wrong.

A common ecological case where the textbook answer is wrong: testing whether count data are overdispersed. We fit a Poisson generalised linear model and a negative binomial alternative, and compare them with a likelihood-ratio test. The trouble is that Poisson is negative binomial with the dispersion pushed to its limit, and that limit sits on the boundary of the parameter space. When a null value lies on a boundary, the likelihood-ratio statistic does not follow the chi-squared distribution everyone reaches for. The parametric bootstrap fixes it without any special theory.

Why the chi-squared is wrong here

The negative binomial has a dispersion parameter; write it so that the Poisson corresponds to one over that parameter being zero. Poisson is the negative binomial with the extra variance switched off. That zero is the smallest value the parameter can take, so under the null the truth sits at the edge of the allowed range. Standard likelihood theory assumes the null value is in the interior, with wiggle room on both sides. On a boundary that assumption fails, and the likelihood-ratio statistic no longer has an asymptotic chi-squared distribution on one degree of freedom. The correct large-sample answer is a fifty-fifty mixture of a point mass at zero and a chi-squared on one degree of freedom. Half the time the fitted dispersion wants to go negative, is held at zero, and the statistic is exactly zero.

Using the naive chi-squared ignores that spike at zero, so it places its critical value too far out and rejects too rarely. In overdispersion terms: the naive test is conservative and will tell you to keep Poisson when you should not.

Overdispersed counts along a gradient

We simulate abundance counts with modest overdispersion along an environmental gradient, then fit both models.

Setup: packages, colours and plot theme

suppressMessages({library(MASS); library(ggplot2)})
te_ink <- "#275139"; te_green <- "#4e7a5a"; te_gold <- "#c98a1b"
te_rose <- "#a5524a"; te_slate <- "#3f5c7a"; te_grey <- "#8a8a80"
theme_te <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = "#2f3e35"),
          plot.title = element_text(face = "bold", colour = "#233a2c", size = rel(1.02)),
          plot.subtitle = element_text(colour = "#54655b", size = rel(0.86)),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e3d8", linewidth = 0.3),
          plot.background = element_rect(fill = "#f5f4ee", colour = NA),
          panel.background = element_rect(fill = "#f5f4ee", colour = NA),
          axis.text = element_text(colour = "#54655b"),
          legend.position = "bottom", plot.margin = margin(8, 12, 8, 12))
}

set.seed(243)
n <- 90
x <- as.numeric(scale(runif(n, 0, 10)))     # standardised environmental gradient
mu <- exp(1.15 + 0.5 * x)
y <- rnbinom(n, size = 6, mu = mu)           # counts with mild overdispersion
disp_obs <- var(y) / mean(y)                 # crude variance-to-mean ratio

The counts have a variance-to-mean ratio of 2.13, but that crude ratio is inflated by the gradient itself: Poisson counts whose mean rises along the gradient already give a ratio well above one, so it says little about overdispersion. The question is whether there is real extra spread around the fitted mean, or only noise.

The likelihood-ratio statistic and its wrong p-value

The statistic is twice the gain in log-likelihood from Poisson to negative binomial. We clamp it at zero, since the negative binomial can never fit worse than the Poisson it contains.

lrt_stat <- function(y, x) {
  m_pois <- glm(y ~ x, family = poisson)
  m_nb <- suppressWarnings(tryCatch(glm.nb(y ~ x), error = function(e) NULL))
  if (is.null(m_nb)) return(0)
  max(0, 2 * (as.numeric(logLik(m_nb)) - as.numeric(logLik(m_pois))))
}
LRT_obs <- lrt_stat(y, x)
p_naive <- 1 - pchisq(LRT_obs, df = 1)         # naive chi-squared(1)
p_mix   <- 0.5 * (1 - pchisq(LRT_obs, df = 1)) # asymptotic 50:50 mixture

The observed statistic is 2.99. Fed through the naive chi-squared on one degree of freedom it gives a p-value of 0.084, while the boundary-aware mixture halves that to 0.042. One statistic, two p-values, and with a threshold at 0.05 the pair can land on opposite sides of it. The naive p-value is exactly twice the mixture one, because it forgets that half the null mass sits at zero.

The parametric bootstrap null

Rather than trust either asymptotic formula, we build the null distribution directly. Under the null the data are Poisson, so we simulate from the fitted Poisson model, refit both models to each simulated dataset, and collect the likelihood-ratio statistics.

m_pois <- glm(y ~ x, family = poisson)
mu_hat <- fitted(m_pois)
set.seed(909)
Bpb <- 1500
LRT_star <- replicate(Bpb, lrt_stat(rpois(n, mu_hat), x))
p_boot <- (1 + sum(LRT_star >= LRT_obs)) / (Bpb + 1)
frac0  <- mean(LRT_star < 1e-6)                # share of the null exactly at zero

The parametric bootstrap p-value is 0.022, agreeing with the mixture that the overdispersion is real, and if anything more decisive. It also shows why the boundary matters: 59% of the simulated null statistics are exactly zero, a point mass the chi-squared curve has no way to represent. That fraction is even a little above the asymptotic fifty per cent, because in samples this size the fitted dispersion is pulled to the boundary slightly more than half the time.

A histogram of simulated null likelihood-ratio statistics over the positive range, with a solid half chi-squared curve running a little above the bars across most of the range and pulling further above them to the right, and a dashed full chi-squared curve well above both. A vertical line marks the observed statistic between two marked critical values, and text notes that about 59 per cent of the null is a spike at zero. — Figure 1: The parametric bootstrap null distribution of the likelihood-ratio statistic (1500 simulations from the fitted Poisson). A large share of the null is a spike at zero (annotated), and the positive part is close to half a chi-squared on one degree of freedom (solid) but lighter at this sample size, and nowhere near the full chi-squared the naive test assumes (dashed), which sits far higher and so overstates the tail badly. The observed statistic (vertical line) lies beyond the mixture critical value but short of the naive one.

The naive test is conservative

Because the simulated statistics come from the null, they let us read off the naive test’s real error rate. A test with a nominal five per cent level should reject five per cent of the time when the null holds.

naive_size <- mean(LRT_star > qchisq(0.95, 1))   # actual size of the naive 5% test
alpha_grid <- seq(0.01, 0.20, by = 0.01)
size_naive <- sapply(alpha_grid, function(a) mean(LRT_star > qchisq(1 - a, 1)))

At its nominal five per cent, the naive chi-squared test actually rejects only 1.1% of the time under the null: it throws away most of its power to detect overdispersion in exchange for a false sense of caution. The parametric bootstrap reads its critical value from these same simulated statistics, so on them it rejects at the nominal rate automatically; that is not a check of its real size, which is only approximately nominal because the Poisson mean it simulates from is itself estimated.

Two panels. Left: a bar chart of three p-values with a dashed line at 0.05; the naive bar is above the line and the mixture and bootstrap bars are below it. Right: a line of actual size versus nominal size lying below the dashed identity line, showing the naive test rejects less often than advertised. — Figure 2: Left: the three p-values for the observed test. The naive chi-squared sits above 0.05 (keep Poisson) while the mixture and the parametric bootstrap sit below it (adopt the negative binomial). Right: the naive test’s actual rejection rate under the null against its nominal level; the curve stays below the diagonal everywhere, so the naive test is uniformly conservative. The parametric bootstrap is not drawn: it takes its critical values from these same simulations, so it would sit on the diagonal automatically.

What to take away

The parametric bootstrap earns its place whenever you can simulate from the model under the null but cannot trust the asymptotic reference distribution. Boundary parameters are the classic trap: testing a variance component against zero, testing overdispersion, testing a mixture component, all put the null value on an edge, and all break the usual chi-squared. Simulating from the fitted null model sidesteps the theory and gives an approximately calibrated p-value, here rescuing a real overdispersion signal that the naive test would have dismissed. The cost is that the parametric bootstrap trusts the null model completely: if the Poisson mean structure is itself misspecified, simulating from it inherits that error. Resample from the model you are willing to defend, and the next post shows what goes wrong when even the data-resampling bootstrap is applied to data that are not independent.

References

Efron 1979 Annals of Statistics 7(1):1-26 (10.1214/aos/1176344552)

Self and Liang 1987 Journal of the American Statistical Association 82(398):605-610 (10.1080/01621459.1987.10478472)

DiCiccio and Efron 1996 Statistical Science 11(3):189-228 (10.1214/ss/1032280214)

Davison and Hinkley 1997 Bootstrap Methods and their Application. ISBN 978-0-521-57391-7

Manly 2007 Randomization, Bootstrap and Monte Carlo Methods in Biology, third edition. ISBN 978-1-58488-541-2

--- title: "The parametric bootstrap and a boundary problem" description: "Simulate from a fitted model to build a null distribution the asymptotic chi-squared gets wrong, and test Poisson against negative binomial overdispersion." date: "2026-06-07 10:00" date-modified: "2026-09-30" categories: [R, resampling, GLM, ecology tutorial] image: thumbnail.png image-alt: "The null distribution of a likelihood-ratio statistic with a large spike at zero, matching a half chi-squared curve rather than a full one." --- *Corrected 30 September 2026: the post said the parametric bootstrap hits the nominal level by construction and gives the right p-value; it now explains that the match on the same simulations is automatic and that the real level is only approximately nominal, in the text, the figure caption and the closing section. It also no longer reads the raw variance-to-mean ratio as a hint of overdispersion, since the gradient alone inflates that ratio.* The previous post resampled the data. The parametric bootstrap resamples from a *model* instead: fit the model you believe holds under the null hypothesis, simulate fresh datasets from it, and recompute your statistic on each. The spread of those simulated statistics is the null distribution you need. This is the right tool when the data-resampling bootstrap does not apply, and above all when a textbook asymptotic distribution is quietly wrong. A common ecological case where the textbook answer is wrong: testing whether count data are overdispersed. We fit a Poisson generalised linear model and a negative binomial alternative, and compare them with a likelihood-ratio test. The trouble is that Poisson is negative binomial with the dispersion pushed to its limit, and that limit sits on the *boundary* of the parameter space. When a null value lies on a boundary, the likelihood-ratio statistic does not follow the chi-squared distribution everyone reaches for. The parametric bootstrap fixes it without any special theory. ## Why the chi-squared is wrong here The negative binomial has a dispersion parameter; write it so that the Poisson corresponds to one over that parameter being zero. Poisson is the negative binomial with the extra variance switched off. That zero is the smallest value the parameter can take, so under the null the truth sits at the edge of the allowed range. Standard likelihood theory assumes the null value is in the interior, with wiggle room on both sides. On a boundary that assumption fails, and the likelihood-ratio statistic no longer has an asymptotic chi-squared distribution on one degree of freedom. The correct large-sample answer is a fifty-fifty mixture of a point mass at zero and a chi-squared on one degree of freedom. Half the time the fitted dispersion wants to go negative, is held at zero, and the statistic is exactly zero. Using the naive chi-squared ignores that spike at zero, so it places its critical value too far out and rejects too rarely. In overdispersion terms: the naive test is conservative and will tell you to keep Poisson when you should not. ## Overdispersed counts along a gradient We simulate abundance counts with modest overdispersion along an environmental gradient, then fit both models. ```{r} #| label: setup #| message: false #| warning: false #| code-fold: true #| code-summary: "Setup: packages, colours and plot theme" suppressMessages({library(MASS); library(ggplot2)}) te_ink <- "#275139"; te_green <- "#4e7a5a"; te_gold <- "#c98a1b" te_rose <- "#a5524a"; te_slate <- "#3f5c7a"; te_grey <- "#8a8a80" theme_te <- function(base_size = 12) { theme_minimal(base_size = base_size) + theme(text = element_text(colour = "#2f3e35"), plot.title = element_text(face = "bold", colour = "#233a2c", size = rel(1.02)), plot.subtitle = element_text(colour = "#54655b", size = rel(0.86)), panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e3d8", linewidth = 0.3), plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA), axis.text = element_text(colour = "#54655b"), legend.position = "bottom", plot.margin = margin(8, 12, 8, 12)) } ``` ```{r} #| label: data set.seed(243) n <- 90 x <- as.numeric(scale(runif(n, 0, 10))) # standardised environmental gradient mu <- exp(1.15 + 0.5 * x) y <- rnbinom(n, size = 6, mu = mu) # counts with mild overdispersion disp_obs <- var(y) / mean(y) # crude variance-to-mean ratio ``` The counts have a variance-to-mean ratio of `r sprintf("%.2f", disp_obs)`, but that crude ratio is inflated by the gradient itself: Poisson counts whose mean rises along the gradient already give a ratio well above one, so it says little about overdispersion. The question is whether there is real extra spread around the fitted mean, or only noise. ## The likelihood-ratio statistic and its wrong p-value The statistic is twice the gain in log-likelihood from Poisson to negative binomial. We clamp it at zero, since the negative binomial can never fit worse than the Poisson it contains. ```{r} #| label: fit-lrt lrt_stat <- function(y, x) { m_pois <- glm(y ~ x, family = poisson) m_nb <- suppressWarnings(tryCatch(glm.nb(y ~ x), error = function(e) NULL)) if (is.null(m_nb)) return(0) max(0, 2 * (as.numeric(logLik(m_nb)) - as.numeric(logLik(m_pois)))) } LRT_obs <- lrt_stat(y, x) p_naive <- 1 - pchisq(LRT_obs, df = 1) # naive chi-squared(1) p_mix <- 0.5 * (1 - pchisq(LRT_obs, df = 1)) # asymptotic 50:50 mixture ``` The observed statistic is `r sprintf("%.2f", LRT_obs)`. Fed through the naive chi-squared on one degree of freedom it gives a p-value of `r sprintf("%.3f", p_naive)`, while the boundary-aware mixture halves that to `r sprintf("%.3f", p_mix)`. One statistic, two p-values, and with a threshold at `0.05` the pair can land on opposite sides of it. The naive p-value is exactly twice the mixture one, because it forgets that half the null mass sits at zero. ## The parametric bootstrap null Rather than trust either asymptotic formula, we build the null distribution directly. Under the null the data are Poisson, so we simulate from the fitted Poisson model, refit both models to each simulated dataset, and collect the likelihood-ratio statistics. ```{r} #| label: pboot #| message: false #| warning: false m_pois <- glm(y ~ x, family = poisson) mu_hat <- fitted(m_pois) set.seed(909) Bpb <- 1500 LRT_star <- replicate(Bpb, lrt_stat(rpois(n, mu_hat), x)) p_boot <- (1 + sum(LRT_star >= LRT_obs)) / (Bpb + 1) frac0 <- mean(LRT_star < 1e-6) # share of the null exactly at zero ``` The parametric bootstrap p-value is `r sprintf("%.3f", p_boot)`, agreeing with the mixture that the overdispersion is real, and if anything more decisive. It also shows why the boundary matters: `r sprintf("%.0f", 100*frac0)`% of the simulated null statistics are exactly zero, a point mass the chi-squared curve has no way to represent. That fraction is even a little above the asymptotic fifty per cent, because in samples this size the fitted dispersion is pulled to the boundary slightly more than half the time. ```{r} #| label: fig-null #| echo: false #| fig-cap: "The parametric bootstrap null distribution of the likelihood-ratio statistic (1500 simulations from the fitted Poisson). A large share of the null is a spike at zero (annotated), and the positive part is close to half a chi-squared on one degree of freedom (solid) but lighter at this sample size, and nowhere near the full chi-squared the naive test assumes (dashed), which sits far higher and so overstates the tail badly. The observed statistic (vertical line) lies beyond the mixture critical value but short of the naive one." #| fig-alt: "A histogram of simulated null likelihood-ratio statistics over the positive range, with a solid half chi-squared curve running a little above the bars across most of the range and pulling further above them to the right, and a dashed full chi-squared curve well above both. A vertical line marks the observed statistic between two marked critical values, and text notes that about 59 per cent of the null is a spike at zero." #| fig-width: 7 #| fig-height: 4.4 pos <- LRT_star[LRT_star > 1e-6] dd <- data.frame(v = pos) crit_naive <- qchisq(0.95, 1); crit_mix <- qchisq(0.90, 1) ggplot(dd, aes(v)) + geom_histogram(aes(y = after_stat(density) * (1 - frac0)), bins = 45, fill = "#dfe8e0", colour = "#c4d3c6", linewidth = 0.2) + stat_function(fun = function(z) 0.5 * dchisq(z, 1), colour = te_green, linewidth = 0.9) + stat_function(fun = function(z) dchisq(z, 1), colour = te_rose, linewidth = 0.7, linetype = "dashed") + geom_vline(xintercept = LRT_obs, colour = te_ink, linewidth = 0.8) + geom_vline(xintercept = crit_mix, colour = te_green, linewidth = 0.5, linetype = "dotted") + geom_vline(xintercept = crit_naive, colour = te_rose, linewidth = 0.5, linetype = "dotted") + annotate("text", x = 6.2, y = 0.34, label = sprintf("spike at zero: %.0f%%", 100 * frac0), hjust = 0, size = 3.1, colour = "#233a2c") + annotate("text", x = 6.2, y = 0.30, label = "solid: half chi-squared (correct)", hjust = 0, size = 3, colour = te_green) + annotate("text", x = 6.2, y = 0.265, label = "dashed: full chi-squared (naive)", hjust = 0, size = 3, colour = te_rose) + annotate("text", x = LRT_obs + 0.15, y = 0.11, label = "observed", hjust = 0, size = 3, colour = te_ink) + coord_cartesian(xlim = c(0, 12), ylim = c(0, 0.42)) + labs(title = "The null is a mixture, not a chi-squared", subtitle = "Positive part of 1500 simulated null statistics, with the two asymptotic curves", x = "likelihood-ratio statistic", y = "density (scaled to positive mass)") + theme_te() + theme(legend.position = "none") ``` ## The naive test is conservative Because the simulated statistics come from the null, they let us read off the naive test's real error rate. A test with a nominal five per cent level should reject five per cent of the time when the null holds. ```{r} #| label: calib naive_size <- mean(LRT_star > qchisq(0.95, 1)) # actual size of the naive 5% test alpha_grid <- seq(0.01, 0.20, by = 0.01) size_naive <- sapply(alpha_grid, function(a) mean(LRT_star > qchisq(1 - a, 1))) ``` At its nominal five per cent, the naive chi-squared test actually rejects only `r sprintf("%.1f", 100*naive_size)`% of the time under the null: it throws away most of its power to detect overdispersion in exchange for a false sense of caution. The parametric bootstrap reads its critical value from these same simulated statistics, so on them it rejects at the nominal rate automatically; that is not a check of its real size, which is only approximately nominal because the Poisson mean it simulates from is itself estimated. ```{r} #| label: fig-calib #| echo: false #| fig-cap: "Left: the three p-values for the observed test. The naive chi-squared sits above 0.05 (keep Poisson) while the mixture and the parametric bootstrap sit below it (adopt the negative binomial). Right: the naive test's actual rejection rate under the null against its nominal level; the curve stays below the diagonal everywhere, so the naive test is uniformly conservative. The parametric bootstrap is not drawn: it takes its critical values from these same simulations, so it would sit on the diagonal automatically." #| fig-alt: "Two panels. Left: a bar chart of three p-values with a dashed line at 0.05; the naive bar is above the line and the mixture and bootstrap bars are below it. Right: a line of actual size versus nominal size lying below the dashed identity line, showing the naive test rejects less often than advertised." #| fig-width: 7 #| fig-height: 3.9 dp <- data.frame(method = factor(c("naive\nchi-squared", "0.5:0.5\nmixture", "parametric\nbootstrap"), levels = c("naive\nchi-squared", "0.5:0.5\nmixture", "parametric\nbootstrap")), p = c(p_naive, p_mix, p_boot), col = c(te_rose, te_gold, te_ink)) pA <- ggplot(dp, aes(method, p, fill = col)) + geom_col(width = 0.66) + geom_hline(yintercept = 0.05, linetype = "dashed", colour = "#233a2c", linewidth = 0.6) + geom_text(aes(label = sprintf("%.3f", p)), vjust = -0.5, size = 3, colour = "#233a2c") + scale_fill_identity() + coord_cartesian(ylim = c(0, 0.095)) + labs(title = "One statistic, three verdicts", x = NULL, y = "p-value") + theme_te() + theme(axis.text.x = element_text(size = 7.5)) dcal <- data.frame(nominal = alpha_grid, actual = size_naive) pB <- ggplot(dcal, aes(nominal, actual)) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "#8a8a80", linewidth = 0.5) + geom_line(colour = te_rose, linewidth = 0.9) + geom_point(colour = te_rose, size = 1.3) + labs(title = "Naive test is conservative", subtitle = "Actual vs nominal rejection rate", x = "nominal level", y = "actual size under null") + theme_te() + theme(legend.position = "none") grid::grid.newpage() grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(pA, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## What to take away The parametric bootstrap earns its place whenever you can simulate from the model under the null but cannot trust the asymptotic reference distribution. Boundary parameters are the classic trap: testing a variance component against zero, testing overdispersion, testing a mixture component, all put the null value on an edge, and all break the usual chi-squared. Simulating from the fitted null model sidesteps the theory and gives an approximately calibrated p-value, here rescuing a real overdispersion signal that the naive test would have dismissed. The cost is that the parametric bootstrap trusts the null model completely: if the Poisson mean structure is itself misspecified, simulating from it inherits that error. Resample from the model you are willing to defend, and the next post shows what goes wrong when even the data-resampling bootstrap is applied to data that are not independent. ## References Efron 1979 Annals of Statistics 7(1):1-26 (10.1214/aos/1176344552) Self and Liang 1987 Journal of the American Statistical Association 82(398):605-610 (10.1080/01621459.1987.10478472) DiCiccio and Efron 1996 Statistical Science 11(3):189-228 (10.1214/ss/1032280214) Davison and Hinkley 1997 Bootstrap Methods and their Application. ISBN 978-0-521-57391-7 Manly 2007 Randomization, Bootstrap and Monte Carlo Methods in Biology, third edition. ISBN 978-1-58488-541-2 ## Related tutorials - [The nonparametric bootstrap](../the-nonparametric-bootstrap/) - [Vuong test for zero inflation, and what to use instead](../vuong-test-for-zero-inflation/) - [Checking a bootstrap](../checking-a-bootstrap/) - [Testing a fitted distribution](../goodness-of-fit-with-estimated-parameters/)