【转】R语言 Logistic回归~变量筛选

发布于：2023/10/17 作者：XIAO油菜花 浏览量：999+ 评论：0

标签： 生物信息 R语言 工具/方法/算法



原文链接：https://zhuanlan.zhihu.com/p/535482448

前排提示：这篇干货只能作为参考，具体实际应用还需根据具体情况进行调参。

 

Logistic回归预测模型思路：

  1. 模型构建
  2. 模型评价
  3. 模型验证

最优模型

  1. 模型能够反映自变量与因变量之间的真实联系
  2. 模型能使用的自变量数目要尽可能的少

模型构建中第一个问题就是变量筛选：最常规的办法是先单后多（先进行单因素分析，单因素有意义的再一起纳入多因素模型中）；但如果变量数目过多，变量间存在共线性或存在较多缺失值而又不愿舍弃含缺失值的样本时，先单后多就存在很多局限性。 此时可以使用具有变量筛选功能的方法： -共线性问题：岭回归(Ridge Regression)，LASSO回归，弹性网络(Elastic Net Regression) -缺失值情况：随机森林模型

常见的筛选变量方法

  1. 正则技术（岭回归、LASSO回归、弹性网络）
  2. 支持向量机
  3. 逐步回归（向前法、向后法、向前向后法）
  4. 最优子集（Best Subset Select）
  5. 树模型（使用的较少）
  6. 随机森林模型
  7. 主成分分析（提取多个自变量的主成分，将主成分得分作为最终的自变量）

常见的模型评价指标

  1. 拟合优度检验（卡方值&P值）
  2. ROC（AUC、sen、spe、accuracy等）
  3. calibration（C-index）
  4. MSE

模型验证是为了防止过拟合情况（所构建的模型对于本次数据有很好的效果，但是对全新的数据效果不理想）

常见的模型验证方法

  1. cross validation（简单交叉，K-fold cross validation，N-fold cross validation）
  2. bootstrap
  3. cross validation + bootstrap（目前最常用）

 

变量筛选~~先单后多

如果某个变量单因素分析时P < 0.05（这个标准可以根据实际案例设置成0.1或者更严格的0.01），就纳入多因素模型。 单因素分析可以用t检验，卡方检验，秩和检验，Logistic回归等。
因变量：Group，1表示疾病，0表示对照 自变量：连续变量和分类变量都有（将分类变量处理成factor--2分类是否处理成因子不影响结果，多分类必须处理成因子）
library(tidyverse) # 加载数据处理包
# 读取数据
data <- read_csv('Clinical/Clinical.RocData.Modified.csv',show_col_types = F)
str(data) # 查看数据类型
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) # 查看各变量的基本统计信息

# 连续型自变量
x1 <- colnames(data)[4:60]
# 分类自变量
x2 <- colnames(data)[2:3]

# t检验--数据符合正态分布
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,
                 # 指定非参数检验变量，exact选项可以指定确切概率检验的变量
                 nonnormal = x1)
write.csv(result2, 'Clinical/单因素分析-秩和检验.csv')

# 单因素Logistic
model <- glm(Group~TP, data = data, family = binomial())
# 查看模型结果
summary(model)$coefficients
# 计算OR值及可信区间
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值及可信区间
exp(cbind('OR' = coef(model), confint(model)))

# 非常规先单后多-筛选协变量
# 指定自变量X是FN，剩下的都是协变量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#'))] #去除一些变量名命名不规范的变量，以免引起错误
lapply(var, covar_method)
# 将阳性的协变量与自变量一起进行多因素回归

 

变量筛选~~LASSO

对于高维数据，普通的变量筛选方法并不见效或者需要消耗大量计算成本和时间成本；且难以避免模型的过拟合及多重共线性问题。此时需要在模型拟合的RSS最小化过程中加入一个正则化项，称之为收缩惩罚；这个惩罚项包含了一个希腊字母λ和对系数β的权重规范化。（RSS+收缩惩罚最小化） 正则化可以对高维数据的系数进行限制，甚至将其缩减到0，避免多重共线性，也可以有效避免过拟合。（岭回归，LASSO，弹性回归） 岭回归中，正则化项是所有变量系数的平方和(L2-norm)，当λ增加时，系数βj缩小，趋向于0，但是永不为0. LASSO回归中，正则化项是变量系数的绝对值和(L1-norm)，这个收缩惩罚项可以使βj收缩到0，因此LASSO具有变量筛选的功能。但是当自变量存在高度共线性或高度两两相关时，LASSO可能会将某个自变量强制删除，这会损失模型预测能力！ 弹性网络中，当α等于0时，弹性网络等价于岭回归；当α等于1时，等价于LASSO；弹性网络技能做到岭回归不能做的变量筛选，又能实现LASSO不能做的变量分组。
library(corrplot) # 相关系数分析用
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数','LYMPH数','MONO数','EO数','BASO数') 
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值的数据
corr <- cor(as.matrix(data))
write.csv(corr, 'Clinical/Correlation.csv')
corrplot.mixed(corr) # 简单进行可视化,查看是存在多重共线性，存在才进行LASSO

