Answer

How do you simulate a prediction interval from a fitted regression model?

Inspired by a question on Cross Validated · November 6, 2026

regressionconfidence intervals

Two sources of uncertainty

A prediction interval answers the question: where is a future observation likely to fall? Suppose the model is a normal linear regression, possibly on a log scale (for example log survival time, as in a log-normal AFT survival model). A new value at predictor value x₀ differs from your fitted prediction for two reasons.

• Estimation error. The intercept, slope and residual SD were estimated from a finite sample, so the true mean at x₀ could be somewhat higher or lower than the fitted line says.
• Individual scatter. Even if you knew the true mean exactly, a single new case would still land somewhere around it, with spread given by the residual SD.

A confidence interval for the mean covers only the first source. A prediction interval must cover both. Simulation handles both naturally: each draw picks one plausible version of the model, then generates a new outcome from that version.

The recipe, for one case and for a group

1. Fit the model and keep the estimates and their estimated covariance matrix.
2. Draw one set of parameters: the residual SD and the coefficients. In a likelihood-based fit, a common choice is a multivariate normal draw centred on the estimates of the coefficients and log SD, using the estimated covariance matrix.
3. Compute the model's mean at the new predictor value with those drawn coefficients.
4. For one new case, add a single noise term: the drawn SD times a standard normal value. For a group of m new cases with the same predictor value, add m separate noise terms, one per case, and average the m simulated outcomes.
5. Repeat steps 2 to 4 many thousands of times and keep the middle 95% of the results as the interval.

The detail that matters most for the group version is that the parameters are drawn once per repetition and shared by everyone in the group. All m new cases live in the same world: if the true slope is steeper than you estimated, it is steeper for all of them. The noise, by contrast, is independent from case to case, so it partly averages out. That is why the interval for a group average shrinks as m grows, but never below the confidence interval for the mean: the estimation error does not average away.

If you skip step 2 and plug in the point estimates, you get a plug-in interval that ignores estimation error. For a single case the damage is usually modest, because individual scatter dominates. For a group average it can be serious, because the scatter has been divided by the square root of m while the estimation error has not.

See it in R and Python

This example fits a straight line to 40 simulated points and predicts at x = 8, both for one new case and for the mean of 10 new cases. For an ordinary linear model the exact intervals have a closed form based on the t distribution, so we can check the simulation against them. Here the residual SD is drawn from its scaled inverse chi-square distribution and the coefficients are then drawn given that SD; for this model, that reproduces the exact t-based interval apart from simulation noise. The plug-in version is shown for contrast.

set.seed(2026)
n <- 40
x <- runif(n, 0, 10)
y <- 1 + 0.3 * x + rnorm(n, sd = 0.8)      # e.g. a log-scale outcome
fit <- lm(y ~ x)
x0 <- 8; m <- 10                            # new x value, size of new group
s  <- summary(fit)$sigma; df <- fit$df.residual
pf <- predict(fit, data.frame(x = x0), se.fit = TRUE)
tq <- qt(0.975, df)

# Exact (t-based) 95% intervals
exact_one <- pf$fit + c(-1, 1) * tq * sqrt(pf$se.fit^2 + s^2)
exact_bar <- pf$fit + c(-1, 1) * tq * sqrt(pf$se.fit^2 + s^2 / m)

# Simulation: draw sigma, then coefficients given sigma, then new outcomes
S <- 20000
XtXi  <- solve(crossprod(model.matrix(fit)))
sig   <- sqrt(df * s^2 / rchisq(S, df))
z     <- matrix(rnorm(2 * S), S, 2) %*% chol(XtXi)
beta  <- sweep(z * sig, 2, coef(fit), "+")
mu0   <- beta[, 1] + beta[, 2] * x0                  # one mean per draw
y_one <- mu0 + sig * rnorm(S)                        # one new person
y_bar <- mu0 + sig * rowMeans(matrix(rnorm(S * m), S, m))  # mean of m people

# Plug-in: ignores uncertainty in the estimates
plug_one <- pf$fit + c(-1, 1) * qnorm(0.975) * s
plug_bar <- pf$fit + c(-1, 1) * qnorm(0.975) * s / sqrt(m)

res <- rbind(exact_one = exact_one,
             sim_one   = quantile(y_one, c(.025, .975)),
             plug_one  = plug_one,
             exact_bar = exact_bar,
             sim_bar   = quantile(y_bar, c(.025, .975)),
             plug_bar  = plug_bar)
colnames(res) <- c("lower", "upper")
round(cbind(res, width = res[, 2] - res[, 1]), 3)
round(c(fit = unname(pf$fit), se_fit = pf$se.fit, sigma = s), 3)

import numpy as np
from scipy import stats

