为猿包 Moran's I 函数在 R 中解决“if (obs <= ei) 2 * pv else 2 * (1 - pv) : 需要 TRUE/FALSE 的缺失值”中的错误

如何解决为猿包 Moran's I 函数在 R 中解决“if (obs <= ei) 2 * pv else 2 * (1 - pv) : 需要 TRUE/FALSE 的缺失值”中的错误

开发人员！

我遇到了一条错误消息（“错误 in if (obs

library(ape)

nrstp <- data.frame(
     X = c(300226.9,300224.6,300226.4,300226.1,300224.0,300225.7,300226.3,300227.1),Y = c(5057949,5057952,5057950,5057956,5057949),V3 = c(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 米缓冲区内其他线串最近的点）？

提前感谢您的帮助！

解决方法

控制流选择 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,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) 是否总是正确的解决方案。