Практическая работа №7

9. Анализ данных типа времени жизни

Глушков Егор Александрович, гр. 20.М04-мм

---

Данные SurMart.xls о дожитии и длительности безрецидивного периоды больных с глиобластомами.

• time – выживаемость (дней)
• time.c – индикатор цензурирования выживаемости (0 – полные данные, 1 – цензурированные)
• relapse_free – безрецедивный период (дней)
• relapse_free.c – индикатор цензурирования безрецедивного периода
• oper – вид оперативного вмешательства (0 – нет, 1 – открытое, 2 – стереотакс. биопсия, 3 – стереотакс. криотомия)
• speech_disorders – речевые нарушения (0 – нет, 1 – ?, 2 – способность общения не нарушена, 3 – нарушена способность общения, 4 – неспособность к общению) – категориальная переменная
• gen – пол (1 – мужской, 2 – женский) – бинарная переменная
• karn – общее состояние по Карновскому – метрическая переменная (ковариата)

1. Построить кривые Каплана-Мейера

   1. по выживанию,

   2. по длительности безрецедивного периода.

   для групп, отличающихся:

   1. значениями категориальной переменной,

   2. бинарной переменной.

   Привести средние и медианы дожития, проверить значимость различия групп.

2. Внести в модель влияние ковариаты (метрические переменные) и проверить значимость её влияния.

---

Кривые Каплана-Мейера по выживанию

library(readxl)

X <- read_excel("SurMart.xls")

surv_data <- na.omit(X[ , c(4, 3, 6, 5, 17, 1, 7)])
colnames(surv_data) <- c('time', 'time.c', 'relapse_free', 'relapse_free.c', 'speech_disorders', 'gen', 'karn')
surv_data$speech_disorders <- as.factor(surv_data$speech_disorders)
surv_data$gen <- as.factor(surv_data$gen)

summary(surv_data)

      time            time.c        relapse_free    relapse_free.c   speech_disorders gen          karn      
 Min.   :  9.00   Min.   :0.0000   Min.   :   1.0   Min.   :0.0000   0:481            1:342   Min.   :20.00  
 1st Qu.: 59.00   1st Qu.:1.0000   1st Qu.:   8.0   1st Qu.:0.0000   1:  9            2:262   1st Qu.:60.00  
 Median : 69.00   Median :1.0000   Median :  66.5   Median :1.0000   2: 56                    Median :70.00  
 Mean   : 70.18   Mean   :0.8046   Mean   : 440.9   Mean   :0.7053   3: 50                    Mean   :71.84  
 3rd Qu.: 79.00   3rd Qu.:1.0000   3rd Qu.: 476.5   3rd Qu.:1.0000   4:  8                    3rd Qu.:80.00  
 Max.   :336.00   Max.   :1.0000   Max.   :6923.0   Max.   :1.0000                            Max.   :90.00  

Строим кривую дожития (в качестве времени – выживаемость в днях)

library(survival)
library(ggfortify)

time_surv <- Surv(surv_data$time, event = surv_data$time.c)
km_fit <- survfit(time_surv ~ 1, data=surv_data)

# summary(km_fit)
# plot(km_fit)
autoplot(km_fit)

hist(surv_data$time)

print(km_fit, print.rmean = TRUE)

Call: survfit(formula = time_surv ~ 1, data = surv_data)

         n     events     *rmean *se(rmean)     median    0.95LCL    0.95UCL 
    604.00     486.00      74.91       1.68      78.00      78.00      78.00 
    * restricted mean with upper limit =  336 

Согласно усеченному средему и медиане, выживаемость составляет 75-78 суток

Кривые для групп, отличающиеся по категориальной переменной “речевые нарушения”

km_speech <- survfit(time_surv ~ speech_disorders, data=surv_data)
# plot(km_speech, lty=seq(0, 4))
autoplot(km_speech)

print(km_speech, print.rmean = TRUE)

Call: survfit(formula = time_surv ~ speech_disorders, data = surv_data)

                     n events *rmean *se(rmean) median 0.95LCL 0.95UCL
speech_disorders=0 481    403   74.7      0.611     78      78      78
speech_disorders=1   9      9   66.7      5.200     69      58      NA
speech_disorders=2  56     40   81.2     11.877     68      68      69
speech_disorders=3  50     32   87.7     18.321     58      58      68
speech_disorders=4   8      2  150.7     75.662     58      58      NA
    * restricted mean with upper limit =  336 

Заметно, как сильно отличается усеченное среднее для группы 4 и для всех остальных: оно больше почти в два раза (однако мало объектов, чтобы всерьез говорить об этом). Среднее группы 3 также превосходит усеченные средние других групп. В то же время медианное значение наоборот больше у группы 0 и с ростом номера группы лишь уменьшается.

