使用ape软件包在R中计算Moran's I时出现“Error in if (obs <= ei) 2 * pv else 2 * (1 - pv) : missing value where TRUE/FALSE needed”错误的解决方法

Question

开发者们！

我遇到了一个错误信息

如果(obs <= ei)，则出现错误：2 * pv，否则出现错误：2 * (1 - pv) ：需要TRUE/FALSE的缺失值。

这个错误阻止我从ape包中获取Moran's 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's I以检查空间自相关性。到目前为止，这似乎是R中唯一有效的Moran's I函数。经过几次测试（我有数千组点），发现此错误只会发生在输入向量仅具有一个值时（我尝试了其他数字而不是0，它仍然引发此错误）。

请问有人可以帮助我改进这个代码吗？或者有更好的建议来计算Moran's I或者测试线串的空间自相关性（这些点组是线串的原点以及离这样的原点10米缓冲区内最近的其他线串上的点）？

非常感谢您提前的任何帮助！

- VeggieB

---

你的 V3 变量是什么？ - Istrel

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)是否总是正确的解决方案。

- Abdur Rohman

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的值应该不同。

- Triss