模型的作用，往小了说就是发文章，往大了说就是给临床医生的诊疗提供决策参考，最终模型的理想出路就是回归临床，变成商业项目。
比如，乳腺癌的21基因检测，70基因检测。
一开始，每个患者临床信息就那么多，20来个吧，只够构建一种预测模型。
后来，每个患者的样本可以测序了，
测了RNAseq，每个患者就多了2万个特征，这时候根据前期删选变量的依据不同，可以作出各种模型。
测了miRNA，又来一波模型
再测甲基化，还能来一波。
RNAseq和miRNA叠加，又来。
都做了一遍，怎么办？
用RNAseq的数据推测免疫浸润，把每个患者的免疫浸润量化一遍，又能预测模型了。
自噬火了，把自噬的基因挑选出来做模型
m6A火了，把他们家族成员拿出来也可以做。
反正，没有穷尽。
做模型，就两个条件，第一，你得有分类指标，比如化疗敏感或者耐药，生存周期，第二，你得有变量，越多越好。
模型发文虽好，但是也有大坑，尤其当一些临床都没去过的人来构建模型，坑就更多了。
比如，在训练集中产生的cutoff值，跟测试集中的cutoff值，应该是一样的，而不是各自产生不同的cutoff，就好比，如果AFP大于400诊断肝癌，那么这个400就不应该跟着患者变化。
因为这个模型最终是要用在临床的，cutoff值在变，那么还有什么意义呢？
比如，在RNAseq数据中训练，去芯片数据中测试，两个数据产生的途径和原理都不一样，我觉得是不妥的，除非你的模型泛化能力极强。
又比如，模型构建的文章，最好能把原始数据展示出来，这样，所有人都可以重复，而很多文章不给数据。
那，模型构建还需要学习么？作为拓宽眼界的一种方法，值得学习。
技术不是万能，但是没有技术万万不能。
学习已经普及的技术，是为了让自己以后更好的创造。
比如，我看到有的临床医生，用机器学习提取图像特征，那么这个特征可以是几百个，然后是不是又可以做模型了呢，当然！文章都已经出来了。
回到今天包子的帖子，我十分推荐。
因为，写得很细致，事无巨细，生怕你学不会。
跟着教程就可以学习一个新的技能。
如果你还不知道如何去学习别人的帖子，这里有一个示范的视频可以参考。
果子小视频: 给初学者的诚意GEO数据分析教程。
事无巨细这种优秀的品质，会随着技术的提高慢慢消失，写的清楚是为了检测自己是否掌握，当你真的掌握了之后，脑子会自动过滤，写都不愿意写出来的。
请珍惜此刻的包子。
之前他也写了几篇质量很高的帖子，我放在了文末。
大家好，我是那个努力健身想瘦成馒头的包子，真的是好久不见啊，不知道大家疫情期间过得怎么样，我是被困在家整整两个月，虽然工作效率会差很多，但是还是坚持在学习一些东西，毕竟看到自己的前辈们奋斗在一线，虽然暂时帮不上什么忙，但是至少不能虚度自己的时光吧。而且最近又有了新的健身心得，自己试了之后感觉还不错，我再归纳总结以下后面跟大家分享，哈哈哈。
对了，今天要给大家介绍的这个Logistic回归模型真的让我很激动，因为这应该也是大家最感兴趣的一种算法。之前,我们已经介绍了简单线性回归模型，但是这种模型只能预测连续性的因变量，并且对临床医生而言接触最多而且最熟悉的还是要算Logistic回归模型。
因为我们在现实生活中，我们要预测的模型的因变量大多数都是二分类变量，例如要预测病人用药后的疗效（好/坏），预测肿瘤是否发生转移（是/否）等等，我们接下来介绍的Logistic回归就是可以用来预测定性结果的算法。
虽然计算机技术的革新和很多机器学习方法的崛起，可以用于预测定性结果的算法还有很多，例如决策树分类模型，随机森林分类模型和支持向量机等。但是在分类问题上Logistic回归模型仍然可以和其它算法平分秋色。
其实我们只需要弄懂下面三个部分，基本上就可以掌握Logistic回归的大部分应用：
1.什么是Logistic回归？
2.模型构建与评价
3.绘制ROC曲线图
因为只要是一本关于统计或者机器学校的书都会狠狠的介绍Logistic回归的概念，因此在这里我们没有必要在概念上花太多时间。Logistic回归的方程可以有很多种展示形式，这里我们就只介绍一种：
log[p/(1-p)] = b0 + b1*x1 + b2*x2 + ... + bn*xn
等式右边是不是看起来像极了我们熟悉的线性回归方程，这是因为Logistic回归就是以对数发生比为响应变量进行线性拟合，这里的系数是通过极大似然估计得到的，而不是我们之前熟悉的最小二乘法。同我们之前讲的简单线性回归方程一样，b0和b1是回归的β系数。当b1大于0时，表明x的增加会增加事件发生的概率p；相反,b1小于0时，表明增加x将与减少概率p。
等式左边的p表示事件发生的概率，log[p/(1-p)]表示odds的对数，odds（称为几率、比值、比数），是指某事件发生的可能性(概率)与不发生的可能性（概率）之比。因此log[p/(1-p)]的值表示为odd的对数，称为log-odd值或logit值。
例如，如果某个饮食习惯会导致糖尿病阳性的概率是0.8，那么“非糖尿病”的概率是1-0.8 = 0.2，而几率odd = p/(1-p) = 0.8/ 0.2 = 4,这个值就是临床上经常使用的OR（Odd Ratio）值。
下面终于要开始我们今天最令人兴奋的一个部分，我们要用Logistic回归来帮助我们解决真正的临床问题。我们使用的是MASS包中的biopsy数据，这个数据集收集了患者组织样本的检测数据，其中有11个变量分别为：
> library(MASS)
> data(biopsy)
> str(biopsy)
'data.frame': 699 obs. of 11 variables:
$ ID : chr '1000025' '1002945' '1015425' '1016277' ...
$ V1 : int 5 5 3 6 4 8 1 2 2 4 ...
$ V2 : int 1 4 1 8 1 10 1 1 1 2 ...
$ V3 : int 1 4 1 8 1 10 1 2 1 1 ...
$ V4 : int 1 5 1 1 3 8 1 1 1 1 ...
$ V5 : int 2 7 2 3 2 7 2 2 2 2 ...
$ V6 : int 1 10 2 4 1 10 10 1 1 1 ...
$ V7 : int 3 3 3 3 3 9 3 3 1 2 ...
$ V8 : int 1 2 1 7 1 7 1 1 1 1 ...
$ V9 : int 1 1 1 1 1 1 1 1 5 1 ...
$ class: Factor w/ 2 levels 'benign','malignant': 1 1 1 1 1 2 1 1 1 1 ...
ID:样本编码
V1-V9: 患者9个水平的细胞检测数据，医疗团队对其进行了评分和编码，因此全部为整数型变量，今天我们主要讨论的是Logistic回归预测分类的能力，因此我们不展开讨论每个变量的细节。
class: 肿瘤诊断的结果，良性或者恶性，肿瘤诊断的结果。
#去除缺失值
biopsy.dat <- na.omit(biopsy)
#删除第一列
biopsy.dat <- biopsy.dat[,-1]
#相关性检测
library(corrplot)
bc <- cor(biopsy.dat[,1:9])
corrplot.mixed(bc)
#训练集与测试集
set.seed(333)
ind <- sample(2,nrow(biopsy.dat),replace = TRUE,prob = c(0.7,0.3))
train <- biopsy.dat[ind == 1,]
test <- biopsy.dat[ind == 2,]
其实可以用函数createDataPartition(p = 0.7, list = FALSE)来快速进行分类，但是sample函数可以更好的帮助我们理解按照70/30比例分类的过程。所有准备工作结束，接下来我们就可以愉快的开始建模了。
其实构建模型的代码很简单：
#family = binomial这个参数表示我们使用的是Logistic回归。
full.fit <- glm(class~.,family = binomial,data = train)
使用summary函数查看预测变量的系数及P值：
summary(full.fit)
Call:
glm(formula = class ~ ., family = binomial, data = train)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.4125 -0.1137 -0.0647 0.0265 2.4561
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.76818 1.27651 -7.652 1.97e-14 ***
V1 0.44497 0.14624 3.043 0.002345 **
V2 -0.06525 0.23325 -0.280 0.779660
V3 0.30168 0.25944 1.163 0.244911
V4 0.30947 0.13039 2.373 0.017623 *
V5 0.01449 0.20706 0.070 0.944217
V6 0.39716 0.10775 3.686 0.000228 ***
V7 0.56501 0.19481 2.900 0.003729 **
V8 0.21455 0.12131 1.769 0.076954 .
V9 0.51266 0.30162 1.700 0.089187 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 617.608 on 482 degrees of freedom
Residual deviance: 79.509 on 473 degrees of freedom
AIC: 99.509
Number of Fisher Scoring iterations: 8
使用confint函数查看95%置信区间：
confint(full.fit)
2.5 % 97.5 %
(Intercept) -12.689180688 -7.5962990
V1 0.177112712 0.7595758
V2 -0.514940729 0.4109586
V3 -0.215593581 0.8164297
V4 0.051170515 0.5742843
V5 -0.386684463 0.4316285
V6 0.195255645 0.6224505
V7 0.204588354 0.9753586
V8 -0.017441016 0.4643443
V9 0.002451551 1.0714645
计算优势比，就是我们临床中常说的OR值，使用exp(coef())函数:
exp(coef(full.fit))
(Intercept) V1 V2 V3 V4 V5 V6 V7
5.724458e-05 1.560443e+00 9.368295e-01 1.352129e+00 1.362697e+00 1.014593e+00 1.487594e+00 1.759460e+00
V8 V9
1.239306e+00 1.669730e+00
通过上面的代码，我们似乎就得到了很多可以用来向别人展示的数据，但是这只是表面的现象，要挖掘这些数据背后真正的含义，我们要做到远远不止如此。
首先我们讨论模型潜在的多重共线性问题，我们之前的帖子已经讨论过，直接用vif()函数就可以解决：
#多重共线性问题
library(car)
vif(full.fit)
V1 V2 V3 V4 V5 V6 V7 V8 V9
1.239979 3.388806 3.339807 1.168101 1.480804 1.227726 1.233238 1.207534 1.099869
根据VIF经验法则，一般VIF值大于5为存在共线性可能，因此我们的数据不存在共线性问题，接下来要进行的就是模型验证。
首先我们查看一下模型在训练集和测试集的表现分别怎么样：
#训练集
library(glmnet)
y <- ifelse(biopsy.dat$class == 'malignant',1,0) #0表示良性，1表示恶性
train.probs <- predict(full.fit,type = 'response') #得到训练集预测概率
train.probs <- ifelse(train.probs>= 0.05,1,0) #预测概率大于0.05默认为肿瘤发生，同样设置为1
train_class <- y[ind ==1] #得到实际肿瘤发生情况
table(train_class,train.probs) #混淆矩阵比较实际发生与预测发生的关系
得到结果如下：
train.probs
train_class 0 1
0 297 23
1 1 162
预测正确率为：（162+297）/ 483 * 100% = 95%
#测试集
test.probs <- predict(full.fit,newdata = test,type = 'response') #注意这里要将数据改为测试集数据
test.probs <- ifelse(test.probs>= 0.05,1,0) #预测概率大于0.05默认为肿瘤发生，同样设置为1
test_class <- y[ind ==2] #得到实际肿瘤发生情况
table(test_class,test.probs) #混淆矩阵比较实际发生与预测发生的关系
得到结果如下：
test.probs
test_class 0 1
0 120 4
1 0 76
预测正确率为：（120+76）/ 200 * 100% = 98%
从结果来看模型的准确率还是很高的，但是我们在面临实际的临床问题的时候，结果经常差强人意，以上是我们通过所有变量得到的模型，模型显得十分的厚重，并且并不是纳入变量越多模型越可靠，并且如果我们想要构建自己的预测模型时，太多的预测变量纳入不利于模型的推广，因此我们要想办法找出真正对模型预测其最大作用的相关变量。
这里我们使用bestglm包进行特征变量的简单筛选，并且使用BIC的最优子集方法：
1.赤池信息准则（Akaike Information Criterion，AIC）：AIC是衡量统计模型拟合优良性的一种标准，由日本统计学家赤池弘次在1974年提出。AIC的基本思想是惩罚包含额外变量的模型，因此我们只要记住AIC越低，模型越好。
2.贝叶斯信息准则（Bayesian Information Criterion，BIC）:BIC贝叶斯信息准则与AIC相似，用于模型选择，1978年由Schwarz提出。BIC的惩罚项比AIC的大，考虑了样本数量，样本数量过多时，可有效防止模型精度过高造成的模型复杂度过高，但是BIC的使用更倾向于选择参数少的简单模型。
3.Mallows’ Cp 准则:是由Colin Mallows发展的一种AIC的变体，同样是Cp越小，模型越好。
验证前准备：
这个程序包使用前需要将结果变量编码为0或者1，并且需要将结果变量放置在数据最后一列：
X <- train[,1:9]
Xy <- data.frame(cbind(X,train_class)) #将编码为0或1的结果变量放置在最后一列
用以下代码进行交叉验证：
library(bestglm)
bestglm(Xy = Xy, #已经格式化后的数据
IC = 'BIC', #表示使用BIC准则
CVArgs = list(Method = 'HTF',K=10,REP=1), #HTF表示方法为K折交叉验证，K指定了均分的份数量
family = binomial) #表示Logistic回归
得到结果如下：
Morgan-Tatar search since family is non-gaussian.
BIC
BICq equivalent for q in (0.0660008195615313, 0.548953346085246)
Best Model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.7431494 1.16892197 -8.335158 7.739474e-17
V1 0.6817560 0.13279418 5.133930 2.837540e-07
V4 0.3840148 0.12023625 3.193836 1.403961e-03
V6 0.4612582 0.09436587 4.887977 1.018773e-06
V7 0.7595874 0.17602832 4.315143 1.594996e-05
经过所有子集的测试之后，以上4个特征提供了最小BIC评分，因此我们选择以上变量进行模型的重新构建，并且查看新模型在测试集上的预测效果：
#模型重新构建，代码和之前一样
BIC.fit <- glm(formula = class ~ V1 + V3 +V4+ V7, family = binomial, data = train)
test.BIC.pro <- predict(BIC.fit,newdata = test,type = 'response')
test.BIC.pro <- ifelse(test.BIC.pro>= 0.05,1,0)
table(test_class,test.BIC.pro)
得到结果如下：
test.BIC.pro
test_class 0 1
0 114 10
1 0 76
虽然测试集上预测正确率有所下降，但是我们预测变量却从9个变为了4个，尽管纯粹的统计学家会推荐简约的模型作为最终的模型，但是最后选择怎么样的模型仍然是需要我们讨论的问题，并且还有很多机器学习的方法可以使用，后面我们都会慢慢接触到，因此我们可以比较各种算法之后，再结合临床经验，选择一个最可靠并且最容易解释的模型作为我们最终的选择。
想象大家肯定不会陌生这个不管是临床还是科研都经常会用到的图形，在ROC图中，Y轴表示真阳性率 True Positive Rate (TPR),X轴表示假阳性率False Positive Rate(FPR),ROC曲线下的面积就是曲线下面积AUC，表示能正确识别阳性案列的概率，因此我们肯定期望越大越好。
为了方便我们进行比较，我们画3个ROC图：全特征模型，BIC选择的简化模型和单变量模型。接下来我们先构建用于对比效果的单变量模型：
#BIC选择之后的模型重新构建
BIC.fit <- glm(formula = class ~ V1 + V3 +V4+ V7, family = binomial, data = train)
test.BIC.pro <- predict(BIC.fit,newdata = test,type = 'response')
test.BIC.pro <- ifelse(test.BIC.pro>= 0.05,1,0)
#单变量模型构建
singel.fit <- glm(class ~ V1, family = binomial, data = train)
test.singel.pro <- predict(singel.fit,newdata = test,type = 'response')
test.singel.pro <- ifelse(test.singel.pro>= 0.05,1,0)
###画ROC图每个都只需3行代码
#单变量模型画图
pred.singel <- prediction(test.singel.pro,test$class)
pred.singel <- performance(pred.singel,'tpr','fpr')
plot(pred.singel,col=3)
#全特征模型画图
library(ROCR)
pred.full <- prediction(test.probs,test$class)
perf.full <- performance(pred.full,'tpr','fpr')
plot(perf.full,main = 'ROC',col= 1，add= TRUE)
#BIC模型画图
pred.BIC <- prediction(test.BIC.pro,test$class)
pred.BIC <- performance(pred.BIC,'tpr','fpr')
plot(pred.BIC,main = 'ROC',col= 2,add= TRUE)
legend(0.6,0.6,c('FULL','BIC','Singel'),1:3)
联系客服
