## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(spBPS)

## ----warning=FALSE, message=F-------------------------------------------------
library(foreach)
library(parallel)
library(doParallel)
library(tictoc)
library(MBA)
library(classInt)
library(RColorBrewer)
library(sp)
library(fields)
library(mvnfast)
library(abind)

## ----results=F----------------------------------------------------------------
# dimensions
n <- 1000
u <- 250
p <- 2
q <- 1

# parameters
B <- c(-0.75, 1.85)
tau2 <- 0.25
sigma2 <- 1
delta <- tau2/sigma2
phi <- 4

set.seed(4-8-15-16-23-42)
# generate sintethic data
crd <- matrix(runif((n+u) * 2), ncol = 2)
X_or <- cbind(rep(1, n+u), matrix(runif((p-1)*(n+u)), ncol = (p-1)))
D <- spBPS:::arma_dist(crd)
Rphi <- exp(-phi * D)
W_or <- matrix(0, n+u) + mniw::rmNorm(1, rep(0, n+u), sigma2*Rphi)
Y_or <- X_or %*% B + W_or + mniw::rmNorm(1, rep(0, n+u), diag(delta*sigma2, n+u))

# train data
crd_s <- crd[1:n, ]
X <- X_or[1:n, ]
W <- W_or[1:n, ]
Y <- Y_or[1:n, ]

# prediction data
crd_u <- crd[-(1:n), ]
X_u <- X_or[-(1:n), ]
W_u <- W_or[-(1:n), ]
Y_u <- Y_or[-(1:n), ]

## ----results=F----------------------------------------------------------------
# priors 
priors <- list(mu_B = matrix(0, nrow = p, ncol = q),
               V_r = diag(10, p),
               Psi = diag(1, q),
               nu = 3)

# hyperparameters values
alfa_seq <- c(0.7, 0.8, 0.9)
phi_seq <- c(3, 4, 5)
hyperpar <- list(alpha = alfa_seq, phi = phi_seq)

## ----results=F----------------------------------------------------------------
# subset dimension
subset_size <- 500

# number of posterior draws
R <- 200

# number of computational cores
n_core <- 1

## ----results=F----------------------------------------------------------------
out <- spBPS(data = list(Y = Y, X = X),
      priors = priors,
      coords = crd_s,
      hyperpar = hyperpar,
      subset_size = subset_size,
      combine_method = "bps",
      draws = R,
      newdata = list(X = X_u, coords = crd_u),
      cores = n_core)

## -----------------------------------------------------------------------------
# statistics computations W
pred_mat_W <- do.call(abind, c(lapply(out$predictive, function(x) x$Wu), along = 3))
post_mean_W <- apply(pred_mat_W, c(1,2), mean)
post_qnt_W <- apply(pred_mat_W, c(1,2), quantile, c(0.025, 0.975))

# Empirical coverage for W
coverage_W <- mean(W_u >= post_qnt_W[1,,1] & W_u <= post_qnt_W[2,,1])
cat("Empirical coverage for Spatial process:", round(coverage_W, 3),"\n")

# statistics computations Y
pred_mat_Y <- do.call(abind, c(lapply(out$predictive, function(x) x$Yu), along = 3))
post_mean_Y <- apply(pred_mat_Y, c(1,2), mean)
post_qnt_Y <- apply(pred_mat_Y, c(1,2), quantile, c(0.025, 0.975))

# Empirical coverage for Y
coverage_Y <- mean(Y_u >= post_qnt_Y[1,,1] & Y_u <= post_qnt_Y[2,,1])
cat("Empirical coverage for Response:", round(coverage_Y, 3),"\n")

# Root Mean Square Prediction Error
rmspe_W <- sqrt( mean( (W_u - post_mean_W)^2 ) )
rmspe_Y <- sqrt( mean( (Y_u - post_mean_Y)^2 ) )
cat("RMSPE for Spatial process:", round(rmspe_W, 3), "\n")
cat("RMSPE for Response:", round(rmspe_Y, 3), "\n")

## ----results=F, echo=F--------------------------------------------------------
# True spatial process surface interpolation
h <- 12
surf.W <- MBA::mba.surf(cbind(crd_s, W), no.X = 500, no.Y = 500,
                        exten = TRUE, sp = TRUE, h = h)$xyz.est
surf.brks <- classIntervals(surf.W$z, 100, 'pretty')$brks
col.pal <- colorRampPalette(brewer.pal(11,'RdBu')[11:1])
xlim <- c(0, 1)
zlim <- range(surf.W$z)

# image for plot
iw <- as.image.SpatialGridDataFrame(surf.W)

# BPS surfaces interpolation
h <- 12
surf.Wp <- MBA::mba.surf(cbind(crd_u, post_mean_W), no.X = 500, no.Y = 500,
                         exten = TRUE, sp = TRUE, h = h)$xyz.est
zlimp <- range(surf.Wp$z)

# image for plot
iwp <- as.image.SpatialGridDataFrame(surf.Wp)

## ----echo=F, fig.dim = c(7.25, 4)---------------------------------------------
# Plotting
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))

plot(crd, type="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="Northing", xlab="Easting",
     main="Spatial process")
axis(2, las=1)
axis(1)
image.plot(iw, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)

plot(crd, type="n", cex=0.5, xlim=xlim, axes=F, ylab="Northing", xlab="Easting")
title(main="Spatial process (Prediction)")
mtext(side = 3, paste("RMSPE :", round(rmspe_W, 3)))
axis(2, las=1)
axis(1)
image.plot(iwp, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlimp)

par(oldpar)

## ----results=F, echo=F--------------------------------------------------------
# True response surface interpolation
h <- 12
surf.Y <- MBA::mba.surf(cbind(crd_s, Y), no.X = 500, no.Y = 500,
                        exten = TRUE, sp = TRUE, h = h)$xyz.est
surf.brks <- classIntervals(surf.Y$z, 100, 'pretty')$brks
col.pal <- colorRampPalette(brewer.pal(11,'RdBu')[11:1])
xlim <- c(0, 1)
zlim <- range(surf.Y$z)

# image for plot
iy <- as.image.SpatialGridDataFrame(surf.Y)

# BPS surfaces interpolation
h <- 12
surf.Yp <- MBA::mba.surf(cbind(crd_u, post_mean_Y), no.X = 500, no.Y = 500,
                         exten = TRUE, sp = TRUE, h = h)$xyz.est
zlimp <- range(surf.Yp$z)

# image for plot
iyp <- as.image.SpatialGridDataFrame(surf.Yp)

## ----echo=F, fig.dim = c(7.25, 4)---------------------------------------------
# Plotting
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))

plot(crd, type="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="Northing", xlab="Easting",
     main="Response")
axis(2, las=1)
axis(1)
image.plot(iy, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)

plot(crd, type="n", cex=0.5, xlim=xlim, axes=F, ylab="Northing", xlab="Easting")
title(main="Response (Prediction)")
mtext(side = 3, paste("RMSPE :", round(rmspe_Y, 3)))
axis(2, las=1)
axis(1)
image.plot(iyp, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlimp)

par(oldpar)