Создайте кривую выживаемости с зависимым от времени ковариатом и изменяющимся во времени коэффициентом

Question

Я хочу, чтобы код генерировал кривые выживания в настройках с обоими

• ковариатами, зависящими от времени, и
• коэффициентами, изменяющимися во времени.

Цель состоит в том, чтобы продемонстрировать, как метод выставления счетов влияет на истечение срока действия полиса страхования жизни.Сложность в том, что

1. метод выставления счетов клиентам (счет-фактура или EFT) меняется со временем,
2. влияние метода выставления счетов по истечении времени со временем стирается, и
3. влияние метода выставления счетов на упущение зависит от других ковариат.

После прочтения виньетки по ковариатам, зависящим от времени, я не знаю, как сгенерировать кривые выживания из модели, которая имеет и ковариат, зависящих от времени, и изменяющиеся во времени коэффициенты.

library(survival)

Samp <- data.frame(
  id = c(143,151,680,134),
  time = c(17,16,17,18) ,
  censor= rep(1,4) , 
  covariate = seq(5,20,length.out = 4))
# Lookup provides the values of a tdc
Lookup <- data.frame(
  id =c(rep(134,2),680,143,rep(151,3)) ,
  billing.mode = c("INV",rep("EFT",2),rep("INV",2),"EFT","INV") ,
  switch.time = c(0,3,rep(0,3),2,7))

# create the tdc 
Samp.tdc <- tmerge(data1=Samp,data2=Samp,id=id,
                    lapse=event(time,censor))
Samp.tdc <- tmerge(data1=Samp.tdc,data2=Lookup,id=id,
                    billing.mode=tdc(switch.time,billing.mode))
Samp.tdc$inv = as.numeric(Samp.tdc$billing.mode == "INV")

# the call looks something like this
fit <-coxph(Surv(tstart, tstop, lapse) ~ inv + tt(inv) +
  covariate*inv, data = Samp.tdc, 
            tt = function(x, t, ...) x * t)

Когда я говорю, что хочу создать кривые выживания, я имею в виду предсказанное выживание для фиксированного набора времен и ковариатных значений.Скажем для LpsData ниже.

LpsData <- data.frame(
  tstart = rep(c(0,16,17),times=4),
  tstop = rep(16:18,times=4) ,
  lapse = 0 ,
  covariate = rep(c(10,20),each=3,times=2) ,
  inv = rep(c(0,1),each=6) ,
  curve=rep(c('eft','inv'), each=6)
)

Answer 1

Это довольно сложная проблема, и я лично нахожу возможности пакет survival ограничен в этом отношении. Например, вы должны предварительно укажите функциональную форму изменения времени. Альтернативой является используйте кусочно-экспоненциальные аддитивные модели (ПАММ) , которые можно оценить через mgcv::gam и, следовательно, наследует всю свою гибкость (+ штрафная оценка нелинейных эффектов, включая изменяющиеся во времени эффекты).

В общем, вам нужно решить, какую модель вы хотите установить. Пусть z будет вашим зависящим от времени ковариатом. Чем потенциальные модели могут быть

• линейный ковариатный эффект, линейно изменяющийся во времени, т.е. модель, указанная в Ваш код (mgcv формула: + z * t +)
• нелинейный ковариатный эффект, линейно изменяющийся во времени (формула: + s(z, by = t) +)
• линейный ковариатный эффект, нелинейно изменяющийся во времени (формула: + s(t, by = z) +)
• нелинейный, нелинейно изменяющийся во времени (формула: + te(t, z) +)

Ниже приведен пример использования данных pbc из пакета survival, который также фигурирует в виньетке выживания на зависящих от времени ковариатах (см. также https://adibender.github.io/pammtools/articles/tdcovar.html для сравнение с ПАММ):

library(survival)
library(ggplot2)
theme_set(theme_bw())
library(pammtools)
library(mgcv)

Преобразование данных

Сначала я преобразую данные в формат пошаговых экспоненциальных данных (PED):

pbc <- pbc %>% filter(id <= 312) %>%
  select(id:sex, bili, protime) %>%
  mutate(status = 1L * (status == 2))

## Transform to piece-wise exponential data (PED) format
pbc_ped <- as_ped(
  data = list(pbc, pbcseq),
  formula = Surv(time, status)~. | concurrent(bili, protime, tz_var = "day"),
  id = "id") %>% ungroup()

pbc_ped <- pbc_ped %>%
  mutate(
    log_bili = log(bili),
    log_protime = log(protime))

Fit Поэтапная экспоненциальная аддитивная модель (PAM)

Здесь мне подходит модель с 2 зависящими от времени ковариатами с линейными ковариатами эффекты, нелинейно изменяющиеся во времени (хотя оценки почти линейны из-за штрафа) * ​​1047 *

pbc_pam <- gam(ped_status ~ s(tend, k = 10) + s(tend, by = log_bili) +
  s(tend, by = log_protime),
  data = pbc_ped, family = poisson(), offset = offset)

Прогноз выживания для фиксированных ковариат

Для прогноза I

• создать новый набор данных во все уникальные наблюдаемые моменты времени (для всех неопределенных ковариат будут заданы средние / модусные значения)
• добавить зависящее от времени значение log_bili в каждый момент времени
• добавить прогнозы вероятности выживания + CI, используя add_surv_prob

ndf <- make_newdata(pbc_ped, tend = unique(tend)) %>%
  mutate(log_bili = runif(n(), min(log_bili), max(log_bili))) %>%
  add_surv_prob(pbc_pam) 

График прогнозируемой вероятности выживания

ggplot(ndf, aes(x = tend, y = surv_prob)) +
  geom_surv() +
  geom_ribbon(aes(ymin = surv_lower, ymax = surv_upper), alpha = 0.3) +
  ylim(c(0, 1))

`` `

^{Создано в 2018-12-08 пакетом представ (v0.2.1)}