决策曲线(Decision Curve Analysis)的绘制(R语言)

本文案例数据是NHLBI（美国国家心肺血液研究所）著名的Framingham心脏研究数据集的一个子集。大概长这个样子：

image.png

代码实现：

1.1 绘制决策曲线(Decision Curve)

#install.packages("rmda")
library(rmda)

Data<-read.csv('2.20.Framingham.csv',sep = ',')

a.首先构建一个简单模型

simple<- decision_curve(chdfate~scl,data= Data,
                        family = binomial(link ='logit'),
                        thresholds= seq(0,1, by = 0.01),
                        confidence.intervals = 0.95,
                        study.design = 'case-control',
                        population.prevalence = 0.3)

## Warning in decision_curve(chdfate ~ scl, data = Data, family = binomial(link =
## "logit"), : 33 observation(s) with missing data removed
## Calculating net benefit curves for case-control data. All calculations are done conditional on the outcome prevalence provided.

decision_curve()函数中，threshold设置横坐标阈概率的范围，一般是0-1；但如果有某种具体情况，大家一致认为阈概率达到某个值以上，比如40%，则必须采取干预措施，那么0.4以后的研究就没什么意义了，可以设为0-0.4。by是指每隔多少距离计算一个数据点。Study.design可设置研究类型，是“cohort”还是“case-control”，当研究类型为“case-control”时，还应加上患病率population.prevalance参数，因为在“case-control”研究中无法计算患病率，需要事先提供。

b.再构建一个复杂logistics回归模型complex

complex<-decision_curve(chdfate~scl+sbp+dbp+age+bmi+sex,
                        data = Data,family = binomial(link ='logit'),
                        thresholds = seq(0,1, by = 0.01),
                        confidence.intervals= 0.95,
                        study.design = 'case-control',
                        population.prevalence= 0.3)

## Warning in decision_curve(chdfate ~ scl + sbp + dbp + age + bmi + sex, data =
## Data, : 41 observation(s) with missing data removed
## Calculating net benefit curves for case-control data. All calculations are done conditional on the outcome prevalence provided.

## Note:  The data provided is used to both fit a prediction model and to estimate the respective decision curve. This may cause bias in decision curve estimates leading to over-confidence in model performance.

c.把simple和complex两个模型合成一个list

List<- list(simple,complex)

plot_decision_curve(List,
                    curve.names=c('simple','complex'),
                    cost.benefit.axis =FALSE,col= c('red','blue'),
                    confidence.intervals=FALSE,
                    standardize = FALSE)

## Note: When multiple decision curves are plotted, decision curves for 'All' are calculated using the prevalence from the first DecisionCurve object in the list provided.

image.png

图表解读：从上图可见在阈值在0.1~0.5大致范围内，complex模型的净受益率都比simple模型高。

plot_decision_curve() 函数的对象就是前面定义的List，如果只画一条曲线，直接把List替换成simple或complex即可。curve.names是图例上每条曲线的名字，书写顺序要跟上面合成list时一致。cost.benefit.axis是另外附加的一条横坐标轴，损失收益比，默认值是TRUE。col设置颜色。confidence.intervals设置是否画出曲线的置信区间，standardize设置是否对净受益率（NB）使用患病率进行校正。

# summary(complex,measure= 'NB') #结果很冗长！

Note: 查看complex模型曲线上的各数据点。NB也可以改成sNB，表示经过患病率的标准化。

1.2 绘制临床影响曲线(Clinical Impact Curve)

使用simple模型预测1000人的风险分层，显示“损失:受益”坐标轴，赋以8个刻度，显示置信区间

plot_clinical_impact(simple,population.size= 1000,
                     cost.benefit.axis = T,
                     n.cost.benefits= 8,
                     col =c('red','blue'),
                     confidence.intervals= T,
                     ylim=c(0,1000),
                     legend.position="topright")

image.png

使用complex模型预测1000人的风险分层，显示“损失:受益”坐标轴，赋以8个刻度，显示置信区间

plot_clinical_impact(complex,population.size= 1000,
                     cost.benefit.axis = T,
                     n.cost.benefits= 8,col =c('red','blue'),
                     confidence.intervals=T,
                     ylim=c(0,1000),
                     legend.position="topright")

image.png

图表解读：红色曲线（Number high risk）表示，在各个阈概率下，被simple模型（图3.）或complex模型（图4.）划分为阳性（高风险）的人数；蓝色曲线（Number high risk with outcome）为各个阈概率下真阳性的人数。

DCA算法原理解析

它相当于在回归预测分析的基础上，引入了损失函数。先简单定义几个概念：

P：给真阳性患者施加干预的受益值（比如用某生化指标预测某患者有癌症，实际也有，予以活检，达到了确诊的目的）；
L：给假阳性患者施加干预的损失值（比如预测有癌症，给做了活检，原来只是个增生，白白受了一刀）；
Pi：患者i有癌症的概率，当Pi > Pt时为阳性，给予干预。

所以较为合理的干预的时机是，当且仅当Pi × P >(1 – Pi) × L，即预期的受益高于预期的损失。推导一下可得，Pi > L / ( P + L )即为合理的干预时机，于是把L / ( P + L )定义为Pi的阈值，即Pt。

参考资料

决策曲线(Decision Curve Analysis)的绘制(R语言)