20 votos

Coeficientes dependientes del tiempo en R - ¿cómo hacerlo?

Actualización : Lo siento por otra actualización, pero he encontrado algunas posibles soluciones con polinomios fraccionarios y el riesgo de competencia-paquete que necesito un poco de ayuda con.


El problema

No puedo encontrar una manera fácil de hacer un análisis de coeficiente dependiente del tiempo es en R. Quiero ser capaz de tomar mi coeficiente de variables y convertirlo en un coeficiente dependiente del tiempo (no variable) y luego trazar la variación contra el tiempo:

$\beta_{my\_variable}=\beta_0+\beta_1*t+\beta_2*t^2...$

Posibles soluciones

1) Dividir el conjunto de datos

He mirado este ejemplo (véase la parte 2 de la sesión de laboratorio), pero la creación de un conjunto de datos independiente parece complicada, costosa desde el punto de vista informático y poco intuitiva...

2) Modelos de rango reducido - El paquete coxvc

En paquete coxvc ofrece una forma elegante de resolver el problema. manual . El problema es que el autor ya no está desarrollando el paquete (la última versión es del 23/05/2007), después de algunas conversaciones por correo electrónico he conseguido que el paquete funcione pero una ejecución tardó 5 horas en mi conjunto de datos (140 000 entradas) y da estimaciones extremas al final del periodo. Puede encontrar una versión ligeramente actualizada paquete aquí - Principalmente he actualizado la función de trazado.

Puede que sólo sea cuestión de retoques, pero como el programa no proporciona fácilmente intervalos de confianza y el proceso lleva tanto tiempo, ahora mismo estoy buscando otras soluciones.

3) El paquete timereg

El impresionante paquete timereg también resuelve el problema, pero no estoy seguro de cómo usarlo y no me da una trama fluida.

4) Modelo de tiempo polinómico fraccionario (FPT)

Encontré la excelente disertación de Anika Buchholz sobre "Evaluación de los efectos a largo plazo variables en el tiempo de las terapias y los factores pronósticos" que hace un excelente trabajo cubriendo diferentes modelos. Concluye que FPT propuesto por Sauerbrei et al. parece ser el más apropiado para los coeficientes dependientes del tiempo:

FPT es muy bueno para detectar efectos variables en el tiempo, mientras que el enfoque de Rango Reducido da lugar a modelos demasiado complejos, ya que no incluye la selección de efectos variables en el tiempo.

La investigación parece muy completa pero está un poco fuera de mi alcance. También me pregunto un poco ya que resulta que trabaja con Sauerbrei. Parece sólido sin embargo y supongo que el análisis podría hacerse con el paquete mfp pero no estoy seguro de cómo.

5) El paquete cmprsk

He estado pensando en hacer mi análisis de riesgo competitivo, pero los cálculos me han llevado demasiado tiempo, así que he optado por la regresión de Cox normal. El crr tiene una opción para covariables dependientes del tiempo:

....
cov2        matrix of covariates that will be multiplied 
            by functions of time; if used, often these 
            covariates would also appear in cov1 to give 
            a prop hazards effect plus a time interaction
....

Hay un ejemplo cuadrático, pero no entiendo muy bien dónde aparece el tiempo y no estoy seguro de cómo mostrarlo. También he mirado el archivo test.R pero el ejemplo es básicamente el mismo...

Mi código de ejemplo

He aquí un ejemplo que utilizo para probar las distintas posibilidades

library("survival")
library("timereg")
data(sTRACE)

# Basic cox regression    
surv <- with(sTRACE, Surv(time/365,status==9))
fit1 <- coxph(surv~age+sex+diabetes+chf+vf, data=sTRACE)
check <- cox.zph(fit1)
print(check)
plot(check, resid=F)
# vf seems to be the most time varying

######################################
# Do the analysis with the code from #
# the example that I've found        #
######################################

# Split the dataset according to the splitSurv() from prof. Wesley O. Johnson
# http://anson.ucdavis.edu/~johnson/st222/lab8/splitSurv.ssc
new_split_dataset = splitSuv(sTRACE$time/365, sTRACE$status==9, sTRACE[, grep("(age|sex|diabetes|chf|vf)", names(sTRACE))])

