# x: 自变量的观测值

# y: 因变量的观测值

# weight: 权重矩阵，表示观测点之间的空间关系

# scaled: 是否对观测值和权重进行标准化，默认为 FALSE

# na.rm: 是否移除 x 和 y 中的缺失值，默认为 FALSE

# alternative: 假设检验的备择假设类型，可选 "two.sided", "less", "greater"

，

bivar_Moran.I = function(x, y, weight, scaled = FALSE, na.rm = FALSE, alternative = "two.sided") {

  # 检查权重矩阵是否为方阵

  if (dim(weight)[1] != dim(weight)[2]) 

      stop("'weight' must be a square matrix")

  # 检查权重矩阵的行数和列数是否与观测值 x 和 y 的长度相等

  n <- length(x)

  if (dim(weight)[1] != n || dim(weight)[2] != n) 

      stop("'weight' must have as many rows and columns as observations in 'x' and 'y'")

  # 计算常量 ei

  ei <- -1/(n - 1)

  # 检查 x 和 y 是否有缺失值

  nas <- is.na(x) | is.na(y)

  if (any(nas)) {

    if (na.rm) {

        # 移除缺失值

        x <- x[!nas]

        y <- y[!nas]

        n <- length(x)

        weight <- weight[!nas, !nas]

    } else {

        # 若未移除缺失值，则返回 NA 值

        warning("'x' or 'y' has missing values: maybe you wanted to set na.rm = TRUE?")

        return(list(observed = NA, expected = ei, sd = NA, p.value = NA))

    }

  }

  # 标准化权重矩阵

  ROWSUM <- rowSums(weight)

  ROWSUM[ROWSUM == 0] <- 1

  weight <- weight/ROWSUM

  s <- sum(weight)

  # 中心化 x 和 y

  mx <- mean(x)

  my <- mean(y)

  x <- x - mx

  y <- y - my

  # 计算 Moran's I 统计量

  cv <- sum(weight * x %o% y)

  v.x <- sum(x^2)

  v.y <- sum(y^2)

  obs <- (n/s) * (cv/sqrt(v.x * v.y))

  # 是否进行标准化

  if (scaled) {

    i.max <- (n/s) * (sd(rowSums(weight) * x)/sqrt(v.x/(n - 1)))

    obs <- obs/i.max

  }

  # 计算统计量的方差

  S1 <- 0.5 * sum((weight + t(weight))^2)

  S2 <- sum((apply(weight, 1, sum) + apply(weight, 2, sum))^2)

  s.sq <- s^2

k <- (sum(x^2 * y^2)/n^2)/(v.x/n)^2

  sdi <- sqrt((n * ((n^2 - 3 * n + 3) * S1 - n * S2 + 3 * s.sq) - 

                 k * (n * (n - 1) * S1 - 2 * n * S2 + 6 * s.sq))/((n - 

                                                                     1) * (n - 2) * (n - 3) * s.sq) - 1/((n - 1)^2))

  alternative <- match.arg(alternative, c("two.sided", 

                                          "less", "greater"))

  pv <- pnorm(obs, mean = ei, sd = sdi)

  if (alternative == "two.sided") 

    pv <- if (obs <= ei) 

      2 * pv

  else 2 * (1 - pv)

  if (alternative == "greater") 

    pv <- 1 - pv

  list(observed = obs, expected = ei, sd = sdi, p.value = pv)

}