Title: | Spatial Dependence: Weighting Schemes, Statistics |
---|---|
Description: | A collection of functions to create spatial weights matrix objects from polygon 'contiguities', from point patterns by distance and tessellations, for summarizing these objects, and for permitting their use in spatial data analysis, including regional aggregation by minimum spanning tree; a collection of tests for spatial 'autocorrelation', including global 'Morans I' and 'Gearys C' proposed by 'Cliff' and 'Ord' (1973, ISBN: 0850860369) and (1981, ISBN: 0850860814), 'Hubert/Mantel' general cross product statistic, Empirical Bayes estimates and 'Assunção/Reis' (1999) <doi:10.1002/(SICI)1097-0258(19990830)18:16%3C2147::AID-SIM179%3E3.0.CO;2-I> Index, 'Getis/Ord' G ('Getis' and 'Ord' 1992) <doi:10.1111/j.1538-4632.1992.tb00261.x> and multicoloured join count statistics, 'APLE' ('Li 'et al.' ) <doi:10.1111/j.1538-4632.2007.00708.x>, local 'Moran's I', 'Gearys C' ('Anselin' 1995) <doi:10.1111/j.1538-4632.1995.tb00338.x> and 'Getis/Ord' G ('Ord' and 'Getis' 1995) <doi:10.1111/j.1538-4632.1995.tb00912.x>, 'saddlepoint' approximations ('Tiefelsdorf' 2002) <doi:10.1111/j.1538-4632.2002.tb01084.x> and exact tests for global and local 'Moran's I' ('Bivand et al.' 2009) <doi:10.1016/j.csda.2008.07.021> and 'LOSH' local indicators of spatial heteroscedasticity ('Ord' and 'Getis') <doi:10.1007/s00168-011-0492-y>. The implementation of most of these measures is described in 'Bivand' and 'Wong' (2018) <doi:10.1007/s11749-018-0599-x>, with further extensions in 'Bivand' (2022) <doi:10.1111/gean.12319>. 'Lagrange' multiplier tests for spatial dependence in linear models are provided ('Anselin et al'. 1996) <doi:10.1016/0166-0462(95)02111-6>, as are 'Rao' score tests for hypothesised spatial 'Durbin' models based on linear models ('Koley' and 'Bera' 2023) <doi:10.1080/17421772.2023.2256810>. A local indicators for categorical data (LICD) implementation based on 'Carrer et al.' (2021) <doi:10.1016/j.jas.2020.105306> and 'Bivand et al.' (2017) <doi:10.1016/j.spasta.2017.03.003> was added in 1.3-7. From 'spdep' and 'spatialreg' versions >= 1.2-1, the model fitting functions previously present in this package are defunct in 'spdep' and may be found in 'spatialreg'. |
Authors: | Roger Bivand [cre, aut] , Micah Altman [ctb], Luc Anselin [ctb], Renato Assunção [ctb], Anil Bera [ctb], Olaf Berke [ctb], F. Guillaume Blanchet [ctb], Marilia Carvalho [ctb], Bjarke Christensen [ctb], Yongwan Chun [ctb], Carsten Dormann [ctb], Stéphane Dray [ctb], Dewey Dunnington [ctb] , Virgilio Gómez-Rubio [ctb], Malabika Koley [ctb], Tomasz Kossowski [ctb] , Elias Krainski [ctb], Pierre Legendre [ctb], Nicholas Lewin-Koh [ctb], Angela Li [ctb], Giovanni Millo [ctb], Werner Mueller [ctb], Hisaji Ono [ctb], Josiah Parry [ctb] , Pedro Peres-Neto [ctb], Michał Pietrzak [ctb] , Gianfranco Piras [ctb], Markus Reder [ctb], Jeff Sauer [ctb], Michael Tiefelsdorf [ctb], René Westerholt [ctb], Justyna Wilk [ctb] , Levi Wolf [ctb], Danlin Yu [ctb] |
Maintainer: | Roger Bivand <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.3-11 |
Built: | 2025-01-23 20:30:29 UTC |
Source: | https://github.com/r-spatial/spdep |
The method aggregates a spatial neighbours object, creating a new object listing the neighbours of the aggregates.
## S3 method for class 'nb' aggregate(x, IDs, remove.self = TRUE, ...)
## S3 method for class 'nb' aggregate(x, IDs, remove.self = TRUE, ...)
x |
an nb neighbour object |
IDs |
a character vector of IDs grouping the members of the neighbour object |
remove.self |
default TRUE: remove self-neighbours resulting from aggregation |
... |
unused - arguments passed through |
an nb neighbour object, with empty aggregates dropped.
Method suggested by Roberto Patuelli
Roger Bivand [email protected]
data(used.cars, package="spData") data(state) cont_st <- match(attr(usa48.nb, "region.id"), state.abb) cents <- as.matrix(as.data.frame(state.center))[cont_st,] opar <- par(mfrow=c(2,1)) plot(usa48.nb, cents, xlim=c(-125, -65), ylim=c(25, 50)) IDs <- as.character(state.division[cont_st]) agg_cents <- aggregate(cents, list(IDs), mean) agg_nb <- aggregate(usa48.nb, IDs) plot(agg_nb, agg_cents[, 2:3], xlim=c(-125, -65), ylim=c(25, 50)) text(agg_cents[, 2:3], agg_cents[, 1], cex=0.6) par(opar)
data(used.cars, package="spData") data(state) cont_st <- match(attr(usa48.nb, "region.id"), state.abb) cents <- as.matrix(as.data.frame(state.center))[cont_st,] opar <- par(mfrow=c(2,1)) plot(usa48.nb, cents, xlim=c(-125, -65), ylim=c(25, 50)) IDs <- as.character(state.division[cont_st]) agg_cents <- aggregate(cents, list(IDs), mean) agg_nb <- aggregate(usa48.nb, IDs) plot(agg_nb, agg_cents[, 2:3], xlim=c(-125, -65), ylim=c(25, 50)) text(agg_cents[, 2:3], agg_cents[, 1], cex=0.6) par(opar)
Measure a distance between two points on a plot using locator
; the function checks par("plt")
and par("usr")
to try to ensure that the aspect ratio y/x is 1, that is that the units of measurement in both x and y are equivalent.
airdist(ann=FALSE)
airdist(ann=FALSE)
ann |
annotate the plot with line measured and distance |
a list with members:
dist |
distance measured |
coords |
coordinates between which distance is measured |
Roger Bivand [email protected]
Calculates the autocovariate to be used in autonormal, autopoisson or autologistic regression. Three distance-weighting schemes are available.
autocov_dist(z, xy, nbs = 1, type = "inverse", zero.policy = NULL, style = "B", longlat=NULL)
autocov_dist(z, xy, nbs = 1, type = "inverse", zero.policy = NULL, style = "B", longlat=NULL)
z |
the response variable |
xy |
a matrix of coordinates or a SpatialPoints, |
nbs |
neighbourhood radius; default is 1 |
type |
the weighting scheme: "one" gives equal weight to all data points in the neighbourhood; "inverse" (the default) weights by inverse distance; "inverse.squared" weights by the square of "inverse" |
zero.policy |
default NULL, use global option value; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors |
style |
default “B” (changed from “W” 2015-01-27); style can take values “W”, “B”, “C”, “U”, and “S” |
longlat |
TRUE if point coordinates are longitude-latitude decimal, in which case distances are measured in kilometers; if xy is a SpatialPoints object, the value is taken from the object itself |
A numeric vector of autocovariate values
The validity of this approach strongly hinges on the correct choice of the
neighbourhood scheme! Using style="B" ensures symmetry of the
neighbourhood matrix (i.e. ). Please see Bardos et al. (2015)
for details.
Carsten F. Dormann and Roger Bivand
Augustin N.H., Mugglestone M.A. and Buckland S.T. (1996) An autologistic model for the spatial distribution of wildlife. Journal of Applied Ecology, 33, 339-347; Gumpertz M.L., Graham J.M. and Ristaino J.B. (1997) Autologistic model of spatial pattern of Phytophthora epidemic in bell pepper: effects of soil variables on disease presence. Journal of Agricultural, Biological and Environmental Statistics, 2, 131-156; Bardos, D.C., Guillera-Arroita, G. and Wintle, B.A. (2015) Valid auto-models for spatially autocorrelated occupancy and abundance data. arXiv, 1501.06529.
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) #xy <- cbind(columbus$X, columbus$Y) xy <- st_coordinates(st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)) ac1a <- autocov_dist(columbus$CRIME, xy, nbs=10, style="B", type="one") acinva <- autocov_dist(columbus$CRIME, xy, nbs=10, style="B", type="inverse") acinv2a <- autocov_dist(columbus$CRIME, xy, nbs=10, style="B", type="inverse.squared") plot(ac1a ~ columbus$CRIME, pch=16, ylim=c(0,9000)) points(acinva ~ columbus$CRIME, pch=16, col="red") points(acinv2a ~ columbus$CRIME, pch=16, col="blue") legend("topleft", legend=c("one", "inverse", "inverse.squared"), col=c("black", "red", "blue"), bty="n", pch=16) nb <- dnearneigh(xy, 0, 10) lw <- nb2listw(nb, style="B") ac1b <- lag(lw, columbus$CRIME) all.equal(ac1b, ac1a) nbd <- nbdists(nb, xy) gl <- lapply(nbd, function(x) 1/x) lw <- nb2listw(nb, glist=gl, style="B") acinvb <- lag(lw, columbus$CRIME) all.equal(acinvb, acinva) gl2 <- lapply(nbd, function(x) 1/(x^2)) lw <- nb2listw(nb, glist=gl2, style="B") acinv2b <- lag(lw, columbus$CRIME) all.equal(acinv2b, acinv2a) #xy <- SpatialPoints(xy) #acinva <- autocov_dist(columbus$CRIME, xy, nbs=10, style="W", # type="inverse") #nb <- dnearneigh(xy, 0, 10) #nbd <- nbdists(nb, xy) #gl <- lapply(nbd, function(x) 1/x) #lw <- nb2listw(nb, glist=gl) #acinvb <- lag(lw, columbus$CRIME) #all.equal(acinvb, acinva) acinvc <- autocov_dist(columbus$CRIME, st_centroid(st_geometry(columbus), of_largest_polygon=TRUE), nbs=10, style="W", type="inverse") all.equal(acinvc, acinva)
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) #xy <- cbind(columbus$X, columbus$Y) xy <- st_coordinates(st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)) ac1a <- autocov_dist(columbus$CRIME, xy, nbs=10, style="B", type="one") acinva <- autocov_dist(columbus$CRIME, xy, nbs=10, style="B", type="inverse") acinv2a <- autocov_dist(columbus$CRIME, xy, nbs=10, style="B", type="inverse.squared") plot(ac1a ~ columbus$CRIME, pch=16, ylim=c(0,9000)) points(acinva ~ columbus$CRIME, pch=16, col="red") points(acinv2a ~ columbus$CRIME, pch=16, col="blue") legend("topleft", legend=c("one", "inverse", "inverse.squared"), col=c("black", "red", "blue"), bty="n", pch=16) nb <- dnearneigh(xy, 0, 10) lw <- nb2listw(nb, style="B") ac1b <- lag(lw, columbus$CRIME) all.equal(ac1b, ac1a) nbd <- nbdists(nb, xy) gl <- lapply(nbd, function(x) 1/x) lw <- nb2listw(nb, glist=gl, style="B") acinvb <- lag(lw, columbus$CRIME) all.equal(acinvb, acinva) gl2 <- lapply(nbd, function(x) 1/(x^2)) lw <- nb2listw(nb, glist=gl2, style="B") acinv2b <- lag(lw, columbus$CRIME) all.equal(acinv2b, acinv2a) #xy <- SpatialPoints(xy) #acinva <- autocov_dist(columbus$CRIME, xy, nbs=10, style="W", # type="inverse") #nb <- dnearneigh(xy, 0, 10) #nbd <- nbdists(nb, xy) #gl <- lapply(nbd, function(x) 1/x) #lw <- nb2listw(nb, glist=gl) #acinvb <- lag(lw, columbus$CRIME) #all.equal(acinvb, acinva) acinvc <- autocov_dist(columbus$CRIME, st_centroid(st_geometry(columbus), of_largest_polygon=TRUE), nbs=10, style="W", type="inverse") all.equal(acinvc, acinva)
The data are collected inthe Atlas of condition indices published by the Joao Pinheiro Foundation and UNDP.
A shape polygon object with seven variables:
The identificator
Name of city
The population of city
Health Life Condition Index
Education Life Condition Index
Children Life Condition Index
Economic Life Condition Index
(GDAL37 <- as.numeric_version(unname(sf_extSoftVersion()["GDAL"])) >= "3.7.0") file <- "etc/shapes/bhicv.gpkg.zip" zipfile <- system.file(file, package="spdep") if (GDAL37) { bh <- st_read(zipfile) } else { td <- tempdir() bn <- sub(".zip", "", basename(file), fixed=TRUE) target <- unzip(zipfile, files=bn, exdir=td) bh <- st_read(target) }
(GDAL37 <- as.numeric_version(unname(sf_extSoftVersion()["GDAL"])) >= "3.7.0") file <- "etc/shapes/bhicv.gpkg.zip" zipfile <- system.file(file, package="spdep") if (GDAL37) { bh <- st_read(zipfile) } else { td <- tempdir() bn <- sub(".zip", "", basename(file), fixed=TRUE) target <- unzip(zipfile, files=bn, exdir=td) bh <- st_read(target) }
The function tallies the numbers of neighbours of regions in the neighbours list.
card(nb)
card(nb)
nb |
a neighbours list object of class |
“nb” objects are stored as lists of integer vectors, where the vectors contain either the indices in the range 1:n
for n
as length(nb)
of the neighbours of region i
, or as.integer(0)
to signal no neighbours. The function card(nb)
is used to extract the numbers of neighbours from the “nb” object.
An integer vector of the numbers of neighbours of regions in the neighbours list.
Roger Bivand [email protected]
Bivand R, Pebesma EJ, Gomez-Rubio V, (2008) Applied Spatial Data Analysis with R, Springer, New York, pp. 239-251; Bivand R, Portnov B, (2004) Exploring spatial data analysis techniques using R: the case of observations with no neighbours. In: Anselin L, Florax R, Rey S, (eds.), Advances in Spatial Econometrics, Methodology, Tools and Applications. Berlin: Springer-Verlag, pp. 121-142, doi:10.1007/978-3-662-05617-2_6.
col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) table(card(col.gal.nb))
col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) table(card(col.gal.nb))
The function generates a list of neighbours for a grid of cells. Helper functions are used to convert to and from the vector indices for row and column grid positions, and rook (shared edge) or queen (shared edge or vertex) neighbour definitions are applied by type. If torus is TRUE, the grid is mapped onto a torus, removing edge effects.
cell2nb(nrow, ncol, type="rook", torus=FALSE, legacy=FALSE, x=NULL) vi2mrc(i, nrow, ncol)
cell2nb(nrow, ncol, type="rook", torus=FALSE, legacy=FALSE, x=NULL) vi2mrc(i, nrow, ncol)
nrow |
number of rows in the grid, may also be an object inheriting from class |
ncol |
number of columns in the grid; if |
type |
default rook, may also be queen |
torus |
default FALSE, if TRUE map grid onto torus |
legacy |
default FALSE, nrow/ncol reversed, if TRUE wrong col/row directions (see https://github.com/r-spatial/spdep/issues/20) |
x |
if given, an object inheriting from class |
i |
vector of vector indices corresponding to rowcol, a matrix with two columns of row, column indices |
The function returns an object of class nb
with a list of integer vectors containing neighbour region number ids. See card
for details of “nb” objects.
Roger Bivand [email protected]
nb7rt <- cell2nb(7, 7) summary(nb7rt) xyc <- attr(nb7rt, "region.id") xy <- matrix(as.integer(unlist(strsplit(xyc, ":"))), ncol=2, byrow=TRUE) plot(nb7rt, xy) nb7rt <- cell2nb(7, 7, torus=TRUE) summary(nb7rt) run <- FALSE if (require("sp", quietly=TRUE)) run <- TRUE if (run) { # https://github.com/r-spatial/spdep/issues/20 GT <- GridTopology(c(1, 1), c(1, 1), c(10, 50)) SPix <- as(SpatialGrid(GT), "SpatialPixels") nb_rook_cont <- poly2nb(as(SPix, "SpatialPolygons"), queen=FALSE) nb_rook_dist <- dnearneigh(coordinates(SPix), 0, 1.01) all.equal(nb_rook_cont, nb_rook_dist, check.attributes=FALSE) ## [1] TRUE } if (run) { t.nb <- cell2nb(GT, type='rook', legacy=TRUE) isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] FALSE } if (run) { t.nb <- cell2nb(GT, type='rook') isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] TRUE } if (run) { # https://github.com/r-spatial/spdep/issues/55 # problem reported in issue caused by rep() cycling in unexpected order GT <- GridTopology(c(1, 1), c(1, 1), c(22, 11)) SPix <- as(SpatialGrid(GT), "SpatialPixels") nb_rook_cont <- poly2nb(as(SPix, "SpatialPolygons"), queen=FALSE) nb_rook_dist <- dnearneigh(coordinates(SPix), 0, 1.01) all.equal(nb_rook_cont, nb_rook_dist, check.attributes=FALSE) } if (run) { t.nb <- cell2nb(GT, type='rook', legacy=TRUE) isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] FALSE } if (run) { t.nb <- cell2nb(GT, type='rook', legacy=FALSE) isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] TRUE }
nb7rt <- cell2nb(7, 7) summary(nb7rt) xyc <- attr(nb7rt, "region.id") xy <- matrix(as.integer(unlist(strsplit(xyc, ":"))), ncol=2, byrow=TRUE) plot(nb7rt, xy) nb7rt <- cell2nb(7, 7, torus=TRUE) summary(nb7rt) run <- FALSE if (require("sp", quietly=TRUE)) run <- TRUE if (run) { # https://github.com/r-spatial/spdep/issues/20 GT <- GridTopology(c(1, 1), c(1, 1), c(10, 50)) SPix <- as(SpatialGrid(GT), "SpatialPixels") nb_rook_cont <- poly2nb(as(SPix, "SpatialPolygons"), queen=FALSE) nb_rook_dist <- dnearneigh(coordinates(SPix), 0, 1.01) all.equal(nb_rook_cont, nb_rook_dist, check.attributes=FALSE) ## [1] TRUE } if (run) { t.nb <- cell2nb(GT, type='rook', legacy=TRUE) isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] FALSE } if (run) { t.nb <- cell2nb(GT, type='rook') isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] TRUE } if (run) { # https://github.com/r-spatial/spdep/issues/55 # problem reported in issue caused by rep() cycling in unexpected order GT <- GridTopology(c(1, 1), c(1, 1), c(22, 11)) SPix <- as(SpatialGrid(GT), "SpatialPixels") nb_rook_cont <- poly2nb(as(SPix, "SpatialPolygons"), queen=FALSE) nb_rook_dist <- dnearneigh(coordinates(SPix), 0, 1.01) all.equal(nb_rook_cont, nb_rook_dist, check.attributes=FALSE) } if (run) { t.nb <- cell2nb(GT, type='rook', legacy=TRUE) isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] FALSE } if (run) { t.nb <- cell2nb(GT, type='rook', legacy=FALSE) isTRUE(all.equal(nb_rook_cont, t.nb, check.attributes=FALSE)) ## [1] TRUE }
Calculates Choynowski probability map values.
choynowski(n, x, row.names=NULL, tol = .Machine$double.eps^0.5, legacy=FALSE)
choynowski(n, x, row.names=NULL, tol = .Machine$double.eps^0.5, legacy=FALSE)
n |
a numeric vector of counts of cases |
x |
a numeric vector of populations at risk |
row.names |
row names passed through to output data frame |
tol |
accumulate values for observed counts >= expected until value less than tol |
legacy |
default FALSE using vectorised alternating side |
A data frame with columns:
pmap |
Poisson probability map values: probablility of getting a more “extreme” count than actually observed, one-tailed with less than expected and more than expected folded together |
type |
logical: TRUE if observed count less than expected |
Roger Bivand [email protected]
Choynowski, M (1959) Maps based on probabilities, Journal of the American Statistical Association, 54, 385–388; Cressie, N, Read, TRC (1985), Do sudden infant deaths come in clusters? Statistics and Decisions, Supplement Issue 2, 333–349; Bailey T, Gatrell A (1995) Interactive Spatial Data Analysis, Harlow: Longman, pp. 300–303.
auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE) auckland.nb <- poly2nb(auckland) res <- choynowski(auckland$M77_85, 9*auckland$Und5_81) resl <- choynowski(auckland$M77_85, 9*auckland$Und5_81, legacy=TRUE) all.equal(res, resl) rt <- sum(auckland$M77_85)/sum(9*auckland$Und5_81) ch_ppois_pmap <- numeric(length(auckland$Und5_81)) side <- c("greater", "less") for (i in seq(along=ch_ppois_pmap)) { ch_ppois_pmap[i] <- poisson.test(auckland$M77_85[i], r=rt, T=(9*auckland$Und5_81[i]), alternative=side[(res$type[i]+1)])$p.value } all.equal(ch_ppois_pmap, res$pmap) res1 <- probmap(auckland$M77_85, 9*auckland$Und5_81) table(abs(res$pmap - res1$pmap) < 0.00001, res$type) lt005 <- (res$pmap < 0.05) & (res$type) ge005 <- (res$pmap < 0.05) & (!res$type) cols <- rep("nonsig", length(lt005)) cols[lt005] <- "low" cols[ge005] <- "high" auckland$cols <- factor(cols) plot(auckland[,"cols"], main="Probability map")
auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE) auckland.nb <- poly2nb(auckland) res <- choynowski(auckland$M77_85, 9*auckland$Und5_81) resl <- choynowski(auckland$M77_85, 9*auckland$Und5_81, legacy=TRUE) all.equal(res, resl) rt <- sum(auckland$M77_85)/sum(9*auckland$Und5_81) ch_ppois_pmap <- numeric(length(auckland$Und5_81)) side <- c("greater", "less") for (i in seq(along=ch_ppois_pmap)) { ch_ppois_pmap[i] <- poisson.test(auckland$M77_85[i], r=rt, T=(9*auckland$Und5_81[i]), alternative=side[(res$type[i]+1)])$p.value } all.equal(ch_ppois_pmap, res$pmap) res1 <- probmap(auckland$M77_85, 9*auckland$Und5_81) table(abs(res$pmap - res1$pmap) < 0.00001, res$type) lt005 <- (res$pmap < 0.05) & (res$type) ge005 <- (res$pmap < 0.05) & (!res$type) cols <- rep("nonsig", length(lt005)) cols[lt005] <- "low" cols[ge005] <- "high" auckland$cols <- factor(cols) plot(auckland[,"cols"], main="Probability map")
The data set is now part of the spData package
data(columbus)
data(columbus)
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
The function finds symmetric differences between lists of neighbours (ordering of objects does not matter), returning a nb
neighbour list of those found
diffnb(x, y, verbose=NULL, legacy=TRUE)
diffnb(x, y, verbose=NULL, legacy=TRUE)
x |
an object of class |
y |
an object of class |
verbose |
default NULL, use global option value; report regions ids taken from object attribute "region.id" with differences, ignored when |
legacy |
default TRUE; use legacy code; if FALSE use differences between sparse matrix representations |
A neighbours list with class nb
Roger Bivand [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE) rn <- row.names(columbus) knn1 <- knearneigh(coords, 1) knn2 <- knearneigh(coords, 2) nb1 <- knn2nb(knn1, row.names=rn) nb2 <- knn2nb(knn2, row.names=rn) diffs <- diffnb(nb2, nb1, legacy=TRUE) diffsl <- diffnb(nb2, nb1, legacy=FALSE) all.equal(diffs, diffsl) # call attribute fifference expected diffsd <- union.nb(setdiff.nb(nb1, nb2), setdiff.nb(nb2, nb1)) all.equal(diffs, diffsd) # call attribute fifference expected opar <- par(no.readonly=TRUE) plot(st_geometry(columbus), border="grey", reset=FALSE, main="Plot of first (black) and second (red)\nnearest neighbours") plot(nb1, coords, add=TRUE) plot(diffs, coords, add=TRUE, col="red", lty=2) par(opar)
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE) rn <- row.names(columbus) knn1 <- knearneigh(coords, 1) knn2 <- knearneigh(coords, 2) nb1 <- knn2nb(knn1, row.names=rn) nb2 <- knn2nb(knn2, row.names=rn) diffs <- diffnb(nb2, nb1, legacy=TRUE) diffsl <- diffnb(nb2, nb1, legacy=FALSE) all.equal(diffs, diffsl) # call attribute fifference expected diffsd <- union.nb(setdiff.nb(nb1, nb2), setdiff.nb(nb2, nb1)) all.equal(diffs, diffsd) # call attribute fifference expected opar <- par(no.readonly=TRUE) plot(st_geometry(columbus), border="grey", reset=FALSE, main="Plot of first (black) and second (red)\nnearest neighbours") plot(nb1, coords, add=TRUE) plot(diffs, coords, add=TRUE, col="red", lty=2) par(opar)
The function identifies neighbours of region points by Euclidean distance in the metric of the points between lower (greater than or equal to (changed from version 1.1-7)) and upper (less than or equal to) bounds, or with longlat = TRUE, by Great Circle distance in kilometers. If x
is an "sf"
object and use_s2=
is TRUE
, spherical distances in km are used.
dnearneigh(x, d1, d2, row.names = NULL, longlat = NULL, bounds=c("GE", "LE"), use_kd_tree=TRUE, symtest=FALSE, use_s2=packageVersion("s2") > "1.0.7", k=200, dwithin=TRUE)
dnearneigh(x, d1, d2, row.names = NULL, longlat = NULL, bounds=c("GE", "LE"), use_kd_tree=TRUE, symtest=FALSE, use_s2=packageVersion("s2") > "1.0.7", k=200, dwithin=TRUE)
x |
matrix of point coordinates, an object inheriting from SpatialPoints or an |
d1 |
lower distance bound in the metric of the points if planar coordinates, in km if in geographical coordinates |
d2 |
upper distance boundd in the metric of the points if planar coordinates, in km if in geographical coordinates |
row.names |
character vector of region ids to be added to the neighbours list as attribute |
longlat |
TRUE if point coordinates are geographical longitude-latitude decimal degrees, in which case distances are measured in kilometers; if x is a SpatialPoints object, the value is taken from the object itself, and overrides this argument if not NULL |
bounds |
character vector of length 2, default |
use_kd_tree |
default TRUE, if TRUE, use dbscan |
symtest |
Default FALSE; before release 1.1-7, TRUE - run symmetry check on output object, costly with large numbers of points. |
use_s2 |
default= |
k |
default 200, the number of closest points to consider when searching when using |
dwithin |
default TRUE, if FALSE, use |
The function returns a list of integer vectors giving the region id numbers
for neighbours satisfying the distance criteria. See card
for details of “nb” objects.
Roger Bivand [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE) rn <- row.names(columbus) k1 <- knn2nb(knearneigh(coords)) all.linked <- max(unlist(nbdists(k1, coords))) col.nb.0.all <- dnearneigh(coords, 0, all.linked, row.names=rn) summary(col.nb.0.all, coords) opar <- par(no.readonly=TRUE) plot(st_geometry(columbus), border="grey", reset=FALSE, main=paste("Distance based neighbours 0-", format(all.linked), sep="")) plot(col.nb.0.all, coords, add=TRUE) par(opar) (sfc_obj <- st_centroid(st_geometry(columbus))) col.nb.0.all_sf <- dnearneigh(sfc_obj, 0, all.linked, row.names=rn) all.equal(col.nb.0.all, col.nb.0.all_sf, check.attributes=FALSE) data(state) us48.fipsno <- read.geoda(system.file("etc/weights/us48.txt", package="spdep")[1]) if (as.numeric(paste(version$major, version$minor, sep="")) < 19) { m50.48 <- match(us48.fipsno$"State.name", state.name) } else { m50.48 <- match(us48.fipsno$"State_name", state.name) } xy <- as.matrix(as.data.frame(state.center))[m50.48,] llk1 <- knn2nb(knearneigh(xy, k=1, longlat=FALSE)) (all.linked <- max(unlist(nbdists(llk1, xy, longlat=FALSE)))) ll.nb <- dnearneigh(xy, 0, all.linked, longlat=FALSE) summary(ll.nb, xy, longlat=TRUE, scale=0.5) gck1 <- knn2nb(knearneigh(xy, k=1, longlat=TRUE)) (all.linked <- max(unlist(nbdists(gck1, xy, longlat=TRUE)))) gc.nb <- dnearneigh(xy, 0, all.linked, longlat=TRUE) summary(gc.nb, xy, longlat=TRUE, scale=0.5) plot(ll.nb, xy) plot(diffnb(ll.nb, gc.nb), xy, add=TRUE, col="red", lty=2) title(main="Differences Euclidean/Great Circle") #xy1 <- SpatialPoints((as.data.frame(state.center))[m50.48,], # proj4string=CRS("+proj=longlat +ellps=GRS80")) #gck1a <- knn2nb(knearneigh(xy1, k=1)) #(all.linked <- max(unlist(nbdists(gck1a, xy1)))) #gc.nb <- dnearneigh(xy1, 0, all.linked) #summary(gc.nb, xy1, scale=0.5) xy1 <- st_as_sf((as.data.frame(state.center))[m50.48,], coords=1:2, crs=st_crs("OGC:CRS84")) old_use_s2 <- sf_use_s2() sf_use_s2(TRUE) gck1b <- knn2nb(knearneigh(xy1, k=1)) system.time(o <- nbdists(gck1b, xy1)) (all.linked <- max(unlist(o))) # use s2 brute-force dwithin_matrix approach for s2 <= 1.0.7 system.time(gc.nb.dwithin <- dnearneigh(xy1, 0, all.linked, use_s2=TRUE, dwithin=TRUE)) summary(gc.nb, xy1, scale=0.5) # use s2 closest_edges approach s2 > 1.0.7 if (packageVersion("s2") > "1.0.7") { (system.time(gc.nb.closest <- dnearneigh(xy1, 0, all.linked, dwithin=FALSE))) } if (packageVersion("s2") > "1.0.7") { system.time(gc.nb.dwithin <- dnearneigh(xy1, 0, all.linked, use_s2=TRUE, dwithin=TRUE)) } if (packageVersion("s2") > "1.0.7") { summary(gc.nb.dwithin, xy1, scale=0.5) } if (packageVersion("s2") > "1.0.7") { summary(gc.nb.closest, xy1, scale=0.5) } # use legacy symmetric brute-force approach system.time(gc.nb.legacy <- dnearneigh(xy1, 0, all.linked, use_s2=FALSE)) summary(gc.nb, xy1, scale=0.5) if (packageVersion("s2") > "1.0.7") all.equal(gc.nb.closest, gc.nb.dwithin, check.attributes=FALSE) # legacy is ellipsoidal, s2 spherical, so minor differences expected if (packageVersion("s2") > "1.0.7") all.equal(gc.nb, gc.nb.closest, check.attributes=FALSE) all.equal(gc.nb, gc.nb.dwithin, check.attributes=FALSE) sf_use_s2(old_use_s2) # example of reading points with readr::read_csv() yielding a tibble load(system.file("etc/misc/coords.rda", package="spdep")) class(coords) k1 <- knn2nb(knearneigh(coords, k=1)) all.linked <- max(unlist(nbdists(k1, coords))) dnearneigh(coords, 0, all.linked)
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE) rn <- row.names(columbus) k1 <- knn2nb(knearneigh(coords)) all.linked <- max(unlist(nbdists(k1, coords))) col.nb.0.all <- dnearneigh(coords, 0, all.linked, row.names=rn) summary(col.nb.0.all, coords) opar <- par(no.readonly=TRUE) plot(st_geometry(columbus), border="grey", reset=FALSE, main=paste("Distance based neighbours 0-", format(all.linked), sep="")) plot(col.nb.0.all, coords, add=TRUE) par(opar) (sfc_obj <- st_centroid(st_geometry(columbus))) col.nb.0.all_sf <- dnearneigh(sfc_obj, 0, all.linked, row.names=rn) all.equal(col.nb.0.all, col.nb.0.all_sf, check.attributes=FALSE) data(state) us48.fipsno <- read.geoda(system.file("etc/weights/us48.txt", package="spdep")[1]) if (as.numeric(paste(version$major, version$minor, sep="")) < 19) { m50.48 <- match(us48.fipsno$"State.name", state.name) } else { m50.48 <- match(us48.fipsno$"State_name", state.name) } xy <- as.matrix(as.data.frame(state.center))[m50.48,] llk1 <- knn2nb(knearneigh(xy, k=1, longlat=FALSE)) (all.linked <- max(unlist(nbdists(llk1, xy, longlat=FALSE)))) ll.nb <- dnearneigh(xy, 0, all.linked, longlat=FALSE) summary(ll.nb, xy, longlat=TRUE, scale=0.5) gck1 <- knn2nb(knearneigh(xy, k=1, longlat=TRUE)) (all.linked <- max(unlist(nbdists(gck1, xy, longlat=TRUE)))) gc.nb <- dnearneigh(xy, 0, all.linked, longlat=TRUE) summary(gc.nb, xy, longlat=TRUE, scale=0.5) plot(ll.nb, xy) plot(diffnb(ll.nb, gc.nb), xy, add=TRUE, col="red", lty=2) title(main="Differences Euclidean/Great Circle") #xy1 <- SpatialPoints((as.data.frame(state.center))[m50.48,], # proj4string=CRS("+proj=longlat +ellps=GRS80")) #gck1a <- knn2nb(knearneigh(xy1, k=1)) #(all.linked <- max(unlist(nbdists(gck1a, xy1)))) #gc.nb <- dnearneigh(xy1, 0, all.linked) #summary(gc.nb, xy1, scale=0.5) xy1 <- st_as_sf((as.data.frame(state.center))[m50.48,], coords=1:2, crs=st_crs("OGC:CRS84")) old_use_s2 <- sf_use_s2() sf_use_s2(TRUE) gck1b <- knn2nb(knearneigh(xy1, k=1)) system.time(o <- nbdists(gck1b, xy1)) (all.linked <- max(unlist(o))) # use s2 brute-force dwithin_matrix approach for s2 <= 1.0.7 system.time(gc.nb.dwithin <- dnearneigh(xy1, 0, all.linked, use_s2=TRUE, dwithin=TRUE)) summary(gc.nb, xy1, scale=0.5) # use s2 closest_edges approach s2 > 1.0.7 if (packageVersion("s2") > "1.0.7") { (system.time(gc.nb.closest <- dnearneigh(xy1, 0, all.linked, dwithin=FALSE))) } if (packageVersion("s2") > "1.0.7") { system.time(gc.nb.dwithin <- dnearneigh(xy1, 0, all.linked, use_s2=TRUE, dwithin=TRUE)) } if (packageVersion("s2") > "1.0.7") { summary(gc.nb.dwithin, xy1, scale=0.5) } if (packageVersion("s2") > "1.0.7") { summary(gc.nb.closest, xy1, scale=0.5) } # use legacy symmetric brute-force approach system.time(gc.nb.legacy <- dnearneigh(xy1, 0, all.linked, use_s2=FALSE)) summary(gc.nb, xy1, scale=0.5) if (packageVersion("s2") > "1.0.7") all.equal(gc.nb.closest, gc.nb.dwithin, check.attributes=FALSE) # legacy is ellipsoidal, s2 spherical, so minor differences expected if (packageVersion("s2") > "1.0.7") all.equal(gc.nb, gc.nb.closest, check.attributes=FALSE) all.equal(gc.nb, gc.nb.dwithin, check.attributes=FALSE) sf_use_s2(old_use_s2) # example of reading points with readr::read_csv() yielding a tibble load(system.file("etc/misc/coords.rda", package="spdep")) class(coords) k1 <- knn2nb(knearneigh(coords, k=1)) all.linked <- max(unlist(nbdists(k1, coords))) dnearneigh(coords, 0, all.linked)
droplinks
drops links to and from or just to a region from a neighbours list. The example corresponds to Fingleton's Table 1, (1999) p. 6, for lattices 5 to 19. addlinks1
adds links from a single region to specified regions.
droplinks(nb, drop, sym=TRUE) addlinks1(nb, from, to, sym=TRUE)
droplinks(nb, drop, sym=TRUE) addlinks1(nb, from, to, sym=TRUE)
nb |
a neighbours list object of class |
drop |
either a logical vector the length of |
sym |
TRUE for removal of both "row" and "column" links, FALSE for only "row" links; when adding links, inserts links to the from region from the to regions |
from |
single from region for adding links, either a character vector of length 1 of the named from region corresponding to |
to |
to regions, either a character vector of named from regions corresponding to |
The function returns an object of class nb
with a list of integer vectors containing neighbour region number ids.
Roger Bivand [email protected]
B. Fingleton (1999) Spurious spatial regression: some Monte Carlo results with a spatial unit root and spatial cointegration, Journal of Regional Science 39, pp. 1–19.
rho <- c(0.2, 0.5, 0.95, 0.999, 1.0) ns <- c(5, 7, 9, 11, 13, 15, 17, 19) mns <- matrix(0, nrow=length(ns), ncol=length(rho)) rownames(mns) <- ns colnames(mns) <- rho mxs <- matrix(0, nrow=length(ns), ncol=length(rho)) rownames(mxs) <- ns colnames(mxs) <- rho for (i in 1:length(ns)) { nblist <- cell2nb(ns[i], ns[i]) nbdropped <- droplinks(nblist, ((ns[i]*ns[i])+1)/2, sym=FALSE) listw <- nb2listw(nbdropped, style="W", zero.policy=TRUE) wmat <- listw2mat(listw) for (j in 1:length(rho)) { mat <- diag(ns[i]*ns[i]) - rho[j] * wmat res <- diag(solve(t(mat) %*% mat)) mns[i,j] <- mean(res) mxs[i,j] <- max(res) } } print(mns) print(mxs)
rho <- c(0.2, 0.5, 0.95, 0.999, 1.0) ns <- c(5, 7, 9, 11, 13, 15, 17, 19) mns <- matrix(0, nrow=length(ns), ncol=length(rho)) rownames(mns) <- ns colnames(mns) <- rho mxs <- matrix(0, nrow=length(ns), ncol=length(rho)) rownames(mxs) <- ns colnames(mxs) <- rho for (i in 1:length(ns)) { nblist <- cell2nb(ns[i], ns[i]) nbdropped <- droplinks(nblist, ((ns[i]*ns[i])+1)/2, sym=FALSE) listw <- nb2listw(nbdropped, style="W", zero.policy=TRUE) wmat <- listw2mat(listw) for (j in 1:length(rho)) { mat <- diag(ns[i]*ns[i]) - rho[j] * wmat res <- diag(solve(t(mat) %*% mat)) mns[i,j] <- mean(res) mxs[i,j] <- max(res) } } print(mns) print(mxs)
The function computes global empirical Bayes estimates for rates "shrunk" to the overall mean.
EBest(n, x, family="poisson")
EBest(n, x, family="poisson")
n |
a numeric vector of counts of cases |
x |
a numeric vector of populations at risk |
family |
either "poisson" for rare conditions or "binomial" for non-rare conditions |
Details of the implementation for the "poisson" family are to be found in Marshall, p. 284–5, and Bailey and Gatrell p. 303–306 and exercise 8.2, pp. 328–330. For the "binomial" family, see Martuzzi and Elliott (implementation by Olaf Berke).
A data frame with two columns:
raw |
a numerical vector of raw (crude) rates |
estmm |
a numerical vector of empirical Bayes estimates |
and a parameters
attribute list with components:
a |
global method of moments phi value |
m |
global method of moments gamma value |
Roger Bivand [email protected] and Olaf Berke, Population Medicine, OVC, University of Guelph, CANADA
Marshall R M (1991) Mapping disease and mortality rates using Empirical Bayes Estimators, Applied Statistics, 40, 283–294; Bailey T, Gatrell A (1995) Interactive Spatial Data Analysis, Harlow: Longman, pp. 303–306, Martuzzi M, Elliott P (1996) Empirical Bayes estimation of small area prevalence of non-rare conditions, Statistics in Medicine 15, 1867–1873.
auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE) res <- EBest(auckland$M77_85, 9*auckland$Und5_81) attr(res, "parameters") auckland$estmm000 <- res$estmm*1000 plot(auckland[,"estmm000"], breaks=c(0,2,2.5,3,3.5,5), main="Infant mortality per 1000 per year") data(huddersfield, package="spData") res <- EBest(huddersfield$cases, huddersfield$total, family="binomial") round(res[,1:2],4)*100
auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE) res <- EBest(auckland$M77_85, 9*auckland$Und5_81) attr(res, "parameters") auckland$estmm000 <- res$estmm*1000 plot(auckland[,"estmm000"], breaks=c(0,2,2.5,3,3.5,5), main="Infant mortality per 1000 per year") data(huddersfield, package="spData") res <- EBest(huddersfield$cases, huddersfield$total, family="binomial") round(res[,1:2],4)*100
An empirical Bayes index modification of Moran's I for testing for spatial autocorrelation in a rate, typically the number of observed cases in a population at risk. The index value is tested by using nsim random permutations of the index for the given spatial weighting scheme, to establish the rank of the observed statistic in relation to the nsim simulated values.
EBImoran.mc(n, x, listw, nsim, zero.policy = attr(listw, "zero.policy"), alternative = "greater", spChk=NULL, return_boot=FALSE, subtract_mean_in_numerator=TRUE)
EBImoran.mc(n, x, listw, nsim, zero.policy = attr(listw, "zero.policy"), alternative = "greater", spChk=NULL, return_boot=FALSE, subtract_mean_in_numerator=TRUE)
n |
a numeric vector of counts of cases the same length as the neighbours list in listw |
x |
a numeric vector of populations at risk the same length as the neighbours list in listw |
listw |
a |
nsim |
number of permutations |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "two.sided", or "less" |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
return_boot |
return an object of class |
subtract_mean_in_numerator |
default TRUE, if TRUE subtract mean of z in numerator of EBI equation on p. 2157 in reference (consulted with Renato Assunção 2016-02-19); until February 2016 the default was FALSE agreeing with the printed paper. |
The statistic used is (m is the number of observations):
where:
and:
A list with class htest
and mc.sim
containing the
following components:
statistic |
the value of the observed Moran's I. |
parameter |
the rank of the observed Moran's I. |
p.value |
the pseudo p-value of the test. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data, and the number of simulations. |
res |
nsim simulated values of statistic, final value is observed statistic |
z |
a numerical vector of Empirical Bayes indices as z above |
Roger Bivand [email protected]
Assunção RM, Reis EA 1999 A new proposal to adjust Moran's I for population density. Statistics in Medicine 18, pp. 2147–2162; Bivand RS, Wong DWS 2018 Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748 doi:10.1007/s11749-018-0599-x
nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE) rn <- as.character(nc.sids$FIPS) ncCC89_nb <- read.gal(system.file("weights/ncCC89.gal", package="spData")[1], region.id=rn) EBImoran.mc(nc.sids$SID74, nc.sids$BIR74, nb2listw(ncCC89_nb, style="B", zero.policy=TRUE), nsim=999, alternative="two.sided", zero.policy=TRUE) sids.p <- nc.sids$SID74 / nc.sids$BIR74 moran.mc(sids.p, nb2listw(ncCC89_nb, style="B", zero.policy=TRUE), nsim=999, alternative="two.sided", zero.policy=TRUE)
nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE) rn <- as.character(nc.sids$FIPS) ncCC89_nb <- read.gal(system.file("weights/ncCC89.gal", package="spData")[1], region.id=rn) EBImoran.mc(nc.sids$SID74, nc.sids$BIR74, nb2listw(ncCC89_nb, style="B", zero.policy=TRUE), nsim=999, alternative="two.sided", zero.policy=TRUE) sids.p <- nc.sids$SID74 / nc.sids$BIR74 moran.mc(sids.p, nb2listw(ncCC89_nb, style="B", zero.policy=TRUE), nsim=999, alternative="two.sided", zero.policy=TRUE)
The function computes local empirical Bayes estimates for rates "shrunk" to a neighbourhood mean for neighbourhoods given by the nb
neighbourhood list.
EBlocal(ri, ni, nb, zero.policy = NULL, spChk = NULL, geoda=FALSE)
EBlocal(ri, ni, nb, zero.policy = NULL, spChk = NULL, geoda=FALSE)
ri |
a numeric vector of counts of cases the same length as the neighbours list in nb; if there are many zero counts, some estimates may be affected by division by zero, see https://stat.ethz.ch/pipermail/r-sig-geo/2022-January/028882.html |
ni |
a numeric vector of populations at risk the same length as the neighbours list in nb |
nb |
a |
zero.policy |
default NULL, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
geoda |
default=FALSE, following Marshall's algorithm as interpreted by Bailey and Gatrell, pp. 305-307, and exercise 8.2, pp. 328-330 for the definition of phi; TRUE for the definition of phi used in GeoDa (see discussion on OpenSpace mailing list June 2003: http://agec221.agecon.uiuc.edu/pipermail/openspace/2003-June/thread.html) |
Details of the implementation are to be found in Marshall, p. 286, and Bailey and Gatrell p. 307 and exercise 8.2, pp. 328–330. The example results do not fully correspond to the sources because of slightly differing neighbourhoods, but are generally close. If there are many zero counts, some estimates may be affected by division by zero, see https://stat.ethz.ch/pipermail/r-sig-geo/2022-January/028882.html.
A data frame with two columns:
raw |
a numerical vector of raw (crude) rates |
est |
a numerical vector of local empirical Bayes estimates |
and a parameters
attribute list with components (if both are zero, the estimate will be NaN
, https://stat.ethz.ch/pipermail/r-sig-geo/2022-January/028882.html):
a |
a numerical vector of local phi values |
m |
a numerical vector of local gamma values |
Roger Bivand [email protected], based on contributions by Marilia Carvalho
Marshall R M (1991) Mapping disease and mortality rates using Empirical Bayes Estimators, Applied Statistics, 40, 283–294; Bailey T, Gatrell A (1995) Interactive Spatial Data Analysis, Harlow: Longman, pp. 303–306.
auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE) auckland.nb <- poly2nb(auckland) res <- EBlocal(auckland$M77_85, 9*auckland$Und5_81, auckland.nb) auckland$est000 <- res$est*1000 plot(auckland[,"est000"], breaks=c(0,2,2.5,3,3.5,8), main="Infant mortality per 1000 per year")
auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE) auckland.nb <- poly2nb(auckland) res <- EBlocal(auckland$M77_85, 9*auckland$Und5_81, auckland.nb) auckland$est000 <- res$est*1000 plot(auckland[,"est000"], breaks=c(0,2,2.5,3,3.5,8), main="Infant mortality per 1000 per year")
The function provides simple interactive editing of neighbours lists to
allow unneeded links to be deleted, and missing links to be inserted. It uses identify
to pick the
endpoints of the link to be deleted or added, and asks for confirmation before
committing. If the result is not assigned to a new object, the editing
will be lost - as in edit
.
This method relies on direct contact with the graphics device. Do not use in RStudio.
## S3 method for class 'nb' edit(name, coords, polys=NULL, ..., use_region.id=FALSE)
## S3 method for class 'nb' edit(name, coords, polys=NULL, ..., use_region.id=FALSE)
name |
an object of class |
coords |
matrix of region point coordinates; if missing and polys= inherits from |
polys |
if polygon boundaries supplied, will be used as background; must inherit from |
... |
further arguments passed to or from other methods |
use_region.id |
default |
The function returns an object of class nb
with the edited list of integer vectors containing neighbour region number ids, with added attributes tallying the added and deleted links.
Roger Bivand [email protected]
## Not run: columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) class(columbus) if (FALSE) nnb1 <- edit.nb(col.gal.nb, polys=as(columbus, "Spatial")) ## End(Not run)
## Not run: columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) class(columbus) if (FALSE) nnb1 <- edit.nb(col.gal.nb, polys=as(columbus, "Spatial")) ## End(Not run)
The data set is now part of the spData package
data(eire)
data(eire)
A simple function to compute Geary's C, called by geary.test
and geary.mc
;
geary.intern
is an internal function used to vary the similarity
criterion.
geary(x, listw, n, n1, S0, zero.policy=attr(listw, "zero.policy"), scale=TRUE)
geary(x, listw, n, n1, S0, zero.policy=attr(listw, "zero.policy"), scale=TRUE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
n |
number of zones |
n1 |
n - 1 |
S0 |
global sum of weights |
zero.policy |
default |
scale |
default TRUE, may be FALSE to revert changes made to accommodate |
a list with
C |
Geary's C |
K |
sample kurtosis of x |
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 17.
geary.test
, geary.mc
,
sp.mantel.mc
data(oldcol) col.W <- nb2listw(COL.nb, style="W") str(geary(COL.OLD$CRIME, col.W, length(COL.nb), length(COL.nb)-1, Szero(col.W)))
data(oldcol) col.W <- nb2listw(COL.nb, style="W") str(geary(COL.OLD$CRIME, col.W, length(COL.nb), length(COL.nb)-1, Szero(col.W)))
A permutation test for Geary's C statistic calculated by using nsim random permutations of x for the given spatial weighting scheme, to establish the rank of the observed statistic in relation to the nsim simulated values.
geary.mc(x, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, adjust.n=TRUE, return_boot=FALSE, na.action=na.fail, scale=TRUE)
geary.mc(x, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, adjust.n=TRUE, return_boot=FALSE, na.action=na.fail, scale=TRUE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
nsim |
number of permutations |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), or "less"; this reversal corresponds to that on |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted |
return_boot |
return an object of class |
na.action |
a function (default |
scale |
default TRUE, may be FALSE to revert changes made to accommodate |
A list with class htest
and mc.sim
containing the following components:
statistic |
the value of the observed Geary's C. |
parameter |
the rank of the observed Geary's C. |
p.value |
the pseudo p-value of the test. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data, and the number of simulations. |
res |
nsim simulated values of statistic, final value is observed statistic |
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 63-5.
data(oldcol) set.seed(1) sim1 <- geary.mc(COL.OLD$CRIME, nb2listw(COL.nb, style="W"), nsim=99, alternative="less") sim1 mean(sim1$res) var(sim1$res) summary(sim1$res) colold.lags <- nblag(COL.nb, 3) sim2 <- geary.mc(COL.OLD$CRIME, nb2listw(colold.lags[[2]], style="W"), nsim=99) sim2 summary(sim2$res) sim3 <- geary.mc(COL.OLD$CRIME, nb2listw(colold.lags[[3]], style="W"), nsim=99) sim3 summary(sim3$res) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) try(geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.fail)) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.omit) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.omit) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.exclude) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.exclude) try(geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.pass))
data(oldcol) set.seed(1) sim1 <- geary.mc(COL.OLD$CRIME, nb2listw(COL.nb, style="W"), nsim=99, alternative="less") sim1 mean(sim1$res) var(sim1$res) summary(sim1$res) colold.lags <- nblag(COL.nb, 3) sim2 <- geary.mc(COL.OLD$CRIME, nb2listw(colold.lags[[2]], style="W"), nsim=99) sim2 summary(sim2$res) sim3 <- geary.mc(COL.OLD$CRIME, nb2listw(colold.lags[[3]], style="W"), nsim=99) sim3 summary(sim3$res) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) try(geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.fail)) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.omit) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.omit) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.exclude) geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.exclude) try(geary.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.pass))
Geary's test for spatial autocorrelation using a spatial weights matrix in weights list form. The assumptions underlying the test are sensitive to the form of the graph of neighbour relationships and other factors, and results may be checked against those of geary.mc
permutations.
geary.test(x, listw, randomisation=TRUE, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, adjust.n=TRUE, na.action=na.fail, scale=TRUE)
geary.test(x, listw, randomisation=TRUE, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, adjust.n=TRUE, na.action=na.fail, scale=TRUE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
randomisation |
variance of I calculated under the assumption of randomisation, if FALSE normality |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two.sided". |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted |
na.action |
a function (default |
scale |
default TRUE, may be FALSE to revert changes made to accommodate |
A list with class htest
containing the following components:
statistic |
the value of the standard deviate of Geary's C, in the order given in Cliff and Ord 1973, p. 21, which is (EC - C) / sqrt(VC), that is with the sign reversed with respect to the more usual (C - EC) / sqrt(VC); this means that the “greater” alternative for the Geary C test corresponds to the “greater” alternative for Moran's I test. |
p.value |
the p-value of the test. |
estimate |
the value of the observed Geary's C, its expectation and variance under the method assumption. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the assumption used for calculating the standard deviate. |
data.name |
a character string giving the name(s) of the data. |
The derivation of the test (Cliff and Ord, 1981, p. 18) assumes that the weights matrix is symmetric. For inherently non-symmetric matrices, such as k-nearest neighbour matrices, listw2U()
can be used to make the matrix symmetric. In non-symmetric weights matrix cases, the variance of the test statistic may be negative (thanks to Franz Munoz I for a well documented bug report). Geary's C is affected by non-symmetric weights under normality much more than Moran's I. From 0.4-35, the sign of the standard deviate of C is changed to match Cliff and Ord (1973, p. 21).
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 21, Cliff, A. D., Ord, J. K. 1973 Spatial Autocorrelation, Pion, pp. 15-16, 21; Bivand RS, Wong DWS 2018 Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748 doi:10.1007/s11749-018-0599-x
data(oldcol) geary.test(COL.OLD$CRIME, nb2listw(COL.nb, style="W")) geary.test(COL.OLD$CRIME, nb2listw(COL.nb, style="W"), randomisation=FALSE) colold.lags <- nblag(COL.nb, 3) geary.test(COL.OLD$CRIME, nb2listw(colold.lags[[2]], style="W")) geary.test(COL.OLD$CRIME, nb2listw(colold.lags[[3]], style="W"), alternative="greater") print(is.symmetric.nb(COL.nb)) coords.OLD <- cbind(COL.OLD$X, COL.OLD$Y) COL.k4.nb <- knn2nb(knearneigh(coords.OLD, 4)) print(is.symmetric.nb(COL.k4.nb)) geary.test(COL.OLD$CRIME, nb2listw(COL.k4.nb, style="W")) geary.test(COL.OLD$CRIME, nb2listw(COL.k4.nb, style="W"), randomisation=FALSE) cat("Note non-symmetric weights matrix - use listw2U()\n") geary.test(COL.OLD$CRIME, listw2U(nb2listw(COL.k4.nb, style="W"))) geary.test(COL.OLD$CRIME, listw2U(nb2listw(COL.k4.nb, style="W")), randomisation=FALSE) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) try(geary.test(crime, nb2listw(COL.nb, style="W"), na.action=na.fail)) geary.test(crime, nb2listw(COL.nb, style="W"), zero.policy=TRUE, na.action=na.omit) geary.test(crime, nb2listw(COL.nb, style="W"), zero.policy=TRUE, na.action=na.exclude) try(geary.test(crime, nb2listw(COL.nb, style="W"), na.action=na.pass))
data(oldcol) geary.test(COL.OLD$CRIME, nb2listw(COL.nb, style="W")) geary.test(COL.OLD$CRIME, nb2listw(COL.nb, style="W"), randomisation=FALSE) colold.lags <- nblag(COL.nb, 3) geary.test(COL.OLD$CRIME, nb2listw(colold.lags[[2]], style="W")) geary.test(COL.OLD$CRIME, nb2listw(colold.lags[[3]], style="W"), alternative="greater") print(is.symmetric.nb(COL.nb)) coords.OLD <- cbind(COL.OLD$X, COL.OLD$Y) COL.k4.nb <- knn2nb(knearneigh(coords.OLD, 4)) print(is.symmetric.nb(COL.k4.nb)) geary.test(COL.OLD$CRIME, nb2listw(COL.k4.nb, style="W")) geary.test(COL.OLD$CRIME, nb2listw(COL.k4.nb, style="W"), randomisation=FALSE) cat("Note non-symmetric weights matrix - use listw2U()\n") geary.test(COL.OLD$CRIME, listw2U(nb2listw(COL.k4.nb, style="W"))) geary.test(COL.OLD$CRIME, listw2U(nb2listw(COL.k4.nb, style="W")), randomisation=FALSE) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) try(geary.test(crime, nb2listw(COL.nb, style="W"), na.action=na.fail)) geary.test(crime, nb2listw(COL.nb, style="W"), zero.policy=TRUE, na.action=na.omit) geary.test(crime, nb2listw(COL.nb, style="W"), zero.policy=TRUE, na.action=na.exclude) try(geary.test(crime, nb2listw(COL.nb, style="W"), na.action=na.pass))
The global G statistic for spatial autocorrelation, complementing the local Gi LISA measures: localG
.
globalG.test(x, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, adjust.n=TRUE, B1correct=TRUE, adjust.x=TRUE, Arc_all_x=FALSE, na.action=na.fail)
globalG.test(x, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, adjust.n=TRUE, B1correct=TRUE, adjust.x=TRUE, Arc_all_x=FALSE, na.action=na.fail)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two.sided". |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted |
B1correct |
default TRUE, if TRUE, the erratum referenced below: "On page 195, the coefficient of W2 in B1, (just below center of the page) should be 6, not 3." is applied; if FALSE, 3 is used (as in CrimeStat IV) |
adjust.x |
default TRUE, if TRUE, x values of observations with no neighbours are omitted in the denominator of G |
Arc_all_x |
default FALSE, if Arc_all_x=TRUE and adjust.x=TRUE, use the full x vector in part of the denominator term for G |
na.action |
a function (default |
A list with class htest
containing the following components:
statistic |
the value of the standard deviate of Getis-Ord G. |
p.value |
the p-value of the test. |
estimate |
the value of the observed statistic, its expectation and variance. |
alternative |
a character string describing the alternative hypothesis. |
data.name |
a character string giving the name(s) of the data. |
Hisaji ONO [email protected] and Roger Bivand [email protected]
Getis. A, Ord, J. K. 1992 The analysis of spatial association by use of distance statistics, Geographical Analysis, 24, p. 195; see also Getis. A, Ord, J. K. 1993 Erratum, Geographical Analysis, 25, p. 276; Bivand RS, Wong DWS 2018 Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748 doi:10.1007/s11749-018-0599-x
nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE) sidsrate79 <- (1000*nc.sids$SID79)/nc.sids$BIR79 dists <- c(10, 20, 30, 33, 40, 50, 60, 70, 80, 90, 100) ndists <- length(dists) ZG <- vector(mode="list", length=ndists) names(ZG) <- as.character(dists) milesxy <- cbind(nc.sids$east, nc.sids$north) for (i in 1:ndists) { thisnb <- dnearneigh(milesxy, 0, dists[i]) thislw <- nb2listw(thisnb, style="B", zero.policy=TRUE) ZG[[i]] <- globalG.test(sidsrate79, thislw, zero.policy=TRUE) } t(sapply(ZG, function(x) c(x$estimate[1], x$statistic, p.value=unname(x$p.value)))) for (i in 1:ndists) { thisnb <- dnearneigh(milesxy, 0, dists[i]) thislw <- nb2listw(thisnb, style="B", zero.policy=TRUE) ZG[[i]] <- globalG.test(sidsrate79, thislw, zero.policy=TRUE, alternative="two.sided") } t(sapply(ZG, function(x) c(x$estimate[1], x$statistic, p.value=unname(x$p.value)))) data(oldcol) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) res <- try(globalG.test(crime, nb2listw(COL.nb, style="B"), na.action=na.fail)) globalG.test(crime, nb2listw(COL.nb, style="B"), zero.policy=TRUE, na.action=na.omit) globalG.test(crime, nb2listw(COL.nb, style="B"), zero.policy=TRUE, na.action=na.exclude) try(globalG.test(crime, nb2listw(COL.nb, style="B"), na.action=na.pass))
nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE) sidsrate79 <- (1000*nc.sids$SID79)/nc.sids$BIR79 dists <- c(10, 20, 30, 33, 40, 50, 60, 70, 80, 90, 100) ndists <- length(dists) ZG <- vector(mode="list", length=ndists) names(ZG) <- as.character(dists) milesxy <- cbind(nc.sids$east, nc.sids$north) for (i in 1:ndists) { thisnb <- dnearneigh(milesxy, 0, dists[i]) thislw <- nb2listw(thisnb, style="B", zero.policy=TRUE) ZG[[i]] <- globalG.test(sidsrate79, thislw, zero.policy=TRUE) } t(sapply(ZG, function(x) c(x$estimate[1], x$statistic, p.value=unname(x$p.value)))) for (i in 1:ndists) { thisnb <- dnearneigh(milesxy, 0, dists[i]) thislw <- nb2listw(thisnb, style="B", zero.policy=TRUE) ZG[[i]] <- globalG.test(sidsrate79, thislw, zero.policy=TRUE, alternative="two.sided") } t(sapply(ZG, function(x) c(x$estimate[1], x$statistic, p.value=unname(x$p.value)))) data(oldcol) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) res <- try(globalG.test(crime, nb2listw(COL.nb, style="B"), na.action=na.fail)) globalG.test(crime, nb2listw(COL.nb, style="B"), zero.policy=TRUE, na.action=na.omit) globalG.test(crime, nb2listw(COL.nb, style="B"), zero.policy=TRUE, na.action=na.exclude) try(globalG.test(crime, nb2listw(COL.nb, style="B"), na.action=na.pass))
n.comp.nb()
finds the number of disjoint connected subgraphs in the graph depicted by nb.obj
- a spatial neighbours list object.
n.comp.nb(nb.obj)
n.comp.nb(nb.obj)
nb.obj |
a neighbours list object of class |
If attr(nb.obj, "sym")
is FALSE
and igraph::components
is available, the components of the directed graph will be found by a simple breadth-first search; if igraph::components
is not available, the object will be made symmetric (which may be time-consuming with large numbers of neighbours) and the components found by depth-first search. If attr(nb.obj, "sym")
is TRUE
, the components of the directed graph will be found by depth-first search. The time complexity of algorithms used in native code and through igraph::components
is linear in the sum of the number of nodes and the number of edges in the graph, see https://github.com/r-spatial/spdep/issues/160 for details; very dense neighbour objects will have large numbers of edges.
A list of:
nc |
number of disjoint connected subgraphs |
comp.id |
vector with the indices of the disjoint connected subgraphs that
the nodes in |
Nicholas Lewin-Koh [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) coords <- st_coordinates(st_centroid(st_geometry(columbus))) plot(col.gal.nb, coords, col="grey") col2 <- droplinks(col.gal.nb, 21) res <- n.comp.nb(col2) table(res$comp.id) plot(col2, coords, add=TRUE) points(coords, col=res$comp.id, pch=16) run <- FALSE if (require("igraph", quietly=TRUE) && require("spatialreg", quietly=TRUE)) run <- TRUE if (run) { B <- as(nb2listw(col2, style="B", zero.policy=TRUE), "CsparseMatrix") g1 <- graph_from_adjacency_matrix(B, mode="undirected") c1 <- components(g1) print(c1$no == res$nc) } if (run) { print(all.equal(c1$membership, res$comp.id)) } if (run) { print(all.equal(c1$csize, c(table(res$comp.id)), check.attributes=FALSE)) } if (run) { W <- as(nb2listw(col2, style="W", zero.policy=TRUE), "CsparseMatrix") g1W <- graph_from_adjacency_matrix(W, mode="directed", weighted="W") c1W <- components(g1W, mode="weak") print(all.equal(c1W$membership, res$comp.id, check.attributes=FALSE)) } if (run) { data(house, package="spData") house <- sf::st_as_sf(house) k6 <- knn2nb(knearneigh(house, k=6)) is.symmetric.nb(k6) } if (run) { print(k6) } if (run) { length(k6) + sum(card(k6)) } if (run) { # no pre-computed graph components str(attr(k6, "ncomp")) } if (run) { # raising the subgraph compute ceiling to above |N|+|E| computes and stores the # object in the neighbour object set.SubgraphCeiling(180000L) k6 <- knn2nb(knearneigh(house, k=6)) str(attr(k6, "ncomp")) } if (run) { print(k6) } if (run) { system.time(udir <- n.comp.nb(make.sym.nb(k6))) } if (run) { system.time(dir <- n.comp.nb(k6)) } if (run) { udir$nc } if (run) { dir$nc } if (run) { all.equal(dir, udir) }
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) coords <- st_coordinates(st_centroid(st_geometry(columbus))) plot(col.gal.nb, coords, col="grey") col2 <- droplinks(col.gal.nb, 21) res <- n.comp.nb(col2) table(res$comp.id) plot(col2, coords, add=TRUE) points(coords, col=res$comp.id, pch=16) run <- FALSE if (require("igraph", quietly=TRUE) && require("spatialreg", quietly=TRUE)) run <- TRUE if (run) { B <- as(nb2listw(col2, style="B", zero.policy=TRUE), "CsparseMatrix") g1 <- graph_from_adjacency_matrix(B, mode="undirected") c1 <- components(g1) print(c1$no == res$nc) } if (run) { print(all.equal(c1$membership, res$comp.id)) } if (run) { print(all.equal(c1$csize, c(table(res$comp.id)), check.attributes=FALSE)) } if (run) { W <- as(nb2listw(col2, style="W", zero.policy=TRUE), "CsparseMatrix") g1W <- graph_from_adjacency_matrix(W, mode="directed", weighted="W") c1W <- components(g1W, mode="weak") print(all.equal(c1W$membership, res$comp.id, check.attributes=FALSE)) } if (run) { data(house, package="spData") house <- sf::st_as_sf(house) k6 <- knn2nb(knearneigh(house, k=6)) is.symmetric.nb(k6) } if (run) { print(k6) } if (run) { length(k6) + sum(card(k6)) } if (run) { # no pre-computed graph components str(attr(k6, "ncomp")) } if (run) { # raising the subgraph compute ceiling to above |N|+|E| computes and stores the # object in the neighbour object set.SubgraphCeiling(180000L) k6 <- knn2nb(knearneigh(house, k=6)) str(attr(k6, "ncomp")) } if (run) { print(k6) } if (run) { system.time(udir <- n.comp.nb(make.sym.nb(k6))) } if (run) { system.time(dir <- n.comp.nb(k6)) } if (run) { udir$nc } if (run) { dir$nc } if (run) { all.equal(dir, udir) }
Functions return a graph object containing a list with the vertex
coordinates and the to and from indices defining the edges. Some/all of these functions assume that the coordinates are not exactly regularly spaced. The helper
function graph2nb
converts a graph
object into a neighbour list. The plot functions plot the graph objects.
gabrielneigh(coords, nnmult=3) relativeneigh(coords, nnmult=3) soi.graph(tri.nb, coords, quadsegs=10) graph2nb(gob, row.names=NULL,sym=FALSE) ## S3 method for class 'Gabriel' plot(x, show.points=FALSE, add=FALSE, linecol=par(col), ...) ## S3 method for class 'relative' plot(x, show.points=FALSE, add=FALSE, linecol=par(col),...)
gabrielneigh(coords, nnmult=3) relativeneigh(coords, nnmult=3) soi.graph(tri.nb, coords, quadsegs=10) graph2nb(gob, row.names=NULL,sym=FALSE) ## S3 method for class 'Gabriel' plot(x, show.points=FALSE, add=FALSE, linecol=par(col), ...) ## S3 method for class 'relative' plot(x, show.points=FALSE, add=FALSE, linecol=par(col),...)
coords |
matrix of region point coordinates or SpatialPoints object or |
nnmult |
scaling factor for memory allocation, default 3; if higher values are required, the function will exit with an error; example below thanks to Dan Putler |
tri.nb |
a neighbor list created from tri2nb |
quadsegs |
number of line segments making a quarter circle buffer, see the |
gob |
a graph object created from any of the graph funtions |
row.names |
character vector of region ids to be added to the
neighbours list as attribute |
sym |
a logical argument indicating whether or not neighbors should be symetric (if i->j then j->i) |
x |
object to be plotted |
show.points |
(logical) add points to plot |
add |
(logical) add to existing plot |
linecol |
edge plotting colour |
... |
further graphical parameters as in |
The graph functions produce graphs on a 2d point set that are all subgraphs of the Delaunay triangulation. The relative neighbor graph is defined by the relation, x and y are neighbors if
where d() is the distance, S is the set of points and z is an arbitrary point in S. The Gabriel graph is a subgraph of the delaunay triangulation and has the relative neighbor graph as a sub-graph. The relative neighbor graph is defined by the relation x and y are Gabriel neighbors if
where x,y,z and S are as before. The sphere of influence graph is
defined for a finite point set S, let be the distance from point x
to its nearest neighbor in S, and
is the circle centered on x. Then
x and y are SOI neigbors iff
and
intersect in at
least 2 places. From 2016-05-31, Computational Geometry in C code replaced by calls to functions in dbscan and sf; with a large
quadsegs=
argument, the behaviour of the function is the same, otherwise buffer intersections only closely approximate the original function.
See card
for details of “nb” objects.
A list of class Graph
with the following elements
np |
number of input points |
from |
array of origin ids |
to |
array of destination ids |
nedges |
number of edges in graph |
x |
input x coordinates |
y |
input y coordinates |
The helper functions return an nb
object with a list of integer
vectors containing neighbour region number ids.
Nicholas Lewin-Koh [email protected]
Matula, D. W. and Sokal R. R. 1980, Properties of Gabriel graphs relevant to geographic variation research and the clustering of points in the plane, Geographic Analysis, 12(3), pp. 205-222.
Toussaint, G. T. 1980, The relative neighborhood graph of a finite planar set, Pattern Recognition, 12(4), pp. 261-268.
Kirkpatrick, D. G. and Radke, J. D. 1985, A framework for computational morphology. In Computational Geometry, Ed. G. T. Toussaint, North Holland.
knearneigh
, dnearneigh
,
knn2nb
, card
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) sf_obj <- st_centroid(st_geometry(columbus), of_largest_polygon) sp_obj <- as(sf_obj, "Spatial") coords <- st_coordinates(sf_obj) suppressMessages(col.tri.nb <- tri2nb(coords)) col.gab.nb <- graph2nb(gabrielneigh(coords), sym=TRUE) col.rel.nb <- graph2nb(relativeneigh(coords), sym=TRUE) par(mfrow=c(2,2)) plot(st_geometry(columbus), border="grey") plot(col.tri.nb,coords,add=TRUE) title(main="Delaunay Triangulation", cex.main=0.6) plot(st_geometry(columbus), border="grey") plot(col.gab.nb, coords, add=TRUE) title(main="Gabriel Graph", cex.main=0.6) plot(st_geometry(columbus), border="grey") plot(col.rel.nb, coords, add=TRUE) title(main="Relative Neighbor Graph", cex.main=0.6) plot(st_geometry(columbus), border="grey") if (require("dbscan", quietly=TRUE)) { col.soi.nb <- graph2nb(soi.graph(col.tri.nb,coords), sym=TRUE) plot(col.soi.nb, coords, add=TRUE) title(main="Sphere of Influence Graph", cex.main=0.6) } par(mfrow=c(1,1)) col.tri.nb_sf <- tri2nb(sf_obj) all.equal(col.tri.nb, col.tri.nb_sf, check.attributes=FALSE) col.tri.nb_sp <- tri2nb(sp_obj) all.equal(col.tri.nb, col.tri.nb_sp, check.attributes=FALSE) if (require("dbscan", quietly=TRUE)) { col.soi.nb_sf <- graph2nb(soi.graph(col.tri.nb, sf_obj), sym=TRUE) all.equal(col.soi.nb, col.soi.nb_sf, check.attributes=FALSE) col.soi.nb_sp <- graph2nb(soi.graph(col.tri.nb, sp_obj), sym=TRUE) all.equal(col.soi.nb, col.soi.nb_sp, check.attributes=FALSE) } col.gab.nb_sp <- graph2nb(gabrielneigh(sp_obj), sym=TRUE) all.equal(col.gab.nb, col.gab.nb_sp, check.attributes=FALSE) col.gab.nb_sf <- graph2nb(gabrielneigh(sf_obj), sym=TRUE) all.equal(col.gab.nb, col.gab.nb_sf, check.attributes=FALSE) col.rel.nb_sp <- graph2nb(relativeneigh(sp_obj), sym=TRUE) all.equal(col.rel.nb, col.rel.nb_sp, check.attributes=FALSE) col.rel.nb_sf <- graph2nb(relativeneigh(sf_obj), sym=TRUE) all.equal(col.rel.nb, col.rel.nb_sf, check.attributes=FALSE) dx <- rep(0.25*0:4,5) dy <- c(rep(0,5),rep(0.25,5),rep(0.5,5), rep(0.75,5),rep(1,5)) m <- cbind(c(dx, dx, 3+dx, 3+dx), c(dy, 3+dy, dy, 3+dy)) cat(try(res <- gabrielneigh(m), silent=TRUE), "\n") res <- gabrielneigh(m, nnmult=4) summary(graph2nb(res)) grd <- as.matrix(expand.grid(x=1:5, y=1:5)) #gridded data r2 <- gabrielneigh(grd) set.seed(1) grd1 <- as.matrix(expand.grid(x=1:5, y=1:5)) + matrix(runif(50, .0001, .0006), nrow=25) r3 <- gabrielneigh(grd1) opar <- par(mfrow=c(1,2)) plot(r2, show=TRUE, linecol=2) plot(r3, show=TRUE, linecol=2) par(opar) # example of reading points with readr::read_csv() yielding a tibble load(system.file("etc/misc/coords.rda", package="spdep")) class(coords) graph2nb(gabrielneigh(coords)) graph2nb(relativeneigh(coords))
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) sf_obj <- st_centroid(st_geometry(columbus), of_largest_polygon) sp_obj <- as(sf_obj, "Spatial") coords <- st_coordinates(sf_obj) suppressMessages(col.tri.nb <- tri2nb(coords)) col.gab.nb <- graph2nb(gabrielneigh(coords), sym=TRUE) col.rel.nb <- graph2nb(relativeneigh(coords), sym=TRUE) par(mfrow=c(2,2)) plot(st_geometry(columbus), border="grey") plot(col.tri.nb,coords,add=TRUE) title(main="Delaunay Triangulation", cex.main=0.6) plot(st_geometry(columbus), border="grey") plot(col.gab.nb, coords, add=TRUE) title(main="Gabriel Graph", cex.main=0.6) plot(st_geometry(columbus), border="grey") plot(col.rel.nb, coords, add=TRUE) title(main="Relative Neighbor Graph", cex.main=0.6) plot(st_geometry(columbus), border="grey") if (require("dbscan", quietly=TRUE)) { col.soi.nb <- graph2nb(soi.graph(col.tri.nb,coords), sym=TRUE) plot(col.soi.nb, coords, add=TRUE) title(main="Sphere of Influence Graph", cex.main=0.6) } par(mfrow=c(1,1)) col.tri.nb_sf <- tri2nb(sf_obj) all.equal(col.tri.nb, col.tri.nb_sf, check.attributes=FALSE) col.tri.nb_sp <- tri2nb(sp_obj) all.equal(col.tri.nb, col.tri.nb_sp, check.attributes=FALSE) if (require("dbscan", quietly=TRUE)) { col.soi.nb_sf <- graph2nb(soi.graph(col.tri.nb, sf_obj), sym=TRUE) all.equal(col.soi.nb, col.soi.nb_sf, check.attributes=FALSE) col.soi.nb_sp <- graph2nb(soi.graph(col.tri.nb, sp_obj), sym=TRUE) all.equal(col.soi.nb, col.soi.nb_sp, check.attributes=FALSE) } col.gab.nb_sp <- graph2nb(gabrielneigh(sp_obj), sym=TRUE) all.equal(col.gab.nb, col.gab.nb_sp, check.attributes=FALSE) col.gab.nb_sf <- graph2nb(gabrielneigh(sf_obj), sym=TRUE) all.equal(col.gab.nb, col.gab.nb_sf, check.attributes=FALSE) col.rel.nb_sp <- graph2nb(relativeneigh(sp_obj), sym=TRUE) all.equal(col.rel.nb, col.rel.nb_sp, check.attributes=FALSE) col.rel.nb_sf <- graph2nb(relativeneigh(sf_obj), sym=TRUE) all.equal(col.rel.nb, col.rel.nb_sf, check.attributes=FALSE) dx <- rep(0.25*0:4,5) dy <- c(rep(0,5),rep(0.25,5),rep(0.5,5), rep(0.75,5),rep(1,5)) m <- cbind(c(dx, dx, 3+dx, 3+dx), c(dy, 3+dy, dy, 3+dy)) cat(try(res <- gabrielneigh(m), silent=TRUE), "\n") res <- gabrielneigh(m, nnmult=4) summary(graph2nb(res)) grd <- as.matrix(expand.grid(x=1:5, y=1:5)) #gridded data r2 <- gabrielneigh(grd) set.seed(1) grd1 <- as.matrix(expand.grid(x=1:5, y=1:5)) + matrix(runif(50, .0001, .0006), nrow=25) r3 <- gabrielneigh(grd1) opar <- par(mfrow=c(1,2)) plot(r2, show=TRUE, linecol=2) plot(r3, show=TRUE, linecol=2) par(opar) # example of reading points with readr::read_csv() yielding a tibble load(system.file("etc/misc/coords.rda", package="spdep")) class(coords) graph2nb(gabrielneigh(coords)) graph2nb(relativeneigh(coords))
The function builds a neighbours list for a grid topology. It works for a k-dimentional grid topology, k>=1.
grid2nb(grid, d = [email protected], queen = TRUE, nb = TRUE, self = FALSE)
grid2nb(grid, d = grid@cells.dim, queen = TRUE, nb = TRUE, self = FALSE)
grid |
An object of class |
d |
A scalar (for one dimentional grid) or a length k vector specyfying the number of grid cells in each direction of the k dimentions. |
queen |
Logical. Default is TRUE. To inform if the queen neighbourhood structure should be considered. If FALSE, only a hyper-cube with a common face will be considered neighbour. If TRUE, a single shared coordinate meets the contiguity condition. |
nb |
Default TRUE. If TRUE, return the result as a neighbours
list with class |
self |
Default FALSE, to indicate if the hyper-cube neighbour itself should be considered a neighbour. |
Either a matrix, if “nb” is FALSE or a neighbours list with
class nb
. See card
for details of “nb”
objects.
This applies to a k-dimentional grid topology.
Elias T Krainski [email protected]
nb <- grid2nb(d = c(5L, 5L, 5L)) nb summary(nb) if (require("sp", quietly=TRUE)) { gt <- GridTopology(c(.125,.1), c(.25,.2), c(4L, 5L)) nb1 <- grid2nb(gt, queen = FALSE) nb2 <- grid2nb(gt) sg <- SpatialGrid(gt) plot(sg, lwd=3) plot(nb1, coordinates(sg), add=TRUE, lty=2, col=2, lwd=2) plot(nb2, coordinates(sg), add=TRUE, lty=3, col=4) str(grid2nb(d=5)) }
nb <- grid2nb(d = c(5L, 5L, 5L)) nb summary(nb) if (require("sp", quietly=TRUE)) { gt <- GridTopology(c(.125,.1), c(.25,.2), c(4L, 5L)) nb1 <- grid2nb(gt, queen = FALSE) nb2 <- grid2nb(gt) sg <- SpatialGrid(gt) plot(sg, lwd=3) plot(nb1, coordinates(sg), add=TRUE, lty=2, col=2, lwd=2) plot(nb2, coordinates(sg), add=TRUE, lty=3, col=4) str(grid2nb(d=5)) }
Used to return a factor showing so-called cluster classification for local indicators of spatial association for local Moran's I, local Geary's C (and its multivariate variant) and local Getis-Ord G. This factor vector can be added to a spatial object for mapping. When obj
is of class licd
, a list of up to six factors for measures of local composition (analytical and permutation), local configuration (analytical and permutation), and combined measures, both the interaction of composition and configuration, and a simplified recoding of these.
hotspot(obj, ...) ## Default S3 method: hotspot(obj, ...) ## S3 method for class 'localmoran' hotspot(obj, Prname, cutoff=0.005, quadrant.type="mean", p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'summary.localmoransad' hotspot(obj, Prname, cutoff=0.005, quadrant.type="mean", p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'data.frame.localmoranex' hotspot(obj, Prname, cutoff=0.005, quadrant.type="mean", p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'localG' hotspot(obj, Prname, cutoff=0.005, p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'localC' hotspot(obj, Prname, cutoff=0.005, p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'licd' hotspot(obj, type = "both", cutoff = 0.05, p.adjust = "none", droplevels = TRUE, control = list(), ...)
hotspot(obj, ...) ## Default S3 method: hotspot(obj, ...) ## S3 method for class 'localmoran' hotspot(obj, Prname, cutoff=0.005, quadrant.type="mean", p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'summary.localmoransad' hotspot(obj, Prname, cutoff=0.005, quadrant.type="mean", p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'data.frame.localmoranex' hotspot(obj, Prname, cutoff=0.005, quadrant.type="mean", p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'localG' hotspot(obj, Prname, cutoff=0.005, p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'localC' hotspot(obj, Prname, cutoff=0.005, p.adjust="fdr", droplevels=TRUE, ...) ## S3 method for class 'licd' hotspot(obj, type = "both", cutoff = 0.05, p.adjust = "none", droplevels = TRUE, control = list(), ...)
obj |
An object of class |
Prname |
A character string, the name of the column containing the probability values to be classified by cluster type if found “interesting” |
cutoff |
Default 0.005, the probability value cutoff larger than which the observation is not found “interesting” |
p.adjust |
Default |
droplevels |
Default |
quadrant.type |
Default |
type |
When |
control |
When |
... |
other arguments passed to methods. |
A factor showing so-called cluster classification for local indicators of spatial association. When obj
is of class licd
, a list of up to six factors for measures of local composition (analytical and permutation), local configuration (analytical and permutation), and combined measures, both the interaction of composition and configuration, and a simplified recoding of these.
Roger Bivand
orig <- spData::africa.rook.nb listw <- nb2listw(orig) x <- spData::afcon$totcon set.seed(1) C <- localC_perm(x, listw) Ch <- hotspot(C, Prname="Pr(z != E(Ci)) Sim", cutoff=0.05, p.adjust="none") table(addNA(Ch)) set.seed(1) I <- localmoran_perm(x, listw) Ih <- hotspot(I, Prname="Pr(z != E(Ii)) Sim", cutoff=0.05, p.adjust="none") table(addNA(Ih)) Is <- summary(localmoran.sad(lm(x ~ 1), nb=orig)) Ish <- hotspot(Is, Prname="Pr. (Sad)", cutoff=0.05, p.adjust="none") table(addNA(Ish)) Ie <- as.data.frame(localmoran.exact(lm(x ~ 1), nb=orig)) Ieh <- hotspot(Ie, Prname="Pr. (exact)", cutoff=0.05, p.adjust="none") table(addNA(Ieh)) set.seed(1) G <- localG_perm(x, listw) Gh <- hotspot(G, Prname="Pr(z != E(Gi)) Sim", cutoff=0.05, p.adjust="none") table(addNA(Gh))
orig <- spData::africa.rook.nb listw <- nb2listw(orig) x <- spData::afcon$totcon set.seed(1) C <- localC_perm(x, listw) Ch <- hotspot(C, Prname="Pr(z != E(Ci)) Sim", cutoff=0.05, p.adjust="none") table(addNA(Ch)) set.seed(1) I <- localmoran_perm(x, listw) Ih <- hotspot(I, Prname="Pr(z != E(Ii)) Sim", cutoff=0.05, p.adjust="none") table(addNA(Ih)) Is <- summary(localmoran.sad(lm(x ~ 1), nb=orig)) Ish <- hotspot(Is, Prname="Pr. (Sad)", cutoff=0.05, p.adjust="none") table(addNA(Ish)) Ie <- as.data.frame(localmoran.exact(lm(x ~ 1), nb=orig)) Ieh <- hotspot(Ie, Prname="Pr. (exact)", cutoff=0.05, p.adjust="none") table(addNA(Ieh)) set.seed(1) G <- localG_perm(x, listw) Gh <- hotspot(G, Prname="Pr(z != E(Gi)) Sim", cutoff=0.05, p.adjust="none") table(addNA(Gh))
The function includes the region itself in its own list of neighbours, and sets attribute "self.included" to TRUE; remove.self
reverts the effects of include.self
.
include.self(nb) remove.self(nb)
include.self(nb) remove.self(nb)
nb |
input neighbours list of class |
The function returns an object of class nb
with a list of integer vectors containing neighbour region number ids; attribute "self.included" is set to TRUE.
Roger Bivand [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) coords <- st_coordinates(st_centroid(columbus)) summary(col.gal.nb, coords) summary(include.self(col.gal.nb), coords) summary(remove.self(include.self(col.gal.nb)), coords)
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) coords <- st_coordinates(st_centroid(columbus)) summary(col.gal.nb, coords) summary(include.self(col.gal.nb), coords) summary(remove.self(include.self(col.gal.nb)), coords)
Checks a neighbours list for symmetry/transitivity (if i is a neighbour of j,
then j is a neighbour of i). This holds for distance and contiguity based
neighbours, but not for k-nearest neighbours. The helper function
sym.attr.nb()
calls is.symmetric.nb()
to set the sym
attribute if needed, and make.sym.nb
makes a non-symmetric list symmetric by adding neighbors. is.symmetric.glist
checks a list of general weights corresponding to neighbours for symmetry for symmetric neighbours.
is.symmetric.nb(nb, verbose = NULL, force = FALSE) sym.attr.nb(nb) make.sym.nb(nb) old.make.sym.nb(nb) is.symmetric.glist(nb, glist)
is.symmetric.nb(nb, verbose = NULL, force = FALSE) sym.attr.nb(nb) make.sym.nb(nb) old.make.sym.nb(nb) is.symmetric.glist(nb, glist)
nb |
an object of class |
verbose |
default NULL, use global option value; if TRUE prints non-matching pairs |
force |
do not respect a neighbours list |
glist |
list of general weights corresponding to neighbours |
TRUE if symmetric, FALSE if not; is.symmetric.glist returns a value with an attribute, "d", indicating for failed symmetry the largest failing value.
A new version of make.sym.nb
by Bjarke Christensen is now included. The older version has been renamed old.make.sym.nb
, and their comparison constitutes a nice demonstration of vectorising speedup using sapply
and lapply
rather than loops. When any no-neighbour observations are present, old.make.sym.nb
is used.
Roger Bivand [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) coords <- st_coordinates(st_centroid(columbus)) ind <- row.names(as(columbus, "Spatial")) print(is.symmetric.nb(col.gal.nb, verbose=TRUE, force=TRUE)) k4 <- knn2nb(knearneigh(coords, k=4), row.names=ind) k4 <- sym.attr.nb(k4) print(is.symmetric.nb(k4)) k4.sym <- make.sym.nb(k4) print(is.symmetric.nb(k4.sym))
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) coords <- st_coordinates(st_centroid(columbus)) ind <- row.names(as(columbus, "Spatial")) print(is.symmetric.nb(col.gal.nb, verbose=TRUE, force=TRUE)) k4 <- knn2nb(knearneigh(coords, k=4), row.names=ind) k4 <- sym.attr.nb(k4) print(is.symmetric.nb(k4)) k4.sym <- make.sym.nb(k4) print(is.symmetric.nb(k4.sym))
A permutation test for same colour join count statistics calculated by using nsim random permutations of fx for the given spatial weighting scheme, to establish the ranks of the observed statistics (for each colour) in relation to the nsim simulated values.
joincount.mc(fx, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL)
joincount.mc(fx, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL)
fx |
a factor of the same length as the neighbours and weights objects in listw |
listw |
a |
nsim |
number of permutations |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "two.sided", or "less". |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
A list with class jclist
of lists with class htest
and mc.sim
for each of the k colours containing the following components:
statistic |
the value of the observed statistic. |
parameter |
the rank of the observed statistic. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data. |
p.value |
the pseudo p-value of the test. |
alternative |
a character string describing the alternative hypothesis. |
estimate |
the mean and variance of the simulated distribution. |
res |
nsim simulated values of statistic, the final element is the observed statistic |
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 63-5.
data(oldcol) HICRIME <- cut(COL.OLD$CRIME, breaks=c(0,35,80), labels=c("low","high")) names(HICRIME) <- rownames(COL.OLD) joincount.mc(HICRIME, nb2listw(COL.nb, style="B"), nsim=99, alternative="two.sided") joincount.test(HICRIME, nb2listw(COL.nb, style="B"), alternative="two.sided")
data(oldcol) HICRIME <- cut(COL.OLD$CRIME, breaks=c(0,35,80), labels=c("low","high")) names(HICRIME) <- rownames(COL.OLD) joincount.mc(HICRIME, nb2listw(COL.nb, style="B"), nsim=99, alternative="two.sided") joincount.test(HICRIME, nb2listw(COL.nb, style="B"), alternative="two.sided")
A function for tallying join counts between same-colour and different colour spatial objects, where neighbour relations are defined by a weights list. Given the global counts in each colour, expected counts and variances are calculated under non-free sampling, and a z-value reported. Since multiple tests are reported, no p-values are given, allowing the user to adjust the significance level applied. Jtot is the count of all different-colour joins.
joincount.multi(fx, listw, zero.policy = attr(listw, "zero.policy"), spChk = NULL, adjust.n=TRUE) ## S3 method for class 'jcmulti' print(x, ...)
joincount.multi(fx, listw, zero.policy = attr(listw, "zero.policy"), spChk = NULL, adjust.n=TRUE) ## S3 method for class 'jcmulti' print(x, ...)
fx |
a factor of the same length as the neighbours and weights objects in listw |
listw |
a |
zero.policy |
default |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted consistently (up to and including spdep 0.3-28 the adjustment was inconsistent - thanks to Tomoki NAKAYA for a careful bug report) |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
x |
object to be printed |
... |
arguments to be passed through for printing |
A matrix with class jcmulti
with row and column names for observed and expected counts, variance, and z-value.
The derivation of the test (Cliff and Ord, 1981, p. 18) assumes that the weights matrix is symmetric. For inherently non-symmetric matrices, such as k-nearest neighbour matrices, listw2U()
can be used to make the matrix symmetric. In non-symmetric weights matrix cases, the variance of the test statistic may be negative.
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 20; Upton, G., Fingleton, B. 1985 Spatial data analysis by example: point pattern and qualitative data, Wiley, pp. 158–170.
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B") joincount.multi(HICRIME, lw) col_geoms <- st_geometry(columbus) col_geoms[21] <- st_buffer(col_geoms[21], dist=-0.05) st_geometry(columbus) <- col_geoms (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B", zero.policy=TRUE) joincount.multi(HICRIME, lw) ## Not run: data(oldcol) HICRIME <- cut(COL.OLD$CRIME, breaks=c(0,35,80), labels=c("low","high")) names(HICRIME) <- rownames(COL.OLD) joincount.multi(HICRIME, nb2listw(COL.nb, style="B")) data(hopkins, package="spData") image(1:32, 1:32, hopkins[5:36,36:5], breaks=c(-0.5, 3.5, 20), col=c("white", "black")) box() hopkins.rook.nb <- cell2nb(32, 32, type="rook") unlist(spweights.constants(nb2listw(hopkins.rook.nb, style="B"))) hopkins.queen.nb <- cell2nb(32, 32, type="queen") hopkins.bishop.nb <- diffnb(hopkins.rook.nb, hopkins.queen.nb, verbose=FALSE) hopkins4 <- hopkins[5:36,36:5] hopkins4[which(hopkins4 > 3, arr.ind=TRUE)] <- 4 hopkins4.f <- factor(hopkins4) table(hopkins4.f) joincount.multi(hopkins4.f, nb2listw(hopkins.rook.nb, style="B")) cat("replicates Upton & Fingleton table 3.4 (p. 166)\n") joincount.multi(hopkins4.f, nb2listw(hopkins.bishop.nb, style="B")) cat("replicates Upton & Fingleton table 3.6 (p. 168)\n") joincount.multi(hopkins4.f, nb2listw(hopkins.queen.nb, style="B")) cat("replicates Upton & Fingleton table 3.7 (p. 169)\n") ## End(Not run) GDAL37 <- as.numeric_version(unname(sf_extSoftVersion()["GDAL"])) >= "3.7.0" file <- "etc/shapes/GB_2024_southcoast_50m.gpkg.zip" zipfile <- system.file(file, package="spdep") if (GDAL37) { sc50m <- st_read(zipfile) } else { td <- tempdir() bn <- sub(".zip", "", basename(file), fixed=TRUE) target <- unzip(zipfile, files=bn, exdir=td) sc50m <- st_read(target) } sc50m$Winner <- factor(sc50m$Winner, levels=c("Con", "Green", "Lab", "LD")) plot(sc50m[,"Winner"], pal=c("#2297E6", "#61D04F", "#DF536B", "#F5C710")) nb_sc_50m <- poly2nb(sc50m, row.names=as.character(sc50m$Constituency)) sub2 <- attr(nb_sc_50m, "region.id")[attr(nb_sc_50m, "ncomp")$comp.id == 2L] iowe <- match(sub2[1], attr(nb_sc_50m, "region.id")) diowe <- c(st_distance(sc50m[iowe,], sc50m)) meet_criterion <- sum(diowe <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(diowe)[1:meet_criterion]] nb_sc_50m_iowe <- addlinks1(nb_sc_50m, from = cands[1], to = cands[3:meet_criterion]) ioww <- match(sub2[2], attr(nb_sc_50m, "region.id")) dioww <- c(st_distance(sc50m[ioww,], sc50m)) meet_criterion <- sum(dioww <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(dioww)[1:meet_criterion]] nb_sc_50m_iow <- addlinks1(nb_sc_50m_iowe, from = cands[2], to = cands[3:meet_criterion]) nb_sc_1_2 <- nblag_cumul(nblag(nb_sc_50m_iow, 2)) joincount.multi(sc50m$Winner, nb2listw(nb_sc_1_2, style="B"))
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B") joincount.multi(HICRIME, lw) col_geoms <- st_geometry(columbus) col_geoms[21] <- st_buffer(col_geoms[21], dist=-0.05) st_geometry(columbus) <- col_geoms (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B", zero.policy=TRUE) joincount.multi(HICRIME, lw) ## Not run: data(oldcol) HICRIME <- cut(COL.OLD$CRIME, breaks=c(0,35,80), labels=c("low","high")) names(HICRIME) <- rownames(COL.OLD) joincount.multi(HICRIME, nb2listw(COL.nb, style="B")) data(hopkins, package="spData") image(1:32, 1:32, hopkins[5:36,36:5], breaks=c(-0.5, 3.5, 20), col=c("white", "black")) box() hopkins.rook.nb <- cell2nb(32, 32, type="rook") unlist(spweights.constants(nb2listw(hopkins.rook.nb, style="B"))) hopkins.queen.nb <- cell2nb(32, 32, type="queen") hopkins.bishop.nb <- diffnb(hopkins.rook.nb, hopkins.queen.nb, verbose=FALSE) hopkins4 <- hopkins[5:36,36:5] hopkins4[which(hopkins4 > 3, arr.ind=TRUE)] <- 4 hopkins4.f <- factor(hopkins4) table(hopkins4.f) joincount.multi(hopkins4.f, nb2listw(hopkins.rook.nb, style="B")) cat("replicates Upton & Fingleton table 3.4 (p. 166)\n") joincount.multi(hopkins4.f, nb2listw(hopkins.bishop.nb, style="B")) cat("replicates Upton & Fingleton table 3.6 (p. 168)\n") joincount.multi(hopkins4.f, nb2listw(hopkins.queen.nb, style="B")) cat("replicates Upton & Fingleton table 3.7 (p. 169)\n") ## End(Not run) GDAL37 <- as.numeric_version(unname(sf_extSoftVersion()["GDAL"])) >= "3.7.0" file <- "etc/shapes/GB_2024_southcoast_50m.gpkg.zip" zipfile <- system.file(file, package="spdep") if (GDAL37) { sc50m <- st_read(zipfile) } else { td <- tempdir() bn <- sub(".zip", "", basename(file), fixed=TRUE) target <- unzip(zipfile, files=bn, exdir=td) sc50m <- st_read(target) } sc50m$Winner <- factor(sc50m$Winner, levels=c("Con", "Green", "Lab", "LD")) plot(sc50m[,"Winner"], pal=c("#2297E6", "#61D04F", "#DF536B", "#F5C710")) nb_sc_50m <- poly2nb(sc50m, row.names=as.character(sc50m$Constituency)) sub2 <- attr(nb_sc_50m, "region.id")[attr(nb_sc_50m, "ncomp")$comp.id == 2L] iowe <- match(sub2[1], attr(nb_sc_50m, "region.id")) diowe <- c(st_distance(sc50m[iowe,], sc50m)) meet_criterion <- sum(diowe <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(diowe)[1:meet_criterion]] nb_sc_50m_iowe <- addlinks1(nb_sc_50m, from = cands[1], to = cands[3:meet_criterion]) ioww <- match(sub2[2], attr(nb_sc_50m, "region.id")) dioww <- c(st_distance(sc50m[ioww,], sc50m)) meet_criterion <- sum(dioww <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(dioww)[1:meet_criterion]] nb_sc_50m_iow <- addlinks1(nb_sc_50m_iowe, from = cands[2], to = cands[3:meet_criterion]) nb_sc_1_2 <- nblag_cumul(nblag(nb_sc_50m_iow, 2)) joincount.multi(sc50m$Winner, nb2listw(nb_sc_1_2, style="B"))
The BB join count test for spatial autocorrelation using a spatial
weights matrix in weights list form for testing whether same-colour joins
occur more frequently than would be expected if the zones were labelled
in a spatially random way. The assumptions underlying the test are
sensitive to the form of the graph of neighbour relationships and other
factors, and results may be checked against those of joincount.mc
permutations.
joincount.test(fx, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", sampling="nonfree", spChk=NULL, adjust.n=TRUE) ## S3 method for class 'jclist' print(x, ...)
joincount.test(fx, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", sampling="nonfree", spChk=NULL, adjust.n=TRUE) ## S3 method for class 'jclist' print(x, ...)
fx |
a factor of the same length as the neighbours and weights objects in listw |
listw |
a |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two.sided". |
sampling |
default “nonfree”, may be “free” |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted consistently (up to and including spdep 0.3-28 the adjustment was inconsistent - thanks to Tomoki NAKAYA for a careful bug report) |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
x |
object to be printed |
... |
arguments to be passed through for printing |
A list with class jclist
of lists with class htest
for each of the k colours containing the following components:
statistic |
the value of the standard deviate of the join count statistic. |
p.value |
the p-value of the test. |
estimate |
the value of the observed statistic, its expectation and variance under non-free sampling. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data. |
The derivation of the test (Cliff and Ord, 1981, p. 18) assumes that the weights matrix is symmetric. For inherently non-symmetric matrices, such as k-nearest neighbour matrices, listw2U()
can be used to make the matrix symmetric. In non-symmetric weights matrix cases, the variance of the test statistic may be negative.
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, pp. 19-20.
joincount.mc
, joincount.multi
, listw2U
data(oldcol) HICRIME <- cut(COL.OLD$CRIME, breaks=c(0,35,80), labels=c("low","high")) names(HICRIME) <- rownames(COL.OLD) joincount.test(HICRIME, nb2listw(COL.nb, style="B")) joincount.test(HICRIME, nb2listw(COL.nb, style="B"), sampling="free") joincount.test(HICRIME, nb2listw(COL.nb, style="C")) joincount.test(HICRIME, nb2listw(COL.nb, style="S")) joincount.test(HICRIME, nb2listw(COL.nb, style="W")) by(card(COL.nb), HICRIME, summary) print(is.symmetric.nb(COL.nb)) coords.OLD <- cbind(COL.OLD$X, COL.OLD$Y) COL.k4.nb <- knn2nb(knearneigh(coords.OLD, 4)) print(is.symmetric.nb(COL.k4.nb)) joincount.test(HICRIME, nb2listw(COL.k4.nb, style="B")) cat("Note non-symmetric weights matrix - use listw2U()\n") joincount.test(HICRIME, listw2U(nb2listw(COL.k4.nb, style="B"))) columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B") joincount.test(HICRIME, lw) col_geoms <- st_geometry(columbus) col_geoms[21] <- st_buffer(col_geoms[21], dist=-0.05) st_geometry(columbus) <- col_geoms (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B", zero.policy=TRUE) joincount.test(HICRIME, lw)
data(oldcol) HICRIME <- cut(COL.OLD$CRIME, breaks=c(0,35,80), labels=c("low","high")) names(HICRIME) <- rownames(COL.OLD) joincount.test(HICRIME, nb2listw(COL.nb, style="B")) joincount.test(HICRIME, nb2listw(COL.nb, style="B"), sampling="free") joincount.test(HICRIME, nb2listw(COL.nb, style="C")) joincount.test(HICRIME, nb2listw(COL.nb, style="S")) joincount.test(HICRIME, nb2listw(COL.nb, style="W")) by(card(COL.nb), HICRIME, summary) print(is.symmetric.nb(COL.nb)) coords.OLD <- cbind(COL.OLD$X, COL.OLD$Y) COL.k4.nb <- knn2nb(knearneigh(coords.OLD, 4)) print(is.symmetric.nb(COL.k4.nb)) joincount.test(HICRIME, nb2listw(COL.k4.nb, style="B")) cat("Note non-symmetric weights matrix - use listw2U()\n") joincount.test(HICRIME, listw2U(nb2listw(COL.k4.nb, style="B"))) columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B") joincount.test(HICRIME, lw) col_geoms <- st_geometry(columbus) col_geoms[21] <- st_buffer(col_geoms[21], dist=-0.05) st_geometry(columbus) <- col_geoms (nb <- poly2nb(columbus)) lw <- nb2listw(nb, style="B", zero.policy=TRUE) joincount.test(HICRIME, lw)
The function returns a matrix with the indices of points belonging to the set of the k nearest neighbours of each other. If longlat = TRUE, Great Circle distances are used. A warning will be given if identical points are found.
knearneigh(x, k=1, longlat = NULL, use_kd_tree=TRUE)
knearneigh(x, k=1, longlat = NULL, use_kd_tree=TRUE)
x |
matrix of point coordinates, an object inheriting from SpatialPoints or an |
k |
number of nearest neighbours to be returned; where identical points are present, |
longlat |
TRUE if point coordinates are longitude-latitude decimal degrees, in which case distances are measured in kilometers; if x is a SpatialPoints object or an |
use_kd_tree |
logical value, if the dbscan package is available, use for finding k nearest neighbours when coordinates are planar, and when there are no identical points; from https://github.com/r-spatial/spdep/issues/38, the input data may have more than two columns if dbscan is used |
The underlying legacy C code is based on the knn
function in the class package.
A list of class knn
nn |
integer matrix of region number ids |
np |
number of input points |
k |
input required k |
dimension |
number of columns of x |
x |
input coordinates |
Roger Bivand [email protected]
knn
, dnearneigh
,
knn2nb
, kNN
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE) col.knn <- knearneigh(coords, k=4) plot(st_geometry(columbus), border="grey") plot(knn2nb(col.knn), coords, add=TRUE) title(main="K nearest neighbours, k = 4") data(state) us48.fipsno <- read.geoda(system.file("etc/weights/us48.txt", package="spdep")[1]) if (as.numeric(paste(version$major, version$minor, sep="")) < 19) { m50.48 <- match(us48.fipsno$"State.name", state.name) } else { m50.48 <- match(us48.fipsno$"State_name", state.name) } xy <- as.matrix(as.data.frame(state.center))[m50.48,] llk4.nb <- knn2nb(knearneigh(xy, k=4, longlat=FALSE)) gck4.nb <- knn2nb(knearneigh(xy, k=4, longlat=TRUE)) plot(llk4.nb, xy) plot(diffnb(llk4.nb, gck4.nb), xy, add=TRUE, col="red", lty=2) title(main="Differences between Euclidean and Great Circle k=4 neighbours") summary(llk4.nb, xy, longlat=TRUE, scale=0.5) summary(gck4.nb, xy, longlat=TRUE, scale=0.5) #xy1 <- SpatialPoints((as.data.frame(state.center))[m50.48,], # proj4string=CRS("+proj=longlat +ellps=GRS80")) #gck4a.nb <- knn2nb(knearneigh(xy1, k=4)) #summary(gck4a.nb, xy1, scale=0.5) xy1 <- st_as_sf((as.data.frame(state.center))[m50.48,], coords=1:2, crs=st_crs("OGC:CRS84")) old_use_s2 <- sf_use_s2() sf_use_s2(TRUE) system.time(gck4a.nb <- knn2nb(knearneigh(xy1, k=4))) summary(gck4a.nb, xy1, scale=0.5) sf_use_s2(FALSE) system.time(gck4a.nb <- knn2nb(knearneigh(xy1, k=4))) summary(gck4a.nb, xy1, scale=0.5) sf_use_s2(old_use_s2) # https://github.com/r-spatial/spdep/issues/38 run <- FALSE if (require("dbscan", quietly=TRUE)) run <- TRUE if (run) { set.seed(1) x <- cbind(runif(50), runif(50), runif(50)) out <- knearneigh(x, k=5) knn2nb(out) } if (run) { # fails because dbscan works with > 2 dimensions but does # not work with duplicate points try(out <- knearneigh(rbind(x, x[1:10,]), k=5)) }
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE) col.knn <- knearneigh(coords, k=4) plot(st_geometry(columbus), border="grey") plot(knn2nb(col.knn), coords, add=TRUE) title(main="K nearest neighbours, k = 4") data(state) us48.fipsno <- read.geoda(system.file("etc/weights/us48.txt", package="spdep")[1]) if (as.numeric(paste(version$major, version$minor, sep="")) < 19) { m50.48 <- match(us48.fipsno$"State.name", state.name) } else { m50.48 <- match(us48.fipsno$"State_name", state.name) } xy <- as.matrix(as.data.frame(state.center))[m50.48,] llk4.nb <- knn2nb(knearneigh(xy, k=4, longlat=FALSE)) gck4.nb <- knn2nb(knearneigh(xy, k=4, longlat=TRUE)) plot(llk4.nb, xy) plot(diffnb(llk4.nb, gck4.nb), xy, add=TRUE, col="red", lty=2) title(main="Differences between Euclidean and Great Circle k=4 neighbours") summary(llk4.nb, xy, longlat=TRUE, scale=0.5) summary(gck4.nb, xy, longlat=TRUE, scale=0.5) #xy1 <- SpatialPoints((as.data.frame(state.center))[m50.48,], # proj4string=CRS("+proj=longlat +ellps=GRS80")) #gck4a.nb <- knn2nb(knearneigh(xy1, k=4)) #summary(gck4a.nb, xy1, scale=0.5) xy1 <- st_as_sf((as.data.frame(state.center))[m50.48,], coords=1:2, crs=st_crs("OGC:CRS84")) old_use_s2 <- sf_use_s2() sf_use_s2(TRUE) system.time(gck4a.nb <- knn2nb(knearneigh(xy1, k=4))) summary(gck4a.nb, xy1, scale=0.5) sf_use_s2(FALSE) system.time(gck4a.nb <- knn2nb(knearneigh(xy1, k=4))) summary(gck4a.nb, xy1, scale=0.5) sf_use_s2(old_use_s2) # https://github.com/r-spatial/spdep/issues/38 run <- FALSE if (require("dbscan", quietly=TRUE)) run <- TRUE if (run) { set.seed(1) x <- cbind(runif(50), runif(50), runif(50)) out <- knearneigh(x, k=5) knn2nb(out) } if (run) { # fails because dbscan works with > 2 dimensions but does # not work with duplicate points try(out <- knearneigh(rbind(x, x[1:10,]), k=5)) }
The function converts a knn
object returned by knearneigh
into a neighbours list of class nb
with a list of integer vectors
containing neighbour region number ids.
knn2nb(knn, row.names = NULL, sym = FALSE)
knn2nb(knn, row.names = NULL, sym = FALSE)
knn |
A knn object returned by |
row.names |
character vector of region ids to be added to the neighbours list as attribute |
sym |
force the output neighbours list to symmetry |
The function returns an object of class nb
with a list of integer vectors containing neighbour region number ids. See card
for details of “nb” objects.
Roger Bivand [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_coordinates(st_centroid(columbus)) col.knn <- knearneigh(coords, k=4) plot(st_geometry(columbus), border="grey") plot(knn2nb(col.knn), coords, add=TRUE) title(main="K nearest neighbours, k = 4") # example of reading points with readr::read_csv() yielding a tibble load(system.file("etc/misc/coords.rda", package="spdep")) class(coords) knn2nb(knearneigh(coords, k=4))
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) coords <- st_coordinates(st_centroid(columbus)) col.knn <- knearneigh(coords, k=4) plot(st_geometry(columbus), border="grey") plot(knn2nb(col.knn), coords, add=TRUE) title(main="K nearest neighbours, k = 4") # example of reading points with readr::read_csv() yielding a tibble load(system.file("etc/misc/coords.rda", package="spdep")) class(coords) knn2nb(knearneigh(coords, k=4))
Using a listw
sparse representation of a spatial weights matrix, compute the lag vector
## S3 method for class 'listw' lag(x, var, zero.policy=attr(listw, "zero.policy"), NAOK=FALSE, ...)
## S3 method for class 'listw' lag(x, var, zero.policy=attr(listw, "zero.policy"), NAOK=FALSE, ...)
x |
a |
var |
a numeric vector the same length as the neighbours list in listw |
zero.policy |
default |
NAOK |
If 'FALSE', the presence of 'NA' values is regarded as an error; if 'TRUE' then any 'NA' or 'NaN' or 'Inf' values in var are represented as an NA lagged value. |
... |
additional arguments |
a numeric vector the same length as var
Roger Bivand [email protected]
data(oldcol) Vx <- lag.listw(nb2listw(COL.nb, style="W"), COL.OLD$CRIME) plot(Vx, COL.OLD$CRIME) plot(ecdf(COL.OLD$CRIME)) plot(ecdf(Vx), add=TRUE, col.points="red", col.hor="red") is.na(COL.OLD$CRIME[5]) <- TRUE VxNA <- lag.listw(nb2listw(COL.nb, style="W"), COL.OLD$CRIME, NAOK=TRUE)
data(oldcol) Vx <- lag.listw(nb2listw(COL.nb, style="W"), COL.OLD$CRIME) plot(Vx, COL.OLD$CRIME) plot(ecdf(COL.OLD$CRIME)) plot(ecdf(Vx), add=TRUE, col.points="red", col.hor="red") is.na(COL.OLD$CRIME[5]) <- TRUE VxNA <- lag.listw(nb2listw(COL.nb, style="W"), COL.OLD$CRIME, NAOK=TRUE)
A simple function to compute Lee's L statistic for bivariate spatial data;
lee(x, y, listw, n, S2, zero.policy=attr(listw, "zero.policy"), NAOK=FALSE)
lee(x, y, listw, n, S2, zero.policy=attr(listw, "zero.policy"), NAOK=FALSE)
x |
a numeric vector the same length as the neighbours list in listw |
y |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
n |
number of zones |
S2 |
Sum of squared sum of weights by rows. |
zero.policy |
default |
NAOK |
if 'TRUE' then any 'NA' or 'NaN' or 'Inf' values in x are passed on to the foreign function. If 'FALSE', the presence of 'NA' or 'NaN' or 'Inf' values is regarded as an error. |
a list of
L |
Lee's L statistic |
local L |
Lee's local L statistic |
Roger Bivand and Virgiio Gómez-Rubio [email protected]
Lee (2001). Developing a bivariate spatial association measure: An integration of Pearson's r and Moran's I. J Geograph Syst 3: 369-385
data(boston, package="spData") lw<-nb2listw(boston.soi) x<-boston.c$CMEDV y<-boston.c$CRIM z<-boston.c$RAD Lxy<-lee(x, y, lw, length(x), zero.policy=TRUE) Lxz<-lee(x, z, lw, length(x), zero.policy=TRUE)
data(boston, package="spData") lw<-nb2listw(boston.soi) x<-boston.c$CMEDV y<-boston.c$CRIM z<-boston.c$RAD Lxy<-lee(x, y, lw, length(x), zero.policy=TRUE) Lxz<-lee(x, z, lw, length(x), zero.policy=TRUE)
A permutation test for Lee's L statistic calculated by using nsim random permutations of x and y for the given spatial weighting scheme, to establish the rank of the observed statistic in relation to the nsim simulated values.
lee.mc(x, y, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", na.action=na.fail, spChk=NULL, return_boot=FALSE)
lee.mc(x, y, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", na.action=na.fail, spChk=NULL, return_boot=FALSE)
x |
a numeric vector the same length as the neighbours list in listw |
y |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
nsim |
number of permutations |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "two.sided", or "less". |
na.action |
a function (default |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
return_boot |
return an object of class |
A list with class htest
and mc.sim
containing the following components:
statistic |
the value of the observed Lee's L. |
parameter |
the rank of the observed Lee's L. |
p.value |
the pseudo p-value of the test. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data, and the number of simulations. |
res |
nsim simulated values of statistic, final value is observed statistic |
Roger Bivand, Virgilio Gómez-Rubio [email protected]
Lee (2001). Developing a bivariate spatial association measure: An integration of Pearson's r and Moran's I. J Geograph Syst 3: 369-385
data(boston, package="spData") lw<-nb2listw(boston.soi) x<-boston.c$CMEDV y<-boston.c$CRIM lee.mc(x, y, nsim=99, lw, zero.policy=TRUE, alternative="two.sided") #Test with missing values x[1:5]<-NA y[3:7]<-NA lee.mc(x, y, nsim=99, lw, zero.policy=TRUE, alternative="two.sided", na.action=na.omit)
data(boston, package="spData") lw<-nb2listw(boston.soi) x<-boston.c$CMEDV y<-boston.c$CRIM lee.mc(x, y, nsim=99, lw, zero.policy=TRUE, alternative="two.sided") #Test with missing values x[1:5]<-NA y[3:7]<-NA lee.mc(x, y, nsim=99, lw, zero.policy=TRUE, alternative="two.sided", na.action=na.omit)
Lee's L test for spatial autocorrelation using a spatial weights matrix in weights list form. The assumptions underlying the test are sensitive to the form of the graph of neighbour relationships and other factors, and results may be checked against those of lee.mc
permutations.
lee.test(x, y, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", na.action=na.fail, spChk=NULL)
lee.test(x, y, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", na.action=na.fail, spChk=NULL)
x |
a numeric vector the same length as the neighbours list in listw |
y |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of greater (default), less or two.sided. |
na.action |
a function (default |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
A list with class htest
containing the following components:
statistic |
the value of the standard deviate of Lee's L. |
p.value |
the p-value of the test. |
estimate |
the value of the observed Lee's L, its expectation and variance under the method assumption. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the assumption used for calculating the standard deviate. |
data.name |
a character string giving the name(s) of the data. |
See Lee (2004) for details on how the asymptotic expectation and variance of Lee's L is computed. In particular, check Lee (2004), table 1, page 1690.
This test may fail for large datasets as the computation of the asymptotic expectation and variance requires the use of dense matrices.
Roger Bivand and Virgilio Gómez-Rubio [email protected]
Lee (2004). A generalized significance testing method for global measures of spatial association: an extension of the Mantel test. Environment and Planning A 2004, volume 36, pages 1687 - 1703
data(oldcol) col.W <- nb2listw(COL.nb, style="W") crime <- COL.OLD$CRIME lee.test(crime, crime, col.W, zero.policy=TRUE) #Test with missing values x<-crime y<-crime x[1:5]<-NA y[3:7]<-NA lee.test(x, y, col.W, zero.policy=TRUE, na.action=na.omit) # lee.test(x, y, col.W, zero.policy=TRUE)#Stops with an error data(boston, package="spData") lw<-nb2listw(boston.soi) x<-boston.c$CMEDV y<-boston.c$CRIM lee.test(x, y, lw, zero.policy=TRUE, alternative="less") #Test with missing values x[1:5]<-NA y[3:7]<-NA lee.test(x, y, lw, zero.policy=TRUE, alternative="less", na.action=na.omit)
data(oldcol) col.W <- nb2listw(COL.nb, style="W") crime <- COL.OLD$CRIME lee.test(crime, crime, col.W, zero.policy=TRUE) #Test with missing values x<-crime y<-crime x[1:5]<-NA y[3:7]<-NA lee.test(x, y, col.W, zero.policy=TRUE, na.action=na.omit) # lee.test(x, y, col.W, zero.policy=TRUE)#Stops with an error data(boston, package="spData") lw<-nb2listw(boston.soi) x<-boston.c$CMEDV y<-boston.c$CRIM lee.test(x, y, lw, zero.policy=TRUE, alternative="less") #Test with missing values x[1:5]<-NA y[3:7]<-NA lee.test(x, y, lw, zero.policy=TRUE, alternative="less", na.action=na.omit)
Local indicators for categorical data combine a measure of local composition in a window given by the per-observation set of neighbouring observations, with a local multi-category joincount test simplified to neighbours with the same or different categories compared to the focal observation
licd_multi(fx, listw, zero.policy = attr(listw, "zero.policy"), adjust.n = TRUE, nsim = 0L, iseed = NULL, no_repeat_in_row = FALSE, control = list())
licd_multi(fx, listw, zero.policy = attr(listw, "zero.policy"), adjust.n = TRUE, nsim = 0L, iseed = NULL, no_repeat_in_row = FALSE, control = list())
fx |
a factor with two or more categories, of the same length as the neighbours and weights objects in listw |
listw |
a |
zero.policy |
default |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted |
nsim |
default 0, number of conditonal permutation simulations |
iseed |
default NULL, used to set the seed; the output will only be reproducible if the count of CPU cores across which computation is distributed is the same |
no_repeat_in_row |
default |
control |
comp_binary= |
The original code may be found at doi:10.5281/zenodo.4283766
local_comp |
data.frame object with LICD local composition columns: ID, category_i, count_like_i, prop_i, count_nbs_i, pbinom_like_BW, pbinom_unlike_BW, pbinom_unlike_BW_alt, chi_BW_i, chi_K_i, anscombe_BW |
local_config |
data.frame object with LICD local configuration columns: ID, jcm_chi_obs, jcm_count_BB_obs, jcm_count_BW_obs, jcm_count_WW_obs, pval_jcm_obs_BB, pval_jcm_obs_WW, pval_jcm_obs_BW |
local_comp_sim |
data.frame object with permutation-based LICD local composition columns: ID, pbinom_like_BW, pbinom_unlike_BW, pbinom_unlike_BW_alt, rank_sim_chi_BW, rank_sim_chi_K, rank_sim_anscombe |
local_config_sim |
data.frame object with permutation-based LICD local configuration columns: ID, jcm_chi_sim_rank, pval_jcm_obs_BB, pval_jcm_obs_BW, pval_jcm_obs_WW |
In order to increase the numbers of neighbours using nblag
and nblag_cumul
is advisable; use of binary weights is advised and are in any case used for the composition measure
Roger Bivand [email protected] based on earlier code by Tomasz M. Kossowski, Justyna Wilk and Michał B. Pietrzak
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 20;
Upton, G., Fingleton, B. 1985 Spatial data analysis by example: point pattern and qualitative data, Wiley, pp. 158–170;
Boots, B., 2003. Developing local measures of spatial association for categorical data. Journal of Geographical Systems 5, 139–160;
Boots, B., 2006. Local configuration measures for categorical spatial data: binary regular lattices. Journal of Geographical Systems 8 (1), 1–24;
Pietrzak, M.B., Wilk, J., Kossowski, T., Bivand, R.S., 2014. The application of local indicators for categorical data (LICD) in the spatial analysis of economic development. Comparative Economic Research 17 (4), 203–220 doi:10.2478/cer-2014-0041;
Bivand, R.S., Wilk, J., Kossowski, T., 2017. Spatial association of population pyramids across Europe: The application of symbolic data, cluster analysis and join-count tests. Spatial Statistics 21 (B), 339–361 doi:10.1016/j.spasta.2017.03.003;
Francesco Carrer, Tomasz M. Kossowski, Justyna Wilk, Michał B. Pietrzak, Roger S. Bivand, The application of Local Indicators for Categorical Data (LICD) to explore spatial dependence in archaeological spaces, Journal of Archaeological Science, 126, 2021, doi:10.1016/j.jas.2020.105306
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) (nb <- poly2nb(columbus)) lw <- nb2listw(nblag_cumul(nblag(nb, 2)), style="B") obj <- licd_multi(HICRIME, lw) str(obj) h_obj <- hotspot(obj) str(h_obj) table(h_obj$both_recode) columbus$both <- h_obj$both_recode plot(columbus[, "both"]) GDAL37 <- as.numeric_version(unname(sf_extSoftVersion()["GDAL"])) >= "3.7.0" file <- "etc/shapes/GB_2024_southcoast_50m.gpkg.zip" zipfile <- system.file(file, package="spdep") if (GDAL37) { sc50m <- st_read(zipfile) } else { td <- tempdir() bn <- sub(".zip", "", basename(file), fixed=TRUE) target <- unzip(zipfile, files=bn, exdir=td) sc50m <- st_read(target) } sc50m$Winner <- factor(sc50m$Winner, levels=c("Con", "Green", "Lab", "LD")) plot(sc50m[,"Winner"], pal=c("#2297E6", "#61D04F", "#DF536B", "#F5C710")) nb_sc_50m <- poly2nb(sc50m, row.names=as.character(sc50m$Constituency)) sub2 <- attr(nb_sc_50m, "region.id")[attr(nb_sc_50m, "ncomp")$comp.id == 2L] iowe <- match(sub2[1], attr(nb_sc_50m, "region.id")) diowe <- c(st_distance(sc50m[iowe,], sc50m)) meet_criterion <- sum(diowe <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(diowe)[1:meet_criterion]] nb_sc_50m_iowe <- addlinks1(nb_sc_50m, from = cands[1], to = cands[3:meet_criterion]) ioww <- match(sub2[2], attr(nb_sc_50m, "region.id")) dioww <- c(st_distance(sc50m[ioww,], sc50m)) meet_criterion <- sum(dioww <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(dioww)[1:meet_criterion]] nb_sc_50m_iow <- addlinks1(nb_sc_50m_iowe, from = cands[2], to = cands[3:meet_criterion]) nb_sc_1_2 <- nblag_cumul(nblag(nb_sc_50m_iow, 2)) lw <- nb2listw(nb_sc_1_2, style="B") licd_obj <- licd_multi(sc50m$Winner, lw) h_obj <- hotspot(licd_obj) sc50m$both <- h_obj$both_recode plot(sc50m[, "both"]) ljc <- local_joincount_uni(factor(sc50m$Winner == "LD"), chosen="TRUE", lw) sc50m$LD_pv <- ljc[, 2] plot(sc50m[, "LD_pv"], breaks=c(0, 0.025, 0.05, 0.1, 0.5, 1))
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) (nb <- poly2nb(columbus)) lw <- nb2listw(nblag_cumul(nblag(nb, 2)), style="B") obj <- licd_multi(HICRIME, lw) str(obj) h_obj <- hotspot(obj) str(h_obj) table(h_obj$both_recode) columbus$both <- h_obj$both_recode plot(columbus[, "both"]) GDAL37 <- as.numeric_version(unname(sf_extSoftVersion()["GDAL"])) >= "3.7.0" file <- "etc/shapes/GB_2024_southcoast_50m.gpkg.zip" zipfile <- system.file(file, package="spdep") if (GDAL37) { sc50m <- st_read(zipfile) } else { td <- tempdir() bn <- sub(".zip", "", basename(file), fixed=TRUE) target <- unzip(zipfile, files=bn, exdir=td) sc50m <- st_read(target) } sc50m$Winner <- factor(sc50m$Winner, levels=c("Con", "Green", "Lab", "LD")) plot(sc50m[,"Winner"], pal=c("#2297E6", "#61D04F", "#DF536B", "#F5C710")) nb_sc_50m <- poly2nb(sc50m, row.names=as.character(sc50m$Constituency)) sub2 <- attr(nb_sc_50m, "region.id")[attr(nb_sc_50m, "ncomp")$comp.id == 2L] iowe <- match(sub2[1], attr(nb_sc_50m, "region.id")) diowe <- c(st_distance(sc50m[iowe,], sc50m)) meet_criterion <- sum(diowe <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(diowe)[1:meet_criterion]] nb_sc_50m_iowe <- addlinks1(nb_sc_50m, from = cands[1], to = cands[3:meet_criterion]) ioww <- match(sub2[2], attr(nb_sc_50m, "region.id")) dioww <- c(st_distance(sc50m[ioww,], sc50m)) meet_criterion <- sum(dioww <= units::set_units(5000, "m")) cands <- attr(nb_sc_50m, "region.id")[order(dioww)[1:meet_criterion]] nb_sc_50m_iow <- addlinks1(nb_sc_50m_iowe, from = cands[2], to = cands[3:meet_criterion]) nb_sc_1_2 <- nblag_cumul(nblag(nb_sc_50m_iow, 2)) lw <- nb2listw(nb_sc_1_2, style="B") licd_obj <- licd_multi(sc50m$Winner, lw) h_obj <- hotspot(licd_obj) sc50m$both <- h_obj$both_recode plot(sc50m[, "both"]) ljc <- local_joincount_uni(factor(sc50m$Winner == "LD"), chosen="TRUE", lw) sc50m$LD_pv <- ljc[, 2] plot(sc50m[, "LD_pv"], breaks=c(0, 0.025, 0.05, 0.1, 0.5, 1))
The function makes a "spatial neighbour"
object representation
(similar to the S-PLUS spatial statististics module representation of a
"listw"
spatial weights object. sn2listw()
is the inverse function to listw2sn()
, creating a "listw"
object from a "spatial neighbour"
object.
listw2sn(listw) sn2listw(sn, style = NULL, zero.policy = NULL, from_mat2listw=FALSE)
listw2sn(listw) sn2listw(sn, style = NULL, zero.policy = NULL, from_mat2listw=FALSE)
listw |
a |
sn |
a |
style |
default NULL, missing, set to "M" and warning given; if not "M", passed to |
zero.policy |
default NULL, use global option value; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors |
from_mat2listw |
default FALSE, set TRUE if called from |
listw2sn()
returns a data frame with three columns, and with class spatial.neighbour
:
from |
region number id for the start of the link (S-PLUS row.id) |
to |
region number id for the end of the link (S-PLUS col.id) |
weights |
weight for this link |
Roger Bivand [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) col.listw <- nb2listw(col.gal.nb) col.listw$neighbours[[1]] col.listw$weights[[1]] col.sn <- listw2sn(col.listw) str(col.sn)
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1]) col.listw <- nb2listw(col.gal.nb) col.listw$neighbours[[1]] col.listw$weights[[1]] col.sn <- listw2sn(col.listw) str(col.sn)
The function reports the estimates of tests chosen among five statistics for
testing for spatial dependence in linear models. The statistics are
the simple RS test for error dependence (“RSerr”), the simple RS test
for a missing spatially lagged dependent variable (“RSlag”), variants
of these adjusted for the presence of the other (“adjRSerr”
tests for error dependence in the possible presence of a missing lagged
dependent variable, “adjRSlag” the other way round), and a portmanteau test
(“SARMA”, in fact “RSerr” + “adjRSlag”). Note: from spdep 1.3-2, the tests are re-named “RS” - Rao's score tests, rather than “LM” - Lagrange multiplier tests to match the naming of tests from the same family in SDM.RStests
.
lm.RStests(model, listw, zero.policy=attr(listw, "zero.policy"), test="RSerr", spChk=NULL, naSubset=TRUE) lm.LMtests(model, listw, zero.policy=attr(listw, "zero.policy"), test="LMerr", spChk=NULL, naSubset=TRUE) ## S3 method for class 'RStestlist' print(x, ...) ## S3 method for class 'RStestlist' summary(object, p.adjust.method="none", ...) ## S3 method for class 'RStestlist.summary' print(x, digits=max(3, getOption("digits") - 2), ...)
lm.RStests(model, listw, zero.policy=attr(listw, "zero.policy"), test="RSerr", spChk=NULL, naSubset=TRUE) lm.LMtests(model, listw, zero.policy=attr(listw, "zero.policy"), test="LMerr", spChk=NULL, naSubset=TRUE) ## S3 method for class 'RStestlist' print(x, ...) ## S3 method for class 'RStestlist' summary(object, p.adjust.method="none", ...) ## S3 method for class 'RStestlist.summary' print(x, digits=max(3, getOption("digits") - 2), ...)
model |
an object of class |
listw |
a |
zero.policy |
default |
test |
a character vector of tests requested chosen from RSerr, RSlag, adjRSerr, adjRSlag, SARMA; test="all" computes all the tests. |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
naSubset |
default TRUE to subset listw object for omitted observations in model object (this is a change from earlier behaviour, when the |
x , object
|
object to be printed |
p.adjust.method |
a character string specifying the probability value adjustment (see |
digits |
minimum number of significant digits to be used for most numbers |
... |
printing arguments to be passed through |
The two types of dependence are for spatial lag and spatial error
:
where is a well-behaved, uncorrelated error
term. Tests for a missing spatially lagged dependent variable test
that
, tests for spatial autocorrelation of
the error
test whether
.
is a spatial weights matrix; for the tests used
here they are identical.
A list of class RStestlist
of htest
objects, each with:
statistic |
the value of the Rao's score (a.k.a Lagrange multiplier) test. |
parameter |
number of degrees of freedom |
p.value |
the p-value of the test. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data. |
Roger Bivand [email protected] and Andrew Bernat
Anselin, L. 1988 Spatial econometrics: methods and models. (Dordrecht: Kluwer); Anselin, L., Bera, A. K., Florax, R. and Yoon, M. J. 1996 Simple diagnostic tests for spatial dependence. Regional Science and Urban Economics, 26, 77–104 doi:10.1016/0166-0462(95)02111-6; Malabika Koley (2024) Specification Testing under General Nesting Spatial Model, https://sites.google.com/view/malabikakoley/research.
data(oldcol) oldcrime.lm <- lm(CRIME ~ HOVAL + INC, data = COL.OLD) summary(oldcrime.lm) lw <- nb2listw(COL.nb) res <- lm.RStests(oldcrime.lm, listw=lw, test="all") summary(res) if (require("spatialreg", quietly=TRUE)) { oldcrime.slx <- lmSLX(CRIME ~ HOVAL + INC, data = COL.OLD, listw=lw) summary(lm.RStests(oldcrime.slx, listw=lw, test=c("adjRSerr", "adjRSlag"))) }
data(oldcol) oldcrime.lm <- lm(CRIME ~ HOVAL + INC, data = COL.OLD) summary(oldcrime.lm) lw <- nb2listw(COL.nb) res <- lm.RStests(oldcrime.lm, listw=lw, test="all") summary(res) if (require("spatialreg", quietly=TRUE)) { oldcrime.slx <- lmSLX(CRIME ~ HOVAL + INC, data = COL.OLD, listw=lw) summary(lm.RStests(oldcrime.slx, listw=lw, test=c("adjRSerr", "adjRSlag"))) }
Moran's I test for spatial autocorrelation in residuals from an estimated linear model (lm()
).
lm.morantest(model, listw, zero.policy=attr(listw, "zero.policy"), alternative = "greater", spChk=NULL, resfun=weighted.residuals, naSubset=TRUE)
lm.morantest(model, listw, zero.policy=attr(listw, "zero.policy"), alternative = "greater", spChk=NULL, resfun=weighted.residuals, naSubset=TRUE)
model |
an object of class |
listw |
a |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two.sided". |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
resfun |
default: weighted.residuals; the function to be used to extract residuals from the |
naSubset |
default TRUE to subset listw object for omitted observations in model object (this is a change from earlier behaviour, when the |
A list with class htest
containing the following components:
statistic |
the value of the standard deviate of Moran's I. |
p.value |
the p-value of the test. |
estimate |
the value of the observed Moran's I, its expectation and variance under the method assumption. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data. |
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 203,
data(oldcol) oldcrime1.lm <- lm(CRIME ~ 1, data = COL.OLD) oldcrime.lm <- lm(CRIME ~ HOVAL + INC, data = COL.OLD) lm.morantest(oldcrime.lm, nb2listw(COL.nb, style="W")) lm.LMtests(oldcrime.lm, nb2listw(COL.nb, style="W")) lm.morantest(oldcrime.lm, nb2listw(COL.nb, style="S")) lm.morantest(oldcrime1.lm, nb2listw(COL.nb, style="W")) moran.test(COL.OLD$CRIME, nb2listw(COL.nb, style="W"), randomisation=FALSE) oldcrime.wlm <- lm(CRIME ~ HOVAL + INC, data = COL.OLD, weights = I(1/AREA_PL)) lm.morantest(oldcrime.wlm, nb2listw(COL.nb, style="W"), resfun=weighted.residuals) lm.morantest(oldcrime.wlm, nb2listw(COL.nb, style="W"), resfun=rstudent)
data(oldcol) oldcrime1.lm <- lm(CRIME ~ 1, data = COL.OLD) oldcrime.lm <- lm(CRIME ~ HOVAL + INC, data = COL.OLD) lm.morantest(oldcrime.lm, nb2listw(COL.nb, style="W")) lm.LMtests(oldcrime.lm, nb2listw(COL.nb, style="W")) lm.morantest(oldcrime.lm, nb2listw(COL.nb, style="S")) lm.morantest(oldcrime1.lm, nb2listw(COL.nb, style="W")) moran.test(COL.OLD$CRIME, nb2listw(COL.nb, style="W"), randomisation=FALSE) oldcrime.wlm <- lm(CRIME ~ HOVAL + INC, data = COL.OLD, weights = I(1/AREA_PL)) lm.morantest(oldcrime.wlm, nb2listw(COL.nb, style="W"), resfun=weighted.residuals) lm.morantest(oldcrime.wlm, nb2listw(COL.nb, style="W"), resfun=rstudent)
The function implements Tiefelsdorf's exact global Moran's I test.
lm.morantest.exact(model, listw, zero.policy = attr(listw, "zero.policy"), alternative = "greater", spChk = NULL, resfun = weighted.residuals, zero.tol = 1e-07, Omega=NULL, save.M=NULL, save.U=NULL, useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) ## S3 method for class 'moranex' print(x, ...)
lm.morantest.exact(model, listw, zero.policy = attr(listw, "zero.policy"), alternative = "greater", spChk = NULL, resfun = weighted.residuals, zero.tol = 1e-07, Omega=NULL, save.M=NULL, save.U=NULL, useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) ## S3 method for class 'moranex' print(x, ...)
model |
an object of class |
listw |
a |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of greater (default), less or two.sided. |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
resfun |
default: weighted.residuals; the function to be used to extract residuals from the |
zero.tol |
tolerance used to find eigenvalues close to absolute zero |
Omega |
A SAR process matrix may be passed in to test an alternative hypothesis, for example |
save.M |
return the full M matrix for use in |
save.U |
return the full U matrix for use in |
useTP |
default FALSE, if TRUE, use truncation point in integration rather than upper=Inf, see Tiefelsdorf (2000), eq. 6.7, p.69 |
truncErr |
when useTP=TRUE, pass truncation error to truncation point function |
zeroTreat |
when useTP=TRUE, pass zero adjustment to truncation point function |
x |
a moranex object |
... |
arguments to be passed through |
A list of class moranex
with the following components:
statistic |
the value of the exact standard deviate of global Moran's I. |
p.value |
the p-value of the test. |
estimate |
the value of the observed global Moran's I. |
method |
a character string giving the method used. |
alternative |
a character string describing the alternative hypothesis. |
gamma |
eigenvalues (excluding zero values) |
oType |
usually set to "E" |
data.name |
a character string giving the name(s) of the data. |
df |
degrees of freedom |
Markus Reder and Roger Bivand
Roger Bivand, Werner G. Müller and Markus Reder (2009) "Power calculations for global and local Moran's I." Computational Statistics & Data Analysis 53, 2859-2872.
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) lm.morantest(e.lm, nb2listw(eire.nb)) lm.morantest.sad(e.lm, nb2listw(eire.nb)) lm.morantest.exact(e.lm, nb2listw(eire.nb)) lm.morantest.exact(e.lm, nb2listw(eire.nb), useTP=TRUE)
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) lm.morantest(e.lm, nb2listw(eire.nb)) lm.morantest.sad(e.lm, nb2listw(eire.nb)) lm.morantest.exact(e.lm, nb2listw(eire.nb)) lm.morantest.exact(e.lm, nb2listw(eire.nb), useTP=TRUE)
The function implements Tiefelsdorf's application of the Saddlepoint approximation to global Moran's I's reference distribution.
lm.morantest.sad(model, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, resfun=weighted.residuals, tol=.Machine$double.eps^0.5, maxiter=1000, tol.bounds=0.0001, zero.tol = 1e-07, Omega=NULL, save.M=NULL, save.U=NULL) ## S3 method for class 'moransad' print(x, ...) ## S3 method for class 'moransad' summary(object, ...) ## S3 method for class 'summary.moransad' print(x, ...)
lm.morantest.sad(model, listw, zero.policy=attr(listw, "zero.policy"), alternative="greater", spChk=NULL, resfun=weighted.residuals, tol=.Machine$double.eps^0.5, maxiter=1000, tol.bounds=0.0001, zero.tol = 1e-07, Omega=NULL, save.M=NULL, save.U=NULL) ## S3 method for class 'moransad' print(x, ...) ## S3 method for class 'moransad' summary(object, ...) ## S3 method for class 'summary.moransad' print(x, ...)
model |
an object of class |
listw |
a |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of greater (default), less or two.sided. |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
resfun |
default: weighted.residuals; the function to be used to extract residuals from the |
tol |
the desired accuracy (convergence tolerance) for |
maxiter |
the maximum number of iterations for |
tol.bounds |
offset from bounds for |
zero.tol |
tolerance used to find eigenvalues close to absolute zero |
Omega |
A SAR process matrix may be passed in to test an alternative hypothesis, for example |
save.M |
return the full M matrix for use in |
save.U |
return the full U matrix for use in |
x |
object to be printed |
object |
object to be summarised |
... |
arguments to be passed through |
The function involves finding the eigenvalues of an n by n matrix, and numerically finding the root for the Saddlepoint approximation, and should therefore only be used with care when n is large.
A list of class moransad
with the following components:
statistic |
the value of the saddlepoint approximation of the standard deviate of global Moran's I. |
p.value |
the p-value of the test. |
estimate |
the value of the observed global Moran's I. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data. |
internal1 |
Saddlepoint omega, r and u |
internal2 |
f.root, iter and estim.prec from |
df |
degrees of freedom |
tau |
eigenvalues (excluding zero values) |
Roger Bivand [email protected]
Tiefelsdorf, M. 2002 The Saddlepoint approximation of Moran's I and local Moran's Ii reference distributions and their numerical evaluation. Geographical Analysis, 34, pp. 187–206; Bivand RS, Wong DWS 2018 Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748 doi:10.1007/s11749-018-0599-x
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) lm.morantest(e.lm, nb2listw(eire.nb)) lm.morantest.sad(e.lm, nb2listw(eire.nb)) summary(lm.morantest.sad(e.lm, nb2listw(eire.nb))) e.wlm <- lm(OWNCONS ~ ROADACC, data=eire, weights=RETSALE) lm.morantest(e.wlm, nb2listw(eire.nb), resfun=rstudent) lm.morantest.sad(e.wlm, nb2listw(eire.nb), resfun=rstudent)
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) lm.morantest(e.lm, nb2listw(eire.nb)) lm.morantest.sad(e.lm, nb2listw(eire.nb)) summary(lm.morantest.sad(e.lm, nb2listw(eire.nb))) e.wlm <- lm(OWNCONS ~ ROADACC, data=eire, weights=RETSALE) lm.morantest(e.wlm, nb2listw(eire.nb), resfun=rstudent) lm.morantest.sad(e.wlm, nb2listw(eire.nb), resfun=rstudent)
The bivariate join count (BJC) evaluates event occurrences in predefined regions and tests if the co-occurrence of events deviates from complete spatial randomness.
local_joincount_bv( x, z, listw, nsim = 199, alternative = "two.sided" )
local_joincount_bv( x, z, listw, nsim = 199, alternative = "two.sided" )
x |
a binary variable either numeric or logical |
z |
a binary variable either numeric or logical with the same length as |
listw |
a listw object containing binary weights created, for example, with |
nsim |
the number of conditional permutation simulations |
alternative |
default |
There are two cases that are evaluated in the bivariate join count. The first being in-situ colocation (CLC) where xi = 1 and zi = 1. The second is the general form of the bivariate join count (BJC) that is used when there is no in-situ colocation.
The BJC case "is useful when x and z cannot occur in the same location, such as when x and z correspond to two different values of a single categorical variable" or "when x and z can co-locate, but do not" (Anselin and Li, 2019). Whereas the CLC case is useful in evaluating simultaneous occurrences of events.
The local bivariate join count statistic requires a binary weights list which can be generated with nb2listw(nb, style = "B")
.
P-values are only reported for those regions that match the CLC or BJC criteria. Others will not have an associated p-value.
P-values are estimated using a conditional permutation approach. This creates a reference distribution from which the observed statistic is compared.
a data.frame
with two columns join_count
and p_sim
and number of rows equal to the length of arguments x
.
Josiah Parry [email protected]
Anselin, L., & Li, X. (2019). Operational Local Join Count Statistics for Cluster Detection. Journal of geographical systems, 21(2), 189–210. doi:10.1007/s10109-019-00299-x
data("oldcol") listw <- nb2listw(COL.nb, style = "B") # Colocation case x <- COL.OLD[["CP"]] z <- COL.OLD[["EW"]] set.seed(1) res <- local_joincount_bv(x, z, listw) na.omit(res) # no colocation case z <- 1 - x set.seed(1) res <- local_joincount_bv(x, z, listw) na.omit(res)
data("oldcol") listw <- nb2listw(COL.nb, style = "B") # Colocation case x <- COL.OLD[["CP"]] z <- COL.OLD[["EW"]] set.seed(1) res <- local_joincount_bv(x, z, listw) na.omit(res) # no colocation case z <- 1 - x set.seed(1) res <- local_joincount_bv(x, z, listw) na.omit(res)
The univariate local join count statistic is used to identify clusters of rarely occurring binary variables. The binary variable of interest should occur less than half of the time.
local_joincount_uni( fx, chosen, listw, alternative = "two.sided", nsim = 199, iseed = NULL, no_repeat_in_row=FALSE )
local_joincount_uni( fx, chosen, listw, alternative = "two.sided", nsim = 199, iseed = NULL, no_repeat_in_row=FALSE )
fx |
a binary variable either numeric or logical |
chosen |
a scalar character containing the level of |
listw |
a listw object containing binary weights created, for example, with |
alternative |
default |
nsim |
the number of conditional permutation simulations |
iseed |
default NULL, used to set the seed; the output will only be reproducible if the count of CPU cores across which computation is distributed is the same |
no_repeat_in_row |
default |
The local join count statistic requires a binary weights list which can be generated with nb2listw(nb, style = "B")
. Additionally, ensure that the binary variable of interest is rarely occurring in no more than half of observations.
P-values are estimated using a conditional permutation approach. This creates a reference distribution from which the observed statistic is compared. For more see Geoda Glossary.
a data.frame
with two columns BB
and Pr()
and number of rows equal to the length of x
.
Josiah Parry [email protected]
Anselin, L., & Li, X. (2019). Operational Local Join Count Statistics for Cluster Detection. Journal of geographical systems, 21(2), 189–210. doi:10.1007/s10109-019-00299-x
data(oldcol) fx <- as.factor(ifelse(COL.OLD$CRIME < 35, "low-crime", "high-crime")) listw <- nb2listw(COL.nb, style = "B") set.seed(1) (res <- local_joincount_uni(fx, chosen = "high-crime", listw))
data(oldcol) fx <- as.factor(ifelse(COL.OLD$CRIME < 35, "low-crime", "high-crime")) listw <- nb2listw(COL.nb, style = "B") set.seed(1) (res <- local_joincount_uni(fx, chosen = "high-crime", listw))
The Local Geary is a local adaptation of Geary's C statistic of spatial autocorrelation. The Local Geary uses squared differences to measure dissimilarity unlike the Local Moran. Low values of the Local Geary indicate positive spatial autocorrelation and large refers to negative spatial autocorrelation.
Inference for the Local Geary is based on a permutation approach which compares the observed value to the reference distribution under spatial randomness. localC_perm()
returns a pseudo p-value. This is not an analytical p-value and is based on the number of permutations and as such should be used with care.
localC(x, ..., zero.policy=NULL) ## Default S3 method: localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'formula' localC(formula, data, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'list' localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'matrix' localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'data.frame' localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) localC_perm(x, ..., zero.policy=NULL, iseed=NULL, no_repeat_in_row=FALSE) ## Default S3 method: localC_perm(x, listw, nsim = 499, alternative = "two.sided", ..., zero.policy=attr(listw, "zero.policy"), iseed=NULL, no_repeat_in_row=FALSE) ## S3 method for class 'formula' localC_perm(formula, data, listw, nsim = 499, alternative = "two.sided", ..., zero.policy=attr(listw, "zero.policy"), iseed=NULL, no_repeat_in_row=FALSE)
localC(x, ..., zero.policy=NULL) ## Default S3 method: localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'formula' localC(formula, data, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'list' localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'matrix' localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) ## S3 method for class 'data.frame' localC(x, listw, ..., zero.policy=attr(listw, "zero.policy")) localC_perm(x, ..., zero.policy=NULL, iseed=NULL, no_repeat_in_row=FALSE) ## Default S3 method: localC_perm(x, listw, nsim = 499, alternative = "two.sided", ..., zero.policy=attr(listw, "zero.policy"), iseed=NULL, no_repeat_in_row=FALSE) ## S3 method for class 'formula' localC_perm(formula, data, listw, nsim = 499, alternative = "two.sided", ..., zero.policy=attr(listw, "zero.policy"), iseed=NULL, no_repeat_in_row=FALSE)
x |
a numeric vector, numeric matrix, or list. See details for more. |
formula |
A one-sided formula determining which variables to be used. |
listw |
a |
data |
Used when a formula is provided. A matrix or data frame containing the variables in the formula |
nsim |
The number of simulations to be used for permutation test. |
alternative |
A character defining the alternative hypothesis. Must be one of |
... |
other arguments passed to methods. |
zero.policy |
default |
iseed |
default NULL, used to set the seed;the output will only be reproducible if the count of CPU cores across which computation is distributed is the same |
no_repeat_in_row |
default |
The Local Geary can be extended to a multivariate context. When x
is a numeric vector, the univariate Local Geary will be calculated. To calculate the multivariate Local Moran provide either a list or a matrix. When x
is a list, each element must be a numeric vector of the same length and of the same length as the neighbours in listw
. In the case that x
is a matrix the number of rows must be the same as the length of the neighbours in listw
.
While not required in the univariate context, the standardized Local Geary is calculated. The multivariate Local Geary is always standardized.
The univariate Local Geary is calculated as and the multivariate Local Geary is calculated as
as described in Anselin (2019).
A numeric vector containing Local Geary statistic with attribute pseudo-p
when localC_perm()
is used. pseudo-p
is an 8 column matrix containing
E.Ci |
expectation of the Local Geary statistic based on permutation sample |
Var.Ci |
variance of Local Geary based on permutation sample |
Z.Ci |
standard deviate of Local Geary based on permutation sample |
Pr() |
p-value of Local Geary statistic using |
Pr() Sim |
|
Pr(folded) Sim |
the simulation folded [0, 0.5] range ranked p-value (based on https://github.com/pysal/esda/blob/4a63e0b5df1e754b17b5f1205b8cadcbecc5e061/esda/crand.py#L211-L213) |
Skewness |
the output of |
Kurtosis |
the output of |
Josiah Parry [email protected] and Roger Bivand
Anselin, L. (1995), Local Indicators of Spatial Association—LISA. Geographical Analysis, 27: 93-115. doi:10.1111/j.1538-4632.1995.tb00338.x
Anselin, L. (2019), A Local Indicator of Multivariate Spatial Association: Extending Gearys c. Geogr Anal, 51: 133-150. doi:10.1111/gean.12164
orig <- spData::africa.rook.nb listw <- nb2listw(orig) x <- spData::afcon$totcon (A <- localC(x, listw)) listw1 <- nb2listw(droplinks(sym.attr.nb(orig), 3, sym=TRUE), zero.policy=TRUE) (A1 <- localC(x, listw1, zero.policy=FALSE)) (A2 <- localC(x, listw1, zero.policy=TRUE)) run <- FALSE if (require(rgeoda, quietly=TRUE)) run <- TRUE if (run) { W <- create_weights(as.numeric(length(x))) for (i in 1:length(listw$neighbours)) { set_neighbors_with_weights(W, i, listw$neighbours[[i]], listw$weights[[i]]) update_weights(W) } set.seed(1) B <- local_geary(W, data.frame(x)) all.equal(A, lisa_values(B)) } if (run) { set.seed(1) C <- localC_perm(x, listw, nsim = 499, conditional=TRUE, alternative="two.sided") cor(ifelse(lisa_pvalues(B) < 0.5, lisa_pvalues(B), 1-lisa_pvalues(B)), attr(C, "pseudo-p")[,6]) } # pseudo-p values probably wrongly folded https://github.com/GeoDaCenter/rgeoda/issues/28 tmap_ok <- FALSE if (require(tmap, quietly=TRUE)) tmap_ok <- TRUE if (run) { # doi: 10.1111/gean.12164 guerry_path <- system.file("extdata", "Guerry.shp", package = "rgeoda") g <- st_read(guerry_path)[, 7:12] cor(st_drop_geometry(g)) #(Tab. 1) lw <- nb2listw(poly2nb(g)) moran(g$Crm_prs, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Crm_prp, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Litercy, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Donatns, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Infants, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Suicids, lw, n=nrow(g), S0=Szero(lw))$I } if (run) { o <- prcomp(st_drop_geometry(g), scale.=TRUE) cor(st_drop_geometry(g), o$x[,1:2])^2 #(Tab. 2) } if (run) { g$PC1 <- o$x[, "PC1"] brks <- c(min(g$PC1), natural_breaks(k=6, g["PC1"]), max(g$PC1)) if (tmap_ok) { tmap4 <- packageVersion("tmap") >= "3.99" if (tmap4) { tm_shape(g) + tm_polygons(fill="PC1", fill.scale=tm_scale(values="brewer.rd_yl_gn", breaks=brks, midpoint=0), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("PC1", breaks=brks, midpoint=0) + tm_borders() # Fig. 1 } } else { pplot(g["PC1"], breaks=brks) } } if (run) { g$PC2 <- -1*o$x[, "PC2"] # eigenvalue sign arbitrary brks <- c(min(g$PC2), natural_breaks(k=6, g["PC2"]), max(g$PC2)) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="PC2", fill.scale=tm_scale(values="brewer.rd_yl_gn", breaks=brks, midpoint=0), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("PC2", breaks=brks, midpoint=0) + tm_borders() # Fig. 2 } } else { plot(g["PC2"], breaks=brks) } } if (run) { w <- queen_weights(g) lm_PC1 <- local_moran(w, g["PC1"], significance_cutoff=0.01, permutations=99999) g$lm_PC1 <- factor(lisa_clusters(lm_PC1), levels=0:4, labels=lisa_labels(lm_PC1)[1:5]) is.na(g$lm_PC1) <- g$lm_PC1 == "Not significant" g$lm_PC1 <- droplevels(g$lm_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lm_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lm_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 3 } } else { plot(g["lm_PC1"]) } } if (run) { set.seed(1) lm_PC1_spdep <- localmoran_perm(g$PC1, lw, nsim=9999) q <- attr(lm_PC1_spdep, "quadr")$pysal g$lm_PC1_spdep <- q is.na(g$lm_PC1_spdep) <- lm_PC1_spdep[,6] > 0.02 # note folded p-values g$lm_PC1_spdep <- droplevels(g$lm_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lm_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lm_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 3 } } else { plot(g["lm_PC1_spdep"]) } } if (run) { lg_PC1 <- local_g(w, g["PC1"], significance_cutoff=0.01, permutations=99999) g$lg_PC1 <- factor(lisa_clusters(lg_PC1), levels=0:2, labels=lisa_labels(lg_PC1)[0:3]) is.na(g$lg_PC1) <- g$lg_PC1 == "Not significant" g$lg_PC1 <- droplevels(g$lg_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lg_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lg_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 4 (wrong) } } else { plot(g["lg_PC1"]) } g$lg_PC1a <- cut(g$PC1, c(-Inf, mean(g$PC1), Inf), labels=c("Low", "High")) is.na(g$lg_PC1a) <- lisa_pvalues(lg_PC1) >= 0.01 g$lg_PC1a <- droplevels(g$lg_PC1a) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lg_PC1a", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lg_PC1a", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 4 (guess) } } else { plot(g["lg_PC1"]) } } if (run) { lc_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.01, permutations=99999) g$lc_PC1 <- factor(lisa_clusters(lc_PC1), levels=0:4, labels=lisa_labels(lc_PC1)[1:5]) is.na(g$lc_PC1) <- g$lc_PC1 == "Not significant" g$lc_PC1 <- droplevels(g$lc_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 5 } } else { plot(g["lc_PC1"]) } } if (run) { set.seed(1) system.time(lc_PC1_spdep <- localC_perm(g$PC1, lw, nsim=9999, alternative="two.sided")) } if (run) { if (require(parallel, quietly=TRUE)) { ncpus <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L # test with single core if (ncpus > 1L) ncpus <- 1L cores <- get.coresOption() set.coresOption(ncpus) system.time(lmc_PC1_spdep1 <- localC_perm(g$PC1, lw, nsim=9999, alternative="two.sided", iseed=1)) set.coresOption(cores) } } if (run) { g$lc_PC1_spdep <- attr(lc_PC1_spdep, "cluster") is.na(g$lc_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.01 g$lc_PC1_spdep <- droplevels(g$lc_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 5 } } else { plot(g["lc_PC1_spdep"]) } } if (run) { g$both_PC1 <- interaction(g$lc_PC1, g$lm_PC1) g$both_PC1 <- droplevels(g$both_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="both_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("both_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 6 } } else { plot(g["both_PC1"]) } } if (run) { lc005_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.005, permutations=99999) g$lc005_PC1 <- factor(lisa_clusters(lc005_PC1), levels=0:4, labels=lisa_labels(lc005_PC1)[1:5]) is.na(g$lc005_PC1) <- g$lc005_PC1 == "Not significant" g$lc005_PC1 <- droplevels(g$lc005_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc005_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 7 } } else { plot(g["lc005_PC1"]) } if (run) { g$lc005_PC1_spdep <- attr(lc_PC1_spdep, "cluster") is.na(g$lc005_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.005 g$lc005_PC1_spdep <- droplevels(g$lc005_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc005_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 7 } } else { plot(g["lc005_PC1_spdep"]) } } if (run) { lc001_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.001, permutations=99999) g$lc001_PC1 <- factor(lisa_clusters(lc001_PC1), levels=0:4, labels=lisa_labels(lc001_PC1)[1:5]) is.na(g$lc001_PC1) <- g$lc001_PC1 == "Not significant" g$lc001_PC1 <- droplevels(g$lc001_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc001_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 8 } } else { plot(g["lc001_PC1"]) } } if (run) { g$lc001_PC1_spdep <- attr(lc_PC1_spdep, "cluster") is.na(g$lc001_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.001 g$lc001_PC1_spdep <- droplevels(g$lc001_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc001_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 8 } } else { plot(g["lc001_PC1_spdep"]) } } } if (run) { lc_PC2 <- local_geary(w, g["PC2"], significance_cutoff=0.01, permutations=99999) g$lc_PC2 <- factor(lisa_clusters(lc_PC2), levels=0:4, labels=lisa_labels(lc_PC2)[1:5]) is.na(g$lc_PC2) <- g$lc_PC2 == "Not significant" g$lc_PC2 <- droplevels(g$lc_PC2) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc_PC2", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc_PC2", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 9 } } else { plot(g["lc_PC2"]) } } if (run) { lmc_PC <- local_multigeary(w, g[c("PC1","PC2")], significance_cutoff=0.00247, permutations=99999) g$lmc_PC <- factor(lisa_clusters(lmc_PC), levels=0:1, labels=lisa_labels(lmc_PC)[1:2]) is.na(g$lmc_PC) <- g$lmc_PC == "Not significant" g$lmc_PC <- droplevels(g$lmc_PC) table(interaction((p.adjust(lisa_pvalues(lmc_PC), "fdr") < 0.01), g$lmc_PC)) } if (run) { if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_PC", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_PC", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 10 } } else { plot(g["lmc_PC"]) } } if (run) { set.seed(1) lmc_PC_spdep <- localC_perm(g[c("PC1","PC2")], lw, nsim=9999, alternative="two.sided") all.equal(lisa_values(lmc_PC), c(lmc_PC_spdep)) } if (run) { cor(attr(lmc_PC_spdep, "pseudo-p")[,6], lisa_pvalues(lmc_PC)) } if (run) { g$lmc_PC_spdep <- attr(lmc_PC_spdep, "cluster") is.na(g$lmc_PC_spdep) <- p.adjust(attr(lmc_PC_spdep, "pseudo-p")[,6], "fdr") > 0.01 g$lmc_PC_spdep <- droplevels(g$lmc_PC_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_PC_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_PC_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 10 } } else { plot(g["lmc_PC_spdep"]) } } if (run) { lmc_vars <- local_multigeary(w, st_drop_geometry(g)[, 1:6], significance_cutoff=0.00247, permutations=99999) g$lmc_vars <- factor(lisa_clusters(lmc_vars), levels=0:1, labels=lisa_labels(lmc_vars)[1:2]) is.na(g$lmc_vars) <- g$lmc_vars == "Not significant" g$lmc_vars <- droplevels(g$lmc_vars) table(interaction((p.adjust(lisa_pvalues(lmc_vars), "fdr") < 0.01), g$lmc_vars)) } if (run) { if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_vars", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_vars", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 11 } } else { plot(g["lmc_vars"]) } } if (run) { set.seed(1) system.time(lmc_vars_spdep <- localC_perm(st_drop_geometry(g)[, 1:6], lw, nsim=9999, alternative="two.sided")) } if (run) { all.equal(lisa_values(lmc_vars), c(lmc_vars_spdep)) } if (run) { cor(attr(lmc_vars_spdep, "pseudo-p")[,6], lisa_pvalues(lmc_vars)) } if (run) { if (require(parallel, quietly=TRUE)) { ncpus <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L # test with single core if (ncpus > 1L) ncpus <- 1L cores <- get.coresOption() set.coresOption(ncpus) system.time(lmc_vars_spdep1 <- localC_perm(st_drop_geometry(g)[, 1:6], lw, nsim=9999, alternative="two.sided", iseed=1)) set.coresOption(cores) } } if (run) { all.equal(lisa_values(lmc_vars), c(lmc_vars_spdep1)) } if (run) { cor(attr(lmc_vars_spdep1, "pseudo-p")[,6], lisa_pvalues(lmc_vars)) } if (run) { g$lmc_vars_spdep <- attr(lmc_vars_spdep1, "cluster") is.na(g$lmc_vars_spdep) <- p.adjust(attr(lmc_vars_spdep1, "pseudo-p")[,6], "fdr") > 0.01 g$lmc_vars_spdep <- droplevels(g$lmc_vars_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_vars_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_vars_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 11 } } else { plot(g["lmc_vars_spdep"]) } } ## Not run: library(reticulate) use_python("/usr/bin/python", required = TRUE) gp <- import("geopandas") ps <- import("libpysal") W <- listw2mat(listw) w <- ps$weights$full2W(W, rownames(W)) w$transform <- "R" esda <- import("esda") lM <- esda$Moran_Local(x, w) all.equal(unname(localmoran(x, listw, mlvar=FALSE)[,1]), c(lM$Is)) # confirm x and w the same lC <- esda$Geary_Local(connectivity=w)$fit(scale(x)) # np$std missing ddof=1 n <- length(x) D0 <- spdep:::geary.intern((x - mean(x)) / sqrt(var(x)*(n-1)/n), listw, n=n) # lC components probably wrongly ordered https://github.com/pysal/esda/issues/192 o <- match(round(D0, 6), round(lC$localG, 6)) all.equal(c(lC$localG)[o], D0) # simulation order not retained lC$p_sim[o] attr(C, "pseudo-p")[,6] ## End(Not run)
orig <- spData::africa.rook.nb listw <- nb2listw(orig) x <- spData::afcon$totcon (A <- localC(x, listw)) listw1 <- nb2listw(droplinks(sym.attr.nb(orig), 3, sym=TRUE), zero.policy=TRUE) (A1 <- localC(x, listw1, zero.policy=FALSE)) (A2 <- localC(x, listw1, zero.policy=TRUE)) run <- FALSE if (require(rgeoda, quietly=TRUE)) run <- TRUE if (run) { W <- create_weights(as.numeric(length(x))) for (i in 1:length(listw$neighbours)) { set_neighbors_with_weights(W, i, listw$neighbours[[i]], listw$weights[[i]]) update_weights(W) } set.seed(1) B <- local_geary(W, data.frame(x)) all.equal(A, lisa_values(B)) } if (run) { set.seed(1) C <- localC_perm(x, listw, nsim = 499, conditional=TRUE, alternative="two.sided") cor(ifelse(lisa_pvalues(B) < 0.5, lisa_pvalues(B), 1-lisa_pvalues(B)), attr(C, "pseudo-p")[,6]) } # pseudo-p values probably wrongly folded https://github.com/GeoDaCenter/rgeoda/issues/28 tmap_ok <- FALSE if (require(tmap, quietly=TRUE)) tmap_ok <- TRUE if (run) { # doi: 10.1111/gean.12164 guerry_path <- system.file("extdata", "Guerry.shp", package = "rgeoda") g <- st_read(guerry_path)[, 7:12] cor(st_drop_geometry(g)) #(Tab. 1) lw <- nb2listw(poly2nb(g)) moran(g$Crm_prs, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Crm_prp, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Litercy, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Donatns, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Infants, lw, n=nrow(g), S0=Szero(lw))$I moran(g$Suicids, lw, n=nrow(g), S0=Szero(lw))$I } if (run) { o <- prcomp(st_drop_geometry(g), scale.=TRUE) cor(st_drop_geometry(g), o$x[,1:2])^2 #(Tab. 2) } if (run) { g$PC1 <- o$x[, "PC1"] brks <- c(min(g$PC1), natural_breaks(k=6, g["PC1"]), max(g$PC1)) if (tmap_ok) { tmap4 <- packageVersion("tmap") >= "3.99" if (tmap4) { tm_shape(g) + tm_polygons(fill="PC1", fill.scale=tm_scale(values="brewer.rd_yl_gn", breaks=brks, midpoint=0), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("PC1", breaks=brks, midpoint=0) + tm_borders() # Fig. 1 } } else { pplot(g["PC1"], breaks=brks) } } if (run) { g$PC2 <- -1*o$x[, "PC2"] # eigenvalue sign arbitrary brks <- c(min(g$PC2), natural_breaks(k=6, g["PC2"]), max(g$PC2)) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="PC2", fill.scale=tm_scale(values="brewer.rd_yl_gn", breaks=brks, midpoint=0), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("PC2", breaks=brks, midpoint=0) + tm_borders() # Fig. 2 } } else { plot(g["PC2"], breaks=brks) } } if (run) { w <- queen_weights(g) lm_PC1 <- local_moran(w, g["PC1"], significance_cutoff=0.01, permutations=99999) g$lm_PC1 <- factor(lisa_clusters(lm_PC1), levels=0:4, labels=lisa_labels(lm_PC1)[1:5]) is.na(g$lm_PC1) <- g$lm_PC1 == "Not significant" g$lm_PC1 <- droplevels(g$lm_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lm_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lm_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 3 } } else { plot(g["lm_PC1"]) } } if (run) { set.seed(1) lm_PC1_spdep <- localmoran_perm(g$PC1, lw, nsim=9999) q <- attr(lm_PC1_spdep, "quadr")$pysal g$lm_PC1_spdep <- q is.na(g$lm_PC1_spdep) <- lm_PC1_spdep[,6] > 0.02 # note folded p-values g$lm_PC1_spdep <- droplevels(g$lm_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lm_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lm_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 3 } } else { plot(g["lm_PC1_spdep"]) } } if (run) { lg_PC1 <- local_g(w, g["PC1"], significance_cutoff=0.01, permutations=99999) g$lg_PC1 <- factor(lisa_clusters(lg_PC1), levels=0:2, labels=lisa_labels(lg_PC1)[0:3]) is.na(g$lg_PC1) <- g$lg_PC1 == "Not significant" g$lg_PC1 <- droplevels(g$lg_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lg_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lg_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 4 (wrong) } } else { plot(g["lg_PC1"]) } g$lg_PC1a <- cut(g$PC1, c(-Inf, mean(g$PC1), Inf), labels=c("Low", "High")) is.na(g$lg_PC1a) <- lisa_pvalues(lg_PC1) >= 0.01 g$lg_PC1a <- droplevels(g$lg_PC1a) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lg_PC1a", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lg_PC1a", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 4 (guess) } } else { plot(g["lg_PC1"]) } } if (run) { lc_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.01, permutations=99999) g$lc_PC1 <- factor(lisa_clusters(lc_PC1), levels=0:4, labels=lisa_labels(lc_PC1)[1:5]) is.na(g$lc_PC1) <- g$lc_PC1 == "Not significant" g$lc_PC1 <- droplevels(g$lc_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 5 } } else { plot(g["lc_PC1"]) } } if (run) { set.seed(1) system.time(lc_PC1_spdep <- localC_perm(g$PC1, lw, nsim=9999, alternative="two.sided")) } if (run) { if (require(parallel, quietly=TRUE)) { ncpus <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L # test with single core if (ncpus > 1L) ncpus <- 1L cores <- get.coresOption() set.coresOption(ncpus) system.time(lmc_PC1_spdep1 <- localC_perm(g$PC1, lw, nsim=9999, alternative="two.sided", iseed=1)) set.coresOption(cores) } } if (run) { g$lc_PC1_spdep <- attr(lc_PC1_spdep, "cluster") is.na(g$lc_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.01 g$lc_PC1_spdep <- droplevels(g$lc_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 5 } } else { plot(g["lc_PC1_spdep"]) } } if (run) { g$both_PC1 <- interaction(g$lc_PC1, g$lm_PC1) g$both_PC1 <- droplevels(g$both_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="both_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("both_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 6 } } else { plot(g["both_PC1"]) } } if (run) { lc005_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.005, permutations=99999) g$lc005_PC1 <- factor(lisa_clusters(lc005_PC1), levels=0:4, labels=lisa_labels(lc005_PC1)[1:5]) is.na(g$lc005_PC1) <- g$lc005_PC1 == "Not significant" g$lc005_PC1 <- droplevels(g$lc005_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc005_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 7 } } else { plot(g["lc005_PC1"]) } if (run) { g$lc005_PC1_spdep <- attr(lc_PC1_spdep, "cluster") is.na(g$lc005_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.005 g$lc005_PC1_spdep <- droplevels(g$lc005_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc005_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 7 } } else { plot(g["lc005_PC1_spdep"]) } } if (run) { lc001_PC1 <- local_geary(w, g["PC1"], significance_cutoff=0.001, permutations=99999) g$lc001_PC1 <- factor(lisa_clusters(lc001_PC1), levels=0:4, labels=lisa_labels(lc001_PC1)[1:5]) is.na(g$lc001_PC1) <- g$lc001_PC1 == "Not significant" g$lc001_PC1 <- droplevels(g$lc001_PC1) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc001_PC1", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 8 } } else { plot(g["lc001_PC1"]) } } if (run) { g$lc001_PC1_spdep <- attr(lc_PC1_spdep, "cluster") is.na(g$lc001_PC1_spdep) <- attr(lc_PC1_spdep, "pseudo-p")[,6] > 0.001 g$lc001_PC1_spdep <- droplevels(g$lc001_PC1_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc005_PC1_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc001_PC1_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 8 } } else { plot(g["lc001_PC1_spdep"]) } } } if (run) { lc_PC2 <- local_geary(w, g["PC2"], significance_cutoff=0.01, permutations=99999) g$lc_PC2 <- factor(lisa_clusters(lc_PC2), levels=0:4, labels=lisa_labels(lc_PC2)[1:5]) is.na(g$lc_PC2) <- g$lc_PC2 == "Not significant" g$lc_PC2 <- droplevels(g$lc_PC2) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lc_PC2", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lc_PC2", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 9 } } else { plot(g["lc_PC2"]) } } if (run) { lmc_PC <- local_multigeary(w, g[c("PC1","PC2")], significance_cutoff=0.00247, permutations=99999) g$lmc_PC <- factor(lisa_clusters(lmc_PC), levels=0:1, labels=lisa_labels(lmc_PC)[1:2]) is.na(g$lmc_PC) <- g$lmc_PC == "Not significant" g$lmc_PC <- droplevels(g$lmc_PC) table(interaction((p.adjust(lisa_pvalues(lmc_PC), "fdr") < 0.01), g$lmc_PC)) } if (run) { if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_PC", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_PC", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 10 } } else { plot(g["lmc_PC"]) } } if (run) { set.seed(1) lmc_PC_spdep <- localC_perm(g[c("PC1","PC2")], lw, nsim=9999, alternative="two.sided") all.equal(lisa_values(lmc_PC), c(lmc_PC_spdep)) } if (run) { cor(attr(lmc_PC_spdep, "pseudo-p")[,6], lisa_pvalues(lmc_PC)) } if (run) { g$lmc_PC_spdep <- attr(lmc_PC_spdep, "cluster") is.na(g$lmc_PC_spdep) <- p.adjust(attr(lmc_PC_spdep, "pseudo-p")[,6], "fdr") > 0.01 g$lmc_PC_spdep <- droplevels(g$lmc_PC_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_PC_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_PC_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 10 } } else { plot(g["lmc_PC_spdep"]) } } if (run) { lmc_vars <- local_multigeary(w, st_drop_geometry(g)[, 1:6], significance_cutoff=0.00247, permutations=99999) g$lmc_vars <- factor(lisa_clusters(lmc_vars), levels=0:1, labels=lisa_labels(lmc_vars)[1:2]) is.na(g$lmc_vars) <- g$lmc_vars == "Not significant" g$lmc_vars <- droplevels(g$lmc_vars) table(interaction((p.adjust(lisa_pvalues(lmc_vars), "fdr") < 0.01), g$lmc_vars)) } if (run) { if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_vars", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_vars", textNA="Insignificant", colorNA="gray95") + tm_borders() # Fig. 11 } } else { plot(g["lmc_vars"]) } } if (run) { set.seed(1) system.time(lmc_vars_spdep <- localC_perm(st_drop_geometry(g)[, 1:6], lw, nsim=9999, alternative="two.sided")) } if (run) { all.equal(lisa_values(lmc_vars), c(lmc_vars_spdep)) } if (run) { cor(attr(lmc_vars_spdep, "pseudo-p")[,6], lisa_pvalues(lmc_vars)) } if (run) { if (require(parallel, quietly=TRUE)) { ncpus <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L # test with single core if (ncpus > 1L) ncpus <- 1L cores <- get.coresOption() set.coresOption(ncpus) system.time(lmc_vars_spdep1 <- localC_perm(st_drop_geometry(g)[, 1:6], lw, nsim=9999, alternative="two.sided", iseed=1)) set.coresOption(cores) } } if (run) { all.equal(lisa_values(lmc_vars), c(lmc_vars_spdep1)) } if (run) { cor(attr(lmc_vars_spdep1, "pseudo-p")[,6], lisa_pvalues(lmc_vars)) } if (run) { g$lmc_vars_spdep <- attr(lmc_vars_spdep1, "cluster") is.na(g$lmc_vars_spdep) <- p.adjust(attr(lmc_vars_spdep1, "pseudo-p")[,6], "fdr") > 0.01 g$lmc_vars_spdep <- droplevels(g$lmc_vars_spdep) if (tmap_ok) { if (tmap4) { tm_shape(g) + tm_polygons(fill="lmc_vars_spdep", fill.scale=tm_scale(values="brewer.set3", value.na="gray95", label.na="Insignificant"), fill.legend=tm_legend(position=tm_pos_in("left", "bottom"), frame=FALSE, item.r=0)) } else { tm_shape(g) + tm_fill("lmc_vars_spdep", textNA="Insignificant", colorNA="gray95") + tm_borders() # rep. Fig. 11 } } else { plot(g["lmc_vars_spdep"]) } } ## Not run: library(reticulate) use_python("/usr/bin/python", required = TRUE) gp <- import("geopandas") ps <- import("libpysal") W <- listw2mat(listw) w <- ps$weights$full2W(W, rownames(W)) w$transform <- "R" esda <- import("esda") lM <- esda$Moran_Local(x, w) all.equal(unname(localmoran(x, listw, mlvar=FALSE)[,1]), c(lM$Is)) # confirm x and w the same lC <- esda$Geary_Local(connectivity=w)$fit(scale(x)) # np$std missing ddof=1 n <- length(x) D0 <- spdep:::geary.intern((x - mean(x)) / sqrt(var(x)*(n-1)/n), listw, n=n) # lC components probably wrongly ordered https://github.com/pysal/esda/issues/192 o <- match(round(D0, 6), round(lC$localG, 6)) all.equal(c(lC$localG)[o], D0) # simulation order not retained lC$p_sim[o] attr(C, "pseudo-p")[,6] ## End(Not run)
The local spatial statistic G is calculated for each zone based on the spatial weights object used. The value returned is a Z-value, and may be used as a diagnostic tool. High positive values indicate the posibility of a local cluster of high values of the variable being analysed, very low relative values a similar cluster of low values. For inference, a Bonferroni-type test is suggested in the references, where tables of critical values may be found (see also details below).
localG(x, listw, zero.policy=attr(listw, "zero.policy"), spChk=NULL, GeoDa=FALSE, alternative = "two.sided", return_internals=TRUE) localG_perm(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy"), spChk=NULL, alternative = "two.sided", iseed=NULL, fix_i_in_Gstar_permutations=TRUE, no_repeat_in_row=FALSE)
localG(x, listw, zero.policy=attr(listw, "zero.policy"), spChk=NULL, GeoDa=FALSE, alternative = "two.sided", return_internals=TRUE) localG_perm(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy"), spChk=NULL, alternative = "two.sided", iseed=NULL, fix_i_in_Gstar_permutations=TRUE, no_repeat_in_row=FALSE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
zero.policy |
default NULL, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
GeoDa |
default FALSE, if TRUE, drop x values for no-neighbour and self-neighbour only observations from all summations |
nsim |
default 499, number of conditonal permutation simulations |
alternative |
a character string specifying the alternative hypothesis, must be one of |
return_internals |
default |
iseed |
default NULL, used to set the seed; the output will only be reproducible if the count of CPU cores across which computation is distributed is the same |
fix_i_in_Gstar_permutations |
default |
no_repeat_in_row |
default |
If the neighbours member of listw has a "self.included" attribute set
to TRUE, the Gstar variant, including the self-weight ,
is calculated and returned. The returned vector will have a "gstari"
attribute set to TRUE. Self-weights must be included by using the
include.self
function before converting
the neighbour list to a spatial weights list with nb2listw
as
shown below in the example.
The critical values of the statistic under assumptions given in the references for the 95th percentile are for n=1: 1.645, n=50: 3.083, n=100: 3.289, n=1000: 3.886.
A vector of G or Gstar standard deviate values, with attributes "gstari" set to TRUE or FALSE, "call" set to the function call, and class "localG". For conditional permutation, the returned value is the same as for localG()
, and the simulated standard deviate is returned as column "StdDev.Gi"
in attr(., "internals")
.
Conditional permutations added for comparative purposes; permutations are over the whole data vector omitting the observation itself, and from 1.2-8 fixing the observation itself as its own neighbour for local G-star.
Roger Bivand [email protected]
Ord, J. K. and Getis, A. 1995 Local spatial autocorrelation statistics: distributional issues and an application. Geographical Analysis, 27, 286–306; Getis, A. and Ord, J. K. 1996 Local spatial statistics: an overview. In P. Longley and M. Batty (eds) Spatial analysis: modelling in a GIS environment (Cambridge: Geoinformation International), 261–277; Bivand RS, Wong DWS 2018 Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748 doi:10.1007/s11749-018-0599-x
data(getisord, package="spData") # spData 0.3.2 changes x, y, xyz object names to go_x, go_y, go_xyz to # avoid putting these objects into the global environment via lazy loading if (exists("go_xyz") && packageVersion("spData") >= "0.3.2") { xyz <- go_xyz x <- go_x y <- go_y } xycoords <- cbind(xyz$x, xyz$y) nb30 <- dnearneigh(xycoords, 0, 30) G30 <- localG(xyz$val, nb2listw(nb30, style="B")) G30[length(xyz$val)-136] set.seed(1) G30_sim <- localG_perm(xyz$val, nb2listw(nb30, style="B")) G30_sim[length(xyz$val)-136] nb60 <- dnearneigh(xycoords, 0, 60) G60 <- localG(xyz$val, nb2listw(nb60, style="B")) G60[length(xyz$val)-136] nb90 <- dnearneigh(xycoords, 0, 90) G90 <- localG(xyz$val, nb2listw(nb90, style="B")) G90[length(xyz$val)-136] nb120 <- dnearneigh(xycoords, 0, 120) G120 <- localG(xyz$val, nb2listw(nb120, style="B")) G120[length(xyz$val)-136] nb150 <- dnearneigh(xycoords, 0, 150) G150 <- localG(xyz$val, nb2listw(nb150, style="B")) G150[length(xyz$val)-136] brks <- seq(-5,5,1) cm.col <- cm.colors(length(brks)-1) image(x, y, t(matrix(G30, nrow=16, ncol=16, byrow=TRUE)), breaks=brks, col=cm.col, asp=1) text(xyz$x, xyz$y, round(G30, digits=1), cex=0.7) polygon(c(195,225,225,195), c(195,195,225,225), lwd=2) title(main=expression(paste("Values of the ", G[i], " statistic"))) G30s <- localG(xyz$val, nb2listw(include.self(nb30), style="B")) cat("value according to Getis and Ord's eq. 14.2, p. 263 (1996)\n") G30s[length(xyz$val)-136] cat(paste("value given by Getis and Ord (1996), p. 267", "(division by n-1 rather than n \n in variance)\n")) G30s[length(xyz$val)-136] * (sqrt(sum(scale(xyz$val, scale=FALSE)^2)/length(xyz$val)) / sqrt(var(xyz$val))) image(x, y, t(matrix(G30s, nrow=16, ncol=16, byrow=TRUE)), breaks=brks, col=cm.col, asp=1) text(xyz$x, xyz$y, round(G30s, digits=1), cex=0.7) polygon(c(195,225,225,195), c(195,195,225,225), lwd=2) title(main=expression(paste("Values of the ", G[i]^"*", " statistic")))
data(getisord, package="spData") # spData 0.3.2 changes x, y, xyz object names to go_x, go_y, go_xyz to # avoid putting these objects into the global environment via lazy loading if (exists("go_xyz") && packageVersion("spData") >= "0.3.2") { xyz <- go_xyz x <- go_x y <- go_y } xycoords <- cbind(xyz$x, xyz$y) nb30 <- dnearneigh(xycoords, 0, 30) G30 <- localG(xyz$val, nb2listw(nb30, style="B")) G30[length(xyz$val)-136] set.seed(1) G30_sim <- localG_perm(xyz$val, nb2listw(nb30, style="B")) G30_sim[length(xyz$val)-136] nb60 <- dnearneigh(xycoords, 0, 60) G60 <- localG(xyz$val, nb2listw(nb60, style="B")) G60[length(xyz$val)-136] nb90 <- dnearneigh(xycoords, 0, 90) G90 <- localG(xyz$val, nb2listw(nb90, style="B")) G90[length(xyz$val)-136] nb120 <- dnearneigh(xycoords, 0, 120) G120 <- localG(xyz$val, nb2listw(nb120, style="B")) G120[length(xyz$val)-136] nb150 <- dnearneigh(xycoords, 0, 150) G150 <- localG(xyz$val, nb2listw(nb150, style="B")) G150[length(xyz$val)-136] brks <- seq(-5,5,1) cm.col <- cm.colors(length(brks)-1) image(x, y, t(matrix(G30, nrow=16, ncol=16, byrow=TRUE)), breaks=brks, col=cm.col, asp=1) text(xyz$x, xyz$y, round(G30, digits=1), cex=0.7) polygon(c(195,225,225,195), c(195,195,225,225), lwd=2) title(main=expression(paste("Values of the ", G[i], " statistic"))) G30s <- localG(xyz$val, nb2listw(include.self(nb30), style="B")) cat("value according to Getis and Ord's eq. 14.2, p. 263 (1996)\n") G30s[length(xyz$val)-136] cat(paste("value given by Getis and Ord (1996), p. 267", "(division by n-1 rather than n \n in variance)\n")) G30s[length(xyz$val)-136] * (sqrt(sum(scale(xyz$val, scale=FALSE)^2)/length(xyz$val)) / sqrt(var(xyz$val))) image(x, y, t(matrix(G30s, nrow=16, ncol=16, byrow=TRUE)), breaks=brks, col=cm.col, asp=1) text(xyz$x, xyz$y, round(G30s, digits=1), cex=0.7) polygon(c(195,225,225,195), c(195,195,225,225), lwd=2) title(main=expression(paste("Values of the ", G[i]^"*", " statistic")))
The function implements the test statistic for local hotspots on specific pairwise evaluated distance bands, as proposed by Westerholt et al. (2015). Like the hotspot estimator
, the
statistic is given as z-scores that can be evaluated accordingly. The idea of the method is to identify hotspots in datasets that comprise several, difficult-to-separate processes operating at different scales. This is often the case in complex user-generated datasets such as those from Twitter feeds. For example, a football match could be reflected in tweets from pubs, homes, and the stadium vicinity. These exemplified phenomena represent different processes that may be detected at different geometric scales. The
method enables this identification by specifying a geometric scale band and strictly calculating all statistical quantities such as mean and variance solely from respective relevant observations that interact on the range of the adjusted scale band. In addition, in each neighbourhood not only the relationships to the respective central unit, but all scale-relevant relationships are considered. In this way, hotspots can be detected on specific scale ranges independently of other scales. The statistic is given as:
with
localGS(x, listw, dmin, dmax, attr, longlat = NULL)
localGS(x, listw, dmin, dmax, attr, longlat = NULL)
x |
a |
listw |
a |
dmin |
a lower distance bound (greater than or equal) |
dmax |
an upper distance bound (less than or equal) |
attr |
the name of the attribute of interest |
longlat |
default NULL; TRUE if point coordinates are longitude-latitude decimal degrees, in which case distances are measured in kilometres; if x is a SpatialPoints object, the value is taken from the object itself, and overrides this argument if not NULL; distances are measured in map units if FALSE or NULL |
Only pairs of observations with a shared distance (in map units) on the interval [dmin
, dmax
] that are within a maximum radius of dmax
around a corresponding output observation are considered. Thereby, also the mean values and variance terms estimated within the measure are adjusted to the scale range under consideration. For application examples of the method see Westerholt et al. (2015) (applied to tweets) and Sonea & Westerholt (2021) (applied in an access to banking scenario).
A vector of values that are given as z-scores.
René Westerholt [email protected]
Westerholt, R., Resch, B. & Zipf, A. 2015. A local scale-sensitive indicator of spatial autocorrelation for assessing high-and low-value clusters in multiscale datasets. International Journal of Geographical Information Science, 29(5), 868–887, doi:10.1080/13658816.2014.1002499.
Sonea, A. and Westerholt, R. (2021): Geographic and temporal access to basic banking services in Wales. Applied Spatial Analysis and Policy, 14 (4), 879–905, doi:10.1007/s12061-021-09386-3.
boston.tr <- sf::st_read(system.file("shapes/boston_tracts.gpkg", package="spData")[1]) boston.tr_utm <- st_transform(boston.tr, 32619) #26786 boston_listw1 <- nb2listwdist(dnearneigh(st_centroid(boston.tr_utm), 1, 2000), boston.tr_utm, type = "dpd", alpha = 2, zero.policy = TRUE, dmax = 9500) boston_listw2 <- nb2listwdist(dnearneigh(st_centroid(boston.tr_utm), 5000, 9500), boston.tr_utm, type = "dpd", alpha = 2, zero.policy = TRUE, dmax = 9500) boston_RM_gsi_1 <- localGS(boston.tr_utm, boston_listw1, 1, 2000, "RM", FALSE) boston_RM_gsi_2 <- localGS(boston.tr_utm, boston_listw2, 2000, 9500, "RM", FALSE)
boston.tr <- sf::st_read(system.file("shapes/boston_tracts.gpkg", package="spData")[1]) boston.tr_utm <- st_transform(boston.tr, 32619) #26786 boston_listw1 <- nb2listwdist(dnearneigh(st_centroid(boston.tr_utm), 1, 2000), boston.tr_utm, type = "dpd", alpha = 2, zero.policy = TRUE, dmax = 9500) boston_listw2 <- nb2listwdist(dnearneigh(st_centroid(boston.tr_utm), 5000, 9500), boston.tr_utm, type = "dpd", alpha = 2, zero.policy = TRUE, dmax = 9500) boston_RM_gsi_1 <- localGS(boston.tr_utm, boston_listw1, 1, 2000, "RM", FALSE) boston_RM_gsi_2 <- localGS(boston.tr_utm, boston_listw2, 2000, 9500, "RM", FALSE)
The local spatial statistic Moran's I is calculated for each zone based on the spatial weights object used. The values returned include a Z-value, and may be used as a diagnostic tool. The statistic is:
, and its expectation and variance were given in Anselin (1995), but those from Sokal et al. (1998) are implemented here.
localmoran(x, listw, zero.policy=attr(listw, "zero.policy"), na.action=na.fail, conditional=TRUE, alternative = "two.sided", mlvar=TRUE, spChk=NULL, adjust.x=FALSE) localmoran_perm(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy"), na.action=na.fail, alternative = "two.sided", mlvar=TRUE, spChk=NULL, adjust.x=FALSE, sample_Ei=TRUE, iseed=NULL, no_repeat_in_row=FALSE)
localmoran(x, listw, zero.policy=attr(listw, "zero.policy"), na.action=na.fail, conditional=TRUE, alternative = "two.sided", mlvar=TRUE, spChk=NULL, adjust.x=FALSE) localmoran_perm(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy"), na.action=na.fail, alternative = "two.sided", mlvar=TRUE, spChk=NULL, adjust.x=FALSE, sample_Ei=TRUE, iseed=NULL, no_repeat_in_row=FALSE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
zero.policy |
default default |
na.action |
a function (default |
conditional |
default TRUE: expectation and variance are calculated using the conditional randomization null (Sokal 1998 Eqs. A7 & A8). Elaboration of these changes available in Sauer et al. (2021). If FALSE: expectation and variance are calculated using the total randomization null (Sokal 1998 Eqs. A3 & A4). |
alternative |
a character string specifying the alternative hypothesis, must be one of greater, less or two.sided (default). |
mlvar |
default TRUE: values of local Moran's I are reported using the variance of the variable of interest (sum of squared deviances over n), but can be reported as the sample variance, dividing by (n-1) instead; both are used in other implementations. |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
adjust.x |
default FALSE, if TRUE, x values of observations with no neighbours are omitted in the mean of x |
nsim |
default 499, number of conditonal permutation simulations |
sample_Ei |
default TRUE; if conditional permutation, use the sample $E_i$ values, or the analytical values, leaving only variances calculated by simulation. |
iseed |
default NULL, used to set the seed; the output will only be reproducible if the count of CPU cores across which computation is distributed is the same |
no_repeat_in_row |
default |
The values of local Moran's I are divided by the variance (or sample variance) of the variable of interest to accord with Table 1, p. 103, and formula (12), p. 99, in Anselin (1995), rather than his formula (7), p. 98. The variance of the local Moran statistic is taken from Sokal et al. (1998) p. 334, equations 4 & 5 or equations 7 & 8 located depending on user specification. By default, the implementation divides by n, not (n-1) in calculating the variance and higher moments. Conditional code contributed by Jeff Sauer and Levi Wolf.
Ii |
local moran statistic |
E.Ii |
expectation of local moran statistic; for |
Var.Ii |
variance of local moran statistic; for |
Z.Ii |
standard deviate of local moran statistic; for |
Pr() |
p-value of local moran statistic using |
Pr() Sim |
For |
Pr(folded) Sim |
the simulation folded [0, 0.5] range ranked p-value (based on https://github.com/pysal/esda/blob/4a63e0b5df1e754b17b5f1205b8cadcbecc5e061/esda/crand.py#L211-L213) |
Skewness |
For |
Kurtosis |
For |
In addition, an attribute data frame "quadr"
with mean and median quadrant columns, and a column splitting on the demeaned variable and lagged demeaned variable at zero.
Conditional permutations added for comparative purposes; permutations are over the whole data vector omitting the observation itself. For p-value adjustment, use p.adjust()
or p.adjustSP()
on the output vector.
Roger Bivand [email protected]
Anselin, L. 1995. Local indicators of spatial association, Geographical Analysis, 27, 93–115; Getis, A. and Ord, J. K. 1996 Local spatial statistics: an overview. In P. Longley and M. Batty (eds) Spatial analysis: modelling in a GIS environment (Cambridge: Geoinformation International), 261–277; Sokal, R. R, Oden, N. L. and Thomson, B. A. 1998. Local Spatial Autocorrelation in a Biological Model. Geographical Analysis, 30. 331–354; Bivand RS, Wong DWS 2018 Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748 doi:10.1007/s11749-018-0599-x; Sauer, J., Oshan, T. M., Rey, S., & Wolf, L. J. 2021. The Importance of Null Hypotheses: Understanding Differences in Local Moran’s under Heteroskedasticity. Geographical Analysis. doi:10.1111/gean.12304
Bivand, R. (2022), R Packages for Analyzing Spatial Data: A Comparative Case Study with Areal Data. Geographical Analysis, 54(3), 488-518. doi:10.1111/gean.12319
data(afcon, package="spData") oid <- order(afcon$id) resI <- localmoran(afcon$totcon, nb2listw(paper.nb)) printCoefmat(data.frame(resI[oid,], row.names=afcon$name[oid]), check.names=FALSE) hist(resI[,5]) mean(resI[,1]) sum(resI[,1])/Szero(nb2listw(paper.nb)) moran.test(afcon$totcon, nb2listw(paper.nb)) # note equality for mean() only when the sum of weights equals # the number of observations (thanks to Juergen Symanzik) resI <- localmoran(afcon$totcon, nb2listw(paper.nb)) printCoefmat(data.frame(resI[oid,], row.names=afcon$name[oid]), check.names=FALSE) hist(p.adjust(resI[,5], method="bonferroni")) totcon <-afcon$totcon is.na(totcon) <- sample(1:length(totcon), 5) totcon resI.na <- localmoran(totcon, nb2listw(paper.nb), na.action=na.exclude, zero.policy=TRUE) if (class(attr(resI.na, "na.action")) == "exclude") { print(data.frame(resI.na[oid,], row.names=afcon$name[oid]), digits=2) } else print(resI.na, digits=2) resG <- localG(afcon$totcon, nb2listw(include.self(paper.nb))) print(data.frame(resG[oid], row.names=afcon$name[oid]), digits=2) set.seed(1) resI_p <- localmoran_perm(afcon$totcon, nb2listw(paper.nb)) printCoefmat(data.frame(resI_p[oid,], row.names=afcon$name[oid]), check.names=FALSE)
data(afcon, package="spData") oid <- order(afcon$id) resI <- localmoran(afcon$totcon, nb2listw(paper.nb)) printCoefmat(data.frame(resI[oid,], row.names=afcon$name[oid]), check.names=FALSE) hist(resI[,5]) mean(resI[,1]) sum(resI[,1])/Szero(nb2listw(paper.nb)) moran.test(afcon$totcon, nb2listw(paper.nb)) # note equality for mean() only when the sum of weights equals # the number of observations (thanks to Juergen Symanzik) resI <- localmoran(afcon$totcon, nb2listw(paper.nb)) printCoefmat(data.frame(resI[oid,], row.names=afcon$name[oid]), check.names=FALSE) hist(p.adjust(resI[,5], method="bonferroni")) totcon <-afcon$totcon is.na(totcon) <- sample(1:length(totcon), 5) totcon resI.na <- localmoran(totcon, nb2listw(paper.nb), na.action=na.exclude, zero.policy=TRUE) if (class(attr(resI.na, "na.action")) == "exclude") { print(data.frame(resI.na[oid,], row.names=afcon$name[oid]), digits=2) } else print(resI.na, digits=2) resG <- localG(afcon$totcon, nb2listw(include.self(paper.nb))) print(data.frame(resG[oid], row.names=afcon$name[oid]), digits=2) set.seed(1) resI_p <- localmoran_perm(afcon$totcon, nb2listw(paper.nb)) printCoefmat(data.frame(resI_p[oid,], row.names=afcon$name[oid]), check.names=FALSE)
Given two continuous numeric variables, calculate the bivariate Local Moran's I.
localmoran_bv(x, y, listw, nsim = 199, scale = TRUE, alternative="two.sided", iseed=1L, no_repeat_in_row=FALSE, zero.policy=attr(listw, "zero.policy"))
localmoran_bv(x, y, listw, nsim = 199, scale = TRUE, alternative="two.sided", iseed=1L, no_repeat_in_row=FALSE, zero.policy=attr(listw, "zero.policy"))
x |
a numeric vector of same length as |
y |
a numeric vector of same length as |
listw |
a listw object for example as created by |
nsim |
the number of simulations to run. |
scale |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "two.sided", or "less". |
iseed |
default NULL, used to set the seed; the output will only be reproducible if the count of CPU cores across which computation is distributed is the same. |
no_repeat_in_row |
default |
zero.policy |
default default |
The Bivariate Local Moran, like its global counterpart, evaluates the value
of x at observation i with its spatial neighbors' value of y. The value of is xi * Wyi. Or, in simpler words, the local bivariate Moran is the result of multiplying x by the spatial lag of y. Formally it is defined as
a data.frame
containing two columns Ib
and p_sim
containing the local bivariate Moran's I and simulated p-values respectively.
Josiah Parry [email protected]
Anselin, Luc, Ibnu Syabri, and Oleg Smirnov. 2002. “Visualizing Multivariate Spatial Correlation with Dynamically Linked Windows.” In New Tools for Spatial Data Analysis: Proceedings of the Specialist Meeting, edited by Luc Anselin and Sergio Rey. University of California, Santa Barbara: Center for Spatially Integrated Social Science (CSISS).
# load columbus datay columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")) nb <- poly2nb(columbus) listw <- nb2listw(nb) set.seed(1) (res <- localmoran_bv(columbus$CRIME, columbus$INC, listw, nsim = 499)) columbus$hs <- hotspot(res, Prname="Pr(folded) Sim", cutoff=0.05, quadrant.type="pysal", p.adjust="none") if (require("tmap", quietly=TRUE)) { tmap4 <- packageVersion("tmap") >= "3.99" if (tmap4) { tm_shape(columbus) + tm_polygons(fill="hs", fill.scale=tm_scale(values="brewer.set3"), fill.legend=tm_legend(position=tm_pos_in("left", "top"), frame=FALSE, item.r=0), lwd=0.01) } else { tm_shape(columbus) + tm_fill("hs") } } moran.plot(x=columbus$CRIME, y=columbus$INC, listw=listw)
# load columbus datay columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")) nb <- poly2nb(columbus) listw <- nb2listw(nb) set.seed(1) (res <- localmoran_bv(columbus$CRIME, columbus$INC, listw, nsim = 499)) columbus$hs <- hotspot(res, Prname="Pr(folded) Sim", cutoff=0.05, quadrant.type="pysal", p.adjust="none") if (require("tmap", quietly=TRUE)) { tmap4 <- packageVersion("tmap") >= "3.99" if (tmap4) { tm_shape(columbus) + tm_polygons(fill="hs", fill.scale=tm_scale(values="brewer.set3"), fill.legend=tm_legend(position=tm_pos_in("left", "top"), frame=FALSE, item.r=0), lwd=0.01) } else { tm_shape(columbus) + tm_fill("hs") } } moran.plot(x=columbus$CRIME, y=columbus$INC, listw=listw)
localmoran.exact
provides exact local Moran's Ii tests under the null hypothesis, while localmoran.exact.alt
provides exact local Moran's Ii tests under the alternative hypothesis. In this case, the model may be a fitted model derived from a model fitted by spatialreg::errorsarlm
, with the covariance matrix can be passed through the Omega=
argument.
localmoran.exact(model, select, nb, glist = NULL, style = "W", zero.policy = NULL, alternative = "two.sided", spChk = NULL, resfun = weighted.residuals, save.Vi = FALSE, useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) localmoran.exact.alt(model, select, nb, glist = NULL, style = "W", zero.policy = NULL, alternative = "two.sided", spChk = NULL, resfun = weighted.residuals, Omega = NULL, save.Vi = FALSE, save.M = FALSE, useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) ## S3 method for class 'localmoranex' print(x, ...) ## S3 method for class 'localmoranex' as.data.frame(x, row.names=NULL, optional=FALSE, ...)
localmoran.exact(model, select, nb, glist = NULL, style = "W", zero.policy = NULL, alternative = "two.sided", spChk = NULL, resfun = weighted.residuals, save.Vi = FALSE, useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) localmoran.exact.alt(model, select, nb, glist = NULL, style = "W", zero.policy = NULL, alternative = "two.sided", spChk = NULL, resfun = weighted.residuals, Omega = NULL, save.Vi = FALSE, save.M = FALSE, useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) ## S3 method for class 'localmoranex' print(x, ...) ## S3 method for class 'localmoranex' as.data.frame(x, row.names=NULL, optional=FALSE, ...)
model |
an object of class |
select |
an integer vector of the id. numbers of zones to be tested; if missing, all zones |
nb |
a list of neighbours of class |
glist |
a list of general weights corresponding to neighbours |
style |
can take values W, B, C, and S |
zero.policy |
default NULL, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA |
alternative |
a character string specifying the alternative hypothesis, must be one of greater (default), less or two.sided. |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
resfun |
default: weighted.residuals; the function to be used to extract residuals from the |
Omega |
A SAR process matrix may be passed in to test an alternative hypothesis, for example |
save.Vi |
if TRUE, return the star-shaped weights lists for each zone tested |
save.M |
if TRUE, save a list of left and right M products |
useTP |
default FALSE, if TRUE, use truncation point in integration rather than upper=Inf, see Tiefelsdorf (2000), eq. 6.7, p.69 |
truncErr |
when useTP=TRUE, pass truncation error to truncation point function |
zeroTreat |
when useTP=TRUE, pass zero adjustment to truncation point function |
x |
object to be printed |
row.names |
ignored argument to |
optional |
ignored argument to |
... |
arguments to be passed through |
A list with class localmoranex
containing "select" lists, each with
class moranex
with the following components:
statistic |
the value of the exact standard deviate of global Moran's I. |
p.value |
the p-value of the test. |
estimate |
the value of the observed local Moran's Ii. |
method |
a character string giving the method used. |
alternative |
a character string describing the alternative hypothesis. |
gamma |
eigenvalues (two extreme values for null, vector for alternative) |
oType |
usually set to "E", but set to "N" if the integration leads to an out of domain value for |
data.name |
a character string giving the name(s) of the data. |
df |
degrees of freedom |
i |
zone tested |
Vi |
zone tested |
When the alternative is being tested, a list of left and right M products in attribute M.
Markus Reder and Roger Bivand
Bivand RS, Müller W, Reder M (2009) Power calculations for global and local Moran’s I. Comput Stat Data Anal 53:2859–2872; Bivand RS, Wong DWS 2018 Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748 doi:10.1007/s11749-018-0599-x
lm.morantest.exact
, localmoran.sad
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) localmoran.sad(e.lm, nb=eire.nb) localmoran.exact(e.lm, nb=eire.nb) localmoran.exact(e.lm, nb=eire.nb, useTP=TRUE) run <- FALSE if (requireNamespace("spatialreg", quietly=TRUE)) run <- TRUE if (run) { e.errorsar <- spatialreg::errorsarlm(OWNCONS ~ ROADACC, data=eire, listw=nb2listw(eire.nb)) lm.target <- lm(e.errorsar$tary ~ e.errorsar$tarX - 1) localmoran.exact.alt(lm.target, nb=eire.nb) } if (run) { Omega <- spatialreg::invIrW(nb2listw(eire.nb), rho=e.errorsar$lambda) Omega1 <- tcrossprod(Omega) localmoran.exact.alt(lm.target, nb=eire.nb, Omega=Omega1) } if (run) { localmoran.exact.alt(lm.target, nb=eire.nb, Omega=Omega1, useTP=TRUE) }
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) localmoran.sad(e.lm, nb=eire.nb) localmoran.exact(e.lm, nb=eire.nb) localmoran.exact(e.lm, nb=eire.nb, useTP=TRUE) run <- FALSE if (requireNamespace("spatialreg", quietly=TRUE)) run <- TRUE if (run) { e.errorsar <- spatialreg::errorsarlm(OWNCONS ~ ROADACC, data=eire, listw=nb2listw(eire.nb)) lm.target <- lm(e.errorsar$tary ~ e.errorsar$tarX - 1) localmoran.exact.alt(lm.target, nb=eire.nb) } if (run) { Omega <- spatialreg::invIrW(nb2listw(eire.nb), rho=e.errorsar$lambda) Omega1 <- tcrossprod(Omega) localmoran.exact.alt(lm.target, nb=eire.nb, Omega=Omega1) } if (run) { localmoran.exact.alt(lm.target, nb=eire.nb, Omega=Omega1, useTP=TRUE) }
The function implements Tiefelsdorf's application of the Saddlepoint
approximation to local Moran's Ii's reference distribution. If the model
object is of class "lm", global independence is assumed; if of class
"sarlm", global dependence is assumed to be represented by the spatial
parameter of that model. Tests are reported separately for each zone
selected, and may be summarised using summary.localmoransad
.
Values of local Moran's Ii agree with those from localmoran()
,
but in that function, the standard deviate - here the Saddlepoint
approximation - is based on the randomisation assumption.
localmoran.sad(model, select, nb, glist=NULL, style="W", zero.policy=NULL, alternative="two.sided", spChk=NULL, resfun=weighted.residuals, save.Vi=FALSE, tol = .Machine$double.eps^0.5, maxiter = 1000, tol.bounds=0.0001, save.M=FALSE, Omega = NULL) ## S3 method for class 'localmoransad' print(x, ...) ## S3 method for class 'localmoransad' summary(object, ...) ## S3 method for class 'summary.localmoransad' print(x, ...) listw2star(listw, ireg, style, n, D, a, zero.policy=attr(listw, "zero.policy"))
localmoran.sad(model, select, nb, glist=NULL, style="W", zero.policy=NULL, alternative="two.sided", spChk=NULL, resfun=weighted.residuals, save.Vi=FALSE, tol = .Machine$double.eps^0.5, maxiter = 1000, tol.bounds=0.0001, save.M=FALSE, Omega = NULL) ## S3 method for class 'localmoransad' print(x, ...) ## S3 method for class 'localmoransad' summary(object, ...) ## S3 method for class 'summary.localmoransad' print(x, ...) listw2star(listw, ireg, style, n, D, a, zero.policy=attr(listw, "zero.policy"))
model |
an object of class |
select |
an integer vector of the id. numbers of zones to be tested; if missing, all zones |
nb |
a list of neighbours of class |
glist |
a list of general weights corresponding to neighbours |
style |
can take values W, B, C, and S |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of greater (default), less or two.sided. |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
resfun |
default: weighted.residuals; the function to be used to extract residuals from the |
save.Vi |
if TRUE, return the star-shaped weights lists for each zone tested |
tol |
the desired accuracy (convergence tolerance) for |
maxiter |
the maximum number of iterations for |
tol.bounds |
offset from bounds for |
save.M |
if TRUE, save a list of left and right M products in a list for the conditional tests, or a list of the regression model matrix components |
Omega |
A SAR process matrix may be passed in to test an alternative hypothesis, for example |
x |
object to be printed |
object |
object to be summarised |
... |
arguments to be passed through |
listw |
a |
ireg |
a zone number |
n |
internal value depending on listw and style |
D |
internal value depending on listw and style |
a |
internal value depending on listw and style |
The function implements the analytical eigenvalue calculation together with trace shortcuts given or suggested in Tiefelsdorf (2002), partly following remarks by J. Keith Ord, and uses the Saddlepoint analytical solution from Tiefelsdorf's SPSS code.
If a histogram of the probability values of the saddlepoint estimate for the assumption of global independence is not approximately flat, the assumption is probably unjustified, and re-estimation with global dependence is recommended.
No n by n matrices are needed at any point for the test assuming no global dependence, the star-shaped weights matrices being handled as listw lists. When the test is made on residuals from a spatial regression, taking a global process into account. n by n matrices are necessary, and memory constraints may be reached for large lattices.
A list with class localmoransad
containing "select" lists, each with
class moransad
with the following components:
statistic |
the value of the saddlepoint approximation of the standard deviate of local Moran's Ii. |
p.value |
the p-value of the test. |
estimate |
the value of the observed local Moran's Ii. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data. |
internal1 |
Saddlepoint omega, r and u |
df |
degrees of freedom |
tau |
maximum and minimum analytical eigenvalues |
i |
zone tested |
Roger Bivand [email protected]
Tiefelsdorf, M. 2002 The Saddlepoint approximation of Moran's I and local Moran's Ii reference distributions and their numerical evaluation. Geographical Analysis, 34, pp. 187–206.
localmoran
, lm.morantest
,
lm.morantest.sad
, errorsarlm
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) lw <- nb2listw(eire.nb) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) e.locmor <- summary(localmoran.sad(e.lm, nb=eire.nb)) e.locmor mean(e.locmor[,1]) sum(e.locmor[,1])/Szero(lw) lm.morantest(e.lm, lw) # note equality for mean() only when the sum of weights equals # the number of observations (thanks to Juergen Symanzik) hist(e.locmor[,"Pr. (Sad)"]) e.wlm <- lm(OWNCONS ~ ROADACC, data=eire, weights=RETSALE) e.locmorw1 <- summary(localmoran.sad(e.wlm, nb=eire.nb, resfun=weighted.residuals)) e.locmorw1 e.locmorw2 <- summary(localmoran.sad(e.wlm, nb=eire.nb, resfun=rstudent)) e.locmorw2 run <- FALSE if (requireNamespace("spatialreg", quietly=TRUE)) run <- TRUE if (run) { e.errorsar <- spatialreg::errorsarlm(OWNCONS ~ ROADACC, data=eire, listw=lw) if (packageVersion("spatialreg") < "1.1.7") spatialreg::print.sarlm(e.errorsar) else print(e.errorsar) } if (run) { lm.target <- lm(e.errorsar$tary ~ e.errorsar$tarX - 1) Omega <- tcrossprod(spatialreg::invIrW(lw, rho=e.errorsar$lambda)) e.clocmor <- summary(localmoran.sad(lm.target, nb=eire.nb, Omega=Omega)) e.clocmor } if (run) { hist(e.clocmor[,"Pr. (Sad)"]) }
eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1]) row.names(eire) <- as.character(eire$names) eire.nb <- poly2nb(eire) lw <- nb2listw(eire.nb) e.lm <- lm(OWNCONS ~ ROADACC, data=eire) e.locmor <- summary(localmoran.sad(e.lm, nb=eire.nb)) e.locmor mean(e.locmor[,1]) sum(e.locmor[,1])/Szero(lw) lm.morantest(e.lm, lw) # note equality for mean() only when the sum of weights equals # the number of observations (thanks to Juergen Symanzik) hist(e.locmor[,"Pr. (Sad)"]) e.wlm <- lm(OWNCONS ~ ROADACC, data=eire, weights=RETSALE) e.locmorw1 <- summary(localmoran.sad(e.wlm, nb=eire.nb, resfun=weighted.residuals)) e.locmorw1 e.locmorw2 <- summary(localmoran.sad(e.wlm, nb=eire.nb, resfun=rstudent)) e.locmorw2 run <- FALSE if (requireNamespace("spatialreg", quietly=TRUE)) run <- TRUE if (run) { e.errorsar <- spatialreg::errorsarlm(OWNCONS ~ ROADACC, data=eire, listw=lw) if (packageVersion("spatialreg") < "1.1.7") spatialreg::print.sarlm(e.errorsar) else print(e.errorsar) } if (run) { lm.target <- lm(e.errorsar$tary ~ e.errorsar$tarX - 1) Omega <- tcrossprod(spatialreg::invIrW(lw, rho=e.errorsar$lambda)) e.clocmor <- summary(localmoran.sad(lm.target, nb=eire.nb, Omega=Omega)) e.clocmor } if (run) { hist(e.clocmor[,"Pr. (Sad)"]) }
Local spatial heteroscedasticity is calculated for each location based on the spatial weights object used. The statistic is:
with
and
Its expectation and variance are given in Ord & Getis (2012). The exponent a allows for investigating different types of mean dispersal.
LOSH(x, listw, a=2, var_hi=TRUE, zero.policy=attr(listw, "zero.policy"), na.action=na.fail, spChk=NULL)
LOSH(x, listw, a=2, var_hi=TRUE, zero.policy=attr(listw, "zero.policy"), na.action=na.fail, spChk=NULL)
x |
a numeric vector of the same length as the neighbours list in listw |
listw |
a |
a |
the exponent applied to the local residuals; the default value of 2 leads to a measure of heterogeneity in the spatial variance |
var_hi |
default TRUE, the moments and the test statistics are calculated for each location; if FALSE, only the plain LOSH measures, |
zero.policy |
default |
na.action |
a function (default |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
In addition to the LOSH measure, the values returned include local spatially weighted mean values and local residuals
estimated about these means. These values facilitate the interpretation of LOSH values. Further, if specified through
var_hi
, the statistical moments and the test statistics as proposed by Ord & Getis (2012) are also calculated and returned.
Hi |
LOSH statistic |
E.Hi |
(optional) expectation of LOSH |
Var.Hi |
(optional) variance of LOSH |
Z.Hi |
(optional) the approximately Chi-square distributed test statistics |
x_bar_i |
local spatially weighted mean values |
ei |
residuals about local spatially weighted mean values |
René Westerholt [email protected]
Ord, J. K., & Getis, A. 2012. Local spatial heteroscedasticity (LOSH), The Annals of Regional Science, 48 (2), 529–539.
data(boston, package="spData") resLOSH <- LOSH(boston.c$NOX, nb2listw(boston.soi)) hist(resLOSH[,"Hi"]) mean(resLOSH[,"Hi"])
data(boston, package="spData") resLOSH <- LOSH(boston.c$NOX, nb2listw(boston.soi)) hist(resLOSH[,"Hi"]) mean(resLOSH[,"Hi"])
The function implements the chi-square based test statistic for local spatial heteroscedasticity (LOSH) as proposed by Ord & Getis (2012).
LOSH.cs(x, listw, zero.policy = attr(listw, "zero.policy"), na.action = na.fail, p.adjust.method = "none", spChk = NULL)
LOSH.cs(x, listw, zero.policy = attr(listw, "zero.policy"), na.action = na.fail, p.adjust.method = "none", spChk = NULL)
x |
a numeric vector of the same length as the neighbours list in listw |
listw |
a |
zero.policy |
default |
na.action |
a function (default |
p.adjust.method |
a character string specifying the probability value adjustment for multiple tests, default "none"; see |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
The test uses a = 2 (see LOSH
) because chi-square based inference is not applicable with other exponents. The function makes use of LOSH
in its calculations.
Hi |
LOSH statistic |
E.Hi |
expectation of LOSH |
Var.Hi |
variance of LOSH |
Z.Hi |
the approximately chi-square distributed test statistics |
x_bar_i |
local spatially weighted mean values |
ei |
residuals about local spatially weighted mean values |
Pr() |
p-values for |
René Westerholt [email protected]
Ord, J. K., & Getis, A. 2012. Local spatial heteroscedasticity (LOSH), The Annals of Regional Science, 48 (2), 529–539.
data(boston, package="spData") resLOSH <- LOSH.cs(boston.c$NOX, nb2listw(boston.soi)) hist(resLOSH[,"Hi"]) mean(resLOSH[,"Hi"])
data(boston, package="spData") resLOSH <- LOSH.cs(boston.c$NOX, nb2listw(boston.soi)) hist(resLOSH[,"Hi"]) mean(resLOSH[,"Hi"])
The function draws inferences about local spatial heteroscedasticity (LOSH) by means of the randomisation-based Monte-Carlo bootstrap proposed by Xu et al. (2014).
LOSH.mc(x, listw, a = 2, nsim = 99, zero.policy = attr(listw, "zero.policy"), na.action = na.fail, spChk = NULL, adjust.n = TRUE, p.adjust.method = "none")
LOSH.mc(x, listw, a = 2, nsim = 99, zero.policy = attr(listw, "zero.policy"), na.action = na.fail, spChk = NULL, adjust.n = TRUE, p.adjust.method = "none")
x |
a numeric vector of the same length as the neighbours list in listw |
listw |
a |
a |
the exponent applied to the local residuals; the default value of 2 leads to a measure of heterogeneity in the spatial variance |
nsim |
the number of randomisations used in the bootstrap |
zero.policy |
default |
na.action |
a function (default |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted |
p.adjust.method |
a character string specifying the probability value adjustment for multiple tests, default "none"; see |
The test calculates LOSH (see LOSH
) and estimates pseudo p-values from a conditional bootstrap. Thereby, the i-th value in each location is held fixed, whereas all other values are permuted nsim
times over all other spatial units.
Hi |
LOSH statistic |
E.Hi |
expectation of LOSH |
Var.Hi |
variance of LOSH |
Z.Hi |
the approximately chi-square distributed test statistics |
x_bar_i |
local spatially weighted mean values |
ei |
residuals about local spatially weighted mean values |
Pr() |
p-values for |
René Westerholt [email protected]
Ord, J. K., & Getis, A. 2012. Local spatial heteroscedasticity (LOSH), The Annals of Regional Science, 48 (2), 529–539; Xu, M., Mei, C. L., & Yan, N. 2014. A note on the null distribution of the local spatial heteroscedasticity (LOSH) statistic. The Annals of Regional Science, 52 (3), 697–710.
data(columbus, package="spData") resLOSH_mc <- LOSH.mc(columbus$CRIME, nb2listw(col.gal.nb), 2, 100) summary(resLOSH_mc) resLOSH_cs <- LOSH.cs(columbus$CRIME, nb2listw(col.gal.nb)) summary(resLOSH_cs) plot(resLOSH_mc[,"Pr()"], resLOSH_cs[,"Pr()"])
data(columbus, package="spData") resLOSH_mc <- LOSH.mc(columbus$CRIME, nb2listw(col.gal.nb), 2, 100) summary(resLOSH_mc) resLOSH_cs <- LOSH.cs(columbus$CRIME, nb2listw(col.gal.nb)) summary(resLOSH_cs) plot(resLOSH_mc[,"Pr()"], resLOSH_cs[,"Pr()"])
The function converts a square spatial weights matrix, optionally a sparse matrix to a weights list object, optionally adding region IDs from the row names of the matrix, as a sequence of numbers 1:nrow(x), or as given as an argument. The style can be imposed by rebuilting the weights list object internally.
mat2listw(x, row.names = NULL, style=NULL, zero.policy = NULL)
mat2listw(x, row.names = NULL, style=NULL, zero.policy = NULL)
x |
A square non-negative matrix with no NAs representing spatial weights; may be a matrix of class “sparseMatrix” |
row.names |
row names to use for region IDs |
style |
default NULL, missing, set to "M" and warning given; if not "M", passed to |
zero.policy |
default NULL, use global option value; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors |
A listw
object with the following members:
style |
"M", meaning matrix style, underlying style unknown, or assigned style argument in rebuilt object |
neighbours |
the derived neighbours list |
weights |
the weights for the neighbours derived from the matrix |
Roger Bivand [email protected]
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col005 <- dnearneigh(st_coordinates(st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)), 0, 0.5, as.character(columbus$NEIGNO)) summary(col005) col005.w.mat <- nb2mat(col005, style="W", zero.policy=TRUE) try(col005.w.b <- mat2listw(col005.w.mat, style="W")) col005.w.b <- mat2listw(col005.w.mat, style="W", zero.policy=TRUE) summary(col005.w.b$neighbours) diffnb(col005, col005.w.b$neighbours) col005.w.mat.3T <- kronecker(diag(3), col005.w.mat) col005.w.b.3T <- mat2listw(col005.w.mat.3T, style="W", zero.policy=TRUE) summary(col005.w.b.3T$neighbours) run <- FALSE if (require("spatialreg", quiet=TRUE)) run <- TRUE if (run) { W <- as(nb2listw(col005, style="W", zero.policy=TRUE), "CsparseMatrix") try(col005.spM <- mat2listw(W)) col005.spM <- mat2listw(W, style="W", zero.policy=TRUE) summary(col005.spM$neighbours) } if (run) { diffnb(col005, col005.spM$neighbours) } if (run && require("Matrix", quiet=TRUE)) { IW <- kronecker(Diagonal(3), W) col005.spM.3T <- mat2listw(as(IW, "CsparseMatrix"), style="W", zero.policy=TRUE) summary(col005.spM.3T$neighbours) }
columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE) col005 <- dnearneigh(st_coordinates(st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)), 0, 0.5, as.character(columbus$NEIGNO)) summary(col005) col005.w.mat <- nb2mat(col005, style="W", zero.policy=TRUE) try(col005.w.b <- mat2listw(col005.w.mat, style="W")) col005.w.b <- mat2listw(col005.w.mat, style="W", zero.policy=TRUE) summary(col005.w.b$neighbours) diffnb(col005, col005.w.b$neighbours) col005.w.mat.3T <- kronecker(diag(3), col005.w.mat) col005.w.b.3T <- mat2listw(col005.w.mat.3T, style="W", zero.policy=TRUE) summary(col005.w.b.3T$neighbours) run <- FALSE if (require("spatialreg", quiet=TRUE)) run <- TRUE if (run) { W <- as(nb2listw(col005, style="W", zero.policy=TRUE), "CsparseMatrix") try(col005.spM <- mat2listw(W)) col005.spM <- mat2listw(W, style="W", zero.policy=TRUE) summary(col005.spM$neighbours) } if (run) { diffnb(col005, col005.spM$neighbours) } if (run && require("Matrix", quiet=TRUE)) { IW <- kronecker(Diagonal(3), W) col005.spM.3T <- mat2listw(as(IW, "CsparseMatrix"), style="W", zero.policy=TRUE) summary(col005.spM.3T$neighbours) }
A simple function to compute Moran's I, called by moran.test
and moran.mc
;
moran(x, listw, n, S0, zero.policy=attr(listw, "zero.policy"), NAOK=FALSE)
moran(x, listw, n, S0, zero.policy=attr(listw, "zero.policy"), NAOK=FALSE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
n |
number of zones |
S0 |
global sum of weights |
zero.policy |
default |
NAOK |
if 'TRUE' then any 'NA' or 'NaN' or 'Inf' values in x are passed on to the foreign function. If 'FALSE', the presence of 'NA' or 'NaN' or 'Inf' values is regarded as an error. |
a list of
I |
Moran's I |
K |
sample kurtosis of x |
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 17.
data(oldcol) col.W <- nb2listw(COL.nb, style="W") crime <- COL.OLD$CRIME str(moran(crime, col.W, length(COL.nb), Szero(col.W))) is.na(crime) <- sample(1:length(crime), 10) str(moran(crime, col.W, length(COL.nb), Szero(col.W), NAOK=TRUE))
data(oldcol) col.W <- nb2listw(COL.nb, style="W") crime <- COL.OLD$CRIME str(moran(crime, col.W, length(COL.nb), Szero(col.W))) is.na(crime) <- sample(1:length(crime), 10) str(moran(crime, col.W, length(COL.nb), Szero(col.W), NAOK=TRUE))
Given two continuous numeric variables, calculate the bivariate Moran's I. See details for more.
moran_bv(x, y, listw, nsim = 499, scale = TRUE)
moran_bv(x, y, listw, nsim = 499, scale = TRUE)
x |
a numeric vector of same length as |
y |
a numeric vector of same length as |
listw |
a listw object for example as created by |
nsim |
the number of simulations to run. |
scale |
default |
The Global Bivariate Moran is defined as
It is important to note that this is a measure of autocorrelation of X
with the spatial lag of Y. As such, the resultant measure may overestimate the amount of
spatial autocorrelation which may be a product of the inherent correlation of X and Y. The output object is of class "boot"
, so that plots and confidence intervals are available using appropriate methods.
An object of class "boot"
, with the observed statistic in component t0
.
Josiah Parry [email protected]
Wartenberg, D. (1985), Multivariate Spatial Correlation: A Method for Exploratory Geographical Analysis. Geographical Analysis, 17: 263-283. doi:10.1111/j.1538-4632.1985.tb00849.x
data(boston, package = "spData") x <- boston.c$CRIM y <- boston.c$NOX listw <- nb2listw(boston.soi) set.seed(1) res_xy <- moran_bv(x, y, listw, nsim=499) res_xy$t0 boot::boot.ci(res_xy, conf=c(0.99, 0.95, 0.9), type="basic") plot(res_xy) set.seed(1) lee_xy <- lee.mc(x, y, listw, nsim=499, return_boot=TRUE) lee_xy$t0 boot::boot.ci(lee_xy, conf=c(0.99, 0.95, 0.9), type="basic") plot(lee_xy) set.seed(1) res_yx <- moran_bv(y, x, listw, nsim=499) res_yx$t0 boot::boot.ci(res_yx, conf=c(0.99, 0.95, 0.9), type="basic") plot(res_yx) set.seed(1) lee_yx <- lee.mc(y, x, listw, nsim=499, return_boot=TRUE) lee_yx$t0 boot::boot.ci(lee_yx, conf=c(0.99, 0.95, 0.9), type="basic") plot(lee_yx)
data(boston, package = "spData") x <- boston.c$CRIM y <- boston.c$NOX listw <- nb2listw(boston.soi) set.seed(1) res_xy <- moran_bv(x, y, listw, nsim=499) res_xy$t0 boot::boot.ci(res_xy, conf=c(0.99, 0.95, 0.9), type="basic") plot(res_xy) set.seed(1) lee_xy <- lee.mc(x, y, listw, nsim=499, return_boot=TRUE) lee_xy$t0 boot::boot.ci(lee_xy, conf=c(0.99, 0.95, 0.9), type="basic") plot(lee_xy) set.seed(1) res_yx <- moran_bv(y, x, listw, nsim=499) res_yx$t0 boot::boot.ci(res_yx, conf=c(0.99, 0.95, 0.9), type="basic") plot(res_yx) set.seed(1) lee_yx <- lee.mc(y, x, listw, nsim=499, return_boot=TRUE) lee_yx$t0 boot::boot.ci(lee_yx, conf=c(0.99, 0.95, 0.9), type="basic") plot(lee_yx)
A permutation test for Moran's I statistic calculated by using nsim random permutations of x for the given spatial weighting scheme, to establish the rank of the observed statistic in relation to the nsim simulated values.
moran.mc(x, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", na.action=na.fail, spChk=NULL, return_boot=FALSE, adjust.n=TRUE)
moran.mc(x, listw, nsim, zero.policy=attr(listw, "zero.policy"), alternative="greater", na.action=na.fail, spChk=NULL, return_boot=FALSE, adjust.n=TRUE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
nsim |
number of permutations |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of "greater" (default), "two.sided", or "less". |
na.action |
a function (default |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
return_boot |
return an object of class |
adjust.n |
default TRUE, if FALSE the number of observations is not adjusted for no-neighbour observations, if TRUE, the number of observations is adjusted |
A list with class htest
and mc.sim
containing the following components:
statistic |
the value of the observed Moran's I. |
parameter |
the rank of the observed Moran's I. |
p.value |
the pseudo p-value of the test. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string giving the method used. |
data.name |
a character string giving the name(s) of the data, and the number of simulations. |
res |
nsim simulated values of statistic, final value is observed statistic |
Roger Bivand [email protected]
Cliff, A. D., Ord, J. K. 1981 Spatial processes, Pion, p. 63-5.
data(oldcol) colw <- nb2listw(COL.nb, style="W") nsim <- 99 set.seed(1234) sim1 <- moran.mc(COL.OLD$CRIME, listw=colw, nsim=nsim) sim1 mean(sim1$res[1:nsim]) var(sim1$res[1:nsim]) summary(sim1$res[1:nsim]) colold.lags <- nblag(COL.nb, 3) set.seed(1234) sim2 <- moran.mc(COL.OLD$CRIME, nb2listw(colold.lags[[2]], style="W"), nsim=nsim) summary(sim2$res[1:nsim]) sim3 <- moran.mc(COL.OLD$CRIME, nb2listw(colold.lags[[3]], style="W"), nsim=nsim) summary(sim3$res[1:nsim]) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) try(moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.fail)) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.omit) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.omit) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.exclude) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.exclude) try(moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.pass))
data(oldcol) colw <- nb2listw(COL.nb, style="W") nsim <- 99 set.seed(1234) sim1 <- moran.mc(COL.OLD$CRIME, listw=colw, nsim=nsim) sim1 mean(sim1$res[1:nsim]) var(sim1$res[1:nsim]) summary(sim1$res[1:nsim]) colold.lags <- nblag(COL.nb, 3) set.seed(1234) sim2 <- moran.mc(COL.OLD$CRIME, nb2listw(colold.lags[[2]], style="W"), nsim=nsim) summary(sim2$res[1:nsim]) sim3 <- moran.mc(COL.OLD$CRIME, nb2listw(colold.lags[[3]], style="W"), nsim=nsim) summary(sim3$res[1:nsim]) crime <- COL.OLD$CRIME is.na(crime) <- sample(1:length(crime), 10) try(moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.fail)) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.omit) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.omit) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, na.action=na.exclude) moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, zero.policy=TRUE, return_boot=TRUE, na.action=na.exclude) try(moran.mc(crime, nb2listw(COL.nb, style="W"), nsim=99, na.action=na.pass))
A plot of spatial data against its spatially lagged values, augmented by reporting the summary of influence measures for the linear relationship between the data and the lag. If zero policy is TRUE, such observations are also marked if they occur.
moran.plot(x, listw, y=NULL, zero.policy=attr(listw, "zero.policy"), spChk=NULL, labels=NULL, xlab=NULL, ylab=NULL, quiet=NULL, plot=TRUE, return_df=TRUE, ...)
moran.plot(x, listw, y=NULL, zero.policy=attr(listw, "zero.policy"), spChk=NULL, labels=NULL, xlab=NULL, ylab=NULL, quiet=NULL, plot=TRUE, return_df=TRUE, ...)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
y |
an optional numeric vector the same length as the neighbours list in listw for a bi-variate plot |
zero.policy |
default |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |
labels |
character labels for points with high influence measures, if set to FALSE, no labels are plotted for points with large influence |
xlab |
label for x axis |
ylab |
label for x axis |
quiet |
default NULL, use !verbose global option value; if TRUE, output of summary of influence object suppressed |
plot |
default TRUE, if false, plotting is suppressed |
return_df |
default TRUE, invisibly return a data.frame object; if FALSE invisibly return an influence measures object |
... |
further graphical parameters as in |
The function returns a data.frame object with coordinates and influence measures if return_df
is TRUE, or an influence object from influence.measures
.
Roger Bivand [email protected]
Anselin, L. 1996. The Moran scatterplot as an ESDA tool to assess local instability in spatial association. pp. 111–125 in M. M. Fischer, H. J. Scholten and D. Unwin (eds) Spatial analytical perspectives on GIS, London, Taylor and Francis; Anselin, L. 1995. Local indicators of spatial association, Geographical Analysis, 27, 93–115
localmoran
, influence.measures
data(afcon, package="spData") mp <- moran.plot(afcon$totcon, nb2listw(paper.nb), labels=as.character(afcon$name), pch=19) moran.plot(as.vector(scale(afcon$totcon)), nb2listw(paper.nb), labels=as.character(afcon$name), xlim=c(-2, 4), ylim=c(-2,4), pch=19) if (require(ggplot2, quietly=TRUE)) { xname <- attr(mp, "xname") ggplot(mp, aes(x=x, y=wx)) + geom_point(shape=1) + geom_smooth(formula=y ~ x, method="lm") + geom_hline(yintercept=mean(mp$wx), lty=2) + geom_vline(xintercept=mean(mp$x), lty=2) + theme_minimal() + geom_point(data=mp[mp$is_inf,], aes(x=x, y=wx), shape=9) + geom_text(data=mp[mp$is_inf,], aes(x=x, y=wx, label=labels, vjust=1.5)) + xlab(xname) + ylab(paste0("Spatially lagged ", xname)) } columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")) nb <- poly2nb(columbus) listw <- nb2listw(nb) moran.plot(x=columbus$CRIME, y=columbus$INC, listw=listw) moran.plot(x=columbus$INC, y=columbus$CRIME, listw=listw)
data(afcon, package="spData") mp <- moran.plot(afcon$totcon, nb2listw(paper.nb), labels=as.character(afcon$name), pch=19) moran.plot(as.vector(scale(afcon$totcon)), nb2listw(paper.nb), labels=as.character(afcon$name), xlim=c(-2, 4), ylim=c(-2,4), pch=19) if (require(ggplot2, quietly=TRUE)) { xname <- attr(mp, "xname") ggplot(mp, aes(x=x, y=wx)) + geom_point(shape=1) + geom_smooth(formula=y ~ x, method="lm") + geom_hline(yintercept=mean(mp$wx), lty=2) + geom_vline(xintercept=mean(mp$x), lty=2) + theme_minimal() + geom_point(data=mp[mp$is_inf,], aes(x=x, y=wx), shape=9) + geom_text(data=mp[mp$is_inf,], aes(x=x, y=wx, label=labels, vjust=1.5)) + xlab(xname) + ylab(paste0("Spatially lagged ", xname)) } columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")) nb <- poly2nb(columbus) listw <- nb2listw(nb) moran.plot(x=columbus$CRIME, y=columbus$INC, listw=listw) moran.plot(x=columbus$INC, y=columbus$CRIME, listw=listw)
Moran's test for spatial autocorrelation using a spatial weights matrix in weights list form. The assumptions underlying the test are sensitive to the form of the graph of neighbour relationships and other factors, and results may be checked against those of moran.mc
permutations.
moran.test(x, listw, randomisation=TRUE, zero.policy=attr(listw, "zero.policy"), alternative="greater", rank = FALSE, na.action=na.fail, spChk=NULL, adjust.n=TRUE, drop.EI2=FALSE)
moran.test(x, listw, randomisation=TRUE, zero.policy=attr(listw, "zero.policy"), alternative="greater", rank = FALSE, na.action=na.fail, spChk=NULL, adjust.n=TRUE, drop.EI2=FALSE)
x |
a numeric vector the same length as the neighbours list in listw |
listw |
a |
randomisation |
variance of I calculated under the assumption of randomisation, if FALSE normality |
zero.policy |
default |
alternative |
a character string specifying the alternative hypothesis, must be one of greater (default), less or two.sided. |
rank |
logical value - default FALSE for continuous variables, if TRUE, uses the adaptation of Moran's I for ranks suggested by Cliff and Ord (1981, p. 46) |
na.action |
a function (default |
spChk |
should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use |