差分
このページの2つのバージョン間の差分を表示します。
両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン | ||
r:survival [2018/12/09 21:53] – watalu | r:survival [2018/12/11 19:42] (現在) – watalu | ||
---|---|---|---|
行 1: | 行 1: | ||
==== Rによる生存時間解析 ==== | ==== Rによる生存時間解析 ==== | ||
+ | |||
+ | [[r: | ||
+ | |||
+ | * [[r: | ||
+ | * [[r: | ||
+ | * [[r: | ||
+ | |||
+ | ただし、実験の内容を実施するためには、上の3つのページは読まずとも問題なく、下の説明に従って進めれば十分である。 | ||
=== データの書き方 === | === データの書き方 === | ||
行 94: | 行 102: | ||
信頼性データを図示する関数draw.reliability.data.diagramを用意した。 | 信頼性データを図示する関数draw.reliability.data.diagramを用意した。 | ||
- | |||
- | < | ||
- | draw.reliability.data.diagram <- function(x, | ||
- | if(is.numeric(x)==TRUE) { | ||
- | n <- length(x) | ||
- | if(sort==TRUE) { | ||
- | x.time <- x[sort.list(x, | ||
- | x.status <- rep(1,n) | ||
- | id <- c(1:n) | ||
- | } else { | ||
- | x.time <- x | ||
- | x.status <- rep(1,n) | ||
- | id <- c(1:n) | ||
- | } | ||
- | } else if (is.data.frame(x)==TRUE) { | ||
- | n <- dim(x)[1] | ||
- | if(sort==TRUE) { | ||
- | x.time <- x[sort.list(x[, | ||
- | x.status <- x[sort.list(x[, | ||
- | id <- row.names(x)[sort.list(x[, | ||
- | } else { | ||
- | x.time <- x[,1] | ||
- | x.status <- x[,2] | ||
- | id <- row.names(x) | ||
- | } | ||
- | } else if (is.matrix(x)==TRUE) { | ||
- | n <- dim(x)[1] | ||
- | if(sort==TRUE) { | ||
- | x.time <- x[sort.list(x[, | ||
- | x.status <- x[sort.list(x[, | ||
- | id <- rawnames(x)[sort.list(x[, | ||
- | } else { | ||
- | x.time <- x[,1] | ||
- | x.status <- x[,2] | ||
- | id <- rawnames(x) | ||
- | } | ||
- | } | ||
- | plot(x.time, | ||
- | | ||
- | | ||
- | | ||
- | | ||
- | points(x.time, | ||
- | for( i in c(1:n) ) { | ||
- | lines(c(0, | ||
- | } | ||
- | } | ||
- | </ | ||
- | |||
上に与えたデータの図を描くには次の1行を実行する。 | 上に与えたデータの図を描くには次の1行を実行する。 | ||
< | < | ||
- | draw.reliability.data.diagram(X) | + | plot.reliability.data.diagram(X) |
</ | </ | ||
行 157: | 行 116: | ||
</ | </ | ||
- | これを実行しておかないと、すぐ次のコマンドの実行時に赤い文字でエラーが表示される。 | + | これを実行しておかないと、すぐ次のコマンドの実行時に青い文字でエラーが表示されるかもしれない。 |
=== 生存曲線と故障率曲線のノンパラメトリック推定 === | === 生存曲線と故障率曲線のノンパラメトリック推定 === | ||
行 298: | 行 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(" | + | |
- | | + | |
- | } | + | |
- | </ | + | |
これで、どの分布に最も近いかを検討できる。 | これで、どの分布に最も近いかを検討できる。 | ||
行 398: | 行 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) | ||
</ | </ |