本文案例数据是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。
参考资料
网友评论