Logistic回歸預測模型思路：
1.模型構建
2.模型評價
3.模型驗證

最優(yōu)模型
1.模型能夠反映自變量與因變量之間的真實聯(lián)系
2.模型能使用的自變量數(shù)目要盡可能的少

模型構建中第一個問題就是變量篩選：最常規(guī)的辦法是先單后多（先進行單因素分析赌躺，單因素有意義的再一起納入多因素模型中）狼牺；但如果變量數(shù)目過多，變量間存在共線性或存在較多缺失值而又不愿舍棄含缺失值的樣本時礼患，先單后多就存在很多局限性缅叠。
此時可以使用具有變量篩選功能的方法：
-共線性問題：嶺回歸(Ridge Regression)，LASSO回歸弹囚，彈性網(wǎng)絡(Elastic Net Regression)
-缺失值情況：隨機森林模型

常見的篩選變量方法
1.正則技術（嶺回歸鸥鹉、LASSO回歸、彈性網(wǎng)絡）
2.支持向量機
3.逐步回歸（向前法绪撵、向后法祝蝠、向前向后法）
4.最優(yōu)子集（Best Subset Select）
5.樹模型（使用的較少）
6.隨機森林模型
7.主成分分析（提取多個自變量的主成分，將主成分得分作為最終的自變量）

常見的模型評價指標
1.擬合優(yōu)度檢驗（卡方值&P值）
2.ROC（AUC细溅、sen喇聊、spe蹦狂、accuracy等）
3.calibration（C-index）
4.MSE

模型驗證是為了防止過擬合情況（所構建的模型對于本次數(shù)據(jù)有很好的效果凯楔，但是對全新的數(shù)據(jù)效果不理想）

常見的模型驗證方法
1.cross validation（簡單交叉，K-fold cross validation邻遏，N-fold cross validation）
2.bootstrap
3.cross validation + bootstrap（目前最常用）

變量篩選~~先單后多

如果某個變量單因素分析時P < 0.05（這個標準可以根據(jù)實際案例設置成0.1或者更嚴格的0.01）虐骑，就納入多因素模型。
單因素分析可以用t檢驗糊饱，卡方檢驗颠黎，秩和檢驗盏缤，Logistic回歸等。

因變量：Group台舱，1表示疾病，0表示對照
自變量：連續(xù)變量和分類變量都有（將分類變量處理成factor--2分類是否處理成因子不影響結果柜去，多分類必須處理成因子）

library(tidyverse) # 加載數(shù)據(jù)處理包
# 讀取數(shù)據(jù)
data <- read_csv('Clinical/Clinical.RocData.Modified.csv',show_col_types = F)
str(data) # 查看數(shù)據(jù)類型
names(data) # 查看變量名稱
# 將分類變量處理成因子
data$Group <- factor(data$Group, levels = c('0','1'), labels = c('Con', 'AD'))
data$EducationLevel <- factor(data$EducationLevel, levels = c('1','2','3','4','5'),
                              labels = c('小學', '初中', '高中', '大學', '研究生'))
data$Gender <- factor(data$Gender, levels = c('0','1'), labels = c('Female', 'Male'))
summary(data) # 查看各變量的基本統(tǒng)計信息

# 連續(xù)型自變量
x1 <- colnames(data)[4:60]
# 分類自變量
x2 <- colnames(data)[2:3]

# t檢驗--數(shù)據(jù)符合正態(tài)分布
library(tableone)
table1 <- CreateTableOne(vars = c(x1, x2), # 指定對哪些變量進行分析
                         data = data,
                         factorVars = x2, # 指定分類變量
                         strata = 'Group', # 指定分組
                         # 是否對總樣本進行分析
                         addOverall = T)
result1 <- print(table1, 
                 # 是否對分類變量全部展示
                 showAllLevels = T) 
write.csv(result1, 'Clinical/單因素分析-t檢驗.csv')

# 秩和檢驗
library(tableone)
table2 <- CreateTableOne(vars = c(x1, x2), # 指定對哪些變量進行分析
                         data = data,
                         factorVars = x2, # 指定分類變量
                         strata = 'Group', # 指定分組
                         # 是否對總樣本進行分析
                         addOverall = F)
