差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン
前のリビジョン
r:survival [2018/12/09 17:53] watalur:survival [2018/12/11 19:42] (現在) watalu
行 1: 行 1:
 ==== Rによる生存時間解析 ==== ==== Rによる生存時間解析 ====
 +
 +[[r:survival:code]]の方に記したコードを、先にRに貼り付けて読み込ませておくこと。いくつかのグラフを描く関数を用意しておいたものである。参考までに、用いているRの関数を個別に詳細に説明したページも用意している。
 +
 +  * [[r:survival:data.frame]]は寿命データの記し方を説明している
 +  * [[r:mass:fitdistr]]は完全データからの確率分布のパラメータの最尤法に基づいて推定する
 +  * [[r:survival:survfit]]は打ち切りを含む寿命データから、生存関数をノンパラメトリックに推定する
 +
 +ただし、実験の内容を実施するためには、上の3つのページは読まずとも問題なく、下の説明に従って進めれば十分である。
  
 === データの書き方 === === データの書き方 ===
行 89: 行 97:
 <code> <code>
 X$status X$status
 +</code>
 +
 +=== 信頼性データ図 ===
 +
 +信頼性データを図示する関数draw.reliability.data.diagramを用意した。
 +上に与えたデータの図を描くには次の1行を実行する。
 +
 +<code>
 +plot.reliability.data.diagram(X)
 </code> </code>
  
行 99: 行 116:
 </code> </code>
  
-これを実行しておかないと、すぐ次のコマンドの実行時にい文字でエラーが表示される。+これを実行しておかないと、すぐ次のコマンドの実行時にい文字でエラーが表示されるかもしれない 
 + 
 +=== 生存曲線と故障率曲線のノンパラメトリック推定 ===
  
-=== 生存曲線と故障率曲線の推定 ===+特定確率分布を仮定せずに母数や分布を推定する方法を、ノンパラメトリックな推定法と呼ぶ。
  
-カプランさんとマイヤーさんが提案した、特定の確率分布を仮定せずに生存関数を推定する方法カプラン・マイヤー推定量、もしくはカプラン・マイヤーいう。生存時間データから寿命分布を推定する際に、まずは何も仮定せずに生存関数を推定して眺めてみることが、第一歩となる。+カプランさんとマイヤーさんが提案した、生存関数をノンパラメトリックに推定する方法カプラン・マイヤー法である。カプラン・マイヤー推定量も呼ばれる。生存時間データから寿命分布を推定する際に、まずは何も仮定せずに生存関数を推定して眺めてみることが、第一歩となる。
  
 <code> <code>
