このパッケージは、カーネル平滑化法を使用して、右打ち切りデータからハザード関数muhaz
を推定します。私の質問は、計算するハザード関数の信頼区間を取得する方法はありますか?muhaz
options(scipen=999)
library(muhaz)
data(ovarian, package="survival")
attach(ovarian)
fit1 <- muhaz(futime, fustat)
plot(fit1, lwd=3, ylim=c(0,0.002))
上記の例では、muhaz.object
fit
にいくつかのエントリfit1$msemin
,がありますがfit1$var.min
、fit1$haz.est
それらの長さは の半分ですfit1$haz.est
。
ハザード関数の信頼区間を抽出することが可能であれば、何か考えはありますか?
編集:@ user20650が提案したことに基づいて、次のことを試しました
options(scipen=999)
library(muhaz)
data(ovarian, package="survival")
fit1 <- muhaz(ovarian$futime, ovarian$fustat,min.time=0, max.time=744)
h.df<-data.frame(est=fit1$est.grid, h.orig=fit1$haz.est)
for (i in 1:10000){
d.s.onarian<-ovarian[sample(1:nrow(ovarian), nrow(ovarian), replace = T),]
d.s.muhaz<-muhaz(d.s.onarian$futime, d.s.onarian$fustat, min.time=0, max.time=744 )
h.df<-cbind(h.df, d.s.muhaz$haz.est)
}
h.df$upper.ci<-apply(h.df[,c(-1,-2)], 1, FUN=function(x) quantile(x, probs = 0.975))
h.df$lower.ci<-apply(h.df[,c(-1,-2)], 1, FUN=function(x) quantile(x, probs = 0.025))
plot(h.df$est, h.df$h.orig, type="l", ylim=c(0,0.003), lwd=3)
lines(h.df$est, h.df$upper.ci, lty=3, lwd=3)
lines(h.df$est, h.df$lower.ci, lty=3, lwd=3)
max.time の設定は機能しているようで、すべてのブートストラップ サンプルには同じ推定グリッド ポイントがあります。ただし、得られたCIはほとんど意味がありません。通常、間隔は t=0 では狭く、時間とともに広くなる (情報が少なくなり、不確実性が増す) と予想されますが、得られた間隔は、時間とともに多かれ少なかれ一定であるように見えます。