# glgm: Generalized Linear Geostatistical Models In geostatsp: Geostatistical Modelling with Likelihood and Bayes

## Description

Fits a generalized linear geostatistical model or a log-Gaussian Cox process using `inla`

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21``` ```## S4 method for signature 'NULL,ANY,ANY,ANY' glgm(formula, data, grid, covariates=NULL, ...) ## S4 method for signature 'numeric,ANY,ANY,ANY' glgm(formula, data, grid, covariates=NULL, ...) ## S4 method for signature 'character,ANY,ANY,ANY' glgm(formula, data, grid, covariates=NULL, ...) ## S4 method for signature 'formula,ANY,numeric,ANY' glgm(formula, data, grid, covariates=NULL, ...) ## S4 method for signature 'formula,Raster,Raster,ANY' glgm(formula, data, grid, covariates=NULL, buffer=0, ...) ## S4 method for signature 'formula,Spatial,Raster,ANY' glgm(formula, data, grid, covariates=NULL, buffer=0, ...) ## S4 method for signature 'formula,data.frame,Raster,data.frame' glgm(formula, data, grid, covariates=NULL, shape=1, priorCI=NULL, mesh=FALSE,...) lgcp(formula=NULL, data, grid, covariates=NULL, border, ...) ```

## Arguments

 `data` An object of class ` SpatialPointsDataFrame` containing the data. `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 `NULL` for an intercept-only model. `formula` Model formula, defaults to a linear combination of each of the layers in the `covariates` object. The spatial random effect should not be supplied but the default can be overridden with a ` f(space,..)` term. For `glgm` the response variable defaults to the first variable in the `data` object, and `formula` can be an integer or character string specifying the response variable. For `lgcp`, the formula should be one-sided. `priorCI` list with elements named `range`, `sd`, `sdNugget`, `intercept`, or `betas`. See Details. `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. `mesh` Currently unimplemented options for using a mesh rather than a grid for the Markov random field approximation. `border` boundary of the region on which an LGCP is defined, passed to `mask` `...` Additional options passed to `inla`

## 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

cov[U(s+h), U(s)] = (2^{1-ν}/Gamma(ν)) d^ν besselK(d, ν)

d= |h| √{8 ν}/range

where ν is the shape parameter. The range parameter produced by `glgm` multiplies the range parameter from `INLA` by the cell size.

Elements of `priorCI` can be named:

• `range` for the range parameter (not the scale)

• `sd` for the standard deviation (not the variance or precision)

• `sdNugget` for the standard deviation of the observation error for Gaussian data

• `gammaShape` for the shape parameter of a Gamma distributed observations

• `intercept` for the intercept

• `betas` for the regression coefficients

Each element is either a vector of length 2 or a two-column matrix.

• `c(lower=a, upper=b)` specifies 0.025 and 0.975 quantiles of prior distributions

• `c(u=a, alpha=b)` gives a penalized complexity prior for standard deviation parameters. pr(σ > a) = b

• A matrix has a first column of parameter values and a second column of densities

## 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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102``` ```# use the 'safe' version of INLA on linux systems if(Sys.info()['sysname'] =='Linux' & requireNamespace("INLA")) { INLA::inla.setOption(inla.call = system.file(paste("bin/linux/", ifelse(.Machine\$sizeof.pointer == 4, "32", "64"), 'bit/inla.static', sep=''), package="INLA")) } # geostatistical model for the swiss rainfall data require("geostatsp") data("swissRain") swissRain\$lograin = log(swissRain\$rain) swissFit = glgm(formula="lograin", data=swissRain, grid=30, covariates=swissAltitude, family="gaussian", buffer=20000, priorCI=list(sd=c(0.01, 5), range=c(50000,500000), sdNugget = c(0.01, 5)), control.mode=list(theta=c(1.6,-0.25,2.9),restart=TRUE) ) 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 plot(swissFit\$parameters\$range\$prior,type="l", ylim=c(0,max(swissFit\$parameters\$range\$posterior[,"y"])), xlim=c(0, 500000)) abline(v=swissFit\$parameters\$range\$userPriorCI,col="yellow") abline(v=swissFit\$parameters\$range\$priorCI,col="orange") lines(swissFit\$parameters\$range\$posterior, col='blue') } if(interactive() | Sys.info()['user'] =='patrick') { 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)) } ## Not run: load(url("http://www.filefactory.com/file/frd1mhownd9/n/CHE_adm0_RData")) thenames = GNcities(bbox(gadm),max=12) swissTiles = openmap(bbox(gadm),zoom=8,type="nps") par(mar=c(0,0,0,0)) plot(gadm) plot(swissTiles, add=TRUE) library('RColorBrewer') mycol=rev(brewer.pal(4,"RdYlGn")) plot(mask( projectRaster(swissExc, crs=proj4string(gadm)), gadm), breaks = c(0, 0.2, 0.8, 0.95, 1.00001), col=mycol, alpha=0.5,add=TRUE) plot(gadm, add=TRUE, lwd=2, border='blue') points(thenames,cex=0.5) text(thenames, labels=thenames\$name,pos=3, vfont=c("gothic german","plain"),cex=1.5) ## End(Not run) # a log-Gaussian Cox process example if(interactive() | Sys.info()['user'] =='patrick') { myPoints = SpatialPoints(cbind(rbeta(100,2,2), rbeta(100,3,4))) myPoints@bbox = cbind(c(0,0), c(1,1)) mycov = raster(matrix(rbinom(100, 1, 0.5), 10, 10), 0, 1, 0, 1) names(mycov)="x1" res = lgcp(data=myPoints, grid=20, covariates=mycov, formula=~factor(x1), priorCI=list(sd=c(0.9, 1.1), range=c(0.4, 0.41)) ) plot(res\$raster[["predict.exp"]]) plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5) } ```