行 166: 行 185:
 が故障率曲線のノンパラメトリックな推定値となる。 が故障率曲線のノンパラメトリックな推定値となる。
  
 +=== ノンパラメトリックな推定量の長所と短所 ===
 +
 +ノンパラメトリックな推定量は、確率分布を仮定しないので、当初から色眼鏡で見ることがなく、データをそのまま吟味できるという長所を持つ。実際に観測された時点から先の予測ができない、という短所も合わせ持つ。また、精度よく確率分布を推定するには、とても大きいサンプルを必要とすることも、短所である。
 +
 +信頼性は時間に関する品質である。販売あるいは投入したすべてのシステムが故障した後に、寿命分布を正しく推定できても、意味がない。可能ならば市場に投入したり現場で運用を開始する前や大量生産に踏み切る前の、設計開発段階や試作段階で、寿命分布を推定しておきたい。そうして、品質保証制度に基づく補償請求に応じる費用や、運用期間中の総保守費用を予測しておきたい。更に可能ならば、それらの費用を削減するように、設計や開発を行いたいのである。
 +
 +そこでシステムや製品の信頼性や保全性を設計・計画する際には、具体的に寿命が従う確率分布(寿命分布)をパラメトリックに推定して、それに基づいて対応することが通例となっている。
 +
 +=== パラメトリックな寿命分布の例 ===
 +== ワイブル分布 ==
 +
 +
 +== 対数正規分布 ==
 +
 +
 +== ガンマ分布 ==
 +
 +
 +== グンベル分布 ==
 +
 +
 +=== 故障時間データからのパラメトリックな寿命分布の推定法 ===
 +== すべてが故障しているデータ ==
 +
 +すべてが故障しているデータの場合は、MASSパッケージの中の関数fitdistrを用いて、確率分布の未知母数の最尤推定値と標準誤差を得ることができる。まず、MASSパッケージを読み込む。
 +
 +<code>
 +library(MASS)
 +</code>
 +
 +次に、関数fitdistrを用いて、ワイブル分布を仮定し、その2つの母数を推定してみる。
 +<code>
 +fitdistr(X$time, densfun="weibull")
 +</code>
 +その結果、次の数値を得る。
 +<code>
 +     shape       scale  
 +  1.6822280   6.8219220 
 + (0.4279334) (1.3515357)
 +</code>
 +ワイブル分布には形状母数shapeと尺度母数scaleがある。それぞれの点推定値は1.6822280と6.8219220であり、標準誤差は括弧の中に表示される。
 +
 +これらは、次のようにしても表示される。
 +<code>
 +X.fit.weibull <- fitdistr(X$time, densfun="weibull")
 +names(X.fit.weibull)
 +X.fit.weibull$estimate
 +X.fit.weibull$sd
 +X.fit.weibull$loglik
 +</code>
 +最後の数値は、最尤推定法で最大化した対数尤度関数の最大値である。この値を参考に、寿命分布を選ぶこともある。
 +
 +対数正規分布も推定しておく。
 +<code>
 +X.fit.lognormal <- fitdistr(X$time, densfun="lognormal")
 +names(X.fit.lognormal)
 +X.fit.lognormal$estimate
 +X.fit.lognormal$sd
 +X.fit.lognormal$loglik
 +</code>
 +
 +同様に、ガンマ分布。
 +<code>
 +X.fit.gamma <- fitdistr(X$time, densfun="gamma")
 +names(X.fit.gamma)
 +X.fit.gamma$estimate
 +X.fit.gamma$sd
 +X.fit.gamma$loglik
 +</code>
 +
 +=== 関数 ===
 +
 +完全データについて、これら三つの確率分布を推定し、ヒストグラムと比較するための関数compare.histgram.and.densitiesを用意した。
 +
 +オプションのhistをhist=FALSEとすると、ヒストグラムは描かず、密度関数のみを描くようになる。
 +
 +またカプラン・マイヤー法で推定した生存関数と共に描くための関数compare.KM.and.distributionsも用意した。
 +これで、どの分布に最も近いかを検討できる。
 +
 +<code>
 +X.fits <- compare.histgram.and.densities(X)
 +compare.KM.and.distributions(X.fits,data=X)
 +</code>
 +
 +== 一部が打ち切られているデータ ==
 +
 +一部の寿命の観測が、故障まで継続されず、途中で打ち切られているデータを、打ち切りデータという。打ち切りデータの場合は、MASSパッケージの中の関数fitdistrは用いることができない。代わりに、survivalパッケージの中の関数survregを用いる。この関数も確率分布の未知母数の最尤推定値と標準誤差を得ることができる。まずsurvivalパッケージを読み込む。
 +
 +<code>
 +library(survival)
 +</code>
 +
 +関数survregは、ワイブル分布weibull、対数正規分布lognormal、ロジスティック分布logistic、対数ロジスティック分布loglogisticの4種類の確率分布の最尤推定法を提供する。ガンマ分布には対応していない。
 +そのため以下では、ワイブル分布と対数正規分布のみを対象として、話を進める。
 +
 +<code>
 +X.survreg.weibull <- survreg(Surv(time, status)~1, data=X, dist="weibull")
 +names(X.survreg.weibull)
 +c(1/X.survreg.weibull$scale,exp(X.survreg.weibull$coefficients))
 +summary(X.survreg.weibull)
 +X.fit.weibull$loglik
 +</code>
 +
 +<code>
 +X.survreg.lognormal <- survreg(Surv(time, status)~1, data=X, dist="lognormal")
 +names(X.survreg.lognormal)
 +c(X.survreg.lognormal$coefficients, X.survreg.lognormal$scale)
 +summary(X.survreg.lognormal)
 +X.fit.lognormal$loglik
 +</code>
 +
 +これらをカプラン・マイヤー曲線に重ねて描く関数compare.KM.and.distributions.2 も用意した。
 +
 +<code>
 +compare.KM.and.distributions.2(X)
 +</code>