diseasemapping-package: Disease Mapping

Description Author(s) Examples

Description

Functions for calculating observed and expected counts by region, and manipulating posterior samples from Bayesian models produced by glmmBUGS.

Author(s)

Patrick Brown

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
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
# creating SMR's
data('kentucky')

if(FALSE) {
	# must have an internet connection to do the following
	larynxRates= cancerRates("USA", year=1998:2002,site="Larynx")
	# get rid of under 10's
	larynxRates = larynxRates[-grep("_(0|5)$",names(larynxRates))]
	dput(larynxRates)
} else {
larynxRates = structure(c(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, 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_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_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"))

}

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

if(require('mapmisc', quietly=TRUE)) {
	mycol = colourScale(kentucky2$SMR, breaks=9, 
	dec=-log10(0.5), style='equal', transform='sqrt')
	plot(kentucky2, col=mycol$plot)
	legendBreaks('topleft', mycol)
}

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


kBYM = bym(observed ~ offset(logExpected) + poverty, kentucky2, 
	 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

}

# an example of a spatial model with glmmBUGS

## Not run: 
# run the model
library('spdep')
popDataAdjMat = poly2nb(ontario,row.names=as.character(ontario[["CSDUID"]]))

library('glmmBUGS')
forBugs = glmmBUGS(formula=observed + logExpected ~ 1,
  effects="CSDUID", family="poisson", spatial=popDataAdjMat,
  data=ontario@data)
startingValues = forBugs$startingValues
source("getInits.R")
library('R2WinBUGS')
ontarioResult = bugs(forBugs$ragged, getInits, parameters.to.save = names(getInits()),
    model.file="model.bug", n.chain=3, n.iter=100, n.burnin=10, n.thin=2,
      program="winbugs", debug=TRUE)

ontarioParams = restoreParams(ontarioResult, forBugs$ragged)
ontarioSummary = summaryChain(ontarioParams)

# merge results back in to popdata
ontario = mergeBugsData(ontario, ontarioSummary)

## End(Not run)


# running the same thing with INLA
## Not run: 
library('INLA')
# get rid of regions with no neighbours
ontario2 = ontario[! as.character(ontario$CSDUID) 
c("3510005" ,"3501007", "3537001", "3551031", "3560065", "3560062"),]
popDataAdjMat2 = poly2nb(ontario2,
	row.names=as.character(ontario2[["CSDUID"]]))
nb2INLA("nb.graph",popDataAdjMat2)

ontario2$CSDUID = as.character(ontario2$CSDUID)
 
prior.iid=prior.besag=c(.42,.00015)
formula.bym = observed ~  f(CSDUID, 
    model = "bym", graph = "nb.graph", values=CSDUID ,
    param = c(prior.iid, prior.besag))
 
result1.bym = inla(formula.bym,family="poisson",data=ontario2@data,
	offset=ontario2@data$logExpected,
    verbose=FALSE, keep = TRUE, 
    control.predictor=list(compute=TRUE))
    

tomerge = result1.bym$summary.random$CSDUID
rownames(tomerge) = tomerge$ID
 
ontario2@data =   cbind(ontario2@data, 
	tomerge[as.character(ontario2$CSDUID),] )   
	
require('mapmisc') 
mycol = colourScale(ontario2$mean, breaks=9, 
	dec=1, style='equal', transform='sqrt')
plot(ontario2, col=mycol$plot)
legendBreaks('topleft', mycol)
	

## End(Not run)

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