R语言绘制生存曲线95%区间

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了R语言绘制生存曲线95%区间相关的知识,希望对你有一定的参考价值。

参考技术A 1. 安装和加载包

绘制Kaplan-Meier生存曲线需要用到的R包:survminer和survival。

library(survminer) # 加载包

library(survival) # 加载包

2 拟合曲线

R中使用survfit()函数来拟合生存曲线。

fit.3<-survfit(Surv(住院天数+病程,组别)~cd1656,data=data)

3. 绘制曲线函数

ggsurvplot(fit, data = NULL, fun = NULL, color = NULL,

          palette = NULL, linetype = 1, conf.int = FALSE,

          pval = FALSE, pval.method = FALSE,

          test.for.trend = FALSE, surv.median.line = "none",

          risk.table = FALSE, cumevents = FALSE,

          cumcensor = FALSE, tables.height = 0.25,

          group.by = NULL, facet.by = NULL, add.all = FALSE,

          combine = FALSE, ggtheme = theme_survminer(),

          tables.theme = ggtheme, ...)

# 参数解释

fit # 拟合的生存曲线对象

data # 用来拟合生存曲线的数据集

fun  # 常用三个字符参数;

# "event"绘制累积事件(f(y)=1-y),

# "cumhaz"绘制累积危害函数(f(y)=-log(y));

# "pct"绘制生存概率(百分比)。

color # 设置生存曲线的颜色。

# 如果只有1条曲线,则直接设置color="blue";

# 如果有多条曲线,默认color="strata",按分组为生存曲线着色;

# 也可以自定义调色板来设置曲线颜色。

palette # 调色板,默认"hue"。

# 可选调色板有"grey","npg","aaas","lancet",

# "jco", "ucscgb","uchicago","simpsons"和"rickandmorty".

linetype = 1 # 设置曲线线型。可以按"strata"设置线型;

# 或按数字向量c(1, 2)或按字符向量c("solid", "dashed")设置

conf.int # 逻辑词;默认FASLE;为TRUE则绘制曲线置信区间

pval = FALSE # 逻辑词;为TRUE则将统计检验计算的p值添加到图上;

# 为数字,则直接指定P值大小,如pval = 0.03;

# 为字符串,则添加字符串到图上,如pval = "p-value: 0.031"

pval.method  # 逻辑词,是否添加计算p值的统计方法的文本;

# 只有当 pval = TRUE时, 才会在图上添加检验方法文本

test.for.trend # 逻辑词,默认为FALSE;

# 为TRUE则返回趋势p值的检验,趋势检验旨在检验生存曲线的有序差异

surv.median.line # 在中位生存时间点处绘制水平或垂直线的字符向量;

# 可用值有"none"、"hv"、"h"、"v";其中v绘制垂直线,h绘制水平线。

risk.table = FALSE  # 逻辑词,图上是否添加风险表;

# "absolute" 显示处于风险中的绝对数量;

# "percentage" 显示处于风险中的百分比数量

# "abs_pct" 显示处于风险中的绝对数量和百分比

cumevents # 逻辑词,是否添加累计事件表

cumcensor # 逻辑词,是否添加累计删失表

tables.height = 0.25 # 生存曲线图下所有生存表的高度,数值0-1之间

group.by  # 包含分组变量名称的字符向量,向量长度≤2

facet.by # 字符向量,指定绘制分面生存曲线的分组变量(应≤2)的名称

ggtheme=theme_survminer() # 设置ggplot2主题,如theme_bw()

tables.theme # 作用于生存表的ggplot2主题名称

# 有theme_survminer、theme_cleantable()

add.all = FALSE # 逻辑词;是否添加总患者生存曲线到主生存图中

MTTF 的置信区间 - R 中的 Weibull 生存曲线

【中文标题】MTTF 的置信区间 - R 中的 Weibull 生存曲线【英文标题】:confidence interval for MTTF - Weibull survival curve in R 【发布时间】:2015-07-27 09:59:38 【问题描述】:

我正在尝试在 R 中实施 Delta 方法来计算 Weibull 生存曲线的 MTTF 方差。形状参数是alpha,比例参数是delta。方差 = var;协方差 = cov.

等式是:

var(mttf) = var(alpha)*[d(mttf)/d(alpha)]^2 + 
2*cov(alpha,delta)*d(mttf)/d(alpha)*d(mttf)/d(delta)
 + var(delta)*[d(mttf/d(delta)]^2.    

地点:

d(mttf)/d(alpha) = gamma(1+1/delta)

d(mttf)/d(delta) = -alpha/delta^2 * gamma(1+1/delta) * digamma(1+1/delta)

所以方程变为:

var(mttf) = var(alpha)*[gamma(1+1/delta)]^2 +
 2*cov(alpha,delta)*gamma(1+1/delta)*(-alpha/delta^2 * gamma(1+1/delta) * digamma(1+1/delta))
 + var(delta)*[-alpha/delta^2 * gamma(1+1/delta) * digamma(1+1/delta)]^2

我可以从方差-协方差矩阵中获取 var(alpha)var(delta)cov(alpha,delta)

拟合的weibull模型称为ajust

vcov(ajust)
a=ajust$var[2,2]*ajust$scale^2
b=ajust$var[1,2]*ajust$scale
matriz=matrix(c(ajust$var[1,1],b,b,a),ncol=2,nrow=2)

var(alpha) = matriz[2,2]
var(delta) = matriz[1,1]
cov(alpha,delta) = matriz[1,2] or matriz[2,1]

还有更多

alpha=coef[2]
delta=coef[1]

其中 coef 是一个变量,其中包含来自 survreg 调整的参数 alpha 和 delta。

所以,计算 MTTF:

mttf<-coef[2]*(gamma((1+(1/coef[1]))))

并计算 mttf 方差:

var_mttf=matriz[2,2]*(gamma(1+1/coef[1]))^2+
2*matriz[1,2]*((-coef[2]/(coef[1]^2))*gamma(1+1/coef[1])*digamma(1+1/coef[1]))+
matriz[1,1]*((-coef[2]/(coef[1]^2))*gamma(1+1/coef[1])*digamma(1+1/coef[1]))^2

但不幸的是,我的 mttf 方差与我从互联网论文中获取的任何示例都不匹配。我修改了太多次了...

整个代码是:

require(survival)
require(stats)
require(gnlm)

time<-c(0.22,  0.5, 0.88,   1.00,   1.32,   1.33,   1.54,   1.76,   2.50,   3.00,   3.00,   3.00,   3.00)
cens<-c(1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  0,  0,  0)

#Weibull adjust with survreg
ajust<-survreg(Surv(time,cens)~1,dist='weibull')
alpha<-exp(ajust$coefficients[1])
beta<-1/ajust$scale

#Weibull coefficients
coef<-cbind(beta,alpha)

#MTTF
mttf<-coef[2]*(gamma((1+(1/coef[1]))))

#Data from variance-covariance matrix:
vcov(ajust)
a=ajust$var[2,2]*ajust$scale^2
b=ajust$var[1,2]*ajust$scale
matriz=matrix(c(ajust$var[1,1],b,b,a),ncol=2,nrow=2)

#MTTF variance - delta method
var_mttf=matriz[2,2]*(gamma(1+1/coef[1]))^2+
  2*matriz[1,2]*((-coef[2]/(coef[1]^2))*gamma(1+1/coef[1])*digamma(1+1/coef[1]))+
  matriz[1,1]*((-coef[2]/(coef[1]^2))*gamma(1+1/coef[1])*digamma(1+1/coef[1]))^2

#standard error - MTTF
se_mttf=sqrt(var_mttf)

#MTTF confidence intervall (95% confidence) 
upper=mttf+1.960*sqrt(var_mttf)
lower=mttf-1.960*sqrt(var_mttf)

所以,从我获取这些数据的论文中得出的结果是:

MTTF standard error = 0.47
MTTF upper = 2.98
MTTF lower = 1.15 

这与我的代码结果相差甚远。

但是纸上的 alphadelta 和 MTTF 与我的代码具有相同的值:

alpha = 2.273151
delta = 1.417457
MTTF = 2.067864

拜托,我想和你们分享这个困难,他们在 R 方面比我更有经验。

问候,维尼修斯。

【问题讨论】:

【参考方案1】:

我建议beta需要大于-1,但根据我自己的计算;贝塔 =2。

【讨论】:

虽然我计算出 beta 为 -4。因为它的 (delta) 系数与 beta vcov(ajust) alpha[2,1] delta[1,1] x^3-4/x^2+5 是向量 T-1 的潜在向量吗。 f(x,y)=x^3-4/x^2+5。设 f(x,y)=(1,1), y=x^3-4/x^2+5, y=2;从矩阵 var(-4)=matrix [1,2], var(x^2)=matrix[1,1], beta

以上是关于R语言绘制生存曲线95%区间的主要内容,如果未能解决你的问题,请参考以下文章

R语言绘制生存曲线估计|生存分析|如何R作生存曲线图

R语言生存分析之组间生存曲线的对比: Log-Rank检验绘制漂亮的生存曲线

R语言使用aov函数进行单因素方差分析(One-way ANOVA)R语言使用gplots包中的plotmeans函数可用于生成组均值及其置信区间f的图(95%置信区间的治疗方法的曲线图)

R语言cox回归模型案例(绘制列线图校正曲线):放疗是否会延长胰脏癌手术患者的生存时间

R语言Kaplan-Meier绘制生存分析Log-rank假设检验Cox回归曲线实战案例:恶性黑色素瘤的术后数据生存分析

R语言生存分析之COX比例风险模型构建及列线图(nomogram)校准曲线(calibration curve)绘制示例