[R] need 95% confidence interval bands on cubic extrapolation

[Marc Schwartz (via MN)]

On Tue, 2005-12-20 at 13:04 -0800, James Salsman wrote:
> Dear R experts:
> 
> I need to get this plot, but also with 95% confidence interval bands:
> 
>    hour <- c(1, 2, 3, 4, 5, 6)
>    millivolts <- c(3.5, 5, 7.5, 13, 40, 58)
> 
>    plot(hour, millivolts, xlim=c(1,10), ylim=c(0,1000))
> 
>    pm <- lm(millivolts ~ poly(hour, 3))
> 
>    curve(predict(pm, data.frame(hour=x)), add=TRUE)
> 
> How can the 95% confidence interval band curves be plotted too?
> 
> Sincerely,
> James Salsman
> 
> P.S.  I know I should be using data frames instead of parallel lists.
> This is just a simple example.

There is an example in ?predict.lm.

Given your data, something like the following will work:

hour <- c(1, 2, 3, 4, 5, 6)
millivolts <- c(3.5, 5, 7.5, 13, 40, 58)

pm <- lm(millivolts ~ poly(hour, 3))

# Now create a new dataset with an interval
# of hours that fits your data above
# This is then used in predict.lm() below
# Smaller increments will create smoother lines in the plot
new <- data.frame(hour = seq(1, 6, 0.5))


# Use the new data and generate confidence intervals
# based upon the model
clim <- predict(pm, new, interval = "confidence")


> clim
         fit        lwr      upr
1   4.400794 -17.659582 26.46117
2   2.879712 -12.954245 18.71367
3   2.817460 -14.317443 19.95236
4   4.252232 -12.822969 21.32743
5   7.222222  -8.051125 22.49557
6  11.765625  -2.374270 25.90552
7  17.920635   2.647288 33.19398
8  25.725446   8.650246 42.80065
9  35.218254  18.083351 52.35316
10 46.437252  30.603295 62.27121
11 59.420635  37.360259 81.48101


# Now use matplot to draw the fitted line (black)
# and the CI's (red)
matplot(new$hour, clim,
        lty = c(1, 2, 2), 
        col = c("black", "red", "red"),
        type = "l", ylab = "predicted y")


See ?predict.lm and ?matplot for more information.

HTH,

Marc Schwartz