6

SAS LIFEREG に、プロットしたい加速故障時間モデルがあります。SAS はグラフ作成が非常に苦手なので、実際に R で曲線のデータを再生成し、そこにプロットしたいと思います。SAS は、尺度 (指数分布が 1 に固定されている場合)、切片、および曝露された集団または曝露されていない集団に含まれる回帰係数を出力します。

2 つの曲線があり、1 つは曝露された集団用で、もう 1 つは曝露されていない集団用です。モデルの 1 つは指数分布であり、次のようにデータとグラフを作成しました。

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 番目のデータ セットがあり、そのデータを生成する試みが完全に失敗したことです。それは私のRスキルを超えています。

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

4

1 に答える 1

8

対数正規モデルの時間 t における生存関数は、R で表すことができます1 - plnorm()。ここ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")

そして結果の数字：

指数モデル

対数正規モデル

ご了承ください

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

と同じです

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

これは、対数正規生存関数の正準方程式の Φ です: S(t) = 1 − Φ((ln(t) − µ) / σ)

ご存じのとおり、サバイバル タスク ビューにリストされているように、左打ち切り、右打ち打ち、または間隔打ち切りで加速故障時間モデルを処理できる R パッケージが多数あります。 SAS。

于 2012-06-19T16:25:17.253 に答える