计算有缺失值时的 Mahalanabois 距离

Question

在 R 中，我试图计算马氏距离来检查数据集中是否存在异常值，以测试多元方差分析的假设之一。我的数据集中缺少值。我最初尝试过 mahalanabois 函数，但这似乎不适用于缺失值，所以我尝试了 modi 包中的 MDmiss 函数。这适用于我的两个变量（DO 和 chla）都缺少值的情况。但是，如果我仅缺少 chla 或 DO 中的数据，则不会计算距离。当我缺少缺失值时，MDmiss 和 mahalanobis 函数都不会返回距离。

我还尝试在原始马哈拉诺比斯距离函数中使用 is.na 和 na.omit 参数，但这也不起作用。我已经提供了一个示例数据集。 感谢您的帮助。 谢谢。

envdata <- data.frame(WaterTemp = c(56.7, 56.4, 60.8,60.6, 59.3, 57.5, 57.9, 65.8,59.2, 59), SPC = c(46600, 47520, 47821, 47801, 47999, 47418, 47646, 49156, 46350, 46260), Salinity = c(30.28, 30.92, 31.54, 31.34, 31.24, 30.87, 31.03, 32.17, 30.12, 30.05), DO = c(NA, NA, 96, NA, NA, NA, NA, 101, 99, 103), Chla = c(7.045, NA, 8.358, NA, NA, NA, 6.306, 26.84, NA, NA))


#Check for outliers using the Mahalanobis distance
#https://www.statology.org/mahalanobis-distance-r/

#Mahalanobis only works on numeric data. Make new data frame with only numeric variables 
#Convert integers to numeric
envdata <- envdata %>% mutate(SPC = as.numeric(envdata$SPC), DO = as.numeric(envdata$DO))
envdata_numeric <- envdata %>% dplyr::select(WaterTemp, SPC, Salinity, DO, Chla)

#create new column in data frame to hold Mahalanobis distances
envdata_numeric$mahal <- mahalanobis(envdata_numeric, colMeans(envdata_numeric, na.rm = TRUE), cov(envdata_numeric))

#create new column in data frame to hold p-value for each Mahalanobis distance
envdata_numeric$p <- pchisq(envdata_numeric$mahal, df = 4, lower.tail = FALSE)
#Df = (c-1)
#DF = 5-1

envdata_numeric

#***#error with calculating distances. Possibly because of NA values. Try this other package. https://search.r-project.org/CRAN/refmans/modi/html/MDmiss.html
devtools::install_github("martinSter/modi")
library(modi)

#create new column in data frame to hold Mahalanobis distances
envdata_numeric$mahal <- MDmiss(envdata_numeric, colMeans(envdata_numeric), cov(envdata_numeric))

Answer 1

您显示的数据有问题，列

DO

和

Chal

共线。 也就是说，您只有两个完整的观察结果（参见下面

envdata_numeric

的第3行和第8行）：

envdata_numeric <- structure(list(WaterTemp = c(56.7, 56.4, 60.8, 60.6, 59.3, 57.5, 
57.9, 65.8, 59.2, 59), SPC = c(46600, 47520, 47821, 47801, 47999, 
47418, 47646, 49156, 46350, 46260), Salinity = c(30.28, 30.92, 
31.54, 31.34, 31.24, 30.87, 31.03, 32.17, 30.12, 30.05), DO = c(NA, 
NA, 96, NA, NA, NA, NA, 101, 99, 103), Chla = c(7.045, NA, 8.358, 
NA, NA, NA, 6.306, 26.84, NA, NA)), class = "data.frame", row.names = c(NA, 
-10L))

# WaterTemp   SPC Salinity  DO   Chla
# 1       56.7 46600    30.28  NA  7.045
# 2       56.4 47520    30.92  NA     NA
# 3       60.8 47821    31.54  96  8.358
# 4       60.6 47801    31.34  NA     NA
# 5       59.3 47999    31.24  NA     NA
# 6       57.5 47418    30.87  NA     NA
# 7       57.9 47646    31.03  NA  6.306
# 8       65.8 49156    32.17 101 26.840
# 9       59.2 46350    30.12  99     NA
# 10      59.0 46260    30.05 103     NA

粗略地说，您正在尝试查找异常值或计算距离，但是您没有足够的信息来围绕点云“绘制椭圆体”。这就是几何上

mahalanobis

正在做的事情。我概述了下面的情况：白色圆圈是没有

NA

的列，大红色表示位于更高维度的两个点（第3行和第8行）。通过 2 个点和中心可以绘制无限多个椭球体（我画了 2 个）。

无论如何，如果我将一些数据点添加到

DO

列中，例如到第 1 行

100

然后继续进行插补（我使用了

mice

包）我可以正式 计算距离。正如您将看到的 p 值 将 > 0.1。这意味着无论算法如何工作，即使根据 3 个观察值来判断异常值也是不够的。太多了

NA

。

library(mice)
envdata_numeric[1, "DO"] <- 100
envdata_numeric_imp <- complete(mice(envdata_numeric))
envdata_numeric_imp$maha <- mahalanobis(envdata_numeric_imp, 
                                        colMeans(envdata_numeric_imp), 
                                        cov(envdata_numeric_imp))

envdata_numeric_imp$p = pchisq(envdata_numeric_imp$maha, df = 4, 
                               lower.tail = FALSE)


envdata_numeric_imp

输出：

Water   Temp   SPC Salinity  DO   Chla     maha         p
1       56.7 46600    30.28 100  7.045 1.274517 0.8656841
2       56.4 47520    30.92 103  7.045 3.554027 0.4697112
3       60.8 47821    31.54  96  8.358 7.201919 0.1255948
4       60.6 47801    31.34 103  6.306 3.968202 0.4103263
5       59.3 47999    31.24  96  6.306 5.790871 0.2153200
6       57.5 47418    30.87 101 26.840 6.985705 0.1366456
7       57.9 47646    31.03 101  6.306 1.523915 0.8223970
8       65.8 49156    32.17 101 26.840 7.254101 0.1230542
9       59.2 46350    30.12  99  6.306 3.556350 0.4693622
10      59.0 46260    30.05 103  7.045 3.890395 0.4210425

Answer 2

我使用了一种新方法， I 计算矩阵 S 的行列式。 该行列式取决于未知值。 然后我最小化正行列式（二次形式），然后得到预测值