diseasemapping (version 2.0.6)

bym-methods: Fit the BYM model

Description

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

Usage

# S4 method for formula,ANY,ANY,missing
bym(formula,data,adjMat,region.id,...)
# S4 method for formula,ANY,missing,missing
bym(formula,data,adjMat,region.id,...)
# S4 method for formula,SpatVector,NULL,character
bym(formula, data, adjMat, region.id, ...)
# S4 method for formula,SpatVector,missing,character
bym(formula, data, adjMat, region.id, ...)
# S4 method for formula,SpatVector,matrix,character
bym(formula,data,adjMat,region.id,...)
# S4 method for 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,...)

Value

A list containing

inla

results from the call to 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.

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 c( '\\code{inla} in the \\code{INLA} package', '\\code{\\link[INLA]{inla}}' )[1+requireNamespace('INLA', quietly=TRUE)] , such as c( '\\code{control.inla}', '\\code{\\link[INLA]{control.inla}}' )[1+requireNamespace('INLA', quietly=TRUE)]

Author

Patrick Brown

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

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