17
$\begingroup$

I have some data that I fitted using a LOESS model in R, giving me this:

enter image description here

The data has one predictor and one response, and it is heteroscedastic.

I also added confidence intervals. The problem is that the intervals are confidence intervals for the line, whereas I am interested in the prediction intervals. For example, the bottom panel is more variable then the top panel, but this is not captured in the intervals.

This question is slightly related: Understanding the confidence band from a polynomial regression, especially the answer by @AndyW, however in his example he uses the relatively straightforward interval="predict" argument that exists in predict.lm, but it is absent from predict.loess.

So I have two very related questions:

  1. How do I get the pointwise prediction intervals for LOESS?
  2. How can I predict values that will capture that interval, i.e. generate a bunch of random numbers that will eventually look somewhat like the original data?

It is possible that I don't need LOESS and should use something else, but I am unfamiliar with my options. Basically it should fit the line using local regression or multiple linear regression, giving me error estimates for the lines, and in addition also different variances for different explanatory variables, so I can predict the distribution of the response variable (y) at certain x values.

$\endgroup$
  • $\begingroup$ Is this a pointwise prediction interval? $\endgroup$ – Glen_b Mar 13 '15 at 5:31
  • $\begingroup$ What do you mean by "this"? And I'm not sure if it's pointwise or not. My question 2 is what I'm looking for - unfortunately I am not familiar with the nomenclature. $\endgroup$ – Gimelist Mar 13 '15 at 5:34
  • $\begingroup$ By 'this' I mean "the thing the question is asking about in the title" $\endgroup$ – Glen_b Mar 13 '15 at 5:38
  • $\begingroup$ So I'm not sure - see my previous comment. Basically I'm looking for an interval that will capture the actual variance in the data points, as described in my question. $\endgroup$ – Gimelist Mar 13 '15 at 5:41
  • 1
    $\begingroup$ The spread may be variable (that's why I opted for local regression in the first place). Single predictor. $\endgroup$ – Gimelist Mar 13 '15 at 7:52
12
$\begingroup$

I don't know how to do prediction bands with the original loess function but there is a function loess.sd in the msir package that does just that! Almost verbatim from the msir documentation:

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)

enter image description here

Your second question is a bit trickier since loess.sd doesn't come with a prediction function, but you can hack it together by linearly interpolating the predicted means and SDs you get out of loess.sd (using approx). These can, in turn, be used to simulate data using a normal distribution with the predicted means and SDs:

# 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)

enter image description here

$\endgroup$
  • $\begingroup$ Exactly what I was looking for. When looking at the method he used by seeing the code with loess.sd, it's not too different from what @rnso suggested in a comment to another question of mine. Thanks! $\endgroup$ – Gimelist Apr 9 '15 at 12:30
  • $\begingroup$ Bootstrap to generate the intervals? $\endgroup$ – SmallChess Jan 30 '17 at 14:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.