R/localC.R

In spdep: Spatial Dependence: Weighting Schemes, Statistics

localC <- function(x, ..., zero.policy=NULL) {
  UseMethod("localC")
}


localC.default <- function(x, listw, ..., zero.policy=attr(listw, "zero.policy")) {
  # check listw object
  if (!inherits(listw, "listw"))
    stop(paste(deparse(substitute(listw)), "is not a listw object"))

  # check missing values
  if (any(is.na(x))) stop(paste("NA in ", deparse(substitute(x))))

  localC_calc(scale(x), listw, zero.policy=zero.policy)
}

localC.formula <- function(formula, data, listw, ..., zero.policy=attr(listw, "zero.policy")) {
  # check listw object
  if (!inherits(listw, "listw"))
    stop(paste(deparse(substitute(listw)), "is not a listw object."))

  # check if sf object in data
  if (inherits(data, "sf")) {
    data[[attr(data, "sf_column")]] <- NULL
    data <- as.data.frame(data)
  }

  if (any(sapply(data, class) == "list"))
    stop("`data` cannot contain list columns or elements.")

  form_terms <- terms(formula, data = data)

  if (!attr(form_terms, "response") == 0)
    stop("`formula` must be one-sided with no response variable.")

  df <- model.frame(formula, data = data)

  char_cols <- colnames(df)[sapply(df, class) == "character"]

  if (length(char_cols) > 0)
    stop(paste("Formula contains character vectors:", char_cols))

  rowSums(apply(scale(df), 2, localC_calc, listw, zero.policy=zero.policy)) / ncol(df)

}

localC.list <- function(x, listw, ..., zero.policy=attr(listw, "zero.policy")) {

  if (!inherits(listw, "listw"))
    stop(paste(deparse(substitute(listw)), "is not a listw object,"))

  if (!length(unique(lengths(x))) == 1) {
    stop("Elements of x must be of equal length.")
  }

  x <- scale(Reduce(cbind, x))
  rowSums(apply(x, 2, localC_calc, listw, zero.policy=zero.policy)) / ncol(x)

}


localC.matrix <- function(x, listw, ..., zero.policy=attr(listw, "zero.policy")) {

  if (!inherits(listw, "listw"))
    stop(paste(deparse(substitute(listw)), "is not a listw object"))

  if (inherits(x, "character")) stop("x must be a numeric matrix.")

  rowSums(apply(scale(x), 2, localC_calc, listw, zero.policy=zero.policy)) / ncol(x)
}

localC.data.frame <- function(x, listw, ..., zero.policy=attr(listw, "zero.policy")) {

  if (inherits(x, "sf")) {
    x[[attr(x, "sf_column")]] <- NULL
    x <- as.data.frame(x)
  }

  if (!all(sapply(x, class) %in% c("numeric", "integer")))
    stop("Columns of x must be numeric.")

  rowSums(apply(scale(x), 2, localC_calc, listw, zero.policy=zero.policy)) / ncol(x)

}


localC_perm <- function(x, ..., zero.policy=NULL, iseed=NULL,
    no_repeat_in_row=FALSE) {
  UseMethod("localC_perm")
}

