Tutorial for spVarBayes

Overview

spVarBayes provides scalable Bayesian inference for spatial data using Variational Bayes (VB) and Nearest Neighbor Gaussian Processes (NNGP). All methods are designed to work efficiently even with 100,000 spatial locations, offering a practical alternative to traditional MCMC. It includes:

Mean Field Approximation (MFA) and its Linear Response (LR) correction
NNGP variational distributions for spatial random effects
Joint variational distributions for regression coefficients and spatial effects

Installation

The package is available on the github. To install this package in your R, you can use the code: devtools::install_github("jfsong100/spVarBayes"). Once installed, you can

library(spVarBayes)

This tutorial walks through generating data, fitting each method, and making predictions.

Simulate Data

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)))
}

set.seed(1010)
n <- 1500

coords <- cbind(runif(n,0,5), runif(n,0,5))

# Remove close points
remove_close_points <- function(x, y, threshold) {

  points <- cbind(x, y)
  
  repeat {
    neighbors <- nn2(data = points, k = 2, searchtype = "radius", radius = threshold)
    distances <- neighbors$nn.dists[, 2]
    if (all(is.na(distances) | distances >= threshold)) {
      break 
    }
    closest_index <- which.min(distances)
    points <- points[-closest_index, ]
  }
  
  list(x = points[, 1], y = points[, 2])
}

remove_coords <- remove_close_points(coords[,1], coords[,2], 0.015)
cleaned_x <- remove_coords$x
cleaned_y <- remove_coords$y

coords <- cbind(cleaned_x,cleaned_y)

n = nrow(coords)

x <- cbind(rnorm(n), rnorm(n))
B <- as.matrix(c(1,5))

sigma2_true <- 5
tau2_true <- 1
phi_true <- 6

D <- as.matrix(dist(coords))
R <- exp(-phi_true*D)
w <- rmvn(1, rep(0,n), sigma2_true*R)
y <- rnorm(n, x%*%B + w, sqrt(tau2_true))

# Split into training set and test set
n_train <- 1000
train_index <- sample(1:n, n_train)
y_train <- y[train_index]
x_train <- x[train_index,]
w_train <- w[train_index,]
coords_train <- coords[train_index,]

y_test <- y[-train_index]
x_test <- x[-train_index,]
w_test <- w[-train_index,]
coords_test <- coords[-train_index,]

Fit NNGP

NNGP <- spVB_NNGP(y = y_train,X = x_train,coords=coords_train, n.neighbors = 15, 
                       n.neighbors.vi = 3,
                       rho = 0.85, max_iter = 1500, covariates = TRUE)
#> Warning in spVB_NNGP(y = y_train, X = x_train, coords = coords_train,
#> n.neighbors = 15, : We recommend using higher n.neighbors especially for small
#> Phi.
#> ---------------------------------------- 
#>           Model description 
#> ---------------------------------------- 
#> Using NNGP Gaussian family for Variational Approximation. 
#> 
#> Using 15 nearest neighbors for the prior. 
#> 
#> Using 3 nearest neighbors for the variational family. 
#> 
#> Data is ordered using Sum_coords ordering. 
#> 
#> Using BRISC estimation for diagonals elements in the variational distribuiton for spatial random effects. 
#> 
#> BRISC estimation for var is 0.927327 . 
#> 
#> Using BRISC estimation for phi as initial value. 
#> 
#> BRISC estimation for phi is 4.524489 . 
#> 
#> Using BRISC estimation for tau.sq as initial value. 
#> 
#> BRISC estimation for tau.sq is 1.116363 . 
#> 
#> Using BRISC estimation for sigma.sq as initial value. 
#> 
#> BRISC estimation for sigma.sq is 5.476393 . 
#> 
#> ---------------------------------------- 
#>           Model fitting 
#> ---------------------------------------- 
#> spVB-NNGP model fit with covariates X and using full batch 
#> 
#> Assume block independence for coefficients beta and w 
#> [*---------] 10%
#> [**--------] 20%
#> [***-------] 30%
#> [****------] 40%
#> [*****-----] 50%
#> [******----] 60%
#> [*******---] 70%
#> [********--] 80%
#> [*********-] 90%
#> [**********] 100%
#> Warning in spVB_NNGP(y = y_train, X = x_train, coords = coords_train,
#> n.neighbors = 15, : The output data (y, X, coords, ...) is reordered, use ord
#> to get the original ordered data.

