Merge pull request #167 from r-spatial/boots_licd

Boots licd
r-spatial · Aug 30, 2024 · 393a6c3 · 393a6c3
2 parents 277b48a + 5177dc8
commit 393a6c3
Show file tree

Hide file tree

Showing 18 changed files with 568 additions and 12 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -21,6 +21,8 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"),
            comment = c(ORCID = "0000-0002-9415-4582")),
 	person("Virgilio", "Gómez-Rubio", role = "ctb"),
 	person("Malabika", "Koley", role = "ctb"),
+	person("Tomasz", "Kossowski", role = "ctb",
+           comment = c(ORCID = "0000-0002-9976-4398")),
 	person("Elias", "Krainski", role = "ctb"),
 	person("Pierre", "Legendre", role = "ctb"),
 	person("Nicholas", "Lewin-Koh", role = "ctb"),
@@ -31,11 +33,15 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"),
 	person("Josiah", "Parry", role = "ctb",
            comment = c(ORCID = "0000-0001-9910-865X")),
 	person("Pedro", "Peres-Neto", role = "ctb"),
+	person("Michał", "Pietrzak", role = "ctb",
+           comment = c(ORCID = 0000-0002-9263-4478)),
 	person("Gianfranco", "Piras", role = "ctb"),
 	person("Markus", "Reder", role = "ctb"),
 	person("Jeff", "Sauer", role = "ctb"),
 	person("Michael", "Tiefelsdorf", role = "ctb"),
 	person("René", "Westerholt", role="ctb"),
+	person("Justyna", "Wilk", role = "ctb",
+           comment = c(ORCID = "0000-0003-1495-2910")),
 	person("Levi", "Wolf", role = "ctb"),
 	person("Danlin", "Yu", role = "ctb"))
 Depends: R (>= 3.3.0), methods, spData (>= 2.3.1), sf 

diff --git a/NAMESPACE b/NAMESPACE
@@ -11,7 +11,8 @@ importFrom(stats, influence.measures, lag, punif, lm, var, integrate,
              aggregate, "coefficients", "fitted", "gaussian",
              "model.frame", "model.matrix", "model.response", "na.fail",
              "na.omit", "naresid", "optimise", "printCoefmat", "resid",
-             "terms", "uniroot", "weighted.residuals", "weights", "median")
+             "terms", "uniroot", "weighted.residuals", "weights", "median",
+             "pbinom")
 
 importFrom(deldir, deldir)
 importFrom(boot, boot)
@@ -92,6 +93,7 @@ export(set.mcOption, get.mcOption, set.coresOption, get.coresOption,
 
 export(localmoran_bv, moran_bv,local_joincount_uni, local_joincount_bv)
 
+export(licd_multi)
 
 S3method(print, nb)
 S3method(summary, nb)
@@ -151,3 +153,4 @@ S3method(hotspot, summary.localmoransad)
 S3method(hotspot, data.frame.localmoranex)
 S3method(hotspot, localC)
 S3method(hotspot, localG)
+S3method(hotspot, licd)
diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,7 @@
 # Version 1.3-6 (development)
 
+* adding prototype of LICD ESDA function `licd_multi` and `hotspot` method
+
 * #162 add option for no-neighbour checking for `poly2nb` - default report whether no-neighbour observations are present
 
 * #162 change the default `snap=` argument to `poly2nb` to 10mm 

diff --git a/R/licd_boots.R b/R/licd_boots.R
diff --git a/R/lisa_perm.R b/R/lisa_perm.R
@@ -63,6 +63,7 @@ run_perm <- function(fun, idx, env, iseed, varlist) {
         oo <- lapply(idx, function(i) fun(i, env))
         out <- do.call("rbind", oo)
     }
+    attr(out, "ncpus") <- ncpus
     out
 }
 