result2 <- print(table2, 
                 # 是否對分類變量全部展示
                 showAllLevels = F,
                 # 指定非參數(shù)檢驗變量嗓奢，exact選項可以指定確切概率檢驗的變量
                 nonnormal = x1)
write.csv(result2, 'Clinical/單因素分析-秩和檢驗.csv')

# 單因素Logistic
model <- glm(Group~TP, data = data, family = binomial())
# 查看模型結果
summary(model)$coefficients
# 計算OR值及可信區(qū)間
exp(cbind('OR' = coef(model), confint(model)))

# 多因素模型
model <- glm(Group~TP+`LYMPH#`+SBP+DBP+NEUT+FN+`APOA1/APOB`+APOA1+ALB+GLB, data = data, family = binomial())
# 查看模型結果
summary(model)$coefficients
# 計算OR值及可信區(qū)間
exp(cbind('OR' = coef(model), confint(model)))

# 非常規(guī)先單后多-篩選協(xié)變量
# 指定自變量X是FN股耽，剩下的都是協(xié)變量Z
covar_method <- function(var){
  model <- glm(Group~FN, data = data, family = binomial())
  coef <- coef(model)[2]
  form <- as.formula(paste0('Group~FN+',var))
  model2 <- glm(form, data = data, family = binomial())
  coef2 <- coef(model2)[2]
  ratio <- abs(coef2-coef)/coef
  if (ratio > 0.1) {
    return(var)
  }
}

var <- c(x1,x2)
var <- var[-which(var %in% c('FN', 'MoCA-B','P-LCR','HDL-c','LDL-c','RDE-SD','APOA1/APOB','A/G','NEUT#',
                             'LYMPH#','MONO#','EO#','BASO#'))] #去除一些變量名命名不規(guī)范的變量物蝙，以免引起錯誤
lapply(var, covar_method)
# 將陽性的協(xié)變量與自變量一起進行多因素回歸

變量篩選~~LASSO

對于高維數(shù)據(jù)敢艰，普通的變量篩選方法并不見效或者需要消耗大量計算成本和時間成本钠导；且難以避免模型的過擬合及多重共線性問題。此時需要在模型擬合的RSS最小化過程中加入一個正則化項票堵，稱之為收縮懲罰换衬；這個懲罰項包含了一個希臘字母λ和對系數(shù)β的權重規(guī)范化证芭。（RSS+收縮懲罰最小化）
正則化可以對高維數(shù)據(jù)的系數(shù)進行限制担映，甚至將其縮減到0蝇完，避免多重共線性，也可以有效避免過擬合氢架。（嶺回歸岖研，LASSO，彈性回歸）
嶺回歸中害淤，正則化項是所有變量系數(shù)的平方和(L2-norm)拓售，當λ增加時础淤，系數(shù)βj縮小，趨向于0莹菱，但是永不為0.
LASSO回歸中吱瘩，正則化項是變量系數(shù)的絕對值和(L1-norm)使碾，這個收縮懲罰項可以使βj收縮到0，因此LASSO具有變量篩選的功能拘鞋。但是當自變量存在高度共線性或高度兩兩相關時矢门，LASSO可能會將某個自變量強制刪除祟剔，這會損失模型預測能力！
彈性網(wǎng)絡中宣旱，當α等于0時浑吟，彈性網(wǎng)絡等價于嶺回歸耗溜；當α等于1時抖拴，等價于LASSO；彈性網(wǎng)絡技能做到嶺回歸不能做的變量篩選轩触，又能實現(xiàn)LASSO不能做的變量分組。

