diff --git a/DESCRIPTION b/DESCRIPTION
index 71c7c5c7..0b1a5ec2 100644
--- 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
index cca7a696..7d84fd7f 100644
--- 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
index 925c4ec6..52a4be81 100644
--- 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
new file mode 100644
index 00000000..d0fe98a5
--- /dev/null
+++ b/R/licd_boots.R
@@ -0,0 +1,450 @@
+licd_multi <- function(fx, listw, zero.policy=attr(listw, "zero.policy"),
+    adjust.n=TRUE, nsim = 0L, iseed = NULL, no_repeat_in_row=FALSE,
+    control=list()) {
+    con <- list(comp_binary=TRUE, binomial_punif_alternative="greater",
+        jcm_same_punif_alternative="less", jcm_diff_punif_alternative="greater")
+    nmsC <- names(con)
+    con[(namc <- names(control))] <- control
+    if (length(noNms <- namc[!namc %in% nmsC])) 
+        warning("unknown names in control: ", paste(noNms, collapse = ", "))
+     if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
+        "is not a listw object"))
+    if(!is.factor(fx)) stop(paste(deparse(substitute(fx)),
+        "is not a factor"))
+    if (any(is.na(fx))) stop("NA in factor")
+    if (is.null(zero.policy))
+        zero.policy <- get("zeroPolicy", envir = .spdepOptions)
+    stopifnot(is.logical(zero.policy))
+    n <- length(listw$neighbours)
+    if (n != length(fx)) stop("objects of different length")
+    nb <- listw$neighbours
+    crd <- card(nb)
+    if (any(crd == 0L) && !zero.policy)
+        stop("regions with no neighbours found - zero.policy not supported")
+    ifx <- as.integer(fx)
+    k <- length(levels(fx))
+    tifx <- table(ifx)/n
+    if (k < 2) stop("must be at least two levels in factor")
+    wt <- listw$weights
+    almost.1 <- 1 - 64 * .Machine$double.eps
+
+    env <- new.env()
+    assign("crd", card(nb), envir = env) # cardinality
+    assign("nb", nb, envir = env) # nbs
+    assign("lww", wt, envir = env) # weights
+    assign("nsim", nsim, envir=env) # sims
+    assign("fxi", fx, envir = env) # x col
+    assign("xi", ifx, envir = env) # x col
+    assign("tifx", tifx, envir = env) # ifx props
+    assign("n", n, envir=env)
+    assign("k", k, envir=env)
+    assign("lk", levels(fx), envir=env)
+    assign("listw", listw, envir=env)
+    assign("no_repeat_in_row", no_repeat_in_row, envir=env)
+    assign("almost.1", almost.1, envir=env)
+    assign("comp_binary", con$comp_binary, envir=env)
+    varlist = ls(envir = env)
+
+    permLICD_int <- function(i, env) {
+
+        crdi <- get("crd", envir = env)[i] # get the ith cardinality
+        x <- get("xi", envir = env) # get xs
+        fx <- get("fxi", envir = env) # get xs
+        x_i <- x[-i] # remove the observed x value for cond. permutations
+        xi <- x[i]
+        nb_i <- get("nb", envir = env)[[i]] # nbs for ith element
+        w_i <- get("lww", envir = env)[[i]] # weights for ith element
+        nsim <- get("nsim", envir = env) # no. simulations
+        permutation <- nsim > 0
+        n_i <- get("n", envir=env) - 1L
+        tifx <- get("tifx", envir=env)
+        listw <- get("listw", envir=env)
+        almost.1 <- get("almost.1", envir=env)
+        k <- get("k", envir=env)
+        lk <- get("lk", envir=env)
+        comp_binary <- get("comp_binary", envir=env)
+
+        local_jcm_BW <- function(xi, sn, wc) {
+            xi <- factor(as.numeric(!xi)+1, levels=c(1L, 2L),
+                labels=c("1", "2"))
+            tab <- table(xi)
+            kBW <- length(tab)
+            y <- factor(paste(xi[sn[,1]], ":", xi[sn[,2]], sep=""),
+                levels=as.vector(outer(1:kBW, 1:kBW, 
+        	FUN=function(X,Y) paste(X,Y,sep=":"))))
+            res <- matrix(tapply(sn[,3], y, sum), ncol=kBW)/2
+            res[is.na(res)] <- 0
+            if (all(dim(res) == c(1L, 1L))) res <- rbind(cbind(res, 0), c(0, 0))
+            BWlk <- c("B", "W")
+            ldiag <- res[2,1] + res[1,2]
+            ntab <- as.numeric(as.vector(tab))
+            Ejc <- (wc$S0*(ntab*(ntab-1))) / (2*wc$nwcn1)
+            Vjc <- (wc$S1*(ntab*(ntab-1))) / (wc$nwcn1)
+            Vjc <- Vjc + (((wc$S2 - 2*wc$S1)*ntab*(ntab-1)*(ntab-2)) /
+                (wc$nwcn2))
+            Vjc <- Vjc + (((wc$S02 + wc$S1 - wc$S2)*ntab*(ntab-1)*(ntab-2)*
+                (ntab-3)) / (wc$nwcn2*wc$n3))
+            Vjc <- (0.25 * Vjc) - Ejc^2
+            nrns <- ntab[2]*ntab[1]
+            pntab <- (ntab*(ntab-1))
+            nrns1 <- pntab[2]*pntab[1]
+            Exp <- (wc$S0*(nrns)) / (wc$nwcn1)
+            Var <- (2*wc$S1*nrns)/(wc$nwcn1)
+            Var <- Var + (((wc$S2 - 2*wc$S1) * nrns * (ntab[2]+ntab[1]-2))/ 
+                (wc$nwcn2))
+            Var <- Var + ((4*(wc$S02 + wc$S1 - wc$S2)*nrns1) / (wc$nwcn2*wc$n3))
+            Var <- (0.25 * Var) - Exp^2
+            JtotExp <- sum(Exp)
+            O_jc_BW <- c(diag(res), ldiag)
+            E_jc_BW <- c(Ejc, Exp)
+            V_jc_BW <- c(Vjc, Var)
+            jc_conf_chi_BW <- local_chi(O_jc_BW, E_jc_BW)
+            list(jc_conf_chi_BW, O_jc_BW, E_jc_BW, V_jc_BW)
+        }
+        local_chi <- function(O, E) {
+            sum((O - E)^2 / E, na.rm=TRUE)
+        }
+        local_anscombe <- function(s, n) {
+            asin(sqrt((s+(3/8))/(n+(3/4))))
+        }
+
+        if (crdi == 0L) return(rep(NA_real_, 28))
+
+        x_nb_i <- c(xi, x[nb_i])
+        x_nb_i_xi <- x_nb_i == xi
+        w_i_i <- c(0, w_i)
+        if (comp_binary) c1_comp_obs_i <- sum(x_nb_i_xi)
+        else c1_comp_obs_i <- sum(w_i_i * x_nb_i_xi) + 1
+        c2_prop_level_i <- tifx[xi]
+        if (comp_binary) c3_crdip1 <- crdi + 1
+        else c3_crdip1 <- sum(w_i_i) + 1
+        c4_comp_bin_BW_i <- pbinom(c1_comp_obs_i, c3_crdip1,
+            c2_prop_level_i, lower.tail=TRUE)
+        c5_comp_bin_BW_i <- pbinom(c1_comp_obs_i, c3_crdip1,
+            c2_prop_level_i, lower.tail=FALSE)
+        c5a_comp_bin_BW_i <- pbinom(c1_comp_obs_i-1, c3_crdip1,
+            c2_prop_level_i, lower.tail=FALSE)
+        O_BW <- c(c1_comp_obs_i, (c3_crdip1-c1_comp_obs_i))
+        E_BW <- c3_crdip1 * c(c2_prop_level_i, 1-c2_prop_level_i)
+        c6_comp_chi_BW_i <- local_chi(O_BW, E_BW)
+        O_k <- table(factor(x_nb_i, levels=1:k))
+        E_k <- c3_crdip1 * tifx
+        c7_comp_chi_K_i <- local_chi(O_k, E_k)
+        c8_comp_ans_BW_i <- local_anscombe(c1_comp_obs_i, c3_crdip1)
+
+        sub_i <- sort(c(i, nb_i))
+        sub_i <- 1:n %in% sub_i
+        sub_xi <- x[sub_i] == xi
+        i_here <- which(i %in% which(sub_i))
+        listw_subi <- spdep::subset.listw(listw, sub_i)
+ 	sn <- spdep::listw2sn(listw_subi)
+	wc <- spdep::spweights.constants(listw_subi, zero.policy=zero.policy, 
+		adjust.n=adjust.n)
+        wc$S02 <- wc$S0*wc$S0
+        wc$nwcn1 <- wc$n*wc$n1
+        wc$nwcn2 <- wc$nwcn1*wc$n2
+
+        local_jcm_obs <- local_jcm_BW(sub_xi, sn, wc)
+        local_jcm_obs_chi_BW <- local_jcm_obs[[1]]
+        local_jcm_obs_count_BB <- local_jcm_obs[[2]][1]
+        local_jcm_obs_count_WW <- local_jcm_obs[[2]][2]
+        local_jcm_obs_count_BW <- local_jcm_obs[[2]][3]
+        zs <- (local_jcm_obs[[2]] - local_jcm_obs[[3]]) / 
+            sqrt(local_jcm_obs[[4]])
+        local_jcm_obs_z_BB <- zs[1]
+        local_jcm_obs_z_WW <- zs[2]
+        local_jcm_obs_z_BW <- zs[3]
+        local_jcm_all_BB <- all(sub_xi)
+
+        if (permutation) {
+
+            no_repeat_in_row <- get("no_repeat_in_row", envir=env)
+        # create matrix of replicates for composition
+            if (no_repeat_in_row) {
+                samples <- .Call("perm_no_replace", as.integer(nsim),
+                  as.integer(n_i), as.integer(crdi), PACKAGE="spdep")
+                sx_i <- matrix(x_i[samples], ncol=crdi, nrow=nsim)
+            } else {
+                sx_i <- matrix(sample(x_i, size = crdi * nsim, replace = TRUE),
+                    ncol = crdi, nrow = nsim)
+            }
+
+            sx_i <- cbind(rep(xi, times=nsim), sx_i)
+        
+            c1_comp_sim_i <- apply(sx_i, 1,
+               function(y) ifelse(comp_binary, sum(y == xi),
+                   sum(w_i_i * (y == xi)) + 1))
+            c1_comp_sim_i_rank <- rank(c(c1_comp_sim_i,
+                c1_comp_obs_i))[(nsim + 1L)]
+            c4_comp_bin_BW_sim_i <- sapply(c1_comp_sim_i, function(y) {
+                pbinom(y, c3_crdip1, c2_prop_level_i, lower.tail=TRUE)
+            })
+            c4_comp_bin_BW_sim_i_rank <- rank(c(c4_comp_bin_BW_sim_i,
+                almost.1 * c4_comp_bin_BW_i))[(nsim + 1L)]
+            c5_comp_bin_BW_sim_i <- sapply(c1_comp_sim_i, function(y) {
+                pbinom(y, c3_crdip1, c2_prop_level_i, lower.tail=FALSE)
+            })
+            c5_comp_bin_BW_sim_i_rank <- rank(c(c5_comp_bin_BW_sim_i,
+                almost.1 * c5_comp_bin_BW_i))[(nsim + 1L)]
+            c5a_comp_bin_BW_sim_i <- sapply(c1_comp_sim_i, function(y) {
+                pbinom(y-1, c3_crdip1, c2_prop_level_i, lower.tail=FALSE)
+            })
+            c5a_comp_bin_BW_sim_i_rank <- rank(c(c5a_comp_bin_BW_sim_i,
+                almost.1 * c5a_comp_bin_BW_i))[(nsim + 1L)]
+            c6_comp_chi_BW_sim_i <- sapply(c1_comp_sim_i, function(y) {
+                local_chi(c(y, c3_crdip1-y), E_BW)
+            })
+            c6_comp_chi_BW_sim_i_rank <- rank(c(c6_comp_chi_BW_sim_i,
+                almost.1 * c6_comp_chi_BW_i))[(nsim + 1L)]
+            c7_comp_chi_K_sim_i <- apply(sx_i, 1, function(y) {
+                O_k <- table(factor(y, levels=1:k))
+                local_chi(O_k, E_k)
+            })
+            c7_comp_chi_K_sim_i_rank <- rank(c(c7_comp_chi_K_sim_i,
+                almost.1 * c7_comp_chi_K_i))[(nsim + 1L)]
+            c8_comp_ans_BW_sim_i <- sapply(c1_comp_sim_i, function(s) {
+                local_anscombe(s, c3_crdip1)
+            })
+            c8_comp_ans_BW_sim_i_rank <- rank(c(c8_comp_ans_BW_sim_i,
+                almost.1 * c8_comp_ans_BW_i))[(nsim + 1L)]
+
+        # create matrix of replicates for configuration
+            x_nb_iBW <- x[nb_i] == xi
+            if (no_repeat_in_row) {
+                samples <- .Call("perm_no_replace", as.integer(nsim),
+                    as.integer(crdi), as.integer(crdi), PACKAGE="spdep")
+                sx_i_i <- matrix(x_nb_iBW[samples], ncol=crdi, nrow=nsim)
+            } else {
+                sx_i_i <- matrix(sample(x_nb_iBW, size = crdi * nsim,
+                    replace = TRUE), ncol = crdi, nrow = nsim)
+            }
+
+            xi1 <- rep(TRUE, times=nsim)
+            if (i_here == 1L) sx_i_i <- cbind(xi1, sx_i_i)
+            else if (i_here == (crdi+1)) sx_i_i <- cbind(sx_i_i, xi1)
+            else sx_i_i <- cbind(sx_i_i[, 1:(i_here-1)], xi1, sx_i_i[, 
+                i_here:crdi])
+
+            local_jcm_sim <- apply(sx_i_i, 1, function(y) {
+                local_jcm_BW(y, sn, wc)
+            }, simplify=FALSE)
+
+            local_jcm_chi_BW_sim <- sapply(local_jcm_sim, "[[", 1)
+            local_jcm_chi_BW_sim_rank <- rank(c(local_jcm_chi_BW_sim,
+                almost.1 * local_jcm_obs_chi_BW))[(nsim + 1L)]
+
+            zs <- t(sapply(local_jcm_sim, function(y) {
+                (y[[2]] - y[[3]]) / sqrt(y[[4]])
+            }))
+
+            local_jcm_z_BB_sim_rank <- rank(c(zs[,1],
+                almost.1 * local_jcm_obs_z_BB))[(nsim + 1L)]
+            local_jcm_z_WW_sim_rank <- rank(c(zs[,2],
+                almost.1 * local_jcm_obs_z_WW))[(nsim + 1L)]
+            local_jcm_z_BW_sim_rank <- rank(c(zs[,3],
+                almost.1 * local_jcm_obs_z_BW))[(nsim + 1L)]
+        } else {
+            c1_comp_sim_i_rank <- c4_comp_bin_BW_sim_i_rank <-
+            c5_comp_bin_BW_sim_i_rank <- c5a_comp_bin_BW_sim_i_rank <- 
+            c6_comp_chi_BW_sim_i_rank <- c7_comp_chi_K_sim_i_rank <- 
+            c8_comp_ans_BW_sim_i_rank <- local_jcm_chi_BW_sim_rank <- 
+            local_jcm_z_BB_sim_rank <- local_jcm_z_BW_sim_rank <- 
+            local_jcm_z_WW_sim_rank <- NA_real_
+        }
+       
+        res_i <- c(xi, c1_comp_obs_i, unname(c2_prop_level_i), c3_crdip1,
+            c4_comp_bin_BW_i, c5_comp_bin_BW_i, c5a_comp_bin_BW_i,
+            c6_comp_chi_BW_i, c7_comp_chi_K_i, c8_comp_ans_BW_i,
+            local_jcm_obs_chi_BW, local_jcm_obs_count_BB,
+            local_jcm_obs_count_BW, local_jcm_obs_count_WW, 
+            local_jcm_obs_z_BB, local_jcm_obs_z_WW, local_jcm_obs_z_BW,
+            c1_comp_sim_i_rank, c4_comp_bin_BW_sim_i_rank,
+            c5_comp_bin_BW_sim_i_rank, c5a_comp_bin_BW_sim_i_rank, 
+            c6_comp_chi_BW_sim_i_rank, c7_comp_chi_K_sim_i_rank,
+            c8_comp_ans_BW_sim_i_rank, local_jcm_chi_BW_sim_rank,
+            local_jcm_z_BB_sim_rank, local_jcm_z_BW_sim_rank,
+            local_jcm_z_WW_sim_rank, local_jcm_all_BB)
+        res_i
+    }
+    out <- run_perm(fun=permLICD_int, idx=1:n, env=env, iseed=iseed,
+        varlist=varlist)
+
+    local_comp <- data.frame(ID=1:n, category_i=out[,1], count_like_i=out[,2],
+        prop_i=out[,3], count_nbs_i=out[,4], pbinom_like_BW=out[,5],
+        pbinom_unlike_BW=out[,6], pbinom_unlike_BW_alt=out[,7],
+        chi_BW_i=out[,8], chi_K_i=out[,9], anscombe_BW=out[,10])
+
+    pval_jcm_obs_BB <- pnorm(out[,15], lower.tail=FALSE)
+    pval_jcm_obs_WW <- pnorm(out[,16], lower.tail=FALSE)
+    pval_jcm_obs_BW <- pnorm(out[,17], lower.tail=TRUE)
+    sameB <- as.logical(out[,29])
+    if (any(sameB)) {
+        pval_jcm_obs_BB[sameB] <- 0
+        pval_jcm_obs_WW[sameB] <- 1
+        pval_jcm_obs_BW[sameB] <- 1
+    }
+    pval_jcm_obs_BB[is.nan(pval_jcm_obs_BB)] <- 1
+    pval_jcm_obs_WW[is.nan(pval_jcm_obs_WW)] <- 1
+    pval_jcm_obs_BW[is.nan(pval_jcm_obs_BW)] <- 1
+    local_config <- data.frame(ID=1:n, jcm_chi_obs=out[,11],
+        jcm_count_BB_obs=out[,12], jcm_count_BW_obs=out[,13],
+        jcm_count_WW_obs=out[,14], pval_jcm_obs_BB=pval_jcm_obs_BB,
+        pval_jcm_obs_WW=pval_jcm_obs_WW, pval_jcm_obs_BW=pval_jcm_obs_BW)
+    local_comp_sim <- local_config_sim <- NULL
+    if (nsim > 0) {
+        pr_bpnsim <- probs_lut("pbinom_like_BW", nsim,
+            alternative=con$binomial_punif_alternative)
+        local_comp_sim <- data.frame(ID=1:n, rank_sim_count_like_i=out[,18],
+            pbinom_like_BW=pr_bpnsim[out[,19]],
+            pbinom_unlike_BW=pr_bpnsim[out[,20]],
+            pbinom_unlike_BW_alt=pr_bpnsim[out[,21]],
+            rank_sim_chi_BW=out[,22], rank_sim_chi_K=out[,23],
+            rank_sim_anscombe_BW=out[,24])
+        pr_jcmnsim <- probs_lut("jcm_same", nsim,
+            alternative=con$jcm_same_punif_alternative)
+        pr_jcmnsim1 <- probs_lut("jcm_diff", nsim,
+            alternative=con$jcm_diff_punif_alternative)
+        pval_jcm_obs_BB_sim <- pr_jcmnsim[out[,26]]
+        pval_jcm_obs_BW_sim <- pr_jcmnsim[out[,27]]
+        pval_jcm_obs_WW_sim <- pr_jcmnsim1[out[,28]]
+        if (any(sameB)) {
+            pval_jcm_obs_BB_sim[sameB] <- 0
+            pval_jcm_obs_WW_sim[sameB] <- 1
+            pval_jcm_obs_BW_sim[sameB] <- 1
+        }
+        pval_jcm_obs_BB_sim[is.nan(pval_jcm_obs_BB)] <- 1
+        pval_jcm_obs_WW_sim[is.nan(pval_jcm_obs_WW)] <- 1
+        pval_jcm_obs_BW_sim[is.nan(pval_jcm_obs_BW)] <- 1
+        
+        local_config_sim <- data.frame(ID=1:n, jcm_chi_sim_rank=out[,25],
+        pval_jcm_obs_BB_sim=pval_jcm_obs_BB_sim,
+        pval_jcm_obs_BW_sim=pval_jcm_obs_BW_sim,
+        pval_jcm_obs_WW_sim=pval_jcm_obs_WW_sim)
+    }
+
+    colnames(out) <- c("category_i", "count_like_i", "prop_i", "count_nbs_i",
+        "bin_like_BW", "bin_unlike_BW", "bin_unlike_BW_alt", "chi_BW_i",
+        "chi_K_i", "anscombe_BW", "jcm_chi_obs", "jcm_count_BB_obs", 
+        "jcm_count_BW_obs", "jcm_count_WW_obs", "local_jcm_obs_z_BB",
+        "local_jcm_obs_z_WW", "local_jcm_obs_z_BW",
+
+        "rank_sim_count_like_i", "rank_sim_bin_like_BW",
+        "rank_sim_bin_unlike_BW", "rank_sim_bin_unlike_BW_alt",
+        "rank_sim_chi_BW", "rank_sim_chi_K", "rank_sim_anscombe_BW",
+        "jcm_chi_sim_rank", "jcm_z_BB_sim_rank", "jcm_z_BW_sim_rank",
+        "jcm_z_WW_sim_rank", "local_jcm_all_BB")
+    res <- list(local_comp=local_comp, local_config=local_config, local_comp_sim=local_comp_sim, local_config_sim=local_config_sim)
+    attr(res, "out") <- out
+    attr(res, "ncpus") <- attr(out, "ncpus")
+    attr(res, "nsim") <- nsim
+    attr(res, "con") <- con
+    class(res) <- c("licd", class(res))
+    res
+}
+
+hotspot.licd <- function(obj, type="both", cutoff=0.05, p.adjust="none", 
+    droplevels=TRUE, control=list(), ...) {
+    con <- list(binomial_sidak=2, binomial_overlap=TRUE, jcm_sidak=3)
+    nmsC <- names(con)
+    con[(namc <- names(control))] <- control
+    if (length(noNms <- namc[!namc %in% nmsC])) 
+        warning("unknown names in control: ", paste(noNms, collapse = ", "))
+    comp <- config <- FALSE
+    type <- match.arg(type, c("both", "comp", "config"))
+    if (type == "both") comp <- config <- TRUE
+    else if (type == "comp") comp <- TRUE
+    else config = TRUE
+    binom_sc <- 1-(1-cutoff)^(1/con$binomial_sidak)
+    jcm_sc <- 1-(1-cutoff)^(1/con$jcm_sidak)
+    local_comp <- NULL
+    local_comp_sim <- NULL
+    if (comp) {
+        if (!con$binomial_overlap) unlike <- obj$local_comp$pbinom_unlike_BW
+        else unlike <- obj$local_comp$pbinom_unlike_BW_alt
+        like <- obj$local_comp$pbinom_like_BW
+        local_comp <- factor(
+            ifelse(p.adjust(unlike, p.adjust) < binom_sc, "Cluster", 
+            ifelse(p.adjust(like, p.adjust) < binom_sc, "Outlier", "Dispersed")))
+        local_comp_sim <- NULL
+        if (attr(obj, "nsim") > 0) {
+            if (!con$binomial_overlap)
+                unlike <- obj$local_comp_sim$pbinom_unlike_BW
+            else unlike <- obj$local_comp_sim$pbinom_unlike_BW_alt
+            like <- obj$local_comp_sim$pbinom_like_BW
+            local_comp_sim <- factor(
+            ifelse(p.adjust(unlike, p.adjust) < binom_sc, "Cluster", 
+            ifelse(p.adjust(like, p.adjust) < binom_sc, "Outlier", "Dispersed")))
+        }
+    }
+    local_config <- NULL
+    local_config_sim <- NULL
+    if (config) {
+        BB <- obj$local_config$pval_jcm_obs_BB
+        WW <- obj$local_config$pval_jcm_obs_WW
+        BW <- obj$local_config$pval_jcm_obs_BW
+        BB <- p.adjust(BB, p.adjust)
+        WW <- p.adjust(WW, p.adjust)
+        BW <- p.adjust(BW, p.adjust)
+        JC.pvalue <- cbind(BB, WW, BW)
+        crit <- suppressWarnings(apply(JC.pvalue, 1,
+            function(y) min(y, na.rm=TRUE) < jcm_sc))
+        wh13 <- apply(JC.pvalue, 1, function(y) {
+            o <- which.min(y)
+            if (length(o) == 0L) o <- 2L
+            o
+        })
+        local_config <- rep("No cluster", nrow(JC.pvalue))
+        local_config[crit & wh13 == 3] <- "Dispersed"
+        local_config[crit & wh13 == 1] <- "Cluster"
+        local_config <- factor(local_config)
+        local_config_sim <- NULL
+        if (attr(obj, "nsim") > 0) {
+            BB <- obj$local_config_sim$pval_jcm_obs_BB_sim
+            WW <- obj$local_config_sim$pval_jcm_obs_WW_sim
+            BW <- obj$local_config_sim$pval_jcm_obs_BW_sim
+            BB <- p.adjust(BB, p.adjust)
+            WW <- p.adjust(WW, p.adjust)
+            BW <- p.adjust(BW, p.adjust)
+            JC.pvalue <- cbind(BB, WW, BW)
+            crit <- suppressWarnings(apply(JC.pvalue, 1,
+                function(y) min(y, na.rm=TRUE) < jcm_sc))
+            wh13 <- apply(JC.pvalue, 1, function(y) {
+                o <- which.min(y)
+                if (length(o) == 0L) o <- 2L
+                o
+            })
+            local_config_sim <- rep("No cluster", nrow(JC.pvalue))
+            local_config_sim[crit & wh13 == 3] <- "Dispersed"
+            local_config_sim[crit & wh13 == 1] <- "Cluster"
+            local_config_sim <- factor(local_config_sim)
+        }
+    }
+    both <- NULL
+    both_sim <- NULL
+    both_recode <- NULL
+    both_recode_sim <- NULL
+    if (type == "both") {
+        both <- interaction(local_comp, local_config)
+        both_recode <- rep("No cluster", length(local_comp))
+        both_recode[both == "Cluster.Cluster"] <- "Cluster"
+        both_recode[both == "Cluster.No cluster"] <- "Clump"
+        both_recode[both == "Outlier.No cluster"] <- "Outlier"
+        both_recode[both == "Dispersed.Dispersed"] <- "Dispersed"
+        both_recode[both == "Outlier.Dispersed"] <- "Outlier in dispersion area"
+        both_recode <- factor(both_recode)
+        if (attr(obj, "nsim") > 0) {
+            both_sim <- interaction(local_comp_sim, local_config_sim)
+            both_recode_sim <- rep("No cluster", length(local_comp))
+            both_recode_sim[both_sim == "Cluster.Cluster"] <- "Cluster"
+            both_recode_sim[both_sim == "Cluster.No cluster"] <- "Clump"
+            both_recode_sim[both_sim == "Outlier.No cluster"] <- "Outlier"
+            both_recode_sim[both_sim == "Dispersed.Dispersed"] <- "Dispersed"
+            both_recode_sim[both_sim == "Outlier.Dispersed"] <- "Outlier in dispersion area"
+            both_recode_sim <- factor(both_recode_sim)
+        }
+    }
+    list(ID=obj$local_comp$ID, local_comp=local_comp,
+        local_comp_sim=local_comp_sim, local_config=local_config,
+        local_config_sim=local_config_sim, both=both, both_sim=both_sim,
+        both_recode=both_recode, both_recode_sim=both_recode_sim)
+}
diff --git a/R/lisa_perm.R b/R/lisa_perm.R
index bedd48fd..e4e7605a 100644
--- 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
index 687ce90a..74bfcd5d 100644
--- 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
index 54646f67..6fd681b3 100644
--- 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
index 378cb5e7..e25cbdd7 100644
--- 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
index 9bb3963c..9c523333 100644
--- 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
index 1f555c01..6bd7866c 100644
--- 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
index 66a29c54..0a40debf 100644
--- 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
index ce3bf632..327da19b 100644
--- 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
new file mode 100644
index 00000000..30e190cf
--- /dev/null
+++ 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
index daef1608..33843e7c 100644
--- 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
index a7012d9b..260e86ec 100644
--- 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
index f4029387..70699fb4 100644
--- 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
index 70e7ada4..767deec1 100644
--- 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}
 
diff --git a/man/localmoran_bv.Rd b/man/localmoran_bv.Rd
index afab28c4..3b0092c0 100644
--- a/man/localmoran_bv.Rd
+++ b/man/localmoran_bv.Rd
@@ -17,7 +17,7 @@ localmoran_bv(x, y, listw, nsim = 199, scale = TRUE, alternative="two.sided",
 \item{scale}{default \code{TRUE}.}
 
 \item{alternative}{a character string specifying the alternative hypothesis, must be one of "greater" (default), "two.sided", or "less".}
-\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}}
   \item{zero.policy}{default 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}
 }