预测模型样本量-R实战

原创聂志强实战医学统计2022-03-07 13:02:06

本章思维导图：

前言

笔者一直以来认为临床预测模型样本量不是采用10*EPV(event per variable)，就是ROC95%CI计算，当检索依据时发现之前理解太过片面。BMJ文献(PMID: 32188600)提供了样本量范例，计算十分严谨。

结论如下：1.构建预测模型最好将所有可用数据用于模型开发（即避免数据拆分），并使用重采样方法（例如bootstrap自抽样）进行内部验证；2.在为二分类（binary）或生存数据（time-to-event）结局开发预测模型时，所需样本量计算方法有①经验法10*EPV；② 10*PPV，10倍候选预测变量；③研究人群总数、事件率、模型预期性能；3.四步法计算预测模型样本量，推荐R包pmsampsize。

正文

10EPV → 10EPP：开发预测模型需要一个训练数据集（development dataset），训练集的样本量必须足够大才能开发出在应用于目标人群中的新个体时的可靠的预测模型。经验法10EPV常常被误解，真正的EPV是确保对于考虑包含在预测模型方程中的每个预测变量至少有10个结局事件（即回归方程中的每个 β 项对应10个Y=1阳性事件，β0不算，Y=1为占比<Y=0的阳性事件）。一般来说，Y=β0+β1X1…βnXn， β1~βn对应n个预测变量。但是，如果某些预测变量是无序多分类转换为哑变量，例如血型ABO，此时三水平血型变量需要2个哑变量β，即β(A vs. O) 和 β(A vs. B)。另外，如果某些预测变量是连续变量根据Q4等非线性转化可能需要2~3个β来拟合非线性效应，例如age、血压。最后交互作用也会增加参数β个数。目前关于EPV的倍数还存在一些争议，不同的模拟研究结果并不一致，例如Vittinghoff 和McCulloch认为在某些情况下10倍EPV过于严格甚至可以去到5倍，Van Smeden等人却认为目前支持logistic回归EPV规则的证据是薄弱的。

因此，BMJ提倡最好应用每候选预测参数事件数（events per candidatepredictor parameter，EPP）来替代EPV。Candidate这个词很重要，因为模型过度拟合的程度受所考虑的参数总量影响，而不仅仅是包含在最终的final模型方程中的参数数量。然而尽管如此，10 EPP 规则的充分性也经常受到争论，有的意见是5EPP，有的是15~50EPP。EPP 实际上是具有场景特异性的，它不仅与候选预测参数相关事件的数量有关，也与参与人群总数、研究人群中结局事件的比例（发生率）以及模型的预期预测性能有关。所以，简而言之，10EPP规则使用需要经过更严谨的4步法来实现，4步中每一步都会进行样本量计算，其中最大的样本量就是所需的样本量。也就是N= Max (Nstep1,2,3,4)。4步法涉及结局为二分类、时间相关结局、连续变量结局。

Step1：实现准确估计结局整体风险或平均结果值需要的样本量？可以理解为率或均数或或OR/RR/HR 95%CI。

案例1：半数个体会发生的二分类型结局事件，需要至少385例的样本量才能将总体结局事件比例的置信区间限制到0.45到0.55，此时在真实值0.5附近的误差最多为0.05。为了在结局事件比例为0.1和0.2的情况下得到相同的误差幅度，分别需要至少139和246例样本量。二分类结局样本量公式如下：

对于生存数据类型结局事件，需要确定一个关键时间点以及预期的结局事件率。例如，在整个随访中每 100 人年的预期事件发生率为 10 次，如果10 年的结局事件概率为 0.63，为确保估计误差≤0.05（置信区间为0.58到0.68），样本量必须包括总共2366人年的随访。对于连续型结局事件，必须预先指定结局值的预期均值和方差，以及模型解释的预期变异百分比。

Step2：使所有个体的预测值具有较小平均误差需要的样本量？Van Smeden等人研究发现预测变量数目、总样本量和结局事件比例是导致模型平均预测准确性的三个主要因素。平均绝对预测误差（MAPE）表示预测模型对目标人群的结局概率的预测误差。

案例2：对于结局事件发生率为0.3的二分类结局，预测变量数为10，想要实现观测和真实结局概率之间的平均绝对误差为0.05样本量至少为461例，每个预测变量对应的结局数（EPP）为13.8。二分类结局样本量公式如下：

对于连续型结局事件，覆盖预测值范围的准确预测需要精确估计残差的标准偏差。要在估计的残差标准偏差中实现小于 10% 的乘法误差（multiplicative error），所需的样本量就是 234+P，其中 P 是所考虑的预测参数的数量。

