R语言统计学

1 R语言求平均数、总和、中位数等

read.table("d:\\data.csv",header=TRUE,sep=",")->a
age<-a[,1]
mean(age)
sum(age)
median(age)
var(age) 
sd(age) 
max(age)
quantile(age) quantile(age,0.5)  quantile(age,0.51)

批量运算时，如批量求平均数，标准差、四分位数：

read.table("e:\\data.csv",header=TRUE,sep=",")->a
box<-array(0,dim=c(变量数,2))
for (i in 1:变量数)
{x<-a[1:分组所在行数行数,i+之前的列数]
box[i,1]<-mean (x, na.rm=TRUE)
box[i,2]<- sd (x, na.rm=TRUE)
}
box
View(box)
write.csv(pvalue,file="e:\\mean+sd.csv")

求发病率及95%置信区间

library(Hmisc)
binom.test(发病数, 总人数, con = 0.95)$conf  #binom.test方法
binconf(发病数, 总人数, alpha = 0.05, method = "all") #binconf方法

绘制直方图代码 histogram(直方图)

read.table("d:\\data.csv",header=TRUE,sep=",")->a
age<-a[,1]
hist(age,breaks=12,col="lightgreen",border="pink")

还可以补充：

hist(age,breaks=5,col="lightgreen",main="lightgreen",border="pink")

代码解释： hist为直方图 histogram，breaks为组间距，col为设置填充颜色，border为边界颜色。

常用：

hist(age, breaks = "Sturges", freq = FALSE, probability = !freq, right = TRUE, density = c(10,50), angle=120,angleborder="pink",col = "lightgray", border="pink", main = paste("Histogram of" , xname), xlim = range(breaks), ylim = NULL, xlab = xname, ylab, axes = TRUE, plot = TRUE, labels = FALSE, nclass = NULL, warn.unused = TRUE, ...)

1.breaks

1.breaks为数字向量时（常用）：指定直方图在哪些点截断，一定要包含数据的极值，否则会报错 2.breaks为单个数字时：指定直方图有多少个直方（会有偏差） 3.breaks可以为函数或者字符串值，用来计算直方图的断点

2.freq

freq：设置直方图纵轴时表示频数还是概率密度，FALSE展示概率密度，默认为FALSE

3.probability

和freq意义相反

4. right 直方位于左还是右

4.density，angle

设置直方图中填充斜线的密度和角度。

read.table("d:\\data.csv",header=TRUE,sep=",")->a

age<-a[,2]

hist(age,breaks=5,col="lightgreen",main="lightgreen", density=c(10,50), angle=120,angleborder="pink")

5.labels

1.labels=TRUE，这时候每一个直方图上面的数字代表了，频数（左上图） 2.labels为一个字符向量时，相当于在直方图上加注释

6.main xlab ylab xlim ylim axes plot

分别设置标题，x轴、y轴的轴标题和限制、是否显示坐标轴（这里是x轴和y轴）

密度函数曲线

read.table("d:\\data.csv",header=TRUE,sep=",")->a
B1<-a[1:20,16]
B2<-a[121:703,16]
y1=density(B1) 
plot(y1, col="red",lwd=3)
y2=density(B2)
plot(y2, col="red",lwd=3)

B1为第1-20行，第16列；B2位第121行，第16列；y1、y2为B1/B2的密度；plot为密度曲线，col为颜色，lwd为曲线宽度。

两条曲线绘制于同一图片：

read.table("d:\\data.csv",header=TRUE,sep=",")->a
B1<-a[1:20,16]
B2<-a[121:703,16]
y1=density(B1) 
y2=density(B2)
plot(y1, xlim=c(10,100), ylim=c(0,1),col="red",lwd=3)
par(new=TRUE)
plot(y2, xlim=c(10,100), ylim=c(0,1),col="green",lwd=3)

xlim为横轴的数值范围，ylim为纵轴数值范围（注：不全部纳入范围有可能y1、y2差距过大跑出曲线图）；par(new=TRUE)为将两个曲线置于1个图表。

卡方检验

卡方检验是对两个或两个以上样本率（构成比）进行差别比较的统计方法。例如：病例组、对照组男女比例、生死比例的等 二分类定类变量 检验。

SNP，是否突变

SNP	case	control	Total
0	60	32	92
1	23	15	38
Total	83	47	130

