Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Applied Spatial Data Analysis with R Roger Bivand1
12 June 2014
1
Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway;
[email protected] Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Overview 1
2 3
4
Spatial data in R Study area Representing spatial data Coordinate reference systems Raster data Exporting data Handling spatial data in R: GRASS interface Handling spatial data in R: methods Topology predicates and operations Overlay Worked examples in spatial statistics Disease mapping Spatial autocorrelation Regression Bayesian spatial regression with INLA Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
R refresher
Many objects of interest in data analysis can be expressed as lists of numbers R sees the world this way too, and almost everything is expressed as vectors or lists of one kind or another R at its simplest behaves like an overgrown calculator, so that: > 2 + 2 [1] 4 > 2 * pi * 10 [1] 62.83185
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
R refresher Examples are made up of the R command line prompt > followed by the instruction in italics, and then the output The output of a vector printed to the screen is preceeded by the number in sequence of the vector value in square brackets; here the vector only has a length of one R lets us assign values to symbolic variables, so we can also do: > a a [1] 2 > a + a [1] 4 Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
R refresher Assignment does not print any output (usually), but we can see what is inside a symbolic variable (object) by entering its name R is built up using functions, some of which we can use on our symbolic variable (or object) a str() shows the structure of a in summary, it is numeric — here num, and has the value 2; length() tells us how long it is > str(a) num 2 > length(a) [1] 1
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
R refresher
class() tells us what class it belongs to
The values returned by functions can be assigned to new symbolic variables, here ca, which is a character vector > ca ca [1] "numeric" > str(ca) chr "numeric"
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
R refresher
R can express arithmetic and other operations directly, or through functions, which most often are vectorised. Let’s generate some data: > a b z a[1] + a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9] + + a[10] [1] 15 > sum(a) [1] 15 > sum(a^2) [1] 25 > sum(a * b) [1] 23
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
R refresher When we talk of arguments to a function, we mean the (possibly named) arguments passed to the function for evaluation We can access the arguments using the function help pages, or briefly using args(): > args(sum) function (..., na.rm = FALSE) NULL
Now we will look at why data-analytical software should be used to analyse data Problems with the data are very common, respondents who do not answer a particular question, equipment that fails occasionally, rain smudging field data collection sheets, etc., and the software needs to be able to handle them Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
R refresher As we see, sum() has an argument called na.rm, which is set to FALSE by default This means that NA — not available — values will not be removed by default, and the summation will fail, returning an NA, and alerting the analyst to problems with the data > a[3] a [1]
1
1 NA
1
1
2
2
2
2
2
> sum(a) [1] NA > sum(a, na.rm = TRUE) [1] 14 Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Olinda Hansen’s Disease data set The data are discussed in Lara, T. et al. (2001) Leprosy surveillance in Olinda, Brazil, using spatial analysis techniques. Cadernos de Sa´ ude P´ ublica, 17(5): 1153–1162. They use 243 1990 census districts, and include the 1991-1996 counts of new cases (CASES), 1993 year-end population figures (POP) and a deprivation index (DEPRIV) The data have been made available in cooperation with Marilia Sa Carvalho, and their use here is partly based on teaching notes that we wrote together some years ago The data now take the form of a shapefile, but were initially assembled from a Mapinfo file of census district borders and text files of variable values Roger Bivand
Applied Spatial Data Analysis with R
Olinda Hansen’s Disease data set So that we’ve got something to go on, let’s read in the data set, using an interface function in the rgdal package (packages are kept in libraries and extend R’s basic functions and classes): > library(rgdal)
The function as we will see later takes at least two arguments, a data source name, here the current directory, and a layer name, here the name of the shapefile without extension. Once we have checked the names in the imported object, we can use the spplot method in sp to make a map: > olinda names(olinda) [1] "AREA" [7] "SET"
"PERIMETER" "SETOR_" "CASES" "POP"
"SETOR_ID" "DEPRIV"
"VAR5"
"DENS_DEMO"
> spplot(olinda, "DEPRIV", col.regions = grey.colors(20, 0.9, 0.3))
Map of deprivation variable
0.8
0.6
0.4
0.2
0.0
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Data frames
Very often data are stored in a data frame with rows representing observations, here 1990 census districts in Olinda In data frames, the columns can be of different classes — in a matrix, all the columns have to have the same class The columns are accessed using the $ operator; the columns are list elements of the data frame We can also use square brackets to access data frame values, like look-up in a data base
Roger Bivand
Applied Spatial Data Analysis with R
Olinda data
> str(as(olinda, "data.frame"))
'data.frame': 243 obs. of 10 variables: $ AREA : num 79139 151153 265037 137696 121873 ... $ PERIMETER: num 1270 2189 2818 2098 3353 ... $ SETOR_ : num 2 3 4 5 6 7 8 9 10 11 ... $ SETOR_ID : num 1 2 3 4 5 6 9 10 11 12 ... $ VAR5 : int 242 224 223 229 228 222 218 225 227 221 ... $ DENS_DEMO: num 4380 10563 6642 13174 13812 ... $ SET : num 242 224 223 229 228 222 218 225 227 221 ... $ CASES : num 1 6 5 1 1 4 6 2 1 6 ... $ POP : num 337 1550 1711 1767 1638 ... $ DEPRIV : num 0.412 0.168 0.192 0.472 0.302 0.097 0.141 0.199 0.228 0.223 .
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Object framework To begin with, all contributed packages for handling spatial data in R had different representations of the data. This made it difficult to exchange data both within R between packages, and between R and external file formats and applications. The result has been an attempt to develop shared classes to represent spatial data in R, allowing some shared methods and many-to-one, one-to-many conversions. We chose to use new-style classes to represent spatial data, and are confident that this choice was justified.
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Spatial objects The foundation object is the Spatial class, with just two slots (new-style class objects have pre-defined components called slots) The first is a bounding box, and is mostly used for setting up plots The second is a CRS class object defining the coordinate reference system, and may be set to CRS(as.character(NA)), its default value. Operations on Spatial* objects should update or copy these values to the new Spatial* objects being created
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Coordinate reference systems The EPSG list among other sources is used in the workhorse PROJ.4 library, which as implemented by Frank Warmerdam handles transformation of spatial positions between different CRS This library is interfaced with R in the rgdal package, and the CRS class is defined partly in sp, partly in rgdal A CRS object is defined as a character NA string or a valid PROJ.4 CRS definition The validity of the definition can only be checked if rgdal is loaded
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Spatial points
The most basic spatial data object is a point, which may have 2 or 3 dimensions A single coordinate, or a set of such coordinates, may be used to define a SpatialPoints object; coordinates should be of mode double and will be promoted if not The points in a SpatialPoints object may be associated with a row of attributes to create a SpatialPointsDataFrame object The coordinates and attributes may, but do not have to be keyed to each other using ID values
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Spatial points classes and their slots
SpatialPointsDataFrame
Spatial
SpatialPoints
bbox
coords.nrs
proj4string
data SpatialPoints data.frame
coords Spatial
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Spatial*DataFrames
Spatial*DataFrames are constructed as Janus-like objects, looking to mapping and GIS people as “their” kind of object, as a collection of geometric freatures, with data associated with each feature. But to data analysts, the object “is” a data frame, because it behaves like one.
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Spatial lines and polygons A Line object is just a spaghetti collection of 2D coordinates; a Polygon object is a Line object with equal first and last coordinates A Lines object is a list of Line objects, such as all the contours at a single elevation; the same relationship holds between a Polygons object and a list of Polygon objects, such as islands belonging to the same county A comment is added to each Polygons object to indicate which interior ring belongs to which exterior ring (SFS) SpatialLines and SpatialPolygons objects are made using lists of Lines or Polygons objects respectively SpatialLinesDataFrame and SpatialPolygonsDataFrame objects are defined using SpatialLines and SpatialPolygons objects and standard data frames, and the ID fields are here required to match the data frame row names Roger Bivand
Applied Spatial Data Analysis with R
Spatial Polygons classes and slots SpatialLines
Lines
Line
lines
Lines
coords
Spatial
ID
SpatialPolygons
Polygons
Polygon
polygons
Polygons
labpt
plotOrder
plotOrder
area
Spatial
labpt ID
hole
area Spatial bbox proj4string
comment
ringDir coords
Back to Olinda So we can treat olinda as a data frame, for example adding an expected cases variable. The object also works within model fitting functions, like glm. > olinda$Expected str(olinda, max.level = 2) Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots ..@ data :'data.frame': 243 obs. of 11 variables: ..@ polygons :List of 243 .. .. [list output truncated] ..@ plotOrder : int [1:243] 243 56 39 231 224 60 58 232 145 66 ... ..@ bbox : num [1:2, 1:2] 288955 9111044 298446 9120528 .. ..- attr(*, "dimnames")=List of 2 ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots > str(as(olinda, "data.frame")) 'data.frame': 243 obs. of 11 variables: $ AREA : num 79139 151153 265037 137696 121873 ... $ PERIMETER: num 1270 2189 2818 2098 3353 ... $ SETOR_ : num 2 3 4 5 6 7 8 9 10 11 ... $ SETOR_ID : num 1 2 3 4 5 6 9 10 11 12 ... $ VAR5 : int 242 224 223 229 228 222 218 225 227 221 ... $ DENS_DEMO: num 4380 10563 6642 13174 13812 ... $ SET : num 242 224 223 229 228 222 218 225 227 221 ... $ CASES : num 1 6 5 1 1 4 6 2 1 6 ... $ POP : num 337 1550 1711 1767 1638 ... $ DEPRIV : num 0.412 0.168 0.192 0.472 0.302 0.097 0.141 0.199 0.228 0.223 ... $ Expected : num 1.4 6.43 7.1 7.33 6.8 ...
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Back to Olinda
The object also works within model fitting functions, like glm; note the number of rows. > str(model.frame(CASES ~ DEPRIV + offset(log(POP)), olinda), give.attr = FALSE) 'data.frame': 241 obs. of 3 variables: $ CASES : num 1 6 5 1 1 4 6 2 1 6 ... $ DEPRIV : num 0.412 0.168 0.192 0.472 0.302 0.097 0.141 0.199 0.228 0.223 ... $ offset(log(POP)): num 5.82 7.35 7.44 7.48 7.4 ...
Roger Bivand
Applied Spatial Data Analysis with R
The GridTopology class
The GridTopology class is the key constituent of raster representations, giving the coordinates of the south-west raster cell, the two cell sizes in the metric of the coordinates, giving the step to successive centres, and the numbers of cells in each dimension. We’ll return to examples after importing some data:
> library(sp) > getClass("GridTopology") Class "GridTopology" [package "sp"] Slots: Name: cellcentre.offset Class: numeric Name: Class:
cellsize numeric
Name: Class:
cells.dim integer
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Spatial grids and pixels There are two representations for data on regular rectangular grids (oriented N-S, E-W): SpatialPixels and SpatialGrid SpatialPixels are like SpatialPoints objects, but the coordinates have to be regularly spaced; the coordinates are stored, as are grid indices SpatialPixelsDataFrame objects only store attribute data
where it is present, but need to store the coordinates and grid indices of those grid cells SpatialGridDataFrame objects do not need to store
coordinates, because they fill the entire defined grid, but they need to store NA values where attribute values are missing Roger Bivand
Applied Spatial Data Analysis with R
Spatial grid and pixels classes and their slots
Spatial
SpatialGridDataFrame SpatialGrid
SpatialGrid
bbox
data
Spatial grid
proj4string GridTopology cellcentre.offset
SpatialPixelsDataFrame
cellsize
SpatialPixels
SpatialPixels
data
grid
cells.dim
grid.index
SpatialPoints
SpatialPoints
coords Spatial
data.frame
A SpatialGridDataFrame object The space shuttle flew a radar topography mission in 2000, giving 90m resolution elevation data for most of the world. The data here have been warped to a UTM projection, but for the WGS84 datum - we’ll see below how to project and if need be datum transform Spatial objects: > DEM summary(DEM) Object of class SpatialGridDataFrame Coordinates: min max x 288307.3 299442.7 y 9109607.8 9120835.2 Is projected: TRUE proj4string : [+proj=utm +zone=25 +south +ellps=WGS84 +units=m +no_defs] Grid attributes: cellcentre.offset cellsize cells.dim x 288353.4 92.02797 121 y 9109653.8 92.02797 122 Data attributes: Min. 1st Qu. Median Mean 3rd Qu. Max. -1.00 3.00 11.00 19.45 31.00 88.00 > is.na(DEM$band1) image(DEM, "band1", axes = TRUE, col = terrain.colors(20))
9110000
9112000
9114000
9116000
9118000
9120000
Elevation at Olinda
288000
290000
292000
294000
296000
298000
The original data set CRS Very often, a good deal of work is needed to find the actual coordinate reference system of a data set. Typically, researchers use what is to hand, without thinking that incomplete CRS specifications limit their ability to integrate other data. Here we were told that the map was from the 1990 census, that it was most likely UTM zone 25 south, but not its datum. Investigation showed that until very recently, the most used datum in Brazil has been Corrego Alegre: > proj4string(olinda) [1] NA > EPSG EPSG[grep("Corrego Alegre", EPSG$note), 1:2]
154 460 2741 2742 2743 2744 3140 3141 3142 3143 3144
code 4225 5524 5536 5537 5538 5539 22521 22522 22523 22524 22525
# # # # #
note # Corrego Alegre 1970-72 # Corrego Alegre 1961 # Corrego Alegre 1961 / UTM zone 21S # Corrego Alegre 1961 / UTM zone 22S # Corrego Alegre 1961 / UTM zone 23S # Corrego Alegre 1961 / UTM zone 24S Corrego Alegre 1970-72 / UTM zone 21S Corrego Alegre 1970-72 / UTM zone 22S Corrego Alegre 1970-72 / UTM zone 23S Corrego Alegre 1970-72 / UTM zone 24S Corrego Alegre 1970-72 / UTM zone 25S
The original data set CRS
In order to insert the parameters for the transformation to the WGS84 datum, we have to override the initial imported values: > set_ReplCRS_warn(FALSE) [1] FALSE > proj4string(olinda) proj4string(olinda)
[1] "+init=epsg:22525 +proj=utm +zone=25 +south +ellps=intl +towgs84=-206,172,-6,0,0,0,0 +units=m +no_defs > proj4string(DEM) [1] "+proj=utm +zone=25 +south +ellps=WGS84 +units=m +no_defs" > proj4string(DEM) proj4string(DEM) [1] "+init=epsg:31985 +proj=utm +zone=25 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" > set_ReplCRS_warn(TRUE) [1] TRUE
Getting olinda to WGS84
As we see, although both olinda and DEM are UTM zone 25 south, they differ in their ellipsoid and datum. Using spTransform methods in rgdal we can undertake a datum shift for olinda, making it possible to overplot in register: > olinda1 image(DEM, "band1", axes = TRUE, col = terrain.colors(20)) > plot(olinda1, add = TRUE, lwd = 0.5)
9110000
9112000
9114000
9116000
9118000
9120000
Getting olinda to WGS84
288000
290000
292000
294000
296000
298000
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
Reading rasters There are very many raster and image formats; some allow only one band of data, others think data bands are RGB, while yet others are flexible There is a simple readAsciiGrid function in maptools that reads ESRI Arc ASCII grids into SpatialGridDataFrame objects; it does not handle CRS and has a single band Much more support is available in rgdal in the readGDAL function, which — like readOGR — finds a usable driver if available and proceeds from there Using arguments to readGDAL, subregions or bands may be selected, which helps handle large rasters Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Study area Representing spatial data Coordinate reference systems Raster data Exporting data
But let’s keep things local Using Landsat 7 panchromatic imagery, we can get a similar contextual effect, here with histogram equalised 15m resolution: > pan > > > + > +
bb0 > > + > > >
bb0 + + + > > +
writeRAST6(DEM, "dem", flags = "o") execGRASS("g.region", rast = "dem") execGRASS("r.resamp.rst", input = "dem", ew_res = 14.25, ns_res = 14.25, elev = "DEM_resamp", slope = "slope", aspect = "aspect") execGRASS("g.region", rast = "DEM_resamp") DEM_out execGRASS("r.watershed", elevation = "DEM_resamp", + stream = "stream", threshold = 1000L, + convergence = 5L, memory = 300L) > execGRASS("r.thin", input = "stream", + output = "stream1", iterations = 200L)
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
9116000 9114000 9110000
> execGRASS("r.to.vect", input = "stream1", + output = "stream", feature = "line") > stream1 execGRASS("r.stream.extract", + elevation = "DEM_resamp", + stream_vect = "v_stream", + threshold = 1000) > stream1 library(rgeos)
predicates and operations for geometries.
> getScale()
The simplest are measures for planar
[1] 1e+08
geometries — only these are supported. A
> > > + >
complication is that computational
pols res0[as.integer(row.names(lnsInt)) + + 1] all.equal(res, res0) [1] TRUE
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Topology predicates and operations Overlay
Intersections - channels and tracts > system.time(tree tree[[2]]
of doing this increase with the number of
[1] 46 47 48 50
objects, as the predicate output matrix is
> res1 system.time(for (i in seq(along = res1)) { + if (!is.null(tree[[i]])) { + gi > + > > + + + + + +
buf50m + > > + > +
o_ndvi ch pm names(pm) [1] "raw" [4] "pmap"
"expCount" "relRisk"
> olinda2$SMR table(findInterval(pm$pmap, + seq(0, 1, 1/10))) 4 5 6 9 12 16
7 8 8 11
40
9 10 8 57 0
1 2 3 77 17 26
Frequency
20
If the underlying distribution of the data does not agree with our assumption, we may get several possible processes mixed up, overdispersion with spatial dependence:
60
80
Poisson probability map
0.0
0.2
0.4
0.6 pm$pmap
Roger Bivand
Applied Spatial Data Analysis with R
0.8
1.0
Some Empirical Bayes rates?
EB_ml
6
5
We can compare the relative rate and the EB relative rate: > library(DCluster) > olinda2$RR olinda2$EB_ml spplot(olinda2, c("RR", "EB_ml"), + col.regions = brewer.pal(10, + "RdBu"), at = c(seq(0, + 1, 0.25), seq(1, 6, + 1)))
4
RR
3
2
1
0
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Disease mapping Spatial autocorrelation Regression Bayesian spatial regression with INLA
Neighbours
●
In order to investigate spatial dependence, we need a list of neighbours. We take Queen contiguities of the tract polygons: > nb nb
●
●●
●
●
● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
●
● ● ● ●
● ● ● ●
● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
●● ●
●
●
●
●
●
●
●
● ● ●●
● ● ● ● ●
●
●
●
●
● ●
● ●
● ●
● ●
●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
Neighbour list object: Number of regions: 241 Number of nonzero links: 1324 Percentage nonzero weights: 2.279575 Average number of links: 5.493776
●
●
Roger Bivand
Applied Spatial Data Analysis with R
Spatial data in R Handling spatial data in R: GRASS interface Handling spatial data in R: methods Worked examples in spatial statistics
Disease mapping Spatial autocorrelation Regression Bayesian spatial regression with INLA
Local Empirical Bayes smoothing 6
EB_mm_local
If instead of shrinking to a global rate, we shrink to a local rate, we may be able to take unobserved heterogeneity into account; here we use the list of neighbours:
5
4
RR
EB_ml
3
2
> olinda2$Observed eb2 olinda2$EB_mm_local lw set.seed(130709) > moran.boot moran.pgboot EBImoran.mc(olinda2$CASES, olinda2$Expected, lw, nsim = 999) Monte-Carlo simulation of Empirical Bayes Index data: cases: olinda2$CASES, risk population: olinda2$Expected weights: lw number of simulations + 1: 1000 statistic = 0.3941, observed rank = 1000, p-value = 0.001 alternative hypothesis: greater
Roger Bivand
Applied Spatial Data Analysis with R
Maybe some modelling? To show how modelling and cartography can be inter-woven, let us fit some quasipoisson models, after seeing that — once again — spatial dependence and overdispersion seem to be entwined: > m0p DeanB(m0p) Dean's P_B test for overdispersion data: m0p P_B = 50.2478, p-value < 2.2e-16 alternative hypothesis: greater > > > > >
m0 > > > + + + + +
if (!("INLA" %in% .packages(all = TRUE))) source("http://www.math.ntnu.no/inla/givemeINLA.R") library(INLA) olinda2$INLA_ID
