---
title: "Double Bayesian Predictive Stacking for Spatial Analysis - Tutotial"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Double Bayesian Predictive Stacking for Spatial Analysis - Tutotial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

We provide a brief tutorial of the `spBPS` package. Here we shows the implementation of the Double Bayesian Predictive Stacking on synthetically univariate generated data. In particular, we focus on parallel computing using the packages `parallel`, `doParallel`; but it is not mandatory: it suffices to make it sequential. For any further details please refer to [@presicce2024]. More examples, for multivariate data, are available in documentation, and functions help.

```{r setup}
library(spBPS)
```

### Working packages
```{r, 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)
```


### Data generation
We generate data from the model detailed in Equation (2.4) [@presicce2024], over a unit square.
```{r, 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), ]
```

### Setting priors and hyperparameters
```{r, 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)
```

### Setting dimensions
```{r, results=F}
# subset dimension
subset_size <- 500

# number of posterior draws
R <- 200

# number of computational cores
n_core <- 1
```

### Double BPS parallel fit
Parallel implementation, exploiting 1 computing core.
```{r, 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)
```

### Results collection
```{r}
# 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")
```

### Plot results
```{r, 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)
```

```{r, 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)
```

```{r, 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)
```

```{r, 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)
```
------