MTTF 的置信区间 - R 中的 Weibull 生存曲线
Posted
技术标签:
【中文标题】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
这与我的代码结果相差甚远。
但是纸上的 alpha、delta 和 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以上是关于MTTF 的置信区间 - R 中的 Weibull 生存曲线的主要内容,如果未能解决你的问题,请参考以下文章
R语言使用aov函数进行单因素方差分析(One-way ANOVA)R语言使用gplots包中的plotmeans函数可用于生成组均值及其置信区间f的图(95%置信区间的治疗方法的曲线图)