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,SpatialPolygonsDataFrame,NULL,character
bym(formula, data, adjMat, region.id, ...)
# S4 method for formula,SpatialPolygonsDataFrame,missing,character
bym(formula, data, adjMat, region.id, ...)
# S4 method for formula,SpatialPolygonsDataFrame,nb,character
bym(formula,data,adjMat,region.id,...)
# S4 method for formula,data.frame,nb,character
bym(
formula,data,adjMat,region.id,
priorCI=list(sdSpatial=c(0.01,2),sdIndep=c(0.01,2)),
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.

priorCI

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 inla in the INLA package , such as control.inla

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.

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 priorCI argument can be a list containing elements named sdSpatial and sdIndep, each being a vector of length 2 with 2.5pct and 97.5pct quantiles for the prior distributions of the standard deviations $\sigma_1$ and $\sigma_2$ respectively. Gamma prior distributions for the precision parameters $1/\sigma_1^2$ and $1/\sigma_2^2$ yielding quantiles specified for the standard deviations are computed, and used with the model="bym" option to f .

The other possible format for priorCI is to have 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 priorCI$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 priorCI$propSpatial = c(0.5, 0.9) implies $pr(\phi < 0.5) = 0.9$.

Examples

Run this code

# NOT RUN {
data('kentucky')

# must have an internet connection to do the following
# }
# NOT RUN {
	larynxRates= cancerRates("USA", year=1998:2002,site="Larynx")
	dput(larynxRates)
# }
# NOT RUN {
larynxRates = structure(c(0, 0, 0, 0, 0, 0, 1e-06, 6e-06, 2.3e-05, 4.5e-05, 
9.9e-05, 0.000163, 0.000243, 0.000299, 0.000343, 0.000308, 0.000291, 
0.000217, 0, 0, 0, 0, 0, 1e-06, 1e-06, 3e-06, 8e-06, 1.3e-05, 
2.3e-05, 3.5e-05, 5.8e-05, 6.8e-05, 7.5e-05, 5.5e-05, 4.1e-05, 
3e-05), .Names = c("M_0", "M_5", "M_10", "M_15", "M_20", "M_25", 
"M_30", "M_35", "M_40", "M_45", "M_50", "M_55", "M_60", "M_65", 
"M_70", "M_75", "M_80", "M_85", "F_0", "F_5", "F_10", "F_15", 
"F_20", "F_25", "F_30", "F_35", "F_40", "F_45", "F_50", "F_55", 
"F_60", "F_65", "F_70", "F_75", "F_80", "F_85"), 
site = "Larynx", area = "USA, SEER", year = "1998-2002")

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

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

if( require("spdep", quietly=TRUE)) {

kBYM = bym(observed ~ offset(logExpected) + poverty, kentucky, 
	 priorCI = list(sdSpatial=c(0.1, 5), sdIndep=c(0.1, 5)),
	 control.mode=list(theta=c(3.52, 3.35),restart=TRUE))

	kBYM$par$summary
	
	if(requireNamespace('geostatsp', quietly=TRUE))
		kBYM$data$exc1 = geostatsp::excProb(
			kBYM$inla$marginals.fitted.bym, log(1.2)
			)
}  else {
	kBYM = list()
}




if(require('mapmisc', quietly=TRUE) & length(kBYM$data$fitted.exp)){

thecol = colourScale(kBYM$data$fitted.exp, 
	breaks=5, dec=1, opacity = 0.7)

map.new(kBYM$data)

# }
# NOT RUN {
kmap = openmap(kBYM$data)
plot(kmap,add=TRUE)
# }
# NOT RUN {
plot(kBYM$data, col=thecol$plot,add=TRUE)
legendBreaks("topleft", thecol)

}

# }

Run the code above in your browser using DataLab