基于mclust初始聚类模型分配新数据聚类的异常问题咨询

背景 #

我们通过在线调研获取初始数据，使用library(mclust)开展聚类分析，之后又获取了规模更小的第二轮调研数据（seconddata），希望用基于初始数据（firstdata）构建的模型为新数据分配聚类标签。

初始流程：先对firstdata做标准化（scale）再聚类；为新数据分配标签时，同样对seconddata做标准化，但结果出现异常（比如原占比最高的聚类出现空簇），怀疑是数据集分布差异导致，但不确定问题所在，也不确定当前方法是否可行，需要正确匹配新数据到已有聚类的方法。

用户提供的示例代码（流程与真实数据一致）：

# library -----------------------------------------------------------------
library(tidyverse)
library(ggplot2)
library(scatterplot3d)
library(hopkins)
library(factoextra)
library(NbClust)
library(mclust)
library(ggpubr)
library(writexl)
library(reshape2)
library(hms)
library("lubridate") 

set.seed(1591593)

# first dataset -----------------------------------------------------------

n1 = 100
#number of observations in first data frame

firstdata <- expand.grid(n1 = 1:n1,
                         x1 = NA,
                         x2 = NA,
                         x3 = NA)


firstdata$x1 <- sample(runif(n = n1, min = 50, max = 100), replace = TRUE)
firstdata$x2 <- sample(runif(n = n1, min = 20, max = 25), replace = TRUE)
firstdata$x3 <- sample(runif(n = n1, min = 30, max = 80), replace = TRUE)
#creating data

firstdata = firstdata %>%
  select(contains("x")) %>%
  mutate_all(scale)
#scale data

par(mfrow = c(1, 1)) # display 1 plot
scatterplot3d(firstdata$x1, firstdata$x2, firstdata$x3)


## model-based clustering ---------------------------------------------

# Fit model
c_model <- Mclust(firstdata)

# Show optimal model
summary(c_model) 

# Plot all models
fviz_mclust(c_model, "BIC", palette = "jco")

# Print cluster sizes
table(c_model$classification)

# Cluster plot
fviz_mclust(c_model, "classification", geom = "point", 
            pointsize = 1.5, palette = "jco")

# Cluster no. (classification) as variable (new dataframe)
firstdata_cluster = firstdata
firstdata_cluster$cluster = as.factor(c_model$classification)

# Compute probability to belong to resp. cluster and not to other cluster
firstdata_cluster$probability = 1-(c_model$uncertainty)

# Cluster label in classifications
clusters = MclustDA(firstdata, class = firstdata_cluster$cluster)
summary(clusters)


# second dataset ----------------------------------------------------------

n2 = 10
#number of observations in second data frame

seconddata <- expand.grid(n2= 1:n2,
                            x1 = NA,
                            x2 = NA,
                            x3 = NA)

seconddata$x1 <- sample(runif(n = n2, min = 50, max = 100), replace = TRUE)
seconddata$x2 <- sample(runif(n = n2, min = 20, max = 25), replace = TRUE)
seconddata$x3 <- sample(runif(n = n2, min = 30, max = 80), replace = TRUE)
# create second data frame


## assigning clusters to the new data ------------------------------------

seconddata = seconddata %>%
  select(contains("x")) %>%
  mutate_all(scale)
# scale second data frame

prediction_seconddata = predict(clusters, seconddata)
round(prediction_seconddata$z, 0)
#assign clusters to the new second data frame/new data using the first model

prediction_seconddata_cluster<-as.data.frame(round(prediction_seconddata$z, 0))
prediction_seconddata_cluster
#assigned clusters of the second data frame

核心问题与解决方法 #

1. 标准化方式错误（最关键） #

问题根源：对seconddata单独调用scale()，会使用新数据自身的均值和标准差做标准化，而非初始firstdata的统计量，导致新数据分布与训练模型时的数据分布不一致，进而引发聚类分配异常。

