Package 'spdep' reference manual

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>. 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], 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], Gianfranco Piras [ctb], Markus Reder [ctb], Jeff Sauer [ctb], Michael Tiefelsdorf [ctb], René Westerholt [ctb], Levi Wolf [ctb], Danlin Yu [ctb]
Maintainer:	Roger Bivand <[email protected]>
License:	GPL (>= 2)
Version:	1.3-4
Built:	2024-04-29 14:32:27 UTC
Source:	https://github.com/r-spatial/spdep

Compute Local Geary statistic

Description

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.

Usage

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)

Arguments

`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 `listw` object created for example by `nb2listw`.
`data`	Used when a formula is provided. A matrix or data frame containing the variables in the formula `formula`.
`nsim`	The number of simulations to be used for permutation test.
`alternative`	A character defining the alternative hypothesis. Must be one of `"two.sided"`, `"less"` or `"greater"`.
`...`	other arguments passed to methods.
`zero.policy`	default `attr(listw, "zero.policy")` as set when `listw` was created, if attribute not set, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA.
`iseed`	default NULL, used to set the seed for possible parallel RNGs
`no_repeat_in_row`	default `FALSE`, if `TRUE`, sample conditionally in each row without replacements to avoid duplicate values, https://github.com/r-spatial/spdep/issues/124

Details

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 $c_i = \sum_j w_{ij}(x_i - x_j)^2$ and the multivariate Local Geary is calculated as $c_{k,i} = \sum_{v=1}^{k} c_{v,i}$ as described in Anselin (2019).

Value

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 `pnorm()` using standard deviates based on permutation sample means and standard deviations
`Pr() Sim`	`rank()` and `punif()` of observed statistic rank for [0, 1] p-values using `alternative=`
`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 `e1071::skewness()` for the permutation samples underlying the standard deviates
`Kurtosis`	the output of `e1071::kurtosis()` for the permutation samples underlying the standard deviates

Author(s)

Josiah Parry [email protected] and Roger Bivand

References

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 Geary's c. Geogr Anal, 51: 133-150. doi:10.1111/gean.12164

Examples

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
## Not run: 
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) 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) 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) 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) 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) 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) tm_shape(g) + tm_fill("lg_PC1", 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) 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) 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) 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) 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) 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) 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) 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) 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) 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) 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) 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) tm_shape(g) + tm_fill("lmc_vars_spdep", textNA="Insignificant",
    colorNA="gray95") + tm_borders() # rep. Fig. 11
  else plot(g["lmc_vars_spdep"])
}


## End(Not run)
## 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)

Saddlepoint approximation of local Moran's Ii tests

Description

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.

Usage

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"))

Arguments

`model`	an object of class `lm` returned by `lm` (assuming no global spatial autocorrelation), or an object of class `sarlm` returned by a spatial simultaneous autoregressive model fit (assuming global spatial autocorrelation represented by the model spatial coefficient); weights may be specified in the `lm` fit, but offsets should not be used
`select`	an integer vector of the id. numbers of zones to be tested; if missing, all zones
`nb`	a list of neighbours of class `nb`
`glist`	a list of general weights corresponding to neighbours
`style`	can take values W, B, C, and S
`zero.policy`	default `attr(listw, "zero.policy")` as set when `listw` was created, if attribute not set, 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 `get.spChkOption()`
`resfun`	default: weighted.residuals; the function to be used to extract residuals from the `lm` object, may be `residuals`, `weighted.residuals`, `rstandard`, or `rstudent`
`save.Vi`	if TRUE, return the star-shaped weights lists for each zone tested
`tol`	the desired accuracy (convergence tolerance) for `uniroot`
`maxiter`	the maximum number of iterations for `uniroot`
`tol.bounds`	offset from bounds for `uniroot`
`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 `Omega <- invIrW(listw, rho=0.1); Omega <- tcrossprod(Omega)`, `chol()` is taken internally
`x`	object to be printed
`object`	object to be summarised
`...`	arguments to be passed through
`listw`	a `listw` object created for example by `nb2listw`
`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

Details

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.

Value

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

Author(s)

Roger Bivand [email protected]

References

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.

Examples

eire <- st_read(system.file("shapes/eire.shp", package="spData")[1])
row.names(eire) <- as.character(eire$names)
st_crs(eire) <- "+proj=utm +zone=30 +ellps=airy +units=km"
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)"])
}

`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

`dist`	distance measured
`coords`	coordinates between which distance is measured

`z`	the response variable
`xy`	a matrix of coordinates or a SpatialPoints, `sf` or `sfc` points object
`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