library(corrplot) # 相關系數(shù)分析用
rm(list = ls())
gc()
data <- read_csv('Clinical/Clinical.RocData.Modified.csv',show_col_types = F)
str(data)
names(data) # 查看變量名稱
names(data)[9] <- 'MoCA_B' 
names(data)[22:26] <- c('NEUT數(shù)','LYMPH數(shù)','MONO數(shù)','EO數(shù)','BASO數(shù)') 
names(data)[32:33] <- c('RDE_SD','P_LCR') 
names(data)[37] <- 'A與G的比值' 
names(data)[54:55] <- c('HDL_c','LDL_c') 
names(data)[58] <- 'APOA1與APOB的比值' 
data <- na.omit(data) # 進行NA的行刪除，因為LASSO無法處理含有NA值的數(shù)據(jù)
corr <- cor(as.matrix(data))
write.csv(corr, 'Clinical/Correlation.csv')
corrplot.mixed(corr) # 簡單進行可視化,查看是存在多重共線性榨为，存在才進行LASSO

# LASSO(正則化技術)不要將分類變量處理成factor随闺，若涉及多分類變量，手動設置啞變量（可用ifelse函數(shù)設置）
# 3分類舉例：若data中存在一個變量X有A/B/C三個水平龄句，此時需要將X處理成兩個啞變量
data$X_B <- ifelse(data$X == 'B', 1, 0)
data$X_C <- ifelse(data$X == 'C', 1, 0)

library(glmnet) # 通過glmnet函數(shù)進行嶺回歸分歇，LASSO回歸，彈性網(wǎng)絡
library(caret) # 幫助鑒定合適的參數(shù)
# 正則化技術需要將數(shù)據(jù)儲存在矩陣里面职抡，而不能是數(shù)據(jù)框
# 將因變量和自變量處理成矩陣误甚，自變量類型不能是double，否則報錯
X <- as.matrix(data[2:60])
Y <- as.matrix(data[1])

# glmnet()語法中alpha=0表示嶺回歸擅威，1表示LASSO回歸
myRidge <- glmnet(X, Y, alpha = 0, family = 'binomial', nlambda = 1000)
myRidge$lambda[1000] # 選最優(yōu)模型的lambda值
plot(myRidge) # L1范數(shù)與系數(shù)值之間的關系
plot(myRidge, xvar = 'lambda') # lambda值減小奕翔，壓縮參數(shù)也減小派继，而系數(shù)絕對值增大
plot(myRidge, xvar = 'dev') # 隨著偏差百分比增加捻艳，系數(shù)絕對值增加
# 查看系數(shù)
myCoef <- coef(myRidge, s = 4.525281)
write.csv(as.matrix(myCoef),'Clinical/嶺回歸系數(shù).csv')

# glmnet()語法中alpha=0表示嶺回歸认轨，1表示LASSO回歸
myLasso <- glmnet(X, Y, alpha = 1, family = 'binomial', nlambda = 500) # glmnet默認運行100次
NROW(myLasso$lambda)
myLasso$lambda[500] # 選最優(yōu)模型的lambda值
myLasso$df[500] # 最優(yōu)模型留下的變量數(shù)
# lambda = 0.004525281 收斂于最優(yōu)解，有7個變量留下來
# 繪制圖形
plot(myLasso, xvar = 'lambda', label = T)
lasso.coef <- coef(myLasso, s = 0.004525281) # 最優(yōu)解下的回歸系數(shù)
# 只有篩選出的自變量才有回歸系數(shù)
write.csv(as.matrix(lasso.coef),'Clinical/LASSO回歸系數(shù).csv')

# 通過交叉驗證進行LASSO回歸
lambdas <- seq(0,0.5,length.out = 200)
set.seed(20220629)
# nfolds = 3杉畜，表示3折交叉驗證
cv.lasso <- cv.glmnet(X, Y, alpha = 1, lambda = lambdas, nfolds = 5, family = 'binomial')
plot(cv.lasso) # 縱坐標是MSE
# 兩條虛線：均方誤差最小時對應的lambda對數(shù)值衷恭；距離最小均方誤差1個標準誤時對應的lambda對數(shù)值
# 一般第二條虛線對應的是我們的最優(yōu)解
plot(cv.lasso$glmnet.fit, xvar = 'lambda', label = T)
plot(cv.lasso$glmnet.fit, xvar = 'dev', label = T)
# 如何找到最優(yōu)lambda：通常距離最小均方誤差（MSE）1個標準誤時對應的lambda
lasso_lse <- cv.lasso$lambda.1se #提取最優(yōu)lambda
lasso.coef <- coef(cv.lasso$glmnet.fit, s = lasso_lse, exact = F)
# 沒有回歸系數(shù)的變量即為已剔除變量