w_var_NNGP <- spVB_get_Vw(NNGP)

Sampling and prediction

NNGP_w_samples <- spVB_w_sampling(NNGP, n.samples = 5000)$p.w.samples
#> Sampling from NNGP variational distribution.

NNGP_predict <- predict(NNGP, coords.0 = coords_test, X.0 = x_test, covariates = TRUE, n.samples = 5000)
#> Sampling from NNGP variational distribution. 
#> ----------------------------------------
#>       Prediction description
#> ----------------------------------------
#> Model fit with 1000 observations.
#> 
#> Using the exponential spatial correlation model.
#> 
#> Using 15 nearest neighbors.
#> 
#> Predicting at 470 locations.
#> 
#> 
#> 
#> Source not compiled with OpenMP support.
#> ----------------------------------------
#>               Predicting
#> ----------------------------------------
#> Location: 100 of 470, 21.28%
#> Location: 200 of 470, 42.55%
#> Location: 300 of 470, 63.83%
#> Location: 400 of 470, 85.11%
#> Location: 470 of 470, 100.00%

Fit NNGP joint model

NNGP_joint <- spVB_NNGP(y = y_train,X = x_train,coords=coords_train, n.neighbors = 15, 
                        n.neighbors.vi = 3,
                        rho = 0.85, max_iter = 1500, covariates = TRUE, joint = TRUE)
#> Warning in spVB_NNGP(y = y_train, X = x_train, coords = coords_train,
#> n.neighbors = 15, : We recommend using higher n.neighbors especially for small
#> Phi.
#> ---------------------------------------- 
#>           Model description 
#> ---------------------------------------- 
#> Using NNGP Gaussian family for Variational Approximation. 
#> 
#> Using 15 nearest neighbors for the prior. 
#> 
#> Using 3 nearest neighbors for the variational family. 
#> 
#> Data is ordered using Sum_coords ordering. 
#> 
#> Using BRISC estimation for diagonals elements in the variational distribuiton for spatial random effects. 
#> 
#> BRISC estimation for var is 0.927327 . 
#> 
#> Using BRISC estimation for phi as initial value. 
#> 
#> BRISC estimation for phi is 4.524489 . 
#> 
#> Using BRISC estimation for tau.sq as initial value. 
#> 
#> BRISC estimation for tau.sq is 1.116363 . 
#> 
#> Using BRISC estimation for sigma.sq as initial value. 
#> 
#> BRISC estimation for sigma.sq is 5.476393 . 
#> 
#> ---------------------------------------- 
#>           Model fitting 
#> ---------------------------------------- 
#> spVB-NNGP model fit with covariates X and using full batch 
#> 
#> Joint modeling for coefficients beta and w 
#> [*---------] 10%
#> [**--------] 20%
#> [***-------] 30%
#> [****------] 40%
#> [*****-----] 50%
#> [******----] 60%
#> [*******---] 70%
#> [********--] 80%
#> [*********-] 90%
#> [**********] 100%
#> Warning in spVB_NNGP(y = y_train, X = x_train, coords = coords_train,
#> n.neighbors = 15, : The output data (y, X, coords, ...) is reordered, use ord
#> to get the original ordered data.


w_var_NNGP_joint <- spVB_get_Vw(NNGP_joint)

Sampling and prediction

NNGP_joint_w_samples <- spVB_joint_sampling(NNGP_joint, n.samples = 5000)$p.w.samples
#> Joint Sampling from NNGP variational distribution for beta and w.

NNGP_joint_predict <- predict(NNGP_joint, coords.0 = coords_test, X.0 = x_test, covariates = TRUE, n.samples = 5000)
#> Joint Sampling from NNGP variational distribution for beta and w. 
#> ----------------------------------------
#>       Prediction description
#> ----------------------------------------
#> Model fit with 1000 observations.
#> 
#> Using the exponential spatial correlation model.
#> 
#> Using 15 nearest neighbors.
#> 
#> Predicting at 470 locations.
#> 
#> 
#> 
#> Source not compiled with OpenMP support.
#> ----------------------------------------
#>               Predicting
#> ----------------------------------------
#> Location: 100 of 470, 21.28%
#> Location: 200 of 470, 42.55%
#> Location: 300 of 470, 63.83%
#> Location: 400 of 470, 85.11%
#> Location: 470 of 470, 100.00%

