I am attempting to better understand how to use simulation techniques when working with linear regression models. For context, if I have a model that, say, predicts the margin of victory between two opponents, how can I simulate the contest. While the regression model using predict() in R give me the fitted value, I understand there is still variance in this result -- a prediction interval.

How can I simulate the actual result, using the model, to find how often a player scores more than his opponent?

My initial thinking: finding the prediction interval based on the model. Then, treating the interval as a uniform distribution between the lower and upper bound and then selecting random values within that distribution as the final outcomes. This feels wrong. Help is much appreciated and if you do know R any code would be further appreciated.

  • 4
    $\begingroup$ You're close. The uniform distribution is wrong. You want to take draws from the normal distribution with a mean of the fitted value and the variance of the residual variance. Note that you're assuming lots with OLS (homoskedasticity, etc). $\endgroup$ Commented Mar 21, 2017 at 19:38

1 Answer 1


Start with writing down the model

$$ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i \\ \varepsilon_i \sim \mathcal{N}(0, \sigma^2) $$

what is equivalent to

$$ \mu_i = \beta_0 + \beta_1 x_i \\ y_i \sim \mathcal{N}(\mu_i, \sigma^2) $$

Regression model predicts conditional mean of $Y$ and the "noise" around the mean is normally distributed. So the correct simulation would take the predictions from the model and add normally distributed noise to them, where $\sigma^2$ is the residual variance.

This is exactly what simulate.lm method does:

simulate.lm <- function(object, nsim = 1, seed = NULL, ...)
# -- snip --
    ftd <- fitted(object)   # == napredict(*, object$fitted)
# -- snip --
    n <- length(ftd)
    ntot <- n * nsim
# -- snip --
    vars <- deviance(object)/ df.residual(object)
# -- snip --
    ftd + rnorm(ntot, sd = sqrt(vars))
# -- snip --
  • $\begingroup$ Thanks for the great solution. Is your "object" the fitted predicted value? $\endgroup$
    – a.powell
    Commented Mar 23, 2017 at 0:52
  • $\begingroup$ Object is an R lm class object. This is the code taken from the R sources for the S3 simulate.lm method. $\endgroup$
    – Tim
    Commented Mar 23, 2017 at 5:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.