|
# 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)
}
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 08:36
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社