\(H_0: h_1(t) = h_2(t) = ... = h_n(t)\)
\(H_1: h_i(t_0) \neq h_j(t_0)\) хотя бы для одной пары i, j и \(t_0\)

survdiff(time_surv ~ speech_disorders, data=surv_data)

Call:
survdiff(formula = time_surv ~ speech_disorders, data = surv_data)

                     N Observed Expected (O-E)^2/E (O-E)^2/V
speech_disorders=0 481      403  443.795      3.75     53.63
speech_disorders=1   9        9    5.140      2.90      3.51
speech_disorders=2  56       40   25.680      7.99     10.57
speech_disorders=3  50       32   10.654     42.77     49.49
speech_disorders=4   8        2    0.732      2.20      2.26

 Chisq= 70.3  on 4 degrees of freedom, p= 2e-14 

p = 2e-14, различия кривых дожития между группами (по нарушениям речи) значимы

Кривые для групп, отличающиеся по бинарной переменной “пол”

km_gen <- survfit(time_surv ~ gen, data=surv_data)
# plot(km_gen, lty=c(1, 2))
autoplot(km_gen)

print(km_gen, print.rmean = TRUE)

Call: survfit(formula = time_surv ~ gen, data = surv_data)

        n events *rmean *se(rmean) median 0.95LCL 0.95UCL
gen=1 342    263   76.5       2.75     78    77.6      78
gen=2 262    223   73.7       1.58     78    70.0      78
    * restricted mean with upper limit =  336 

Для мужского пола среднее время выживания (в днях) равно 76.5 (усеченное среднее), медианное – 78. Для женского пола среднее – 73.7, медианное – также 78.

\(H_0: h_1(t) = h_2(t) = ... = h_n(t)\)
\(H_1: h_i(t_0) \neq h_j(t_0)\) хотя бы для одной пары i, j и \(t_0\)

survdiff(time_surv ~ gen, data=surv_data)

Call:
survdiff(formula = time_surv ~ gen, data = surv_data)

        N Observed Expected (O-E)^2/E (O-E)^2/V
gen=1 342      263      278     0.839      2.64
gen=2 262      223      208     1.124      2.64

 Chisq= 2.6  on 1 degrees of freedom, p= 0.1 

p = 0.1, различия кривых дожития между полами не значимы на уровне 0.1 > 0.05

Кривые Каплана-Мейера по длительности безрецидивного периода

rel_surv <- Surv(surv_data$relapse_free, event = surv_data$relapse_free.c)
km_fit.2 <- survfit(rel_surv ~ 1, data=surv_data)

# summary(km_fit.2)
# plot(km_fit.2)
autoplot(km_fit.2)

hist(surv_data$relapse_free)

print(km_fit.2, print.rmean = TRUE)

Call: survfit(formula = rel_surv ~ 1, data = surv_data)

         n     events     *rmean *se(rmean)     median    0.95LCL    0.95UCL 
     604.0      426.0      955.4       94.1      150.0       79.0      259.0 
    * restricted mean with upper limit =  6923 

Медианное значение безрецидивного периода в днях составляет 150, тогда как усечённое среднее – 955.4 дня.

Кривые для групп, отличающиеся по категориальной переменной “речевые нарушения”

km_speech.2 <- survfit(rel_surv ~ speech_disorders, data=surv_data)
# plot(km_speech.2, lty=seq(0, 4))
autoplot(km_speech.2)

print(km_speech.2, print.rmean = TRUE)

Call: survfit(formula = rel_surv ~ speech_disorders, data = surv_data)

                     n events  *rmean *se(rmean) median 0.95LCL 0.95UCL
speech_disorders=0 481    350 1042.90     105.73    228      96     335
speech_disorders=1   9      7   71.37      34.91      9       7      NA
speech_disorders=2  56     42  374.21      95.34     53      16     320
speech_disorders=3  50     25 1549.75     677.40     58      19      NA
speech_disorders=4   8      2    8.88       1.05     10      NA      NA
    * restricted mean with upper limit =  6923 

Если не обращать внимания на те группы, где мало наблюдений (1, 4) и потому сложно всерьез брать это в расчет, то для групп 0, 2, 3 усеченные средние сильно отличаются: разница между некоторыми из них 3-5 раз. Медиана группы 1 также намного больше примерно равных медиан групп 2 и 3.

\(H_0: h_1(t) = h_2(t) = ... = h_n(t)\)
\(H_1: h_i(t_0) \neq h_j(t_0)\) хотя бы для одной пары i, j и \(t_0\)

survdiff(rel_surv ~ speech_disorders, data=surv_data)

