1# Demo for prediction and confidence intervals 2 3par(ask = TRUE) 4# Classical Gaussian Version 5sub <- "Classical Gaussian Version" 6plot(cars$speed,cars$dist,xlim=c(0,27),ylim=c(-10,130), main = "Ezekiel Data", 7 sub = sub, xlab = "speed (mph)",ylab = "stopping distance (ft)") 8ss = seq(0,27,len=100) 9flm <- lm(dist ~ speed + I(speed^2),data = cars) 10cflm <- predict(flm,list(speed=ss), data = cars, interval = "confidence") 11pflm <- predict(flm,list(speed=ss), data = cars, interval = "prediction") 12matlines(ss, pflm,lty = c(1,2,2),col=c("black","red","red")) 13matlines(ss, cflm,lty = c(1,2,2),col=c("black","blue","blue")) 14 15# Now try several different ways to compute alternative bands based on QR 16fl=rq(dist~speed+I(speed^2),data=cars,tau= 0.25) 17fu=rq(dist~speed+I(speed^2),data=cars,tau= 0.75) 18 19# 1. Covariance matrix estimation using "nid" 20sub <- "Covariance matrix estimation using nid" 21plot(cars$speed,cars$dist,xlim=c(0,27),ylim=c(-10,130), main = "Ezekiel Data", 22 sub = sub, xlab = "speed (mph)",ylab = "stopping distance (ft)") 23pl <- predict(fl,list(speed=ss),interval = "confidence", se = "nid") 24pu <- predict(fu,list(speed=ss),interval = "confidence", se = "nid") 25matlines(ss,pl,lty = c(1,2,0), col=c("black","red","red")) 26matlines(ss,pu,lty = c(1,0,2), col=c("black","red","red")) 27# 2. Covariance matrix estimation using xy-bootstrap 28sub <- "Covariance matrix estimation using xy-bootstrap" 29plot(cars$speed,cars$dist,xlim=c(0,27),ylim=c(-10,130), main = "Ezekiel Data", 30 sub = sub, xlab = "speed (mph)",ylab = "stopping distance (ft)") 31pl <- predict(fl,list(speed=ss),interval = "confidence", se = "boot",bsmethod = "xy") 32pu <- predict(fu,list(speed=ss),interval = "confidence", se = "boot",bsmethod = "xy") 33matlines(ss,pl,lty = c(1,2,0), col=c("black","red","red")) 34matlines(ss,pu,lty = c(1,0,2), col=c("black","red","red")) 35# 3. Covariance matrix estimation using weighted xy-bootstrap 36sub <- "Covariance matrix estimation using weighted xy-bootstrap" 37plot(cars$speed,cars$dist,xlim=c(0,27),ylim=c(-10,130), main = "Ezekiel Data", 38 sub = sub, xlab = "speed (mph)",ylab = "stopping distance (ft)") 39pl <- predict(fl,list(speed=ss),interval = "confidence", se = "boot",bsmethod = "wxy") 40pu <- predict(fu,list(speed=ss),interval = "confidence", se = "boot",bsmethod = "wxy") 41matlines(ss,pl,lty = c(1,2,0), col=c("black","red","red")) 42matlines(ss,pu,lty = c(1,0,2), col=c("black","red","red")) 43# 4. Percentile method using xy-bootstrap 44sub <- "Percentile method using xy-bootstrap" 45plot(cars$speed,cars$dist,xlim=c(0,27),ylim=c(-10,130), main = "Ezekiel Data", 46 sub = sub, xlab = "speed (mph)",ylab = "stopping distance (ft)") 47pl <- predict(fl,list(speed=ss),interval = "confidence", 48 type="percentile", se = "boot", bsmethod = "xy") 49pu <- predict(fu,list(speed=ss),interval = "confidence", 50 type="percentile", se = "boot", bsmethod = "xy") 51matlines(ss,pl,lty = c(1,2,0), col=c("black","red","red")) 52matlines(ss,pu,lty = c(1,0,2), col=c("black","red","red")) 53# 5. Direct method using "iid" covariance matrix 54sub <- "Direct method using iid covariance matrix" 55plot(cars$speed,cars$dist,xlim=c(0,27),ylim=c(-10,130), main = "Ezekiel Data", 56 sub = sub, xlab = "speed (mph)",ylab = "stopping distance (ft)") 57pl <- predict(fl,list(speed=ss),interval = "confidence", type="direct", se = "iid") 58pu <- predict(fu,list(speed=ss),interval = "confidence", type="direct", se = "iid") 59matlines(ss,pl,lty = c(1,2,0), col=c("black","red","red")) 60matlines(ss,pu,lty = c(1,0,2), col=c("black","red","red")) 61# 6. Direct method using "nid" covariance matrix 62sub <- "Direct method using nid covariance matrix" 63plot(cars$speed,cars$dist,xlim=c(0,27),ylim=c(-10,130), main = "Ezekiel Data", 64 sub = sub, xlab = "speed (mph)",ylab = "stopping distance (ft)") 65pl <- predict(fl,list(speed=ss),interval = "confidence", type="direct", se = "nid") 66pu <- predict(fu,list(speed=ss),interval = "confidence", type="direct", se = "nid") 67matlines(ss,pl,lty = c(1,2,0), col=c("black","red","red")) 68matlines(ss,pu,lty = c(1,0,2), col=c("black","red","red")) 69 70par(ask = FALSE) 71