¿Existe una función o paquete que simule predicciones para un objeto devuelto desde lm ()?

¿Existe una función única, similar a "runif", "rnorm" y similares que produzca predicciones simuladas para un modelo lineal? Puedo codificarlo por mi cuenta, pero el código es feo y asumo que esto es algo que alguien ha hecho antes.

slope = 1.5
intercept = 0
x = as.numeric(1:10)
e = rnorm(10, mean=0, sd = 1)
y = slope * x + intercept + e
fit = lm(y ~ x, data = df)
newX = data.frame(x = as.numeric(11:15))

Lo que me interesa es una función que se parece a la siguiente línea:

sims = rlm(1000, fit, newX)

Esa función devolvería 1000 simulaciones de valores y, basadas en las nuevas variables x.

Respuestas a la pregunta(1)

Su respuesta a la pregunta