Fit MFA

MFA <- spVB_MFA(y = y_train,X = x_train,coords=coords_train, covariates = TRUE, 
                n.neighbors = 15, rho = 0.85, max_iter = 1000, LR = FALSE)
#> Warning in spVB_MFA(y = y_train, X = x_train, coords = coords_train, covariates
#> = TRUE, : We recommend using higher n.neighbors especially for small Phi.
#> ---------------------------------------- 
#>           Model description 
#> ---------------------------------------- 
#> Using Mean-field Approximation family for Variational Approximation 
#> 
#> Data is ordered using Sum_coords ordering. 
#> 
#> Using BRISC estimation for diagonals elements in the variational distribuiton for spatial random effects. 
#> 
#> BRISC estimation for var is 0.927327 . 
#> 
#> Using BRISC estimation for phi as initial value. 
#> 
#> BRISC estimation for phi is 4.524489 . 
#> 
#> Using BRISC estimation for tau.sq as initial value. 
#> 
#> BRISC estimation for tau.sq is 1.116363 . 
#> 
#> Using BRISC estimation for sigma.sq as initial value. 
#> 
#> BRISC estimation for sigma.sq is 5.476393 . 
#> 
#> ---------------------------------------- 
#>           Model fitting 
#> ---------------------------------------- 
#> spVB-MFA model fit with covariates X and using full batch 
#> 
#> [*---------] 10%
#> [**--------] 20%
#> [***-------] 30%
#> [****------] 40%
#> [*****-----] 50%
#> [******----] 60%
#> [*******---] 70%
#> [********--] 80%
#> [*********-] 90%
#> [**********] 100%
#> Warning in spVB_MFA(y = y_train, X = x_train, coords = coords_train, covariates
#> = TRUE, : The output data (y, X, coords, ...) is reordered, use ord to get the
#> original ordered data.

Sampling and prediction

MFA_w_samples <- spVB_w_sampling(MFA, n.samples = 5000)$p.w.samples
#> Sampling from mean-field variational distribution.

MFA_predict <- predict(MFA, coords.0 = coords_test, X.0 = x_test, covariates = TRUE, n.samples = 5000)
#> Sampling from mean-field variational distribution. 
#> ----------------------------------------
#>       Prediction description
#> ----------------------------------------
#> Model fit with 1000 observations.
#> 
#> Using the exponential spatial correlation model.
#> 
#> Using 15 nearest neighbors.
#> 
#> Predicting at 470 locations.
#> 
#> 
#> 
#> Source not compiled with OpenMP support.
#> ----------------------------------------
#>               Predicting
#> ----------------------------------------
#> Location: 100 of 470, 21.28%
#> Location: 200 of 470, 42.55%
#> Location: 300 of 470, 63.83%
#> Location: 400 of 470, 85.11%
#> Location: 470 of 470, 100.00%

Fit MFA with linear response

MFA_pre_fit <- spVB_MFA(y = y_train,X = x_train,coords=coords_train, covariates = TRUE, 
                   n.neighbors = 15, rho = 0.85, max_iter = 1000, LR = TRUE)
#> Warning in spVB_MFA(y = y_train, X = x_train, coords = coords_train, covariates
#> = TRUE, : We recommend using higher n.neighbors especially for small Phi.
#> ---------------------------------------- 
#>           Model description 
#> ---------------------------------------- 
#> Using Mean-field Approximation family for Variational Approximation 
#> 
#> Data is ordered using Sum_coords ordering. 
#> 
#> Using BRISC estimation for diagonals elements in the variational distribuiton for spatial random effects. 
#> 
#> BRISC estimation for var is 0.927327 . 
#> 
#> Using BRISC estimation for phi as initial value. 
#> 
#> BRISC estimation for phi is 4.524489 . 
#> 
#> Using BRISC estimation for tau.sq as initial value. 
#> 
#> BRISC estimation for tau.sq is 1.116363 . 
#> 
#> Using BRISC estimation for sigma.sq as initial value. 
#> 
#> BRISC estimation for sigma.sq is 5.476393 . 
#> 
#> ---------------------------------------- 
#>           Model fitting 
#> ---------------------------------------- 
#> spVB-MFA model fit with covariates X and using full batch 
#> 
#> [*---------] 10%
#> [**--------] 20%
#> [***-------] 30%
#> [****------] 40%
#> [*****-----] 50%
#> [******----] 60%
#> [*******---] 70%
#> [********--] 80%
#> [*********-] 90%
#> [**********] 100%
#> Warning in spVB_MFA(y = y_train, X = x_train, coords = coords_train, covariates
#> = TRUE, : The output data (y, X, coords, ...) is reordered, use ord to get the
#> original ordered data.


