対数正規生存関数の生成/プロット

Question

私はプロットしたいと思っているSAS LIFEREGの加速故障時間モデルを持っています。 SASはグラフ作成に深刻な悪影響を与えるので、実際にRの曲線のデータを再生成してそこにプロットしたいと思います。 SASは、曝露されたまたは未曝露の集団に存在するためのスケール（1に固定された指数分布の場合）、切片、および回帰係数を記載する。

露出されたものと未露光のものの2つの曲線があります。モデルの一つは、指数分布である、と私はそうのようなデータとグラフを作成しました：

私はこれを与え

intercept <- 5.00 
effect<- -0.500 
data<- data.frame(time=seq(0:180)-1) 
data$s_unexposed <- apply(data,1,function(row) exp(-(exp(-intercept))*row[1])) 
data$s_exposed <- apply(data,1,function(row) exp(-(exp(-(intercept+effect))*row[1]))) 

plot(data$time,data$s_unexposed, type="l", ylim=c(0,1) ,xaxt='n', 
    xlab="Days since Infection", ylab="Percent Surviving", lwd=2) 
axis(1, at=c(0, 20, 40, 60, 80, 100, 120, 140, 160, 180)) 
lines(data$time,data$s_exposed, col="red",lwd=2) 
legend("topright", c("ICU Patients", "Non-ICU Patients"), lwd=2, col=c("red","black")) 

：

ない、これまできれいなグラフが、私ggplot2の周りに私のやり方を本当に分かっていないでください。しかし、もっと重要なのは、指数関数ではなく、ログ正規分布から得られる第2のデータセットがあり、そのデータを生成しようとする試みがまったく失敗したことです。正規分布などでcdfを組み込むと失敗しましたそれは私のRスキルを超えています。

同じ数字を使用し、1のスケールパラメータを使用して、誰かが正しい方向に向けることができますか？

Answer 1

対数正規モデルの時刻tにおける生存関数は、1 - plnorm()でRで表すことができます。ここで、plnorm()は対数正規累積分布関数です。説明するために、我々は利便性のための機能にあなたのプロットを最初に出してあげる：

## Function to plot aft data 
plot.aft <- function(x, legend = c("ICU Patients", "Non-ICU Patients"), 
    xlab = "Days since Infection", ylab="Percent Surviving", lwd = 2, 
    col = c("red", "black"), at = c(0, 20, 40, 60, 80, 100, 120, 140, 160, 180), 
     ...) 
{ 
    plot(x[, 1], x[, 2], type = "l", ylim = c(0, 1), xaxt = "n", 
      xlab = xlab, ylab = ylab, col = col[2], lwd = 2, ...) 
    axis(1, at = at) 
    lines(x[, 1], x[, 3], col = col[1], lwd=2) 
    legend("topright", legend = legend, lwd = lwd, col = col) 
} 

次に、我々は、係数、変数、およびモデルを指定し、指数の生存確率を生成し、対数正規ますモデル：

## Specify coefficients, variables, and linear models 
beta0 <- 5.00 
beta1 <- -0.500 
icu <- c(0, 1) 
t <- seq(0, 180) 
linmod <- beta0 + (beta1 * icu) 
names(linmod) <- c("unexposed", "exposed") 

## Generate s(t) from exponential AFT model 
s0.exp <- dexp(exp(-linmod["unexposed"]) * t) 
s1.exp <- dexp(exp(-linmod["exposed"]) * t) 

## Generate s(t) from lognormal AFT model 
s0.lnorm <- 1 - plnorm(t, meanlog = linmod["unexposed"]) 
s1.lnorm <- 1 - plnorm(t, meanlog = linmod["exposed"]) 

最後に、我々は生存確率をプロットすることができます

## Plot survival 
plot.aft(data.frame(t, s0.exp, s1.exp), main = "Exponential model") 
plot.aft(data.frame(t, s0.lnorm, s1.lnorm), main = "Log-normal model") 

、そして得られた数値：

S：

plnorm(t, meanlog = linmod["exposed"]) 

は対数正規生存関数の正規方程式におけるΦある

pnorm((log(t) - linmod["exposed"])/1) 

と同じであること0

注（t）= 1-Φ（（ln（t）-μ）/σ）

survival task viewにリストされているように、左、右、または区間の打ち切りで加速された失敗時間モデルを処理できる多くのRパッケージがあります。 SAS経由。