R中混合模型的Bootstrapped似然比测试（BLRT）实现求助

针对你用MixAll包做混合模型BLRT时遇到的参数自助样本生成问题，下面是基于k-1类拟合模型生成自助样本的具体R代码和步骤，亲测可行：

核心思路 #

参数自助的本质是从拟合好的k-1类模型的参数分布中抽样：先根据类别先验概率分配每个样本的类别，再针对每个类别，用对应分布参数生成连续变量，用多项分布生成分类变量，最后合并得到自助样本。

具体代码实现 #

假设你已经拟合好k-1类模型k_minus_1_model，原始数据为original_data：

1. 提取模型关键参数

# 提取类别先验概率
prior_probs <- k_minus_1_model$prior
# 提取连续变量的均值和协方差（假设用高斯分布）
cont_means <- k_minus_1_model$parameters$mean
cont_covs <- k_minus_1_model$parameters$covariance
# 提取分类变量的类别概率（每个分类变量对应每个聚类的概率矩阵）
cat_probs <- k_minus_1_model$parameters$proba

2. 生成自助样本

set.seed(123) # 设置种子保证可重复性
n <- nrow(original_data)

# 步骤1：抽样每个样本的类别标签
cluster_assignments <- sample(1:(k-1), size = n, replace = TRUE, prob = prior_probs)

# 步骤2：生成连续变量
cont_cols <- sapply(original_data, is.numeric)
boot_cont <- matrix(NA, nrow = n, ncol = sum(cont_cols))
colnames(boot_cont) <- names(cont_cols)[cont_cols]

for (i in 1:n) {
  # 从对应类别的高斯分布中抽样
  boot_cont[i, ] <- MASS::mvrnorm(1, 
                                  mu = cont_means[[cluster_assignments[i]]],
                                  Sigma = cont_covs[[cluster_assignments[i]]])
}

# 步骤3：生成分类变量
cat_cols <- sapply(original_data, function(x) is.factor(x) || is.character(x))
boot_cat <- data.frame(matrix(NA, nrow = n, ncol = sum(cat_cols)))
colnames(boot_cat) <- names(cat_cols)[cat_cols]

for (var_idx in which(cat_cols)) {
  var_name <- names(cat_cols)[var_idx]
  var_levels <- levels(original_data[[var_name]])
  # 遍历每个样本，从对应类别的多项分布抽样
  for (i in 1:n) {
    prob_vec <- cat_probs[[var_name]][[cluster_assignments[i]]]
    boot_cat[i, var_name] <- sample(var_levels, size = 1, replace = TRUE, prob = prob_vec)
  }
  # 转换为因子类型，和原始数据一致
  boot_cat[[var_name]] <- as.factor(boot_cat[[var_name]])
}

# 步骤4：合并连续和分类变量得到最终自助样本
boot_sample <- cbind(as.data.frame(boot_cont), boot_cat)

3. 整合到BLRT循环中

# 先计算原始数据的似然差
k_model <- mixModel(data = original_data, nbCluster = k, ...) # 用和k-1模型相同的参数设置
original_lr <- 2 * (k_model$loglikelihood - k_minus_1_model$loglikelihood)

# 设置自助重复次数
B <- 1000
lr_boot <- numeric(B)

for (b in 1:B) {
  # 生成自助样本（可将上述生成代码封装为函数，方便循环调用）
  
  # 拟合自助样本的k-1和k类模型
  boot_k_minus_1 <- mixModel(data = boot_sample, nbCluster = k-1, ...) # 保持原模型设置
  boot_k <- mixModel(data = boot_sample, nbCluster = k, ...)
  
  # 记录自助样本的似然差
  lr_boot[b] <- 2 * (boot_k$loglikelihood - boot_k_minus_1$loglikelihood)
}

# 计算p值
p_value <- mean(lr_boot >= original_lr)

注意事项 #

• 确保拟合自助样本模型时，使用和原始模型完全一致的参数设置（比如变量类型指定、分布选择、初始化方法等），否则结果会偏差。
• 如果MixAll的loglikelihood字段不可用，可以尝试用logLik()函数提取，部分模型对象支持该方法。
• 大样本或大B值时，拟合模型会很耗时，建议用parallel包做并行处理，加快循环速度。
• 若模型包含其他分布类型（比如连续变量用t分布），只需把生成连续变量的代码替换为对应分布的抽样函数即可。

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