修正方法：保存firstdata标准化时的均值和标准差，用这些参数对seconddata做标准化：

# 处理初始数据，保存标准化参数
firstdata_raw <- firstdata %>% select(contains("x"))
# 提取每个变量的均值和标准差
scale_params <- lapply(firstdata_raw, function(col) {
  list(mean = mean(col, na.rm = TRUE), sd = sd(col, na.rm = TRUE))
})

# 用保存的参数标准化firstdata
firstdata_scaled <- as.data.frame(
  mapply(function(col, params) (col - params$mean)/params$sd,
         firstdata_raw, scale_params, SIMPLIFY = FALSE)
)
colnames(firstdata_scaled) <- names(firstdata_raw)

# 用相同参数标准化seconddata
seconddata_raw <- seconddata %>% select(contains("x"))
seconddata_scaled <- as.data.frame(
  mapply(function(col, params) (col - params$mean)/params$sd,
         seconddata_raw, scale_params, SIMPLIFY = FALSE)
)
colnames(seconddata_scaled) <- names(seconddata_raw)

2. 直接使用Mclust模型预测，无需MclustDA #

问题根源：MclustDA是用于判别分析的工具，而我们只需要基于已有的聚类模型分配标签，直接使用原始的c_model（Mclust对象）做预测更直接且准确。

修正方法：

# 用初始聚类模型直接预测新数据
pred_result <- predict(c_model, seconddata_scaled)

# 查看分配的聚类标签
pred_result$classification

# 查看每个样本属于各聚类的概率
pred_result$z

3. 验证数据分布差异（可选） #

如果修正后仍有异常，可先验证两轮数据的分布差异：

• 绘制变量的直方图/箱线图对比firstdata_raw和seconddata_raw的分布
• 计算变量的均值、标准差、分位数等统计量，确认差异是否在可接受范围内

完整修正后的核心流程 #

# 1. 加载库与设置随机种子
library(mclust)
library(tidyverse)
set.seed(1591593)

# 2. 处理初始数据，保存标准化参数
n1 = 100
firstdata_raw <- tibble(
  x1 = sample(runif(n = n1, min = 50, max = 100), replace = TRUE),
  x2 = sample(runif(n = n1, min = 20, max = 25), replace = TRUE),
  x3 = sample(runif(n = n1, min = 30, max = 80), replace = TRUE)
)

# 提取标准化参数
scale_params <- lapply(firstdata_raw, function(col) {
  list(mean = mean(col), sd = sd(col))
})

# 标准化初始数据
firstdata_scaled <- firstdata_raw %>%
  mutate(
    x1 = (x1 - scale_params$x1$mean)/scale_params$x1$sd,
    x2 = (x2 - scale_params$x2$mean)/scale_params$x2$sd,
    x3 = (x3 - scale_params$x3$mean)/scale_params$x3$sd
  )

# 3. 构建聚类模型
c_model <- Mclust(firstdata_scaled)
summary(c_model)

# 4. 处理新数据，用相同参数标准化
n2 = 10
seconddata_raw <- tibble(
  x1 = sample(runif(n = n2, min = 50, max = 100), replace = TRUE),
  x2 = sample(runif(n = n2, min = 20, max = 25), replace = TRUE),
  x3 = sample(runif(n = n2, min = 30, max = 80), replace = TRUE)
)

seconddata_scaled <- seconddata_raw %>%
  mutate(
    x1 = (x1 - scale_params$x1$mean)/scale_params$x1$sd,
    x2 = (x2 - scale_params$x2$mean)/scale_params$x2$sd,
    x3 = (x3 - scale_params$x3$mean)/scale_params$x3$sd
  )

# 5. 预测新数据的聚类标签
pred_result <- predict(c_model, seconddata_scaled)
# 输出分配的聚类标签
cat("新数据聚类标签：\n")
print(pred_result$classification)
# 输出各聚类概率
cat("\n新数据聚类概率：\n")
print(pred_result$z)

内容的提问来源于stack exchange，提问作者Linus