# 彈性網(wǎng)絡：尋找α和λ的最優(yōu)組合
grid <- expand.grid(.alpha = seq(0, 1, by = 0.1), .lambda = seq(0, 0.2, by = 0.01))
table(grid)
as.matrix(head(grid))
# trainControl函數(shù)設定重抽樣的方法LOOCV(留一法)随珠，cv(簡單交叉驗證)窗看，bootstrp
con <- trainControl(method = 'LOOCV') # 留一法消耗的時間是最多的
# 彈性網(wǎng)絡中必須將因變量處理成因子
data$Group <- factor(data$Group)
set.seed(20220629)
enet.train <- train(Group ~ ., data = data, method = 'glmnet',
                    trControl = con, tuneGrid = grid)
enet.train # 選擇原則是Accuracy最大，最優(yōu)參數(shù)是α = 0.7软瞎；λ = 0.2

# 用最優(yōu)組合來擬合模型
enet <- glmnet(X, Y, family = 'binomial', alpha = 0.7, lambda = 0.2)
enet.coef <- coef(enet, s = 0.2, exact = T)
enet.coef

變量篩選~~支持向量機

支持向量機（Support Vector Machine涤浇，SVM）是一類有監(jiān)督學習方式對數(shù)據(jù)進行二分類的廣義線性分類器遂唧，其決策邊界是對學習樣本求解的最大邊距超平面：求解能夠正確劃分數(shù)據(jù)集并且?guī)缀伍g隔最大的分離超平面盖彭，利用該超平面使任何一類的數(shù)據(jù)均勻劃分。對于線性可分的訓練數(shù)據(jù)而言铺呵，線性可分離的超平面有無窮多個片挂，但幾何間隔最大的分離超平面是唯一的贞盯。
在使用SVM時，并不是所有問題都線性可分闷愤，所以需要使用不同的核函數(shù)讥脐。正因為核函數(shù)的使用，使得超平面可以是任意形狀俱萍。
常用的核函數(shù)有：線性核函數(shù)枪蘑，多項式核函數(shù)芋齿，徑向基核函數(shù)觅捆，sigmoid核函數(shù)。
SVM實際應用中掂摔，很少會有一個超平面能將不同類別數(shù)據(jù)完全分開乙漓，所以對劃分邊界近似線性的數(shù)據(jù)使用軟間隔的方法释移，允許數(shù)據(jù)跨過超平面玩讳，這樣會使一些樣本分類錯誤；通過對分類錯誤樣本加以懲罰同诫，可以在最大間隔和確保劃分超平面正確分類之間尋找一個平衡误窖。

rm(list = ls())
gc()
data <- read_csv('Clinical/Clinical.RocData.Modified.csv',show_col_types = F) %>% as.data.frame()
data <- na.omit(data) # 進行NA的行刪除
str(data)
names(data) # 查看變量名稱
names(data)[9] <- 'MoCA_B' 
names(data)[22:26] <- c('NEUT數(shù)','LYMPH數(shù)','MONO數(shù)','EO數(shù)','BASO數(shù)') 
names(data)[32:33] <- c('RDE_SD','P_LCR') 
names(data)[37] <- 'A與G的比值' 
names(data)[54:55] <- c('HDL_c','LDL_c') 
names(data)[58] <- 'APOA1與APOB的比值' 

library(e1071)
# 將因變量處理成因子型
data$Group <- factor(data$Group)
# 線性核函數(shù)
# tune.svm函數(shù)進行交叉驗證選擇最優(yōu)cost成本函數(shù)
set.seed(20220629)
linner.tune <- tune.svm(Group ~ ., data = data, kernal = 'linner', cost = c(0.001,0.01,0.1,1,5,10))
summary(linner.tune)
# best parameters: cost = 0.01
# best performance: 0.1583333
best.linner <- linner.tune$best.model
best.linner
# 對擬合的最佳模型進行檢驗
linner.pred <- predict(best.linner, newdata = data)
table(linner.pred, data$Group)

