CI上限直奔250，仍可发北大核心？——小狐解密 附R代码

阳性样本不足10人，多因素建模OR置信区间上限直奔250，居然也能发北大核心（《中国组织工程研究 》为中国康复医学会主办的北大核心期刊）？

——自学统计的艺术史小狐（用的对象账号发帖）为您解密！

本文内含相关问题的解决方案及R代码教学示范。欢迎点赞、收藏。

前一段时间，本人在本论坛吐槽了一篇中文期刊文章，阳性样本不足10人还构建了多因素的回归模型。

一些网友发现自己构建模型也存在相似的问题，给我发了私信进行咨询。

有类似问题的朋友，可观察Cox或逻辑回归模型结果，主要查看各个变量的风险比HR/OR及置信区间是否正常。如果没有大于10或者数据特别小（0.0几），那就说明数据是OK的。可以正常分析，按常规的套路走完流程就好。

近期我系统整理了对这篇文章的意见，形成了给编辑部的投诉信。

详细内容如下。如果回归模型系数不正常，或没有外部验证集，可在文章中找到相关问题的解决方案及参考文献。

您好，我是一名毕业于XX大学艺术与设计学院的研究生，也是贵刊的读者。近日在细致阅读您刊2020年第24卷第18期的文章《全膝关节置换治疗类风湿性膝关节炎的早期急性并发症》（doi: 10.3969/j.issn.2095-4344.2648，时为XX大学附属医院在读研究生XX为第一作者）时，我对这篇文章临床模型的构建产生了很大的疑问。

根据原文的表2，本研究共纳入了300例全膝关节置换术患者，但只有不到10例出现了术后30天内的并发症。数据中阳性结局占比过低，并不满足建立二分类结局的临床预测模型的“10EPV原则”，即每个预测变量至少对应10个阳性结局，否则会导致混杂因素的影响过大。作者最终筛选了3个独立预测因素，那么至少需要30个阳性结局才能建立有效的模型。此外，全膝关节置换术后30天内发生并发症在临床上并不算是罕见事件。因此，利用阳性样本极少的数据建立的逻辑回归模型，不仅存在统计学上的缺陷，也缺乏临床实践的指导意义。作者没有在研究开始前进行合理的样本量估算，并在论文中报告相关的过程和结果，说明其对统计学基本逻辑概念的认识不够系统。

原文的表3显示，有脑血管疾病史的患者相比于无脑血管疾病史的患者，发生并发症的风险高达30.736倍（31.736-1=30.736）。我咨询了流行病学专家，他们认为，当逻辑回归的OR或Cox回归的HR大于3时，就很难在临床上进行合理解释。从统计学角度来看，这可能是由于结局过于不均衡，导致了脑血管疾病对结局影响的高估。同时这样的模型也可能低估“假阳性率”。换言之，这样的模型容易产生“第二类错误”，极易造成“误诊”。

此外，“脑血管疾病”变量在多因素逻辑回归中的OR的95%置信区间（95% CI）过于宽泛（4.053-248.517），CI的上限更是直奔“二百五”（ 能得出250如此之大的数据，说明作者本人也🤦‍，嘿嘿嘿），这说明模型的估计误差过大。作者应采用一种稳健的误差方差估计法——夹心方差估算量（sandwich variance estimators）对95% CI进行矫正，以让95% CI回归于一个相对更合理的范围，并根据上述稳健估计的95% CI重新计算p值。实现这一步骤的R代码请原文作者参考文末附录。

变量在模型中的系数过高，亦反映出模型存在“过拟合”的严重问题。原文的图2A和图2B显示，“年龄”和“病程”两个变量各自单因素逻辑回归的ROC曲线下面积不到0.7，但“脑血管疾病”的ROC曲线下面积达到了0.865，趋近于0.9，进一步说明了“脑血管疾病”单因素逻辑回归模型及多因素逻辑回归存在过拟合的问题。过拟合的模型在训练数据上的表现可能非常好，但在未见过的数据（测试集或验证集）上的表现通常较差，即表现出较低的泛化能力。因此，在建模之后，作者应当充分验证模型。具体来说，有以下几点。

首先，作者仅提供未经重抽样算法内部验证的逻辑回归单因素的ROC曲线及ROC曲线下面积（原文图2），根本不足以证明模型具有良好的区分度。作者应进一步提供多因素逻辑回归模型的ROC曲线可视化及曲线下面积数据。其次，作者还应提供精准率-召回率（PR）曲线可视化及PR曲线线下面积数据，以证明模型未低估假阳性。第三，作者在绘制多因素回归模型的ROC曲线/PR曲线及计算曲线下面积时，应采用重抽样算法，如1000次bootstrap自助法，以进行有效的内部验证，并提高模型的证据等级。

有关于临床预测模型开发与验证的相关学术规范，可详参2015年发表于《BMJ》期刊上的文章《个体预后与诊断预测模型研究报告规范——TRIPOD声明 [Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): the TRIPOD statement]》（10.1136/bmj.g7594）。

若作者仍坚持在不进一步搜集阳性结局受试者患病资料的基础上，采用机器学习构建多因素模型（逻辑回归属于广义上的机器学习），即便采用了类似于SMOTE-NC的过采样或欠采样数据均衡手段，仍极有可能造成模型过拟合。因此，作者应考虑使用特定的不均衡数据处理学习器，如不均衡随机森林，而不是常规的逻辑回归。不均衡随机森林能够通过对结局情况为占比较少的那组样本加大权重，从而在保证高准确率的前提下，尽可能识别出结局阳性者。也就是说，这样的预测模型能够有效避免误诊。

同时，我注意到，作者在讨论章节中提到了ASA（美国麻醉医师协会）麻醉等级，并指出这个变量对研究的结局有影响。那么，作者理应在模型中也考虑这一变量。忽略这样一个重要变量可能会导致模型存在偏倚，从而影响模型的准确性和可靠性。

……最后是结束语，就不展示了。

接下来是实现压缩置信区间的R代码

library(survival)

library(gtsummary)

#lung数据集为R包survival自带的肺癌数据集

lung=na.omit(lung)

table(lung$status)

# 1  2

# 47 120

47/nrow(lung)

#[1] 0.2814371

#阳性结局占比约三成

mod <- glm(status-1~sex, data = lung, family = binomial())

summary(mod)

# 使用gtsummary包输出模型的系数

# 将系数转换为OR，使用exp()函数即可

mod%>%

 tbl_regression(estimate_fun = function(x) style_number(x, digits = 3))

#未处理的置信区间：-1.684~-0.295

#未处理的p value：0.005

vcovHC(mod, type = "HC")

sandwich_se <- diag(vcovHC(mod, type = "HC"))^0.5

sandwich_se

#lower bound

coef(mod)+qnorm(0.975)*sandwich_se

coef(mod)+qnorm(0.025)*sandwich_se

#夹心估计的置信区间：-1.675~-0.291

#可以看到置信区间被压缩

#但压缩不多，因为测试数据结局阳性比例没有特别偏低，此代码仅作为演示

z_stat <- coef(mod)/sandwich_se

p_values <- pchisq(z_stat^2, 1, lower.tail=FALSE)

p_values

 # (Intercept)         sex

 # 2.112155e-05 5.368168e-03

#p值几乎没有变化，这是因为测试数据结局阳性比例没有特别偏低