弹性网络回归(Elastic Net Regression)是一种用于处理回归问题的机器学习算法，它结合了岭回归(Ridge Regression)和Lasso回归(Least Absolute Shrinkage and Selection Operator Regression)的优点。弹性网络回归在应对高维数据和具有多重共线性的特征时表现出色，它可以有效地减少模型的过拟合，并在存在高度相关的特征时保持稳定性。自适应弹性网络模型则结合了弹性网络回归和自适应正则化技术。
优点：解决多重共线性、稳定、适用于大数据集。
例：用自适应弹性网络分析大气污染物对精子向前活动率的影响
1.数据导入和预处理
#加载包library(dplyr);library(boot);library(glmnet);library(gcdnet);library(SGL);library(ROCR);library(pROC)library(e1071);library(caret);library(DMwR2);library(matrixStats);library(gtools);library(readxl)#读入数据dataall<-read_excel("file path")#指定自变量vocter1<c('ma_age','ma_edu','ma_season','P3PM2.5','P3PM10','P3SO2','P3NO2','P3CO','P3O3')names(dataall)#提取自变量到ind.dataind.data<-dataall[vocter1]#指定因变量vocter2<-c('Lnumber','Lvolume','Ldensity','active','PR')#提取其中一个因变量PR到depen.datadepen.data<-dataall$PR# 检查自变量和因变量数据集中是否有缺失值any(is.na(ind.data))any(is.na(depen.data))#设置受惩罚自变量、不受惩罚协变量个数。（#受惩罚自变量个数6个为'P3PM2.5','P3PM10','P3SO2','P3NO2','P3CO','P3O3'，不受惩罚协变量个数3个为'ma_age','ma_edu','ma_season'）preds<-6covars<-3a <- 1:3 #a是ind.data中未惩罚协变量数量的范围b <- 2:6 #b从aENET结果中提取的行范围一般只排除第一行(截距)####################################################################
✳受惩罚自变量：指在建模过程中受到某种惩罚或约束的自变量。这种惩罚或约束通常是为了控制模型的复杂度、防止过拟合或进行特征选择。它们可以使得部分自变量的系数趋向于零，从而实现变量选择和降维的效果。在本例中为各个污染物。
✳不受惩罚的协变量：不受惩罚的协变量通常指的是在正则化方法中，这些协变量的系数没有受到惩罚，或者受到的惩罚非常小。这些变量在模型拟合过程中不会被强制趋于零，它们的系数可以根据数据的拟合情况来进行估计。在本例中为年龄、教育水平、取精季节。
2.设置五折交叉验证
aENET_bin = function(ind.data, depen.data, preds, covars, a, b, out){#设置随机数种子set.seed(5435672)#为lambda2调优参数创建一个五折交叉验证列表num.obs=nrow(ind.data)#行数foldid=sample(rep(seq(5),length=nrow(ind.data))) #得到一个随机化的1到5之间的整数序列，其长度与 ind.data 的行数相同lambda2.list=seq(0,1,by=0.01) #创建候选lambda2值num.lambda2=length(lambda2.list)#min.cv用于存储每个lambda2值的交叉验证误差min.cv=numeric(num.lambda2)pf_en<-c(rep(1,preds),rep(0,covars))pf_en_2<-c(rep(1,preds),rep(0,covars))length(pf_en)length(pf_en_2)
3.构建函数并运用
#这个函数运行一个连续的结果，如果想运行分类结果需要改变method="logit"#最小二乘方法（Least Squares）进行回归分析，在回归分析中使用损失函数（loss function）来评估预测的准确性for(i in 1:num.lambda2){cv.en=cv.gcdnet(ind.data, depen.data, foldid=foldid, lambda2=lambda2.list[i],method="ls", pred.loss="loss", pf=pf_en, pf2=pf_en_2, standardize=T)min.cv[i]=min(cv.en$cvm)}lambda2.opt=lambda2.list[which.min(min.cv)]lambda.list=cv.en$lambdapar(mar=c(5,5,4,2))cv.en=cv.gcdnet(ind.data, depen.data, foldid=foldid, lambda2=lambda2.opt,method="ls", pred.loss="loss", pf=pf_en, pf2=pf_en_2, standardize=T)plot(cv.en)
#运行函数时会输出以上图片，上图最低处对应具有最小交叉验证误差的最佳惩罚参数值lambda2的log值#将值输入到下面的lambda=exp()的括号中，上图：lambda=exp(-1.8)#初始拟合fit.en = gcdnet(ind.data, depen.data, lambda=exp(), standardize=T,lambda2=lambda2.opt, method="ls", pf=pf_en, pf2=pf_en_2)beta.en = coef(fit.en)v <- log(ncol(ind.data))/log(num.obs)gamma <- ceiling(2*v/(1-v))+1x.sd <- colSds(as.matrix(ind.data))*((num.obs-1)/num.obs)^0.5beta.enet.star <- beta.en[-1,]*x.sdada.wts <- (abs(beta.enet.star)+1/num.obs)^(-gamma)ada.wts[a] <- 0min.cv <- numeric(num.lambda2)#选择最优的lambda2参数for(i in 1:num.lambda2){cv.adenet = cv.gcdnet(ind.data, depen.data, foldid=foldid, lambda2=lambda2.list[i],method="ls", pf=ada.wts, pred.loss="loss", standardize=T)min.cv[i] = min(cv.adenet$cvm)}lambda2.adenet.opt=lambda2.list[which.min(min.cv)]#使用最优lambda2参数进行的最终拟合cv.adenet <- cv.gcdnet(ind.data, depen.data, foldid=foldid, lambda2=lambda2.adenet.opt,method="ls", pf=ada.wts, pred.loss="loss", standardize=T)plot(cv.adenet)
#上一步输入完成后，运行函数时会输出以上图片，上图最低处对应具有最小交叉验证误差的最佳惩罚参数值lambda2的log值#将值输入到下面的lambda=exp()的括号中，上图 lambda=exp(-6.4)fit.all <- gcdnet(ind.data, depen.data, lambda=exp(),lambda2=lambda2.adenet.opt, pf=ada.wts, method="ls", standardize=T )#对数据进行汇总beta.all <- as.matrix(coef(fit.all))#回归系数矩阵化beta.allindex.a=which(beta.all[-1]!=0)#计算除截距项外非0系数的索引p.ne0=length(index.a)x.wi=as.matrix(cbind(1,ind.data))sigma.2=t(depen.data-x.wi%*%beta.all)%*%(depen.data-x.wi%*%beta.all)/(num.obs-p.ne0-1)x.mean=colMeans(ind.data)#x.stdi=matrix(0,num.obs,ncol(ind.data))#class(x.stdi)# 检查矩阵维度print(nrow(x.stdi)) # 输出行数print(ncol(x.stdi)) # 输出列数class(ind.data)class(x.mean)ind.dataa<-as.matrix(ind.data)x.stdi <- matrix(0, nrow = 5503, ncol = 9)#循环将自变量数据标准化for(i in 1:9) {x.stdi[ ,i]<-(ind.dataa[ ,i]-x.mean[i]) / x.sd[i]}class(x.stdi)x.stdi.a=x.stdi[,index.a]#提取非0系数对应的列Sig.a.stdi=t(x.stdi.a)%*%x.stdi.avar.a.stdi=as.numeric(sigma.2)*(1+lambda2.adenet.opt/2)^2*solve(Sig.a.stdi+diag(rep(num.obs*lambda2.adenet.opt,p.ne0))+ (num.obs*lambda2.adenet.opt/2)^2*solve(Sig.a.stdi))x.sd.a.inv=diag(1/x.sd[index.a])var.a.ori=x.sd.a.inv%*%var.a.stdi%*%t(x.sd.a.inv)var.int=t(as.matrix(x.mean[index.a]))%*%var.a.ori%*%as.matrix(x.mean[index.a])name.x.a=c("(Intercept)",colnames(ind.data)[index.a])se.a=data.frame(se=sqrt(c(var.int,diag(var.a.ori))),row.names=name.x.a)pvalue.a=data.frame(pvalue=2-2*apply(abs(beta.all[c(1,index.a+1)])/se.a,1,pnorm),row.names=name.x.a)result.adenet=cbind(numeric(ncol(ind.data)+1),matrix(NA,ncol(ind.data)+1,2))colnames(result.adenet)=c("beta","SE","p-value")rownames(result.adenet)=c("(Intercept)",colnames(ind.data))#构建数据框，包含回归系数、标准误和p值的数据框result.adenet[name.x.a,1]=beta.all[name.x.a,]result.adenet[name.x.a,2]=se.a[name.x.a,]result.adenet[name.x.a,3]=pvalue.a[name.x.a,]result.adenet = as.data.frame(result.adenet)result.adenet$loci = result.adenet$beta - 1.96*result.adenet$SEresult.adenet$hici = result.adenet$beta + 1.96*result.adenet$SE#计算（环境风险评分ERS）权重部分weights = (beta.all[,1])[b]matri = as.matrix(as.data.frame(ind.data[, names(ind.data) %in% names(weights)]))ers = matri%*%as.matrix(weights, ncol = 1)enetlist = list(ers, result.adenet)return(enetlist)}#运用函数依次需要输入变量：自变量数据集、因变量数据集、受惩罚自变量、不受惩罚协变量个数、不受惩罚协变量数量范围和从结果中需要提取得行temp1 = aENET_bin(ind.data, depen.data , 6, 3, 1:3, 2:6, out)temp1[2]
4.得出结果
输出结果：模型帮我们筛选出了NO2、CO、O3这三个自变量，它们的系数β都为负，其中只有O3的p值小于0.05,说明在P3暴露期O3对精子向前活动率有着显著的负向影响。
5.输出结果为excel
library(openxlsx)# 创建一个新的Excel工作簿wb <- createWorkbook()# 在工作簿中创建一个新的工作表addWorksheet(wb, "Sheet1")# 将数据表格写入工作表中writeData(wb, sheet = "Sheet1", x = temp1[2])# 保存工作簿为Excel文件saveWorkbook(wb, file = "C:\\Users\\段甜甜\\Desktop\\P1PR.xlsx")
