glgm-methods: Generalized Linear Geostatistical Models

Description

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

Usage

## S3 method for class 'NULL,ANY,ANY,ANY':
glgm(formula=NULL, data,  grid, 
				covariates=NULL, 
				...)
## S3 method for class 'numeric,ANY,ANY,ANY':
glgm(formula, data,  grid, 
				covariates=NULL, 
				...)
## S3 method for class 'character,ANY,ANY,ANY':
glgm(formula, data,  grid, 
				covariates=NULL, 
				...)
## S3 method for class 'formula,ANY,numeric,ANY':
glgm(formula, data,  grid, 
				covariates=NULL, 
				...)
## S3 method for class 'formula,Raster,Raster,ANY':
glgm(formula, data,  grid, 
				covariates=NULL, buffer=0,
				...)
## S3 method for class 'formula,Spatial,Raster,ANY':
glgm(
	formula, data,  grid, 
				covariates=NULL, buffer=0,
				...)
## S3 method for class '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

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

priorCI

list with named elements of 0.025 and 0.975 quantiles of prior distributions, or a single vector giving the prior quantiles for the range parameter. List elements can be named: range for the range parameter (not the scale); sd

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

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, random.sd,random.X0.025quant, random.X0.5quant, random.X0.975quant, random.kld
predict.mean, predict.sd,predict.X0.025quant, predict.X0.5quant, predict.X0.975quant, predict.kld
predict.exp
predict.invlogitOnly supplied if a binomial response variable was used.
parameters contains a list with element
summary
and range and sd elements containing, for the range and standard deviation parameters respectively,
prior
posterior

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-shape)/gamma(shape) )*d^shape*besselK(d, shape)$$ $$d= |h|*sqrt(8*shape)/range$$ The range parameter produced by glgm multiplies the range parameter from INLA by the cell size.

Examples

Run this code

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

		

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)

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