生成/绘制对数正态生存函数

Question

生成/绘制对数正态生存函数

7

我在SAS LIFEREG中有一个加速失效时间模型，我想绘制它。由于 SAS 在绘图方面极其糟糕，我想在 R 中重新生成曲线的数据并在那里绘制它们。SAS 输出了一个比例尺（在指数分布固定为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") )

这给我带来了以下结果:

enter image description here

虽然不是最美丽的图表，但我并不太了解如何使用ggplot2来美化它。更重要的是，我有第二组数据，其来源于对数正态分布，而不是指数分布，我的尝试生成该数据已经彻底失败了 - 对于正态分布的累积分布函数等的融合将其超出了我的R技能范围。

有没有人能够指导我朝正确的方向前进，使用相同的数字和1的比例参数？

- EpiGrad

当您使用ODS时，SAS通常会提供非常漂亮的曲线。如果不使用SAS Graph，是否有在SAS中绘制生存曲线的选项？也许有一个默认图形看起来不错。 - Michael R. Chernick

1

在我看来，这个问题涉及SO-CV重叠领域，但比起SO更适合CV。虽然它是一个编程问题，但回答需要一些统计专业知识，因此根据CV的faq，应该放在CV上。 - jthetzel

据我所知，LIFEREG可以生成风险图和一些诊断图，但不能生成生存函数。公正地说，大多数人通常希望LIFETEST生成生存函数，但在这种特殊情况下我不需要。 - Fomite

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jthetzel · Accepted Answer

对于一个对数正态模型来说，时间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软件包可以处理左、右或区间截尾的加速失效时间模型，如survival task view中所列，如果您喜欢R而不是SAS。