bym: Fit the BYM model

In diseasemapping: Modelling Spatial Variation in Disease Risk for Areal Data

Fit the BYM model

Description

Uses inla to fit a Besag, York and Mollie disease mapping model

Usage

## S4 method for signature 'formula,ANY,ANY,missing'
bym(formula,data,adjMat,region.id,...)
## S4 method for signature 'formula,ANY,missing,missing'
bym(formula,data,adjMat,region.id,...)
## S4 method for signature 'formula,SpatVector,NULL,character'
bym(formula, data, adjMat, region.id, ...)
## S4 method for signature 'formula,SpatVector,missing,character'
bym(formula, data, adjMat, region.id, ...)
## S4 method for signature 'formula,SpatVector,matrix,character'
bym(formula,data,adjMat,region.id,...)
## S4 method for signature 'formula,data.frame,matrix,character'
bym(
formula,data,adjMat,region.id,
prior=list(sd=c(0.1,0.5),propSpatial=c(0.5,0.5)),
family="poisson",formula.fitted=formula,...)

Arguments

formula	model formula, defaults to intercept-only model suitable for output from getSMR if data is a SpatialPolygonsDataFrame.
data	The observations and covariates for the model, can be output from getSMR.
adjMat	An object of class nb containing the adjacency matrix. If not supplied it will be computed from data, but is required if data is a SpatialPolygonDataFrame
region.id	Variable in data giving identifiers for the spatial regions. If not supplied, row numbers will be used.
prior	named list of vectors specifying priors, see Details
family	distribution of the observations, defaults to poisson
formula.fitted	formula to use to compute the fitted values, defaults to the model formula but may, for example, exclude individual-level covariates.
...	Additional arguments passed to \Sexpr[results=rd]{c( '\\\code{inla} in the \\\code{INLA} package', '\\\code{\\\\link[INLA]{inla}}' )[1+requireNamespace('INLA', quietly=TRUE)] }, such as \Sexpr[results=rd]{c( '\\\code{control.inla}', '\\\code{\\\\link[INLA]{control.inla}}' )[1+requireNamespace('INLA', quietly=TRUE)]}

Details

The Besag, York and Mollie model for Poisson distributed case counts is:

Y_i \sim Poisson(O_i \lambda_i)

\log(\mu_i) = X_i \beta + U_i

U_i \sim BYM(\sigma_1^2 , \sigma_2^2)

• Y_i is the response variable for region i, on the left side of the formula argument.

• O_i is the 'baseline' expected count, which is specified in formula on the log scale with \log(O_i) an offset variable.

• X_i are covariates, on the right side of formula

• U_i is a spatial random effect, with a spatially structured variance parameter \sigma_1^2 and a spatially independent variance \sigma_2^2.

The prior has elements named sd and propSpatial, which specifies model="bym2" should be used with penalized complexity priors. The sd element gives a prior for the marginal standard deviation \sigma_0 =\sqrt{\sigma_1^2+\sigma_2^2}. This prior is approximately exponential, and prior$sd = c(1, 0.01) specifies a prior probability pr(\sigma_0 > 1) = 0.01. The propSpatial element gives the prior for the ratio \phi = \sigma_1/\sigma_0. Having prior$propSpatial = c(0.5, 0.9) implies pr(\phi < 0.5) = 0.9.

Value

A list containing

inla	results from the call to \Sexpr[results=rd]{c( '\\\code{inla}', '\\\code{\\\\link[INLA]{inla}}' )[1+requireNamespace('INLA', quietly=TRUE)]}. Two additional elements are added: marginals.bym for the marginal distributions of the spatial random effects, and marginals.fitted.bym for the marginals of the fitted values.
data	A data.frame or SpatialPolygonsDataFrame containing posterior means and quantiles of the spatial random effect and fitted values.
parameters	Prior and posterior distributions of the two covariance parameters, and a table summary with posterior quantiles of all model parameters.

Author(s)

Patrick Brown

See Also

getSMR

Examples

data('kentucky')
kentucky = terra::unwrap(kentucky)

# get rid of under 10s
larynxRates = larynxRates[-grep("_(0|5)$",names(larynxRates))]

kentucky = getSMR(kentucky, larynxRates, larynx, regionCode="County")

if(requireNamespace('INLA')) {
  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)
}

kBYM <- try(
  bym(
    observed ~ offset(logExpected) + poverty, kentucky,
	  prior= list(sd=c(0.1, 0.5), propSpatial=c(0.5, 0.5))
    ), silent=TRUE)

if(length(grep("parameters", names(kBYM)))) {
  kBYM$parameters$summary
}



if( require("mapmisc", quietly=TRUE) && length(grep("parameters", names(kBYM))) ) {
  thecol = colourScale(kBYM$data$fitted.exp, breaks=5, dec=1, style="equal")
  terra::plot(kBYM$data, col=thecol$plot)
  legendBreaks("topleft", thecol)
}

diseasemapping documentation built on Sept. 22, 2023, 1:07 a.m.