a<-c(60,32,23,15)
dim(a)<-c(2,2)
chisq.test(a,correct=FALSE)
	Pearson's Chi-squared test
data:  a
X-squared = 0.25638, df = 1, p-value = 0.6126
注：样本总数≥40时用chisq.test(a,correct=FALSE)产生卡方检验；样本总数＜40时用
chisq.test(a,correct=FALSE)产生矫正的卡方检验;如果有n＜5时改用Fisher精确检验

做多个（n个）卡方检验时，在exel中如下排列：

Factor1   Factor2   Factor3   Factor4
case1
contro1
case0
contro0

然后读取

read.table("d:\\chisq.csv",header=TRUE,sep=",")->a
box<-array(0,dim=c(n,2))
for (i in 1:n)
{x<-a[,i]
dim(x)<-c(2,2)
box[i,1]<-chisq.test(x,correct=FALSE)$p.value
box[i,2]<-chisq.test(x,correct=FALSE)$statistic
}
box
View(box)
write.csv(box,file="d:\\chisq.testoutcome.csv")

快速制作二联表，group1为分类1，group2为分类2，数据集为data
dataforfisher=data
dataforfisher=rename(dataforfisher,c("group1"=))
dataforfisher=rename(dataforfisher,c("group2"=))
condition1=""
condition2=""

a=nrow(filter(dataforfisher,group1==condition1,group2==condition2))
b=nrow(filter(dataforfisher,group1==condition1,group2!=condition2))
c=nrow(filter(dataforfisher,group1!=condition1,group2==condition2))
d=nrow(filter(dataforfisher,group1!=condition1,group2!=condition2))
dim=c(a,b,c,d)
dim(dim)<-c(2,2)
chisq.test(dim,correct=FALSE)

独立样本的非参数检验

MHO	age	lum	          HIP	         L/Tr	glu0	glu2	HBA1C
YES	34	123.000 	126.000 	0.976 	5.920 	8.230 	6.200 
YES	25	133.000 	137.000 	0.971 	4.460 	5.540 	5.200 
YES	41	119.000 	115.000 	1.035 	4.820 	6.410 	5.400 
YES	21	147.000 	137.000 	1.073 	5.530 	3.320 	6.000 
YES	30	163.000 	145.000 	1.124 	4.810 	6.690 	5.900 
YES	23	104.000 	109.000 	0.954 	3.840 	8.060 	5.100 
YES	30	146.000 	145.000 	1.007 	5.140 	5.730 	5.300 
NO	34	99.000 	        108.000 	0.917 	6.950 	17.710 	6.800 
NO	29	128.000 	132.000 	0.970 	7.430 	12.380 	7.200 
NO	42	113.000 	130.000 	0.869 	5.150 	9.480 	5.500 
NO	31	109.000 	118.000 	0.924 	9.700 	15.200 	8.500 
NO	33	141.000 	135.000 	1.044 	6.220 	6.930 	5.800 
NO	53	100.000 	111.000 	0.901 	6.510 	18.290 	8.200

x个病例组，y个对照组，n个观察因素factors。作疾病和正常的wilcoxon秩和检验。

1.1 单因素的wilcoxon秩和检验，x数值变量，y为二分类变量，数据集data
wilcox.test(x~ y, data=data) #x可为age，y为MHO


1.2 批量检验
read.table("d:\\data.csv",header=TRUE,sep=",")->a
pvalue<-array(0,dim=c(n,1))
for (i in 1:n){test<-wilcox.test(a[,i+1][1:x],a[,i+1][x+1:x+y])
pvalue[i,1]<-test$p.value
}
pvalue
write.csv(pvalue,file="d:\\wilcoxon.csv")

i+1为因子所在列数，前面有几列不需要的，就加几

配对T检验

T检验

T检验（T-test）主要是为了比较数据样本之间是否具有显著性的差异。如病例组、对照组 身高、体重等连续数值的定量型变量，T检验的前提条件是假设样本服从或者近似服从正态分布。

独立样本T检验：不同个体、同一变量的2个样本，如病人和正常人群的身高水平差异；正常组织和肿瘤组织同一基因的表达水平（定量值）

配对样本T检验：同一个体同一变量前后的两个样本，两个如同一病人术前术后血色素水平；小明1月前和1月后考试成绩。