Step3:最小化过度拟合所需的样本量。Riley等人建议在开发模型时注意在较小收缩(shrinkage，又称惩罚或正则化)率（≤0.1，即预期收缩因子S≥0.9）下确定样本量和预测变量的数量。例如对于二分类结局，R²cs为评价模型性能的保守指标Cox-Snell R²，它反映模型信噪比(signal/noise)，可影响多参数估计和模型潜在过拟合（接近0时，考虑参数估计不可靠出现过拟合；接近1时，反之欠拟合）。二分类结局R²cs 范围是0~max(R²cs)，连续结局R²cs 范围是0~1，R²cs设定需要参考外部文献证据和专家意见。因此，即使对于潜在预测性能优良的模型，预期的R²cs也可能很小。二分类结局样本量公式如下：

案例3：结局事件发生的比例为 0.3，在开发具有多达20个候选预测参数和至少0.1的预期R²cs的logistic回归模型时，需要1698例样本量以确保预期收缩为10% (S=0.9)，此时对应的EPP为25.5。所需的样本量和EPP对R²cs的选择很敏感，R²cs的预期值越低，所需的样本量越大。故建议保守选择R²cs！

Step4：保证开发模型和R²cs的乐观优化调整值之间差异较小？听起来非常拗口，其实乐观优化调整值就是R²cs加了个限制，R²_Nagelkerke=R²cs/max(R²cs)。乐观调整的R²_Nagelkerke值是模型在目标人群中拟合的更现实（近似无偏）估计。二分类结局样本量公式如下，对于二分类结局，建议δ≤0.05：

案例4:在开发具有预期R²cs为0.2的logistic回归模型时，并且在结局事件比例为0.05（因此max(R²cs) 为0.33），需要1079 例样本量才能确保乐观R²_Nagelkerke的预期乐观度仅为0.05。计算过程如下，结局事件发生率为0.05，假设结局事件发生的期望值E=5，n=100，对于Logistic模型，结局事件发生率分别为0.5、0.4、0.3、0.2、0.1、0.05和0.01时，对应的max(R²cs)分别为0.75、0.74、0.71、0.63、0.48、0.33和0.11。

综上所述，BMJ文中对三大类outcome进行了详细4步法样本量阐述并提供了细分代码，分别是二分类结局B1~B4, 时间关联结局T1~T3，连续结局C1~C4，并在pmsampsize包中对应了4step代码。

PASS、R 实战

案例1，率的95%CI估计，文献计算为简单估计N=385,更严格应该是exact法。其他OR/RR/均数区间均类似。

pmsampsize包可以针对预测模型进行4步法的样本量计算，取最大值为最终所需样本量。

案例2，shiny包制作网页工具https://mvansmeden.shinyapps.io/BeyondEPV，推荐参数：n为模型开发数据集样本量；φ为预期结局事件比例（≤0.5）；P为预测变量数目（≤30）。注意shiny工具用的是rMPSE (square root of the mean squared predictionerror)，原文用的是MAPE (mean absoluteprediction error)。案例2中MAPE设为0.05，网页版近似结果需要设置rMAPE0.07。

分类结局B1~B4, 时间关联结局T1~T3，连续结局C1~C4，并在pmsampsize包中对应了4step代码如下，方法学部分step3-4计算复杂直接用pmsampsize是否方便，且一次性给出4step过程所需样本量！

注意pmsampsize包输出的criteria1~3，分别对应的是B3、B4、B1! B1固定绝对单侧宽度±0.05，也就是事件率0.05±0.05，明显不合理，故还是推荐PASS计算B1样本量。

eg1 二分类型结局事件: North 等人开发了一个根据妊娠15 周时测量的临床指标预测孕妇发生子痫前期风险的模型。这些临床指标包括阴道流血、年龄、既往流产史、家族史、吸烟和饮酒。该模型包括 13 个预测参数，C-index为0.71。本例中的新研究旨在通过纳入额外的预测因子（例如，生物标志物和超声测量）来改进该模型或其他预测子痫前期的模型。由于二分类型结局，新预测模型的样本量计算需要包括标准B1到B4。输入参数：女性中发展为子痫前期的总体比例(0.05)、候选预测参数的数目（在本例中假设新增预测因子到30）。对于0.05 的结果比例，max(R²cs)值为 0.33。如果我们保守地假设新模型能解释15%的变异性，则预期的R²cs值为0.15×0.33=0.05。现在我们可以检查标准B1、B3和B4：

