# Simulation Using Linear Regression

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.

• 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). Commented Mar 21, 2017 at 19:38

$$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.

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 --
}

• Thanks for the great solution. Is your "object" the fitted predicted value? Commented Mar 23, 2017 at 0:52
• Object is an R lm class object. This is the code taken from the R sources for the S3 simulate.lm method.
– Tim
Commented Mar 23, 2017 at 5:35