快来看ls是什么意思(lso是什么意思)

（1）lasso筛选变量，logistic回归建模；（2）lasso回归参数建模；（3）lasso中多分类变量的处理；（4）logistic回归模型列线图绘制与内外部数据验证。

 

接着上次的笔记，此次笔记的任务是利用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(training

set[,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 lassoprint(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 curvenfolds表示进行几折验证，默认是10lsocv ##print(lsocv)。

，glmnet模型交叉验证的结果

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

[1] 0.02179716plot(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/nulldevdevicance函数返回的是(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 variablesin 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=Xt。

rain,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) testdata: trainingset$class, fitted(lgrfit2)

X-squared = 16.432, df = 8, p-value = 0.0366统计学的结果显示校准度一般一个优秀的模型不仅应该对建模的内部数据表现良好，也应该对外的的测试集数据有着不错的预测能力。

文章一开始已经生成了验证集数据：testingset<-biopsy[-n,]，以下命令可以分别获得验证数据的校准度统计学校验、校准度曲线及C-indextestingset<-as.data.frame(testing。

set)pre.test<- predict(lgrfit2,testingset,type = "response")library(ResourceSelection)hoslem.test(testing

set$class,pre.test,g=10)Hosmer and Lemeshow goodness of fit (GOF) testdata: testingset$class, pre.test

X-squared = 5.9867, df = 8, p-value = 0.6487library(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 < casesArea 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$coefficients0]#用lasso的估计系数替换掉lrm的估计系数lsoNom

g<-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$coefficients0]

#用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) testdata: 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 < casesArea under the curve: 0.996logistic回归与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=24 & diabetes$BMI<30,1,diabetes$BMIc

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

$BMI>=35 & diabetes$BMI<40,3,diabetes$BMIc)diabetes$BMIc=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] #去掉变量Insulinlsoc.fit<-glmnet(x,y,family="binomial",alpha=1,dfmax=10)lsoc.fitassess.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$coefficients0]#用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。— END —原创不易，欢迎“在看”关注“一统浆糊”获取更多信息