生存分析之R包survival的参数生存回归模型

生存分析之R包survival的参数生存回归模型 by 小白鱼的生统笔记

生存分析之R包survival的参数生存回归模型

续前文生存分析（Survival Analysis）。

前文已经分别简介了非参数（如Kaplan-Meier分析）和半参数（如Cox回归模型）生存分析方法，本篇继续以R包survival中的方法为例展示参数生存回归模型。

与非参数或半参数方法相比，参数生存回归模型其假设存在特定形式的生存时间分布，并根据已知的生存时间分布模型估计模型参数，最后以分布模型计算生存率。

参数模型可分为以下几种：（1）参数比例风险模型，采用Cox模型的形式，但在基线风险上采用参数形式；（2）附加危害模型，其中预测变量以累加而非替代的方式贡献风险；（3）与传统线性回归最相似的AFT模型。

由于事件时间是正值，并且通常偏离正态分布，因此对称正态分布成为紧密拟合数据的较差选择，需要其它一些假定的参数分布模型代替正态分布，如指数分布（exponential）、Weibull分布、广义伽玛分布（generalized gamma）、对数正态分布（log-normal）、log-logistic等。在实践中，通过比较适合各种不同分布的模型来选择要使用的参数分布，考虑的候选分布的选择取决于先验假设或对数据的已知了解程度。但需要注意的是，没有分布可以提供完美的拟合，并且对于同一数据而言可能有多种合适的参数分布可供选择。

相比非参数（如Kaplan-Meier分析）和半参数（如Cox回归模型）生存分析方法，参数生存回归模型使用并不广泛（主要归因于生存时间分布难以描述，若分布模型选择不当将导致较差的统计功效，相比之下Kaplan-Meier分析或Cox回归无需考虑这一点），因此本篇仅简单演示下其在R语言中的过程，大致对它作个了解。

示例数据

例如，这里模拟一个服从Weibull分布的数据集作为演示。

#模拟数据set.seed(123)N      <- 180                  # number of observationsP      <- 3                    # number of groupssex    <- factor(sample(c("f", "m"), N, replace=TRUE))  # stratification factorX      <- rnorm(N, 0, 1)       # continuous covariateIV     <- factor(rep(LETTERS[1:P], each=N/P))  # factor covariateIVeff  <- c(0, -1, 1.5)        # effects of factor levels (1 -> reference level)Xbeta  <- 0.7*X + IVeff[unclass(IV)] + rnorm(N, 0, 2)weibA  <- 1.5                  # Weibull shape parameterweibB  <- 100                  # Weibull scale parameterU      <- runif(N, 0, 1)       # uniformly distributed RVeventT <- ceiling((-log(U)*weibB*exp(-Xbeta))^(1/weibA))   # simulated event time
# censoring due to study end after 120 daysobsLen <- 120                  # length of observation timecensT  <- rep(obsLen, N)       # censoring time = end of studyobsT   <- pmin(eventT, censT)  # observed censored event timesstatus <- eventT <= censT      # has event occured?dfSurv <- data.frame(obsT, status, sex, X, IV)          # data frame
head(dfSurv)

该模拟数据集中，obsT代表了患者的生存时间（单位为月）、status为生存状态（TRUE代表去世，FALSE代表生存）、性别（sex）、X为一组连续变量，IV为一组分类变量。

参数生存回归模型示例

survival包中，可使用函数survreg()拟合参数生存回归模型。

拟合患者生存时间与自变量的关系，同时考虑一组连续变量和一组分类变量。并且由于上述模拟数据集中的生存时间为Weibull分布，因此参数模型选择使用“weibull”。

library(survival)
#参数回归模型，详情 ?survreg#dist 指定假定的 y 变量的分布，这里以 weibull 为例fitWeib <- survreg(Surv(obsT, status) ~ X + IV, dist = 'weibull', data = dfSurv)fitWeib
summary(fitWeib)

显示全模型的p=8e-13<<0.01，是非常显著的，说明在多个自变量中包含了对生存时间具有影响的因素。上方统计量则显示了每个自变量的效应。

根据特定观察值对新数据进行拟合。

假设这里给定两组个体，使用上述拟合的参数回归模型预测生存能力。

#预测生存能力#模拟两组新个体dfNew <- data.frame(sex = factor(c('m', 'm'), levels = levels(dfSurv$sex)),    X = c(0, 0), IV = factor(c('A', 'C'), levels = levels(dfSurv$IV)))
#估计生存率的百分位数percs <- (1:99)/100FWeib <- predict(fitWeib, newdata=dfNew, type = 'quantile', p = percs, se = TRUE)FWeib

结果中，给出了随时间的生存率分布，以及置信区间。

绘制两组新观测值的估计的生存曲线。

#生存曲线绘制matplot(cbind(FWeib$fit[1, ], FWeib$fit[1, ] - 2*FWeib$se.fit[1, ], FWeib$fit[1, ] + 2*FWeib$se.fit[1, ]),     1-percs, type = 'l', main = expression(paste('Weibull-Fit ', hat(S)(t), ' mit SE')),    xlab = 't', ylab = 'Survival', lty = c(1, 2, 2), lwd = 2, col = 'blue')matlines(cbind(FWeib$fit[2, ], FWeib$fit[2, ] - 2*FWeib$se.fit[2, ], FWeib$fit[2, ] + 2*FWeib$se.fit[2, ]),     1-percs, col = 'red', lwd = 2)legend(x = 'topright', lwd = 2, lty = c(1, 2, 1, 2), col = c('blue', 'blue', 'red', 'red'),    legend = c('sex=m, X=0, IV=A', '+- 2*SE', 'sex=m, X=0, IV=C', '+- 2*SE'))