R生存分析AFT
γ = 1/scale =1/0.902
α = exp(−(Intercept)γ)=exp(-(7.111)*γ)
> library(survival)
> myfit=survreg(Surv(futime, fustat)~1 , ovarian, dist="weibull",scale=0)
> summary(myfit)
Call:
survreg(formula = Surv(futime, fustat) ~ 1, data = ovarian, dist = "weibull",
scale = 0)
Value Std. Error z p
(Intercept) 7.111 0.293 24.292 2.36e-130
Log(scale) -0.103 0.254 -0.405 6.86e-01
Scale= 0.902
Weibull distribution
Loglik(model)= -98 Loglik(intercept only)= -98
Number of Newton-Raphson Iterations: 5
n= 26
画生存函数图
d<- seq(0, 2000, length.out=10000)
h<-1-pweibull(d,shape=1/0.902,scale=exp(7.111))
df<-data.frame(t=d,s=h)
library(ggplot2)
ggplot(df,aes(x=t,y=s))+
geom_line(colour="green")+
ggtitle("s(t) \n 生存函数")
1. Surv
Description
创建一个生存对象，通常用作模型公式中的响应变量。 参数匹配是此功能的特殊功能，请参阅下面的详细信息。
Surv(time, time2, event,
type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'),
origin=0)
is.Surv(x)
Arguments
time
对于右删失数据，这是一个跟踪时间。对于区间数据，第一个参数是区间的开始时间。
event
状态指示，通常，0=活着，1=死亡。其他选择是TRUE/FALSE (TRUE = 死亡) or 1/2 (2=死亡)。对于区间删失数据，状态指示，0=右删失, 1=事件时间, 2=左删失, 3=区间删失.
右删失(Right Censoring)：只知道实际寿命大于某数；
左删失(Left Censoring)：只知道实际寿命小于某数；
区间删失(Interval Censoring)：只知道实际寿命在一个时间区间内。
time2
区间删失区间的结束时间或仅对过程数据进行计数。
type
指定删失类型。 "right", "left", "counting", "interval", "interval2" or "mstate".
如果event变量是一个因子，假定type="mstate"。如果没有指定参数time2，type="right"；如果指定参数time2，type="counting"
Surv使用示例
> str(lung) 'data.frame': 228 obs. of 10 variables: $ inst : num 3 3 3 5 1 12 7 11 1 7 ... $ time : num 306 455 1010 210 883 ... $ status : num 2 2 1 2 2 1 2 2 2 2 ... $ age : num 74 68 56 57 60 74 68 71 53 61 ... $ sex : num 1 1 1 1 1 1 2 2 1 1 ... $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ... $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ... $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ... $ meal.cal : num 1175 1225 NA 1150 NA ... $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ... > with(lung, Surv(time, status)) [1] 306 455 1010+ 210 883 1022+ 310 361 218 [10] 166 170 654 728 71 567 144 613 707 [19] 61 88 301 81 624 371 394 520 574 [28] 118 390 12 473 26 533 107 53 122 [37] 814 965+ 93 731 460 153 433 145 583 [46] 95 303 519 643 765 735 189 53 246 [55] 689 65 5 132 687 345 444 223 175 [64] 60 163 65 208 821+ 428 230 840+ 305 [73] 11 132 226 426 705 363 11 176 791 [82] 95 196+ 167 806+ 284 641 147 740+ 163 [91] 655 239 88 245 588+ 30 179 310 477 [100] 166 559+ 450 364 107 177 156 529+ 11 [109] 429 351 15 181 283 201 524 13 212 [118] 524 288 363 442 199 550 54 558 207 [127] 92 60 551+ 543+ 293 202 353 511+ 267 [136] 511+ 371 387 457 337 201 404+ 222 62 [145] 458+ 356+ 353 163 31 340 229 444+ 315+ [154] 182 156 329 364+ 291 179 376+ 384+ 268 [163] 292+ 142 413+ 266+ 194 320 181 285 301+ [172] 348 197 382+ 303+ 296+ 180 186 145 269+ [181] 300+ 284+ 350 272+ 292+ 332+ 285 259+ 110 [190] 286 270 81 131 225+ 269 225+ 243+ 279+ [199] 276+ 135 79 59 240+ 202+ 235+ 105 224+ [208] 239 237+ 173+ 252+ 221+ 185+ 92+ 13 222+ [217] 192+ 183 211+ 175+ 197+ 203+ 116 188+ 191+ [226] 105+ 174+ 177+ > str(heart) 'data.frame': 172 obs. of 8 variables: $ start : num 0 0 0 1 0 36 0 0 0 51 ... $ stop : num 50 6 1 16 36 39 18 3 51 675 ... $ event : num 1 1 0 1 0 1 1 1 0 1 ... $ age : num -17.16 3.84 6.3 6.3 -7.74 ... $ year : num 0.123 0.255 0.266 0.266 0.49 ... $ surgery : num 0 0 0 0 0 0 0 0 0 0 ... $ transplant: Factor w/ 2 levels "0","1": 1 1 1 2 1 2 1 1 1 2 ... $ id : num 1 2 3 3 4 4 5 6 7 7 ... > Surv(heart$start, heart$stop, heart$event) [1] ( 0.0, 50.0] ( 0.0, 6.0] ( 0.0, 1.0+] [4] ( 1.0, 16.0] ( 0.0, 36.0+] ( 36.0, 39.0] [7] ( 0.0, 18.0] ( 0.0, 3.0] ( 0.0, 51.0+] [10] ( 51.0, 675.0] ( 0.0, 40.0] ( 0.0, 85.0] [13] ( 0.0, 12.0+] ( 12.0, 58.0] ( 0.0, 26.0+] [16] ( 26.0, 153.0] ( 0.0, 8.0] ( 0.0, 17.0+] [19] ( 17.0, 81.0] ( 0.0, 37.0+] ( 37.0,1387.0] [22] ( 0.0, 1.0] ( 0.0, 28.0+] ( 28.0, 308.0] [25] ( 0.0, 36.0] ( 0.0, 20.0+] ( 20.0, 43.0] [28] ( 0.0, 37.0] ( 0.0, 18.0+] ( 18.0, 28.0] [31] ( 0.0, 8.0+] ( 8.0,1032.0] ( 0.0, 12.0+] [34] ( 12.0, 51.0] ( 0.0, 3.0+] ( 3.0, 733.0] [37] ( 0.0, 83.0+] ( 83.0, 219.0] ( 0.0, 25.0+] [40] ( 25.0,1800.0+] ( 0.0,1401.0+] ( 0.0, 263.0] [43] ( 0.0, 71.0+] ( 71.0, 72.0] ( 0.0, 35.0] [46] ( 0.0, 16.0+] ( 16.0, 852.0] ( 0.0, 16.0] [49] ( 0.0, 17.0+] ( 17.0, 77.0] ( 0.0, 51.0+] [52] ( 51.0,1587.0+] ( 0.0, 23.0+] ( 23.0,1572.0+] [55] ( 0.0, 12.0] ( 0.0, 46.0+] ( 46.0, 100.0] [58] ( 0.0, 19.0+] ( 19.0, 66.0] ( 0.0, 4.5+] [61] ( 4.5, 5.0] ( 0.0, 2.0+] ( 2.0, 53.0] [64] ( 0.0, 41.0+] ( 41.0,1408.0+] ( 0.0, 58.0+] [67] ( 58.0,1322.0+] ( 0.0, 3.0] ( 0.0, 2.0] [70] ( 0.0, 40.0] ( 0.0, 1.0+] ( 1.0, 45.0] [73] ( 0.0, 2.0+] ( 2.0, 996.0] ( 0.0, 21.0+] [76] ( 21.0, 72.0] ( 0.0, 9.0] ( 0.0, 36.0+] [79] ( 36.0,1142.0+] ( 0.0, 83.0+] ( 83.0, 980.0] [82] ( 0.0, 32.0+] ( 32.0, 285.0] ( 0.0, 102.0] [85] ( 0.0, 41.0+] ( 41.0, 188.0] ( 0.0, 3.0] [88] ( 0.0, 10.0+] ( 10.0, 61.0] ( 0.0, 67.0+] [91] ( 67.0, 942.0+] ( 0.0, 149.0] ( 0.0, 21.0+] [94] ( 21.0, 343.0] ( 0.0, 78.0+] ( 78.0, 916.0+] [97] ( 0.0, 3.0+] ( 3.0, 68.0] ( 0.0, 2.0] [100] ( 0.0, 69.0] ( 0.0, 27.0+] ( 27.0, 842.0+] [103] ( 0.0, 33.0+] ( 33.0, 584.0] ( 0.0, 12.0+] [106] ( 12.0, 78.0] ( 0.0, 32.0] ( 0.0, 57.0+] [109] ( 57.0, 285.0] ( 0.0, 3.0+] ( 3.0, 68.0] [112] ( 0.0, 10.0+] ( 10.0, 670.0+] ( 0.0, 5.0+] [115] ( 5.0, 30.0] ( 0.0, 31.0+] ( 31.0, 620.0+] [118] ( 0.0, 4.0+] ( 4.0, 596.0+] ( 0.0, 27.0+] [121] ( 27.0, 90.0] ( 0.0, 5.0+] ( 5.0, 17.0] [124] ( 0.0, 2.0] ( 0.0, 46.0+] ( 46.0, 545.0+] [127] ( 0.0, 21.0] ( 0.0, 210.0+] (210.0, 515.0+] [130] ( 0.0, 67.0+] ( 67.0, 96.0] ( 0.0, 26.0+] [133] ( 26.0, 482.0+] ( 0.0, 6.0+] ( 6.0, 445.0+] [136] ( 0.0, 428.0+] ( 0.0, 32.0+] ( 32.0, 80.0] [139] ( 0.0, 37.0+] ( 37.0, 334.0] ( 0.0, 5.0] [142] ( 0.0, 8.0+] ( 8.0, 397.0+] ( 0.0, 60.0+] [145] ( 60.0, 110.0] ( 0.0, 31.0+] ( 31.0, 370.0+] [148] ( 0.0, 139.0+] (139.0, 207.0] ( 0.0, 160.0+] [151] (160.0, 186.0] ( 0.0, 340.0] ( 0.0, 310.0+] [154] (310.0, 340.0+] ( 0.0, 28.0+] ( 28.0, 265.0+] [157] ( 0.0, 4.0+] ( 4.0, 165.0] ( 0.0, 2.0+] [160] ( 2.0, 16.0] ( 0.0, 13.0+] ( 13.0, 180.0+] [163] ( 0.0, 21.0+] ( 21.0, 131.0+] ( 0.0, 96.0+] [166] ( 96.0, 109.0+] ( 0.0, 21.0] ( 0.0, 38.0+] [169] ( 38.0, 39.0+] ( 0.0, 31.0+] ( 0.0, 11.0+] [172] ( 0.0, 6.0]
2.survreg
拟合参数生存回归模型。 这些是用于时间变量的任意变换的位置尺度模型; 最常见的情况使用对数变换，建立加速失效时间模型。
survreg(formula, data, weights, subset,
na.action, dist="weibull", init=NULL, scale=0,
control,parms=NULL,model=FALSE, x=FALSE,
y=TRUE, robust=FALSE, score=FALSE, ...)
dist
y变量的假设分布。"weibull", "exponential", "gaussian", "logistic","lognormal" and "loglogistic"。
scale
可选的固定值。如果设置<=0，scale将被估计
# survreg's scale = 1/(rweibull shape) # survreg's intercept = log(rweibull scale) # survreg结果中输出的scale与“rweibull scale”不同
control
控制值列表，参考survreg.control()
survreg.control
survreg.control(maxiter=30, rel.tolerance=1e-09,
toler.chol=1e-10, iter.max, debug=0, outer.max=10)
maxiter
最大迭代次数
rel.tolerance
“相对容忍度”来声明收敛
toler.chol
Tolerance to declare Cholesky decomposition singular
iter.max
与maxiter相同
debug
打印调试信息
outer.max
用于选择惩罚参数的外部迭代的最大数目
convergence tolerance
/x(k+1)/=/x(k)/*Tr+Ta
其中：
k 迭代次数；
x(k+1) x的k次迭代计算值；
x(k) x的k次迭代初始值；
Tr 相对误差
Ta绝对误差
// 表示绝对值
因此，从上述公式我们可以得到一个重要结论：设定计算迭代误差的时候，要全面权衡物理量的绝对值大小，同时要衡量收敛迭代值的相对误差，这样才能正确满意的设定自己需要的计算误差！
不能简单的以为绝对误差和相对误差越小越好，也要杜绝tolerance设定时的随意性（按照公式进行合理的设定是一个优秀的过程工程师的基本素质）
survreg使用示例
> library(survival)
> str(ovarian)
'data.frame': 26 obs. of 6 variables:
$ futime : num 59 115 156 421 431 448 464 475 477 563 ...
$ fustat : num 1 1 1 0 1 0 1 1 0 1 ...
$ age : num 72.3 74.5 66.5 53.4 50.3 ...
$ resid.ds: num 2 2 2 2 2 1 2 2 2 1 ...
$ rx : num 1 1 1 2 1 1 2 2 1 2 ...
$ ecog.ps : num 1 1 2 1 1 2 2 2 1 2 ...
> # 拟合指数模型，以下2个模型是一样的
> survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist='weibull',
+ scale=1)
Call:
survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian,
dist = "weibull", scale = 1)
Coefficients:
(Intercept) ecog.ps rx
6.9618376 -0.4331347 0.5815027
Scale fixed at 1 # 注意这里
Loglik(model)= -97.2 Loglik(intercept only)= -98
Chisq= 1.67 on 2 degrees of freedom, p= 0.43
n= 26
> survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian,
+ dist="exponential")
Call:
survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian,
dist = "exponential")
Coefficients:
(Intercept) ecog.ps rx
6.9618376 -0.4331347 0.5815027
Scale fixed at 1 # 注意这里
Loglik(model)= -97.2 Loglik(intercept only)= -98
Chisq= 1.67 on 2 degrees of freedom, p= 0.43
n= 26
> ##########################################
>
> str(lung)
'data.frame': 228 obs. of 10 variables:
$ inst : num 3 3 3 5 1 12 7 11 1 7 ...
$ time : num 306 455 1010 210 883 ...
$ status : num 2 2 1 2 2 1 2 2 2 2 ...
$ age : num 74 68 56 57 60 74 68 71 53 61 ...
$ sex : num 1 1 1 1 1 1 2 2 1 1 ...
$ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
$ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
$ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
$ meal.cal : num 1175 1225 NA 1150 NA ...
$ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...
> # weibull分布,如果设置scale<=0，模型将使用最大似然估计，估计最优的scale
> myfit<-survreg(Surv(time, status) ~ ph.ecog + age + sex, data=lung,dist = "weibull")
> myfit
Call:
survreg(formula = Surv(time, status) ~ ph.ecog + age + sex, data = lung,
dist = "weibull")
Coefficients:
(Intercept) ph.ecog age sex
6.273435252 -0.339638098 -0.007475439 0.401090541
Scale= 0.731109
Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4
Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06
n=227 (1 observation deleted due to missingness)
> summary(myfit)
Call:
survreg(formula = Surv(time, status) ~ ph.ecog + age + sex, data = lung,
dist = "weibull")
Value Std. Error z p
(Intercept) 6.27344 0.45358 13.83 1.66e-43
ph.ecog -0.33964 0.08348 -4.07 4.73e-05
age -0.00748 0.00676 -1.11 2.69e-01
sex 0.40109 0.12373 3.24 1.19e-03
Log(scale) -0.31319 0.06135 -5.11 3.30e-07
Scale= 0.731
Weibull distribution
Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4
Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06
Number of Newton-Raphson Iterations: 5
n=227 (1 observation deleted due to missingness)
# survreg's scale = 1/(rweibull shape)
# survreg's intercept = log(rweibull scale)
# survreg结果中输出的scale与“rweibull scale”不同
##############################################
> # weibull分布，以sex对数据分组，产生2个不同的scale
> myfit1<-survreg(Surv(time, status) ~ ph.ecog + age + strata(sex), data=lung,dist = "weibull")
> summary(myfit1)
Call:
survreg(formula = Surv(time, status) ~ ph.ecog + age + strata(sex),
data = lung, dist = "weibull")
Value Std. Error z p
(Intercept) 6.73235 0.42396 15.880 8.75e-57
ph.ecog -0.32443 0.08649 -3.751 1.76e-04
age -0.00581 0.00693 -0.838 4.02e-01
sex=1 -0.24408 0.07920 -3.082 2.06e-03
sex=2 -0.42345 0.10669 -3.969 7.22e-05
Scale:
sex=1 sex=2
0.783 0.655
Weibull distribution
Loglik(model)= -1137.3 Loglik(intercept only)= -1146.2
Chisq= 17.8 on 2 degrees of freedom, p= 0.00014
Number of Newton-Raphson Iterations: 5
n=227 (1 observation deleted due to missingness)
################################################
# 具有高斯误差的线性回归，以及左删失数据
> survreg(Surv(durable, durable>0, type='left') ~ age + quant,
+ data=tobin, dist='gaussian', scale = 0)
Call:
survreg(formula = Surv(durable, durable > 0, type = "left") ~
age + quant, data = tobin, dist = "gaussian", scale = 0)
Coefficients:
(Intercept) age quant
15.14486636 -0.12905928 -0.04554166
Scale= 5.57254
Loglik(model)= -28.9 Loglik(intercept only)= -29.5
Chisq= 1.1 on 2 degrees of freedom, p= 0.58
n= 20
3.survdiff
测试在两条或更多条生存曲线之间是否存在差异
survdiff(formula, data, subset, na.action, rho=0)
rho = 0 表示log-rank or Mantel-Haenszel检验；
rho = 1 表示Wilcoxon检验
log-rank检验(时序检验)--生存时间分布近似weibull分布或者属于比例风险模型时效率较高
似然比检验(likelihood ratio test)--生存时间分布近似指数分布时效率较高
wilcoxon检验--生存时间分布近似对数正态分布效率较高
survdiff使用示例
> survdiff(Surv(futime, fustat) ~ rx,data=ovarian)
Call:
survdiff(formula = Surv(futime, fustat) ~ rx, data = ovarian)
N Observed Expected (O-E)^2/E (O-E)^2/V
rx=1 13 7 5.23 0.596 1.06
rx=2 13 5 6.77 0.461 1.06
Chisq= 1.1 on 1 degrees of freedom, p= 0.303
> # 用strata来控制协变量的影响
> survdiff(Surv(time, status) ~ pat.karno + strata(inst), data=lung)
Call:
survdiff(formula = Surv(time, status) ~ pat.karno + strata(inst),
data = lung)
n=224, 4 observations deleted due to missingness.
N Observed Expected (O-E)^2/E (O-E)^2/V
pat.karno=30 2 1 0.692 0.13720 0.15752
pat.karno=40 2 1 1.099 0.00889 0.00973
pat.karno=50 4 4 1.166 6.88314 7.45359
pat.karno=60 30 27 16.298 7.02790 9.57333
pat.karno=70 41 31 26.358 0.81742 1.14774
pat.karno=80 50 38 41.938 0.36978 0.60032
pat.karno=90 60 38 47.242 1.80800 3.23078
pat.karno=100 35 21 26.207 1.03451 1.44067
Chisq= 21.4 on 7 degrees of freedom, p= 0.00326
> survdiff(Surv(time, status) ~ pat.karno + sex, data=lung)
Call:
survdiff(formula = Surv(time, status) ~ pat.karno + sex, data = lung)
n=225, 3 observations deleted due to missingness.
N Observed Expected (O-E)^2/E (O-E)^2/V
pat.karno=30, sex=1 1 1 0.246 2.31e+00 2.33e+00
pat.karno=30, sex=2 1 0 0.412 4.12e-01 4.16e-01
pat.karno=40, sex=1 1 1 0.132 5.68e+00 5.72e+00
pat.karno=40, sex=2 1 0 1.204 1.20e+00 1.22e+00
pat.karno=50, sex=1 2 2 0.111 3.23e+01 3.25e+01
pat.karno=50, sex=2 2 2 0.968 1.10e+00 1.11e+00
pat.karno=60, sex=1 18 17 9.173 6.68e+00 7.14e+00
pat.karno=60, sex=2 12 10 6.064 2.56e+00 2.68e+00
pat.karno=70, sex=1 30 23 16.737 2.34e+00 2.65e+00
pat.karno=70, sex=2 11 8 9.527 2.45e-01 2.62e-01
pat.karno=80, sex=1 32 26 24.570 8.32e-02 9.91e-02
pat.karno=80, sex=2 19 13 16.311 6.72e-01 7.55e-01
pat.karno=90, sex=1 31 25 24.992 2.66e-06 3.21e-06
pat.karno=90, sex=2 29 13 24.420 5.34e+00 6.33e+00
pat.karno=100, sex=1 21 15 14.002 7.11e-02 7.91e-02
pat.karno=100, sex=2 14 6 13.131 3.87e+00 4.25e+00
Chisq= 66.1 on 15 degrees of freedom, p= 2.17e-08
anova
计算一个或多个拟合模型的方差（或偏差）表的分析
mod1<-survreg(Surv(time, status) ~ ph.ecog + age + sex, data=lung,dist = "weibull")
mod2<-survreg(Surv(time, status) ~ ph.ecog + age + sex+ph.karno*pat.karno, data=lung,dist = "weibull")
anova(mod1,mod2,test="Chisq")
Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
1 ph.ecog + age + sex 222 2264.877 NA NA NA
2 ph.ecog + age + sex + ph.karno * pat.karno 215 2209.372 = 7 55.50527 1.183848e-09
library(MASS)
# 简洁模型更好
stepAIC(mod1)
stepAIC(mod2)
