二分类结局资料的lasso回归与列线图绘制

接着上次的笔记，此次笔记的任务是利用lasso回归建立预测模型并绘制列线图。

在目前发表的论文中，lasso回归大都只是作为一种变量的筛选方法。首先通过lasso回归获得系数不为0的解释变量，再利用这些筛选到的变量进行多重回归建立预测模型，实际上这是relaxed lasso的一种特殊情况（γ=0）。这种做法用于预测问题不大，但一般不用于因果推断。我们也可以直接利用lasso回归的参数来建模预测，但这样做的人并不多，教程也比较少，但笔者觉得这种方法也是一种不错的选择。

在开始本次笔记的查阅之旅前，建议先看一下《Lasso基础与JMP示例》一文中关于Lasso回归的一些基础介绍。

本次笔记内容：

【1】Lasso回归筛选变量与模型评估

【2】利用lasso筛选的变量进行logistic建模、绘制列线图并进行内外部评估

【3】Lasso回归参数建模、绘制列线图并进行内外部评估

【4】lasso回归中分类变量的处理

示例：MASS程序包中的Biopsy（Biopsy Data on Breast Cancer Patients），数据来自美国University of Wisconsin Hospitals的乳腺癌数据库，数据库含有截至到1992年7月15日的699名乳腺肿瘤患者的活检检查结果，另外保存了患者的9个特征，每个属性都有1到10分的评分。

【0】数据准备，切分数据集

##导入数据集

data(biopsy,package="MASS")

biopsy<-na.omit(biopsy)

biopsy$class<-ifelse(biopsy$class=="benign",0,1)

##按7:3将数据集分割为训练集和测试集合

set.seed(20220612)

ss<-sample(nrow(biopsy),nrow(biopsy)*0.7)

trainingset<-biopsy[ss,]

testingset<-biopsy[-ss,]

很多函数都可以切割数据集，比如sample.split{caTools}、createDataPartition{caret}等。

library(caTools)

set.seed(20220612)

split<-sample.split(biopsy$class, SplitRatio = 0.70)

trainingset2<-subset(biopsy, split==TRUE)

testingset2<-subset(biopsy, split == FALSE) 

library(caret)

library(ggplot2);library(lattice)

set.seed(20220612)

split<-createDataPartition(biopsy$class,p=0.7,list=FALSE)

trainingset3<-biopsy[split,]

testingset3<-biopsy[-split,]

##生成解释变量和结局变量的矩阵格式，glmnet数据格式是矩阵

Xtrain<-as.matrix(trainingset[,2:10])

Ytrain<-as.matrix(trainingset[,11])

Xtest<-as.matrix(testingset[,2:10])

Ytest<-as.matrix(testingset[,11])

【1】Lasso回归筛选变量与模型评估

library(glmnet)

library(Matrix)

##lasso回归

lsofit<-glmnet(Xtrain,Ytrain,family="binomial",alpha=1)

#解释变量须是矩阵形式；family=c("gaussian", "binomial", "poisson", "multinomial", "cox", "mgaussian")，gaussian适用于连续型因变量, mgaussian适用于连续型的多元因变量，poisson适用于计数因变量,binomial适用于二分类因变量, multinomial适用于多分类的因变量, cox适用于生存资料；弹性网络alpha取值在0和1之间，0≤alpha≤1，取值1时拟合lasso回归，取值0时拟合领回归；nlambda为λ值的数量，默认等于100；dfmax和pmax可以分别限定模型中的变量的数量；relax=TRUE表示将进行Relaxed lasso

print(lsofit)

#结果会产生三列结果，分别是Df (非0系数的变量数量), %dev (模型解释的偏差的百分数) 和 Lambda (相应的λ值)。偏差（deviance）即-2倍的Log-likelihood

随着λ值的减少，纳入的变量越来越多，模型能够解释的偏差也越来越大。模型构建第77步就结束了，意味着可解释的偏差不再随着λ值的减少而增加了。模型自由度为4时（截距+4个变量），模型可解释的偏差百分比在41.48%~53.14%，λ值在0.16880~0.22310之间。

##lasso回归系数

coef.apprx<-coef(lsofit, s=0.2) 

