variance - confidence interval for MTTF - Weibull survival curve in R -

i trying implement delta method in r calculate mttf variance of weibull survival curve. shape parameter alpha , scale parameter 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 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) , cov(alpha,delta) variance-covariance matrix.

the fitted weibull model 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 variable contains parameters alpha , delta 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 mttf variance not match example took internet papers. revised 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 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 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, paper took these data results are:

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

which far results of code.

but alpha, delta , mttf paper has same values of code:

alpha = 2.273151 delta = 1.417457 mttf = 2.067864 

please, share difficulty guys, have more experience in r me.

regards, vinícius.