rng = np.random.default_rng(2026)
n = 40
x = rng.uniform(0, 10, n)
y = 1 + 0.3 * x + rng.normal(0, 0.8, n)        # e.g. a log-scale outcome
X = np.column_stack([np.ones(n), x])
b, *_ = np.linalg.lstsq(X, y, rcond=None)
df = n - 2
s = np.sqrt(np.sum((y - X @ b) ** 2) / df)
XtXi = np.linalg.inv(X.T @ X)
x0, m = 8, 10                                  # new x value, size of new group
v0 = np.array([1, x0])
fit = v0 @ b
se_fit = s * np.sqrt(v0 @ XtXi @ v0)
tq = stats.t.ppf(0.975, df)
pm = np.array([-1, 1])

# Exact (t-based) 95% intervals
exact_one = fit + pm * tq * np.sqrt(se_fit**2 + s**2)
exact_bar = fit + pm * tq * np.sqrt(se_fit**2 + s**2 / m)

# Simulation: draw sigma, then coefficients given sigma, then new outcomes
S = 20000
sig = np.sqrt(df * s**2 / rng.chisquare(df, S))
z = rng.normal(size=(S, 2)) @ np.linalg.cholesky(XtXi).T
beta = b + z * sig[:, None]
mu0 = beta[:, 0] + beta[:, 1] * x0                 # one mean per draw
y_one = mu0 + sig * rng.normal(size=S)             # one new person
y_bar = mu0 + sig * rng.normal(size=(S, m)).mean(axis=1)  # mean of m people

# Plug-in: ignores uncertainty in the estimates
z975 = stats.norm.ppf(0.975)
plug_one = fit + pm * z975 * s
plug_bar = fit + pm * z975 * s / np.sqrt(m)

rows = {
    "exact_one": exact_one, "sim_one": np.percentile(y_one, [2.5, 97.5]),
    "plug_one": plug_one, "exact_bar": exact_bar,
    "sim_bar": np.percentile(y_bar, [2.5, 97.5]), "plug_bar": plug_bar,
}
for name, (lo, hi) in rows.items():
    print(f"{name:10s} lower {lo:.3f} upper {hi:.3f} width {hi - lo:.3f}")
print(f"fit {fit:.3f} se_fit {se_fit:.3f} sigma {s:.3f}")

          lower upper width
exact_one 1.805 4.922 3.118
sim_one   1.835 4.923 3.088
plug_one  1.919 4.808 2.890
exact_bar 2.711 4.016 1.305
sim_bar   2.704 4.022 1.319
plug_bar  2.907 3.820 0.914
   fit se_fit  sigma 
 3.364  0.223  0.737

The figures above come from one seeded R run. The predicted mean at x = 8 is 3.364 with a standard error of 0.223, and the residual SD is 0.737. The simulated intervals match the exact ones to within simulation noise: [1.835, 4.923] against [1.805, 4.922] for one case, and [2.704, 4.022] against [2.711, 4.016] for the mean of 10. The plug-in interval is only about 7% narrower for one case (width 2.890 against 3.118), but about 30% narrower for the group mean (0.914 against 1.305), so it would cover the future group average far less often than 95% of the time. Python's random numbers differ from R's, so its figures differ slightly, but they show the same pattern.

Practical points

• Draw the SD too. Holding the residual SD at its estimate makes the interval too narrow, mainly in small samples. A normal draw for log SD, the scale survreg() reports it on, is a reasonable approximation unless the sample is small; for a plain linear model the chi-square draw used above is exact.
• Censored survival data. With censoring there is no closed-form interval, so simulation from the fitted model (for example a log-normal model fitted with survreg() in the survival package) is the practical route. The same recipe applies; just use the coefficients, log scale parameter and covariance matrix the fit reports.
• Going back to the original scale. Quantiles survive any increasing transformation, so exponentiating the limits for one case gives a valid interval for the time itself. For a group, the exponentiated average of log times is the geometric mean time; for the ordinary mean, exponentiate each simulated case first, then average.
• Enough draws. Tail percentiles are noisier than the median. With 20,000 draws, the simulated limits in the example land within about 0.03 of the exact ones; rerun with another seed to see how much yours move.
• The model is still an assumption. The interval is only as good as the distributional assumption behind the noise. Check the residuals before trusting the tails.

If you are unsure whether a prediction interval or a confidence interval is what your question needs, the DASS blog covers interval estimates in more depth, and our test chooser can help with the surrounding analysis.

How to report a prediction interval in APA style (7th edition)

APA has no standard abbreviation for a prediction interval, so write it out the first time and say whether it covers one new case or a group average. Using the example above:

"At x = 8, the predicted outcome was 3.364, 95% prediction interval (PI) [1.805, 4.922] for a single new case. For the mean of 10 new cases with the same predictor value, the 95% PI was [2.711, 4.016]."

If you simulated the interval, say so in the method section, for example: "Prediction intervals were obtained from 20,000 simulations that drew the model parameters from their estimated sampling distribution and added residual noise for each new case." State the scale (log or original units) next to the numbers.

Related tools and guides

• APA 7 formatter for regression
• What does 95% confidence mean?
• Normality of residuals vs the outcome variable
• Prediction interval (Wikipedia)
• survival package (CRAN)

More answered questions

• What does it mean to control for a variable in regression?
• When can you drop the intercept from a linear regression?
• Where does a curve level off? Estimating a plateau point

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.