library(pmsampsize)
pmsampsize( 
  type="b", #b对应binary二分类型
  rsquared = 0.05, #预期的R2cs值
  parameters = 30, #候选预测参数的数目
  prevalence = 0.05, #结局事件发生的总体比例
  seed = 123456 #设置的随机种子
  )

标准B3而言，至少需要5249例女性，对应263个结局事件和8.75的 EPP。这个样本量才足以确保所需的预期收缩率仅为10%（以最大限度地减少过度拟合的可能）。要检查标准B2，需要应用案例2公式或Shiny，结果表明需要544例女性才能将预测概率的平均绝对误差控制为≤0.05。这远低于标准B3所需的5249例女性的要求。如果招募5249例女性时不切实际的话（由于时间、成本或数据收集的实际限制等原因），可以通过减少候选预测因子的数目来减少所需的样本量。例如，如果是20个而不是30个候选预测变量，则满足所有四个标准的样本量要求为3500例女性和175个结局事件（仍为8.75的EPP）。pmsamplesize也可以利用C-indx来计算，本来cstatistic = 0.71，计算B3结果为9493例，更加可怕。

综上所述pmsamplesize计算预测模型样本量不光有统计理论依据，还比10*EPV更省事件数，比如本例B3计算263例事件<< 10EPV的300例。

eg2 生存数据类型结局事件:许多预后模型可预测首次复发性静脉血栓栓塞 (VTE)治疗停止后复发的风险。例如，Ensor 等人的模型包括了年龄、性别、第一凝块部位、D-二聚体水平，以及从停止治疗到测量D-二聚体的滞后时间（通常约为30天）。模型的C-index为0.69，调整后的R²cs为0.051（对应于8%的总变异）。本例中的新研究旨在通过包含额外的预测变量来改进此类模型。新模型所需的样本量必须至少满足标准T1到T3。这要求我们输入预测 VTE 复发风险的关键时间点（例如两年），以及候选预测参数的数量（n= 30）、预期的平均随访时长（2.07 年）和结局事件的发生比例（0.065，或每1000人年随访发生65次VTE复发），以及 R²cs的保守值 (0.051)。这些设置的值均基于Ensor等人发表的VTE研究。

library(pmsampsize)
pmsampsize(
    type = "s", #s对应survival生存数据类型
    rsquared = 0.051, #预期的R2cs值
    parameters = 30, #候选预测参数的数目
    rate = 0.065,#结局事件发生的总体比例
    timepoint = 2, #进行预测的关键时间点
    meanfup = 2.07, #预期的平均随访时长，请注意timepoint和meanfup的单位要一致
    seed = 123456#设置的随机种子
  )

结果表明至少需要5143例样本，对应692个结局事件和23.1的EPP。这远远超过10 EPP，并且由仅10%的收缩因子（标准T2）产生，能解释仅8%的总变异。如果候选预测参数的数量减少到20，则所需的样本量减少到3429（EPP仍为23.1）。

eg3 连续型结局事件: Hudda等人开发了预测4至15岁儿童和青少年的无脂肪体重的模型。模型包括了基于身高、体重、年龄、性别和种族的10个预测参数。该模型可以估计具体个体的脂肪量（体重减去预测的无脂肪质量）。在外部验证中，该模型的R2cs为0.90。让我们假设该模型需要更新（例如，由于种群行为和环境在10年内的变化），并且需要在新模型中考虑额外的 10 个预测参数（因此总共20个参数）。模型的训练集的样本量大小必须满足方框1中C1到C4四个标准。这要求我们指定预期的R2cs (0.90)、候选预测参数的数量 (n=20) 、总体结局的平均值 (26.7 kg)、目标人群中无脂肪质量的标准差 (8.7 kg)。

library(pmsampsize)
pmsampsize(
  type = "c",  #c对应continuous连续型
  rsquared = 0.9, #预期的R2cs值
  parameters = 20, #候选预测参数的数目
  intercept = 26.7, #总体结局的平均值
  sd = 8.7,#总体结局的标准差
  seed = 123456 #设置的随机种子
)

需要254例样本，因此每个预测参数需要12.7例样本。254的样本量取决于精确估计模型标准偏差（标准 C3）所需的数量，因为只需要68例样本即可最大限度地减少过度拟合（标准C1和C2）。

发散思维