# 多項式核函數(shù)
set.seed(20220629)
poly.tune <- tune.svm(Group ~ ., data = data, kernal = 'polynomial', 
                      degree = c(3,4,5), 
                      coef0 = c(0.1,0.5,1,2,3,4))
summary(poly.tune)
# best parameters: degree = 3, coef0 = 0.1
# best performance: 0.08333333 
best.poly <- poly.tune$best.model
best.poly
# 對擬合的最佳模型進行檢驗
poly.pred <- predict(best.poly, newdata = data)
table(poly.pred, data$Group)

# 徑向基核函數(shù)
set.seed(20220629)
rbf.tune <- tune.svm(Group ~ ., data = data, kernal = 'radial', gamma = c(0.5,1,2,3,4,5,6))
summary(rbf.tune)
# best parameters: gamma = 3
# best performance: 0.425 
best.rbf <- rbf.tune$best.model
best.rbf
# 對擬合的最佳模型進行檢驗
rbf.pred <- predict(best.rbf, newdata = data)
table(rbf.pred, data$Group)

# sigmoid核函數(shù)
set.seed(20220629)
sigmoid.tune <- tune.svm(Group ~ ., data = data, kernal = 'sigmoid', 
                         gamma = c(0.1,0.5,1,2,3,4,5),
                         coef0 = c(0.1,0.5,1,2,3,4,5))
summary(sigmoid.tune)
# best parameters: gamma = 3, coef0 = 0.1
# best performance: 0.425 
best.sigmoid <- sigmoid.tune$best.model
best.sigmoid
# 對擬合的最佳模型進行檢驗
sigmoid.pred <- predict(best.sigmoid, newdata = data)
table(sigmoid.pred, data$Group)

# 四種核函數(shù)最佳模型
library(caret)
confusionMatrix(poly.pred, data$Group, positive = '1')

# 支持向量機的變量篩選
set.seed(20220629)
# 設置篩選方法
selecMeth <- rfeControl(functions = lrFuncs, 
                        method = 'cv', # 指定是交叉驗證
                        number = 10) # 指定nfold數(shù)
# rfe函數(shù)中有3種變量篩選方法：svmLinner; svmPoly; svmRadial
# size指定自變量個數(shù)
# 指定自變量和因變量
X <- data[,2:60]
Y <- data[,1]
svm.feature <- rfe(X,Y,sizes = 30:1, rfeControl = selecMeth, method = 'svmPoly')
svm.feature # TP, MMSE, MoCA_B, DBP
# 利用篩選出的自變量進行分析
svm.4 <- svm(Group ~ TP+MMSE+MoCA_B+DBP, data = data,
             kernel = 'polynomial', degree = 3, coef0 = 0.1)
names(data)
# 對模型進行預測
svm.4.pred <- predict(svm.4, newdata = data[,c(34,8,9,6)])
table(svm.4.pred, data$Group)
# 圖形繪制
plot(svm.4, data = data, 
     TP ~ MoCA_B, # 只展現(xiàn)2維的，維度過高不易觀察
     svSymbol = 2, # 指定支持向量的形狀為三角形
     dataSymbol = 1) # 支持向量之外的數(shù)據(jù)為圓圈

變量篩選~~逐步回歸(stepwise)

逐步回歸的基本思想：將變量一個一個引入或刪除吭服，引入的條件是其偏回歸平方和檢驗呈顯著性艇棕。每引入一個變量串塑，對已入選回歸模型的舊變量逐個進行檢驗桩匪，將不顯著的變量刪除，以保證回歸模型中的每一個自變量都顯著闺骚。直到不再引入新變量或刪除舊變量為止僻爽。包括3種方法：向前法贾惦，向后法，向前向后法碰镜。