`nrow`	number of rows in the grid, may also be an object inheriting from class `"SpatialGrid"` or `"GridTopology"` only in `cell2nb`
`ncol`	number of columns in the grid; if `nrow` or `x` is an object inheriting from class `"SpatialGrid"` or `"GridTopology"`, it may be omitted
`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 `"SpatialGrid"` or `"GridTopology"`, and replaces `nrow` and `ncol`
`i`	vector of vector indices corresponding to rowcol, a matrix with two columns of row, column indices

`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 `ppois` version, if true use original version written from sources and iterating down to `tol`

`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

`x`	an object of class `nb`
`y`	an object of class `nb`
`verbose`	default NULL, use global option value; report regions ids taken from object attribute "region.id" with differences

`x`	matrix of point coordinates, an object inheriting from SpatialPoints or an `"sf"` or `"sfc"` object; if the `"sf"` or `"sfc"` object geometries are in geographical coordinates (`use_s2=FALSE`, `sf::st_is_longlat(x) == TRUE` and `sf::sf_use_s2() == TRUE`), s2 will be used to find the neighbours because it will (we hope) use spatial indexing https://github.com/r-spatial/s2/issues/125 as opposed to the legacy method which uses brute-force (at present s2 also uses brute-force)
`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 `region.id`, default `seq(1, nrow(x))`
`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 `c("GE", "LE")`, (GE: greater than or equal to, LE: less than or equal to) that is the finite and closed interval `[d1, d2]`, `d1 <= x <= d2`. The first element may also be `"GT"` (GT: greater than), the second `"LT"` (LT: less than) for finite, open intervals excluding the bounds; the first bound default was changed from `"GT"` to `"GE"` in release 1.1-7. When creating multiple distance bands, finite, half-open right-closed intervals may be used until the final interval to avoid overlapping on bounds: `"GE", "LT"`, that is `[d1, d2)`, `d1 <= x < d2`
`use_kd_tree`	default TRUE, if TRUE, use dbscan `frNN` if available (permitting 3D distances).
`symtest`	Default FALSE; before release 1.1-7, TRUE - run symmetry check on output object, costly with large numbers of points.
`use_s2`	default=`packageVersion("s2") > "1.0.7"`, as of s2 > 1.0-7, distance bound computations use spatial indexing so when `sf::sf_use_s2()` is TRUE, `s2::s2_dwithin_matrix()` will be used for distances on the sphere for `"sf"` or `"sfc"` objects if s2 > 1.0-7.
`k`	default 200, the number of closest points to consider when searching when using `s2::s2_closest_edges()`
`dwithin`	default TRUE, if FALSE, use `s2::s2_closest_edges()`, both if `use_s2=TRUE`, `sf::st_is_longlat(x) == TRUE` and `sf::sf_use_s2() == TRUE`; `s2::s2_dwithin_matrix()` yields the same lists of neighbours as `s2::s2_closest_edges()` is `k=` is set correctly.

`nb`	a neighbours list object of class `nb`
`drop`	either a logical vector the length of `nb`, or a character vector of named regions corresponding to `nb`'s region.id attribute, or an integer vector of region numbers
`sym`	TRUE for removal of both "row" and "column" links, FALSE for only "row" links

`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

`raw`	a numerical vector of raw (crude) rates
`estmm`	a numerical vector of empirical Bayes estimates

`a`	global method of moments phi value
`m`	global method of moments gamma value

`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 `listw` object created for example by `nb2listw`
`nsim`	number of permutations
`zero.policy`	default `attr(listw, "zero.policy")` as set when `listw` was created, if attribute not set, 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), "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 `get.spChkOption()`
`return_boot`	return an object of class `boot` from the equivalent permutation bootstrap rather than an object of class `htest`
`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.

`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

`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 `nb` object of neighbour relationships
`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 `get.spChkOption()`
`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)

`C`	Geary's C
`K`	sample kurtosis of x

`nc`	number of disjoint connected subgraphs
`comp.id`	vector with the indices of the disjoint connected subgraphs that the nodes in `nb.obj` belong to

`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

`obj`	An object of class `localmoran`, `localC` or `localG`
`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 `"fdr"`, the `p.adjust()` methood used, one of `c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")`
`droplevels`	Default `TRUE`, should empty levels of the input cluster factor be dropped
`quadrant.type`	Default `"mean"`, for `"localmoran"` objects only, can be `c("mean", "median", "pysal")` to partition the Moran scatterplot; `"mean"` partitions on the means of the variable and its spatial lag, `"median"` on medians of the variable and its spatial lag, `"pysal"` at zero for the centred variable and its spatial lag
`...`	other arguments passed to methods.

