差分
このページの2つのバージョン間の差分を表示します。
両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン | ||
r:survival [2018/12/09 21:58] – watalu | r:survival [2018/12/11 19:42] (現在) – watalu | ||
---|---|---|---|
行 1: | 行 1: | ||
==== Rによる生存時間解析 ==== | ==== Rによる生存時間解析 ==== | ||
- | [[r: | + | [[r: |
+ | |||
+ | * [[r: | ||
+ | * [[r: | ||
+ | * [[r: | ||
+ | |||
+ | ただし、実験の内容を実施するためには、上の3つのページは読まずとも問題なく、下の説明に従って進めれば十分である。 | ||
=== データの書き方 === | === データの書き方 === | ||
行 99: | 行 105: | ||
< | < | ||
- | draw.reliability.data.diagram(X) | + | plot.reliability.data.diagram(X) |
</ | </ | ||
行 110: | 行 116: | ||
</ | </ | ||
- | これを実行しておかないと、すぐ次のコマンドの実行時に赤い文字でエラーが表示される。 | + | これを実行しておかないと、すぐ次のコマンドの実行時に青い文字でエラーが表示されるかもしれない。 |
=== 生存曲線と故障率曲線のノンパラメトリック推定 === | === 生存曲線と故障率曲線のノンパラメトリック推定 === | ||
行 251: | 行 257: | ||
=== 関数 === | === 関数 === | ||
- | 完全データについて、これら三つの確率分布を推定し、ヒストグラムと比較するための関数を用意した。 | + | 完全データについて、これら三つの確率分布を推定し、ヒストグラムと比較するための関数compare.histgram.and.densitiesを用意した。 |
- | < | + | |
- | compare.histgram.and.densities | + | |
- | require(MASS) | + | |
- | est.weibull <- fitdistr(X$time, | + | |
- | est.lnorm <- fitdistr(X$time, | + | |
- | est.gamma <- fitdistr(X$time, | + | |
- | if(hist==TRUE) { | + | |
- | hist(X$time, | + | |
- | | + | |
- | plot(function(x){return(dlnorm(x, | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | } else { | + | |
- | plot(function(x){return(dlnorm(x, | + | |
- | | + | |
- | | + | |
- | | + | |
- | } | + | |
- | plot(function(x){return(dweibull(x, | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | plot(function(x){return(dgamma(x, | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | legend(" | + | |
- | | + | |
- | | + | |
- | | + | |
- | return(list(weibull=est.weibull, | + | |
- | } | + | |
- | </ | + | |
オプションのhistをhist=FALSEとすると、ヒストグラムは描かず、密度関数のみを描くようになる。 | オプションのhistをhist=FALSEとすると、ヒストグラムは描かず、密度関数のみを描くようになる。 | ||
- | またカプラン・マイヤー法で推定した生存関数と共に描くための関数も用意した。 | + | またカプラン・マイヤー法で推定した生存関数と共に描くための関数compare.KM.and.distributionsも用意した。 |
- | < | + | |
- | compare.KM.and.distributions | + | |
- | require(survival) | + | |
- | km <- survfit(Surv(time, | + | |
- | plot(km, | + | |
- | | + | |
- | | + | |
- | plot(function(x) pweibull(x, | + | |
- | lower.tail=FALSE), | + | |
- | | + | |
- | plot(function(x) plnorm(x, | + | |
- | lower.tail=FALSE), | + | |
- | | + | |
- | plot(function(x) pgamma(x, | + | |
- | lower.tail=FALSE), | + | |
- | | + | |
- | legend(" | + | |
- | | + | |
- | paste(" | + | |
- | paste(" | + | |
- | paste(" | + | |
- | | + | |
- | } | + | |
- | </ | + | |
これで、どの分布に最も近いかを検討できる。 | これで、どの分布に最も近いかを検討できる。 | ||
行 351: | 行 296: | ||
</ | </ | ||
- | これらをカプラン・マイヤー曲線に重ねて描く関数も用意した。 | + | これらをカプラン・マイヤー曲線に重ねて描く関数compare.KM.and.distributions.2 |
- | + | ||
- | < | + | |
- | compare.KM.and.distributions.2 | + | |
- | require(survival) | + | |
- | km <- survfit(Surv(time, | + | |
- | fits <- compare.histgram.and.densities(data, | + | |
- | survreg.weibull <- survreg(Surv(time, | + | |
- | fits$weibull$estimate[1] <- 1/ | + | |
- | fits$weibull$estimate[2] <- exp(survreg.weibull$coefficients[1]) | + | |
- | survreg.lognormal <- survreg(Surv(time, | + | |
- | fits$lnorm$estimate[1] <- survreg.lognormal$coefficients[1] | + | |
- | fits$lnorm$estimate[2] <- survreg.lognormal$scale | + | |
- | plot(km, | + | |
- | | + | |
- | | + | |
- | plot(function(x) pweibull(x, | + | |
- | lower.tail=FALSE), | + | |
- | | + | |
- | plot(function(x) plnorm(x, | + | |
- | lower.tail=FALSE), | + | |
- | | + | |
- | legend(" | + | |
- | | + | |
- | paste(" | + | |
- | paste(" | + | |
- | | + | |
- | } | + | |
- | </ | + | |
< | < | ||
compare.KM.and.distributions.2(X) | compare.KM.and.distributions.2(X) | ||
</ | </ |