2012-02-05 13 views
23

私はWeibullモデルを生存データに合わせてプロットしようとしています。データには、2006年から2010年にかけての共変量、コホートが1つしかありません。したがって、2010年のコホートの生存曲線をプロットするために、次の2行のコードに何を追加するかについてのアイデアはありますか?survreg(Rのパッケージサバイバル)によって生成された生存曲線をプロットするにはどうすればよいですか?

library(survival) 
s <- Surv(subSetCdm$dur,subSetCdm$event) 
sWei <- survreg(s ~ cohort,dist='weibull',data=subSetCdm) 

Cox PHモデルで同じことを達成するのは簡単で、次のような行があります。問題は、survfit()が型survivのオブジェクトを受け入れないことです。

sCox <- coxph(s ~ cohort,data=subSetCdm) 
cohort <- factor(c(2010),levels=2006:2010) 
sfCox <- survfit(sCox,newdata=data.frame(cohort)) 
plot(sfCox,col='green') 

データ肺(生存パッケージから)を使用して、私が達成しようとしているのはここです。

#create a Surv object 
s <- with(lung,Surv(time,status)) 

#plot kaplan-meier estimate, per sex 
fKM <- survfit(s ~ sex,data=lung) 
plot(fKM) 

#plot Cox PH survival curves, per sex 
sCox <- coxph(s ~ as.factor(sex),data=lung) 
lines(survfit(sCox,newdata=data.frame(sex=1)),col='green') 
lines(survfit(sCox,newdata=data.frame(sex=2)),col='green') 

#plot weibull survival curves, per sex, DOES NOT RUN 
sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung) 
lines(survfit(sWei,newdata=data.frame(sex=1)),col='red') 
lines(survfit(sWei,newdata=data.frame(sex=2)),col='red') 
+2

完全な例を投稿した場合、私はそれを理解しようとします。 subSetCdmオブジェクトが必要です。 try dput(subSetCdm) –

+3

'?predict.survreg'には例があります。 –

答えて

18

・ホープ、このことができますし、私はいくつかの誤解を招くようなミスを犯していない:上からコピー

:ワイブル、使用する

#create a Surv object 
    s <- with(lung,Surv(time,status)) 

    #plot kaplan-meier estimate, per sex 
    fKM <- survfit(s ~ sex,data=lung) 
    plot(fKM) 

    #plot Cox PH survival curves, per sex 
    sCox <- coxph(s ~ as.factor(sex),data=lung) 
    lines(survfit(sCox,newdata=data.frame(sex=1)),col='green') 
    lines(survfit(sCox,newdata=data.frame(sex=2)),col='green') 

はヴィンセントからのコメント再予測

#plot weibull survival curves, per sex, 
    sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung) 

    lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red") 
    lines(predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red") 

plot output

プロットと予測の分位数の順序が逆転していました。これを行うより良い方法がありますが、ここで動作します。がんばろう!

+0

ティム、クイック質問です。 genetで上記の部分集合ではない部分を再作成したい場合、例えばwith_sWei < - survreg(s〜1、dist = 'weibull'、data = lung)を開始します。予測関数のnewdata部分の仕様をどのように変更しますか?私はあなたが上で指定した方法をgrokしようとしています... – Chris

+0

こんにちはクリス、私は困っています、申し訳ありませんが、おそらく他の回答者の1つは知っています。そうでなければ、新しい質問かもしれません。 –

13

代替オプションは、パッケージflexsurvを使用することです。これは、survivalパッケージにいくつかの追加機能を提供します - パラメトリック回帰関数flexsurvreg()には、あなたが求めるものを行う素敵なプロット方法が含まれています。

上記のように肺を使用する。

#create a Surv object 
s <- with(lung,Surv(time,status)) 

require(flexsurv) 
sWei <- flexsurvreg(s ~ as.factor(sex),dist='weibull',data=lung) 
sLno <- flexsurvreg(s ~ as.factor(sex),dist='lnorm',data=lung) 

plot(sWei) 
lines(sLno, col="blue") 

output from plot.flexsurvreg

あなたはtype引数を使用して累積ハザードやハザードスケールでプロットし、ci引数で信頼区間を追加することができます。

6

これは以下のコードを使用Tim Riffe's answerを明らかにするだけノートれる:予測()メソッドのための変位値を与えているので

lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red") 
lines(predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red") 

2鏡画像シーケンス、seq(.01,.99,by=.01)seq(.99,.01,by=-.01)理由でありますイベント分布f(t)、すなわちf(t)の逆CDFの値であるが、生存曲線は1(CDFのf)対tをプロットしている。言い換えれば、p対predict(p)をプロットするとCDFが得られ、1-pとpredict(p)をプロットすると生存曲線(1-CDF)が得られます。次のコードは、より透明でp値の任意のベクトルに一般化されます:

pct <- seq(.01,.99,by=.01) 
lines(predict(sWei, newdata=list(sex=1),type="quantile",p=pct),1-pct,col="red") 
lines(predict(sWei, newdata=list(sex=2),type="quantile",p=pct),1-pct,col="red") 
関連する問題