AIC赤池信息準則：衡量統(tǒng)計模型擬合優(yōu)良性（Goodness of fit）的一種標準绪颖，通常最小AIC值相對應的模型可作為最優(yōu)對象甜奄。

逐步回歸存在的爭議：可能得到一個好的模型贺嫂，但是不能保證該模型是最佳模型第喳，因為不是每一個可能的模型都被評價了。（可以使用自由子集回歸的方法）

rm(list = ls())
gc()
data <- read_csv('Clinical/Clinical.RocData.Modified.csv',show_col_types = F)
str(data)
names(data) # 查看變量名稱
names(data)[9] <- 'MoCA_B' 
names(data)[22:26] <- c('NEUT數(shù)','LYMPH數(shù)','MONO數(shù)','EO數(shù)','BASO數(shù)') 
names(data)[32:33] <- c('RDE_SD','P_LCR') 
names(data)[37] <- 'A與G的比值' 
names(data)[54:55] <- c('HDL_c','LDL_c') 
names(data)[58] <- 'APOA1與APOB的比值' 
# 因變量：Group悠抹，1表示疾病楔敌，0表示對照
# 自變量：連續(xù)變量和分類變量都有（將分類變量處理成factor--2分類是否處理成因子不影響結果驻谆，
# 多分類必須處理成因子）
# 將分類變量處理成因子
data <- na.omit(data) # 進行NA的行刪除
data$Group <- factor(data$Group, levels = c('0','1'), labels = c('Con', 'AD'))
data$EducationLevel <- factor(data$EducationLevel, levels = c('1','2','3','4','5'),
                              labels = c('小學', '初中', '高中', '大學', '研究生'))
data$Gender <- factor(data$Gender, levels = c('0','1'), labels = c('Female', 'Male'))
summary(data) # 查看各變量的基本統(tǒng)計信息

model.old <- glm(Group~., data = data, family = binomial())
summary(model.old)$coefficients

# 逐步回歸stepAIC函數(shù)在MASS包中
library(MASS)
model.both <- stepAIC(model.old, direction = 'both') # 'forward', 'backward'
# 查看模型結果
summary(model.both)$coefficients
# 計算OR值及可信區(qū)間
exp(cbind('OR' = coef(model.both), confint(model.both)))

變量篩選~~最優(yōu)子集

進行回歸分析時，通常我們獲取的自變量并不是全部有用伙判，一般可以人為根據(jù)經(jīng)驗判斷篩選對因變量有影響的自變量宴抚；如果不是該領域的專家甫煞，則需要利用算法獲得最貼近真實模型的回歸模型抚吠，比如最優(yōu)子集回歸。
最優(yōu)子集回歸：對p個預測變量的所有可能組合分別使用回歸進行擬合蕊玷。若有p個解釋變量垃帅，總共存在2^p個可用于建模的變量子集剪勿，根據(jù)RSS和R方等評級模型指標的改善情況厕吉，從中選擇一個最優(yōu)模型。

rm(list = ls())
gc()
data <- read_csv('Clinical/Clinical.RocData.Modified.csv',show_col_types = F)
str(data)
names(data) # 查看變量名稱
names(data)[9] <- 'MoCA_B' 
names(data)[22:26] <- c('NEUT數(shù)','LYMPH數(shù)','MONO數(shù)','EO數(shù)','BASO數(shù)') 
names(data)[32:33] <- c('RDE_SD','P_LCR') 
names(data)[37] <- 'A與G的比值' 
names(data)[54:55] <- c('HDL_c','LDL_c') 
names(data)[58] <- 'APOA1與APOB的比值' 
# 因變量：Group运悲，1表示疾病班眯，0表示對照
# 自變量：連續(xù)變量和分類變量都有（將分類變量處理成factor--2分類是否處理成因子不影響結果署隘，
# 多分類必須處理成因子）
# 將分類變量處理成因子
data <- na.omit(data) # 進行NA的行刪除

