Add normalized Moran's I statistic and tests

.Rbuildignore

@@ -7,3 +7,5 @@ README.md
 .github
 .git
 _pkgdown.yml
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,7 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+src/*.o
+src/*.so
+src/*.dll

DESCRIPTION

@@ -43,10 +43,33 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"),
 	person("Justyna", "Wilk", role = "ctb",
            comment = c(ORCID = "0000-0003-1495-2910")),
 	person("Levi", "Wolf", role = "ctb"),
-	person("Danlin", "Yu", role = "ctb"))
+	person("Danlin", "Yu", role = "ctb"),
+	person("Giuliani", "Diego", role = "ctb", 
+	         comment = c(ORCID = "0000-0002-7198-6714")))
 Depends: R (>= 3.3.0), methods, spData (>= 2.3.1), sf 
 Imports: stats, deldir, boot (>= 1.3-1), graphics, utils, grDevices, units, s2, e1071, sp (>= 1.0)
-Suggests: spatialreg (>= 1.2-1), Matrix, parallel, dbscan, RColorBrewer, lattice, xtable, foreign, igraph, RSpectra, knitr, classInt, tmap, spam, ggplot2, rmarkdown, tinytest, rgeoda (>= 0.0.11.1), mipfp, Guerry, codingMatrices
+Suggests: 
+    spatialreg (>= 1.2-1),
+    Matrix,
+    parallel,
+    dbscan,
+    RColorBrewer,
+    lattice,
+    xtable,
+    foreign,
+    igraph,
+    RSpectra,
+    knitr,
+    classInt,
+    tmap,
+    spam,
+    ggplot2,
+    rmarkdown,
+    tinytest,
+    rgeoda (>= 0.0.11.1),
+    mipfp,
+    Guerry,
+    codingMatrices
 URL: https://github.com/r-spatial/spdep/, https://r-spatial.github.io/spdep/
 BugReports: https://github.com/r-spatial/spdep/issues/
 Description: A collection of functions to create spatial weights matrix
@@ -95,3 +118,4 @@ Description: A collection of functions to create spatial weights matrix
 License: GPL (>= 2)
 VignetteBuilder: knitr
 RoxygenNote: RoxygenNote: 6.1.1
+Config/testthat/edition: 3

NAMESPACE

@@ -105,6 +105,10 @@ export(plot_spatialcoords, plot_moran, plot_spatialscree, localdelta,
         factorial_coordinates, plot_factorialcoords, plot_factorialscree,
         cornish_fisher)
 
+export(normalized.moran)
+export(normalized.moran.mc)
+export(normalized.moran.test)
+
 S3method(print, nb)
 S3method(summary, nb)
 S3method(print, summary.nb)

NEWS.md

@@ -4,6 +4,8 @@
 
 * add `NO_LEGACY_NONAP`
 
+* add `normalized.moran.test` and `normalized.moran.mc` to compute the I_6 normalized Moran's I index and perform the associated tests by Tillé et al. (2026)
+
 # Version 1.4-2 (2026-02-13)
 
 * robustness checks on `licd_multi`, update `hotspot` method for `licd_multi` output

R/normalized.moran.R

@@ -0,0 +1,22 @@
+normalized.moran <- function (x, listw, zero.policy = attr(listw, "zero.policy"),
+                              NAOK = FALSE)
+{
+  if (is.null(zero.policy))
+    zero.policy <- get.ZeroPolicyOption()
+  stopifnot(is.logical(zero.policy))
+  n1 <- length(listw$neighbours)
+  x <- c(x)
+  if (n1 != length(x))
+    stop("objects of different length")
+  xx <- mean(x, na.rm = NAOK)
+  z <- x - xx
+  lz <- spdep::lag.listw(listw, z, zero.policy = zero.policy, NAOK = NAOK)
+  if(NAOK) {
+    z_lz <- cbind(z,lz)
+    z_lz <- na.omit(z_lz)
+    z <- z_lz[,1]
+    lz <- z_lz[,2]
+  }
+  I <- cor(z, lz)
+  I
+}

R/normalized.moran.mc.R

@@ -0,0 +1,101 @@
+normalized.moran.mc <- function(x, listw, nsim, zero.policy = attr(listw, "zero.policy"),
+                                alternative = "greater", na.action = na.fail, spChk = NULL,
+                                return_boot = FALSE, adjust.n = TRUE)
+{
+  alternative <- match.arg(alternative, c("greater", "less",
+                                          "two.sided"))
+  wname <- deparse(substitute(listw))
+  if (!inherits(listw, "listw"))
+    stop(wname, "is not a listw object")
+  xname <- deparse(substitute(x))
+  if (!is.numeric(x))
+    stop(xname, "is not a numeric vector")
+  if (is.null(attr(listw$weights, "W")))
+    message(wname, " is not a listw object of style W, then it has been row-standardized")
+  if (is.null(zero.policy))
+    zero.policy <- get.ZeroPolicyOption()
+  stopifnot(is.logical(zero.policy))
+  stopifnot(length(zero.policy) == 1L)
+  if (missing(nsim))
+    stop("nsim must be given")
+  if (is.null(spChk))
+    spChk <- get.spChkOption()
+  stopifnot(is.logical(spChk))
+  stopifnot(length(spChk) == 1L)
+  if (spChk && !chkIDs(x, listw))
+    stop("Check of data and weights ID integrity failed")
+  cards <- card(listw$neighbours)
+  if (!zero.policy && any(cards == 0))
+    stop("regions with no neighbours found")
+  stopifnot(length(na.action) == 1L)
+  if (deparse(substitute(na.action)) == "na.pass")
+    stop("na.pass not permitted")
+  x <- na.action(x)
+  na.act <- attr(x, "na.action")
+  if (!is.null(na.act)) {
+    subset <- !(1:length(listw$neighbours) %in% na.act)
+    listw <- subset(listw, subset, zero.policy = zero.policy)
+    if (return_boot)
+      message("NA observations omitted: ", paste(na.act,
+                                                 collapse = ", "))
+  }
+  n <- length(listw$neighbours)
+  if (n != length(x))
+    stop("objects of different length")
+  gamres <- suppressWarnings(nsim > gamma(n + 1))
+  if (gamres)
+    stop("nsim too large for this number of observations")
+  if (nsim < 1)
+    stop("nsim too small")
+  if (adjust.n)
+    n <- n - sum(cards == 0L)
+  if (is.null(attr(listw$weights, "W"))) {
+    W <- spdep::listw2mat(listw)
+    listw <- spdep::mat2listw(W, style="W", zero.policy = zero.policy)
+  }
+  if (return_boot) {
+    moran_boot <- function(var, i, ...) {
+      var <- var[i]
+      return(normalized.moran(x = var, ...))
+    }
+    p_setup <- spdep:::parallel_setup(NULL)
+    parallel <- p_setup$parallel
+    ncpus <- p_setup$ncpus
+    cl <- p_setup$cl
+    res <- boot::boot(x, statistic = moran_boot, R = nsim, sim = "permutation",
+                      listw = listw, zero.policy = zero.policy,
+                      parallel = parallel, ncpus = ncpus, cl = cl)
+    return(res)
+  }
+  res <- numeric(length = nsim + 1)
+  for (i in 1:nsim) res[i] <- normalized.moran(sample(x), listw,
+                                               zero.policy)
+  res[nsim + 1] <- normalized.moran(x, listw, zero.policy)
+  rankres <- rank(res)
+  xrank <- rankres[length(res)]
+  diff <- nsim - xrank
+  diff <- ifelse(diff > 0, diff, 0)
+  if (alternative == "less")
+    pval <- punif((diff + 1)/(nsim + 1), lower.tail = FALSE)
+  else if (alternative == "greater")
+    pval <- punif((diff + 1)/(nsim + 1))
+  else pval <- punif(abs(xrank - (nsim + 1)/2)/(nsim + 1),
+                     0, 0.5, lower.tail = FALSE)
+  if (!is.finite(pval) || pval < 0 || pval > 1)
+    warning("Out-of-range p-value: reconsider test arguments")
+  statistic <- res[nsim + 1]
+  names(statistic) <- "statistic"
+  parameter <- xrank
+  names(parameter) <- "observed rank"
+  method <- "Monte-Carlo simulation of Normalized Moran I"
+  data.name <- paste(xname, "\nweights:", wname, ifelse(is.null(na.act),
+                                                        "", paste("\nomitted:", paste(na.act, collapse = ", "))),
+                     "\nnumber of simulations + 1:", nsim + 1, "\n")
+  lres <- list(statistic = statistic, parameter = parameter,
+               p.value = pval, alternative = alternative, method = method,
+               data.name = data.name, res = res)
+  if (!is.null(na.act))
+    attr(lres, "na.action") <- na.act
+  class(lres) <- c("htest", "mc.sim")
+  lres
+}

R/normalized.moran.test.R

@@ -0,0 +1,91 @@
+normalized.moran.test <- function (x, listw, zero.policy = attr(listw, "zero.policy"), alternative = "greater",
+                                   na.action = na.fail,
+                                   spChk = NULL, adjust.n = TRUE)
+{
+  alternative <- match.arg(alternative, c("greater", "less",
+                                          "two.sided"))
+  wname <- deparse(substitute(listw))
+  if (!inherits(listw, "listw"))
+    stop(wname, "is not a listw object")
+  xname <- deparse(substitute(x))
+  if (!is.numeric(x))
+    stop(xname, " is not a numeric vector")
+  if (is.null(attr(listw$weights, "W")))
+    message(wname, " is not a listw object of style W, then it has been row-standardized")
+  if (is.null(zero.policy))
+    zero.policy <- get.ZeroPolicyOption()
+  stopifnot(is.logical(zero.policy))
+  stopifnot(length(zero.policy) == 1L)
+  if (is.null(spChk))
+    spChk <- get.spChkOption()
+  stopifnot(is.logical(spChk))
+  stopifnot(length(spChk) == 1L)
+  if (spChk && !chkIDs(x, listw))
+    stop("Check of data and weights ID integrity failed")
+  stopifnot(length(na.action) == 1L)
+  NAOK <- deparse(substitute(na.action)) == "na.pass"
+  x <- na.action(x)
+  na.act <- attr(x, "na.action")
+  if (!is.null(na.act)) {
+    subset <- !(1:length(listw$neighbours) %in% na.act)
+    listw <- subset(listw, subset, zero.policy = zero.policy)
+  }
+  n <- length(listw$neighbours)
+  if (n != length(x))
+    stop("objects of different length")
+  W <- listw2mat(listw)
+  if (is.null(attr(listw$weights, "W")))
+    listw <- spdep::mat2listw(W, style="W", zero.policy = zero.policy)
+  I <- normalized.moran(x, listw, zero.policy = zero.policy,
+                        NAOK = NAOK)
+  one <- rep(1,n)
+  In <- diag(rep(1,n))
+  P <- In - one %*% t(one)/n
+  D <- diag(rowSums(W))
+  invD <- solve(D)
+  B <- P %*% invD %*% W
+  C <- t(W) %*% invD %*% P %*% invD %*% W
+  B_tilde <- (B + t(B)) / 2
+
+  tr <- function(M) sum(diag(M))
+  trP <- tr(P)
+  trC <- tr(C)
+  trB <- tr(B)
+  trB2_tilde <- tr(B_tilde %*% B_tilde)
+  trP2 <- tr(P %*% P)
+  trC2 <- tr(C %*% C)
+  trBP <- tr(B_tilde %*% P)
+  trBC <- tr(B_tilde %*% C)
+  trPC <- tr(P %*% C)
+
+  EI <- -1/sqrt((n-1)*trC)
+  VI <- (2/(trP*trC)) * (trB2_tilde + trB^2/(4*trP^2)*trP2 + trB^2/(4*trC^2)*trC2 -
+                           trB/trP*trBP - trB/trC*trBC +
+                           trB^2/(2*trP*trC)*trPC)
+  if (VI < 0)
+    warning("Negative variance,\ndistribution of variable does not meet test assumptions")
+  ZI <- (I - EI)/sqrt(VI)
+  statistic <- ZI
+  names(statistic) <- "Normalized Moran I statistic standard deviate"
+  if (alternative == "two.sided")
+    PrI <- 2 * pnorm(abs(ZI), lower.tail = FALSE)
+  else if (alternative == "greater")
+    PrI <- pnorm(ZI, lower.tail = FALSE)
+  else PrI <- pnorm(ZI)
+  if (!is.finite(PrI) || PrI < 0 || PrI > 1)
+    warning("Out-of-range p-value: reconsider test arguments")
+  vec <- c(I, EI, VI)
+  names(vec) <- c("Normalized Moran I statistic", "Expectation", "Variance")
+  method <- "Normalized Moran I test under normality"
+  data.name <- paste(xname, "\nweights:", wname, ifelse(is.null(na.act), "",
+                                                        paste("\nomitted:", paste(na.act, collapse = ", "))),
+                     ifelse(adjust.n && isTRUE(any(sum(card(listw$neighbours) ==
+                                                         0L))), "\nn reduced by no-neighbour observations",
+                            ""), "", "\n")
+  res <- list(statistic = statistic, p.value = PrI, estimate = vec,
+              alternative = alternative, method = method, data.name = data.name)
+  if (!is.null(na.act))
+    attr(res, "na.action") <- na.act
+  class(res) <- "htest"
+  res
+}

inst/tinytest/test_normalized.moran_NA.R

@@ -0,0 +1,7 @@
+library(spdep)
+data(oldcol)
+lw <- nb2listw(COL.nb, style="W")
+is.na(COL.OLD$CRIME) <- 1
+expect_error(o <- normalized.moran.test(COL.OLD$CRIME, lw))
+expect_warning(o <- normalized.moran.test(COL.OLD$CRIME, lw, na.action=na.pass))
+expect_silent(o <- normalized.moran.test(COL.OLD$CRIME, lw, na.action=na.omit))

man/normalized.moran.Rd

@@ -0,0 +1,46 @@
+\name{normalized.moran}
+\alias{normalized.moran}
+\title{Compute Normalized Moran's \eqn{I}}
+\usage{
+normalized.moran(
+  x,
+  listw,
+  zero.policy = attr(listw, "zero.policy"),
+  NAOK = FALSE
+)
+}
+\arguments{
+\item{x}{a numeric vector the same length as the neighbours list in listw}
+
+\item{listw}{a \code{listw} object of style 'W' created for example by \code{nb2listw}}
+
+\item{zero.policy}{default \code{attr(listw, "zero.policy")} as set when \code{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}
+
+\item{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.}
+}
+\value{
+The normalized "\eqn{I_6}" Moran's \eqn{I} of the values in \code{x} is computed, as a numeric vector of length one.
+}
+\description{
+A simple function to compute the normalized "\eqn{I_6}" Moran's \eqn{I} statistic proposed by Tillé et al. (2026), called by \code{normalized.moran.test} and \code{normalized.moran.mc}.
+}
+\section{Author}{
+ Diego Giuliani \email{diego.giuliani@unitn.it}
+}
+
+\section{References}{
+ Tillé Y., Giuliani D., Dickson M.M., Espa G. (2026) Revisiting the Normalization of the Moran's \eqn{I} index: A Correlation-based Approach with Inference, \emph{International Statistical Review}
+}
+
+\examples{
+library(spdep)
+data(oldcol)
+col.W <- nb2listw(COL.nb, style="W")
+crime <- COL.OLD$CRIME
+normalized.moran(crime, col.W)
+is.na(crime) <- sample(1:length(crime), 10)
+normalized.moran(crime, col.W, NAOK=TRUE)
+}
+\seealso{
+\link{normalized.moran.test}, \link{normalized.moran.mc}
+}