Analysing Spatial Data in R: Worked examples: disease ...

1 Analysing Spatial data in R: Worked examples: disease mapping IRoger BivandDepartment of EconomicsNorwegian School of Economics and Business AdministrationBergen, Norway31 August 2007 disease mappingIDisease mapping is one of the areas of applied statistics thatis developing most rapidly and is in most demandIIt involved both Spatial data and methods ranging fromvisualisation to advanced statisticsIWaller and Gotway discuss many approaches, using, amongother data sets, the classic Scottish lip cancer set by countyIBanerjee, Carlin and Gelfand also discuss the same data set,providing code for a WinBUGS exampleScottish lip cancer data setReading the data as downloaded from Waller and Gotway swebsite, we need first to import the county boundaries and assignthe correct CRS.

2 Then transform to the British National Grid:> library(rgdal)> scot_LL <- readOGR(".", "scot")OGR data source with driver: ESRI ShapefileSource: ".", layer: "scot"with 56 rows and 2 columns> proj4string(scot_LL) <- CRS("+proj=longlat ellps=WGS84")> EPSG <- make_EPSG()> EPSG[grep("British National Grid", EPSG$note), 1:2]code note2418 27700 # OSGB 1936 / British National Grid> scot_BNG0 <- spTransform(scot_LL, CRS("+init=epsg:27700"))Scottish lip cancer data setThe data as given by Clayton and Kaldor (1987) use the samecounty ID values, but in a different order, so we need to merge thedata frame so as to output the same order as the county boundarypolygons:> library(maptools)> scot_dat <- (" ", skip = 1)> names(scot_dat) <- c("District", "Observed", "Expected", "PcAFF",+ "Latitude", "Longitude")> (scot_dat) <- formatC(scot_dat$District, width = 2,+ flag = "0")> ID <- formatC(scot_BNG0$ID, width = 2, flag = "0")> scot_BNG1 <- spChFIDs(scot_BNG0, ID)> scot_BNG <- spCbind(scot_BNG1, scot_dat[match(ID, (scot_dat)),+ ])Scottish lip cancer data setThe original data table includes expected numbers of cases, whichare based on age effects using the method of Mantel and Stark(1968), and are proportional to a population at risk after sucheffects have been taken into account.

3 Setting up some palettes, wecan look at maps of observed and expected counts:> bluepal <- colorRampPalette(c("azure1", "steelblue4"))> brks <- c(0, 5, 10, 25, 50, 100)> library(classInt)> O_CI <- classIntervals(scot_BNG$Observed, style = "fixed", fixedBreaks = brks)> E_CI <- classIntervals(scot_BNG$Expected, style = "fixed", fixedBreaks = brks)Lip cancer 1975 1980: observed and expectedObservedunder 55 1010 2525 50over 50 Expectedunder 55 1010 2525 50over 50 Probability maps and smoothingIChoynowski proposed a method for making probability mapsin 1959, still described by contemporary authors it splitsbetween upper and lower tailsIIt is also possible to say justppois(Observed, Expected)toget a similar result. If the observations are distributed asassumed, this will indicate which regions appear unusualIProbability maps typically mirror the Spatial pattern of relativerisks.

4 Using Empirical Bayes smoothing will use the populationat risk to adjust relative risksIEmpirical Bayes smoothing is implemented for themethod-of-moments global and local cases inspdepand forthe ML global case inDClusterChoynowski probability mapThe legacy Choynowski approach isimplemented in spdep> library(spdep)> ch <- choynowski(scot_BNG$Observed,+ scot_BNG$Expected)Choynowski probability maplowN/ShighPoisson probability mapThe Poisson probability map in onetail is identical with the Choynowskiapproach, but shows both tails ratherthan folding them together> pm <- probmap(scot_BNG$Observed,+ scot_BNG$Expected)> names(pm)[1] "raw" "expCount" "relRisk"[4] "pmap"> scot_BNG$SMR <- pm$relRiskPoisson probability mapunder (11) (1) (4) (15) (3) (4)over (18)Weaknesses of probability mapsIf the underlying distribution of thedata does not agree with ourassumption, we may get severalpossible processes mixed up,overdispersion with spatialdependence:> table(findInterval(pm$pmap,+ seq(0, 1, 1/10)))1 2 4 5 6 7 8 9 1012 4 2 4 2 3 4 3 22 Poisson probability mappm$ Bayes smoothingThe method of moments approach is implemented inspdep, whilethe maximum likelihood approach is implemented inDCluster:> eb1 <- EBest(scot_BNG$Observed, scot_BNG$Expected)> unlist(attr(eb1, "parameters"))a > scot_BNG$EB_mm <- eb1$estmm * 100> library(DCluster)> res <- empbaysmooth(scot_BNG$Observed, scot_BNG$Expected)> unlist(res[2.)]

5 3])nu > scot_BNG$EB_ml <- res$smthrr * 100 NeighboursIn order to investigate spatialdependence, we need a list ofneighbours. We take the list given byClayton and Kaldor, re-ordering tosuit out data order:> CK_nb <- (" ",+ = scot_BNG$District)> CK_nbNeighbour list object:Number of regions: 56 Number of nonzero links: 264 Percentage nonzero weights: number of links: Empirical Bayes smoothingIf instead of shrinking to a globalrate, we shrink to a local rate, wemay be able to take unobservedheterogeneity into account; here weuse the list of neighbours:> eb2 <- EBlocal(scot_BNG$Observed,+ scot_BNG$Expected, CK_nb)> scot_BNG$EB_mm_local <- eb2$est *+ 100 SMREB_mmEB_mlEB_mm_local0100200300400500600 Moran sIDCluster provides a permutationbootstrap test for spatialautocorrelation of the differencebetween observed and expectedcounts:> lw <- nb2listw(CK_nb)> (20060614)> <- boot(as(scot_BNG,+ " "), statistic = ,+ R = 999, listw = lw, n = length(CK_nb),+ S0 = Szero(lw))Histogram of tt*Density 3 1123 of Standard Normalt*Moran sIIt also provides parametricbootstraps for variants, including theNegative Binomial.