# bestglm包來進行最優(yōu)子集的選擇
library(bestglm) # 因變量要放在最后一列
data2 <- cbind(data[,2:60], data[,1])
# bestglm包中可用的方法有：AIC; BIC; BICg; BICq; LOOCV; CV
bestAIC <- bestglm(data2, family = binomial(), IC = 'AIC')
# Morgan-Tatar search since family is non-gaussian.
# Error in bestglm(data2, family = binomial(), IC = "AIC") : 
#   p = 59. must be <= 15 for GLM.
data2 <- cbind(data[,2:10], data[,1])
bestAIC <- bestglm(data2, family = binomial(), IC = 'AIC')
# 查看擬合的最優(yōu)模型
attr(bestAIC$BestModel$terms, 'term.labels')
# 模型擬合
glm(Group~MoCA_B, data = data, family = binomial())

# leaps包來進行最優(yōu)子集的選擇; 函數(shù)是regsubsets()
library(leaps)
library(corrplot)
data.cor <- cor(data)
corrplot(data.cor, method = "ellipse") #提示是否存在多重共線性問題
?regsubsets
# 執(zhí)行最優(yōu)子集回歸磁餐，當變量數(shù)大于50時需設置really.big = T阿弃，默認最優(yōu)子集最多選擇8個變量
sub.fit <- regsubsets(Group ~ ., data = data[1:20]) 
best.summary <- summary(sub.fit)
# 按照模型評價標準找到評價指標
which.min(best.summary$cp)#馬洛斯Cp值
which.max(best.summary$adjr2) #調整R2
which.min(best.summary$bic) #貝葉斯信息準則
# 執(zhí)行最優(yōu)子集回歸后返回的是自變量組合的子集回歸方程，以及每個回歸方程對應的評價指標,
# 采用which函數(shù)選取最優(yōu)的回歸方程肴楷。其中調整R2越大越好荠呐，馬洛斯Cp越小越好泥张。

# 將返回結果的調整R2作圖媚创，可以看到在模型變量個數(shù)為8的時候彤恶，調整R2最大声离。
plot(best.summary$adjr2, type = "l", xlab = "numbers of Features", 
     ylab = "adjr2", main = "adjr2 by Feature Inclusion")
plot(sub.fit, scale = "adjr2", main = "Best Subset Features")
coef(sub.fit, 8)
# 多重共線性檢查
# 將篩選的變量建模并進行共線性檢查，方差膨脹系數(shù)大于5說明有嚴重的共線性本刽。
# 對這兩個強相關的變量子寓，我們分別做模型笋除，挑選調整R2大的模型垃它。最終我們保留f4模型嗤瞎。
f2 <- lm(Group ~ Age + BMI + MoCA_B + RBC + HGB + MCH + MPV + MCHC, data = data)
library(car)
vif(f2)
#      Age       BMI    MoCA_B       RBC       HGB       MCH       MPV      MCHC 
# 1.194711  1.950985  1.153166 87.278855 88.169592 26.876267  1.125481  2.322166
# 這兩個強相關的變量分別做模型，挑選R2大的模型
f3 <- lm(Group ~ Age + BMI + MoCA_B + RBC + MPV + MCHC, data = data) 
summary(f3) #調整R2:0.8207
f4 <- lm(Group ~ Age + BMI + MoCA_B + HGB + MPV + MCHC, data = data)
summary(f4) #調整R2:0.8212
f5 <- lm(Group ~ Age + BMI + MoCA_B + MCH + MPV + MCHC, data = data)
summary(f5) #調整R2:0.8193

變量篩選~~隨機森林

為提升模型的預測能力虹菲，我們可以生成多個樹模型毕源，然后將樹模型的結果組合起來。隨機森林的兩個方法：裝袋&變量隨機址愿。
使用數(shù)據(jù)集的約2/3的數(shù)據(jù)建立樹模型响谓，剩下的1/3成為袋外數(shù)據(jù)（驗證前期建立的模型的準確率省艳；在沒有驗證數(shù)據(jù)的情況下跋炕，這是隨機森林的一大優(yōu)勢）辐烂，這個過程重復N次，最后取平均結果胳嘲。每個樹都任其生長分瘾，不進行任何基于測量誤差的剪枝胎围，這意味著每一個樹模型的方差都很大。但是對多個樹模型進行平均化德召，可以降低方差白魂，同時又不增加偏差。
除了對樣本進行隨機選擇上岗，我們對自變量也進行隨機選擇福荸。對于分類問題，每次抽取的自變量數(shù)是自變量總數(shù)的平方根肴掷；對于回歸問題敬锐，每次抽取的自變量數(shù)是自變量總數(shù)的1/3呆瞻。
隨機森林可以對變量的重要性進行評分和排序台夺，并不能實現(xiàn)變量的篩選。