localC_perm.default <- function(x, listw, nsim = 499, alternative = "two.sided",
             ..., zero.policy=attr(listw, "zero.policy"), iseed=NULL, no_repeat_in_row=FALSE) {

  alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
  # checks are inherited from localC no need to implement
  obs <- localC(x, listw, zero.policy=zero.policy)

  # if sf object remove geometry & cast as df
  if (inherits(x, "sf")) {
    x[[attr(x, "sf_column")]] <- NULL
    x <- as.data.frame(x)
  }
  
  if (inherits(x, "list")) {
    xorig <- as.matrix(Reduce(cbind, x))
    x <- scale(xorig)
    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 (inherits(x, c("matrix", "data.frame"))) {
    xorig <- as.matrix(x)
    x <- scale(xorig)
    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 (is.vector(x) & is.numeric(x)) {
    xorig <- as.matrix(x)
    x <- scale(xorig)
    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
    cluster[obs > reps[, 1]] <- 2L
    cluster <- factor(cluster, levels=1:3, labels=c("Positive", "Negative",
      "Undefined"))
  } else {
    a <- scale(c(xorig), scale=FALSE)
    b <- lag(listw, a)
    q <- rep(4L, length(a))
    q[a > 0 & b > 0] <- 1L
    q[a <= 0 & b > 0] <- 3L
    q[a <= 0 & b <= 0] <- 2L
    cluster <- factor(q, levels=1:4, labels=c("High-High", "Low-Low",
      "Other Positive", "Negative"))
  }


  attr(obs, "call") <- match.call()
  attr(obs, "pseudo-p") <- reps
  attr(obs, "cluster") <- cluster
  class(obs) <- c("localC", "numeric")

  obs

}

localC_perm.formula <- function(formula, data, listw,
                                nsim = 499, alternative = "two.sided", ...,
                                zero.policy=attr(listw, "zero.policy"), iseed=NULL, 
                                no_repeat_in_row=FALSE) {

  alternative <- match.arg(alternative, c("less", "two.sided", "greater"))
  # if any data issues the localC formula method will catch it
  obs <- localC(formula, listw, data, zero.policy=zero.policy)

  # check if sf object in data
  if (inherits(data, "sf")) {
    data[[attr(data, "sf_column")]] <- NULL
    data <- as.data.frame(data)
  }

  xorig <- model.frame(formula, data = data)
  x <- scale()

  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
    cluster[obs > reps[, 1]] <- 2L
    cluster <- factor(cluster, levels=1:3, labels=c("Positive", "Negative",
      "Undefined"))
  } else {
    a <- scale(c(xorig), scale=FALSE)
    b <- lag(listw, a)
    q <- rep(4L, length(a))
    q[a > 0 & b > 0] <- 1L
    q[a <= 0 & b > 0] <- 3L
    q[a <= 0 & b <= 0] <- 2L
    cluster <- factor(q, levels=1:4, labels=c("High-High", "Low-Low",
      "Other Positive", "Negative"))
  }

  attr(obs, "call") <- match.call()
  attr(obs, "pseudo-p") <- reps
  attr(obs, "cluster") <- cluster
  class(obs) <- c("localC", "numeric")

  obs

}

# "localC" cluster nlevel==3L (multi) c("Positive", "Negative", "Undefined")

# "localC" cluster nlevel==4L (uni) c("High-High", "Low-Low", "Other Positive", "Negative")

# Local Geary Utils -------------------------------------------------------
localC_calc <- function(x, listw, zero.policy=attr(listw, "zero.policy")) {
  if (any(card(listw$neighbours) == 0L)) {
    res <- geary.intern(x, listw, n=length(listw$neighbours), zero.policy=zero.policy)
  } else {
    xij <- lapply(listw$neighbours, FUN = function(nbs_i) x[nbs_i])
# xij: list of vectors: for each i, x[j] values of its n_i neighbours
    res <- mapply(function(x, j, wi) sum(wi * (j - x)^2),
                x, xij, listw$weights,
                USE.NAMES = FALSE)
# res: numeric vector: for each i, the sum over j of w_{ij} * (x[j] - x[i])^2
  }
  res
}

localC_perm_calc <- function(x, listw, obs, nsim, alternative="two.sided",
  zero.policy=attr(listw, "zero.policy"), iseed=NULL, no_repeat_in_row=FALSE) {
    nc <- ncol(x)
    stopifnot(nc > 0L)
    n <- length(listw$neighbours)

    if (n != nrow(x))stop("Different numbers of observations")
    probs <- probs_lut(stat="C", nsim=nsim, alternative=alternative)
    Prname <- attr(probs, "Prname")

    crd <- card(listw$neighbours)
    z <- scale(x)
    lww <- listw$weights
    env <- new.env()
    assign("z", z, envir=env)
    assign("crd", crd, envir=env)
    assign("lww", lww, envir=env)
    assign("nsim", nsim, envir=env)
    assign("obs", obs, envir=env)
    assign("nc", nc, envir=env)
    assign("n", n, envir=env)
    assign("no_repeat_in_row", no_repeat_in_row, envir=env)
    varlist <- ls(envir = env)

    permC_int <- function(i, env) {
#    permC_int <- function(i, zi, z_i, crdi, wtsi, nsim, Ci, nc) {
      res_i <- rep(as.numeric(NA), 8)
      crdi <- get("crd", envir=env)[i]
      no_repeat_in_row <- get("no_repeat_in_row", envir=env)
      if (crdi > 0) { # if i has neighbours
        nsim <- get("nsim", envir=env)
        n_i <- get("n", envir=env) - 1L
        zi <- get("z", envir=env)[i,,drop=FALSE]
        z_i <- get("z", envir=env)[-i,]
        wtsi <- get("lww", envir=env)[[i]]
        nc <- get("nc", envir=env)
        if (nc == 1L) { # if univariate
          if (no_repeat_in_row) {
            samples <- .Call("perm_no_replace", as.integer(nsim),
              as.integer(n_i), as.integer(crdi), PACKAGE="spdep")
            sz_i <- matrix(z_i[samples], ncol=crdi, nrow=nsim)
          } else {
            sz_i <- matrix(sample(z_i, size=crdi*nsim, replace=TRUE),
              ncol=crdi, nrow=nsim) # permute nsim*#neighbours from z[-i]
          }
          diffs <- (c(zi) - sz_i)^2
          res_p <- c(diffs %*% wtsi)
        } else { # else multivariate
          res_p <- numeric(length=nsim) # for cumulation across columns
          if (no_repeat_in_row) {
            sii <- .Call("perm_no_replace", as.integer(nsim),
              as.integer(n_i), as.integer(crdi), PACKAGE="spdep")
          } else {
            sii <- sample.int(nrow(z_i), size=crdi*nsim, replace=TRUE)
          }
          for (j in 1:nc) { # permute nsim*#neighbours row indices from z[-i]
            # create nsim by crdi matrix of z[-i, j] for j-th column
            sz_i <- matrix(z_i[sii, j], ncol=crdi, nrow=nsim) 
            diffs <- (zi[, j] - sz_i)^2
            res_p <- res_p + c(diffs %*% wtsi) # cumulate across columns
          }
          res_p <- res_p/nc # nsim simulated local Gi
        }
        # res_p length nsim for obs i conditional draws
        res_i[1] <- mean(res_p)
        res_i[2] <- var(res_p)
        Ci <- get("obs", envir=env)[i]
        res_i[5] <- rank(c(res_p, Ci))[(nsim + 1L)]
        res_i[6] <- as.integer(sum(res_p >= Ci))
        res_i[7] <- e1071::skewness(res_p)
        res_i[8] <- e1071::kurtosis(res_p)
      }
      res_i
    }

    res <- run_perm(fun=permC_int, idx=1:n, env=env, iseed=iseed, varlist=varlist)

    res[,3] <- (obs - res[,1])/sqrt(res[,2])
    if (alternative == "two.sided")
      res[,4] <- 2 * pnorm(abs(res[,3]), lower.tail=FALSE)
    else if (alternative == "greater")
      res[,4] <- pnorm(res[,3], lower.tail=FALSE)
    else res[,4] <- pnorm(res[,3])
    res[,5] <- probs[as.integer(res[,5])]
    rnk0 <- as.integer(res[,6])
    drnk0 <- nsim - rnk0
    rnk <- ifelse(drnk0 < rnk0, drnk0, rnk0)
    res[,6] <- (rnk + 1.0) / (nsim + 1.0)
    
    colnames(res) <- c("E.Ci", "Var.Ci", "Z.Ci", Prname,
      paste0(Prname, " Sim"), "Pr(folded) Sim", "Skewness", "Kurtosis")
    res
}