我想绘制一个在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") )
这给了我这个:
这不是有史以来最漂亮的图表，但是我对ggplot2的了解还不够。但更重要的是，我有第二组数据来自对数正态分布，而不是指数分布，而我尝试为此生成数据的尝试完全失败了–将CDF纳入正态分布之类的方法这超出了我的R技能。
有人能够使用相同的数字和1的比例参数将我指向正确的方向吗？
最佳答案
对数正态模型在时间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的情况。
关于r - 生成/绘制对数正态生存函数，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/11103637/