机器学习之随机森林
2020-04-30 本文已影响0人
猪莎爱学习
开篇先看个风险森林图吧~~
image.png
1.准备输入数据
load("TCGA-KIRCsur_model.Rdata")
2.挑选感兴趣的基因构建coxph模型
出自文章Integrated genomic analysis identifies subclasses and prognosis signatures of kidney cancer中，五个miRNA是miR-21,miR-143,miR-10b,miR-192,miR-183
将他们从表达矩阵中取出来，就得到了5个基因在522个肿瘤样本中的表达量，可作为列添加在meta表噶的后面,组成的数据框赋值给dat。
#首先把这5个基因的表达量挑选出来，一行是一个基因，一列是一个样本，根据表达矩阵按行取子集[,]
e=t(exprSet[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
e=log2(e)
#补充一个小小知识点，就是把上面的中划线替换成下划线
#library(stringr)
#str_replace("hsa-mir-21","-","_")
#上面是替换一个，如果是要替换所有，可以加上all
#str_replace_all("hsa-mir-21","-","_")
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183') #改名
dat=cbind(meta,e)
dat$gender=factor(dat$gender)
dat$stage=factor(dat$stage)
colnames(dat)
用survival::coxph()函数构建模型
library(survival)
library(survminer)
s=Surv(time, event) ~ gender + stage + age + miR21+miR143+miR10b+miR192+miR183 #越往前就是越基础，gender,stage,age然后才是miR，不要乱改代码
#s=Surv(time, event) ~ miR21+miR143+miR10b+miR192+miR183
model <- coxph(s, data = dat )
3.模型可视化-森林图
ggforest 风险森林图
options(scipen=1)
ggforest(model, data =dat,
main = "Hazard ratio",
cpositions = c(0.10, 0.22, 0.4),
fontsize = 1.0,
refLabel = "1", noDigits = 4)
4.模型预测
fp <- predict(model)
summary(model)$concordance
library(Hmisc)
options(scipen=200)
with(dat,rcorr.cens(fp,Surv(time, event)))
# 若要找到最佳模型，我们可以进行变量选择，可以采用逐步回归法进行分析
这里只是举个栗子，自己预测自己的C-index是1-with(dat,rcorr.cens(fp,Surv(time, event)))[[1]],实战应该是拿另一个数据集来预测，或者将一个数据集分两半，一半构建模型，一半验证,可以使用机器学习的R包切割数据。
C-index用于计算生存分析中的COX模型预测值与真实之间的区分度（discrimination），也称为Harrell's concordanceindex。C-index在0.5-1之间。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。
5.切割数据构建模型并预测
5.1 切割数据
用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。
library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)
head(sam)
可查看两组一些临床参数切割比例
train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]
prop.table(table(train_meta$stage))
prop.table(table(test_meta$stage))
prop.table(table(test_meta$race))
prop.table(table(train_meta$race))
5.2 切割后的train数据集建模
和上面的建模方法一样。
e=t(train[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
e=log2(e)
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183')
dat=cbind(train_meta,e)
dat$gender=factor(dat$gender)
dat$stage=factor(dat$stage)
colnames(dat)
s=Surv(time, event) ~ gender + stage + age + miR21+miR143+miR10b+miR192+miR183
#s=Surv(time, event) ~ miR21+miR143+miR10b+miR192+miR183
model <- coxph(s, data = dat )
5.3 模型可视化
options(scipen=1)
ggforest(model, data =dat,
main = "Hazard ratio",
cpositions = c(0.10, 0.22, 0.4),
fontsize = 1.0,
refLabel = "1", noDigits = 4)
5.4 用切割后的数据test数据集验证模型
e=t(test[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
e=log2(e)
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183')
test_dat=cbind(test_meta,e)
fp <- predict(model)
summary(model)$concordance
library(Hmisc)
options(scipen=200)
with(test_dat,rcorr.cens(fp,Surv(time, event)))
C-index不到0.6， 模型全靠猜了。
内部评价：训练集
外部评价：测试集
在没有拆分数据的时候我们只有训练集，就自己预测自己
莎莎写于2020.4.30 今天是噶宝生日~🎂
1.准备输入数据
输入数据是肿瘤样本表达矩阵exprSet、临床信息meta
load("sur_model_KIRC.Rdata")
library(randomForest)
library(ROCR)
library(genefilter)
library(Hmisc)
2.构建随机森林模型
输入数据是表达矩阵(仅含tumor样本)和对应的生死。
x=t(log2(exprSet+1))
y=meta$event
#构建模型，一个叫randomForest的函数,运行时间很长,存Rdata跳过
tmp_rf="TCGA_KIRC_miRNA_rf_output.Rdata"
if(!file.exists(tmp_rf)){
rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE ) #这一步运行要20多分钟，所以设置了跳过，如果再运行自己的数据的时候再修改
save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)
#top30的基因
varImpPlot(rf_output, type=2, n.var=30, scale=FALSE,
main="Variable Importance (Gini) for top 30 predictors",cex =.7)
rf_importances=importance(rf_output, scale=FALSE)
head(rf_importances)
#解释量top30的基因,和图上是一样的，从下往上看。
choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30)) #用tail函数取了后30个
image.png
image.png
image.png
这样就实现了倒序
3.模型预测和评估
3.1自己预测自己
rf.prob <- predict(rf_output, x)
re=cbind(y ,rf.prob)
head(re)
3.2 箱线图
对预测结果进行可视化。以实际的生死作为分组，画箱线图整体上查看预测结果。
re=as.data.frame(re)
colnames(re)=c('event','prob')
re$event=as.factor(re$event)
library(ggpubr)
p1 = ggboxplot(re, x = "event", y = "prob",
color = "event", palette = "jco",
add = "jitter")+ stat_compare_means()
p1
对比lasso回归，这个似乎更好一些？
3.3 ROC曲线
library(ROCR)
#library(caret)
pred <- prediction(re[,2], re[,1])
auc = performance(pred,"auc")@y.values[[1]]
auc
自己预测自己，auc值竟然是1。
4.切割数据构建模型并预测
4.1 切割数据
用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。
library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)
train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]
4.2 切割后的train数据集建模
和上面的建模方法一样。
#建模
x = t(log2(train+1))
y = train_meta$event
tmp_rf="TCGA_KIRC_miRNA_t_rf_output.Rdata"
if(!file.exists(tmp_rf)){
rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)
choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))
length(choose_gene)
5.模型预测
用训练集构建模型，预测测试集的生死。
x = t(log2(test+1))
y = test_meta$event
rf.prob <- predict(rf_output, x)
re=cbind(y ,rf.prob)
re=as.data.frame(re)
colnames(re)=c('event','prob')
re$event=as.factor(re$event)
library(ggpubr)
p1 = ggboxplot(re, x = "event", y = "prob",
color = "event", palette = "jco",
add = "jitter")+ stat_compare_means()
p1
再看AUC值。
library(ROCR)
# 训练集模型预测测试集
pred <- prediction(re[,2], re[,1])
auc= performance(pred,"auc")@y.values[[1]]
auc
还行。
心理感受：怎么说呢，这部分内容对我来说太抽象了，我不是很懂，具体哪里不懂也不知道，就是觉得很难懂，很难理解，后面也不知道怎么用，还有一个支持向量机，还有timeROC和风险因子评估的也没有搞懂，代码就不贴了，慢慢学习吧~~~😔
小洁老师画的图，真好看
