《机器学习实用案例解析》((美)DrewConwayJohnMylesWhite)电子书网盘下载免费在线阅读
链接: https://pan.baidu.com/s/1Tlmx4cyoX7iioA6X-rbj1Q
提取码: gh9a书名:机器学习实用案例解析
豆瓣评分:7.6
作者:(美)DrewConwayJohnMylesWhite
出版社: 机械工业出版社
副标题:实用案例解析
原作名: Machine Learning for Hackers
译者: 陈开江 / 刘逸哲 / 孟晓楠 / 罗森林 审校
出版年: 2013-4-1
页数: 320
内容简介
这本书为机器学习技术提供了一些非常棒的案例研究。它并不想成为一本关于机器学习的工具书或者理论书籍,它注重的是一个学习的过程,因而对于任何有一些编程背景和定量思维的人来说,它都是不错的选择。
——Max Shron OkCupid
机器学习是计算机科学和人工智能中非常重要的一个研究领域,近年来,机器学习不但在计算机科学的众多领域中大显身手,而且成为一些交叉学科的重要支撑技术。本书比较全面系统地介绍了机器学习的方法和技术,不仅详细阐述了许多经典的学习方法,还讨论了一些有生命力的新理论、新方法。
全书案例既有分类问题,也有回归问题;既包含监督学习,也涵盖无监督学习。本书讨论的案例从分类讲到回归,然后讨论了聚类、降维、最优化问题等。这些案例包括分类:垃圾邮件识别,排序:智能收件箱,回归模型:预测网页访问量,正则化:文本回归,最优化:密码破解,无监督学习:构建股票市场指数,空间相似度:用投票记录对美国参议员聚类,推荐系统:给用户推荐R语言包,社交网络分析:在Twitter上感兴趣的人,模型比较:给你的问题找到最佳算法。各章对原理的叙述力求概念清晰、表达准确,突出理论联系实际,富有启发性,易于理解。在探索这些案例的过程中用到的基本工具就是R统计编程语言。R语言非常适合用于机器学习的案例研究,因为它是一种用于数据分析的高水平、功能性脚本语言。
本书主要内容:
·开发一个朴素贝叶斯分类器,仅仅根据邮件的文本信息来判断这封邮件是否是垃圾邮件;
·使用线性回归来预测互联网排名前1000网站的PV;
·利用文本回归理解图书中词与词之间的关系;
·通过尝试破译一个简单的密码来学习优化技术;
·利用无监督学习构建股票市场指数,用于衡量整体市场行情的好坏;
·根据美国参议院的投票情况,从统计学的角度对美国参议员聚类;
·通过K近邻算法构建向用户推荐R语言包;
·利用Twitter数据来构建一个“你可能感兴趣的人”的推荐系统;
·模型比较:给你的问题找到最佳算法。
作者简介
【作者介绍】
Drew Conway 机器学习专家,拥有丰富的数据分析与处理工作经验。目前主要利用数学、统计学和计算机技术研究国际关系、冲突和恐怖主义等。他曾作为研究员在美国情报和国防部门供职数年。他拥有纽约大学政治系博士学位,曾为多种杂志撰写文章,是机器学习领域的著名学者。
John Myles White 机器学习专家,拥有丰富的数据分析与处理工作经验。目前主要从理论和实验的角度来研究人类如何做出决定,同时还是几个流行的R语言程序包的主要维护者,包括ProjectTemplate和log4r。他拥有普林斯顿大学哲学系博士学位,曾为多家技术杂志撰稿,发表过许多关于机器学习的论文,并在众多国际会议上发表演讲。
【译者介绍】
罗森林 博士,教授,博导。现任北京理工大学信息系统及安全对抗实验中心主任、专业责任教授。国防科技工业局科学技术委员会成员;《中国医学影像技术杂志》、《中国介入影像与治疗学》编委会委员;全国大学生信息安全技术专题邀请赛专家组副组长;中国人工智能学会智能信息安全专业委员会委员等。主要研究方向为信息安全、数据挖掘、媒体计算、中文信息处理等。负责或参加完成国家自然科学基金、国家科技支撑计划、863计划、国家242计划等省部级以上项目40余项。已发表学术论文90余篇,出版著作8部,出版译著1部,获授权专利3项。
陈开江 新浪微博搜索部研发工程师,曾独立负责微博内容反垃圾系统、微博精选内容挖掘算法、自助客服系统(包括自动回复、主动挖掘、舆情监测)等项目,目前主要从事社交挖掘、推荐算法研究、机器学习、自然语言处理相关工作,研究兴趣是社交网络的个性化推荐。
刘逸哲 阿里巴巴,CBU基础平台部搜索与推荐团队核心技术与query分析方向负责人,机器学习技术领域及圈子负责人。曾任中国雅虎相关性团队、自然语言处理团队算法工程师;AvePoint.inc开发工程师,从事企业级搜索引擎开发。研究兴趣是机器学习、自然语言处理及个性化推荐等算法在大规模数据上的应用。
孟晓楠 一淘广告技术,阿里非搜索广告算法负责人,负责用户行为分析、建模与细分,RTB竞价算法,展示广告CTR预估与SEM优化。曾工作于网易杭州研究院,参与过分布式全文检索系统和网易博客产品的数据挖掘算法开发。研究兴趣是计算广告技术、机器学习、大数据技术、信息检索等。
下游分析
cellranger count 计算的结果只能作为错略观测的结果,如果需要进一步分析聚类细胞,还需要进行下游分析,这里使用官方推荐 R 包(Seurat 3.0)
流程参考官方外周血分析标准流程( https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html )
Rstudio操作过程:
## 安装seurat
install.packages('Seurat')
## 载入seurat包
library(dplyr)
library(Seurat)
## 读入pbmc数据(文件夹路径不能包含中文,注意“/“的方向不能错误,这里读取的是10x处理的文件,也可以处理其它矩阵文件,具体怎样操作现在还不知道,文件夹中的3个文件分别是:barcodes.tsv,genes.tsv,matrix.mtx,文件的名字不能错,否则读取不到)
pbmc.data <- Read10X(data.dir = "D:/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
## 查看稀疏矩阵的维度,即基因数和细胞数
dim(pbmc.data)
pbmc.data[1:10,1:6]
## 创建Seurat对象与数据过滤,除去一些质量差的细胞(这里读取的是单细胞 count 结果中的矩阵目录;在对象生成的过程中,做了初步的过滤;留下所有在>=3 个细胞中表达的基因 min.cells = 3;留下所有检测到>=200 个基因的细胞 min.genes = 200。)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
##计算每个细胞的线粒体基因转录本数的百分比(%),使用[[ ]] 操作符存放到metadata中,mit-开头的为线粒体基因
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
##展示基因及线粒体百分比(这里将其进行标记并统计其分布频率,"nFeature_RNA"为基因数,"nCount_RNA"为细胞数,"percent.mt"为线粒体占比)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
## 过滤细胞:根据上面小提琴图中基因数"nFeature_RNA"和线粒体数"percent.mt",分别设置过滤参数,这里基因数 200-2500,线粒体百分比为小于 5%,保留gene数大于200小于2500的细胞;目的是去掉空GEMs和1个GEMs包含2个以上细胞的数据;而保留线粒体基因的转录本数低于5%的细胞,为了过滤掉死细胞等低质量的细胞数据。
pbmc <- subset(pbmc, subset = nFeature_RNA >200 &nFeature_RNA <2500 &percent.mt <5)
## 表达量数据标准化,LogNormalize的算法:A = log( 1 + ( UMIA ÷ UMITotal ) × 10000
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#pbmc <- NormalizeData(pbmc) 或者用默认的
## 鉴定表达高变基因(2000个),用于下游分析,如PCA;
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
## 提取表达量变化最高的10个基因;
top10 <- head(VariableFeatures(pbmc), 10)
top10
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10)
CombinePlots(plots = list(plot1, plot2))
plot1<-VariableFeaturePlot(object=pbmc)
plot2<-LabelPoints(plot=plot1,points=top10,repel=TRUE)
CombinePlots(plots=list(plot1,plot2))
## PCA分析:
# PCA分析数据准备,使用ScaleData()进行数据归一化;默认只是标准化高变基因(2000个),速度更快,不影响PCA和分群,但影响热图的绘制。
#pbmc <- ScaleData(pbmc,vars.to.regress ="percent.mt")
## 而对所有基因进行标准化的方法如下:
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
## 线性降维(PCA),默认用高变基因集,但也可通过features参数自己指定;
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
## 展示 pca 结果(最简单的方法)
DimPlot(object=pbmc,reduction="pca")
## 检查PCA分群结果, 这里只展示前5个PC,每个PC只显示5个基因;
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
##PC_ 1
##Positive: RPS27, MALAT1, RPS6, RPS12, RPL13
##Negative: CSTA, FCN1, CST3, LYZ, LGALS2
##PC_ 2
##Positive: NKG7, GZMA, CST7, KLRD1, CCL5
##Negative: RPL34, RPL32, RPL13, RPL39, LTB
##PC_ 3
##Positive: MS4A1, CD79A, BANK1, IGHD, CD79B
##Negative: IL7R, RPL34, S100A12, VCAN, AIF1
##PC_ 4
##Positive: RPS18, RPL39, RPS27, MALAT1, RPS8
##Negative: PPBP, PF4, GNG11, SDPR, TUBB1
##PC_ 5
##Positive: PLD4, FCER1A, LILRA4, SERPINF1, LRRC26
##Negative: MS4A1, CD79A, LINC00926, IGHD, FCER2
## 展示主成分基因分值
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
## 绘制pca散点图
DimPlot(pbmc, reduction = "pca")
## 画第1个或15个主成分的热图;
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
## 确定数据集的分群个数
# 鉴定数据集的可用维度,方法1:Jackstraw置换检验算法;重复取样(原数据的1%),重跑PCA,鉴定p-value较小的PC;计算‘null distribution’(即零假设成立时)时的基因scores。虚线以上的为可用维度,也可以调整 dims 参数,画出所有 pca 查看。
#pbmc <- JackStraw(pbmc, num.replicate = 100)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
#JackStrawPlot(pbmc, dims = 1:15)
# 方法2:肘部图(碎石图),基于每个主成分对方差解释率的排名。
ElbowPlot(pbmc)
## 细胞聚类:分群个数这里选择10,建议尝试选择多个主成分个数做下游分析,对整体影响不大;在选择此参数时,建议选择偏高的数字(为了获取更多的稀有分群,“宁滥勿缺”);有些亚群很罕见,如果没有先验知识,很难将这种大小的数据集与背景噪声区分开来。
## 非线性降维(UMAP/tSNE)基于PCA空间中的欧氏距离计算nearest neighbor graph,优化任意两个细胞间的距离权重(输入上一步得到的PC维数) 。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
## 接着优化模型,resolution参数决定下游聚类分析得到的分群数,对于3K左右的细胞,设为0.4-1.2 能得到较好的结果(官方说明);如果数据量增大,该参数也应该适当增大。
pbmc <- FindClusters(pbmc, resolution = 0.5)
## 使用Idents()函数可查看不同细胞的分群;
head(Idents(pbmc), 5)
## 结果:AAACCTGAGGTGCTAG AAACCTGCAGGTCCAC AAACCTGCATGGAATA AAACCTGCATGGTAGG AAACCTGCATTGGCGC
1 3 0 10 2
Levels: 0 1 2 3 4 5 6 7 8 9 10 11
## Seurat提供了几种非线性降维的方法进行数据可视化(在低维空间把相似的细胞聚在一起),比如UMAP和t-SNE,运行UMAP需要先安装'umap-learn'包,这里不做介绍,两种方法都可以使用,但不要混用,如果混用,后面的结算结果会将先前的聚类覆盖掉,只能保留一个。
## 这里采用基于TSNE的聚类方法。
pbmc <- RunTSNE(pbmc, dims = 1:10)
## 用DimPlot()函数绘制散点图,reduction = "tsne",指定绘制类型;如果不指定,默认先从搜索 umap,然后 tsne, 再然后 pca;也可以直接使用这3个函数PCAPlot()、TSNEPlot()、UMAPPlot(); cols,pt.size分别调整分组颜色和点的大小;
DimPlot(pbmc,reduction = "tsne",label = TRUE,pt.size = 1.5)
## 这里采用基于图论的聚类方法
pbmc<-RunUMAP(object=pbmc,dims=1:10)
DimPlot(object=pbmc,reduction="umap")
## 细胞周期归类
pbmc<- CellCycleScoring(object = pbmc, g2m.features = cc.genes$g2m.genes, s.features = cc.genes$s.genes)
head(x = pbmc@meta.data)
DimPlot(pbmc,reduction = "tsne",label = TRUE,group.by="Phase",pt.size = 1.5)
## 存储结果
saveRDS(pbmc, file = "D:/pbmc_tutorial.rds")
save(pbmc,file="D:/res0.5.Robj")
## 寻找cluster 1的marker
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
## 结果: p_val avg_logFC pct.1 pct.2 p_val_adj
MT-CO1 0.000000e+00 -0.6977083 0.985 0.996 0.000000e+00
RPS27 2.182766e-282 0.3076454 1.000 0.999 3.480202e-278
MT-CO3 2.146399e-274 -0.4866429 0.995 0.997 3.422218e-270
DUSP1 2.080878e-247 -1.7621662 0.376 0.745 3.317752e-243
RPL34 8.647733e-244 0.3367755 1.000 0.997 1.378795e-239
##寻找每一cluster的marker
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
# A tibble: 24 x 7
# Groups: cluster [12]
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl><dbl><dbl> <dbl><fct> <chr>
1 2.29e-123 0.636 0.344 0.097 3.65e-119 0 CD8B
2 7.62e-113 0.487 0.632 0.305 1.22e-108 0 LEF1
3 2.04e- 74 0.483 0.562 0.328 3.25e- 70 1 LEF1
4 1.39e- 61 0.462 0.598 0.39 2.22e- 57 1 ITM2A
5 0. 2.69 0.972 0.483 0. 2 GNLY
6 0. 2.40 0.964 0.164 0. 2 GZMB
7 1.31e-121 0.768 0.913 0.671 2.09e-117 3 JUNB
8 2.06e- 94 0.946 0.426 0.155 3.28e- 90 3 RGS1
9 2.05e-255 1.57 0.586 0.09 3.27e-251 4 GZMK
10 2.94e-140 1.57 0.69 0.253 4.68e-136 4 KLRB1
# ... with 14 more rows
## 存储marker
write.table(pbmc.markers,file="D:/allmarker.txt")
## 各种绘图
## 绘制Marker 基因的tsne图
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"),cols = c("gray", "red"))
## 绘制Marker 基因的小提琴图
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
## 绘制分cluster的热图
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
剩下的便是寻找基因 marker 并对细胞类型进行注释(见下回分解)
今天被粉丝发的文章给难住了,又偷偷去学习了一下竞争风险模型,想起之前写的关于竞争风险模型的做法,真的都是皮毛哟,大家见笑了。想着就顺便把所有的生存分析的知识和R语言的做法和论文报告方法都给大家梳理一遍。什么时候用生存分析
当你关心结局和结局发生时间的时候,就要考虑生存分析了,这种既有结局又有时间的数据叫做生存数据,英文叫做Time-to-event data. 只不过因为这个方法医学上用来分析存活情况用的多,所以得名生存分析,反正你就记住一个例子,我要研究汽车发生故障,我也应该用生存分析,因为我既关心是不是有故障,我还关心用了多久(跑了多远)才出故障,就是既有time,又有event,Time-to-event data就用生存分析。
基本概念
首先是删失,对象失访了,脱落了,出现结局之前随访结束了,都叫做删失:
删失又分为左删失,区间删失和右删失,图示如下:
比如我想研究得了A病的人的生存情况,存在的所有可能情形为:
第一种,研究的开始的时候有人已经有A病,这个时候人家已经活了一段时间了,具体多久我不知道,叫做左删失;
第二种,入组随访的时候没病,中途得了A病死了,什么时候得的,没记录下来,叫区间删失;
第三种,得了A病,一直活到了研究结束还没死,叫做右删失。
你看,所有的删失情况造成的后果都是我们没法准确估计发生结局的时间,这也是其名字删失的由来,对于这类数据就需记录为删失数据。
生存分析的种类有哪些
具体的种类是为了回答具体的问题,我们做生存分析常常要回答的问题如下:
一个是描述生存情况,一个是比较,再一个就是探究影响因素。
比如我随访了很多病人,我就想知道随着时间变化这群人的生存概率是如何变化的(描述)?我就可以用简单粗暴的Kaplan-Meier method,又叫乘积极限法,基本思想就是此刻的生存概率等于上一时刻的生存概率乘以此刻的存活率。
比如我手上有如下数据:
time是随访时间,status是结局,我就可以写出如下代码拟合出总体人群的生存曲线:
fit1 <- survfit(Surv(time, status) ~ 1, data = mydata)summary(fit1)
并且得到每个时间点,病人生存概率的估计值和标准误:
如果我还恰好收集了病人的性别,我又想看看男病人和女病人生存情况是不是有不同(比较),我可以对生存曲线分组比较,代码如下:
fit2<- survfit(Surv(time, status) ~ sex, data = mydata)summary(fit2)
输出两组生存曲线比较的log-rank test的统计量:
survdiff(Surv(time, status) ~ sex, data = mydata)
并且我们还可以进行复杂分组的组间比较,大家可以翻看我之前的文章。
但是大家明白,KM法始终是在做单因素分析,而且都是做的分类变量的单因素分析,我们经常还会有的需求是考虑各种混杂的情况下去探讨影响生存时间的因素,这个时候我们就要用到The Cox regression model了。
模型形式如下:
上面的式子把h0移项到左边,等号两边同时取对数就成了一个线性模型,和广义线性模型的理解一样一样的,b就是我们的影响因素的系数,解释为lnHR的改变量,其为正就是危险因素,为负就是保护因素:
The quantities exp(bi) are called hazard ratios (HR). A value of bi greater than zero, or equivalently a hazard ratio greater than one, indicates that as the value of the ith covariate increases, the event hazard increases and thus the length of survival decreases.
再观测一下上面的式子,我们拟合了预测因素对风险比例(t时刻风险比上基础风险)的模型,这个时候暗含的假设就是就是对于每个人在任何时刻风险只有一个常数,就是在所有预测因素的作用下,各个时间的风险是不变的,这个就叫做Proportional Hazards Assumption, 比例风险假设。
做COX比例风险模型的时候都有必要验证这个假设是不是满足,具体方法如下:
我们需要先做一个cox模型,拟合cox模型需要用到coxph函数,代码如下:
res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data = lung)res.cox
模型输出结果里面会有预测变量的β值(coef)和其标准误,有HR(exp(coef))值,还有预测变量的显著性检验结果:
我们将这个模型对象直接喂给cox.zph,就可以得到风险比例假设的检验结果了:
test.ph <- cox.zph(res.cox)test.ph
可以看到,我们的p值都大于0.05,即都满足ph假设。
我们还可以从图上看,看scaled Schoenfeld residuals是不是和时间独立,如果独立则满足ph假设:
ggcoxzph(test.ph)
有竞争事件时
上面写了只有一个截距事件时生存分析的不同目的,描述,比较,和影响因素探讨,接着再来看存在竞争事件的时候该如何描述,比较,和探讨影响因素。
生存分析的另外情况就是竞争风险模型,就是time to event的event有多种,一种发生了另外一种就不可能发生了,这个就是竞争,比如好多文献中都会列举的经典例子:就是在造血干细胞移植人群中,我们关心其疾病复发风险,但是有些人因为移植死了(TRM),死了之后肯定也谈不上复发,如果你把因为移植死了的人都作为删失数据,肯定也是不对的,这个里面就会有两个竞争结局相关性造成的问题,此时我们应该用竞争风险模型来分析。
For example, disease relapse is an event of interest in studies of allogeneic hematopoietic stem cell transplantation (HSCT), as is mortality related to complications of transplantation (transplant-related mortality or TRM). Relapse and TRM are not independent in this setting because these two events are likely related to immunologic effector mechanisms following HSCT, whereby efforts to reduce TRM may adversely affect the risk of relapsemoreover, patients who die from TRM cannot be at further risk of relapse.
在比例风险中我们的结局事件只有一个,我们关心的是哪些因素对事件发生的风险有影响。在竞争风险模型中我们关注的地方又变了,因为我们有好几个结局事件,这个时候我们会关心在竞争事件存在的情况下各个结局事件的累计发病函数是如何随着时间变化的,以及如何来比较不同的累计发病函数,以及如何探讨影响因素(competing risks regression analysis)。
之前写的在非竞争风险模型中累计发病率的组间比较可以用KM法,将纵坐标换一换,用log-rank检验,在竞争风险模型中我们需要用Fine and Gray提到的方法来做(Gray’s test),如果比如说我发现两组(治疗组和对照组)的累计发病风险不一样,我肯定还想探讨哪些因素影响累计发病风险,之前是用COX比例风险模型做的,在竞争结局存在的情况下我们依然是得用Fine and Gray提出的模型(Fine and Gray Model),叫做竞争风险回归模型:
Fine and Gray (6) proposed a method for direct regression modeling of the effect of covariates on the cumulative incidence function for competing risks data. As in any other regression analysis, modeling cumulative incidence functions for competing risks can be used to identify potential prognostic factors for a particular failure in the presence of competing risks, or to assess a prognostic factor of interest after adjusting for other potential risk factors in the model.
首先看竞争风险时候累计发病曲线的比较方法。我手上有数据如下:
其中dis是疾病类型,2分类,ftime是时间,status是结局事件,结局事件有3个水平,多的1个水平是竞争事件。现在我关心不同疾病人群各个事件累积发病曲线有无不同,我可以用cuminc函数结合plot画出各个组的累计发病曲线:
fit=cuminc (ftime, status, dis, cencode = 0) plot(fit)
fit对象的结果中还有每条曲线的比较,从比较结果可以看出两组在事件2的累积发病曲线上是有显著差异的:
上面介绍的方法相当于没有竞争风险的时候的KM法,通过上面的方法我们知道有了不同风险结局累计发病曲线的差异,继续我们会继续看影响因素,要做的就是竞争风险回归模型了,需要用到的函数就是crr。
比如我手上有如下数据,除了时间,结局还包括每个病人的像sex,age等等协变量,我想探讨说这些因素是如何影响病人某个结局的,我就可以写出一个竞争风险回归模型:
mod1 <- crr(ftime,Status,x)summary(mod1)
上面的代码中x是自变量矩阵,运行代码输出竞争风险回归模型的结果如下:
到这儿基本就写完了所有的生存分析的情况,接着再结合一篇论文看看报告方法,论文就是下面这篇,是我随意查的:
S. Chen, H. Sun, M. Heng, X. Tong, P. Geldsetzer, Z. Wang, P. Wu, J. Yang, Y. Hu, C.
Wang, T. Bärnighausen, Factors Predicting Progression to Severe COVID-19: A Competing Risk Survival Analysis of 1753 Patients in Community Isolation in Wuhan, China, Engineering (2021), doi: https://doi.org/10.1016/j.eng.2021.07.021
这个作者用竞争风险模型探讨了重型新冠进程的影响因素,作者报告了每个影响因素的HR和置信区间,p值,这些都很容易做出来。还报告了固定时间点的累积发病的置信区间。在R语言中固定时间点的累计发病置信区间可以用CumIncidence这个函数计算出来:
The CumIncidence() function allows for the pointwise confidence intervals, by simply adding a further argument, level, where we specify the desired confidence level.
比如我想计算第10天各组的累计发病风险的置信区间,我就可以写出如下代码:
CumIncidence (ftime, status, dis, cencode = 0,t=10,level = 0.95)
输出结果如下:
图中就是各个组在时间为10的时候结局累计发病风险的置信区间。
小结
也欢迎大家的意见和建议,大家想了解什么统计方法都可以在文章下留言,说不定我看见了就会给你写教程哦,有疑问欢迎私信。
如果你是一个大学本科生或研究生,如果你正在因为你的统计作业、数据分析、模型构建,科研统计设计等发愁,如果你在使用SPSS, R,Python,Mplus, Excel中遇到任何问题,都可以联系我。因为我可以给您提供最好的,最详细和耐心的数据分析服务。
如果你对Z检验,t检验,方差分析,多元方差分析,回归,卡方检验,相关,多水平模型,结构方程模型,中介调节,量表信效度等等统计技巧有任何问题,请私信我,获取详细和耐心的指导。
If you are a student and you are worried about you statistical #Assignments, #Data #Analysis, #Thesis, #Reports, #Composing, #Quizzes, Exams.. And if you are facing problem in #SPSS, #R-Programming, #Excel, Mplus, then contact me. Because I could provide you the best services for your Data Analysis.
Are you confused with statistical Techniques like z-test, t-test, ANOVA, MANOVA, Regression, Logistic Regression, Chi-Square, Correlation, Association, SEM, multilevel model, mediation and moderation etc. for your Data Analysis...??
Then Contact Me. I will solve your Problem...
加油吧,打工人!
R数据分析:用R语言做潜类别分析LCA
R数据分析:ROC曲线与模型评价实例
R数据分析:用R语言做meta分析
R数据分析:贝叶斯定理的R语言模拟
R数据分析:什么是人群归因分数Population Attributable Fraction
R数据分析:有调节的中介
R数据分析:如何用R做多重插补,实例操练
R数据分析:如何用R做验证性因子分析及画图,实例操练
R数据分析:手把手教你画列线图(Nomogram)及解读结果
R数据分析:倾向性评分匹配完整实例(R实现)
R文本挖掘:文本聚类分析
R数据分析:混合效应模型实例
R文本挖掘:中文文本聚类
R文本挖掘:中文词云生成
R数据分析:临床预测模型的样本量探讨及R实现
R数据分析:多分类逻辑回归
R文本挖掘:社会网络分析
R数据分析:中介效应的做法
R数据分析:随机截距交叉滞后RI-CLPM与传统交叉滞后CLPM
R数据分析:多水平模型详细说明
R数据分析:如何做潜在剖面分析Mplus
R数据分析:生存分析的做法与解释续
R数据分析:竞争风险模型的做法和解释二
R数据分析:多元逻辑斯蒂回归的做法
R数据分析:线性回归的做法和优化实例
R数据分析:双分类变量的交互作用作图
R机器学习:基于树的分类算法的原理及实现
R数据分析:如何做聚类分析,实操解析
R数据分析:潜在剖面分析LPA的做法与解释
欢迎分享,转载请注明来源:夏雨云
评论列表(0条)