术前（1组）和术后（2组）各有8人，每个人有：Age Bmi Weight Height 4个变量。

假设样本：Age1 Bmi1 Weight1 Height1

配对样本：Age2 Bmi2 Weight2 Height2

共产生4个T检验结果

T(Age)           Age1-Age2
T(Bmi)           Bmi1-Bmi2
T(Weight)      Weight1-Weight2
T(Height)      Height1-Height2

A1和B1各有8个样本

A1	A2	A3	A4	B1	B2	B3	B4
7.5	8.2	8.7	9.4	7.3	7.8	7.2	7.3
7.8	8.7	9.5	9.7	8	8.7	7.4	8.4
7.2	7.9	8.3	9.1	7.1	7.2	7.1	7.2
8.1	9	9.6	9.4	7.7	8.4	7.5	7.9
7.3	8.1	8.1	8.8	7.2	7.5	7.2	7.5
7.8	8.9	9.8	9.9	7.2	8.1	7.1	8.5
7.3	8.5	8.5	8.9	7	7.3	7	7.3
8.5	7.9	8.3	9.9	7.6	7.1	7	7.1

读取上述排布的excel

read.table("d:\\pairT.csv",header=TRUE,sep=",")->a

result<-array(0,dim=c(4,2)) #4为产生配对数（根据factors变化），2为列数（不变），产生一个矩阵：

> result<-array(0,dim=c(4,2))
> result
     [,1] [,2]
[1,]    0    0      
[2,]    0    0    
[3,]    0    0    
[4,]    0    0

然后运算配对T检验

> for (i in 1:4){
+     result[i,1]<-t.test(a[,i]-a[,i+4])$statistic
+     result[i,2]<-t.test(a[,i]-a[,i+4])$p.value
+ }
> result
         [,1]         [,2]
[1,] 2.509980 4.039810e-02
[2,] 5.209043 1.240627e-03
[3,] 7.966582 9.361745e-05
[4,] 9.536924 2.921872e-05

导出结果：

write.csv(result,file="d:\\pairToutcome.csv")

多重线性回归、logistic回归与Cox回归是回归分析中最为常见的三类，三者有差异也有共同之处。

当结局变量为连续型资料时可以考虑多重线性回归

结局变量为分类资料时可以考虑logistic回归

结局变量是生存时间和二分类资料时则考虑Cox回归

COX回归模型

COX回归模型，又称“比例风险回归模型(proportional hazards model，简称Cox模型)”，是由英国统计学家D.R.Cox(1972)年提出的一种半参数回归模型。该模型以生存结局和生存时间为因变量，可同时分析众多因素对生存期的影响，能分析带有截尾生存时间的资料，且不要求估计资料的生存分布类型。

共有m个基因，分布excel分布为：time，status，1，2，n

> data
  time status   ABCA12    ABCC4    ABCF1  ABHD14B     ABI3    ACSL5
1 37.2      1 6.469723 8.249256 9.273143 6.723832 6.506526 9.759622
2 60.0      0 6.804131 8.271230 8.515503 6.760753 6.075105 7.746783
3  7.6      1 6.038480 7.874059 8.964774 7.919579 6.232277 9.874736

1. 单变量，批量COX回归：读取excel

read.table("f:\\cox.csv",header=TRUE,sep=",")->data
pvalue<-array(0,dim=c(n,1))
for (i in 1:n){fit<-coxph(Surv(time,status)~data[,i+2],data=data)
pvalue[i,1]<-summary(fit)$coefficients[,'Pr(>|z|)']
}
pvalue
write.csv(pvalue,file="unicoxoutcome.csv")

1. 分析因素
covariates <- colnames(data)[1:189]

2. 构建生存公式
univ_formulas <- sapply(covariates,
                        function(x) as.formula(paste('Surv(time, event=status)~', x)))

3.对每一个特征做cox回归分析
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = data)})

4.提取HR，95%置信区间和p值  转换成数据框，并转置
univ_results <- lapply(univ_models,
                       function(x){ 
                         x <- summary(x)
                         #获取p值
                         p.value<-signif(x$wald["pvalue"], digits=2)
                         #获取HR
                         HR <-signif(x$coef[2], digits=2);
                         #获取95%置信区间
                         HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
                         HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
                         HR <- paste0(HR, " (", 
                                      HR.confint.lower, "-", HR.confint.upper, ")")
                         res<-c(p.value,HR)
                         names(res)<-c("p.value","HR (95% CI for HR)")
                         return(res)
                       })