# LASSO(正则化技术)不要将分类变量处理成factor，若涉及多分类变量，手动设置哑变量（可用ifelse函数设置）
# 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函数进行岭回归，LASSO回归，弹性网络
library(caret) # 帮助鉴定合适的参数
# 正则化技术需要将数据储存在矩阵里面，而不能是数据框
# 将因变量和自变量处理成矩阵，自变量类型不能是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] # 选最优模型的lambda值
plot(myRidge) # L1范数与系数值之间的关系
plot(myRidge, xvar = 'lambda') # lambda值减小，压缩参数也减小，而系数绝对值增大
plot(myRidge, xvar = 'dev') # 随着偏差百分比增加，系数绝对值增加
# 查看系数
myCoef <- coef(myRidge, s = 4.525281)
write.csv(as.matrix(myCoef),'Clinical/岭回归系数.csv')

# glmnet()语法中alpha=0表示岭回归，1表示LASSO回归
myLasso <- glmnet(X, Y, alpha = 1, family = 'binomial', nlambda = 500) # glmnet默认运行100次
NROW(myLasso$lambda)
myLasso$lambda[500] # 选最优模型的lambda值
myLasso$df[500] # 最优模型留下的变量数
# lambda = 0.004525281 收敛于最优解，有7个变量留下来
# 绘制图形
plot(myLasso, xvar = 'lambda', label = T)
lasso.coef <- coef(myLasso, s = 0.004525281) # 最优解下的回归系数
# 只有筛选出的自变量才有回归系数
write.csv(as.matrix(lasso.coef),'Clinical/LASSO回归系数.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对数值；距离最小均方误差1个标准误时对应的lambda对数值
# 一般第二条虚线对应的是我们的最优解
plot(cv.lasso$glmnet.fit, xvar = 'lambda', label = T)
plot(cv.lasso$glmnet.fit, xvar = 'dev', label = T)
# 如何找到最优lambda：通常距离最小均方误差（MSE）1个标准误时对应的lambda
lasso_lse <- cv.lasso$lambda.1se #提取最优lambda
lasso.coef <- coef(cv.lasso$glmnet.fit, s = lasso_lse, exact = F)
# 没有回归系数的变量即为已剔除变量

# 弹性网络：寻找α和λ的最优组合
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函数设定重抽样的方法LOOCV(留一法)，cv(简单交叉验证)，bootstrp
con <- trainControl(method = 'LOOCV') # 留一法消耗的时间是最多的
# 弹性网络中必须将因变量处理成因子
data$Group <- factor(data$Group)
set.seed(20220629)
enet.train <- train(Group ~ ., data = data, method = 'glmnet',
                    trControl = con, tuneGrid = grid)
enet.train # 选择原则是Accuracy最大，最优参数是α = 0.7；λ = 0.2

# 用最优组合来拟合模型
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）是一类有监督学习方式对数据进行二分类的广义线性分类器，其决策边界是对学习样本求解的最大边距超平面：求解能够正确划分数据集并且几何间隔最大的分离超平面，利用该超平面使任何一类的数据均匀划分。对于线性可分的训练数据而言，线性可分离的超平面有无穷多个，但几何间隔最大的分离超平面是唯一的。 在使用SVM时，并不是所有问题都线性可分，所以需要使用不同的核函数。正因为核函数的使用，使得超平面可以是任意形状。 常用的核函数有：线性核函数，多项式核函数，径向基核函数，sigmoid核函数。 SVM实际应用中，很少会有一个超平面能将不同类别数据完全分开，所以对划分边界近似线性的数据使用软间隔的方法，允许数据跨过超平面，这样会使一些样本分类错误；通过对分类错误样本加以惩罚，可以在最大间隔和确保划分超平面正确分类之间寻找一个平衡。
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数','LYMPH数','MONO数','EO数','BASO数') 
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)
# 线性核函数
# tune.svm函数进行交叉验证选择最优cost成本函数
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)

# 多项式核函数
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)

# 径向基核函数
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核函数
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)

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

# 支持向量机的变量筛选
set.seed(20220629)
# 设置筛选方法
selecMeth <- rfeControl(functions = lrFuncs, 
                        method = 'cv', # 指定是交叉验证
                        number = 10) # 指定nfold数
# rfe函数中有3种变量筛选方法：svmLinner; svmPoly; svmRadial
# size指定自变量个数
# 指定自变量和因变量
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, # 只展现2维的，维度过高不易观察
     svSymbol = 2, # 指定支持向量的形状为三角形
     dataSymbol = 1) # 支持向量之外的数据为圆圈

 

变量筛选~~逐步回归(stepwise)