@@ -205,6 +206,7 @@ localmoran_perm <- function(x, listw, nsim=499L, zero.policy=attr(listw, "zero.p
     }
 
     out <- run_perm(fun=permI_int, idx=1:n, env=env, iseed=iseed, varlist=varlist)
+    ncpus <- attr(out, "ncpus")
 
     if (sample_Ei) res[,2] <- out[,1]
     else  res[,2] <- EIc
@@ -213,7 +215,9 @@ localmoran_perm <- function(x, listw, nsim=499L, zero.policy=attr(listw, "zero.p
     if (alternative == "two.sided") 
         res[,5] <- 2 * pnorm(abs(res[,4]), lower.tail=FALSE)
     else if (alternative == "greater") 
-        res[,5] <- pnorm(res[,4], lower.tail=FALSE)
+        res[,5] <- pnorm(res[,4
+
+], lower.tail=FALSE)
     else res[,5] <- pnorm(res[,4])
 # look-up table
     probs <- probs_lut(stat="I", nsim=nsim, alternative=alternative)
@@ -241,6 +245,7 @@ localmoran_perm <- function(x, listw, nsim=499L, zero.policy=attr(listw, "zero.p
     if (!is.null(na.act)) attr(res, "na.action") <- na.act
     attr(res, "quadr") <- data.frame(mean=quadr, median=quadr_med, 
         pysal=quadr_ps)
+    attr(res, "ncpus") <- ncpus
     class(res) <- c("localmoran", class(res))
     res
 }
@@ -398,6 +403,7 @@ localG_perm <- function(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy
     }
 
     out <- run_perm(fun=thisfun, idx=1:n, env=env, iseed=iseed, varlist=varlist)
+    ncpus <- attr(out, "ncpus")
 
 # add simulated standard deviate direct output
     res_sim <- (G - out[,1])
@@ -427,6 +433,7 @@ localG_perm <- function(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy
     attr(res, "cluster") <- cut(x, c(-Inf, mean(x), Inf), labels = c("Low", "High"))
     attr(res, "gstari") <- gstari
     attr(res, "call") <- match.call()
+    attr(res, "ncpus") <- ncpus
     class(res) <- "localG"
     res
 }

diff --git a/R/local-joincount-univariate.R b/R/local-joincount-univariate.R
@@ -106,12 +106,14 @@ local_joincount_uni <- function(fx, chosen, listw,
   # commenting out because this should be handled in `probs_lut()`
   # if (alternative == "two.sided") probs <- probs / 2
   p_ranks <- run_perm(permBB_int, index, env, iseed, varlist)
+  ncpus <- attr(p_ranks, "ncpus")
 
   p_res <- rep(NA_real_, length(x))
   p_res[index] <- probs[floor(p_ranks)]
 
   res <- data.frame(obs, p_res)
   colnames(res) <- c("BB", attr(probs, "Prname"))
+  attr(res, "ncpus") <- ncpus
   res
 }
 

diff --git a/R/local-moran-bv.R b/R/local-moran-bv.R
@@ -101,6 +101,7 @@ localmoran_bv <- function(x, y, listw, nsim = 199, scale = TRUE,
   }
 
   out <- run_perm(fun=permI_bv_int, idx=1:n, env=env, iseed=iseed, varlist=varlist)
+  ncpus <- attr(out, "ncpus")
 
   res <- matrix(nrow=n, ncol=7)
   res[,1] <- obs
@@ -128,6 +129,7 @@ localmoran_bv <- function(x, y, listw, nsim = 199, scale = TRUE,
   attr(res, "quadr") <- data.frame(mean=quadr, median=quadr_med, 
     pysal=quadr_ps)
   class(res) <- c("localmoran", class(res))
+  attr(res, "ncpus") <- ncpus
   res
 }
 
diff --git a/R/localC.R b/R/localC.R
@@ -122,6 +122,7 @@ localC_perm.default <- function(x, listw, nsim = 499, alternative = "two.sided",
     reps <- localC_perm_calc(x, listw, obs, nsim, alternative=alternative,
       zero.policy=zero.policy, iseed=iseed, no_repeat_in_row=no_repeat_in_row)
   }
+
   if (ncol(xorig) > 1L) {
     cluster <- rep(3L, length(obs))
     cluster[obs <= reps[, 1]] <- 1L
@@ -290,6 +291,7 @@ localC_perm_calc <- function(x, listw, obs, nsim, alternative="two.sided",
     }
 
     res <- run_perm(fun=permC_int, idx=1:n, env=env, iseed=iseed, varlist=varlist)
+    ncpus <- attr(res, "ncpus")
 
     res[,3] <- (obs - res[,1])/sqrt(res[,2])
     if (alternative == "two.sided")
@@ -305,6 +307,7 @@ localC_perm_calc <- function(x, listw, obs, nsim, alternative="two.sided",
 
     colnames(res) <- c("E.Ci", "Var.Ci", "Z.Ci", Prname,
       paste0(Prname, " Sim"), "Pr(folded) Sim", "Skewness", "Kurtosis")
+    attr(res, "ncpus") <- ncpus
     res
 }
 

diff --git a/R/subset.nb.R b/R/subset.nb.R
@@ -64,6 +64,8 @@ subset.listw <- function(x, subset, zero.policy=attr(x, "zero.policy"), ...) {
     n <- length(nb)
     if (n != length(subset))
 	stop("neighbours list and subset vector different lengths")
+    if (!is.null(attr(x, "region.id"))) 
+        attr(nb, "region.id") <- attr(x, "region.id")
     subnb <- subset.nb(x=nb, subset=subset)
     if (any(card(subnb) == 0L)) {
         if (!zero.policy) {

diff --git a/R/utils.R b/R/utils.R
@@ -199,6 +199,7 @@ listw2star <- function(listw, ireg, style, n, D, a, zero.policy=attr(listw, "zer
             nb[[jj]] <- ireg
             wts[[jj]] <- iwts[j]
 	}
+        attributes(wts) <- attributes(listw$weights)
     }
     res <- list(style=style, neighbours=nb, weights=wts)
     class(res) <- c("listw", "star")

diff --git a/inst/tinytest/test_lisa_perm.R b/inst/tinytest/test_lisa_perm.R
@@ -20,7 +20,7 @@ expect_silent(no <- system.time(localC_perm(xx, lw, nsim=nsim, iseed=iseed))["el
 expect_silent(no <- system.time(localmoran_bv(x, y, lw, nsim=nsim, iseed=iseed))["elapsed"])
 if (require(parallel, quietly=TRUE)) {
  coresOpt <- get.coresOption()
- nc <- detectCores(logical=FALSE)-1L
+ nc <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L
  nc
  mcOpt <- get.mcOption()
 # set nc to 1L here

diff --git a/man/hotspotmap.Rd b/man/hotspotmap.Rd
@@ -6,8 +6,9 @@
 \alias{hotspot.data.frame.localmoranex}
 \alias{hotspot.localG}
 \alias{hotspot.localC}
+\alias{hotspot.licd}
 
-\title{Cluster classifications for local indicators of spatial association}
+\title{Cluster Classifications for Local Indicators of Spatial Association and Local Indicators for Categorical Data}
 \usage{
 hotspot(obj, ...)
 
@@ -23,6 +24,8 @@ hotspot(obj, ...)
 \method{hotspot}{localG}(obj, Prname, cutoff=0.005, p.adjust="fdr", droplevels=TRUE, ...)
 
 \method{hotspot}{localC}(obj, Prname, cutoff=0.005, p.adjust="fdr", droplevels=TRUE, ...)
+\method{hotspot}{licd}(obj, type = "both", cutoff = 0.05, p.adjust = "none", 
+ droplevels = TRUE, control = list(), ...)
 }
 \arguments{
 \item{obj}{An object of class \code{localmoran}, \code{localC} or \code{localG}}
@@ -37,18 +40,24 @@ hotspot(obj, ...)
 
 \item{quadrant.type}{Default \code{"mean"}, for \code{"localmoran"} objects only, can be \code{c("mean", "median", "pysal")} to partition the Moran scatterplot; \code{"mean"} partitions on the means of the variable and its spatial lag, \code{"median"} on medians of the variable and its spatial lag, \code{"pysal"} at zero for the centred variable and its spatial lag}
 
+\item{type}{When \code{obj} is of class \code{licd}, default \code{both}, may also be \code{comp} for local composition or \code{config} for local configuration}
+
+\item{control}{When \code{obj} is of class \code{licd}, default \code{binomial_sidak} 2, \code{binomial_overlap} TRUE, \code{jcm_sidak} 3. \code{binomial_overlap}  may be set FALSE to avoid the Binomial probability values summing to more than unity - the tests in Boots (2003, p. 141) do overlap (\code{>=} and \code{<=}), and the Šidák exponents may be set to 1 to prevent by-observation correction for 2 Binomial and 3 Normal probability values per observation}
+
 \item{...}{other arguments passed to methods.}
 
 }
 \description{
- 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.
+ 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 \code{obj} is of class \code{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.
 
 }
 
 \value{
- A factor showing so-called cluster classification for local indicators of spatial association.
+ A factor showing so-called cluster classification for local indicators of spatial association. When \code{obj} is of class \code{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.
 }
 
+\seealso{\code{\link{licd_multi}}}
+
 \examples{
 orig <- spData::africa.rook.nb
 listw <- nb2listw(orig)

diff --git a/man/licd_multi.Rd b/man/licd_multi.Rd
@@ -0,0 +1,69 @@
+\name{licd_multi}
+\alias{licd_multi}
+
+\title{Local Indicators for Categorical Data}
+
+\description{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}
+
+\usage{
+licd_multi(fx, listw, zero.policy = attr(listw, "zero.policy"), adjust.n = TRUE,
+ nsim = 0L, iseed = NULL, no_repeat_in_row = FALSE, control = list())
+}
+
+\arguments{
+  \item{fx}{a factor with two or more categories, of the same length as the neighbours and weights objects in listw}
+  \item{listw}{a \code{listw} object 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{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}
+  \item{nsim}{default 0, number of conditonal permutation simulations}
+  \item{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}
+  \item{no_repeat_in_row}{default \code{FALSE}, if \code{TRUE}, sample conditionally in each row without replacements to avoid duplicate values, \url{https://github.com/r-spatial/spdep/issues/124}}
+  \item{control}{comp_binary=TRUE, binomial_punif_alternative="greater",
+        jcm_same_punif_alternative="less", jcm_diff_punif_alternative="greater"}
+}
+\details{The original code may be found at \doi{10.5281/zenodo.4283766}}
+
+\value{
+  \item{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}
+  \item{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}
+  \item{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}
+  \item{local_config_sim}{data.frame object with permutation-based LICD local configuration columns: ID, jcm_chi_sim_rank, pval_jcm_obs_BB_sim, pval_jcm_obs_BW_sim, pval_jcm_obs_WW_sim}
+}
+\references{
+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}
+}
+\author{Roger Bivand \email{Roger.Bivand@nhh.no} based on earlier code by Tomasz M. Kossowski, Justyna Wilk and Michał B. Pietrzak}
+
+\note{In order to increase the numbers of neighbours using \code{\link{nblag}} and \code{\link{nblag_cumul}} is advisable; use of binary weights is advised and are in any case used for the composition measure}
+
+\seealso{\code{\link{joincount.multi}}}
+
+\examples{
+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"])
+}
+
+\keyword{spatial}
diff --git a/man/localC.Rd b/man/localC.Rd
@@ -45,7 +45,7 @@ localC_perm(x, ..., zero.policy=NULL, iseed=NULL, no_repeat_in_row=FALSE)
 \item{alternative}{A character defining the alternative hypothesis. Must be one of \code{"two.sided"}, \code{"less"} or \code{"greater"}.}
 \item{...}{other arguments passed to methods.}
 \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{iseed}{default NULL, used to set the seed for possible parallel RNGs}
+\item{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}
 \item{no_repeat_in_row}{default \code{FALSE}, if \code{TRUE}, sample conditionally in each row without replacements to avoid duplicate values, \url{https://github.com/r-spatial/spdep/issues/124}}
 }
 \description{

diff --git a/man/localG.Rd b/man/localG.Rd
@@ -28,7 +28,7 @@ localG_perm(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy"), spChk=NU
   \item{nsim}{default 499, number of conditonal permutation simulations}
   \item{alternative}{a character string specifying the alternative hypothesis, must be one of \code{"two.sided"} (default), \code{"greater"} or \code{"less"}.}
   \item{return_internals}{default \code{TRUE}, unused}
-  \item{iseed}{default NULL, used to set the seed for possible parallel RNGs}
+  \item{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}
   \item{fix_i_in_Gstar_permutations}{default \code{TRUE} (fix x at self in permutations for local G-star), set \code{FALSE} to use pre-1.2-8 behaviour}
   \item{no_repeat_in_row}{default \code{FALSE}, if \code{TRUE}, sample conditionally in each row without replacements to avoid duplicate values, \url{https://github.com/r-spatial/spdep/issues/124}}
 }

diff --git a/man/local_joincount_uni.Rd b/man/local_joincount_uni.Rd
@@ -23,7 +23,7 @@ local_joincount_uni(
 
 \item{nsim}{the number of conditional permutation simulations}
 
-\item{iseed}{default NULL, used to set the seed for possible parallel RNGs}
+\item{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}
 
 \item{no_repeat_in_row}{default \code{FALSE}, if \code{TRUE}, sample conditionally in each row without replacements to avoid duplicate values, \url{https://github.com/r-spatial/spdep/issues/124}}
 

diff --git a/man/localmoran.Rd b/man/localmoran.Rd
@@ -31,7 +31,7 @@ localmoran_perm(x, listw, nsim=499, zero.policy=attr(listw, "zero.policy"),
   \item{adjust.x}{default FALSE, if TRUE, x values of observations with no neighbours are omitted in the mean of x}
   \item{nsim}{default 499, number of conditonal permutation simulations}
   \item{sample_Ei}{default TRUE; if conditional permutation, use the sample $E_i$ values, or the analytical values, leaving only variances calculated by simulation.}
-  \item{iseed}{default NULL, used to set the seed for possible parallel RNGs}
+  \item{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}
   \item{no_repeat_in_row}{default \code{FALSE}, if \code{TRUE}, sample conditionally in each row without replacements to avoid duplicate values, \url{https://github.com/r-spatial/spdep/issues/124}}
 }
 
@@ -60,7 +60,7 @@ statistics: an overview. In P. Longley and M. Batty (eds) \emph{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{doi:10.1111/gean.12304}
+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}