"not permitted for Cox family"
搞了好久,才发现生存时间里面有0值,所以必须删除0
基于Lasso回归筛选变量构建Cox模型并绘制Nomogram
研究背景
本章是基于Lasso回归筛选变量后,构建Cox回归临床预测模型,并绘制Nomogram图。Cox模型是一种半参数模型,该模型以生存结局和生存时间为因变量,分析多个因素对生存期的影响,常用RR来量化这种结果,绘制Nomogram列线图实现个体预测。有关Lasso回归可见公众号前文章介绍:如何进行高维变量筛选和特征选择(一)?Lasso回归
案例研究
本文数据收集了83例癌症患者的生存资料,包含患者年龄、性别、癌症分期等。研究目的探讨癌症患者生存情况的影响因素并构建预测模型。
临床研究一般有提供多个危险因素,首先做单因素的筛选,具体筛选方法,见公众号之前的文章。本文采用Lasso回归筛选因素。
具体分析步骤是①筛选变量②基于这些变量构建模型③绘制Nomogram图,预测不同时间生存概率。③计算模型c_index(区分度)该步骤用神包rms一步实现。接下来直接上代码。
R代码及解读
##加载包 明确每个包的作用
library(glmnet) ##Lasso回归
library(rms) ## 画列线图;
library(VIM) ## 包中aggr()函数,判断数据缺失情况
library(survival) ## 生存分析包
#读取数据集
dt <-read.csv("cancer.csv")
str(dt) ##查看每个变量结构
aggr(dt,prop=T,numbers=T) #判断数据缺失情况,红色表示有缺失。
dt <- na.omit(dt) 按行删除缺失值
![](https://i-blog.csdnimg.cn/blog_migrate/03cf31cf1794a09e7458174ef22134cc.png)
由图片可看到所有变量都为蓝色,没有缺失值。如果用na.omit()函数按照行删除。
第一步,也是很重要的一步,数据整理。
用for循环语句将数值型变量转为因子变量
for(i in names(dt)[c(4:9)]){dt[,i]<-as.factor(dt[,i])}##筛选变量前,首先将自变量数据(因子变量)转变成矩阵(matrix)## Lasso要求的数据类型
x.factors <- model.matrix(~ dt$sex+dt$trt+dt$bui+dt$ch+dt$p+dt$stage,dt)[,-1]#将矩阵的因子变量与其它定量边量合并成数据框,定义了自变量。
x <-as.matrix(data.frame(x.factors,dt[,3]))#设置应变量,打包生存时间和生存状态(生存数据)
y <- data.matrix(Surv(dt$time,dt$censor))
第二步:Lasso回归筛选变量
#调用glmnet包中的glmnet函数,注意family那里一定要制定是“cox”,如果是做logistic需要换成"binomial"。
fit <-glmnet(x,y,family ="cox",alpha =1)plot(fit,label=T)plot(fit,xvar="lambda",label=T)
#主要在做交叉验证,lasso
fitcv <- cv.glmnet(x,y,family="cox", alpha=1,nfolds=10)plot(fitcv)coef(fitcv, s="lambda.min")
##
#9 x 1 sparse Matrix of class "dgCMatrix"1
##d.sex1 .
##d.trt1 .
##d.bui1 .
##d.ch2 .
##d.ch3 .
##d.ch4 -0.330676
##d.p1 .
##d.stage4 .
##d...3.
![](https://i-blog.csdnimg.cn/blog_migrate/51fe8cd5a74a07b803f71fd3351fdf60.png)
![](https://i-blog.csdnimg.cn/blog_migrate/b0fb282d9ba9dca78d4775dc4e53cab8.png)
该图在之前文章提到,见如何进行高维变量筛选和特征选择(一)?Lasso回归,由上述代码以及图片完成变量筛选,这里只做演示,假设所有的变量都入选了,我们用这些入选的变量构建Cox回归模型。
第三步:构建Cox模型,并检验等比例风险
#拟合cox回归
coxm <-cph(Surv(time,censor==1)~age+sex+trt+bui+ch+p+stage,x=T,y=T,data=dt,surv=T)
cox.zph(coxm)#等比例风险假定
## chisq df p
##age 1.993 1 0.158
##sex 0.363 1 0.547
##trt 3.735 1 0.053
##bui 2.587 1 0.108
##ch 0.296 1 0.587
##p 0.307 1 0.579
##stage 0.395 1 0.530
##GLOBAL 9.802 7 0.200
注意chp()函数的写法,其中因变量需要用Surv()先打包。后面写法同LR。
等比例风险检验:最后面的GLOBAL是整体看,P值大于0.05,全模型整体都是满足的。对于每一个分类来说P值大于0.05,也是满足的。
第四步:绘制nomogram图,注意该函数里面的参数设置。
surv <-Survival(coxm) # 建立生存函数
surv1 <-function(x)surv(1*3,lp=x) # 定义time.inc,3月OS
surv2 <-function(x)surv(1*6,lp=x) # 定义time.inc,6月OS
surv3 <-function(x)surv(1*12,lp=x) # 定义time.inc,1年OS
dd<-datadist(dt) #设置工作环境变量,将数据整合
options(datadist='dd') #设置工作环境变量,将数据整合
plot(nomogram(coxm,fun=list(surv1,surv2,surv3),
lp= F,
funlabel=c('3-Month Survival','6-Month survival','12-Month survival'),
maxscale=100,fun.at=c('0.9','0.85','0.80','0.70','0.6','0.5','0.4','0.3','0.2','0.1')),
xfrac=.45)
#maxscale 参数指定最高分数,一般设置为100或者10分
#fun.at 设置生存率的刻度
#xfrac 设置数值轴与最左边标签的距离,可以调节下数值观察下图片变化情况
plot(nomogram)
![](https://i-blog.csdnimg.cn/blog_migrate/2ec8547266bf3fe172ceddb0bf1c887a.png)
该图的使用,本质上是将Cox回归模型可视化展示,方便临床快速判断。假设有个病人性别为女,trt为0,P期为1,Nomogram用法是在sex变量上找到其值为1的刻度,然后画垂线投影到最上方的points刻度尺上,找到对应的分值为75分,同理找到trt为0的分值约为50分,P为1的对应分值为100,将这三个因素的points值加起来总分225。下一步在下面的Total Points刻度尺上找到225分,向下方的3个轴做垂线,6-Month-survival对应的值在0.6和0.7之间,约为0.65,说明该患者6个月的生存概率值为65%,其他以此类推。
第五步:利用rms包计算模型区分度。
##模型验证
#Concordance index
f<-coxph(Surv(time,censor==1)~age+sex+trt+bui+ch+p+stage,data=d)
sum.surv<-summary(f)
c_index<-sum.surv$concordance
c_index ##
##C se(C)
##0.55396619 0.07664425
该模型的区分度C-index为0.554,其本质同ROC曲线面积。结果显示,该模型的区分度一般。根据前面变量筛选,考虑纳入更多的影响因素和样本。
函数介绍
glmnet()
glmnet(x, y, family = c("gaussian", "binomial", "poisson", "multinomial",
"cox", "mgaussian"), weights = NULL, offset = NULL, alpha = 1,
nlambda = 100, lambda.min.ratio = ifelse(nobs < nvars, 0.01, 1e-04),
lambda = NULL, standardize = TRUE, intercept = TRUE,
thresh = 1e-07, dfmax = nvars + 1, pmax = min(dfmax * 2 + 20,
nvars), exclude = NULL, penalty.factor = rep(1, nvars),
lower.limits = -Inf, upper.limits = Inf, maxit = 1e+05,
type.gaussian = ifelse(nvars < 500, "covariance", "naive"),
type.logistic = c("Newton", "modified.Newton"),
standardize.response = FALSE, type.multinomial = c("ungrouped",
"grouped"), relax = FALSE, trace.it = 0, ...)
cv.glmnet()
cv.glmnet(x, y,family = c("gaussian", "binomial", "poisson", "multinomial",
"cox", "mgaussian"), weights = NULL, offset = NULL, lambda = NULL,
type.measure = c("default", "mse", "deviance", "class", "auc", "mae",
"C"), nfolds = 10, foldid = NULL, alignment = c("lambda",
"fraction"), grouped = TRUE, keep = FALSE, parallel = FALSE,
gamma = c(0, 0.25, 0.5, 0.75, 1), relax = FALSE, trace.it = 0, ...)
Lasso回归复杂度调整的程度由参数lambda来控制,lambda越大模型复杂度的惩罚力度越大,从而获得一个较少变量的模型。Lasso回归和bridge回归都是Elastic Net广义线性模型的特例。除了参数lambda,还有参数alpha,控制对高相关性数据时建模的形状。Lasso回归,alpha=1(R语言glmnet的默认值),brigde回归,alpha=0,一般的elastic net 0<alpha<1.
根据Hastie(斯坦福统计学家), Tibshirani和Wainwright的Statistical Learning with Sparsity(The Lasso and Generalizations),如下五类模型的变量选择可采用R语言的glmnet包来解决。这五类模型分别是:
1. 二分类logistic回归模型
2. 多分类logistic回归模型
3.Possion模型
4.Cox比例风险模型
5.SVM
下面介绍如何使用glmnet包来实现,以二元logistic回归模型为例:
>library("glmnet") #加载该软件包
>cv.fit<-cv.glmnet(x,y,family="binomial")
#x为输入特征,x应该是矩阵格式的,若非矩阵格式,采用as.matrix()转换成矩阵格式,否则,会报如下错误:Error in lognet(x, is.sparse, ix, jx, y, weights, offset, alpha, nobs, :
(串列)对象不能强制改变成'double'种类。
其他模型familiy值不一样,如cox风险比例模型是cox,possion是possion,多分类logistic是multinomial,广义线性模型guassian,family还有一种选择是mgaussian,不知是否是svm模型?
正确答案是:参数family规定了回归模型的类型:
----family="gaussian"适用于一维连续因变量
----family=mgaussian"适用于多维连续因变量
----family="poisson"适用于非负次数因变量(count)
----family="binomial"适用于二元离散因变量(binary)
----family="multinomial"适用于多元离散因变量(category)
>plot(cv.fit)
![](https://i-blog.csdnimg.cn/blog_migrate/b2d7fecbe35d31697a75e26b5862738c.png)