#参数s指定具体的λ值。上述命令是提取λ=0.2时的lasso回归系数，此时模型会有4个非0系数变量。精确系数估算可用coef.exact<-coef(lsofit, s=0.2,exact=TRUE, x=Xtrain, y=Ytrain) 

coef.apprx

predict(lsofit,type="coefficients", s=0.2)

#predict也可以用来提取λ某一个取值时的lasso回归系数

#type=c("link","response","coeffificients","nonzero");当想用新的取值做预测时，可用参数newx来指定预测变量的矩阵，可用于type=c("link","response")时；exact参数同可参见coef函数。link给出的是线性预测值，即进行logit变化前的值；response给出的是概率预测值，即进行logit变换之后的值；coefficients给出的是指定λ值的模型系数；nonzero给出指定的定λ值时系数不为零的模型变量

结果表明当λ=0.2时，V1、V2、V3和V6这4个变量系数不为0。

##系数路径图

plot(lsofit,xvar="lambda",label=TRUE)

#参数xvar=c("norm", "lambda", "dev")，用于指定X轴的变量，norm表示横坐标使用系数的L1范数，lambda表示横坐标使用lnλ，而 "dev"表示横坐标采用解释偏差的百分比

图中每一条曲线代表了相应的一个自变量系数的变化轨迹，纵坐标是系数的值，横坐标可以是L1范数、lnλ或者解释偏差的百分比，有些变量在λ值很大时依旧有非零的系数，然后随着λ值变小不断变大，这些都是容易保留在模型中的变量。

虽然函数glmnet可以获得一系列的结果，但在很多情况下，我们可能更喜欢选择其中一个。交叉验证可能是完成这个工作的最简单和最广泛使用的方法。cv.glmnet可以完成交叉验证的主要功能，plot和predict也可以提供更多的支持。

set.seed(20220612) #设置随机种子，保证K折验证的可重复性

lsocv<-cv.glmnet(Xtrain,Ytrain,family="binomial",nfolds=10)

#family同glmnet函数的family；type.measure用来指定交叉验证选取模型的标准，可取值"default", "mse", "deviance", "class", "auc", "mae", "C"。type.measure的默认值是"deviance"，线性模型是squared-error for gaussian models (type.measure="mse" ), logistic和poisson回归是deviance， Cox模型则是偏似然值（partial-likelihood）。deviance即-2倍的对数似然值，mse是实际值与拟合值的mean squred error，mae即mean absolute error，class是模型分类的错误率(missclassification error)，auc即area under the ROC curve。nfolds表示进行几折验证，默认是10

lsocv ##print(lsocv) ，glmnet模型交叉验证的结果

结果中的min和1se对应的λ值指的是使模型偏差最小或者偏差+一个标准误时的λ值。这两个λ值也可以通过以下命令来获得：

lsocv$lambda.min

[1] 0.003721544

lsocv$lambda.1se

[1] 0.02179716

plot(lsocv) #绘制交叉验证曲线

交叉验证曲线中有两条虚线，一条对应的是偏差的模型，另外一条是偏差+一个标准误的模型。

在确定模型的λ值后，可以回代入glmnet函数建模，系数的估计一般采用函数coef.apprx或者predict，前面已经示例过了。系数的估计也可以在cv.glmnet交叉验证后直接通过函数coef或predict来完成。

coef(lsocv,s="lambda.min") #获取使模型偏差最小时λ值的模型系数

coef(lsocv,s="lambda.1se") #获取使模型偏差最小时λ值+一个标准误时的模型系数

cbind2(coef(lsocv,s="lambda.min"),coef(lsocv,s="lambda.1se")) #合并显示

predict(lsocv,newx=Xtest,type="class",s="lambda.min")

#预测λ=s取值时newx的结果，结果类型type=c("link", "response", "class", "coefficients", "nonzero")。link给出的是线性预测值，即进行logit变化前的值，函数默认值；response给出的是概率预测值，即进行logit变换之后的值；clase给出0/1预测值；coefficients给出的是指定λ值的模型系数；nonzero给出指定的定λ值时系数不为0的模型变量

#s=c("lambda.1se", "lambda.min")

