Say I'm doing a simple multiple regression on the following data (R):

Code:

```
n <- 40 x1 <- rnorm(n, mean=3, sd=1)
x2 <- rnorm(n, mean=4, sd=1.25)
y <- 2*x1 + 3*x2 + rnorm(n, mean=2, sd=1)
mydata <- data.frame(x1, x2, y)
mod <- lm(y ~ x1 + x2, data=mydata)
```

`tmp <- predict(mod) + rnorm(length(predict(mod)), 0, summary(mod)$sigma)`

(R function simulate);`tmp <- rnorm(length(predict(mod)), mean(y), sd(y))`

;