Logistic回归模型组间两两比较：不同数据集营养状态差异分析求助

我希望分析不同数据集间**营养状态（Nutritional.Status）**的差异，为此构建了二项式广义线性模型：

model <- glm(Nutritional.Status ~ Data.origin, family = 'binomial'(link='logit'), data = data)
summary(model)

模型输出结果如下：

> model <- glm(Nutritional.Status ~ Data.origin, family = "binomial", data = data)
> summary(model)

Call:
glm(formula = Nutritional.Status ~ Data.origin, family = "binomial", 
    data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9667  -0.9469  -0.9469   1.4269   1.4269  

Coefficients:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -0.08701    0.41742  -0.208 0.834879    
Data.originIR.recent  1.86478    0.52121   3.578 0.000347 ***
Data.originUK        -0.48261    0.43043  -1.121 0.262185    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 686.54  on 498  degrees of freedom
Residual deviance: 614.61  on 496  degrees of freedom
  (428 observations deleted due to missingness)
AIC: 620.61

Number of Fisher Scoring iterations: 4

我需要进行数据集间的两两比较，判断是否某一数据集中的良好营养状态占比显著高于其他数据集。尝试用emmeans后得到的结果不是我想要的：

pigs.emm.s <- emmeans(model, "Nutritional.Status")
> pairs(pigs.emm.s)
 contrast                    estimate    SE  df z.ratio p.value
 Good - Moderate               -0.489 0.587 Inf  -0.834  0.6818
 Good - (Poor-very poor)        0.785 0.494 Inf   1.588  0.2510
 Moderate - (Poor-very poor)    1.274 0.644 Inf   1.980  0.1172

Results are given on the log odds ratio (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

请问：

1. 如何得到数据集间营养状态（尤其是良好状态占比）的两两比较结果？
2. 现有模型输出是否足以说明不同数据集的营养状态构成存在显著差异，进而判断某一数据集的良好状态占比更高？

解决方案 #

一、先明确模型基准组信息 #

从模型系数输出可知，模型默认以未标注的基准数据集（非IR.recent、UK的组）作为参照：

• Data.originIR.recent系数为1.86且P值<0.001，说明IR.recent组相对基准组，营养状态为“良好”的对数优势比显著更高；
• Data.originUK系数不显著，说明UK组与基准组的良好状态优势比无显著差异。

但这仅为单组与基准组的比较，要完成所有数据集间的两两对比，需调整emmeans的用法。

二、正确用emmeans实现数据集间两两比较 #

你之前的代码是在比较不同营养状态的差异，而非数据集差异。正确做法是针对Data.origin分组计算边际均值，再做两两比较：

# 按Data.origin分组，计算良好营养状态的概率（转换到响应尺度）
origin_emmeans <- emmeans(model, ~ Data.origin, type = "response")
# 执行两两比较，用Tukey方法校正多重比较误差
pairs(origin_emmeans, adjust = "tukey")

• type = "response"将结果从对数优势比转换为概率尺度，即你需要的良好营养状态占比；
• pairs()会输出所有数据集间良好状态占比的对比结果，同时校正多重比较带来的假阳性问题。

三、现有模型输出的结论 #

从模型偏差来看：零偏差686.54，残差偏差614.61，两者差异显著（可通过卡方检验验证：pchisq(686.54 - 614.61, 498 - 496, lower.tail = FALSE)，结果远小于0.05），说明加入Data.origin变量后模型拟合效果显著提升，即不同数据集的营养状态构成确实存在差异。

结合系数可知IR.recent组的良好状态占比显著高于基准组，但要确认IR.recent是否显著高于UK组，必须通过上述两两比较代码的输出结果判断。

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