MFA_LR = spVB_LR(MFA_pre_fit, get_mat = TRUE, get_para = TRUE, n_omp = 5)
#> ------------------------------------------------------- 
#>    Default Number of Threads is 5 
#>    Compute the Linear response for variance 
#> ------------------------------------------------------- 
#>    Update spatial parameters 
#> -------------------------------------------------------

Sampling and prediction

MFA_LR_w_samples <- spVB_LR_sampling(MFA_LR, n.samples = 5000)$p.w.samples
#> Joint Sampling from linear response corrected mean-field variational distribution for beta and w.

MFA_LR_predict <- predict(MFA_LR, coords.0 = coords_test, X.0 = x_test, covariates = TRUE, n.samples = 5000)
#> Joint Sampling from linear response corrected mean-field variational distribution for beta and w. 
#> ----------------------------------------
#>       Prediction description
#> ----------------------------------------
#> Model fit with 1000 observations.
#> 
#> Using the exponential spatial correlation model.
#> 
#> Using 15 nearest neighbors.
#> 
#> Predicting at 470 locations.
#> 
#> 
#> 
#> Source not compiled with OpenMP support.
#> ----------------------------------------
#>               Predicting
#> ----------------------------------------
#> Location: 100 of 470, 21.28%
#> Location: 200 of 470, 42.55%
#> Location: 300 of 470, 63.83%
#> Location: 400 of 470, 85.11%
#> Location: 470 of 470, 100.00%

Compare with spNNGP

n.samples = 5000
starting <- list("phi"= 5, "sigma.sq"=1, "tau.sq"= 0.5)
tuning <- list("phi"=0.15, "sigma.sq"=1.5, "tau.sq"=1.15)
priors <- list("phi.Unif"=c(1/10, 10), "sigma.sq.IG"=c(0.1, 1), "tau.sq.IG"=c(1,1))
cov.model <- "exponential"

library(spNNGP)
intercept = rep(1,n)
m.s <- spNNGP(y_train~x_train-1, coords=coords_train, starting=starting, method="latent",
              n.neighbors=15, priors=priors, tuning=tuning, cov.model=cov.model,
              n.samples=n.samples, n.omp.threads=1, n.report=1000)
#> ----------------------------------------
#>  Building the neighbor list
#> ----------------------------------------
#> ----------------------------------------
#> Building the neighbors of neighbors list
#> ----------------------------------------
#> ----------------------------------------
#>  Model description
#> ----------------------------------------
#> NNGP Latent model fit with 1000 observations.
#> 
#> Number of covariates 2 (including intercept if specified).
#> 
#> Using the exponential spatial correlation model.
#> 
#> Using 15 nearest neighbors.
#> 
#> Number of MCMC samples 5000.
#> 
#> Priors and hyperpriors:
#>  beta flat.
#>  sigma.sq IG hyperpriors shape=0.10000 and scale=1.00000
#>  tau.sq IG hyperpriors shape=1.00000 and scale=1.00000
#>  phi Unif hyperpriors a=0.10000 and b=10.00000
#> 
#> 
#> Source not compiled with OpenMP support.
#> ----------------------------------------
#>      Sampling
#> ----------------------------------------
#> Sampled: 1000 of 5000, 20.00%
#> Report interval Metrop. Acceptance rate: 59.80%
#> Overall Metrop. Acceptance rate: 59.80%
#> -------------------------------------------------
#> Sampled: 2000 of 5000, 40.00%
#> Report interval Metrop. Acceptance rate: 63.70%
#> Overall Metrop. Acceptance rate: 61.75%
#> -------------------------------------------------
#> Sampled: 3000 of 5000, 60.00%
#> Report interval Metrop. Acceptance rate: 62.90%
#> Overall Metrop. Acceptance rate: 62.13%
#> -------------------------------------------------
#> Sampled: 4000 of 5000, 80.00%
#> Report interval Metrop. Acceptance rate: 64.20%
#> Overall Metrop. Acceptance rate: 62.65%
#> -------------------------------------------------
#> Sampled: 5000 of 5000, 100.00%
#> Report interval Metrop. Acceptance rate: 59.00%
#> Overall Metrop. Acceptance rate: 61.92%
#> -------------------------------------------------


