excProb: Exceedance probabilities

Description

Calculate exceedance probabilities pr(X > threshold) from a fitted geostatistical model.

Usage

excProb(x, threshold=0, random=FALSE, template=NULL, templateIdCol=NULL,
nuggetInPrediction=TRUE)

Arguments

Output from either the lgm or glgm functions, or a list of two-column matrices with columns named x and y containing the posterior distributions of random effects, as produced by inla.

threshold

the value which the exceedance probability is calculated with respect to.

random

Calculate exceedances for the random effects, rather than the predicted observations (including fixed effects).

template

A Raster or SpatialPolygonsDataFrame or SpatialPointsDataFrame object which the results will be contained in.

templateIdCol

The data column in template corresponding to names of marginals

nuggetInPrediction

If TRUE, calculate exceedance probabilities of new observations by adding the nugget effect. Otherwise calculate probabilities for the latent process. Ignored if x is output from glgm.

Value

Either a vector of exceedance probabilities or an object of the same class as template.

Details

When x is the output from lgm, pr(Y>threshold) is calculated using the Gaussian distribution using the Kriging mean and conditional variance. When x is from the glgm function, the marginal posteriors are numerically integrated to obtain pr(X > threshold).

Examples

Run this code

# NOT RUN {
	data('swissRain')
	swissFit =  lgm("rain", swissRain, grid=30, 
	boxcox=0.5,fixBoxcox=TRUE,	covariates=swissAltitude)
	swissExc = excProb(swissFit, 20)
	mycol = c("green","yellow","orange","red")
	mybreaks = c(0, 0.2, 0.8, 0.9, 1)
	plot(swissBorder)
	plot(swissExc, breaks=mybreaks, col=mycol,add=TRUE,legend=FALSE)
	plot(swissBorder, add=TRUE)
	legend("topleft",legend=mybreaks, col=c(NA,mycol))
# }
# NOT RUN {
	swissRain$sqrtrain = sqrt(swissRain$rain)
	swissFit2 =  glgm(formula="sqrtrain",data=swissRain, grid=40, 
	covariates=swissAltitude,family="gaussian")
	swissExc = excProb(swissFit2, threshold=sqrt(30))
	swissExc = excProb(swissFit2$inla$marginals.random$space, 0,
		template=swissFit2$raster)
# }
# NOT RUN {
# }