res <- t(as.data.frame(univ_results, check.names = FALSE))
res <-as.data.frame(res,stringsAsFactors=F)

5.森林图
pdf(file="univariate_forest.pdf",width=7)
#左边和右边边距稍微留多一点来写变量名称，pvalue和HR
par(mar=c(5,6,4,13))
#先用小方块画出HR
plot(as.numeric(HR[1,]),1:dim(HR)[2],
     pch=15,cex=2,col="blue",bty='n',yaxt='n',ylab=NA,xlab="Hazard Ratio",
     xlim=range(as.numeric(unlist(HR)))
)
#添加中线
abline(v=1,col="grey",lwd=2,lty=2)
​
for(i in 1:ncol(HR)){
  x=as.numeric(HR[2:3,i])
  #循环画出CI
  lines(x,c(i,i),col="blue")
  #添加变量名
  text(0.2,i,rownames(res)[i],xpd=T,adj = c(0,0))
  #添加p值
  text(2.1,i,as.numeric(res[i,1]),xpd=T,adj = c(0,0))
  #添加HR和CI
  text(2.7,i,as.character(res[i,2]),xpd=T,adj = c(0,0))
}
#添加标题
text(2.1,ncol(HR)+0.5,"pvalue",xpd=T,adj = c(0,0))
text(2.7,ncol(HR)+0.5,"HR(CI)",xpd=T,adj = c(0,0))
dev.off()

2. 多元cox

library(survival)
library(rms)
model<-coxph(Surv(time,status)~gene1+gene2+gene3,data,surv=TRUE)

首先查看likehood ration test , wald test, logrank test三种检验方法的p值，p值小于0.05, 这个回归方程是统计学显著的。说明在这么多自变量中包含了对生存时间具有影响的因素。然后查看每个自变量的p值。最后查看自变量的coef等指标，coef就是偏回归系数，exp(coef)就是HR。

Logistic回归

通过logistic回归，可以解决几个问题：（1） 探索疾病的风险因素；（2） 阐明某一个风险因素与疾病的独立关系；（3） 构建由不同风险因素建立的疾病预测模型；（4） 不同预测模型间的比较；（5） 绘制ROC曲线形象化分析模型的效能；（6） 得出模型的风险因素的风险比值比（OR）及95%可信区间、cutoff值、曲线下面积等参数；

常用的统计软件都可以进行logistic回归分析。本文用R语言解决了几个自己关心的问题：

（1） 不同回归模型的比较；

（2） 在一张图中绘制不同回归模型ROC曲线；

第一步，导入和整理数据（本文省略）；

library(stats)
help(infert) # Description of data
infert <- data.frame(infert)
str(infert) # Check type of variables
summary(infert) # Statistical summary

第二步，根据专业知识确定与疾病相关的“可疑”风险因素，也就是根据文献、研究背景认为和疾病相关的因素。关于多因素回归变量的筛选，是一个专题问题，有时间再细细学习。对于流行病等大型数据库，任何一种变量筛选方法都适用，但“统筹考虑统计学上的单因素分析结果与已知临床专业知识决定纳入”变量是最好的考虑。

model.full <- glm(y~.Gender+Age+BMI+xxx, data=yourdata, family='binomial')
summary(model.full) # 做一个logistic回归，纳入所有“可疑风险因素”

第二步，根据初步模型的统计结果，去除不显著的变量后，重新拟合回归模型。

model.reduced <- glm(y~.Gender+Age+BMI+xxx, data=yourdata, family='binomial')
summary(model.reduced)
model.step <-step(model.full) #逐步回归
summary(model.step)

第三步，比较优化前后的模型model.full与model.reduced，可以使用卡方检验（anova()函数）。

anova(model.full, model.reduced, test="Chisq") #两个模型比较

第四步，绘制ROC曲线

