glgm-methods {geostatsp} | R Documentation |
Generalized Linear Geostatistical Models
Description
Fits a generalized linear geostatistical model or a log-Gaussian Cox process
using inla
Usage
## S4 method for signature 'ANY,ANY,ANY,ANY'
glgm(formula, data, grid, covariates, buffer=0, shape=1, prior, ...)
## S4 method for signature 'formula,SpatRaster,ANY,ANY'
glgm(formula, data, grid, covariates, buffer=0, shape=1, prior, ...)
## S4 method for signature 'formula,SpatVector,ANY,ANY'
glgm(formula, data, grid, covariates, buffer=0, shape=1, prior, ...)
## S4 method for signature 'formula,data.frame,SpatRaster,data.frame'
glgm(formula, data, grid, covariates, buffer=0, shape=1, prior, ...)
lgcp(formula=NULL, data, grid, covariates=NULL, border, ...)
Arguments
data |
An object of class |
grid |
Either an integer giving the number of cells in the x direction, or a raster object which will be used for the spatial random effect. If the cells in the raster are not square, the resolution in the y direction will be adjusted to make it so. |
covariates |
Either a single raster, a list of rasters or a raster stack containing covariate values used when
making spatial predictions. Names of the raster layers or list elements correspond to names in the formula. If
a covariate is missing from the data object it will be extracted from the rasters. Defaults to |
formula |
Model formula, defaults to a linear combination of each of the layers in the |
prior |
list with elements named |
shape |
Shape parameter for the Matern correlation function, must be 1 or 2. |
buffer |
Extra space padded around the data bounding box to reduce edge effects. |
border |
boundary of the region on which an LGCP is defined, passed to |
... |
Additional options passed to
|
Details
This function performs Bayesian inference for generalized linear geostatistical models with INLA. The Markov random field approximation on a regular lattice is used for the spatial random effect. The range parameter is the distance at which the correlation is 0.13, or
where is the shape parameter. The range parameter produced by
glgm
multiplies the range parameter from INLA
by the cell size.
Elements of prior
can be named range
, sd
, or sdObs
. Elements can consist of:
a single value giving the prior median for penalized complexity priors (exponential on the sd or 1/range).
a vector
c(u=a, alpha=b)
giving an quantile and probability for pc priors. For standard deviations alpha is an upper quantile, for the range parameter b = pr(1/range > 1/a).a vector
c(lower=a, upper=b)
giving a 0.025 and 0.975 quantiles for the sd or range.a list of the form
list(prior='loggamma', param=c(1,2))
passed directly to inla.a two-column matrix of prior densities for the sd or range.
Value
A list with two components named inla
, raster
, and parameters
. inla
contains the results of the call to the
inla
function. raster
is a raster stack with the following layers:
random. |
mean, sd, X0.0??quant: Posterior mean, standard deviation, and quantiles of the random effect |
predict. |
mean, sd, X0.0??quant: same for linear predictors, on the link scale |
predict.exp |
posterior mean of the exponential of the linear predictor |
predict.invlogit |
Only supplied if a binomial response variable was used. |
parameters
contains a list with elements:
summary |
a table with parameter estimates and posterior quantiles |
range , sd |
prior and posterior distributions of range and standard deviations |
See Also
inla
in the INLA
package,
https://www.r-inla.org
Examples
# geostatistical model for the swiss rainfall data
if(requireNamespace("INLA", quietly=TRUE) ) {
INLA::inla.setOption(num.threads=2)
# not all versions of INLA support blas.num.threads
try(INLA::inla.setOption(blas.num.threads=2), silent=TRUE)
}
require("geostatsp")
data("swissRain")
swissRain = unwrap(swissRain)
swissAltitude = unwrap(swissAltitude)
swissBorder = unwrap(swissBorder)
swissRain$lograin = log(swissRain$rain)
swissFit = glgm(formula="lograin", data=swissRain,
grid=30,
covariates=swissAltitude, family="gaussian",
buffer=2000,
prior = list(sd=1, range=100*1000, sdObs = 2),
control.inla = list(strategy='gaussian')
)
if(!is.null(swissFit$parameters) ) {
swissExc = excProb(swissFit, threshold=log(25))
swissExcRE = excProb(swissFit$inla$marginals.random$space,
log(1.5),template=swissFit$raster)
swissFit$parameters$summary
matplot(
swissFit$parameters$range$postK[,'x'],
swissFit$parameters$range$postK[,c('y','prior')],
type="l", lty=1, xlim = c(0, 1000),
xlab = 'km', ylab='dens')
legend('topright', lty=1, col=1:2, legend=c('post','prior'))
plot(swissFit$raster[["predict.exp"]])
mycol = c("green","yellow","orange","red")
mybreaks = c(0, 0.2, 0.8, 0.95, 1)
plot(swissBorder)
plot(swissExc, breaks=mybreaks, col=mycol,add=TRUE,legend=FALSE)
plot(swissBorder, add=TRUE)
legend("topleft",legend=mybreaks, fill=c(NA,mycol))
plot(swissBorder)
plot(swissExcRE, breaks=mybreaks, col=mycol,add=TRUE,legend=FALSE)
plot(swissBorder, add=TRUE)
legend("topleft",legend=mybreaks, fill=c(NA,mycol))
}
# a log-Gaussian Cox process example
myPoints = vect(cbind(rbeta(100,2,2), rbeta(100,3,4)))
mycov = rast(matrix(rbinom(100, 1, 0.5), 10, 10), extent=ext(0, 1, 0, 1))
names(mycov)="x1"
if(requireNamespace("INLA", quietly=TRUE) ) {
INLA::inla.setOption(num.threads=2)
# not all versions of INLA support blas.num.threads
try(INLA::inla.setOption(blas.num.threads=2), silent=TRUE)
}
res = lgcp(
formula=~factor(x1),
data=myPoints,
grid=squareRaster(ext(0,1,0,1), 20), covariates=mycov,
prior=list(sd=c(0.9, 1.1), range=c(0.4, 0.41),
control.inla = list(strategy='gaussian'), verbose=TRUE)
)
if(length(res$parameters)) {
plot(res$raster[["predict.exp"]])
plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5)
}