`raw`	a numerical vector of raw (crude) rates
`est`	a numerical vector of local empirical Bayes estimates

`a`	a numerical vector of local phi values
`m`	a numerical vector of local gamma values

`name`	an object of class `nb`
`coords`	matrix of region point coordinates; if missing and polys= inherits from `SpatialPolygons`, the label points of that object are used
`polys`	if polygon boundaries supplied, will be used as background; must inherit from `SpatialPolygons`
`...`	further arguments passed to or from other methods
`use_region.id`	default `FALSE`, in `identify` use 1-based observation numbers, otherwise use the `nb` `region.id` attribute values

`x`	a numeric vector the same length as the neighbours list in listw
`listw`	a `listw` object created for example by `nb2listw`
`n`	number of zones
`n1`	n - 1
`S0`	global sum of weights
`zero.policy`	default `attr(listw, "zero.policy")` as set when `listw` was created, if attribute not set, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA

`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

`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.

`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.

`coords`	matrix of region point coordinates or SpatialPoints object or `sfc` points object
`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 `nQuadSegs` argument in `geos_unary`
`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 `region.id`, default `seq(1, nrow(x))`
`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 `par(..)`

`grid`	An object of class `GridTopology`. One can avoid to supply this by just suplying the dimentions in the `d` argument.
`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 `nb`. If FALSE, the result is a matrix with `3^k` columns if `self = TRUE` or `3^k-1` if `self = FALSE`. Zeros are used for hyper-cubes at boundaries.
`self`	Default FALSE, to indicate if the hyper-cube neighbour itself should be considered a neighbour.

`nb`	an object of class `nb` with a list of integer vectors containing neighbour region number ids.
`verbose`	default NULL, use global option value; if TRUE prints non-matching pairs
`force`	do not respect a neighbours list `sym` attribute and test anyway
`glist`	list of general weights corresponding to neighbours

Package 'spdep'

Help Index

Aggregate a spatial neighbours object

Description

Usage

Arguments

Value

Note

Author(s)

Examples

Measure distance from plot

Description

Usage

Arguments

Value

Author(s)

See Also

Distance-weighted autocovariate

Description

Usage

Arguments

Value

Note

Author(s)

References

See Also

Examples

Data set with 4 life condition indices of Belo Horizonte region

Description

Format

Examples

Cardinalities for neighbours lists

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Generate neighbours list for grid cells

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Choynowski probability map values

Description

Usage

Arguments

Value

Author(s)

References

See Also

Examples

Columbus OH spatial analysis data set

Description

Usage

Examples

Differences between neighbours lists

Description

Usage

Arguments

Value

Author(s)

Examples

Neighbourhood contiguity by distance

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Drop links in a neighbours list

Description

Usage

`fx`	a factor of the same length as the neighbours and weights objects in listw
`listw`	a `listw` object created for example by `nb2listw`
`nsim`	number of permutations
`zero.policy`	default `attr(listw, "zero.policy")` as set when `listw` was created, if attribute not set, 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), "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 `get.spChkOption()`

`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

`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.

`x`	matrix of point coordinates, an object inheriting from SpatialPoints or an `"sf"` or `"sfc"` object; if the `"sf"` or `"sfc"` object geometries are in geographical coordinates (`sf::st_is_longlat(x) == TRUE` and `sf::sf_use_s2() == TRUE`), s2 will be used to find the neighbours because it uses spatial indexing https://github.com/r-spatial/s2/issues/125 as opposed to the legacy method which uses brute-force
`k`	number of nearest neighbours to be returned; where identical points are present, `k` should be at least as large as the largest count of identical points (if `k` is smaller, an error will occur when `s2` is used)
`longlat`	TRUE if point coordinates are 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; longlat will override `kd_tree`
`use_kd_tree`	logical value, if the dbscan package is available, use for finding k nearest neighbours when longlat is FALSE, 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

`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

`knn`	A knn object returned by `knearneigh`
`row.names`	character vector of region ids to be added to the neighbours list as attribute `region.id`, default `seq(1, nrow(x))`
`sym`	force the output neighbours list to symmetry

`x`	a `listw` object created for example by `nb2listw`
`var`	a numeric vector the same length as the neighbours list in listw
`zero.policy`	default `attr(listw, "zero.policy")` as set when `listw` was created, if attribute not set, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA
`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

`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

`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.

`listw`	a `listw` object from for example `nb2listw`
`sn`	a `spatial.neighbour` object
`style`	default NULL, missing, set to "M" and warning given; if not "M", passed to `nb2listw` to re-build the object
`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 `mat2listw`