大家好，我是风间琉璃。最近临床预测模型好像特别火，通过构建临床预测模型对某一个疾病的发生或者是某个疾病并发症的发生进行预测，并对构建相关的评分量表，比如nomogram【列线图，可个性化计算特定肿瘤患者生存率】。根据疾病相关的危险因素作为评分指标，最后将预测模型带入到单独的个体计算出总分，从而评估其患病的风险。 当做到这一步时候，我想问一个问题，那就是对于我们所构建的这个预测模型，怎么进行评价呢？那么今天我给大家介绍怎么对模型进行评价。因为具有较稳定的预测能力和推断自变量权重的特点，logistic在医学论文中广泛应用，那么本次就以logistic方法为例。来对模型进行评价吧。
首先我为大家介绍一下本次使用的数据.这次我选择的疾病是高血压数据。我们的数据一共有ID、Gender、Age、A、hypertension五个变量。该数据的目的是通过对每一个观测（每一行）的性别、年龄和指标A（这个指标可以是分子标志物、临床上的任意指标），以及是否发生高血压疾病进行收集。构建预测模型。明确哪一类人是我们的潜在高血压疾病患者。从而将更多的健康宣传、金钱投入向这一类人群进行更多倾斜。
那么这儿有一个问题，那就是这个模型是否能够准确推断出我们预测的那一类人就是真正会得高血压的患者呢？这个问题推广一下，到医学领域也就是，我们构建的临床模型推测的潜在患病人群是否真的会患病呢？这就是模型的评价。那我们开始吧。第一步，我们需要导入数据 ##setwd(“D:/ROC and Calibration”) ##并且不要ID这一列 dataset = read.csv('disease.csv')
dataset = dataset[2:5]
第二步，在构建模型之前需要处理数据 ##将结局变量转换为因子型变量 ##将数据随机分为训练集和测试集，训练集的作用就是构建模型，测试集的目的就是构建的模型进行评价，我们采取7：3的比例随机将数据进行拆分。###install.packages(“CaTools”) dataset$hypertension = factor(dataset$hypertension, levels = c(0, 1))
library(caTools)
set.seed(123)
split = sample.split(dataset$hypertension, SplitRatio = 0.7)
training_set = subset(dataset, split == TRUE)
test_set = subset(dataset, split == FALSE)
我们来看一看新的数据集的结构。 dim(training_set)
## [1] 280 4
str(training_set)
## 'data.frame': 280 obs. of 4 variables:
## $ Gender : Factor w/ 2 levels "Female","Male": 2 1 2 1 1 1 1 2 2 2 ...
## $ Age : int 19 26 27 27 32 35 26 20 18 29 ...
## $ A : int 19000 43000 58000 84000 150000 65000 80000 86000 82000 80000 ...
## $ hypertension: Factor w/ 2 levels "0","1": 1 1 1 1 2 1 1 1 1 1 ...
好了，现在就得到我们想要的数据模型了。 有两个函数可以构建模型，分别是R本身的stats包里面的glm函数和rms包的lrm函数。我们这儿是使用的glm函数。 classifier = glm(formula = hypertension ~ .,
data = training_set,family ="binomial" )
summary(classifier)
##
## Call:
## glm(formula = hypertension ~ ., family = "binomial", data = training_set)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0629 -0.4945 -0.1071 0.2942 2.5220
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.383e+01 1.811e+00 -7.638 2.20e-14 ***
## GenderMale 2.660e-01 3.792e-01 0.701 0.483
## Age 2.599e-01 3.527e-02 7.368 1.73e-13 ***
## A 3.833e-05 6.833e-06 5.610 2.03e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 364.98 on 279 degrees of freedom
## Residual deviance: 181.48 on 276 degrees of freedom
## AIC: 189.48
##
## Number of Fisher Scoring iterations: 6
我们发现除了性别，年龄和薪资都是结局变量的重要预测因素（P<0.005）。 接下啦我们即将开始对模型进行评价，需要注意的是，我们对模型的评价应该要做三方面，分别是区分度（discrimination）、校正度（Calibration） 和决策曲线分析（ decision curve analysis， DCA）。具体的理论知识我简单介绍一下。区分度就是模型能不能区分是或者否发生某个结局时间。 我们通常以ROC曲线下面积AUC衡量区分能力。0.5-0.7代表模型的区分能力差，这个模型不太适用。0.7-0.9这个模型的区分能力较好，已经能够应用在医学领域。而大于0.9则表明这个模型的区分效果很好。而校正度（calibration）则是，你得出的预测概率是否和真实的结果一致性评价。我们一般用Hosmer-Lemeshow (H-L) test评价，所得统计量卡方值越小、对应P值越大校准度越好。最后的DCA就模型进行预测时是否真正的能够净获益。我们首先进行区分度的检验，而模型的区分能力，具有两个评判指标，分别是ROC曲线下面积和C指数（C-index），对于logistic回归分析来说，我们使用ROC曲线下面积（AUC）这一个指标就可以了。 插上一句，在进行模型的评价时，我们如果用训练集training_set来评价模型显然并不合理，这未免就有王婆卖瓜自卖自夸的嫌疑了。所以，就需要我们的test_set了。首先我们使用训练集所得到的模型对测试集的每一个观测进行预测。得出每一个观测出现结局事件的预测概率，我们展现前面的十个预测值的结果。 y_pred = predict(classifier, newdata = test_set[-4],type = "response")
y_pred[1:10]
## 2 4 5 9 12 14
## 0.023999387 0.009636432 0.003276809 0.002999543 0.006156426 0.010334937
## 18 19 20 22
## 0.293941724 0.368247383 0.438425694 0.564420975
我们可以看到前面十个的预测值的结果其实并不是返回的0或者1，而是返回一个个小数，这是什么呢？其实这代表通过模型进行预测后，预测结局变量为阳性结果的可能性。接下来我们进行ROC曲线的计算，这里我们使用时ROC的pROC包。没有安装的的同学记得安装一下。 ##install.packages(“pROC”)
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
ROC=roc(test_set[,4],y_pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
计算ROC曲线下面积和最佳阈值，并绘图。 auc(ROC)
## Area under the curve: 0.9103
plot(ROC,print.thres=T,print.thres.best.method="youden")
这样就绘制出了我们的ROC曲线，得到曲线下面积是0.9103，根据我们之前所说的这个模型很棒了对不对。并且我们打印出其最佳的阈值，也就是0.2940，其特异度和灵敏度分别是0.844，0.907。这里的最佳阈值即是约登指数最高的值，等于也就是灵敏度与特异度之和减去1。也就是曲线最左上角的那个点。接下来我们稍微美化一下图片 plot(smooth(ROC),col="red",print.auc=T,legacy.axes=T)
legend("bottomright",legend = c("smoothed"),col = "red",lwd = 2)
这样是不是好看多了呢？
|
