R語(yǔ)言中使用非凸懲罰函數(shù)回歸(SCAD、MCP)分析前列腺數(shù)據(jù)|附代碼數(shù)據(jù)
原文鏈接：http://tecdat.cn/?p=20828?
最近我們被客戶要求撰寫(xiě)關(guān)于非凸懲罰函數(shù)回歸的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。
本文使用lasso或非凸懲罰擬合線性回歸，GLM和Cox回歸模型的正則化，特別是最小最大凹度懲罰函數(shù)?(MCP)?和光滑切片絕對(duì)偏差懲罰（SCAD），以及其他L2懲罰的選項(xiàng)（ “彈性網(wǎng)絡(luò)”）。還提供了用于執(zhí)行交叉驗(yàn)證以及擬合后可視化，摘要，推斷和預(yù)測(cè)的實(shí)用程序。
我們研究?前列腺數(shù)據(jù)，它具有8個(gè)變量和一個(gè)連續(xù)因變量，即將進(jìn)行根治性前列腺切除術(shù)的男性的PSA水平（按對(duì)數(shù)尺度）：
X <- data$Xy <- data$y
要將懲罰回歸模型擬合到此數(shù)據(jù)，執(zhí)行以下操作：
reg(X, y)
此處的默認(rèn)懲罰是最小最大凹度懲罰函數(shù)?(MCP)?，但也可以使用SCAD和lasso懲罰。這將產(chǎn)生一個(gè)系數(shù)路徑，我們可以繪制
plot(fit)
注意，變量一次輸入一個(gè)模型，并且在λ的任何給定值下，幾個(gè)系數(shù)均為零。要查看系數(shù)是多少，我們可以使用以下?coef?函數(shù)：
coef(fit, lambda=0.05)# (Intercept) ? ? ?lcavol ? ? lweight ? ? ? ? age ? ? ? ?lbph ? ? ? ? svi # ?0.35121089 ?0.53178994 ?0.60389694 -0.01530917 ?0.08874563 ?0.67256096 # ? ? ? ? lcp ? ? gleason ? ? ? pgg45 # ?0.00000000 ?0.00000000 ?0.00168038
該?summary?方法可用于后選擇推斷：
summary(fit # MCP-penalized linear regression with n=97, p=8# At lambda=0.0500:# -------------------------------------------------# ? Nonzero coefficients ? ? ? ? : ? 6# ? Expected nonzero coefficients: ? 2.54# ? Average mfdr (6 features) ? ?: ? 0.424# # ? ? ? ? Estimate ? ? ?z ? ? mfdr Selected# lcavol ? 0.53179 ?8.880 ?< 1e-04 ? ? ? ?*# svi ? ? ?0.67256 ?3.945 0.010189 ? ? ? ?*# lweight ?0.60390 ?3.666 0.027894 ? ? ? ?*# lbph ? ? 0.08875 ?1.928 0.773014 ? ? ? ?*# age ? ? -0.01531 -1.788 0.815269 ? ? ? ?*# pgg45 ? ?0.00168 ?1.160 0.917570 ? ? ? ?*
在這種情況下，?即使調(diào)整了模型中的其他變量之后，lcavol，?svi以及?lweight?顯然與因變量關(guān)聯(lián)，同時(shí)?lbph，?age和?pgg45?可能只是偶然包括。通常，為了評(píng)估模型在λ的各種值下的預(yù)測(cè)準(zhǔn)確性，將執(zhí)行交叉驗(yàn)證：
plot(cvfit)
使交叉驗(yàn)證誤差最小的λ的值由?cvfit$lambda.min給出，在這種情況下為0.017。將coef?在return的輸出?應(yīng)用于?cv.ncvreg?λ的值的系數(shù)：
coef # ?(Intercept) ? ? ? lcavol ? ? ?lweight ? ? ? ? ?age ? ? ? ? lbph ? ? ? ? ?svi # ?0.494154801 ?0.569546027 ?0.614419811 -0.020913467 ?0.097352536 ?0.752397339 # ? ? ? ? ?lcp ? ? ?gleason ? ? ? ?pgg45 # -0.104959403 ?0.000000000 ?0.005324465
可以通過(guò)predict來(lái)獲得預(yù)測(cè)值?，該選項(xiàng)有多種選擇：
predict(cvfit# 預(yù)測(cè)新觀測(cè)結(jié)果的響應(yīng)# ? ? ? ? 1 ? ? ? ? 2 ? ? ? ? 3 ? ? ? ? 4 ? ? ? ? 5 ? ? ? ? 6 # 0.8304040 0.7650906 0.4262072 0.6230117 1.7449492 0.8449595 # 非零系數(shù)的數(shù)量# 0.01695 # ? ? ? 7 # 非零系數(shù)的特性# ?lcavol lweight ? ? age ? ?lbph ? ? svi ? ? lcp ? pgg45 # ? ? ? 1 ? ? ? 2 ? ? ? 3 ? ? ? 4 ? ? ? 5 ? ? ? 6 ? ? ? 8
請(qǐng)注意，原始擬合（至完整數(shù)據(jù)集）的結(jié)果為?cvfit$fit；不必同時(shí)調(diào)用兩者?ncvreg?和?cv.ncvreg?分析數(shù)據(jù)集。
如，?plot(cvfit$fit)?將產(chǎn)生與上述相同的系數(shù)路徑圖?plot(fit)?。
最受歡迎的見(jiàn)解
1.R語(yǔ)言多元Logistic邏輯回歸 應(yīng)用案例
2.面板平滑轉(zhuǎn)移回歸(PSTR)分析案例實(shí)現(xiàn)
3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）
4.R語(yǔ)言泊松Poisson回歸模型分析案例
5.R語(yǔ)言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)
6.r語(yǔ)言中對(duì)LASSO回歸，Ridge嶺回歸和Elastic Net模型實(shí)現(xiàn)
7.在R語(yǔ)言中實(shí)現(xiàn)Logistic邏輯回歸
8.python用線性回歸預(yù)測(cè)股票價(jià)格
9.R語(yǔ)言如何在生存分析與Cox回歸中計(jì)算IDI，NRI指標(biāo)
