library(car)#scatterplotMatrix；；；；vip(model,....)计算线性模型、广义线性模型和其他模型的方差-通货膨胀和广义方差-通货膨胀因子。 library(glmnet)#用惩罚极大似然拟合广义线性模型。在正则化参数lambda的数值网格上计算套索或elasticnet惩罚的正则化路径。可以处理各种形状的数据，包括非常大的稀疏数据矩阵。拟合线性，logistic和多项式，泊松和Cox回归模型。 library(psych)#corr.tes library(ncvreg)#SCAD 拟合MCP-或scad -惩罚回归模型的系数路径在一个网格值的正则化参数lambda。拟合线性和逻辑回归模型，选择额外的L2惩罚。 library(msgps)#Adaptive Las #=================Heart数据集实例1=================== options(digits=3) Heart#462*9数据集（ncvreg中） attach(Heart) #数据集中各变量指标含义 #sbp: Systolic blood pressure收缩压 #tobacco: Cumulative tobacco consumption, in kg#累积烟草消费（公斤） #ldl: Low-density lipoprotein cholesterol#低密度脂蛋白胆固醇 #adiposity: Adipose tissue concentration#脂肪组织浓度 #famhist: Family history of heart disease (1=Present, 0=Absent)#heart disease家族史 #typea: Score on test designed to measure type-A behavior#测试a型行为的分数 #obesity: Obesity#肥胖 #alcohol: Current consumption of alcohol#目前的酒精消费量 #age: Age of subject#年龄 famhist=as.factor(Heart$famhist)#转换分类变量famhist为因子#stheart=std(Heart$X);函数std接受一个设计矩阵并返回该矩阵的标准化版本(即，每一列的均值为0，平均平方和为1) HeartF=as.data.frame(Heart$X) #is.data.frame(HeartF) #独立性判断，相关系数及相关性检验 corr.test(HeartF,use="complete")#use=c(pairwise,complete) scatterplotMatrix(~sbp+tobacco+ldl+adiposity+famhist+typea+obesity+alcohol+age,data=HeartF, spread=FALSE, lty.smooth=2, main="ScatterplotmydataHeartF",col="grey") X=model.matrix(sbp~.,HeartF)[,-1]#glmnet函数必须要输入一个预测变量X的设计矩阵和#一个响应向量Y,model.matrix可以生成矩阵(将数据框转换），还可将定性变量转变成哑变量 XS=std(X)#对于具有不同量纲预测变量的数据集，需要对数集进行标准化 Y=HeartF$sbp YS=std(Y) #分割数据集 train=sample(1:nrow(X),nrow(X)/2) test=(-train) #-------------训练集--------------------（x*beta)^3+x*beta---- XTr=XS[train,] # X矩阵 YTr=YS[train] #Y矩阵 #dim(XTr) #length(YTr) #-------------测试集--------------------------- XTe=XS[test,] # X矩阵 YTe=YS[test] #dim(XTe) #=================================GLM============================= #library(gvlma) HeartSF=as.data.frame(std(Heart$X)) heart.fit<- lm(sbp~.,data=HeartSF[train,]) summary(heart.fit) vif(heart.fit) opar=par(no.readonly=TRUE) par(mfrow=c(2,2)) plot(heart.fit) par(opar) predmi=predict(heart.fit,newdata=HeartSF[test,]) errGLM=sqrt(mean((predmi-HeartSF[test,]$sbp)^2)) errGLM #=====================elastic-net(alpha=0.5)====================== HeartEN=glmnet(XTr,YTr,alpha=0.5) plot(HeartEN,label = TRUE) coef(HeartEN) cvHeartEN =cv.glmnet(XTr,YTr,type.measure="mse",alpha=0.5)#交叉验证最优？ plot(cvHeartEN) coef(cvHeartEN,s = "lambda.min")# predmi=predict(cvHeartEN,XTe,s="lambda.min") errEN=sqrt(mean((predmi-YTe)^2)) cvHeartEN$lambda.min errEN #============================Lasso=================== HeartLA=glmnet(XTr,YTr) plot(HeartLA,label = TRUE) cvHeartLA =cv.glmnet(XTr,YTr,type.measure="mse")#交叉验证最优？ plot(cvHeartLA)#默认xvar = "lambda",label = TRUE coef(cvHeartLA,s = "lambda.1se") predmi=predict(cvHeartLA,XTe,s="lambda.min") errLA=sqrt(mean((predmi-YTe)^2)) cvHeartLA$lambda.min errLA #============================SCAD=================== HeartSD=ncvreg(XTr,YTr,family="gaussian",penalty="SCAD",gamma=3.7) plot(HeartSD) cvHeartSD=cv.ncvreg(XTr,YTr,family="gaussian",penalty="SCAD",gamma=3.7) plot(cvHeartSD) plot(cvHeartSD$fit) plot(cvHeartSD$lambda,cvHeartSD$cve,lty=1,pch=20,col="blue") opar=par(no.readonly=TRUE) par(mfrow=c(2,2)) plot(cvHeartSD, type='cve') plot(cvHeartSD, type='scale') # sigma hat plot(cvHeartSD, type='snr') plot(cvHeartSD, type='rsq') par(opar) fit=cvHeartSD$fit cvHeartSD$lambda.min#具有最小交叉验证错误的lambda值coef(cvSDH08) betasd <- fit$beta[,cvHeartSD$min]#cvAUTOSD$min=0.151(对应于lambda。min的索引) betasd summary(cvHeartSD) #summary(cvHeartSD)$min/lambda/r.squared#/对应于lambda.min的索引)/lambda序列/模型解释的方差比例，通过交叉验证估计 predmi=predict(cvHeartSD,XTe,s="lambda.min") errSD=sqrt(mean((predmi-YTe)^2)) errSD #===================Adaptive Lasso---=============================== YTrs=as.vector(YTr) YTes=as.vector(YTe) HeartAL=msgps(XTr,YTrs,penalty="alasso",gamma=1)#(penalty="enet",alpha=0---Lasso);(penalty="enet"/gent",0