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