我們?cè)诩韧恼隆妒职咽纸棠闶褂肦語(yǔ)言做LASSO 回歸》中介紹了glmnet包進(jìn)行l(wèi)asso回歸，后臺(tái)不少粉絲發(fā)信息向我問(wèn)到分類變量處理的問(wèn)題，我后面查了一下資料之前文章分類變量沒(méi)有處理，非常抱歉?，F(xiàn)在來(lái)重新聊一聊分類變量的處理。

我們導(dǎo)入glmnet包的時(shí)候可以看到，還需要導(dǎo)入一個(gè)Matrix包，說(shuō)明這個(gè)矩陣包很重要

按照glmnet包的原文如下：

就是告訴我們，除了Cox Model外，其他的表達(dá)都支持矩陣形式，在Cox Model的介紹中，函數(shù)樣式為

說(shuō)明我們應(yīng)該把其他變量變?yōu)榫仃嚨男问?。這樣說(shuō)得不是很明白，下面我們來(lái)舉個(gè)例子說(shuō)明，繼續(xù)使用我們的乳腺癌數(shù)據(jù)（公眾號(hào)回復(fù)：乳腺癌，可以獲得數(shù)據(jù)）我們先導(dǎo)入數(shù)據(jù)和R包

library(glmnet)
library(foreign)
library("survival")
bc <- read.spss("E:/r/Breast cancer survival agec.sav",
                use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)

我們先來(lái)看看數(shù)據(jù)：

age表示年齡，pathsize表示病理腫瘤大?。ɡ迕祝?，lnpos表示腋窩淋巴結(jié)陽(yáng)性，histgrad表示病理組織學(xué)等級(jí)，er表示雌激素受體狀態(tài)，pr表示孕激素受體狀態(tài)，status結(jié)局事件是否死亡，pathscat表示病理腫瘤大小類別（分組變量），ln_yesno表示是否有淋巴結(jié)腫大，time是生存時(shí)間，后面的agec是我們自己設(shè)定的，不用管它。

接下來(lái)刪除缺失變量和把分類變量轉(zhuǎn)成因子

bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)
bc$histgrad<-as.factor(bc$histgrad)
bc$pathscat<-as.factor(bc$pathscat)

我們先來(lái)進(jìn)行一個(gè)lasso的cox模型

glmnet包只能接受矩陣形式的數(shù)據(jù)，我們要分別進(jìn)行轉(zhuǎn)換

先把結(jié)局和時(shí)間提取出來(lái)

y<-bc$status
time<-bc$time

把id，結(jié)局變量，時(shí)間變量和一個(gè)亂七八糟的變量刪掉

data1<-bc[,-c(1,8,11,12)]##把id，結(jié)局變量，時(shí)間變量和一個(gè)亂七八糟的變量刪掉

把分類變量變成啞變量矩陣形式

model_mat <-model.matrix(~ +er+pr+ln_yesno+histgrad+pathscat-1,data1)###把分類變量變成啞變量矩陣形式

重新組成數(shù)據(jù)，也就是我們需要的x

x<-as.matrix(data.frame(age=data1$age,
                        pathsize=data1$pathsize,lnpos=data1$lnpos,model_mat))#重新組合成數(shù)據(jù)

弄好x就可以進(jìn)行分析了，交叉驗(yàn)證最好設(shè)一個(gè)種子，

set.seed(123)
cv.fit <- cv.glmnet(x,Surv(time,y),family="cox", maxit = 1000)
plot(cv.fit)

maxit = 1000是讓它迭代100次的意思，如果迭代沒(méi)到1000次，可能會(huì)出現(xiàn)一次報(bào)錯(cuò)，這在官方說(shuō)明里面也有講到，但我用兩種方法算了一遍，結(jié)果都是一樣的，沒(méi)有錯(cuò)

下圖是官方說(shuō)明

有興趣的可以試一下這樣算，結(jié)果也是一樣的，但也要先設(shè)一個(gè)種子

set.seed(123)
cv.fit1<- cv.glmnet(x,Surv(time,y),family="cox", alpha=1,nfolds=10)
plot(cv.fit1)

取最小值，也都是一樣的

cv.fit$lambda.min
cv.fit1$lambda.min

fit <- glmnet(x, Surv(time,y), family =  "cox", maxit = 1000)
plot(fit)

查看和提取系數(shù)

Coefficients <- coef(fit, s = cv.fit$lambda.min)
Active.Index <- which(Coefficients != 0)
Active.Coefficients <- Coefficients[Active.Index]
Active.Index
Active.Coefficients

上圖標(biāo)出了最后還剩下的變量（指的是它的位置）和變量的系數(shù)，自己對(duì)照x看一下就可以了。值得一提的是我看到官方的示例cox模型只取最小的lambda，這樣大家就不用這么糾結(jié)了，還有一個(gè)是它沒(méi)有預(yù)測(cè)功能，不能進(jìn)行預(yù)測(cè)。

下面來(lái)進(jìn)行Binomial Models，也就是我們的二分類變量模型，其實(shí)就是不用時(shí)間變量就行了，其他都差不多，繼續(xù)拿乳腺癌數(shù)據(jù)演示，懶得找數(shù)據(jù)了，上一篇文章就是拿乳腺癌來(lái)模擬二分類數(shù)據(jù)的（當(dāng)時(shí)沒(méi)找到好的數(shù)據(jù)）。

fit1 = glmnet(x, y, family = "binomial")
plot(fit1, xvar = "dev", label = TRUE)

換成lambda

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

其實(shí)到了這里基本和上一篇差不多了

set.seed(999)
cvfit=cv.glmnet(x,y, family = "binomial")
plot(cvfit)

求出最小值

cvfit$lambda.min#求出最小值
cvfit$lambda.1se#求出最小值一個(gè)標(biāo)準(zhǔn)誤的λ值

求出系數(shù)

coef1<-coef(cvfit, s = "lambda.min")
coef2<-coef(cvfit, s = "lambda.1se")
coef1
coef2

有一個(gè)已經(jīng)被懟沒(méi)有了，只能選coef1了。

到此這篇關(guān)于R語(yǔ)言glmnet包lasso回歸中分類變量處理的文章就介紹到這了,更多相關(guān)R語(yǔ)言lasso回歸分類變量?jī)?nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

最新評(píng)論