bym-methods {diseasemapping} | R Documentation |
Uses inla to fit a Besag, York and Mollie disease mapping model
## 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,...)
formula |
model formula, defaults to intercept-only model suitable for
output from |
data |
The observations and covariates for the model, can be output from
|
adjMat |
An object of class |
region.id |
Variable in |
priorCI |
named list of vectors specifying priors, see Details |
family |
distribution of the observations, defaults to |
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
|
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)
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 σ_1^2 and a spatially independent variance σ_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
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
σ_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.
A list containing
inla |
results from the call to
|
data |
A |
parameters |
Prior and posterior distributions of the two covariance parameters, and a table summary with posterior quantiles of all model parameters. |
Patrick Brown
http://r-inla.org,
inla
glgm
, getSMR
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) }