Simulation Using Linear Regression

Question

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.

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