surv2 <- with(new_split_dataset, Surv(start, stop, event))
fit2 <- coxph(surv2~age+sex+diabetes+chf+I(pspline(stop)*vf), data=new_split_dataset)
print(fit2)

######################################
# Do the analysis by just straifying #
######################################
fit3 <- coxph(surv~age+sex+diabetes+chf+strata(vf), data=sTRACE)
print(fit3)

# High computational cost!
# The price for 259 events
sum((sTRACE$status==9)*1)
# ~240 times larger dataset!
NROW(new_split_dataset)/NROW(sTRACE)

########################################
# Do the analysis with the coxvc and   #
# the timecox from the timereg library #
########################################
Ft_1 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=3))
fit_coxvc1 <- coxvc(surv~vf+sex, Ft_1, rank=2, data=sTRACE)

fit_coxvc2 <- coxvc(surv~vf+sex, Ft_1, rank=1, data=sTRACE)

Ft_3 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=5))
fit_coxvc3 <- coxvc(surv~vf+sex, Ft_3, rank=2, data=sTRACE)

layout(matrix(1:3, ncol=1))
my_plotcoxvc <- function(fit, fun="effects"){
    plotcoxvc(fit,fun=fun,xlab='time in years', ylim=c(-1,1), legend_x=.010)
    abline(0,0, lty=2, col=rgb(.5,.5,.5,.5))
    title(paste("B-spline =", NCOL(fit$Ftime)-1, "df and rank =", fit$rank))
}
my_plotcoxvc(fit_coxvc1)
my_plotcoxvc(fit_coxvc2)
my_plotcoxvc(fit_coxvc3)

# Next group
my_plotcoxvc(fit_coxvc1)

fit_timecox1<-timecox(surv~sex + vf, data=sTRACE)
plot(fit_timecox1, xlab="time in years", specific.comps=c(2,3))

El código da como resultado estos gráficos: Comparación de diferentes configuraciones para coxvc y de la coxvc y el timecox tramas. Supongo que los resultados están bien, pero no creo que sea capaz de explicar el gráfico timecox, parece demasiado complejo...

Mis preguntas (actuales)

  • ¿Cómo hago el análisis FPT en R?
  • ¿Cómo se utiliza la covariable temporal en cmprsk?
  • ¿Cómo puedo representar gráficamente el resultado (preferiblemente con intervalos de confianza)?

8voto

aron Puntos 174

@mpiktas se acercó al ofrecer un modelo factible, sin embargo el término que hay que utilizar para la cuadrática en tiempo=t sería I(t^2)) . Esto es así porque en R la interpretación de fórmulas de "^" crea interacciones y no realiza exponenciación, por lo que la interacción de "t" con "t" es simplemente "t". (¿No debería migrarse esto a SO con una etiqueta [r]?).

Para alternativas a este proceso, que me parece algo dudoso a efectos de inferencia, y que probablemente se ajuste a su interés por utilizar las funciones de apoyo de los paquetes rms/Hmisc de Harrell, véase "Regression Modeling Strategies" de Harrell. Menciona (pero sólo de pasada, aunque cita algunos de sus propios trabajos) la construcción de ajustes spline a la escala temporal para modelar riesgos en forma de bañera. Su capítulo sobre modelos paramétricos de supervivencia describe una serie de técnicas de trazado para comprobar los supuestos de riesgos proporcionales y para examinar la linealidad de los efectos logarítmicos de los riesgos estimados en la escala temporal.

Edición: Una opción adicional es utilizar coxph descrito como una "lista opcional de funciones de transformación temporal".

6voto

Max Gordon Puntos 2361

He cambiado la respuesta a esto ya que ni las respuestas de @DWin ni las de @Zach responden completamente a cómo modelar coeficientes variables en el tiempo. Recientemente he escrito un Correo electrónico: sobre esto. Esto es lo esencial.

El concepto central de la Modelo de regresión de Cox es el función de riesgo , $h(t)$ . Se define como:

$$h(t) = \frac{f(t)}{S(t)}$$

