求解猿包的"Error in if (obs <= ei) 2 * pv else 2 * (1 - pv) : missing value where TRUE/FALSE needed" 莫

Question

developers！

我遇到了一条错误消息

if 中的错误 (obs <= ei) 2 * pv else 2 * (1 - pv) ： 缺少值 其中 需要真/假

阻止我从 Moran 的 I 函数中从猿包中获取值。这是我所做的：

library(ape)
nrstp <- data.frame(
X = c(300226.9, 300224.6, 300226.4, 300226.1, 300224.0, 300226.4, 300225.7, 300226.4, 300226.1, 300226.4, 300226.3, 300226.3, 300227.1),
Y = c(5057949, 5057952, 5057950, 5057950, 5057956, 5057950, 5057950, 5057950, 5057950, 5057950, 5057950, 5057950, 5057949),
V3 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
nrstp = data.frame(nrstp)
dist = as.matrix(dist(cbind(nrstp$X, nrstp$Y)))
invdist = 1/dist
invdist[is.infinite(invdist)] <- 0
moranI = Moran.I(nrstp$V3, invdist)

此代码的目的是从一系列点计算 Moran 的 I，以检查空间自相关。到目前为止，这似乎是 Moran 在 R 中的 I 中唯一起作用的函数。经过几次测试(我有数千组点)，此错误似乎只发生在只有一个值的输入向量上(我尝试了 0 以外的其他数字，它仍然引发此错误)。

有人可以帮助我改进此代码吗？或者他们更好的建议是计算莫兰 I 或从线串测试空间自相关(这些点组是一个线串的原点和该原点 10 米缓冲区内与其他线串最近的点)？

提前感谢您的任何帮助！

Answer 1

控制流选择if(condition) do something要求condition的值不NA。

在您的情况下，obs <= ei会导致NA。这就是生成错误消息missing value where TRUE/FALSE needed的原因。

要了解obs <= ei如何产生NA，您可以在Moran.I函数中检查详细信息：

Moran.I
function (x, weight, scaled = FALSE, na.rm = FALSE, alternative = "two.sided") 
{
if (dim(weight)[1] != dim(weight)[2]) 
stop("'weight' must be a square matrix")
n <- length(x)
if (dim(weight)[1] != n) 
stop("'weight' must have as many rows as observations in 'x'")
ei <- -1/(n - 1)
nas <- is.na(x)
if (any(nas)) {
if (na.rm) {
x <- x[!nas]
n <- length(x)
weight <- weight[!nas, !nas]
}
else {
warning("'x' 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)
m <- mean(x)
y <- x - m
cv <- sum(weight * y %o% y)
v <- sum(y^2)
obs <- (n/s) * (cv/v)
if (scaled) {
i.max <- (n/s) * (sd(rowSums(weight) * y)/sqrt(v/(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(y^4)/n)/(v/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)
}
<bytecode: 0x000001cd5e0715d0>
<environment: namespace:ape>

通过分配x = nrstp$V3和weight = invdist，您将获得mean(x) = 0。这导致y=0、cv = 0、v=0，最后obs = NaN。 因此

obs <= ei
[1] NA

要克服这个问题，您需要确保obs和ei中的每一个都不是NA。在您的情况下，如果mean(x)不为零，则obs <= ei将不会NA。但是，由于我对这个特定主题一无所知，因此我不确定非零mean(x)是否总是正确的解决方案。

Answer 2

问题是你的x都是相同的值。如果你查看Abdur Rohman的代码，函数的计算是

m <- mean(x)
y <- x - m
cv <- sum(weight * y %o% y)
v <- sum(y^2)
obs <- (n/s) * (cv/v)

如果所有x都相同，则m <- mean(x)的平均值显然与所有x值相同，并且y, v, obs为0。 对于obs，你除以cv/v这是NaN

因此，x至少有一个值应该不同