# \donttest{
data(quakes)
# group stations to ensure min 20 observations per factor level
# and reduce number of levels for speed
quakes$stations <- factor( cut( quakes$stations, breaks = c(0,15,19,23,30,39,132)) )
# Artificially split data to create prediction data set
set.seed(1)
quakes.pred <- quakes[ ss <- sample(1:nrow(quakes), 500), ]
quakes <- quakes[ -ss, ]
quakes.geogam <- geoGAM(response = "mag",
covariates = c("stations", "depth"),
coords = c("lat", "long"),
data = quakes,
max.stop = 20,
cores = 1)
## compute model based bootstrap with 10 repetitions (use at least 100)
quakes.boot <- bootstrap(quakes.geogam,
newdata = quakes.pred,
R = 10, cores = 1)
# plot predictive distribution for site in row 9
hist( as.numeric( quakes.boot[ 9, -c(1:2)] ), col = "grey",
main = paste("Predictive distribution at", paste( quakes.boot[9, 1:2], collapse = "/" )),
xlab = "predicted magnitude")
# compute 95 % prediction interval and add to plot
quant95 <- quantile( as.numeric( quakes.boot[ 9, -c(1:2)] ), probs = c(0.025, 0.975) )
abline(v = quant95[1], lty = "dashed")
abline(v = quant95[2], lty = "dashed")
# }
Run the code above in your browser using DataLab