6 > <- boot(as(scot_BNG,+ " "), statistic = ,+ sim = "parametric", = ,+ R = 999, listw = lw, n = length(CK_nb),+ S0 = Szero(lw))Histogram of tt*Density 3 1123 of Standard Normalt*Assun c ao and Reis correctionThe Assun c ao and Reis correction to Moran sIis implemented inspdep:> (scot_BNG$Observed, scot_BNG$Expected, lw, nsim = 999)Monte-Carlo simulation of Empirical Bayes Indexdata: cases: scot_BNG$Observed, risk population: scot_BNG$Expectedweights: lwnumber of simulations + 1: 1000statistic = , observed rank = 1000, p-value = hypothesis: greaterFitting base GLM modelsWe can fit GLMs for the base model with only the intercept, forthe Poisson, quasi-Poisson, and Negative Binomial, to give astarting point with respect to overdispersion:> <- glm(Observed ~ 1 + offset(log(Expected)), data = scot_BNG,+ family = poisson())> <- glm(Observed ~ 1 + offset(log(Expected)), data = scot_BNG,+ family = quasipoisson())> library(MASS)> <- (Observed ~ 1 + offset(log(Expected)), data = scot_BNG)> unlist(summary( )[20:21])NULLT ests for overdispersionTests for overdispersion, based in part on work by Dean, areprovided inDCluster:> ( , )Likelihood ratio test for overdispersiondata: : = , = 1, p-value < estimates:zscore > DeanB( )Dean's P_B test for overdispersiondata.

7 = , p-value < hypothesis: greaterFitting GLMsWe can augment the base model with the percentage occupied inagriculture, forestry and fisheries, as a measure of exposure tosunlight, but this does not seem to alleviate overdispersion much:> <- glm(Observed ~ PcAFF + offset(log(Expected)), data = scot_BNG,+ family = poisson())> <- glm(Observed ~ PcAFF + offset(log(Expected)), data = scot_BNG,+ family = quasipoisson())> <- (Observed ~ PcAFF + offset(log(Expected)), data = scot_BNG)> unlist(summary( )[20:21])NULL> anova( , )Likelihood ratio tests of Negative Binomial ModelsResponse: ObservedModel theta Resid. df 2 x log-lik. Test df1 1 + offset(log(Expected)) 55 PcAFF + offset(log(Expected)) 54 1 vs 2 1LR stat.

8 Pr(Chi)12 of the GLMsIn the same way that we stacked up smoothed rates, we can addthe standardised residulats of the GLM fits to ourSpatialobject,to examine them visually for patterning:> scot_BNG$base_glm_rst <- rstandard( )> scot_BNG$base_glmQ_rst <- rstandard( )> scot_BNG$base_nb_rst <- rstandard( )> scot_BNG$AFF_glm_rst <- rstandard( )> scot_BNG$AFF_glmQ_rst <- rstandard( )> scot_BNG$AFF_nb_rst <- rstandard( )Residuals of the GLMsbase_glm_rstbase_glmQ_rstbase_nb_rst AFF_glm_rstAFF_glmQ_rstAFF_nb_rst 8 6 4 202468 What next?IWe are pretty badly misspecified, but can the spatialdependence be separated from the distributional assumptions?IThere is a forthcoming paper inGeographical Analysisusingpermutation bootstrap on Moran sIof the deviance residualsof GLM fits, but this doesn t help with overdispersionIVirgilio G omez-Rubio has been working on exportingneighbour lists to Brugs/Openbugs and/or WinBUGS, and weare close to having something that can be releasedIWe can already exportSpatialPolygonsto WinBugs, buthere this would require further manual intervention to set linksto islandsExporting the resultsFinally, we write the results so far out as a shapefile, and as a textfile to be read by Mondrian, since dynamic graphics may offerfurther insight.

9 > writePolyShape(scot_BNG, "scot_BNG")> crs <- showWKT(proj4string(scot_BNG), " ")> cat(strwrap(gsub(",", ", ", crs)), sep = "\n")PROJCS["OSGB 1936 / British National Grid", GEOGCS["OSGB 1936",DATUM["D_OSGB_1936", SPHEROID["Airy_1830", , ]],PRIMEM["Greenwich", 0], UNIT["Degree", ]],PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin", 49],PARAMETER["central_meridian", -2], PARAMETER["scale_factor", ], PARAMETER["false_easting", 400000],PARAMETER["false_northing", -100000], UNIT["Meter", 1]]> sp2 Mondrian(scot_BNG, " ")Exploring the results in MondrianSumming upIUsing Spatial classes with shared methods for visualisation,import and export should let researchers get on with what youare good atIWe still need feedback on things that need improving, andcontributions of richer objects to suit the different researchdomainsIWriting small functions to output Spatial objects for othersoftware turns out to be much easier, since we now know howthe Spatial objects are constructedIPlease feel free to use the mailing list to follow up it ispretty likely that somebody else will have seen your problemalready (the archives are good too)