1

我有一个组数据框，我想为其计算每个组的马氏距离。

我正在尝试将 Mahalanobis 函数应用于数百个组，并且一个特定的组由于样本量小（只有两行）而导致问题。

我的数据如下所示：

foo <- data.frame(GRP = c("a","a","b","b","b"),X = c(1,1,15,12,50),
                      Y = c(2.17,12.44,50,70,100))

我从这里借用了一个函数的想法，它看起来如下：

auto.mahalanobis <- function(temp) {
 mahalanobis(temp, 
             center = colMeans(temp, na.rm=T), 
             cov = cov(temp, use="pairwise.complete.obs"),
             tol=1e-20,
             inverted = FALSE
             )
}

根据此处的建议，我将tol参数添加到auto.mahalanobis函数中，以避免在计算小数的协方差矩阵时出现问题。

然后我尝试将此函数与我的数据集一起使用，并收到以下关于奇异矩阵的错误：

 z <- foo %>% group_by(GRP) %>% mutate(mahal = auto.mahalanobis(data.frame(X,Y)))

Error: Problem with `mutate()` input `mahal`.
x Lapack routine dgesv: system is exactly singular: U[1,1] = 0
i Input `mahal` is `auto.mahalanobis(data.frame(X, Y))`.
i The error occurred in group 1: GRP = "a".

相同的功能适用于具有较大样本量的其他组，是否有建议的方法来解决此问题或在样本太小时跳过此类组？

4

1 回答 1

0

最简单的方法可能是：

auto.mahalanobis <- function(temp) {
 m <- try(silent=TRUE,
           mahalanobis(temp, 
             center = colMeans(temp, na.rm=TRUE), 
             cov = cov(temp, use="pairwise.complete.obs"),
             tol=1e-20,
             inverted = FALSE
             ))
 if (!inherits(m,"try-error")) return(m)
 return(rep(NA_real_, length(temp))
}

（未经测试：真正的程序员可能会使用tryCatch()）

如果您认为只有在n==2可以使用if子句时才会出现问题，例如if (length(temp)<min_length) return(rep(NA_real_,length(temp))).

或者，您可以制作一个mahalanobis()使用广义逆 ( MASS::ginv) 而不是常规矩阵求逆 ( ) 的黑客版本solve；我认为这可能（？）可以作为替代品，但尚未检查数学。

于 2021-01-20T22:23:02.930 回答