Reputation: 59
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
Reputation: 1
I suggest that beta needs to be greater than -1 but from my own calculations; beta =2.
Upvotes: 0