Donde el $f(t)$ es el riesgo de que se produzca un suceso en un momento dado mientras el $S(t)$ es la probabilidad de sobrevivir a esa tarifa. El número es, por tanto, una fracción con un rango teórico de $0$ a $\infty$ .

Una característica interesante de la función de riesgo es que podemos incluir observaciones en otros momentos distintos de $time_0$ Por ejemplo, si "Pedro" es operado de una artroplastia de cadera en Inglaterra y llega al cabo de 1 año a Suecia, habrá vivido 1 año cuando decidamos incluirlo. Tenga en cuenta que cualquier paciente que hubiera fallecido antes de eso nunca habría llegado a nuestro conocimiento y, por lo tanto, no podemos incluir a Peter en el $S(t)$ cuando se observa el peligro antes de 1 año. Después de 1 año podemos incluir a Peter.

Al permitir la entrada de sujetos en otros momentos, debemos cambiar la Surv de Surv(time, status) a Surv(start_time, end_time, status) . Mientras que el end_time está fuertemente correlacionada con el resultado start_time está ahora disponible como término de interacción (tal y como se insinuaba en las sugerencias originales). En un entorno normal, el start_time es 0 excepto para unos pocos sujetos que aparecen más tarde, pero si dividimos cada observación en varios periodos, de repente tenemos muchas horas de inicio que no son cero. La única diferencia es que cada observación se produce varias veces y todas, excepto la última, tienen la opción de un resultado no censurado.

La división del tiempo en la práctica

Acabo de publicar en CRAN el Greg que hace que este time-split fácil. Primero comenzaremos con algunas observaciones teóricas:

library(Greg)
test_data <- data.frame(
  id = 1:4,
  time = c(4, 3.5, 1, 5),
  event = c("censored", "dead", "alive", "dead"),
  age = c(62.2, 55.3, 73.7, 46.3),
  date = as.Date(
    c("2003-01-01", 
      "2010-04-01", 
      "2013-09-20",
      "2002-02-23"))
)

Podemos mostrarlo gráficamente con * como indicador de evento:

enter image description here

Si aplicamos el timeSplitter como sigue:

library(dplyr)
split_data <- 
  test_data %>% 
  select(id, event, time, age, date) %>% 
  timeSplitter(by = 2, # The time that we want to split by
               event_var = "event",
               time_var = "time",
               event_start_status = "alive",
               time_related_vars = c("age", "date"))

Obtenemos lo siguiente:

enter image description here

Como puede ver, cada objeto se ha dividido en varios eventos en los que el último intervalo de tiempo contiene el estado real del evento. Esto nos permite ahora construir modelos que tienen simples : términos de interacción (no utilice el * ya que se amplía a time + var + time:var y no nos interesa la hora en sí). No es necesario utilizar el I() aunque si desea comprobar la no linealidad con el tiempo, a menudo creo una variable de interacción temporal independiente a la que añado una spline y luego la muestro utilizando rms::contrast . De todos modos, tu llamada de regresión debería tener ahora este aspecto:

coxp(Surv(start_time, end_time, event) ~ var1 + var2 + var2:time, 
     data = time_split_data)

Utilizando el paquete de supervivencia tt función

También existe una forma de modelar los coeficientes dependientes del tiempo directamente en el paquete de supervivencia utilizando la función tt función. El profesor Therneau ofrece una exhaustiva introducción al tema en su viñeta . Desgraciadamente, en grandes conjuntos de datos esto es difícil de hacer debido a las limitaciones de memoria. Parece que el tt divide el tiempo en trozos muy finos generando en el proceso una matriz enorme.

2voto

Boris Tsirelson Puntos 191

Puede utilizar la función aplicar.rodadura función en PerformanceAnalytics para ejecutar una regresión lineal a través de una ventana móvil, lo que permitirá que sus coeficientes varíen con el tiempo.

Por ejemplo:

library(PerformanceAnalytics)
library(quantmod)
getSymbols(c('AAPL','SPY'), from='01-01-1900')
chart.RollingRegression(Cl(AAPL),Cl(SPY), width=252, attribute='Beta')
#Note: Alpha=y-intercept, Beta=regression coeffient

Esto también funciona con otras funciones.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X