WrapKrigSp() R function from [CircSpaceTime]

Spatial interpolation using wrapped normal model.

WrapKrigSp function computes the spatial prediction for circular spatial data using samples from the posterior distribution of the spatial wrapped normal


WrapKrigSp(WrapSp_out, coords_obs, coords_nobs, x_obs)

Arguments

WrapSp_out: the functions takes the output of WrapSp function
coords_obs: coordinates of observed locations (in UTM)
coords_nobs: coordinates of unobserved locations (in UTM)
x_obs: observed values

Returns

a list of 3 elements

M_out: the mean of the associated linear process on the prediction locations coords_nobs (rows) over all the posterior samples (columns) returned by WrapSp
V_out: the variance of the associated linear process on the prediction locations coords_nobs (rows) over all the posterior samples (columns) returned by WrapSp
Prev_out: the posterior predicted values at the unobserved locations.

Implementation Tips

To facilitate the estimations, the observations x are centered around pi, and the posterior samples of x and posterior mean are changed back to the original scale

Examples


library(CircSpaceTime)
## auxiliary function
rmnorm<-function(n = 1, mean = rep(0, d), varcov){
  d <- if (is.matrix(varcov))
    ncol(varcov)
  else 1
  z <- matrix(rnorm(n * d), n, d) %*% chol(varcov)
  y <- t(mean + t(z))
  return(y)
}

####
# Simulation with exponential spatial covariance function
####
set.seed(1)
n <- 20
coords <- cbind(runif(n,0,100), runif(n,0,100))
Dist <- as.matrix(dist(coords))

rho     <- 0.05
sigma2  <- 0.3
alpha   <- c(0.5)
SIGMA   <- sigma2*exp(-rho*Dist)

Y <- rmnorm(1,rep(alpha,times=n), SIGMA)
theta <- c()
for(i in 1:n) {
  theta[i] <- Y[i]%%(2*pi)
}
rose_diag(theta)

#validation set
val <- sample(1:n,round(n*0.1))

set.seed(12345)
mod <- WrapSp(
  x       = theta[-val],
  coords    = coords[-val,],
  start   = list("alpha"      = c(.36,0.38),
                 "rho"     = c(0.041,0.052),
                 "sigma2"    = c(0.24,0.32),
                 "k"       = rep(0,(n - length(val)))),
  priors   = list("rho"      = c(0.04,0.08), #few observations require to be more informative
                  "sigma2"    = c(2,1),
                  "alpha" =  c(0,10)
  ),
  sd_prop   = list( "sigma2" = 0.1,  "rho" = 0.1),
  iter    = 1000,
  BurninThin    = c(burnin = 500, thin = 5),
  accept_ratio = 0.234,
  adapt_param = c(start = 40000, end = 45000, exp = 0.5),
  corr_fun = "exponential",
  kappa_matern = .5,
  parallel = FALSE,
  #With doParallel, bigger iter (normally around 1e6) and n_cores>=2 it is a lot faster
  n_chains = 2 ,
  n_cores = 1
)
check <- ConvCheck(mod)
check$Rhat ## close to 1 means convergence has been reached
## graphical check
par(mfrow = c(3,1))
coda::traceplot(check$mcmc)
par(mfrow = c(1,1))
##### We move to the spatial interpolation

Krig <- WrapKrigSp(
  WrapSp_out = mod,
  coords_obs =  coords[-val,],
  coords_nobs =  coords[val,],
  x_obs = theta[-val]
)

#### check the quality of the prediction using APE and CRPS
ApeCheck <- APEcirc(theta[val],Krig$Prev_out)
CrpsCheck <- CRPScirc(theta[val],Krig$Prev_out)

References

G. Jona-Lasinio, A .E. Gelfand, M. Jona-Lasinio, "Spatial analysis of wave direction data using wrapped Gaussian processes", The Annals of Applied Statistics, 6 (2012), 1478-1498

WrapKrigSp function