通过glmnet拟合了一系列模型（或者选定其中的一个后），我们还需要对模型进行评估，评估常常使用外部的验证集或者测试集，你也可以用内部数据来评估。同建立模型时一样，我们需要有一个预测矩阵，同时需要确定一个评估指标（如Deviance、Class、AUC、MSE、MAE等）来评估一系列λ值时的模型的性能表现。glmnet提供了多个函数使得这些评估变的更容易，如assess.glmnet、roc.glmnet和confusion.glmnet、deviance.glmnet等。

因为lasso回归一般会有很多个λ值（默认100个，当然可能提前收敛而小于100个），这些函数会评估每一个λ值对应的模型指标，结果就会非常多，这里只罗列命令及命令释义：

Ytrain<-factor(Ytrain)

Ytest<-factor(Ytest)

#下面几个函数有些需要响应变量为因子

assess.glmnet(lsofit,newx=Xtrain,newy=Ytrain)

#生成一系列λ下各种模型评价指标Deviance、Class、AUC、MSE、MAE的值

roc.glmnet(lsofit,newx=Xtrain,newy=Ytrain)

#生成一个包含TPR(true positive rate)和FPR(false positive rate)的向量，可用于绘制ROC曲线

confusion.glmnet(lsofit,newx=Xtrain,newy=Ytrain)

#生成分类交叉表，可以得到预测正确的概率.

deviance(lsofit,newx=Xtrain,newy=Ytrain)

#函数devicance返回的结果与函数assess.glmnet返回的devicance结果并不一致，前面print(lsofit)中结果出现的模型解释的偏差的百分数(%dev)，这里也有一个偏差。在函数deviance中找到了关于偏差和零偏差的解释：偏差(deviance)=2*(loglike_sat - loglike)，其中loglike_sat是饱和模型的对数似然值，饱和模型是每个观测都有一个参数的模型；零偏差( Null deviance)=2*(loglike_sat -loglike(Null))，其中NULL model指的是截距模型（Cox回归中是0模型）。dev.ratio=1-deviance/nulldev。devicance函数返回的是(1-dev.ratio)*nulldev，即偏差。assess.glmnet返回的结果跟交叉验证曲线上的偏差值一致，但显然经过了一些变换，具体怎么转换的没找到，但这不影响对结果的判读。

【2】利用lasso筛选的变量进行logistic建模并绘制列线图

这种做法可以看做是Relaxed lasso的一种特殊情况。The idea of the relaxed lasso is to take a glmnet fifitted object, and then for each lambda, refifit the variables in the active set without any penalization. Relaxed lasso将变量选择过程和参数估计过程分成两步来进行，第1步先利用lasso回归获得重要的变量，第2步再对这些变量进行参数估计。第2步中的参数估计并非OLS，而是引进了一个新的调整参数γ，0≤γ≤1，gamma=1时就是传统的lasso回归（glmnet默认relax=FALSE）；gamma=0则是无惩罚拟合，相当于进行多重回归拟合。利用lasso筛选的变量进行logistic建模就是gamma=0时的Relaxed lasso回归。

set.seed(20220612)

rlsocv<-cv.glmnet(Xtrain,Ytrain,family="binomial",relax=TRUE)

predict(rlsocv,newx=Xtrain,type="coefficients",s="lambda.min", gamma=0)

我们可以看到结果跟后面采用logistic回归的结果非常非常接近了，接下来进行系数提取建模就可以了（参见【3】）。在第【1】步中，我们通过lasso回归筛选到了8个解释变量：V1、V3、V4、V5、V6、V7、V8、V9，我们可以直接利用这8个解释变量进行logistic回归并绘制列线图。

library(rms)

library(Hmisc);library(lattice);library(survival);library(Formula);library(ggplot2);library(SparseM)

NM<-datadist(trainingset)

options(datadist='NM')

lgrfit<-lrm(class~V1+V3+V4+V5+V6+V7+V8+V9,data=trainingset,x=T,y=T) #用lasso回归筛选的变量采用lrm{rm}函数进行logistic回归。列线图绘制需要用到该函数的结果

lgrfit

结果显示C-index：0.996，采用同样变量进行的logistic回归的模型参数跟lasso回归的估计的变量参数是不同的，但与γ=0时的Relaxed lasso结果一致。

Nomg<-nomogram(lgrfit,fun=plogis,fun.at=c(.001,.01,.05,0.1,seq(.2,.8,by=.2),0.9,.95,.99,.999),lp=F,funlabel="malignant Rate")

plot(Nomg)

