16

私は次のsurvregモデルを持っています:

Call:
survreg(formula = Surv(time = (ev.time), event = ev) ~ age, 
    data = my.data, dist = "weib")
             Value Std. Error    z        p
(Intercept) 4.0961     0.5566 7.36 1.86e-13
age         0.0388     0.0133 2.91 3.60e-03
Log(scale)  0.1421     0.1208 1.18 2.39e-01
Scale= 1.15 

Weibull distribution

上記の推定値に基づいて、ハザード関数と生存関数をプロットしたいと思います。or
は使用したくありません(ここで提示されているように、パラメトリックサバイバルまたはここでSO question . predict()pweibull()

機能を利用したいcurve()。これを達成する方法はありますか？survreg のワイブル関数は、通常とは異なるスケールと形状の定義を使用しているようです (たとえば rweibull とは異なります)。

更新：既製の関数を使用せずInterceptにage (+ other potential covariates)、推定値の関数としてハザード/生存率を表現するために本当に必要なものを推測します。Scale*weilbull

4

2 に答える 2

11

あなたが提供した最初のリンクには、実際にこれがどのように機能するかの理論に関する明確な説明と、素敵な例があります. （これをありがとう、それは私が自分の仕事で使用する素晴らしいリソースです。）

関数を使用するには、curveいくつかの関数を引数として渡す必要があります。*weibull関数のファミリが とは異なるワイブルのパラメーター化を使用することは事実ですがsurvreg、最初のリンクで説明したように、簡単に変換できます。また、 のドキュメントからsurvreg:

ワイブル分布をパラメータ化する方法は複数あります。survreg 関数は、rweibull 関数とは異なるパラメーター化である一般的な位置スケール ファミリにそれを組み込み、しばしば混乱を招きます。

  survreg's scale  =    1/(rweibull shape)
  survreg's intercept = log(rweibull scale)

その単純な変換の実装を次に示します。

# The parameters
intercept<-4.0961
scale<-1.15

par(mfrow=c(1,2),mar=c(5.1,5.1,4.1,2.1)) # Make room for the hat.
# S(t), the survival function
curve(pweibull(x, scale=exp(intercept), shape=1/scale, lower.tail=FALSE), 
      from=0, to=100, col='red', lwd=2, ylab=expression(hat(S)(t)), xlab='t',bty='n',ylim=c(0,1))
# h(t), the hazard function
curve(dweibull(x, scale=exp(intercept), shape=1/scale)
      /pweibull(x, scale=exp(intercept), shape=1/scale, lower.tail=FALSE), 
      from=0, to=100, col='blue', lwd=2, ylab=expression(hat(h)(t)), xlab='t',bty='n')
par(mfrow=c(1,1),mar=c(5.1,4.1,4.1,2.1))

生存関数とハザード関数

回答でpweibull関数を使用したくないとおっしゃったことは理解していますが、別のパラメーター化を使用しているため、使用したくないと推測しています。それ以外の場合は、そのパラメーター化pweibullを使用する独自のバージョンを単純に作成できます。survreg

my.weibull.surv<-function(x,intercept,scale) pweibull(x,scale=exp(intercept),shape=1/scale,lower.tail=FALSE)
my.weibull.haz<-function(x,intercept,scale) dweibull(x, scale=exp(intercept), shape=1/scale) / pweibull(x,scale=exp(intercept),shape=1/scale,lower.tail=FALSE)

curve(my.weibull.surv(x,intercept,scale),1,100,lwd=2,col='red',ylim=c(0,1),bty='n')
curve(my.weibull.haz(x,intercept,scale),1,100,lwd=2,col='blue',bty='n')

コメントで述べたように、なぜこれを行うのかわかりませんが (これが宿題でない限り)、必要pweibullに応じてハンドコーディングできdweibullます。

my.dweibull <- function(x,shape,scale) (shape/scale) * (x/scale)^(shape-1) * exp(- (x/scale)^shape)
my.pweibull <- function(x,shape,scale) exp(- (x/scale)^shape)

これらの定義は からそのまま出てき?dweibullます。直接ラップする代わりにpweibull、遅くてテストされていない関数をラップするだけです。dweibull

于 2013-04-26T15:07:54.837 に答える
11

パラメータは次のとおりです。

scale=exp(Intercept+beta*x)あなたの例では、 age=40 としましょう

scale=283.7

形状パラメータは、モデルが出力するスケールの逆数です

shape=1/1.15

次に、危険性は次のとおりです。

curve((shape/scale)*(x/scale)^(shape-1), from=0,to=12,ylab=expression(hat(h)(t)), col="darkblue",xlab="t", lwd=5)

累積ハザード関数は次のとおりです。

curve((x/scale)^(shape), from=0,to=12,ylab=expression(hat(F)(t)), col="darkgreen",xlab="t", lwd=5)

生存関数は 1-累積ハザード関数なので、次のようになります。

curve(1-((x/scale)^(shape)), from=0,to=12,ylab=expression(hat(S)(t)), col="darkred",xlab="t", lwd=5, ylim=c(0,1))

ehaまた、パッケージと機能hweibullをチェックしてください。Hweibull

ワイブル関数

于 2013-04-26T16:51:04.197 に答える