Analysing Spatial Data in R: Worked examples: …

1 Analysing Spatial Data in R: Worked examples:( bayesian ) disease mapping IIVirgilio G omez-RubioDepartment of Epidemiology and Public HeathImperial College LondonLondon, UK31 August 2007 bayesian disease mappingIBayesian Estimation in disease mapping has been one of theleading topics in Spatial statistics in the last 20 yearsIBayesian Hierarchical Models can be used to model complexdata structuresIThe bayesian approach offers an easy approach to theestimation of complex models via Markov Chain Monte CarloISpatial analysis of routinoulsy collected health data isstandard practise nowadaysISpatio-temporal models can be usedIWaller & Gotway (2004) and Banerjee et al. (2003) accountfor a comprehensive summary on Spatial modelsBayesian InferenceIBayesian Inference is based on estimating the probabilitydensity of the parameters in the modelafterobserving thedata, , their posterior distributions:p( |y)Ip( |y) is usually difficult to derive.

2 P( |y) =p(y| )p( ) p(y| )p( ) p(y| )p( )Ip(y| ) is the likelihood of the model, which reflects therelationship between the data and the parametersIp( ) is thepriordistribution of the parameters, which reflectsthe initial information on the parametersIUsually,p( |y) is computed by simulation usingMarkov ChainMonte CarlotechniquesIWinBUGS in a generic software to fit a wide range of uses the Gibbs sampler for of bayesian InferenceISuitable framework to deal with a large number of problemsIPriors can be used to account for initial information (forexample, Spatial dependence)IIf no prior information is available, vague (or non-informative)priors can be used so that the posterior distributionwill onlydependon the data and the models can be used: bayesian Hierarchical ModelsIComplex effects, such as Spatial and/or temporal dependence,can be modeled easilyIWhen the posterior distribution is not in a closed form,different simulation techniques can be used to values are treated similarly as the parameters in themodelMarkov Chain Monte Carlo/Gibbs samplerMCMC aims at simulating a series of values for the parameters inthe model, so that, in the end, these values will be draws from theposterior initial values to every parameter in the model (andmissing values)IAt every step, Gibbs sampler simulates from the fullconditional distribution.

3 P( i| i,y)IAfter aburn inperiod, the simulated values are draws fromthe posteriorp( |y)IConvergence of the simulated values should be assessedWinBUGSIBUGS stands forBayesian inference Using Gibbs SamplerIDeveloped at the MRC and Imperial College LondonIProvides a generic language to bayesian Hierarchical modelsIModels can be specified graphically as wellISeveral utilities to assess the convergence of the chain anddisplay resultsIGeoBUGS is an extension to deal with Spatial models andmapsIPkBUGS is another extension to deal with PharmacokineticsmodelsIA developer interface has been included so that the user canextend the range of functions availableIOpenBUGS is theopen sourceversion of WinBUGSC alling WinBUGS fromRIPackagesR2 WinBUGSandBRUGScan call WinBUGS andOpenBUGS fromRIR2 WinBUGS calls WinBUGS using the scripting language andthen reads the output log fileIBRUGSis an interface to the

4 Actual OpenBUGS (NOTWinBUGS) routinesIR2 WinBUGScan run on several platforms (Windows,Linux/Unix, Mac)IOther alternatives to call WinBUGS externally in differentways are available Cancer Data revisitedWe specification (using the BUGS language)IMortality Data (in alist)ISpatial data describing the neighbourhood structure, in aspecific formatIInitial values of the parametersIOptionally, we may want to export the map information to beused within WinBUGSB ayesian Spatial ModellingOi Poisson( i) i= iEilog( i) = + Xi+ui+viui Normal(0, 2u)vi|v i Normal( j ivj/ni, 2v/ni)f( ) 1f( ) 1 2u Gamma 1(.001,.001) 2v Gamma 1(.001,.001)v1v2v3v4v5viui 2u 2v 2 Xi iEi i= iEiOiJ. Besag, J. York, A. Mollie (1991). bayesian image restoration, with twoapplications in Spatial statistics (with discussion).Annals of the Institute ofStatistical Mathematics43(1), 1-59 Model specification using the BUGS languagemodel{for(i in 1:N){O[i] ~ dpois(mu[i])mu[i]<-theta[i] * E[i]log(theta[i]) <- alpha + beta[1]*PCTAGE65P[i]+beta[2]*PCTOWNHOME[ i]+beta[3]*AVGIDIST[i]+u[i] + v[i]u[i] ~ dnorm(0, precu)SMR[i]<- O[i] / E[i]prob[i]<-step(theta[i]-1)}v[1:N]~ (adj[], weights[], num[], precv)alpha~dflat()for(i in 1.)}

5 3) {beta[i] ~dflat()}precu~dgamma( , )precv~dgamma( , )sigmau<-1/precusigmav<-1/precv}Preparing Read maps> library(maptools)> nymap <- readShapePoly("NY8_utm18") Create list of observed, expected> nymap$EXP <- nymap$POP8 * sum(nymap$Cases)/sum(nymap$POP8) Create adjacency matrix> library(spdep)> nynb <- poly2nb(nymap) Create weights> nyWBweights <- nb2WB(nynb)> d <- c(list(O = nymap$Cases, E = nymap$EXP), N = 281,+ list(PCTAGE65P = nymap$PCTAGE65P, PCTOWNHOME = nymap$PCTOWNHOME,+ AVGIDIST = nymap$AVGIDIST))> inits1 <- list(alpha = 1, beta = c(0, 0, 0), u = rep(0,+ 281), v = rep(0, 281), precu = 1, precv = 1)> inits2 <- list(alpha = 10, beta = c(1, 1, 1), u = rep(1,+ 281), v = rep(1, 281), precu = , precv = )Calling WinBUGS Call WinBUGS> library(R2 WinBUGS)> mfile <- paste(getwd(), " ", sep = "", collapse = "")> tdir <- paste(getwd(), "/NYoutput", sep = "", collapse = "")> (tdir)> res <- bugs(data = c(d, nyWBweights), inits = list(inits1,+ inits2), = c("u", "v", "theta",+ "prob", "sigmau", "sigmav"), = mfile,+ = tdir, = 3, = 2,+ = 6000, = 3000) Add results to map object> nymap$prob <- res$mean$prob> nymap$theta <- res$mean$theta> nymap$u <- res$mean$u> nymap$v <- res$mean$v> logfile <- paste(getwd(), "/ ", sep = "",+ collapse = "")> reslog <- (file = logfile) mapping the resultsSmoothed Relative map 1 Exporting the data to work directly with the maps with spdep> sp2WB(map = nymap, file = " ") map with WB first, and reboot (from R2 WinBUGS) to create the files withdata and initial values> (d)> (" ", " ")

6 > (nyWBweights)> (" ", " ")> (inits1)> (" ", " ")> (inits2)> (" ", " ")Running WinBUGS all needed files in model (health and Spatial ) model (burn in period) of interest and of the statistics of the parameters of the on a mapFurther referencesIS. Banerjee, Carlin and Gelfand (2003).HierarchicalModeling and Analysis for Spatial Data. Chapman & Lawson, Browne and Vidal Rodeiro (2003). disease mapping with WinBUGS and MLwiN. Wiley & : language: Spiegelhalter, Best, Carlin and A. Van derLinde (2002). bayesian Measures of Model Complexity and Fit(with Discussion),Journal of the Royal Statistical Society,Series B64(4), Waller and Gotway (2004).Applied SpatialStatistics for Public Health Data. Wiley & Sons.