接下来我们可以对模型进行一下评估，评估包括内部验证和外部验证，主要指标是区分度（Discrimination）和校准度（Calibration）验证。区分度常用C-index（AUC），校正度常用校准曲线或者Hosmer and Lemeshow goodness of fit (GOF) 检验。前面绘制列线图时已经得出了内部数据验证的C=0.996，校准曲线绘制命令如下：

Cal<-calibrate(lgrfit,method='boot',B=100)

plot(Cal,xlim=c(0,1.0),ylim=c(0,1.0))

结果显示校准曲线与y=x的直线重合性比较不错，表明模型的校准度还不错。当然你可能还是想从统计学上来看一下，可以通过以下命令：

lgrfit2<-glm(class~V1+V3+V4+V5+V6+V7+V8+V9, family=binomial(link = logit), data=trainingset) #用lasso回归筛选的变量采用glm函数进行logistic回归。校准度的验证需要这个函数的结果

library(ResourceSelection)

hoslem.test(trainingset$class,fitted(lgrfit2),g=10)

Hosmer and Lemeshow goodness of fit (GOF) test

data: trainingset$class, fitted(lgrfit2)

X-squared = 16.432, df = 8, p-value = 0.0366

统计学的结果显示校准度一般。

一个优秀的模型不仅应该对建模的内部数据表现良好，也应该对外的的测试集数据有着不错的预测能力。文章一开始已经生成了验证集数据：testingset<-biopsy[-n,]，以下命令可以分别获得验证数据的校准度统计学校验、校准度曲线及C-index。

testingset<-as.data.frame(testingset)

pre.test<- predict(lgrfit2,testingset,type = "response")

library(ResourceSelection)

hoslem.test(testingset$class,pre.test,g=10)

Hosmer and Lemeshow goodness of fit (GOF) test

data: testingset$class, pre.test

X-squared = 5.9867, df = 8, p-value = 0.6487

library(PredictABEL)

plotCalibration(data=testingset,cOutcom=11,predRisk=pre.test,group=10,rangeaxis=c(0,1))

library(pROC)

auc(testingset$class,pre.test)

Setting levels: control = 0, case = 1

Setting direction: controls < cases

Area under the curve: 0.9955

【3】Lasso回归参数建模并绘制列线图

基本思路是提取lasso的回归系数，然后替代logistic回归模型的参数。

lsocoef<-as.numeric(coef(lsocv,s="lambda.min"))

#提取使模型偏差最小时的λ取值对应的lasso的参数（cv.glmnet函数），并将其转化成数值。也可提取glmnet函数的结果参数，命令如下：

#lsocoef<-as.numeric(coef(lsocoef)

#lsocoef<-predict(lsofit,type="coefficients", s=lsocv$lambda.min)

lgrfit_lsocoef<-lrm(class~V1+V3+V4+V5+V6+V7+V8+V9,data=trainingset,x=T,y=T)

lgrfit_lsocoef$coefficients<-lsocoef[abs(lsocoef)>0]

#用lasso的估计系数替换掉lrm的估计系数

lsoNomg<-nomogram(lgrfit_lsocoef,fun=plogis,fun.at=c(.001,.01,.05,0.1,seq(.2,.8,by=.2),0.9,.95,.99,.999),lp=F,funlabel="malignant Rate")

plot(lsoNomg)

lgrfit2_lsocoef<-glm(class~V1+V3+V4+V5+V6+V7+V8+V9, family=binomial(link = logit), data=trainingset)

lgrfit2_lsocoef$coefficients<-lsocoef[abs(lsocoef)>0]

#用lasso的估计系数替换掉glm的估计系数

---------

##外部验证：校准度

lsopre.test<- predict(lgrfit2_lsocoef,testingset,type = "response")

hoslem.test(testingset$class,lsopre.test,g=10)

Hosmer and Lemeshow goodness of fit (GOF) test

data: testingset$class, lsopre.test

X-squared = 2.9674, df = 8, p-value = 0.9364

##外部验证：区分度

auc(testingset$class,lsopre.test)

Setting levels: control = 0, case = 1

Setting direction: controls < cases

Area under the curve: 0.996

logistic回归与lasso回归外部验证的AUC比较（本例差别不大）

roc(testingset$class,pre.test,plot=TRUE,col="blue")

