excProb: Exceedance probabilities

In geostatsp: Geostatistical Modelling with Likelihood and Bayes

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

x	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 SpatRaster or SpatVector 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.

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).

Value

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

Examples

	data('swissRain')
	swissRain = unwrap(swissRain)
	swissAltitude = unwrap(swissAltitude)
	swissBorder = unwrap(swissBorder)
	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))



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)

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

geostatsp documentation built on Dec. 24, 2024, 3 a.m.