install.packages("pROC") #安装pROC包
library(pROC)
yourdata$Pre.full <- predict(model.full,type='response') # 建立预测变量
yourdata$Pre.reduced <- predict(model.reduced,type='response')
# 绘制一条ROC曲线Roc.full
Roc.full <-plot.roc(yourdata$y~, yourdata$Pre.full, percent=TRUE, col="1")
# 添加一条ROC曲线Roc.reduced
Roc.reduced <-lines.roc(yourdata$y~, yourdata$Pre.reduced, percent=TRUE, col="2")
# 可以按此方式继续添加不同模型的ROC曲线
ROC_3 <- lines.roc(yourdata$y~, yourdata$Pre_3, percent=TRUE, col="3")

第五步：计算相关参数

auc(model.full) #计算ROC曲线下面积
auc(model.reduced)
coef(model.full) #计算回归系数
coef(model.reduced)
confint(model.full) # 计算95%置信区间
confint(model.reduced)
exp(coef(model.full))#指数化为优势比
exp(coef(model.reduced))
exp(confint(model.full))# 计算优势比系数的95%置信区间
exp(confint(model.reduced))
predict(model1, type="risk") # predicted values 预测值
residuals(model1, type="deviance") # residuals 残差
#两个ROC曲线比较计算P值
testobj<- roc.test(model.full, model.reduced)

第六步，美化修饰ROC图中元素

#添加文本，可以将两条ROC曲线比较的P值加入图中
text(50, 50, labels=paste("p-value=", format.pval(testobj_1$p.value)), adj=c(0, .5))
#添加图例
legend("bottomright", legend=c("Model-full ", "Model-reduced"), col=c("1", "2"), lwd=2)

第七步，计算Youden系数，寻找cutoff值，R语言和SAS中似乎没有直接给出ROC曲线cutoff值的方法。在R语言中，可以从ROC曲线的数据文件中提取sensitivities和specificities数据进行计算，查找Youden系数最大的预测值（下一节预测模型中将详细介绍）。

Roc.reduced$Youden=Roc.reduced$sensitivities+Roc.reduced$specificities-1

批量

批量做Logistic回归：第1个基因在第k列，共有m个基因/factors（要logitic分析的因素）