plot.roc(roc(testingset$class,lsopre.test),add=TRUE,col="red")

还有件事情需要啰嗦一下，函数glm、lrm的估计参数用lasso回归的估计参数替换掉后，对测试集用新的glm、lrm做预测或直接采用lasso回归后的predict的预测结果是一致的。

用替换系数的glm、lrm做预测：

lsopre.test<- predict(lgrfit2_lsocoef,testingset,type = "response")

lasso回归后的predict命令：

predict(lsofit,newx=Xtest,type="response", s=lsocv$lambda.min)

predict(lsocv,newx=Xtest,type="response",s="lambda.min")

#glmnet或者cv.glmnet的预测函数都可以。lsofit和lsocv的模型表达式见【1】Lasso回归筛选变量部分

【4】分类自变量的处理

我们并不会总是有这样的好运气，刚好解释变量全部像示例那样是连续变量。在进行logistic和Cox的lasso回归时，真是世界中的解释变量很多都是分类变量。二分类变量赋值0/1可以直接按连续变量纳入；但多分类变量尤其是无序多分类变量，按连续变量纳入是不合适的。遗憾的是，lasso回归的经典函数glmnet目前只能处理连续变量。

可行的解决的办法是将多分类变量进行哑变量编码（Dummy Encoding），按哑变量进行处理，哑变量遵循同进同出原则进入模型。也可以考虑进行独热编码（One-Hot Encoding）。与哑变量不同的是，独热编码将多分类变量的每个水平都编译成0/1，而哑变量编码中的参照水平编译成0/0，其他水平编译成0/1。我们在《Lasso回归概述与JMP示例》提到JMP的lasso回归可以直接分析分类变量（因子），但需要注意的是JMP对因子的编码采用的是效应编码（Effect Encoding），这种编码也会把K个水平的分类变量编码成K-1个变量，但参照水平会被编码成-1，JMP也默认高水平为参照水平，效应编码的好处是对因子变量而言，截距不再像SPSS那样是因子变量参照水平的值，而是所有水平的均值，因子的系数表示的是该水平与平均水平的差异，不过这不会影响到模型的效应检验。理解这种因子变量的编码差异，有助于对参数做出正确的解读，笔者曾在《再谈有序多分类的logistic回归：JMP操作及OR的再理解》、《重复测量数据分析系列：方差分析【基于JMP】》等文中对这种效应编码和哑变量编码的结果做出过比较。

示例采用R语言VIM包的diabetes数据集（Indian Prime Diabetes Data），数据集介绍参见《Lasso基础与JMP示例》。原数据集存在大量的缺失数据，是VIM包中用于缺失值分析的的一个数据集。考虑到演示的目的，并未对缺失值进行填补。由于原数据集解释变量全部为连续变量，笔者按界值25、30、35、40将BMI转换成因子变量BMIc：正常(0)、超重(1)、一级肥胖(2)、二级肥胖(3)、三级肥胖(4)。

##导入数据，并进行缺失值删除

data(diabetes,package="VIM")

diabetes<-na.omit(diabetes)

##构造分类变量

diabetes$BMIc<-ifelse (diabetes$BMI<24,0,diabetes$BMI)

diabetes$BMIc<-ifelse (diabetes$BMI>=24 & diabetes$BMI<30,1,diabetes$BMIc)

diabetes$BMIc<-ifelse (diabetes$BMI>=30 & diabetes$BMI<35,2,diabetes$BMIc)

diabetes$BMIc<-ifelse (diabetes$BMI>=35 & diabetes$BMI<40,3,diabetes$BMIc)

diabetes$BMIc<-ifelse (diabetes$BMI>=40,4,diabetes$BMIc)

diabetes$BMIc<-factor(diabetes$BMIc)

##构建响应变量和解释变量矩阵

y<-as.matrix(diabetes[,9])

x1<-as.matrix(diabetes[,c("Pregnancies","Glucose","BloodPressure","SkinThickness","Insulin","DiabetesPedigreeFunction","Age","BMIc")]) #解释变量全部为连续变量的多因素矩阵

head(x1)

x2<-model.matrix(~Pregnancies+Glucose+BloodPressure+SkinThickness+Insulin+DiabetesPedigreeFunction+Age+BMIc,data=diabetes)[,-1] #哑变量编码。构建多因素矩阵，model.matrix(~A+B+…)[,-1] 将为因子编码生成变量0/1哑变量，并保持连续变量不变。最后的[,-1]从 model.matrix 的输出中删除常数项

