【原】KM生存曲线经logRNA检验后也可以计算HR值

号外：中秋节广州3天入门课程报名马上截止：（中秋节一起来学习！）全国巡讲第16站-广州（生信入门课加量不加价）

最近根据基因表达量对病人进行分组后，使用KM生存分析的logRNAK法来检查两个组的病人的生存差异，得到了如下的图：

可以很明显看到，根据基因表达量把病人分成高表达组合低表达组，经过log rank 检验，可以看到两组病人的生存是有统计学显著差异的，而且我们可以看到，基因表达量越高，病人死亡风险越高，所以我们可以把这个基因在表达水平上看做是风险因子，而我们通常判断风险因子或者保护因子是根据hazard ratio来的，可以看官方教材：https://www./topics/medicine-and-dentistry/hazard-ratio

As for the other measures of association, a hazard ratio of 1 means lack of association, a hazard ratio greater than 1 suggests an increased risk, and a hazard ratio below 1 suggests a smaller risk.

所以我就突发奇想，是不是log rank 检验后也可以返回hazard ratio呢？

下面是KM生存分析和log rank检验的具体定义：

KM法即乘积极限法（product-limit method），是现在生存分析最常用的方法，是由Kaplan和Meier于1958年提出，因此称Kaplan-Meier法，通常简称KM法。KM法是这样估计生存曲线：首先计算出活过一定时期的病人再活过下一时期的概率（即生存概率），然后将逐个生存概率相乘，即为相应时段的生存率。

如果有两组病人的的生存曲线有所不同，但生存率差别是否有统计学意义呢，这可以用log rank test(时序检验或对数秩检验)来检验。log rank test是计算不同日期两种（或多种）疗法的暴露人数及出现终点人数，计算不同时期期望人数与实际出现终点的差值，以此可作卡方检验作出判断。当P<0.05，认为两组或多组总体生存曲线差别有统计学意义。

简单的谷歌搜索，就找到了答案，在 https://stats./questions/124489/how-to-calculate-the-hr-and-95ci-using-the-log-rank-test-in-r  教程里面，作者 是根据书本的公式自己写的代码：According to the book Survival Analysis: A Practical Approach, I got two formulas on Page 62 and 66 to do this (as shown below). So I wrote the R code as below, is there anybody know whether I'm right?

有趣的是他也不确定是否正确，所以我找到的链接其实是他在向广大网友求助：

使用的是R语言，代码比较长：

library(survival)
data.survdiff <- survdiff(Surv(time, status) ~ group)
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))


我使用公式计算后发现，他计算的是基因低表达量组的hazard ratio值，也就是说，如果hazard ratio值小于1，表明这个低表达是保护因子，所以这个基因高表达就是风险因子，是不是听起来有点绕口，这个其实就有点类似于我们的差异表达分析，我们使用limma等R包计算差异基因的上调下调也是遇到同样的困境，case相当于control的上调基因，其实就是control相当于case的下调基因。

另外，附上生存分析需理解的的定义：

1. 生存分析，是将终点事件出现与否与对应时间结合起来分析的一种统计方法；

2. 生存时间，是从规定的观察起点到某一特定终点事件出现的时间，如膀胱癌术后5年存活率研究，及膀胱癌手术为观测起点，死亡为事件终点，两点为生存时间；

3. 完全数据，观测起点到终点事件所经历的时间，上述例子即膀胱癌手术到因膀胱癌死亡的时间；

4. 删失数据，因失访、研究结束终点事件未发生或患者死于规定的终点事件以外的原因而终止观察，不能确定具体生存时间的一类数据；

5. 生存概率，表示某时段开始存活的个体到该时段结束仍存活的概率，p=活满某时段的人数/该时段期初有效人口数；

6. 生存率，为观察起点起到研究时间点内各个时段的生存概率的累积概率，S(tk)=p1.p2.pk=S(tk-1).pk;

7. 生存曲线，以生存时间为横轴，将各个时间点的生存率连在一起的曲线图；

8. 中位生存期，又称半数生存期，表示50%的个体存活的时间；

9. PH假定（等比例风险假定），某研究因素对生存的影响不随时间的改变而改变，是COX回归模型建立的前提条件。

实际上，没有必要这样做，因为coxph本身 就可以返回HR值，这里我不得不指出很久以前的批量生存分析教程的一个小瑕疵：

都可以批量做生存分析了，还要网页工具干嘛？

批量COX回归生存分析图，指定挑选lncRNA基因，森林图，ROC曲线打包给你

如果你不理解代码，仅仅是跑所谓高手/大佬的代码，就会出现南辕北辙的后果，无论如何都需要用自己的理解力来判断代码是否适用于你的情况，上下调搞反了，风险因子或者保护因子也会搞混！