bym-methods {diseasemapping}

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 inlain the INLApackage, such as control.inla

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