断更一天，抱歉，因为完全不会，不知道该记录什么，不过，慢慢的我似乎get到了一些点，试着做一下lasso回归分析。但是分析代码大多是copy的，里面的函数其实我并不是很理解，但是本着我不是专业学统计的，我选择放过自己，了解整个数据分析的原理即可。其实作为一个专业种地的，我是不太懂什么是建模的，更何况如何建模？但经过这几天死缠烂打式的恶补，我来说说我的体会。构建模型其实就是用已知来预测未知，我们可以用已知的具有某一特性的完整数据来拟合出一个可以概括该特性的公式，然后再用这个公式去预测其他的数据是否在以后具有产生这一特性，之前也一直在说，线性模型一般是模型中最稳定的，因此被广泛用于预测，而lasso回归其实就是一个带惩罚的线性回归模型。
之前我们已经通过cox回归得到了一些显著相关的基因了，但是基因数量太多了，里面肯定会包含很多很杂的信息，如果直接用这些基因建模那么可能会出现过拟合的现象，在此基础上，是否能找到更具代表性的基因，进而构建一个更加稳定的模型呢？lasso回归来帮你实现。前一篇说过，lasso回归对于所有的样本都带有惩罚，通过减小损失函数中所有的参数值，将不显著相关的样本去掉，然后剩下的就是显著相关性较大的样本了，用这些样本建模肯定更加可靠。在用基因表达模拟生物学特征这件事上，每个相关基因就是一个样本了，以人为例，不同癌症患者的表达谱存在差异，但是他们都是某一种癌症的患者，因此一定会存在该类癌症相关基因的共性表达，不过我们不知道那些是真正的癌症相关基因，因此，要做个筛选，筛选到可能致癌的基因后，每个基因会带有一个参数(回归系数)，将回归系数与表达量相乘，然后再按个体加起来，就是这个人的风险得分(Risk Score)了，粗狂一点，直接用中位数或者均值为cutoff，高于cutoff的就是高风险了，低于cutoff的就是低风险了。对于这些特征基因，我们当然可以用他们来预测其他人的风险得分，如果这些基因的表达模式与高风险的人一致，那么一定要多多关注自己的健康问题了。这就是癌症预测了。那怎么评判这个模型建的到底行不行，首先很直接的一点就是高风险分类和低风险分类的差异要显著，即用KM曲线看看是不是符合高风险的走的早。还有就是利用ROC曲线，这个之前我们好像说过的，关于两类分类问题，原始类为positive、negative，分类后的类别为p'、n'。排列组合后得到4种结果：
于是得到四个指标，分别为：真阳、伪阳、伪阴、真阴。ROC空间将伪阳性率（FPR）定义为 X 轴，真阳性率（TPR）定义为 Y 轴。这两个值由上面四个值计算得到，公式如下： TPR：在所有实际为阳性的样本中，被正确地判断为阳性之比率。TPR=TP/(TP+FN) FPR：在所有实际为阴性的样本中，被错误地判断为阳性之比率。FPR=FP/(FP+TN) 放在具体领域来理解上述两个指标。如在医学诊断中，判断有病的样本。那么尽量把有病的揪出来是主要任务，也就是第一个指标TPR，要越高越好。而把没病的样本误诊为有病的，也就是第二个指标FPR，要越低越好。不难发现，这两个指标之间是相互制约的。如果某个医生对于有病的症状比较敏感，稍微的小症状都判断为有病，那么他的第一个指标应该会很高，但是第二个指标也就相应地变高。最极端的情况下，他把所有的样本都看做有病，那么第一个指标达到1，第二个指标也为1。 以FPR为横轴，TPR为纵轴，得到如下ROC空间：
左上角的点（TPR=1，FPR=0），为完美分类，也就是这个医生医术高明，诊断全对；点A（TPR>FPR），医生A的判断大体是正确的。中线上的点B（TPR=FPR），也就是医生B全都是蒙的，蒙对一半，蒙错一半；下半平面的点C（TPR<FPR），这个医生说你有病，那么你很可能没有病，医生C的话我们要反着听，为真庸医。 上图中一个阈值，得到一个点。现在我们需要一个独立于阈值的评价指标来衡量这个医生的医术如何，也就是遍历所有的阈值，得到ROC曲线。还是一开始的那幅图，假设如下就是某个医生的诊断统计图，直线代表阈值。我们遍历所有的阈值，能够在ROC平面上得到如下的ROC曲线。
曲线距离左上角越近，证明分类器效果越好。 AUC值为ROC曲线所覆盖的区域面积，显然，AUC越大，分类器分类效果越好。 AUC = 1，是完美分类器，采用这个预测模型时，不管设定什么阈值都能得出完美预测。绝大多数预测的场合，不存在完美分类器。 0.5 < AUC < 1，优于随机猜测。这个分类器（模型）妥善设定阈值的话，能有预测价值。 AUC = 0.5，跟随机猜测一样（例：丢铜板），模型没有预测价值。 AUC < 0.5，比随机猜测还差；但只要总是反预测而行，就优于随机猜测。 假设分类器的输出是样本属于正类的socre（置信度），则AUC的物理意义为，任取一对（正、负）样本，正样本的score大于负样本的score的概率。 AUC值的计算：（1）第一种方法：AUC为ROC曲线下的面积，那我们直接计算面积可得。面积为一个个小的梯形面积之和，计算的精度与阈值的精度有关。 （2）第二种方法：根据AUC的物理意义，我们计算正样本score大于负样本的score的概率。取N*M（N为正样本数，M为负样本数）个二元组，比较score，最后得到AUC。时间复杂度为O(N*M)。 （3）第三种方法：与第二种方法相似，直接计算正样本score大于负样本的score的概率。我们首先把所有样本按照score排序，依次用rank表示他们，如最大score的样本，rank=n(n=N+M)，其次为n-1。那么对于正样本中rank最大的样本（rank_max），有M-1个其他正样本比他score小，那么就有(rank_max-1)-(M-1)个负样本比他score小。其次为(rank_second-1)-(M-2)。最后我们得到正样本大于负样本的概率为： 
所以，搞起来吧！首先下载做lasso回归用的包：glmnet install.packages('glmnet')library(glmnet)phe<-read.csv('phe.csv',quote='',header=T,row.names=1,check.names=F) #获取表型数据data<-read.csv('cox_exprdata.csv',quote='',header=T,row.names=1,check.names=F) #获取表达数set.seed(1) #设置种子，为了可以复现结果fit_cv <- cv.glmnet(t(data), as.matrix(phe$event), alpha=1, family = 'binomial', type.measure='auc') #交叉验证，多次建模，找到最合适的参数lambdaplot(fit_cv) #绘制交叉验证结果
两条虚线分别指示了两个特殊的λ值，lambda.min和lambda.1se，lambda.min 是指在所有的 λ 值中，得到最小目标参量均值的那一个。而 lambda.1se 是指在 lambda.min 一个方差范围内得到最简单模型的那一个 λ值。因为λ值到达一定大小之后，继续增加模型自变量个数即缩小λ值，并不能很显著的提高模型性能， lambda.1se 给出的就是一个具备优良性能但是自变量个数最少的模型。但是大家一般都用的lambda.min。查看该lambda下的样本及参数：get_coe <- function(the_fit,the_lamb){ #提取结果的函数Coefficients <- coef(the_fit, s = the_lamb)Active.Index <- which(Coefficients != 0) #至少coef(回归系数)不能为0啊Active.Coefficients <- Coefficients[Active.Index]re <- data.frame(rownames(Coefficients)[Active.Index],Active.Coefficients)re <- data.table('var_names'=rownames(Coefficients)[Active.Index],'coef'=Active.Coefficients)re$expcoef <- exp(re$coef) #这个可要可不要，后期没啥用return(re[order(expcoef)])}
提取该模型： get_coe(fit_cv,fit_cv$lambda.min) # var_names coef expcoef# 1: (Intercept) -1.12352837 0.3251306# 2: X16 -0.20444072 0.8151031# 3: Genotype.0 -0.15247051 0.8585842# 4: Age 0.02762516 1.0280103# 5: X14 0.15049722 1.1624121# 6: X15 0.77624720 2.1733010
然后，你的模就建好了啊 ，算每个人的风险得分吧：
这是我的美女技术顾问给我提供的公式，其实就是sum(expr*coef)。
model<-as.data.frame(get_coe(fit_cv,fit_cv$lambda.min)) #存为数据框model<-model[-1,] #去掉截距，你要看好截距在哪一行啊model<-model[,-3] #去掉expcoefnames(model)<-c('id','coef') #重新命名model_expr<-merge(model,data,by='id') #整合表达rownames(model_expr)<-model$id #设置行名model_expr<-model_expr[,-1]sample<-model$coef*model_expr #对每一个基因乘以对应的系数sample<-apply(sample,2,sum) #对列求和，2是列，1是行summary(sample) #看一下sample的整体统计情况 RS Min. :-11.08343 #最小值 1st Qu.: -1.47933 #一分位点 1/4 Median : 0.04057 #中位点 Mean : 0.61385 #均值 3rd Qu.: 1.59033 #三分为点 3/4 Max. : 18.77088 #最大值
把mean和median保存下来，好吧，我是直接输入的值，不要在意那么多细节： median <- 0.04057mean <- 0.61385
然后按照均值或者中位数进行分组，具体用什么值，开心就好。
sample$group_median<-ifelse(sample$RS>median,'high','low')sample$group_mean<-ifelse(sample$RS>mean,'high','low')
然后这就是做KM曲线的材料了，绘制KM曲线和上次一样，不多说了。
|
