bym: Fit the BYM model

Description Usage Arguments Details Value Author(s) See Also Examples

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
## 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,SpatialPolygonsDataFrame,NULL,character'
bym(formula, data, adjMat, region.id, ...)
## S4 method for signature 'formula,SpatialPolygonsDataFrame,missing,character'
bym(formula, data, adjMat, region.id, ...)
## S4 method for signature 'formula,SpatialPolygonsDataFrame,nb,character'
bym(formula,data,adjMat,region.id,...)
## S4 method for signature '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 \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~Poisson(O_i λ_i)

log(μ_i) = X_i β + U_i

U_i ~ BYM(σ_1^2 , σ_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 σ_1 and σ_2 respectively. Gamma prior distributions for the precision parameters 1/σ_1^2 and 1/σ_2^2 yielding quantiles specified for the standard deviations are computed, and used with the model="bym" option to \Sexpr[results=rd]{c( '\\\code{f}', '\\\code{\\\\link[INLA]{f}}' )[1+requireNamespace('INLA', quietly=TRUE)] }.

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 σ_0 = sqrt(σ_1^2+σ_2^2). This prior is approximately exponential, and priorCI$sd = c(1, 0.01) specifies a prior probability pr(σ_0 > 1) = 0.01. The propSpatial element gives the prior for the ratio φ = σ_1/σ_0. Having priorCI$propSpatial = c(0.5, 0.9) implies pr(φ < 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

https://www.r-inla.org, \Sexpr[results=rd]{c( '\\\code{inla}', '\\\code{\\\\link[INLA]{inla}}' )[1+requireNamespace('INLA', quietly=TRUE)] } glgm, getSMR

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
data('kentucky')

# must have an internet connection to do the following
## Not run: 
	larynxRates= cancerRates("USA", year=1998:2002,site="Larynx")
	dput(larynxRates)

## End(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)

## End(Not run)

plot(kBYM$data, col=thecol$plot,add=TRUE)
legendBreaks("topleft", thecol)

}

diseasemapping documentation built on Sept. 10, 2021, 9:06 a.m.