R でデルタ法を実装して、ワイブル生存曲線の MTTF 分散を計算しようとしています。shape パラメーターはalphaで、scale パラメーターは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)を取得できます。
当てはめたワイブル モデルは 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
しかし、残念ながら、私の mtf 分散は、インターネットの論文から取ったどの例とも一致しません。何度も見直しました…
コード全体は次のとおりです。
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 の経験が豊富な皆さんと共有したいと思います。
よろしく、ヴィニシウス。