临床预测模型按照动态情况分为：1）短期模型：如临床诊断模型，根据研究对象的临床症状和特征，诊断当前患有某种疾病的概率，常常采用横断面设计；2）长期模型：临床预后模型，在当前的疾病状态下，未来某段时间内疾病复发、死亡、伤残以及出现并发症等结局的概率，常常采用RCT或队列；3）动态预测模型：随着时间推移，有关预测因素的信息变得可用，风险也会发生变化。临床预测模型原则上是不推荐用病例对照的，除非基于建立在队列本身的巢式病例对照或者病例队列设计。

常见的外部验证数据集有时间验证或狭义验证（ temporal validationdesign），比如是把队列中2010-2015年入组患者作为训练集，把2016-2020年入组患者作为外部训练集；地理或广泛验证（ spatial validation design），用一个中心/区域数据建模，取其他中心/区域的数据作为外部验证。作为常见的外验证样本量计算上，BMJ文章并没有给出答案。因为作者认为不建议将数据拆分为训练集和测试集，而应使用所有数据搭配bootstrap方法进行内部验证。对于预测模型交叉验证问题，众说纷纭没有统一结论，但是笔者认为作者原意是和title对应的，仅做开发模型即development data，故全文并不适合外验证样本量计算。那是否意味外验证样本量又要回归到10EPV？

TRIPOD申明(模型报告质量)指出预测模型的研究中必须说明样本量是如何确定的，样本数量能体现出预测变量与结局之间的关联是否稳固、预测模型拟合是否合适、模型是否存在潜在过度拟合；样本量是基于统计方法，还是基于现实考虑（比如使用现有数据）。PROBAST申明(模型建立质量)则推荐高质量预测模型研究构建模型需要EPV或EPP≥20，验证模型需要阳性事件≥100。

基于机器学习方法建立的预测模型，样本量并没有严格的定义。一般，机器学习技术可能需要为每个预测变量增加10倍以上的结局事件来确保过度拟合的程度较小，并且即使超过200 EPP也可能表现出不稳定和高度的乐观值。因此，机器学习方法对样本大小的要求应该更高。大多数医学研究数据集的大小更适合使用回归方法建模（包括惩罚和收缩方法），特别是因为回归建模方法产生的模型方程更具有可解释性，便于实施、验证和图形展示，而不是以拿到数据上来就跑各种高阶ML模型。

小结

最近几年临床预测模型各种刷屏，其中有很多样本量很少却结局罕见的研究。此类文章会因为样本量不足导致预测不准确，甚者可能会误导参考临床预测模型做临床决策的医生。然后现实情况是，样本量问题依旧就像缺失值问题一样经常被研究者选择性忽略。笔者越深入挖掘预测模型越发觉，这个科学算命其实存在很多玄学基因在里面，比如随机种子不同模型效果不同，超参数无法穷举等等，甚至高端复杂ML方法不一样超过简单logistic，但是噱头吹逼的时代需要包装成高大上才显得高级。既然建模方法都存在百花齐放，那样本量其实自然就可以自我标榜了。最后笔者小结：4步法统计算法分别是95%CI、MAPE、shrinkage、乐观R²CS，总阳性数可能<10*EPV/EPP 也可能 >>10*EPV/EPP。

参考文献

1.Riley Richard D,Ensor Joie,Snell Kym I E et al. Calculating the sample sizerequired for developing a clinical prediction model.[J] .BMJ, 2020, 368: m441.

2.Vittinghoff Eric,McCulloch Charles E,Relaxing the rule of ten events pervariable in logistic and Cox regression.[J] .Am J Epidemiol, 2007, 165: 710-8.

3.van Smeden Maarten,de Groot Joris A H,Moons Karel G M et al. No rationale for 1variable per 10 events criterion for binary logistic regression analysis.[J].BMC Med Res Methodol, 2016, 16: 163.

4.Riley RD, Snell KIE, Ensor J, Burke DL, Harrell FE, Jr., Moons KG, Collins GS.Minimum sample size required for developing a multivariable prediction model:Part I continuous outcomes. Statistics in Medicine. 2018 (in-press). doi:10.1002/sim.7993.

5.Moons Karel G M,Wolff Robert F,Riley Richard D et al. PROBAST: A Tool to AssessRisk of Bias and Applicability of Prediction Model Studies: Explanation andElaboration.[J] .Ann Intern Med, 2019, 170: W1-W33.

6.van Smeden Maarten,Moons Karel Gm,de Groot Joris Ah et al. Sample size forbinary logistic prediction models: Beyond events per variable criteria.[J].Stat Methods Med Res, 2019, 28: 2455-2474.

原创不易，欢迎转载，转载请标明来自本公众号。