burnin = 2000

summary(m.s)
#> 
#> Call:
#> spNNGP(formula = y_train ~ x_train - 1, coords = coords_train, 
#>     method = "latent", n.neighbors = 15, starting = starting, 
#>     tuning = tuning, priors = priors, cov.model = cov.model, 
#>     n.samples = n.samples, n.omp.threads = 1, n.report = 1000)
#> 
#> Model class is NNGP, method latent, family gaussian.
#> 
#> Model object contains 5000 MCMC samples.
#> 
#> Chain sub.sample:
#> start = 2500
#> end = 5000
#> thin = 1
#> samples size = 2501
#>          2.5%   25%    50%    75%    97.5% 
#> x_train1 0.7740 0.8501 0.8898 0.9278 1.0021
#> x_train2 4.9034 4.9813 5.0237 5.0689 5.1459
#> sigma.sq 4.4661 5.1528 5.5364 5.9301 6.6644
#> tau.sq   0.5474 0.8317 1.0437 1.3002 1.6841
#> phi      3.1174 4.2686 4.7528 5.2848 6.2491

what = apply(m.s$p.w.samples[,((burnin+1):n.samples)], 1, mean)
wvarhat = apply(m.s$p.w.samples[,((burnin+1):n.samples)], 1, var)

spNNGP_predict <- predict(m.s, coords.0 = coords_test, X.0 = x_test, n.samples = 5000)
#> ----------------------------------------
#>  Prediction description
#> ----------------------------------------
#> NNGP Latent model fit with 1000 observations.
#> 
#> Number of covariates 2 (including intercept if specified).
#> 
#> Using the exponential spatial correlation model.
#> 
#> Using 15 nearest neighbors.
#> 
#> Number of MCMC samples 5000.
#> 
#> Predicting at 470 locations.
#> 
#> 
#> 
#> Source not compiled with OpenMP support.
#> -------------------------------------------------
#>      Predicting
#> -------------------------------------------------
#> Location: 100 of 470, 21.28%
#> Location: 200 of 470, 42.55%
#> Location: 300 of 470, 63.83%
#> Location: 400 of 470, 85.11%
#> Location: 470 of 470, 100.00%

Compare the posterior mean with true value

df_means <- data.frame(
  true = rep(w_train, 5),
  est = c(
    NNGP$w_mu[order(NNGP$ord)],
    NNGP_joint$w_mu[order(NNGP_joint$ord)],
    MFA$w_mu[order(MFA$ord)],
    MFA_LR$w_mu[order(MFA_LR$ord)],
    what
  ),
  method = rep(c(
    "spVB-NNGP",
    "spVB-NNGP Joint",
    "spVB-MFA",
    "spVB-MFA Linear Response",
    "spNNGP"
  ), each = length(w_train))
)

ggplot(df_means, aes(x = true, y = est)) +
  geom_point(alpha = 0.7, size = 1.6) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed", linewidth = 1) +
  facet_wrap(~ method, ncol = 3) +
  labs(
    x = "True spatial random effect (w)",
    y = "Estimated posterior mean (w)",
    title = "Posterior Mean vs. True Spatial Random Effect Across Methods"
  ) +
  coord_fixed() +
  theme_minimal(base_size = 12)+
  theme(
    plot.title = element_text(size = 11, hjust = 0.5)
  )

Compare the predicted mean with true value

df_test <- data.frame(
  true = rep(w_test, 5),
  est = c(
    apply(NNGP_predict$p.w.0, 1, mean),
    apply(NNGP_joint_predict$p.w.0, 1, mean),
    apply(MFA_predict$p.w.0, 1, mean),
    apply(MFA_LR_predict$p.w.0, 1, mean),
    apply(spNNGP_predict$p.w.0, 1, mean)
  ),
  method = rep(c(
    "spVB-NNGP",
    "spVB-NNGP Joint",
    "spVB-MFA",
    "spVB-MFA Linear Response",
    "spNNGP"
  ), each = length(w_test))
)

