In addition to confidence intervals, we can also compute a prediction interval which estimate the variability of new data rather than a statistic
In this example, we consider a single variable and compute the frequency each value was covered.
x <- runif(1000)
# Middle 90% of values
xq0 <- quantile(x, probs=c(.05,.95))
bks <- seq(0,1,by=.01)
hist(x, breaks=bks, border=NA,
main='Prediction Interval', font.main=1)
abline(v=xq0)
paste0('we are 90% confident that the a future data point will be between ',
round(xq0[1],2), ' and ', round(xq0[2],2) )
## [1] "we are 90% confident that the a future data point will be between 0.05 and 0.95"In this example, we consider a range for \(y_{i}(x)\) rather than for \(m(x)\). These intervals also take into account the residuals — the variability of individuals around the mean.
# Bivariate Data from USArrests
xy <- USArrests[,c('Murder','UrbanPop')]
colnames(xy) <- c('y','x')
xy0 <- xy[order(xy$x),]For a nice overview of different types of intervals, see https://www.jstor.org/stable/2685212. For an in-depth view, see “Statistical Intervals: A Guide for Practitioners and Researchers” or “Statistical Tolerance Regions: Theory, Applications, and Computation”. See https://robjhyndman.com/hyndsight/intervals/ for constructing intervals for future observations in a time-series context. See Davison and Hinkley, chapters 5 and 6 (also Efron and Tibshirani, or Wehrens et al.)
# From "Basic Regression"
xy0 <- xy[order(xy$x),]
X0 <- unique(xy0$x)
reg_lo <- loess(y~x, data=xy0, span=.8)
preds_lo <- predict(reg_lo, newdata=data.frame(x=X0))
# Jackknife CI
jack_lo <- sapply(1:nrow(xy), function(i){
xy_i <- xy[-i,]
reg_i <- loess(y~x, dat=xy_i, span=.8)
predict(reg_i, newdata=data.frame(x=X0))
})
boot_regs <- lapply(1:399, function(b){
b_id <- sample( nrow(xy), replace=T)
xy_b <- xy[b_id,]
reg_b <- lm(y~x, dat=xy_b)
})
plot(y~x, pch=16, col=grey(0,.5),
dat=xy0, ylim=c(0, 20))
lines(X0, preds_lo,
col=hcl.colors(3,alpha=.75)[2],
type='o', pch=2)
# Estimate Residuals CI at design points
res_lo <- sapply(1:nrow(xy), function(i){
y_i <- xy[i,'y']
preds_i <- jack_lo[,i]
resids_i <- y_i - preds_i
})
res_cb <- apply(res_lo, 1, quantile,
probs=c(.025,.975), na.rm=T)
# Plot
lines( X0, preds_lo +res_cb[1,],
col=hcl.colors(3,alpha=.75)[2], lt=2)
lines( X0, preds_lo +res_cb[2,],
col=hcl.colors(3,alpha=.75)[2], lty=2)
# Smooth estimates
res_lo <- lapply(1:nrow(xy), function(i){
y_i <- xy[i,'y']
x_i <- xy[i,'x']
preds_i <- jack_lo[,i]
resids_i <- y_i - preds_i
cbind(e=resids_i, x=x_i)
})
res_lo <- as.data.frame(do.call(rbind, res_lo))
res_fun <- function(x0, h, res_lo){
# Assign equal weight to observations within h distance to x0
# 0 weight for all other observations
ki <- dunif(res_lo$x, x0-h, x0+h)
ei <- res_lo[ki!=0,'e']
res_i <- quantile(ei, probs=c(.025,.975), na.rm=T)
}
res_lo2 <- sapply(X0, res_fun, h=15, res_lo=res_lo)
lines( X0, preds_lo + res_lo2[1,],
col=hcl.colors(3,alpha=.75)[2], lty=1, lwd=2)
lines( X0, preds_lo + res_lo2[2,],
col=hcl.colors(3,alpha=.75)[2], lty=1, lwd=2)
# Bootstrap Prediction Interval
boot_resids <- lapply(boot_regs, function(reg_b){
e_b <- resid(reg_b)
x_b <- reg_b$model$x
res_b <- cbind(e_b, x_b)
})
boot_resids <- as.data.frame(do.call(rbind, boot_resids))
# Homoskedastic
ehat <- quantile(boot_resids$e_b, probs=c(.025, .975))
x <- quantile(xy$x,probs=seq(0,1,by=.1))
boot_pi <- coef(reg)[1] + x*coef(reg)['x']
boot_pi <- cbind(boot_pi + ehat[1], boot_pi + ehat[2])
# Plot Bootstrap PI
plot(y~x, dat=xy, pch=16, main='Prediction Intervals',
ylim=c(-5,20), font.main=1)
polygon( c(x, rev(x)), c(boot_pi[,1], rev(boot_pi[,2])),
col=grey(0,.2), border=NA)
# Parametric PI (For Comparison)
#pi <- predict(reg, interval='prediction', newdata=data.frame(x))
#lines( x, pi[,'lwr'], lty=2)
#lines( x, pi[,'upr'], lty=2)Perhaps the most common approach to selecting a bandwidth is to minimize error. Leave-one-out Cross-validation minimizes the average “leave-one-out” mean square prediction errors: \[\begin{eqnarray} CV &=& \min_{\mathbf{H}} \quad \frac{1}{n} \sum_{i=1}^{n} \left[ Y_{i} - \widehat{Y_{[i]}}(\mathbf{X},\mathbf{H}) \right]^2, % \widehat{Y_{[i]}}(\mathbf{X},\mathbf{H}) &=& \sum_{j\neq i} k(\mathbf{X}_{j},\mathbf{X}_{i},\mathbf{H}) \left[ \widehat{\alpha}(\mathbf{X}_{j}) + \widehat{\beta}(\mathbf{X}_{j}) \mathbf{X}_{i} \right] \end{eqnarray}\] where \(\widehat{Y_{[i]}}(\mathbf{X},\mathbf{H})\) is the predicted value at \(\mathbf{X}_{i}\) based on a dataset that excludes \(\mathbf{X}_{i}\).
There are many types of cross-validation . For example, one extension is , which splits \(N\) datapoints into \(k=1...K\) groups, each sized \(B\), and predicts values for the left-out group. adjusts for the degrees of freedom, whereas the function in R uses by default. You can refer to extensions on a case by case basis.
Minimizing out-sample prediction error is perhaps the simplest computational approach to choose bandwidths, and it also addresses an issue that plagues observational studies in the social sciences of explanations without predictions. It is a problem if your model explains everything and predicts nothing, but minimizing prediction error is not necessarily “best”.
##################
# Crossvalidated bandwidth for regression
##################
y <- (CASchools$read + CASchools$math) / 2
xy_mat <- data.frame(y=y, x1=CASchools$income)
library(np)
## Grid Search
BWS <- seq(1,10,length.out=20)
BWS_CV <- sapply(BWS, function(bw){
E_bw <- sapply(1:nrow(xy_mat), function(i){
llls <- npreg(y~x1, data=xy_mat[-i,],
bws=bw, regtype="ll",
ckertype='epanechnikov', bandwidth.compute=F)
pred_i <- predict(llls, newdata=xy_mat[i,])
e <- (pred_i- xy_mat[i,'y'])
return(e)
})
return( mean(E_bw^2) )
})
## Plot MSE
par(mfrow=c(1,2))
plot(BWS, BWS_CV, ylab='CV', pch=16,
xlab='bandwidth (h)',)
## Plot Resulting Predictions
bw <- BWS[which.min(BWS_CV)]
llls <- npreg(y~x1, data=xy_mat,
ckertype='epanechnikov',
bws=bw, regtype="ll")
plot(xy_mat$x, predict(llls), pch=16, col=grey(0,.5),
xlab='X', ylab='Predictions')
abline(a=0,b=1, lty=2)
## Built in algorithmic Optimziation
llls2 <- npreg(y~x1, data=xy_mat, ckertype='epanechnikov', regtype="ll")
points(xy_mat$x, predict(llls2), pch=2, col=rgb(1,0,0,.25))
## Add legend
add_legend <- function(...) {
opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0),
mar=c(0, 0, 0, 0), new=TRUE)
on.exit(par(opar))
plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
legend(...)
}
add_legend('topright',
col=c(grey(0,.5),rgb(1,0,0,.25)),
pch=c(16,2),
bty='n', horiz=T,
legend=c('Grid Search', 'NP-algorithm'))##################
# CV Application
##################
## Smoothly Estimate X & Y Density
y <- sort(xy_mat$y)
fy <- npudens(y, bandwidth.compute=TRUE)
x1 <- sort(xy_mat$x1)
fx <- npudens(x1, bandwidth.compute=TRUE)
## Smoothly Estimate How Y changes with X
llls2 <- npreg(y~x1,data=xy_mat,
ckertype='epanechnikov',
regtype="ll", bandwidth.compute=TRUE)
layout( matrix(c(2,0,1,3), ncol=2, byrow=TRUE),
widths=c(4/5,1/5), heights=c(1/5,4/5))
## Joint Distribution
par(mar=c(4,4,1,1))
plot(y~x1, data=xy_mat,
pch=16, col=grey(0,.25),
xlab="District Income (1000$)",
ylab="Test Score")
lines( sort(xy_mat$x), predict(llls2)[order(xy_mat$x1)],
pch=16, col=1)
## Marginal Distribution
par(mar=c(0,4,1,1))
plot(x1, predict(fx),
col=grey(0,1), type='l', axes=F,
xlab='', ylab='')
rug(x1, col=grey(0,.25))
par(mar=c(4,0,1,1))
plot(predict(fy), y,
col=grey(0,1), type='l', axes=F,
xlab='', ylab='')
rug(y, col=grey(0,.25), side=2)