need 95% confidence interval bands on cubic extrapolation
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