library("knitr") ## set global chunk options # opts_chunk$set(fig.path = "figures", fig.align = "center", fig.show = "hold") # options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) opts_chunk$set(prompt = TRUE) options(replace.assign = TRUE, width = 90, prompt = "R> ") ## ----pkgs, echo=TRUE, message=TRUE, warning=TRUE----------- library("spBayes") library("spNNGP") library("fields") library("MBA") library("geoR") library("ggplot2") library("viridis") library("raster") library("sfsmisc") library("extrafont") library("plotROC") library("conflicted") library("coda") library("ggdag") library("dplyr") library("rgdal") library("raster") # font_import(prompt = FALSE) loadfonts() # font_install("fontcm") sys <- local({ I <- Sys.cpuinfo() I[!grepl("^flags", names(I))] }) sys <- sys[!grepl("[.][0-9]+$", names(sys))] conflict_prefer("spDiag", "spNNGP") ## all figures will be held in the figures directory SAVE <- FALSE # save figures if (SAVE) dir.create("figures", showWarnings = FALSE) ### ----simData, echo=TRUE, message=TRUE, warning=TRUE-------- set.seed(1) rmvn <- function(n, mu = 0, V = matrix(1)) { p <- length(mu) if (any(is.na(match(dim(V), p)))) stop("Dimension problem!") D <- chol(V) t(matrix(rnorm(n * p), ncol = p) %*% D + rep(mu, rep(n, p))) } n.mod <- 5000 coords.mod <- cbind(runif(n.mod, 0, 1), runif(n.mod, 0, 1)) coords.ho <- as.matrix(expand.grid(seq(0, 1, length.out = 50), seq(0, 1, length.out = 50))) coords <- rbind(coords.mod, coords.ho) colnames(coords) <- c("x", "y") n <- nrow(coords) beta <- c(1, -0.1) tau.sq <- 0.25 sigma.sq <- 1 phi <- 3 / 0.5 D <- as.matrix(dist(coords)) R <- exp(-phi * D) w <- rmvn(1, rep(0, n), sigma.sq * exp(-phi * D)) x <- rnorm(n) y <- rnorm(n, cbind(1, x) %*% beta + w, sqrt(tau.sq)) ho <- (n.mod + 1):n w.ho <- w[ho] y.ho <- y[ho] x.ho <- x[ho] coords.ho <- coords[ho, ] w <- w[-ho] y <- y[-ho] x <- x[-ho] coords <- coords[-ho, ] ## ----simSeq, cache=TRUE------------------------------------------------------- n.samples <- 2000 starting <- list(phi = 3 / 0.5, sigma.sq = 1, tau.sq = 1) priors <- list(phi.Unif = c(3 / 1, 3 / 0.1), sigma.sq.IG = c(2, 1), tau.sq.IG = c(2, 1)) cov.model <- "exponential" tuning <- list(phi = 0.2) ord <- order(coords[, 1] + coords[, 2]) sim.s <- spNNGP(formula = y ~ x, coords = coords, starting = starting, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.neighbors = 10, method = "latent", ord = ord, n.omp.threads = 4, n.report = 1000, fit.rep = TRUE, sub.sample = list(start = 1000), return.neighbor.info = TRUE) ## ----simRes, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE---- tuning <- list(phi = 0.1, sigma.sq = 0.1, tau.sq = 0.1) sim.r <- spNNGP(formula = y ~ x, coords = coords, starting = starting, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.neighbors = 10, method = "response", ord = ord, n.omp.threads = 4, n.report = 1000, fit.rep = TRUE, sub.sample = list(start = 1000)) ## ----simGP, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE---- set.seed(1) tuning <- list(phi = 0.01, sigma.sq = 0.001, tau.sq = 0.001) sim.gp <- spLM(formula = y ~ x, coords = coords, starting = starting, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples) sim.gp <- spRecover(sim.gp, start = 1000) ## ----simRespSum--------------------------------------------------------------- summary(sim.s) ## ----simConjRMSPE, cache=TRUE------------------------------------------------- theta.alpha <- as.matrix(expand.grid(seq(0.01, 1, length.out = 15), seq(3, 30, length.out = 15))) colnames(theta.alpha) <- c("alpha", "phi") sim.c <- spConjNNGP(y ~ x, coords = as.matrix(coords), cov.model = "exponential", sigma.sq.IG = c(2, 0.5 * var(y)), n.neighbors = 15, ord = ord, theta.alpha = theta.alpha, k.fold = 2, score.rule = "rmspe", fit.rep = TRUE, n.samples = 200, n.omp.threads = 4) ## ----simConjGridHead---------------------------------------------------------- utils::head(sim.c$k.fold.scores, n = 3) ## ----simConjCRPS, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE---- sim.c.crps <- spConjNNGP(y ~ x, coords = as.matrix(coords), cov.model = "exponential", sigma.sq.IG = c(2, 0.5 * var(y)), n.neighbors = 15, theta.alpha = theta.alpha, k.fold = 5, score.rule = "crps", fit.rep = TRUE, n.samples = 200, n.omp.threads = 4) ## ----simConjGrid, echo=TRUE, message=TRUE, warning=TRUE---- k.fold <- as.data.frame(sim.c$k.fold.scores) opt.set <- k.fold[which.min(k.fold[, "rmspe"]), ] pt.size <- 8 text.size <- 15 true.alpha <- tau.sq / sigma.sq true.phi <- phi if (SAVE) png(file = "figures/simConjGridRMSPE.png") ggplot(k.fold, aes(alpha, phi, color = rmspe)) + geom_point(alpha = 1, size = pt.size) + scale_colour_viridis_c(name = "RMSPE", direction = -1, option = "D") + xlab(expression(alpha)) + ylab(expression(phi)) + guides(colour = guide_colorbar(ticks = FALSE)) + theme_bw() + geom_point(aes(x = opt.set[1, "alpha"], y = opt.set[1, "phi"]), color = "black", size = pt.size + 2.5, shape = 1) + geom_point(aes(x = true.alpha, y = true.phi), color = "black", size = pt.size + 2.5, shape = 3) + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() k.fold <- as.data.frame(sim.c.crps$k.fold.scores) opt.set <- k.fold[which.min(k.fold[, "crps"]), ] if (SAVE) png(file = "figures/simConjGridCRPS.png") ggplot(k.fold, aes(alpha, phi, color = rmspe)) + geom_point(alpha = 1, size = pt.size) + scale_colour_viridis_c(name = "CRPS", direction = -1, option = "D") + xlab(expression(alpha)) + ylab(expression(phi)) + guides(colour = guide_colorbar(ticks = FALSE)) + theme_bw() + geom_point(aes(x = opt.set[1, "alpha"], y = opt.set[1, "phi"]), color = "black", size = pt.size + 2.5, shape = 1) + geom_point(aes(x = true.alpha, y = true.phi), color = "black", size = pt.size + 2.5, shape = 3) + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() ## ----simSummary, echo=TRUE, message=TRUE, warning=TRUE----- ## Table 1 values in s.sum, c.sum, r.sum, and gp.sum (see table Latex code at the end of the script). s.sum <- summary(sim.s) c.sum <- summary(mcmc(sim.c$p.beta.theta.samples))$quantiles r.sum <- summary(sim.r) gp.sum <- rbind(summary(sim.gp$p.beta.recover.samples)$quantiles, summary(sim.gp$p.theta.recover.samples)$quantiles) table.1 <- list(beta, sigma.sq, phi, tau.sq, s.sum, c.sum, r.sum, gp.sum) conflict_prefer("spDiag", "spNNGP") ## ----simDiag------------------------------------------------------------------ s.diag <- spDiag(sim.s) r.diag <- spDiag(sim.r) c.diag <- spDiag(sim.c) ## ----simDiagSPLM, echo=TRUE, message=TRUE, warning=TRUE---- ## hack spLM object to work with spNNGP's spDiag conflict_prefer("spDiag", "spNNGP") beta.samples <- sim.gp$p.beta.recover.samples w.samples <- sim.gp$p.w.recover.samples theta.samples <- sim.gp$p.theta.recover.samples y.rep.samples <- apply(rbind(theta.samples[, "tau.sq"], sim.gp$X %*% t(beta.samples) + w.samples), 2, function(x) { rnorm(length(x) - 1, x[2:length(x)], sqrt(x[1])) }) obj <- list(s.indx = 1:nrow(beta.samples), y = sim.gp$Y, X = sim.gp$X, p.beta.samples = beta.samples, p.theta.samples = theta.samples, p.w.samples = w.samples, y.rep.samples = y.rep.samples, type = c("latent", "gaussian")) class(obj) <- c("nngp") gp.diag <- spDiag(obj) ## Table 2 values in s.diag, r.diag, c.diag, and gp.diag (see table Latex code at the end of the script). table.2 <- list(gp.diag, s.diag, r.diag, c.diag) ## ----simPredict, cache=TRUE--------------------------------------------------- s.pred <- predict(sim.s, X.0 = cbind(1, x.ho), coords.0 = coords.ho, sub.sample = list(start = 1000, thin = 10), n.omp.threads = 4, n.report = 1000) ## ----simPredict2, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE---- r.pred <- predict(sim.r, X.0 = cbind(1, x.ho), coords.0 = coords.ho, sub.sample = list(start = 1000, thin = 10), n.omp.threads = 4, n.report = 1000, verbose = FALSE) c.pred <- predict(sim.c, X.0 = cbind(1, x.ho), coords.0 = coords.ho, n.omp.threads = 4, n.report = 1000, verbose = FALSE) ## ----simPredPlot, echo=TRUE, message=TRUE, warning=TRUE---- sim.w.pred <- rbind(data.frame(x = coords.ho[, 1], y = coords.ho[, 2], class = rep("Observed ", nrow(coords.ho)), w = w.ho), data.frame(x = coords.ho[, 1], y = coords.ho[, 2], class = rep("Latent predicted", nrow(coords.ho)), w = apply(s.pred$p.w.0, 1, median))) if (SAVE) png(file = "figures/simWHo.png") ggplot(data = sim.w.pred, aes(y = y, x = x, fill = w)) + geom_raster() + facet_grid(. ~ class) + coord_equal() + xlab("Easting") + ylab("Northing") + scale_fill_viridis_c(name = "", option = "C", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() sim.y.pred <- rbind(data.frame(x = coords.ho[, 1], y = coords.ho[, 2], class = rep("Observed", nrow(coords.ho)), outcome = y.ho), data.frame(x = coords.ho[, 1], y = coords.ho[, 2], class = rep("Latent predicted", nrow(coords.ho)), outcome = apply(s.pred$p.y.0, 1, median)), data.frame(x = coords.ho[, 1], y = coords.ho[, 2], class = rep("Response predicted", nrow(coords.ho)), outcome = apply(r.pred$p.y.0, 1, median)), data.frame(x = coords.ho[, 1], y = coords.ho[, 2], class = rep("Conjugate predicted", nrow(coords.ho)), outcome = c.pred$y.0.hat)) dim(sim.y.pred) if (SAVE) png(file = "figures/simYHo.png") ggplot(data = sim.y.pred, aes(y = y, x = x, fill = outcome)) + geom_raster() + facet_wrap(. ~ class) + coord_equal() + xlab("Easting") + ylab("Northing") + scale_fill_viridis_c(name = "", option = "C", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() ## ----nCPU, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- n <- 1e+05 y <- rnorm(n) coords <- cbind(runif(n, 0, 1), runif(n, 0, 1)) n.cpu <- 1:64 wall.time <- matrix(0, nrow = length(n.cpu), ncol = 2) r.nn <- spNNGP(y ~ 1, coords = coords, starting = starting, method = "response", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = 1, n.omp.threads = 4, return.neighbor.info = TRUE, verbose = FALSE) s.nn <- spNNGP(y ~ 1, coords = coords, starting = starting, method = "latent", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = 1, n.omp.threads = 4, return.neighbor.info = TRUE, verbose = FALSE) n.samples <- 25 for (i in n.cpu) { r.time <- spNNGP(y ~ 1, coords = coords, starting = starting, method = "response", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.omp.threads = i, neighbor.info = r.nn$neighbor.info, verbose = FALSE) s.time <- spNNGP(y ~ 1, coords = coords, starting = starting, method = "latent", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.omp.threads = i, neighbor.info = s.nn$neighbor.info, verbose = FALSE) wall.time[i, 1] <- r.time$run.time[3] wall.time[i, 2] <- s.time$run.time[3] } run.time <- data.frame(time = as.vector(wall.time * 1000 / n.samples / 60), cpu = rep(n.cpu, 2), model = c(rep("Response", nrow(wall.time)), rep("Latent", nrow(wall.time)))) ## ----cpuFig, echo=TRUE, message=TRUE, warning=TRUE--------- if (SAVE) png(file = "figures/cpu_runtime.png") ggplot(run.time, aes(x = cpu, y = time, color = model)) + geom_line() + theme_bw() + xlab("Number of cores") + ylab("Time (minutes per 1000 iterations)") + scale_y_continuous(breaks = round(seq(0, 20, by = 2), 1)) + theme(text = element_text(family = "CM Roman", size = text.size)) + labs(color = "Model") if (SAVE) dev.off() ## ----n, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- set.seed(1) n <- c(1000, 10000, 20000, 5e+05, 1e+05, seq(2e+05, 5e+06, by = 2e+05)) tuning <- list(phi = 0.01, sigma.sq = 0.01, tau.sq = 0.01) n.cpu <- 40 wall.time <- matrix(0, nrow = length(n), ncol = 2) search.time <- matrix(0, nrow = length(n), ncol = 2) n.samples <- 10 for (i in 1:length(n)) { y <- rnorm(n[i]) coords <- cbind(runif(n[i], 0, 1), runif(n[i], 0, 1)) r.time <- spNNGP(y ~ 1, coords = coords, starting = starting, method = "response", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.omp.threads = n.cpu, return.neighbor.info = TRUE, verbose = FALSE) s.time <- spNNGP(y ~ 1, coords = coords, starting = starting, method = "latent", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.omp.threads = n.cpu, return.neighbor.info = TRUE, verbose = FALSE) wall.time[i, 1] <- r.time$run.time[3] wall.time[i, 2] <- s.time$run.time[3] search.time[i, 1] <- r.time$neighbor.info$nn.indx.run.time[3] search.time[i, 2] <- s.time$neighbor.info$nn.indx.run.time[3] + s.time$neighbor.info$u.indx.run.time[3] } run.time <- data.frame(time = as.vector(wall.time * 1000 / n.samples / 60), n = rep(n, 2), model = c(rep("Response", nrow(wall.time)), rep("Latent", nrow(wall.time)))) ## ----nFig, echo=TRUE, message=TRUE, warning=TRUE----------- if (SAVE) png(file = "figures/n_runtime.png") ggplot(run.time, aes(x = n, y = time, color = model)) + geom_line() + theme_bw() + xlab("Number of observations (n)") + ylab("Time (minutes per 1000 iterations)") + theme(text = element_text(family = "CM Roman", size = text.size)) + labs(color = "Model") if (SAVE) dev.off() ## ----bcefLoadPTC, echo=TRUE-------------------------------------------------- load("data/BCEF_PTC.RData") ## ----bcefMap, echo=TRUE-------------------------------------- data(BCEF, package = "spNNGP") set.seed(1) mod <- sample(1:nrow(BCEF), 1e+05) BCEF.mod <- BCEF[mod, ] BCEF.ho <- BCEF[-mod, ] ## subset just for plotting BCEF.mod.plt <- BCEF.mod[seq(1, nrow(BCEF.mod), by = 10), ] BCEF.ho.plt <- BCEF.ho[seq(1, nrow(BCEF.ho), by = 50), ] pts <- rbind(BCEF.mod.plt[, c("x", "y")], BCEF.ho.plt[, c("x", "y")]) pts$type <- as.factor(c(rep("FCH (training)", nrow(BCEF.mod.plt)), rep("FCH (validation)", nrow(BCEF.ho.plt)))) if (SAVE) png(file = "figures/BCEF.png") ggplot(data = BCEF_PTC, aes(y = y, x = x)) + geom_raster(aes(fill = PTC)) + geom_point(data = pts, aes(x = x, y = y, color = type, size = type)) + coord_equal() + xlab("Easting (km)") + ylab("Northing (km)") + scale_color_manual(values = c("#C6417D", "green")) + scale_size_manual(values = c(0.1, 0.1)) + labs(color = "", size = "") + scale_fill_viridis_c(name = "Percent\ncanopy", option = "D", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) + guides(col = guide_legend(override.aes = list(shape = 16, size = 5))) if (SAVE) dev.off() ## ----bcefVario, echo=TRUE------------------------------------ set.seed(1) sub <- sample(1:nrow(BCEF), 25000) m <- lm(BCEF$FCH[sub] ~ BCEF$PTC[sub]) v <- variog(coords = BCEF[sub, c("x", "y")], data = resid(m), uvec = (seq(0, 10, length = 25))) v.fit <- variofit(v, ini.cov.pars = c(50, 4), cov.model = "exponential", minimisation.function = "nls", weights = "equal") cex.lab <- 1.5 cex.axis <- 1.5 if (SAVE) pdf(file = "figures/BCEF-vario.pdf", height = 6, width = 6) plot(1, axes = F, typ = "n", xlim = c(0, 10), ylim = c(0, 55), cex.lab = cex.lab, cex.axis = cex.axis, xlab = "Distance (km)", ylab = "") mtext("Semivariance", side = 2, line = 2.5, cex = cex.lab) axis(1, seq(0, 10, 2), cex.lab = cex.lab, cex.axis = cex.axis) axis(2, seq(0, 55, 10), cex.lab = cex.lab, cex.axis = cex.axis) abline(h = v.fit$nugget, col = "#020212", lw = 2) abline(h = v.fit$cov.pars[1] + v.fit$nugget, col = "#020212", lw = 2) abline(v = 3 / (1 / v.fit$cov.pars[2]), col = "#020212", lw = 2) points(v$u, v$v, pch = 19, cex = 2, col = "#297B8E") lines(v.fit, col = "#020212", lw = 2) if (SAVE) dev.off() ## ----bcefSeq, cache=TRUE, cache.lazy=FALSE------------------------------------ n.samples <- 5000 starting <- list(phi = 3 / 2, sigma.sq = 40, tau.sq = 1) priors <- list(phi.Unif = c(3 / 10, 3 / 0.1), sigma.sq.IG = c(2, 40), tau.sq.IG = c(2, 10)) cov.model <- "exponential" tuning <- list(phi = 0.02) bcef.s <- spNNGP(FCH ~ PTC, coords = c("x", "y"), data = BCEF.mod, starting = starting, method = "latent", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.omp.threads = 40, n.report = 2500, fit.rep = TRUE, sub.sample = list(start = 4000, thin = 10)) ## ----bcefRes, cache=TRUE, cache.lazy=FALSE------------------------------------ tuning <- list(phi = 0.01, sigma.sq = 0.01, tau.sq = 0.005) bcef.r <- spNNGP(FCH ~ PTC, coords = c("x", "y"), data = BCEF.mod, starting = starting, method = "response", n.neighbors = 10, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.omp.threads = 40, n.report = 2500, fit.rep = TRUE, sub.sample = list(start = 4000, thin = 10), verbose = FALSE) ## ----bcefLM, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE---- bcef.lm <- bayesLMRef(lm(FCH ~ PTC, data = BCEF.mod), n.samples = 100) ## hack bayesLMRef object to work with spNNGP's spDiag obj <- list(s.indx = 1:nrow(bcef.lm$p.beta.tauSq.samples), y = bcef.lm$Y, X = bcef.lm$X, p.beta.samples = bcef.lm$p.beta.tauSq.samples[, 1:ncol(bcef.lm$X)], p.theta.samples = bcef.lm$p.beta.tauSq.samples, p.w.samples = matrix(0, nrow(BCEF.mod), nrow(bcef.lm$p.beta.tauSq.samples)), y.rep.samples = apply(bcef.lm$p.beta.tauSq.samples, 1, function(p) { rnorm(nrow(BCEF.mod), bcef.lm$X %*% p[1:ncol(bcef.lm$X)], sqrt(p["tau.sq"])) }), type = c("latent", "gaussian")) class(obj) <- c("nngp") lm.diag <- spNNGP::spDiag(obj) lm.sum <- summary(mcmc(bcef.lm$p.beta.tauSq.samples))$quantiles ## ----bcefConj, cache=TRUE----------------------------------------------------- theta.alpha <- as.matrix(expand.grid(seq(0.1, 1, length.out = 15), seq(3 / 10, 3 / 0.1, length.out = 15))) colnames(theta.alpha) <- c("alpha", "phi") bcef.c <- spConjNNGP(FCH ~ PTC, coords = c("x", "y"), data = BCEF.mod, cov.model = "exponential", sigma.sq.IG = c(2, 40), n.neighbors = 10, theta.alpha = theta.alpha, k.fold = 2, score.rule = "crps", fit.rep = TRUE, n.samples = 200, n.omp.threads = 40, verbose = FALSE) ## ----bcefConjGrid, echo=TRUE, message=TRUE, warning=TRUE---- k.fold <- as.data.frame(bcef.c$k.fold.scores) opt.set <- k.fold[which.min(k.fold[, "crps"]), ] if (SAVE) png(file = "figures/bcefConjGridRMSPE.png") ggplot(k.fold, aes(alpha, phi, color = rmspe)) + geom_point(alpha = 1, size = pt.size) + scale_colour_viridis_c(name = "CRPS", direction = -1, option = "D") + xlab(expression(alpha)) + ylab(expression(phi)) + guides(colour = guide_colorbar(ticks = FALSE)) + theme_bw() + geom_point(aes(x = opt.set[1, "alpha"], y = opt.set[1, "phi"]), color = "black", size = pt.size + 2.5, shape = 1) + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() ## ----bcefSummary, echo=TRUE, message=TRUE, warning=TRUE---- s.sum <- summary(bcef.s) c.sum <- summary(mcmc(bcef.c$p.beta.theta.samples))$quantiles r.sum <- summary(bcef.r) ## Table 3 values in lm.sum, s.sum, c.sum, and r.sum (see table Latex code at the end of the script). table.3 <- list(lm.sum, s.sum, c.sum, r.sum) ## ----bcefDiag, cache=TRUE, echo=TRUE----------------------------------------- s.diag <- spNNGP::spDiag(bcef.s) r.diag <- spNNGP::spDiag(bcef.r) c.diag <- spNNGP::spDiag(bcef.c) ## ----bcefHOPredict, cache=TRUE, echo=TRUE------------------------------------ s.ho.pred <- predict(bcef.s, X.0 = cbind(1, BCEF.ho$PTC), coords.0 = as.matrix(BCEF.ho[, c("x", "y")]), sub.sample = list(start = floor(0.5 * n.samples), thin = 10), n.omp.threads = 40, verbose = FALSE) r.ho.pred <- predict(bcef.r, X.0 = cbind(1, BCEF.ho$PTC), coords.0 = as.matrix(BCEF.ho[, c("x", "y")]), sub.sample = list(start = floor(0.5 * n.samples), thin = 10), n.omp.threads = 40, verbose = FALSE) c.ho.pred <- predict(bcef.c, X.0 = cbind(1, BCEF.ho$PTC), coords.0 = as.matrix(BCEF.ho[, c("x", "y")]), n.omp.threads = 40, verbose = FALSE) ## ----bcefHOPredictSummary, cache=TRUE, echo=TRUE----------------------------- crps <- function(s, y) { var <- apply(s, 1, var) mu <- apply(s, 1, mean) sd <- sqrt(var) y.std <- (y - mu) / sd CRPS <- -sd * (1 / sqrt(pi) - 2 * dnorm(y.std) - y.std * (2 * pnorm(y.std) - 1)) return(mean(CRPS)) } crps.conj <- function(mu, var, y) { sd <- sqrt(var) y.std <- (y - mu) / sd CRPS <- -sd * (1 / sqrt(pi) - 2 * dnorm(y.std) - y.std * (2 * pnorm(y.std) - 1)) return(mean(CRPS)) } rmspe <- function(s, y) { sqrt(mean((apply(s, 1, mean) - y)^2)) } quant <- function(s) { quantile(s, prob = c(0.025, 0.975)) } cover <- function(s, y) { q <- t(apply(s, 1, quant)) sum((q[, 1] < y & q[, 2] > y)) / length(y) * 100 } width <- function(s, y) { q <- t(apply(s, 1, quant)) mean(q[, 2] - q[, 1]) } lm.ho.y.0 <- apply(bcef.lm$p.beta.tauSq.samples, 1, function(p) { rnorm(nrow(BCEF.ho), cbind(1, BCEF.ho$PTC) %*% p[1:ncol(bcef.lm$X)], sqrt(p["tau.sq"])) }) bcef.ho.crps <- c(crps(lm.ho.y.0, BCEF.ho$FCH), crps(s.ho.pred$p.y.0, BCEF.ho$FCH), crps(r.ho.pred$p.y.0, BCEF.ho$FCH), crps.conj(c.ho.pred$y.0.hat, c.ho.pred$y.0.hat.var, BCEF.ho$FCH)) bcef.ho.rmspe <- c(rmspe(lm.ho.y.0, BCEF.ho$FCH), rmspe(s.ho.pred$p.y.0, BCEF.ho$FCH), rmspe(r.ho.pred$p.y.0, BCEF.ho$FCH), rmspe(c.ho.pred$y.0.hat, BCEF.ho$FCH)) bcef.ho.cover <- c(cover(lm.ho.y.0, BCEF.ho$FCH), cover(s.ho.pred$p.y.0, BCEF.ho$FCH), cover(r.ho.pred$p.y.0, BCEF.ho$FCH)) bcef.ho.width <- c(width(lm.ho.y.0, BCEF.ho$FCH), width(s.ho.pred$p.y.0, BCEF.ho$FCH), width(r.ho.pred$p.y.0, BCEF.ho$FCH)) ## ----bcefGridPredict, cache=TRUE, echo=TRUE---------------------------------- s.grd.pred <- predict(bcef.s, X.0 = cbind(1, BCEF_PTC$PTC), coords.0 = as.matrix(BCEF_PTC[, c("x", "y")]), sub.sample = list(start = floor(0.75 * n.samples), thin = 10), n.omp.threads = 40, verbose = FALSE) r.grd.pred <- predict(bcef.r, X.0 = cbind(1, BCEF_PTC$PTC), coords.0 = as.matrix(BCEF_PTC[, c("x", "y")]), sub.sample = list(start = floor(0.75 * n.samples), thin = 10), n.omp.threads = 40, verbose = FALSE) c.grd.pred <- predict(bcef.c, X.0 = cbind(1, BCEF_PTC$PTC), coords.0 = as.matrix(BCEF_PTC[, c("x", "y")]), n.omp.threads = 40, verbose = FALSE) ## ----bcefPredictSurf, echo=TRUE, echo=TRUE, message=TRUE, warning=TRUE---- bcef.mu.y.pred <- rbind(data.frame(x = BCEF_PTC[, "x"], y = BCEF_PTC[, "y"], class = rep("Latent predicted", nrow(BCEF_PTC)), outcome = apply(s.grd.pred$p.y.0, 1, mean)), data.frame(x = BCEF_PTC[, "x"], y = BCEF_PTC[, "y"], class = rep("Response predicted", nrow(BCEF_PTC)), outcome = apply(r.grd.pred$p.y.0, 1, mean)), data.frame(x = BCEF_PTC[, "x"], y = BCEF_PTC[, "y"], class = rep("Conjugate predicted", nrow(BCEF_PTC)), outcome = c.grd.pred$y.0.hat)) bcef.var.y.pred <- rbind(data.frame(x = BCEF_PTC[, "x"], y = BCEF_PTC[, "y"], class = rep("Latent predicted", nrow(BCEF_PTC)), outcome = apply(s.grd.pred$p.y.0, 1, var)), data.frame(x = BCEF_PTC[, "x"], y = BCEF_PTC[, "y"], class = rep("Response predicted", nrow(BCEF_PTC)), outcome = apply(r.grd.pred$p.y.0, 1, var)), data.frame(x = BCEF_PTC[, "x"], y = BCEF_PTC[, "y"], class = rep("Conjugate predicted", nrow(BCEF_PTC)), outcome = c.grd.pred$y.0.hat.var)) if (SAVE) png(file = "figures/BCEFMuPred.png", width = 1000) ggplot(data = bcef.mu.y.pred, aes(y = y, x = x, fill = outcome)) + geom_raster() + facet_wrap(. ~ class) + xlab("Easting (km)") + ylab("Northing (km)") + scale_fill_viridis_c(name = "", option = "C", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() if (SAVE) png(file = "figures/BCEFVarPred.png", width = 1000) ggplot(data = bcef.var.y.pred, aes(y = y, x = x, fill = outcome)) + geom_raster() + facet_wrap(. ~ class) + xlab("Easting (km)") + ylab("Northing (km)") + scale_fill_viridis_c(name = "", option = "C", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() ## Table 4 values in lm.diag, s.diag, r.diag, c.diag, bcef.ho.crps, bcef.ho.rmspe, ## bcef.ho.cover, and bcef.ho.width (see table Latex code at the end of the script). table.4 <- list(lm.diag, s.diag, r.diag, c.diag, bcef.ho.crps, bcef.ho.rmspe, bcef.ho.cover, bcef.ho.width) ## ----miData, cache=TRUE, echo=TRUE, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- data(MI_TSCA, package = "spNNGP") set.seed(1) mod <- sample(1:nrow(MI_TSCA), 15000) mi.mod <- MI_TSCA[mod, ] mi.ho <- MI_TSCA[-mod, ] n.samples <- 10000 ## ----miNS, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- m.ns <- PGLogit(TSCA ~ MIN + MAX + SUP + WIP + AET + DEF, data = mi.mod, fit.rep = TRUE, sub.sample = list(start = 9000), n.samples = n.samples) ## ----miNSPred, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- m.ns.diag <- spNNGP::spDiag(m.ns) m.ns.pred <- cbind(1, as.matrix(mi.ho[, c("MIN", "MAX", "SUP", "WIP", "AET", "DEF")])) %*% t(m.ns$p.beta.samples[m.ns$s.indx, ]) m.ns.pred.prob <- 1 / (1 + exp(-m.ns.pred)) ## ----miNNGP, cache=TRUE------------------------------------------------------- n.samples <- 10000 starting <- list(phi = 3 / 50, sigma.sq = 10) tuning <- list(phi = 0.05) priors <- list(phi.Unif = c(3 / 300, 3 / 10), sigma.sq.IG = c(2, 10)) cov.model <- "exponential" m.s <- spNNGP(TSCA ~ MIN + MAX + SUP + WIP + AET + DEF, coords = c("long", "lat"), data = mi.mod, family = "binomial", method = "latent", n.neighbors = 10, starting = starting, tuning = tuning, priors = priors, cov.model = cov.model, n.samples = n.samples, n.report = 5000, fit.rep = TRUE, sub.sample = list(start = 9000), n.omp.threads = 40) ## ----miNNGPPred, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- m.s.diag <- spNNGP::spDiag(m.s) X.0 <- as.matrix(cbind(1, mi.ho[, c("MIN", "MAX", "SUP", "WIP", "AET", "DEF")])) coords.0 <- as.matrix(mi.ho[, c("long", "lat")]) m.s.pred <- predict(m.s, X.0 = X.0, coords.0 = coords.0, sub.sample = list(start = 9000), n.omp.threads = 40, verbose = FALSE) m.s.pred.prob <- 1 / (1 + exp(-m.s.pred$p.y.0)) ## ----miGPP, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- ## gpp knots <- kmeans(mi.mod[, c("long", "lat")], centers = 25)$centers fit <- glm(TSCA ~ MIN + MAX + SUP + WIP + AET + DEF, data = mi.mod, family = "binomial") beta.starting <- coefficients(fit) beta.tuning <- t(chol(vcov(fit))) n.samples <- 25000 m.gpp <- spGLM(TSCA ~ MIN + MAX + SUP + WIP + AET + DEF, data = mi.mod, family = "binomial", coords = as.matrix(mi.mod[, c("long", "lat")]), knots = as.matrix(knots), starting = list(beta = beta.starting, phi = 3 / 100, sigma.sq = 10, w = 0), tuning = list(beta = beta.tuning, phi = 0.1, sigma.sq = 0.01, w = 0.01), priors = list(phi.Unif = c(3 / 300, 3 / 10), sigma.sq.IG = c(2, 10)), n.samples = n.samples, cov.model = "exponential", verbose = TRUE, n.report = 1000) s.indx <- floor(0.9 * n.samples):n.samples beta.samples <- m.gpp$p.beta.theta.samples[s.indx, 1:ncol(m.gpp$X)] theta.samples <- m.gpp$p.beta.theta.samples[s.indx, c("sigma.sq", "phi")] w.samples <- m.gpp$p.w.samples[, s.indx] y.rep.samples <- apply(1 / (1 + exp(-(m.gpp$X %*% t(beta.samples) + w.samples))), 2, function(x) { rbinom(nrow(m.gpp$X), size = rep(1, nrow(m.gpp$X)), prob = x) }) ## hack spGLM object to work with spNNGP's spDiag obj <- list(s.indx = 1:nrow(beta.samples), y = m.gpp$Y, X = m.gpp$X, weights = rep(1, length(m.gpp$Y)), p.beta.samples = beta.samples, p.theta.samples = theta.samples, p.w.samples = w.samples, y.rep.samples = y.rep.samples, type = c("latent", "binomial")) class(obj) <- c("nngp") gpp.diag <- spNNGP::spDiag(obj) gpp.sum <- summary(mcmc(m.gpp$p.beta.theta.samples[s.indx, ]))$quantiles m.ggp.pred.prob <- spBayes::spPredict(m.gpp, pred.covar = X.0, pred.coords = coords.0, start = min(s.indx), verbose = FALSE)$p.y.predictive.samples ## ROC curve mu.s.pred.prob <- apply(m.s.pred.prob, 1, mean) mu.ns.pred.prob <- apply(m.ns.pred.prob, 1, mean) mu.ggp.pred.prob <- apply(m.ggp.pred.prob, 1, mean) df <- data.frame(d = rep(mi.ho$TSCA, 3), m = c(mu.ns.pred.prob, mu.ggp.pred.prob, mu.s.pred.prob), model = c(rep("Non-spatial", nrow(mi.ho)), rep("GPP", nrow(mi.ho)), rep("NNGP", nrow(mi.ho)))) ## ----miROC, echo=TRUE, message=TRUE, warning=TRUE---------- if (SAVE) png(file = "figures/miROC.png", width = 1000) ggplot(df, aes(d = d, m = m, color = model)) + geom_roc(n.cuts = 0, labels = FALSE) + style_roc(theme = theme_grey) + scale_colour_viridis_d(name = "Model", option = "viridis") + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() ## ----miSummary, cache=TRUE, echo=TRUE, message=TRUE, warning=TRUE---- s.sum <- summary(m.s) ns.sum <- summary(m.ns$p.beta.samples)$quantiles ## Table 5 values in ns.sum, gpp.sum, and s.sum (see table Latex code at the end of the script). table.5 <- list(ns.sum, gpp.sum, s.sum) ## Table 6 values in m.ns.diag, gpp.diag, m.s.diag (see table Latex code at the end of the script). table.6 <- list(m.ns.diag, gpp.diag, m.s.diag) ## ----ndviFigs, echo=TRUE, message=TRUE, warning=TRUE, cache.lazy=FALSE---- elev <- raster("data/kruger-elev.tif") ndvi <- raster("data/kruger-ndvi.tif") aoi <- ggplot2::fortify(readOGR("data", "aoi")) aoi[, c("long", "lat")] <- aoi[, c("long", "lat")] / 1000 ser <- cbind(as.data.frame(ndvi, xy = TRUE), as.data.frame(elev)) colnames(ser) <- c("x", "y", "ndvi", "elev") ser[, c("x", "y")] <- ser[, c("x", "y")] / 1000 ser <- ser[!is.na(ser$ndvi), ] ser.mod <- ser[ser[, 3] != 0, ] ser.ho <- ser[ser[, 3] == 0, ] ser[ser$ndvi == 0, "ndvi"] <- NA m <- lm(ser.mod$ndvi ~ ser.mod$elev) var(resid(m)) set.seed(1) sub <- sample(1:nrow(ser.mod), 25000) v <- variog(coords = ser.mod[sub, c("x", "y")], data = resid(m)[sub], uvec = (seq(0, 200, length = 25))) v.fit <- variofit(v, ini.cov.pars = c(0.015, 100), cov.model = "exponential", minimisation.function = "nls", weights = "equal") cex.lab <- 1.5 cex.axis <- 1.5 if (SAVE) pdf(file = "figures/kruger-vario.pdf", height = 6, width = 6) plot(1, axes = F, typ = "n", xlim = c(0, 200), ylim = c(0, 0.015), cex.lab = cex.lab, cex.axis = cex.axis, xlab = "Distance (km)", ylab = "") mtext("Semivariance", side = 2, line = 2.5, cex = cex.lab) axis(1, seq(0, 200, 50), cex.lab = cex.lab, cex.axis = cex.axis) axis(2, seq(0, 0.015, 0.005), cex.lab = cex.lab, cex.axis = cex.axis) abline(h = v.fit$nugget, col = "#020212", lw = 2) abline(h = v.fit$cov.pars[1] + v.fit$nugget, col = "#020212", lw = 2) abline(v = 3 / (1 / v.fit$cov.pars[2]), col = "#020212", lw = 2) points(v$u, v$v, pch = 19, cex = 2, col = "#297B8E") lines(v.fit, col = "#020212", lw = 2) if (SAVE) dev.off() if (SAVE) png(file = "figures/kruger-elev.png") ggplot(data = ser, aes(y = y, x = x)) + geom_raster(aes(fill = elev)) + coord_equal() + xlab("Easting (km)") + ylab("Northing (km)") + scale_fill_viridis_c(name = "Elev (m)", option = "C", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) + guides(col = guide_legend(override.aes = list(shape = 16, size = 10))) if (SAVE) dev.off() if (SAVE) png(file = "figures/kruger-ndvi.png") ggplot(data = ser, aes(y = y, x = x)) + geom_raster(aes(fill = ndvi)) + geom_polygon(data = aoi, aes(x = long, y = lat), alpha = 0, color = "red") + coord_equal() + xlab("Easting (km)") + ylab("Northing (km)") + scale_fill_viridis_c(name = "NDVI", option = "D", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) + guides(col = guide_legend(override.aes = list(shape = 16, size = 10))) if (SAVE) dev.off() if (SAVE) png(file = "figures/kruger-ndvi-aoi-missing.png") ggplot(data = ser, aes(y = y, x = x)) + xlim(min(aoi$long), max(aoi$long)) + ylim(min(aoi$lat), max(aoi$lat)) + geom_raster(aes(fill = ndvi)) + geom_polygon(data = aoi, aes(x = long, y = lat), alpha = 0, color = "red") + coord_equal() + xlab("Easting (km)") + ylab("Northing (km)") + scale_fill_viridis_c(name = "NDVI", option = "D", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) + guides(col = guide_legend(override.aes = list(shape = 16, size = 10))) if (SAVE) dev.off() ## ----set-seed, echo=TRUE----------------------------------------------------- set.seed(1) ## ----ndviConj, cache=TRUE, cache.lazy=FALSE----------------------------------- theta.alpha <- as.matrix(expand.grid(seq(1.5, 0.01, length.out = 10), seq(3 / 200, 3 / 25, length.out = 10))) colnames(theta.alpha) <- c("alpha", "phi") ser.c <- spConjNNGP(ndvi ~ elev, coords = c("x", "y"), data = ser.mod, cov.model = "exponential", sigma.sq.IG = c(2, 0.01), n.neighbors = 10, theta.alpha = theta.alpha, k.fold = 2, score.rule = "crps", X.0 = cbind(1, ser.ho$elev), coords.0 = as.matrix(ser.ho[, c("x", "y")]), n.omp.threads = 40) ## ----ndviConjSummary, cache=TRUE, cache.lazy=FALSE---------------------------- summary(ser.c, digits = 8) ## ----ndviConjGrid, echo=TRUE, message=TRUE, warning=TRUE---- k.fold <- as.data.frame(ser.c$k.fold.scores) opt.set <- k.fold[which.min(k.fold[, "crps"]), ] if (SAVE) png(file = "figures/serConjGridCRPS.png") ggplot(k.fold, aes(alpha, phi, color = rmspe)) + geom_point(alpha = 1, size = pt.size) + scale_colour_viridis_c(name = "CRPS", direction = -1, option = "D") + xlab(expression(alpha)) + ylab(expression(phi)) + guides(colour = guide_colorbar(ticks = FALSE)) + theme_bw() + geom_point(aes(x = opt.set[1, "alpha"], y = opt.set[1, "phi"]), color = "black", size = pt.size + 2.5, shape = 1) + theme(text = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() ## ----ndviConjSave, cache=TRUE, echo=TRUE------------------------------------- save(ser.c, file = "ser.c") ## ----ndviConjFigs, echo=TRUE, message=TRUE, warning=TRUE---- ser.hat <- ser ser.hat[is.na(ser$ndvi), "ndvi"] <- ser.c$y.0.hat if (SAVE) png(file = "figures/kruger-ndvi-aoi-yhat.png") ggplot(data = ser.hat, aes(y = y, x = x)) + xlim(min(aoi$long), max(aoi$long)) + ylim(min(aoi$lat), max(aoi$lat)) + geom_raster(aes(fill = ndvi)) + geom_polygon(data = aoi, aes(x = long, y = lat), alpha = 0, color = "red") + coord_equal() + xlab("Easting (km)") + ylab("Northing (km)") + scale_fill_viridis_c(name = "NDVI", option = "D", direction = -1) + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) + guides(col = guide_legend(override.aes = list(shape = 16, size = 10))) if (SAVE) dev.off() ser.hat[!is.na(ser$ndvi), "ndvi"] <- NA ser.hat[is.na(ser$ndvi), "ndvi"] <- ser.c$y.0.hat.var if (SAVE) png(file = "figures/kruger-ndvi-aoi-yhatvar.png") ggplot(data = ser.hat, aes(y = y, x = x)) + xlim(min(aoi$long), max(aoi$long)) + ylim(min(aoi$lat), max(aoi$lat)) + geom_raster(aes(fill = ndvi)) + geom_polygon(data = aoi, aes(x = long, y = lat), alpha = 0, color = "red") + coord_equal() + xlab("Easting (km)") + ylab("Northing (km)") + scale_fill_viridis_c(name = "NDVI", option = "D", direction = -1, na.value = "transparent") + theme_bw() + theme(text = element_text(family = "CM Roman", size = text.size)) + guides(col = guide_legend(override.aes = list(shape = 16, size = 10))) if (SAVE) dev.off() ## ----OtherFigs, echo=TRUE, message=TRUE, warning=TRUE------ set.seed(79) ##62 Seed selected to make the locations nicely spaced n <- 10 coords <- cbind(runif(n, 0, 1), runif(n, 0, 1)) ## fit some model to get the neighbor list m.r <- spNNGP(rnorm(n) ~ 1, coords = coords, n.neighbors = 3, starting = list(phi = 1, sigma.sq = 1, tau.sq = 1), tuning = list(phi = 0, sigma.sq = 0, tau.sq = 0), priors = list(phi.Unif = c(0.5, 1), sigma.sq.IG = c(2, 5), tau.sq.IG = c(2, 1)), cov.model = "exponential", n.samples = 1, return.neighbor.info = TRUE) coords <- coords[order(coords[, 1]), ] nn.indx <- m.r$neighbor.info$n.indx m <- matrix(0, n, n) for (i in 1:n) { if (!any(is.na(nn.indx[[i]]))) { m[i, nn.indx[[i]]] <- 1 } } df <- as.data.frame(which(m >= 0, arr.ind = TRUE)) value <- as.vector(m) value[value == 0] <- 0 df <- cbind(df, value) colnames(df) <- c("col", "row", "value") df[, 1] <- 11 - df[, 1] df if (SAVE) png(file = "figures/neighbor-mtrix.png") ggplot(df, aes(row, col, fill = value)) + geom_tile(color = "white") + scale_fill_gradient(low = "gray90", high = "black", na.value = "transparent") + guides(fill = "none") + scale_y_continuous(breaks = 1:10, labels = 10:1) + scale_x_continuous(breaks = 1:10, labels = 1:10) + theme_void() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_text(family = "CM Roman", size = text.size, hjust = 1), axis.text.x = element_text(family = "CM Roman", size = text.size, vjust = 155)) if (SAVE) dev.off() a <- tibble(name = double(), x = double(), y = double(), direction = factor(), to = double(), xend = double(), yend = double(), circular = logical()) a <- add_row(a, name = 1, x = coords[1, 1], y = coords[1, 2], direction = NA, to = NA, xend = NA, yend = NA, circular = FALSE) for (i in 1:n) { if (!any(is.na(nn.indx[[i]]))) { for (j in 1:length(nn.indx[[i]])) { a <- add_row(a, name = i, x = coords[i, 1], y = coords[i, 2], direction = "->", to = nn.indx[[i]][j], xend = coords[nn.indx[[i]][j], 1], yend = coords[nn.indx[[i]][j], 2], circular = FALSE) } } } if (SAVE) png(file = "figures/neighbor-graph.png") a %>% ggplot2::ggplot(ggplot2::aes(x = x, y = y, xend = xend, yend = yend)) + geom_dag_edges() + geom_dag_point(size = 16) + geom_dag_text(col = "white", size = 3.88) + theme_bw() + xlab("Easting") + ylab("Northing") + theme(text = element_text(family = "CM Roman", size = text.size), axis.text.y = element_text(family = "CM Roman", size = text.size), axis.text.x = element_text(family = "CM Roman", size = text.size)) if (SAVE) dev.off() ## Code to make tables frmtCI <- function(x) { apply(x, 1, function(xx) paste0(round(as.numeric(xx["50%"]), 2), " (", paste(round(as.numeric(xx[c("2.5%", "97.5%")]), 2), collapse = ", "), ")")) } q <- c("50%", "2.5%", "97.5%") ix <- c(1:3, 5:4) names(table.1) <- c("beta", "sigma.sq", "phi", "tau.sq", "s.sum", "c.sum", "r.sum", "gp.sum") tab1 <- with(table.1, cbind("Coefficients" = c(beta, sigma.sq, phi, tau.sq), "Full GP" = frmtCI(gp.sum[ix, ]), "Latent" = frmtCI(s.sum[ix, ]), "Response" = frmtCI(r.sum[ix, ]), "Conguate" = c(frmtCI(c.sum[1:3, ]), round(sim.c$theta.alpha[1], 2), frmtCI(c.sum[4,,drop = FALSE]))) ) tab1 names(table.2) <- c("gp.diag", "s.diag", "r.diag", "c.diag") tab2 <- with(table.2, cbind("Full GP" = round(unname(do.call("rbind", gp.diag[c("WAIC", "DIC", "GP", "GRS")])), 2), "Latent" = round(unname(do.call("rbind", s.diag[c("WAIC", "DIC", "GP", "GRS")])), 2), "Response" = round(c(rep(NA, 8), unlist(do.call("rbind", r.diag[c("GP", "GRS")]))), 2), "Conjugate" = round(c(rep(NA, 8), unlist(do.call("rbind", c.diag[c("GP", "GRS")]))), 2)) ) tab2 names(table.3) <- c("lm.sum", "s.sum", "c.sum", "r.sum") tab3 <- with(table.3, cbind("Non-spatial" = c(frmtCI(lm.sum[1:2, q]), rep(NA, 2), frmtCI(lm.sum[3, q, drop = FALSE])), "Latent" = c(frmtCI(s.sum[c(1:3, 5:4), q])), "Response" = frmtCI(r.sum[c(1:3, 5:4), q]), "Conjugate" = c(frmtCI(c.sum[1:3, ]), round(as.numeric(sim.c$theta.alpha[1]), 2), frmtCI(c.sum[4,, drop = FALSE]))) ) tab3 names(table.4) <- c("lm.diag", "s.diag", "r.diag", "c.diag", "bcef.ho.crps", "bcef.ho.rmspe", "bcef.ho.cover", "bcef.ho.width") tab4 <- with(table.4, cbind("Non-Spatial" = round(unname(do.call("rbind", lm.diag[c("WAIC", "DIC", "GP", "GRS")])), 2), "Latent" = round(unname(do.call("rbind", s.diag[c("WAIC", "DIC", "GP", "GRS")])), 2), "Response" = round(c(rep(NA, 8), unlist(do.call("rbind", r.diag[c("GP", "GRS")]))), 2), "Conjugate" = round(c(rep(NA, 8), unlist(do.call("rbind", c.diag[c("GP", "GRS")]))), 2)) ) CRPS <- with(table.4, round(bcef.ho.crps, 2)) RMSPE <- with(table.4, round(bcef.ho.rmspe, 2)) CI_Cover <- c(with(table.4, round(bcef.ho.cover, 2)), NA) CI_Width <- c(with(table.4, round(bcef.ho.width, 2)), NA) tab4.more <- t(data.frame(CRPS = CRPS, RMSPE = RMSPE, CI_Cover = CI_Cover, CI_Width = CI_Width)) colnames(tab4.more) <- colnames(tab4) tab4 <- rbind(tab4, tab4.more) tab4 names(table.5) <- c("ns.sum", "gpp.sum", "s.sum") tab5 <- with(table.5, cbind("Non-Spatial" = c(frmtCI(ns.sum[, q]), rep(NA, 2)), "GPP 25 knots" = frmtCI(gpp.sum[, q]), "NNGP Latent" = frmtCI(s.sum[, q])) ) tab5 names(table.6) <- c("m.ns.diag", "gpp.diag", "m.s.diag") tab6 <- with(table.6, cbind("Non-Spatial" = round(unname(do.call("rbind", m.ns.diag[c("WAIC", "DIC")])), 2), "GPP 25 knots" = round(unname(do.call("rbind", gpp.diag[c("WAIC", "DIC")])), 2), "NNGP Latent" = round(unname(do.call("rbind", m.s.diag[c("WAIC", "DIC")])), 2)) ) tab6 # save.image("done.RData")