Call:
survdiff(formula = rel_surv ~ speech_disorders, data = surv_data)

                     N Observed Expected (O-E)^2/E (O-E)^2/V
speech_disorders=0 481      350   370.62     1.147     9.437
speech_disorders=1   9        7     2.99     5.392     5.773
speech_disorders=2  56       42    30.80     4.071     4.580
speech_disorders=3  50       25    20.38     1.049     1.170
speech_disorders=4   8        2     1.22     0.502     0.567

 Chisq= 13.1  on 4 degrees of freedom, p= 0.01 

p = 0.01, различия кривых дожития между группами значимы, так как нулевая гипотеза о равенстве функций интенсивности для разных групп отвергается.

Кривые для групп, отличающиеся по бинарной переменной “пол”

km_gen.2 <- survfit(rel_surv ~ gen, data=surv_data)
# plot(km_gen.2, lty=c(1, 2))
autoplot(km_gen.2)

print(km_gen.2, print.rmean = TRUE)

Call: survfit(formula = rel_surv ~ gen, data = surv_data)

        n events *rmean *se(rmean) median 0.95LCL 0.95UCL
gen=1 342    227    881        111    124      51     287
gen=2 262    199    958        130    210      82     369
    * restricted mean with upper limit =  6923 

Можно лишь заметить значительную разницу медиан для обоих полов: почти в 2 раза длиннее безрецедивный период у представительниц женского пола, в усеченных средних такой разницы не наблюдается.

\(H_0: h_1(t) = h_2(t) = ... = h_n(t)\)
\(H_1: h_i(t_0) \neq h_j(t_0)\) хотя бы для одной пары i, j и \(t_0\)

survdiff(rel_surv ~ gen, data=surv_data)

Call:
survdiff(formula = rel_surv ~ gen, data = surv_data)

        N Observed Expected (O-E)^2/E (O-E)^2/V
gen=1 342      227      227  4.08e-05   9.2e-05
gen=2 262      199      199  4.65e-05   9.2e-05

 Chisq= 0  on 1 degrees of freedom, p= 1 

Различия кривых дожития между полами не значимы, так как нулевая гипотеза о равенстве интенсивностей для разных групп не отвергается с p = 1.

Регрессионная модель Кокса. Влияние ковариаты

summary(coxph(Surv(time, event = time.c) ~ karn, data=surv_data))

Call:
coxph(formula = Surv(time, event = time.c) ~ karn, data = surv_data)

  n= 604, number of events= 486 

          coef exp(coef)  se(coef)      z Pr(>|z|)    
karn -0.126675  0.881020  0.004224 -29.99   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

     exp(coef) exp(-coef) lower .95 upper .95
karn     0.881      1.135    0.8738    0.8883

Concordance= 0.894  (se = 0.007 )
Likelihood ratio test= 655.1  on 1 df,   p=<2e-16
Wald test            = 899.2  on 1 df,   p=<2e-16
Score (logrank) test = 861.1  on 1 df,   p=<2e-16

Полученные доверительные уровни вероятности порядка <2e-16 позволяют отвергнуть гипотезу о незначимости коэффициента, то есть о равенстве коэффициента нулю. Таким образом, коэффициент отличен от нуля, влияние ковариаты (общее состояние по Карновскому) на выживаемость (в днях) значимо.

summary(coxph(Surv(relapse_free, event = relapse_free.c) ~ karn, data=surv_data))

Call:
coxph(formula = Surv(relapse_free, event = relapse_free.c) ~ 
    karn, data = surv_data)

  n= 604, number of events= 426 

          coef exp(coef)  se(coef)      z Pr(>|z|)    
karn -0.017172  0.982974  0.003554 -4.832 1.35e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

     exp(coef) exp(-coef) lower .95 upper .95
karn     0.983      1.017    0.9762    0.9898

Concordance= 0.577  (se = 0.015 )
Likelihood ratio test= 21.58  on 1 df,   p=3e-06
Wald test            = 23.35  on 1 df,   p=1e-06
Score (logrank) test = 23.39  on 1 df,   p=1e-06

Полученные p.values, имеющие порядок e-06, позволяют отвергнуть гипотезу о незначимости коэффициента, то есть о равенстве коэффициента нулю. Таким образом, коэффициент отличен от нуля, влияние ковариаты на длительность безрецидивного периода значимо.

rho <- sum(surv_data$time.c) / sum(surv_data$time); rho

[1] 0.01146454

fitKM <- survfit(Surv(time, event=time.c)~1, data=surv_data)
plot(fitKM, mark.time = FALSE, col=1, conf.int = FALSE, lty=2)
curve(exp(-rho*x), from=0, to=max(surv_data$time), add=TRUE, col=2)