Vinicius_rr
Vinicius_rr

Reputation: 59

confidence interval for MTTF - Weibull survival curve in R

I am trying to implement the Delta Method in R to calculate the MTTF variance of a Weibull survival curve. The shape parameter is alpha and scale parameter is delta. Variance = var; covariance = cov.

The equation is:

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.    

Where:

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

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

So the equation becomes:

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

I can take var(alpha), var(delta) and cov(alpha,delta) from variance-covariance matrix.

The fitted weibull model is called 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)

And

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

And more

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

Where coef is a variable which contains parameters alpha and delta from survreg adjust.

So, calculating MTTF:

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

And calculating mttf variance:

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

But unfortunatelly my mttf variance does not match to any example I took from internet papers. I revised it too many times...

The whole code is:

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)

So, from paper which I took these data the results are:

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

Which is very far from the results of my code.

But alpha, delta and MTTF from paper has same values of my code:

alpha = 2.273151
delta = 1.417457
MTTF = 2.067864

Please, I would like to share this difficulty with you guys, who have much more experience in R than me.

Regards, Vinícius.

Upvotes: 4

Views: 894

Answers (1)

user243620
user243620

Reputation: 1

I suggest that beta needs to be greater than -1 but from my own calculations; beta =2.

Upvotes: 0

Related Questions