隨機森林的優(yōu)勢：
-可以處理大量輸入變量
-可以評估變量的重要性
-建模時使用無偏估計痴脾，模型泛化能力強
-當數(shù)據(jù)缺失較多時颤介，仍可維持一定精確度
-可以處理混合數(shù)據(jù)（數(shù)值型+因子）

rm(list = ls())
gc()
data <- read_csv('Clinical/Clinical.RocData.Modified.csv',show_col_types = F)
str(data)
names(data) # 查看變量名稱
data <- na.omit(data) # 進行NA的行刪除，防止bug
# 因變量：Group，1表示疾病滚朵，0表示對照
# 自變量：連續(xù)變量和分類變量都有（將分類變量處理成factor--2分類是否處理成因子不影響結果冤灾，
# 多分類必須處理成因子）
# 將分類變量處理成因子
names(data)[9] <- 'MoCA_B' 
names(data)[22:26] <- c('NEUT數(shù)','LYMPH數(shù)','MONO數(shù)','EO數(shù)','BASO數(shù)') 
names(data)[32:33] <- c('RDE_SD','P_LCR') 
names(data)[37] <- 'A與G的比值' 
names(data)[54:55] <- c('HDL_c','LDL_c') 
names(data)[58] <- 'APOA1與APOB的比值' 
# 將因變量處理成factor
data$Group <- factor(data$Group, levels = c('0','1'), labels = c('Con', 'AD'))
data$EducationLevel <- factor(data$EducationLevel, levels = c('1','2','3','4','5'),
                              labels = c('小學', '初中', '高中', '大學', '研究生'))
data$Gender <- factor(data$Gender, levels = c('0','1'), labels = c('Female', 'Male'))
summary(data) # 查看各變量的基本統(tǒng)計信息

# 隨機森林的擬合
library(randomForest)
set.seed(20220627) # 因為方法有隨機選擇，所以需要設置隨機數(shù)辕近，保證代碼的復現(xiàn)性
model1 <- randomForest(Group~., data = data)
model1
# Call:
# randomForest(formula = Group ~ ., data = data) 
# Type of random forest: classification # 使用的是分類隨機森林模型
# Number of trees: 500 # 生成了500棵不同的樹
# No. of variables tried at each split: 7 #每次樹的分枝隨機抽取7個變量
# 
# OOB estimate of  error rate: 2.63% # 袋外數(shù)據(jù)估計的錯誤率是2.63%
# Confusion matrix:
#       Con AD class.error
# Con  20  0  0.00000000 #預測錯誤率是0
# AD    1 17  0.05555556 #預測錯誤率是0.056

# 畫誤差和樹數(shù)量的關系圖
plot(model1)
# 綠色的線表示AD的誤差
# 紅色的線表示Con的誤差
# 黑色的線表示總樣本的誤差

# 找出總樣本最小誤差對應的樹數(shù)量
which.min(model1$err.rate[,1])

# 重新擬合模型韵吨，將ntree = 47
set.seed(20220627) # 因為方法有隨機選擇，所以需要設置隨機數(shù)移宅，保證代碼的復現(xiàn)性
model2 <- randomForest(Group~., data = data, ntree = 47)
model2

# 變量的重要性評分
importance(model2)
# 繪制圖形
varImpPlot(model2) 
# 變量的重要性是指每個變量對基尼指數(shù)平均減少量的貢獻归粉。

如果你關注了我，希望你與我一起學習吞杭，一起成長盏浇！?