问题描述
我使用 ContaminatedMixt package、CNmixt function 对多维数据进行聚类。我的数据维度范围从 3 到 9。大多数数据可以毫无问题地聚类，但我一直在努力对 7D 数据进行聚类。
对于无法包含独立的 reprex 表示歉意，但我知道我自己的数据可靠地发生了错误。您可以通过 GoogleDrive 下载 .csv of the data。
重现错误：
library(ContaminatedMixt)
df = read.csv (file = "210102_7clickcodas.csv",header = TRUE,fileEncoding = "UTF-8-BOM")
df = round(df,3)
clusters = CNmixt (
df,G = 2:10,contamination = TRUE,model = NULL,initialization = "kmeans",parallel = TRUE)
我要求 CNmixt 做的是将 2，然后是 3，然后是 4，...，然后是 10 个集群拟合到数据并选择最佳的集群数量（基于一个标准 - 我使用 ICL，但是BIC、AIC等也是选项）。但是，当我运行此代码时，总是遇到以下错误：
Error in checkForRemoteErrors(val) :
one node produced an error: BLAS/LAPACK routine 'DGEBAL' gave error code -3
有时多个节点产生错误，但错误信息总是相同的。我试图通过将 G 设置为 2:10 范围内的每个值而不是范围本身（G=2，然后 G=3，然后 G=4 等）来解决导致问题的 G 值。看起来 CNmixt 不能将 7D 数据分成 4 个集群。这样做会导致错误抛出并且程序停止运行。
最重要的是，我想： a) 告诉 CNmixt 程序绕过任何导致错误的 G 值并继续。因此，如果 G=4 是问题，程序将拟合 2 个集群、3 个集群，尝试拟合 4 个但失败，然后继续拟合 5 个集群。在拟合 2:9 集群结束时，它将选择优化标准并成功拟合的集群数量。我认为这将涉及直接调整 CNmixt source code（我认为是 .CNmixtG2 函数），但是在不操作 CNmixt 源代码的情况下做到这一点的方法也很棒。
b) 了解为什么会出现此错误也会很棒。我知道其他人之前曾询问过此错误（例如，r msm BLAS/LAPACK routine 'DGEBAL' gave error code -3），并且根据我所阅读的内容，错误代码 -3 似乎意味着第三个 INFO 参数中存在非法值。但事实证明，找出 CNmixt 代码中“DGEBAL”函数被调用的位置、它究竟在做什么以及它为什么失败对我来说真的很有挑战性。
任何有关如何执行此操作的建议将不胜感激！
解决方法
这种方式可能对你有用。如果您查看 CNmixt() 函数的 R 代码，它只是重复调用 CNmixt_main() 函数，尽管该函数不是从命名空间中导出的。我做的第一件事是使用 expand.grid() 设置搜索网格，以获取集群数量和模型的所有组合。
eg <- expand.grid(G=2:10,model=c("EII","VII","EEI","VEI","EVI","VVI","EEE","VEE","EVE","EEV","VVE","VEV","EVV","VVV"),stringsAsFactors=FALSE)
接下来，我初始化一个输出对象 (ICS)，然后循环遍历 eg 的行。对于每一行（模型和集群数量的组合），我估计包裹在 try() 中的模型，这样如果模型抛出错误，它就不会停止循环。对象 out 将只包含错误，它属于 "try-error" 类。如果函数没有因错误而终止，则返回 IC 值，否则返回一个 NA 值向量，该向量与否则将返回的 IC 值向量一样长。最后，将 IC 值添加为 ICS 的新行，该行累积每个模型的值。
ICS <- NULL
for(i in 1:nrow(eg)){
out <- try(ContaminatedMixt:::CNmixt_main (
df,G = eg$G[i],contamination = TRUE,model = eg$model[i],initialization = "kmeans",alphafix = NULL,alphamin = .5,seed = NULL,start.z=NULL,start.v = NULL,start = 0,label = NULL,AICcond = FALSE,iter.max = 1000,threshold = 1e-10,parallel=FALSE,eps=1e-100,verbose=TRUE,doCV = FALSE))
if(inherits(out,"try-error")){
ics <- rep(NA,8)
}else{
ics <- unlist(out$models[[1]]$IC)
}
ICS <- rbind(ICS,ics)
}
接下来，我们可以将 IC 值附加到我们搜索的参数上。
eg <- bind_cols(eg,as_tibble(ICS))
接下来，我们可以根据集群大小找到最佳 ICL 值。这是假设最好的 ICL 是最小的，如果它是最大的，只需将 min(ICL,na.rm=TRUE) 替换为 max(ICL,na.rm=TRUE)
best_per_cluster <- eg %>%
group_by(G) %>%
filter(ICL == min(ICL,na.rm=TRUE)) %>%
arrange(G)
best_per_cluster
# # A tibble: 9 x 10
# # Groups: G [9]
# G model AIC BIC AIC3 CAIC AWE ICL AICc AICu
# <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 2 VII 28960. 28854. 28941. 28835. 28653. 28724. 28960. 28940.
# 2 3 EVI 31190. 30955. 31148. 30913. 30511. 30786. 31188. 31145.
# 3 4 VEI 34230. 33984. 34186. 33940. 33518. 33702. 34228. 34183.
# 4 5 EVI 35896. 35504. 35826. 35434. 34763. 35373. 35891. 35818.
# 5 6 EEI 34959. 34629. 34900. 34570. 34004. 34446. 34955. 34894.
# 6 7 EVI 38211. 37663. 38113. 37565. 36625. 37450. 38201. 38099.
# 7 8 EEI 39301. 38870. 39224. 38793. 38054. 38058. 39294. 39215.
# 8 9 EEI 37999. 37518. 37913. 37432. 36607. 37259. 37991. 37902.
# 9 10 EII 34709. 34205. 34619. 34115. 33252. 33900. 34700. 34607.
最后，我们可以选择最好的模型，再次假设最好的 ICL 是最小的。
best_per_cluster %>%
ungroup %>%
filter(ICL == min(ICL))
# # A tibble: 1 x 10
# G model AIC BIC AIC3 CAIC AWE ICL AICc AICu
# <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 2 VII 28960. 28854. 28941. 28835. 28653. 28724. 28960. 28940.
显然，获得每个集群大小的最佳模型是一个不必要的步骤，但同样的观察可能会很有趣。
此外，我还没有回答您关于错误发生原因的其他问题。当我运行 CNmixt() 函数时，直到 G=10 才出现错误，但是当我使用 CNmixt_main() 运行模型时，我在该过程的早期出现了错误（例如，G=3 和对于 EVE 模型，G=4）。在我的特殊情况下，错误并不像对您来说那么可靠。我的会话信息如下仅供参考。
> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.4
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ContaminatedMixt_1.3.4.1
loaded via a namespace (and not attached):
[1] Rcpp_1.0.5 pillar_1.4.6 compiler_4.0.3
[4] gower_0.2.1 plyr_1.8.6 class_7.3-17
[7] iterators_1.0.12 tools_4.0.3 rpart_4.1-15
[10] ipred_0.9-9 mclust_5.4.6 lubridate_1.7.9
[13] mixture_1.5.1 lifecycle_0.2.0 tibble_3.0.3
[16] gtable_0.3.0 nlme_3.1-149 lattice_0.20-41
[19] pkgconfig_2.0.3 rlang_0.4.7 Matrix_1.2-18
[22] foreach_1.5.0 rstudioapi_0.11 parallel_4.0.3
[25] mvtnorm_1.1-1 prodlim_2019.11.13 stringr_1.4.0
[28] withr_2.2.0 dplyr_1.0.2 pROC_1.16.2
[31] generics_0.0.2 vctrs_0.3.4 recipes_0.1.12
[34] stats4_4.0.3 nnet_7.3-14 grid_4.0.3
[37] caret_6.0-86 tidyselect_1.1.0 data.table_1.13.0
[40] glue_1.4.2 R6_2.4.1 survival_3.2-7
[43] lava_1.6.7 reshape2_1.4.4 ggplot2_3.3.2
[46] purrr_0.3.4 magrittr_1.5 ModelMetrics_1.2.2.2
[49] splines_4.0.3 MASS_7.3-53 scales_1.1.1
[52] codetools_0.2-16 ellipsis_0.3.1 mnormt_2.0.1
[55] timeDate_3043.102 colorspace_1.4-1 stringi_1.5.3
[58] munsell_0.5.0 tmvnsim_1.0-2 crayon_1.3.4