逐步回归的基本思想：将变量一个一个引入或删除，引入的条件是其偏回归平方和检验呈显著性。每引入一个变量，对已入选回归模型的旧变量逐个进行检验，将不显著的变量删除，以保证回归模型中的每一个自变量都显著。直到不再引入新变量或删除旧变量为止。包括3种方法：向前法，向后法，向前向后法。
AIC赤池信息准则：衡量统计模型拟合优良性（Goodness of fit）的一种标准，通常最小AIC值相对应的模型可作为最优对象。
逐步回归存在的争议：可能得到一个好的模型，但是不能保证该模型是最佳模型，因为不是每一个可能的模型都被评价了。（可以使用自由子集回归的方法）
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数','LYMPH数','MONO数','EO数','BASO数') 
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表示对照
# 自变量：连续变量和分类变量都有（将分类变量处理成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) # 查看各变量的基本统计信息

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

# 逐步回归stepAIC函数在MASS包中
library(MASS)
model.both <- stepAIC(model.old, direction = 'both') # 'forward', 'backward'
# 查看模型结果
summary(model.both)$coefficients
# 计算OR值及可信区间
exp(cbind('OR' = coef(model.both), confint(model.both)))

 

变量筛选~~最优子集

进行回归分析时，通常我们获取的自变量并不是全部有用，一般可以人为根据经验判断筛选对因变量有影响的自变量；如果不是该领域的专家，则需要利用算法获得最贴近真实模型的回归模型，比如最优子集回归。 最优子集回归：对p个预测变量的所有可能组合分别使用回归进行拟合。若有p个解释变量，总共存在2^p个可用于建模的变量子集，根据RSS和R方等评级模型指标的改善情况，从中选择一个最优模型。
# 变量筛选--最优子集
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数','LYMPH数','MONO数','EO数','BASO数') 
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表示对照
# 自变量：连续变量和分类变量都有（将分类变量处理成factor--2分类是否处理成因子不影响结果，
# 多分类必须处理成因子）
# 将分类变量处理成因子
data <- na.omit(data) # 进行NA的行删除

# bestglm包来进行最优子集的选择
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')
# 查看拟合的最优模型
attr(bestAIC$BestModel$terms, 'term.labels')
# 模型拟合
glm(Group~MoCA_B, data = data, family = binomial())

# leaps包来进行最优子集的选择; 函数是regsubsets()
library(leaps)
library(corrplot)
data.cor <- cor(data)
corrplot(data.cor, method = "ellipse") #提示是否存在多重共线性问题
?regsubsets
# 执行最优子集回归，当变量数大于50时需设置really.big = T，默认最优子集最多选择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) #贝叶斯信息准则
# 执行最优子集回归后返回的是自变量组合的子集回归方程，以及每个回归方程对应的评价指标,
# 采用which函数选取最优的回归方程。其中调整R2越大越好，马洛斯Cp越小越好。

# 将返回结果的调整R2作图，可以看到在模型变量个数为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)
# 多重共线性检查
# 将筛选的变量建模并进行共线性检查，方差膨胀系数大于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

 

变量筛选~~随机森林

为提升模型的预测能力，我们可以生成多个树模型，然后将树模型的结果组合起来。随机森林的两个方法：装袋&变量随机。 使用数据集的约2/3的数据建立树模型，剩下的1/3成为袋外数据（验证前期建立的模型的准确率），这个过程重复N次，最后取平均结果。每个树都任其生长，不进行任何基于测量误差的剪枝，这意味着每一个树模型的方差都很大。但是对多个树模型进行平均化，可以降低方差，同时又不增加偏差。 除了对样本进行随机选择，我们对自变量也进行随机选择。对于分类问题，每次抽取的自变量数是自变量总数的平方根；对于回归问题，每次抽取的自变量数是自变量总数的1/3。 随机森林可以对变量的重要性进行评分和排序，并不能实现变量的筛选。
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表示对照
# 自变量：连续变量和分类变量都有（将分类变量处理成factor--2分类是否处理成因子不影响结果，
# 多分类必须处理成因子）
# 将分类变量处理成因子
names(data)[9] <- 'MoCA_B' 
names(data)[22:26] <- c('NEUT数','LYMPH数','MONO数','EO数','BASO数') 
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) # 查看各变量的基本统计信息

# 随机森林的拟合
library(randomForest)
library(randomForest)
set.seed(20220627) # 因为方法有随机选择，所以需要设置随机数，保证代码的复现性
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% # 袋外数据估计的错误率是2.63%
# Confusion matrix:
#   Con AD class.error
# Con  20  0  0.00000000 #预测错误率是0
# AD    1 17  0.05555556 #预测错误率是0.056

# 画误差和树数量的关系图
plot(model1)
# 绿色的线表示AD的误差
# 红色的线表示Con的误差
# 黑色的线表示总样本的误差

# 找出总样本最小误差对应的树数量
which.min(model1$err.rate[,1])

# 重新拟合模型，将ntree = 47，mtry设置每次迭代次数
set.seed(20220627) # 因为方法有随机选择，所以需要设置随机数，保证代码的复现性
model2 <- randomForest(Group~., data = data, ntree = 47, mtry = 7)
model2

# 变量的重要性评分
importance(model2)
# 绘制图形
varImpPlot(model2)

 

-END-

发布评论：

登录后方可评论，点击登录或注册。


评论列表：

暂无评论。

苏公网安备 32050602011302号
苏ICP备2020062135号-1
Copyright© Li Shen. All rights reserved.