ggplot(df_test, aes(x = true, y = est)) +
  geom_point(alpha = 0.7, size = 1.6) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed", linewidth = 1) +
  facet_wrap(~ method, ncol = 3) +
  labs(
    x = "True spatial random effect (w)",
    y = "Predicted posterior mean (w)",
    title = "Predicted vs. True Spatial Random Effect (Test Set)"
  ) +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(size = 11, hjust = 0.5)
  )

Compare the posterior variance with spNNGP

df_all <- data.frame(
  ref = rep(wvarhat, 4),
  est = c(
    diag(w_var_NNGP)[order(NNGP$ord)],
    diag(w_var_NNGP_joint)[-(1:ncol(x_train))][order(NNGP_joint$ord)],
    MFA$w_sigma_sq[order(MFA$ord)],
    diag(MFA_LR$updated_mat)[-(1:ncol(x_train))][order(MFA_LR$ord)]
  ),
  method = rep(c("spVB-NNGP", "spVB-NNGP Joint", "spVB-MFA", "spVB-MFA Linear Response"),
               each = length(wvarhat))
)

ggplot(df_all, aes(x = ref, y = est)) +
  geom_point(alpha = 0.7, size = 1.6) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed", linewidth = 1) +
  facet_wrap(~ method, ncol = 2) +
  labs(
    x = "Estimated variance using spNNGP",
    y = "Estimated variance using method",
    title = "Comparison of Estimated Variance Across Methods"
  ) +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(size = 11, hjust = 0.5)
  )

Compare the estimated mean for spatial parameters

theta_est <- matrix(c(
  MFA$theta,
  MFA_LR$theta,
  NNGP$theta,
  NNGP_joint$theta,
  colMeans(m.s$p.theta.samples[-(1:burnin), ])
), nrow = 3)

rownames(theta_est) <- c("sigmasq", "tausq", "phi")
colnames(theta_est) <- c("spVB-MFA", "spVB-MFA LR", "spVB-NNGP", "spVB-NNGP Joint", "spNNGP")

theta_est
#>          spVB-MFA spVB-MFA LR spVB-NNGP spVB-NNGP Joint   spNNGP
#> sigmasq 6.8422530    5.540780  5.688508        5.639810 5.553913
#> tausq   0.5969584    1.115316  1.032520        1.055593 1.058363
#> phi     4.5326399    4.524489  4.532257        4.532795 4.782708

Compare the estimated mean and variance for regression coefficients

beta_est <- matrix(rbind(cbind(
  MFA$beta,
  MFA_LR$beta,
  NNGP$beta,
  NNGP_joint$beta,
  colMeans(m.s$p.beta.samples[-(1:burnin), ])
),cbind(
  diag(MFA$beta_cov),
  diag(MFA_LR$updated_mat)[1:ncol(x_train)],
  diag(NNGP$beta_cov),
  diag(w_var_NNGP_joint)[1:ncol(x_train)],
  apply(m.s$p.beta.samples[-(1:burnin),], 2, var)
)),nrow = 4)

rownames(beta_est) <- c("beta1.mean", "beta2.mean", "beta1.var", "beta2.var")
colnames(beta_est) <- c("spVB-MFA", "spVB-MFA LR", "spVB-NNGP", "spVB-NNGP Joint", "spNNGP")

beta_est
#>                spVB-MFA spVB-MFA LR   spVB-NNGP spVB-NNGP Joint      spNNGP
#> beta1.mean 0.8870426607 0.888654129 0.888590739     0.888604354 0.888493589
#> beta2.mean 5.0161355638 5.020306200 5.021390517     5.021510614 5.022568298
#> beta1.var  0.0005929757 0.003318319 0.001014088     0.001614907 0.003370545
#> beta2.var  0.0006225177 0.003638638 0.001064610     0.001726749 0.003915777

Tutorial for spVarBayes

Jiafang Song

Overview

Installation

Simulate Data

Fit NNGP

Sampling and prediction

Fit NNGP joint model

Sampling and prediction

Fit MFA

Sampling and prediction

Fit MFA with linear response

Sampling and prediction

Compare with spNNGP

Compare the posterior mean with true value

Compare the predicted mean with true value

Compare the posterior variance with spNNGP

Compare the estimated mean for spatial parameters

Compare the estimated mean and variance for regression coefficients