spLMexact {spStack} | R Documentation |
Univariate Bayesian spatial linear model
Description
Fits a Bayesian spatial linear model with spatial process parameters and the noise-to-spatial variance ratio fixed to a value supplied by the user. The output contains posterior samples of the fixed effects, variance parameter, spatial random effects and, if required, leave-one-out predictive densities.
Usage
spLMexact(
formula,
data = parent.frame(),
coords,
cor.fn,
priors,
spParams,
noise_sp_ratio,
n.samples,
loopd = FALSE,
loopd.method = "exact",
verbose = TRUE,
...
)
Arguments
formula |
a symbolic description of the regression model to be fit. See example below. |
data |
an optional data frame containing the variables in the model.
If not found in |
coords |
an |
cor.fn |
a quoted keyword that specifies the correlation function used
to model the spatial dependence structure among the observations. Supported
covariance model key words are: |
priors |
a list with each tag corresponding to a parameter name and containing prior details. |
spParams |
fixed value of spatial process parameters. |
noise_sp_ratio |
noise-to-spatial variance ratio. |
n.samples |
number of posterior samples to be generated. |
loopd |
logical. If |
loopd.method |
character. Ignored if |
verbose |
logical. If |
... |
currently no additional argument. |
Details
Suppose \chi = (s_1, \ldots, s_n)
denotes the n
spatial locations the response y
is observed. With this function, we
fit a conjugate Bayesian hierarchical spatial model
\begin{aligned}
y \mid z, \beta, \sigma^2 &\sim N(X\beta + z, \delta^2 \sigma^2 I_n), \quad
z \mid \sigma^2 \sim N(0, \sigma^2 R(\chi; \phi, \nu)), \\
\beta \mid \sigma^2 &\sim N(\mu_\beta, \sigma^2 V_\beta), \quad
\sigma^2 \sim \mathrm{IG}(a_\sigma, b_\sigma)
\end{aligned}
where we fix the spatial process parameters \phi
and \nu
, the
noise-to-spatial variance ratio \delta^2
and the hyperparameters
\mu_\beta
, V_\beta
, a_\sigma
and b_\sigma
. We utilize
a composition sampling strategy to sample the model parameters from their
joint posterior distribution which can be written as
p(\sigma^2, \beta, z \mid y) = p(\sigma^2 \mid y) \times
p(\beta \mid \sigma^2, y) \times p(z \mid \beta, \sigma^2, y).
We proceed by first sampling \sigma^2
from its marginal posterior,
then given the samples of \sigma^2
, we sample \beta
and
subsequently, we sample z
conditioned on the posterior samples of
\beta
and \sigma^2
(Banerjee 2020).
Value
An object of class spLMexact
, which is a list with the
following tags -
- samples
a list of length 3, containing posterior samples of fixed effects (
beta
), variance parameter (sigmaSq
), spatial effects (z
).- loopd
If
loopd=TRUE
, contains leave-one-out predictive densities.- model.params
Values of the fixed parameters that includes
phi
(spatial decay),nu
(spatial smoothness) andnoise_sp_ratio
(noise-to-spatial variance ratio).
The return object might include additional data used for subsequent prediction and/or model fit evaluation.
Author(s)
Soumyakanti Pan span18@ucla.edu,
Sudipto Banerjee sudipto@ucla.edu
References
Banerjee S (2020). "Modeling massive spatial datasets using a conjugate Bayesian linear modeling framework." Spatial Statistics, 37, 100417. ISSN 2211-6753. doi:10.1016/j.spasta.2020.100417.
Vehtari A, Simpson D, Gelman A, Yao Y, Gabry J (2024). "Pareto Smoothed Importance Sampling." Journal of Machine Learning Research, 25(72), 1-58. URL https://jmlr.org/papers/v25/19-556.html.
See Also
Examples
# load data
data(simGaussian)
dat <- simGaussian[1:100, ]
# setup prior list
muBeta <- c(0, 0)
VBeta <- cbind(c(1.0, 0.0), c(0.0, 1.0))
sigmaSqIGa <- 2
sigmaSqIGb <- 0.1
prior_list <- list(beta.norm = list(muBeta, VBeta),
sigma.sq.ig = c(sigmaSqIGa, sigmaSqIGb))
# supply fixed values of model parameters
phi0 <- 3
nu0 <- 0.75
noise.sp.ratio <- 0.8
mod1 <- spLMexact(y ~ x1, data = dat,
coords = as.matrix(dat[, c("s1", "s2")]),
cor.fn = "matern",
priors = prior_list,
spParams = list(phi = phi0, nu = nu0),
noise_sp_ratio = noise.sp.ratio,
n.samples = 100,
loopd = TRUE, loopd.method = "exact")
beta.post <- mod1$samples$beta
z.post.median <- apply(mod1$samples$z, 1, median)
dat$z.post.median <- z.post.median
plot1 <- surfaceplot(dat, coords_name = c("s1", "s2"),
var_name = "z_true")
plot2 <- surfaceplot(dat, coords_name = c("s1", "s2"),
var_name = "z.post.median")
plot1
plot2