head(x2)

library(cattonum)

diabetes.BMIc.onehot<-catto_onehot(diabetes,10) #独热编码，相当于在x2的基础上增加参照水平的0/1赋值。cattonum包中的函数catto_dummy也可以进行哑变量编码

x3<-as.matrix(diabetes.BMIc.onehot[,c("Pregnancies","Glucose","BloodPressure","SkinThickness","Insulin","DiabetesPedigreeFunction","Age","BMIc0","BMIc1","BMIc2","BMIc3","BMIc4")])

head(x3)

我们在此仅以哑变量编码为例演示一下在lasso回归对多分类变量的处理。这方面教材较少，也请大家批评建议。

library(glmnet);library(Matrix)

set.seed(20220621)

lsocv.dmy<-cv.glmnet(x2,y,family="binomial",nfolds=10)

plot(lsocv.dmy)

lso.dmy<-glmnet(x2,y,family="binomial",alpha=1)

lsodmy.coef<-predict(lso.dmy,type="coefficients", s=lsocv.dmy$lambda.min)

#lsodmy.coef<-predict(lsocv.dmy,newx=x2,type="coefficients",s="lambda.min")

lsodmy.coef

接下来的事情就会很简单，既可以利用这些系数不为0的变量进行logistic建模，只是需要特别注意哑变量应遵循同进同出原则，或者提取lasso回归的系数进行建模。前面已经花了非常多的笔墨来解释这个两个过程，下面只展示一下命令清单和结果。

##lasso筛选的变量进行logistic建模

library(rms)

library(Hmisc);library(lattice);library(survival);library(Formula);library(ggplot2);library(SparseM)

NM<-datadist(diabetes)

options(datadist='NM')

lgsr<-lrm(Outcome~Pregnancies+Glucose+BloodPressure+SkinThickness+DiabetesPedigreeFunction+Age+BMIc, data=diabetes,x=T,y=T)

lgsr

Nom<-nomogram(lgsr,fun=plogis,fun.at=c(.001,.01,.05,0.1,seq(.2,.8,by=.2),0.9,.95,.99,.999),lp=F,funlabel="Diabetes Rate")

plot(Nom)

##lasso回归结果建模

x<-x2[,-5] #去掉变量Insulin

lsoc.fit<-glmnet(x,y,family="binomial",alpha=1,dfmax=10)

lsoc.fit

assess.glmnet(lsoc.fit,newx=x,newy=y)

结果包括不同λ值时的用来评估模型的Deviance、Class、AUC、MSE、MAE等性能指标，此处省略。。。结论是在全部变量都纳入的模型中（哑变量同进同出，相当于10个变量的模型），当λ=0.000272时，模型最佳。

lsoc.coef<-predict(lsoc.fit,type="coefficients", s=0.000272) 

#获得λ=0.000272时的模型系数

#lsoc.fit<-glmnet(x,y,family="binomial",alpha=1,lambda=0.000272)

lsoc.coef<-as.numeric(lsoc.coef) #提取λ=0.000272时lasso的回归系数

library(rms)

library(Hmisc);library(lattice);library(survival);library(Formula);library(ggplot2);library(SparseM)

NM<-datadist(diabetes.BMIc.onehot)

options(datadist='NM')

lgsr.lsocoef<-lrm(Outcome~Pregnancies+Glucose+BloodPressure+SkinThickness+DiabetesPedigreeFunction+Age+BMIc1+BMIc2+BMIc3+BMIc4, data=diabetes.BMIc.onehot,x=T,y=T)

lgsr.lsocoef$coefficients<-lsoc.coef[abs(lsoc.coef)>0]

#用lasso的估计系数替换掉lrm的估计系数

Nomlso<-nomogram(lgsr.lsocoef,fun=plogis,fun.at=c(.001,.01,.05,0.1,seq(.2,.8,by=.2),0.9,.95,.99,.999),lp=F,funlabel="Diabetes Rate")

plot(Nomlso)

因为哑变量的存在，列线图中的BMI对应着4变量，需要注意的是在进行评分预测时，每一个患者只能用其中的一个BMI。