|
这次做的是LASSO回归和多元线性回归的比较，也是大二下一门课的实验，主要用两个数据集，一是RIL.G.Rdata，维度为210*1619，即210个样本，1619个基因型；二是RIL.Phe.Rdata，维度为210*8，即210个样本，8个指标(包含4种表型)。今天做的就是拿基因型数据来预测4种表型中的一个kgw(千粒重)，看看LASSO回归的效果。 library(Matrix)library(glmnet)library(foreach)load('D:/study/sophomore/生信数学基础/RIL.G.Rdata')load('D:/study/sophomore/生信数学基础/RIL.Phe.Rdata')View(G)View(RIL.Phe)phenotype <- as.data.frame(RIL.Phe)View(phenotype)dim(phenotype)#210 8genotype <- as.data.frame(G)View(genotype)dim(genotype)#210 1619#有4个表型，分别是yd(水稻产量)、tp(分蘖数)、gn(单株谷粒数)、kgw(千粒重)，这里选取一个作为因变量genotype$KGW <- phenotype$kgw#数据预处理x <- model.matrix(KGW~.,genotype) #自变量是除了表型的所有变量y <- genotype$KGW #因变量是当前选取的表型set.seed(12) #设置种子数，便于后续复现train <- sample(1:nrow(x),nrow(x)*0.7) #划分训练集和测试集，训练集占0.7x.train <- x[train,]dim(x.train) #查看训练集维数，理论上应该是[nrow(x)*0.7]*ncol(x)#147 1620x.test <- x[-train,]dim(x.test)#63 1620y.train <- y[train] #因变量只有一列length(y.train) #长度与自变量的训练集行数一致#147y.test <- y[-train]length(y.test) #63#构造LASSO模型grid <- 10^seq(10,-2,length=100) #设置一个网格，包含100个从10^(-2)到10^10的数值LASSO.model <- glmnet(x.train,y.train,lambda = grid) #100组lambdaplot(LASSO.model, xvar='lambda', label=TRUE)str(LASSO.model)coef(LASSO.model)[,10] #查看第十组模型的系数上述代码包括数据读入与预处理、数据集划分、LASSO模型构建，下面这个图是LASSO模型的路径图，x轴为log(λ)，y轴是系数，较小的λ值会导致较弱的惩罚和更多的变量被保留在模型中，而较大的λ值会增加惩罚强度，导致更多的变量系数被压缩至零，从而实现变量选择。
接下来查看第十组λ值的拟合结果，可以看到系数几乎都为0，那么可以通过交叉验证得到最优λ值，再进行模型预测。
通过交叉验证可以看到λ最小值为0.0164，也可以通过绘图看一下不同λ取值对MSE的影响。
#接着使用最优λ进行模型预测λ <- cv.LASSO.model$lambda.mincv.LASSO.pred<-predict(cv.LASSO.model,newx=x.test,s=λ)MSE <- mean((cv.LASSO.pred-y.test)^2) #计算均方误差MSE #1.998579SSR <- sum((cv.LASSO.pred-mean(y.test))^2) #计算回归平方和SSR #361.6482SST <- sum((y.test-mean(y.test))^2) #计算总平方和SST #379.6047R_square <- SSR/SST #计算R^2，即决定系数，它是衡量回归模型拟合优度的一个指标R_square #0.952697 # 绘制散点图par(mar=c(5.1, 4.1, 4.1, 4.1)) # 设置图形边界，避免标签超出边界plot(y.test, cv.LASSO.pred, pch=19, col='blue', main='LASSO Regression Predictions', xlab='Actual Values', ylab='Predicted Values', cex.main=1.1, cex.axis=1.1)# 添加R²值到图形usr <- par('usr')#在箱线图上添加p值text(x = (usr[1] (usr[2] - usr[1]) * 0.05), #x坐标在y轴方向偏移一点 y = (usr[3] (usr[4] - usr[3]) * (1-0.05)), #y坐标在x轴方向偏移一点 labels = paste('R² =', round(R_square, 3)), #R²保留2位小数 cex = 1.5, #文本大小 pos = 4) # 添加45度参考线abline(0, 1, col='grey', lty=2)最后就是使用最优λ值进行模型预测，以及计算各个评估指标，绘制拟合图。其中MSE(Mean Squared Error)是均方误差，它能够量化模型预测值与实际观测值之间的差异。MSE越小，表示模型的预测越准确。R²(R-squared)，也称为决定系数，是衡量回归模型拟合优度的一个指标。它表示模型对数据中变异性的解释程度。R²的值介于0和1之间，可以解释为模型对响应变量变异性的百分比。这里的预测结果比较好，R²为0.953，如下图。
|
|
|
