作者热门文章
- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我在 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") )
最佳答案
对数正态模型在时间 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)
关于r - 生成/绘制对数正态生存函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/11103637/
我每天都进行回归分析。就我而言,这通常意味着估计连续和分类预测变量对各种结果的影响。生存分析可能是我执行的最常见的分析。此类分析通常以非常方便的方式出现在期刊中。下面是一个例子: 我想知道是否有人遇到
我是一名优秀的程序员,十分优秀!