单组不均簇聚类数据的拟合优度检验R实现咨询

核心逻辑 #

普通卡方检验不适合聚类数据——同一簇内的观测存在相关性，会导致标准误被低估、假阳性率升高。必须用考虑聚类随机效应的模型调整后，再检验分类变量各水平的频率差异。

(a) 适用的R函数 #

• lme4::glmer()：拟合广义线性混合模型，处理聚类带来的随机效应
• aods3::gof()：针对混合模型做拟合优度检验，可直接验证各水平频率是否一致
• car::Anova()：配合混合模型检验固定效应的显著性
• brms::brm()（可选）：拟合更灵活的贝叶斯多项混合模型，适配复杂场景

(b) 具体模型设定与操作 #

1. 数据预处理 #

先把示例数据整理为规范数据框，并将分类变量转为因子（R中模型默认要求因子类型的分类变量）：

dat <- data.frame(
  var = factor(c("a", "a", "a", "b", "b", "c", "c", "c", "c", "d", "d", "e", "e", "e", "e", "e")),
  clusters = factor(c(1, 1, 2, 2, 3, 3, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12))
)

2. 拟合聚类调整的混合模型 #

因为var是多水平分类变量，我们拟合多项Logistic混合模型，将clusters设为随机截距，控制簇内相关性：

library(lme4)
# 以"a"为参考水平拟合模型，水平数>2时glmer自动处理为多项逻辑回归
model <- glmer(var ~ 1 + (1 | clusters), 
               data = dat, 
               family = binomial(),
               control = glmerControl(optimizer = "bobyqa"))

注：若模型收敛异常，可调整优化器（如上述的bobyqa），或改用brms包的贝叶斯模型。

3. 检验各水平频率的差异 #

方法一：用aods3::gof()做拟合优度检验

该检验直接对比观测频率与模型预测的期望频率，若p值显著，说明各水平频率存在差异：

library(aods3)
gof_result <- gof(model)
print(gof_result)

输出中的Pearson卡方或Deviance卡方的p值<0.05时，拒绝"各水平频率无差异"的原假设。

方法二：似然比检验

通过对比零模型（假设各水平频率相等）和备选模型（允许各水平频率不同）的拟合度，判断差异是否显著：

# 零模型：仅包含聚类随机效应，假设各水平概率一致
null_model <- glmer(var ~ 1 + (1 | clusters), 
                    data = dat, 
                    family = binomial(),
                    control = glmerControl(optimizer = "bobyqa"))

# 备选模型：加入var的固定效应，允许各水平概率不同
alt_model <- glmer(var ~ var + (1 | clusters), 
                   data = dat, 
                   family = binomial(),
                   control = glmerControl(optimizer = "bobyqa"))

# 执行似然比检验
anova(null_model, alt_model)

检验结果的p值<0.05时，说明各水平频率存在统计学显著差异。

关键注意点 #

• 簇大小不均不影响混合模型的有效性，随机截距会自动调整不同簇的权重。
• 若分类变量水平数过多，brms包的多项混合模型稳定性更强，还能输出各水平的概率估计与置信区间。

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