¿Cómo calcular los intervalos de predicción para LOESS?

Tengo algunos datos que ajusté usando un modelo LOESS en R, dándome esto:

Los datos tienen un predictor y una respuesta, y es heteroscedastic.

También agregué intervalos de confianza. El problema es que los intervalos son intervalos de confianza para la línea, mientras que estoy interesado en los intervalos de predicción. Por ejemplo, el panel inferior es más variable que el panel superior, pero esto no se captura en los intervalos.

Esta pregunta está ligeramente relacionada: Comprender la banda de confianza a partir de una regresión polinómica , especialmente la respuesta de @AndyW, sin embargo, en su ejemplo, utiliza el interval="predict"argumento relativamente directo que existe predict.lm, pero está ausente predict.loess.

Entonces tengo dos preguntas muy relacionadas:

1. ¿Cómo obtengo los intervalos de predicción puntuales para LOESS?
2. ¿Cómo puedo predecir valores que capturarán ese intervalo, es decir, generar un montón de números aleatorios que eventualmente se parecerán a los datos originales?

Es posible que no necesite LOESS y deba usar algo más, pero no estoy familiarizado con mis opciones. Básicamente, debe ajustarse a la línea usando regresión local o regresión lineal múltiple, lo que me da estimaciones de error para las líneas, y además también diferentes variaciones para diferentes variables explicativas, por lo que puedo predecir la distribución de la variable de respuesta (y) a ciertos valores de x .

¡No sé cómo hacer bandas de predicción con la loessfunción original , pero hay una función loess.sden el msirpaquete que hace exactamente eso! Casi al pie de la letra de la msirdocumentación:

library(msir)
data(cars)
# Calculates and plots a 1.96 * SD prediction band, that is,
# a 95% prediction band
l <- loess.sd(cars, nsigma = 1.96)
plot(cars, main = "loess.sd(cars)", col="red", pch=19)
lines(l$x, l$y)
lines(l$x, l$upper, lty=2)
lines(l$x, l$lower, lty=2)

Su segunda pregunta es un poco más complicada ya loess.sdque no viene con una función de predicción, pero puede hackearla interpolando linealmente los medios pronosticados y las SD de las que sale loess.sd(usando approx). Estos, a su vez, se pueden utilizar para simular datos utilizando una distribución normal con los medios y las DE pronosticados:

# Simulate x data uniformly and y data acording to the loess fit
sim_x <- runif(100, min(cars[,1]), max(cars[,1]))
pred_mean <- approx(l$x, l$y, xout = sim_x)$y
pred_sd <- approx(l$x, l$sd, xout = sim_x)$y
sim_y <- rnorm(100, pred_mean, pred_sd) 

# Plots 95% prediction bands with simulated data 
plot(cars, main = "loess.sd(cars)", col="red", pch=19)
points(sim_x, sim_y, col="blue")
lines(l$x, l$y)
lines(l$x, l$upper, lty=2)
lines(l$x, l$lower, lty=2)