read.table("d:\\logit.csv",header=TRUE,sep=",")->a
result<-array(0,dim=c(m,4)) ##建立m行，4列的矩阵；设置4列：分别放入OR值，OR值置信区间的上下限，P值
for (i in 1:m){logr<-glm(group~sex+age+duration+a[,i+k-1],data=a, family="binomial")
result[i,1]<-exp(summary(logr)$coefficients[k,1])
for (i in 1:39){logr<-glm(group~sex+age+bmi+duration+a[,i+5],data=a, family="binomial")
result[i,1]<-exp(summary(logr)$coefficients[5,1])
result[i,2]<-exp(summary(logr)$coefficients[5,1]-1.96*summary(logr)$coefficients[5,2])
result[i,3]<-exp(summary(logr)$coefficients[5,1]+1.96*summary(logr)$coefficients[5,2])
result[i,4]<-summary(logr)$coefficients[5,4]
}
result              #result[i,1]后面固定地为[5,1] [5,2] [5,4]
write.csv(result,file="d:\\logisticoutcome.csv")

Lasso回归

在该数据集中，主要包括患者信息，生存时间（futime），生存状态（fustat），以及相关基因的表达情况；当然，在此，我们只选取了数据集中一部分基因的表达。通过筛选，对高维度的数据进行降维，以获得预后相关的基因特征，进而构建预后相关的预测模型。

set.seed(2)   #设定随机种子 
x=as.matrix(rt[,c(3:ncol(rt))]) 
y=data.matrix(Surv(rt$futime,rt$fustat)) 

实操

library(glmnet)
library(foreign)
library(caret)
library(survival)
## 读取数据 目前，glmnet包只能接受矩阵形式的数据，数据框的数据会报错，所以我们先要把数据转换成矩阵形式。
data= 
data=as.matrix(data)
data=na.omit(data)

1.1 初级：构建x表达矩阵和y生存矩阵
x=sur[,4:29]
y=data.matrix(Surv(time=sur$OS.time,event=sur$OS))
fit=glmnet(x,y,family="cox",maxit=10000)
plot(fit,xvar="lambda",label=T)

该图中展示了，在模型建立过程中，在不同的λ状态下入选特征的数量和系数；而且，随着λ的增大，入选的特征参数随之减小，而系数绝对值却随之增大。同时，有相当都的系数一直接近于0，直到最后一个特征加入模

####进行交叉验证,如报错#问题是RealRatingMatrix类扩展了矩阵，而矩阵没有实现包含字符的矩阵。先将你的矩阵转换成数字，
####然后转换，x<-sapply(data.frame(x),as.numeric)
lasso_fit=cv.glmnet(x,y,family="cox",type.measure="deviance")
plot(lasso_fit,label=T)
#其中两条虚线分别指示了两个特殊的λ值 
abline(v = log(c(lasso_fit$lambda.min,lasso_fit$lambda.1se)),lty="dashed")
v = log(c(lasso_fit$lambda.min,lasso_fit$lambda.1se))
v
[1] -2.778627 -1.476155

在该结果中，随着特征数量的不断选择与模拟，最终得到两个模型，分别为最佳模型（左侧虚线）和最简模型（右侧虚线）；其中，上方的数值代表该λ值时纳入模型的特征数量。

###筛选变量##输出预测模型的相关系数与riskScore
coefficient=coef(lasso_fit,s=lasso_fit$lambda.min)
Active.Index=which(as.numeric(coefficient)!=0)
actCoef=as.numeric(coefficient)[Active.Index]
lassoGene=rownames(coefficient)[Active.Index]
lassoGene  #查看模型的基因
geneCoef = cbind(Gene=lassoGene,Coef=actCoef) 
geneCoef   #查看模型的相关系数
      Gene      Coef                 
 [1,] "DCN"     "-0.0574874509423495"
 [2,] "SCARB1"  "0.00479746322116437"
 [3,] "HMOX1"   "-0.172042730652799" 
 [4,] "STC2"    "0.128564983984361"  
 [5,] "FAM162A" "0.227631199155787"  
 [6,] "MXI1"    "0.177346734946575" 


###4.计算riskScore
随后，我们进一步计算模型的最终得分，即通过将每个基因的表达情况与其对应的系数相乘，最好相加得到最终的风险评分（riskScore）。
并且，以所有患者的中位riskScore值为临界值，将患者分成高风险和低风险两组。
FinalGeneExp = data[,lassoGene] 
myFun = function(x){crossprod(as.numeric(x),actCoef)} 
riskScore = apply(FinalGeneExp,1,myFun) 
outCol = c("OS.time", "OS", lassoGene) 
risk = as.vector(ifelse(riskScore > median(riskScore), "high", "low")) 
dat = cbind(x[,outCol], riskScore=as.vector(riskScore), risk)

###5.绘制散点分布图
接下来，通过使用ggpubr包，绘制散点图来对不同生存状态下患者的风险评分进行可视化展示。


library(ggpubr)   #使用ggpubr包绘制散点图
dat=as.data.frame(dat)
p <- ggboxplot(dat, x = "OS", y = "riskScore", 
              color = "OS", palette = "jco", 
              add = "jitter") 
p <- p + stat_compare_means()   #  Add p-value 
p   #得出预测结果

结果展示了不同生存状态下患者的riskScore值，其中0代表存活，1代表死亡。我们可以发现，与存活的患者相比，死亡患者的riskScore值显著升高，且其P值小于0.001。

###6.判断预测结果的准确性
library(ROCR)   #使用ROCR包绘制预测模型的ROC曲线 
library(glmnet) 
library(caret)

pred <- prediction(dat$riskScore, dat$OS) 
perf <- performance(pred,"tpr","fpr") 
performance(pred,"auc")   # shows calculated AUC for model 
plot(perf,colorize=FALSE, col="red")   #绘制ROC曲线 
lines(c(0,1),c(0,1),col = "gray", lty = 4 )























x=model.matrix(data[,4:29])
x=data[,4:29]
y=Surv(data$OS.time,data$OS)
## 切分为训练集核测试集，70%训练，30%测试
set.seed(123)
d_index <- createDataPartition(y,p = 0.7)
train_d <- x[d_index$Resample1,]
test_d <- x[-d_index$Resample1,]

## 数据标准化
scal <- preProcess(train_d,method = c("center","scale"))
train_ds <- predict(scal,train_d)
test_ds <- predict(scal,test_d)

X <- as.matrix(train_ds)
Y =y[d_index$Resample1,]
set.seed(1245)
lasso_model <- cv.glmnet(X,Y,alpha = 1,lambda = lambdas,nfolds =3)
plot(lasso_model)
plot(lasso_model$glmnet.fit, "lambda", label = T)