Hasta ahora se ha considerado el caso más simple de datos de supervivencia: cada individuo tiene un único tiempo de supervivencia, un único tipo de evento final y las observaciones son independientes. Sin embargo, en muchas aplicaciones estas condiciones no se cumplen.
Existen al menos dos situaciones importantes:
Datos agrupados o correlacionados: los tiempos de supervivencia pueden no ser independientes si los individuos pertenecen al mismo grupo, por ejemplo familias, escuelas, pueblos, hospitales o unidades experimentales.
Eventos recurrentes: el evento de interés puede ocurrir más de una vez por individuo, por ejemplo convulsiones, recaídas, hospitalizaciones o infecciones.
En estas notas nos concentraremos en la primera situación: tiempos de supervivencia agrupados y modelos que permiten manejar la dependencia dentro de grupos.
Supongamos que los datos están organizados en grupos o conglomerados. Por ejemplo, en estudios familiares, varias personas pertenecen a la misma familia. En estudios clínicos, puede haber dos observaciones por paciente, como un ojo tratado y un ojo no tratado. En estudios multicéntricos, los pacientes pueden estar agrupados por hospital.
La dificultad principal es que los modelos estándar de Cox asumen independencia entre observaciones. Si esta suposición se viola, los estimadores puntuales pueden seguir siendo razonables en ciertos modelos marginales, pero los errores estándar pueden estar mal calculados.
Hay dos enfoques frecuentes:
Ejemplo: estudio Ashkenazi y cáncer de mama
Struewing et al. estudiaron el efecto de mutaciones del gen BRCA sobre el riesgo de cáncer de mama en una población judía Ashkenazi.
El conjunto de datos original incluía mujeres de ascendencia Ashkenazi. A cada mujer se le determinó si era portadora de una mutación BRCA. Además, se recolectó información sobre sus familiares femeninas de primer grado: edad al diagnóstico de cáncer de mama, si lo tuvieron, o edad actual si nunca fueron diagnosticadas.
Para el ejemplo se construyó un subconjunto con 1960 familias, cada una con dos familiares femeninas. Si una familia tenía tres o más familiares, se seleccionaron dos al azar. El conjunto de datos ashkenazi está disponible en el paquete asaur.
Las variables principales son:
famID: identificador de familia;brcancer: indicador de cáncer de mama, con 1 si la familiar tuvo cáncer y 0 si fue censurada;age: edad de inicio del cáncer o edad al momento de entrevista;mutant: indicador de si la sujeto era portadora de mutación BRCA.La variable mutant es común a las dos familiares de una misma familia, porque se refiere al estado mutacional de la sujeto.
library(asaur)
library(survival)
data(ashkenazi, package = "asaur")
ashkenazi[ashkenazi$famID %in% c(1, 9, 94), ] famID brcancer age mutant
1 1 0 73 0
2 1 0 40 0
7 9 0 89 0
8 9 1 60 0
87 94 1 44 1
88 94 0 45 1La preocupación estadística es que las familiares de una misma familia comparten factores genéticos y ambientales. Por tanto, puede no ser adecuado tratarlas como observaciones independientes.
Ejemplo: retinopatía diabética
El Diabetic Retinopathy Study fue un ensayo clínico aleatorizado para evaluar el efecto de la fotocoagulación láser sobre el tiempo hasta la ceguera.
En cada paciente:
También interesaba evaluar si el tipo de diabetes, de inicio temprano o adulto, afectaba el tiempo hasta la ceguera y modificaba el efecto del tratamiento.
Los datos están en el paquete timereg.
library(timereg)
library(survival)
library(coxme)Loading required package: bdsmatrix
Attaching package: 'bdsmatrix'The following object is masked from 'package:base':
backsolvedata(diabetes, package = "timereg")
head(diabetes) id time status trteye treat adult agedx
1 5 46.24967 0 2 1 2 28
2 5 46.27553 0 2 0 2 28
3 14 42.50684 0 1 1 1 12
4 14 31.34145 1 1 0 1 12
5 16 42.30098 0 1 1 1 9
6 16 42.27406 0 1 0 1 9str(diabetes)'data.frame': 394 obs. of 7 variables:
$ id : int 5 5 14 14 16 16 25 25 29 29 ...
$ time : num 46.2 46.3 42.5 31.3 42.3 ...
$ status: int 0 0 0 1 0 0 0 0 0 1 ...
$ trteye: int 2 2 1 1 1 1 2 2 2 2 ...
$ treat : int 1 0 1 0 1 0 1 0 1 0 ...
$ adult : int 2 2 1 1 1 1 1 1 1 1 ...
$ agedx : int 28 28 12 12 9 9 9 9 13 13 ...Cada paciente contribuye con dos observaciones, una por ojo. Por tanto, las observaciones de los dos ojos de un mismo paciente no son independientes.
En el enfoque marginal se mantiene la estructura usual del modelo de riesgos proporcionales de Cox:
\[ h_i(t)=h_0(t)\exp(z_i\beta), \]
pero se reconoce que las observaciones dentro de un mismo grupo pueden estar correlacionadas.
La idea es:
Supongamos primero que hay una sola covariable y que el estimador del coeficiente es \(\hat\beta\). Denotemos por \(\hat V\) la varianza estimada bajo el modelo de Cox estándar, que asume independencia. Entonces el error estándar usual es:
\[ \operatorname{se}(\hat\beta)=\sqrt{\hat V}. \]
Para corregir la varianza por agrupamiento, se definen residuos tipo score.
Para el sujeto \(j\) en el grupo \(i\), un residuo score puede escribirse esquemáticamente como
\[ s_{ij} = \delta_{ij} \left[ z_{ij}-\bar z(t_{ij}) \right]- \sum_{t_u\leq t_{ij}} \left[ z_i-\bar z(t_u) \right] e^{z_i\beta} \left[ \hat H_0(t_u)-\hat H_0(t_{u-1}) \right]. \]
El primer término coincide con el residuo de Schoenfeld en los tiempos de falla. El segundo término incorpora la contribución acumulada del riesgo esperado.
Si hay \(G\) grupos y el grupo \(i\) contiene \(n_i\) observaciones, se define
\[ C = \sum_{i=1}^{G} \sum_{j=1}^{n_i} \sum_{m=1}^{n_i} s_{ij}s_{im}. \]
Entonces una varianza ajustada por conglomerados puede escribirse como
\[ V^*=\hat V^2 C. \]
El error estándar robusto es
\[ \operatorname{se}_{rob}(\hat\beta)=\sqrt{V^*}. \]
Si hay \(q\) covariables, \(\hat{\boldsymbol\beta}\) es un vector y \(\hat V\) es la matriz de covarianzas del modelo de Cox estándar. Los residuos score son vectores \(1\times q\) y
\[ C = \sum_{i=1}^{G} \sum_{j=1}^{n_i} \sum_{m=1}^{n_i} \mathbf{s}_{ij}^{\top}\mathbf{s}_{im}. \]
La matriz de covarianzas robusta queda
\[ V^*=\hat V C \hat V. \]
Este estimador se conoce como estimador sándwich, porque la matriz empírica \(C\) queda entre dos copias de \(\hat V\).
Los errores estándar ajustados son
\[ \operatorname{se}_{rob}(\hat{\boldsymbol\beta}) = \left[ \operatorname{diag}(V^*) \right]^{1/2}. \]
El segundo enfoque consiste en introducir explícitamente una fuente de variación no observada para cada grupo. Esa fuente se llama fragilidad.
La fragilidad puede interpretarse como un efecto aleatorio multiplicativo sobre la función de riesgo.
Verosimilitud completa para datos de supervivencia independientes
Para una observación independiente \((t_i,\delta_i,z_i)\), la contribución a la verosimilitud puede escribirse como
\[ L_i = f(t_i;\beta)^{\delta_i} S(t_i;\beta)^{1-\delta_i}. \]
Como
\[ h(t)=\frac{f(t)}{S(t)}, \]
también puede escribirse como
\[ L_i = h(t_i;\beta)^{\delta_i} S(t_i;\beta). \]
Bajo riesgos proporcionales,
\[ h(t_i;\beta) = h_0(t_i)\exp(z_i\beta). \]
Además,
\[ S(t_i;\beta) = \exp{(-H_0(t_i)\exp(z_i\beta))}, \]
donde
\[ H_0(t_i) = \int_0^{t_i} h_0(v),dv \]
es el riesgo acumulado basal.
Por tanto,
\[ L(\beta) = \prod_{i=1}^{n} \left[ h_0(t_i)\exp(z_i\beta) \right]^{\delta_i} \exp\left[ -H_0(t_i)\exp(z_i\beta) \right]. \]
Fragilidad compartida por grupo
Supongamos ahora que los individuos están organizados en grupos. Denotemos por \(\omega_i\) la fragilidad común del grupo \(i\).
El riesgo para el sujeto \(j\) del grupo \(i\) se modela como
\[ h_{ij}(t_{ij}) = h_0(t_{ij}) \omega_i \exp(z_{ij}\beta). \]
Si \(\omega_i>1\), los sujetos del grupo \(i\) tienen mayor riesgo que el promedio, después de ajustar por covariables observadas. Si \(\omega_i<1\), tienen menor riesgo.
Una elección común es asumir fragilidades gamma:
\[ \omega_i \sim \operatorname{Gamma}\left(\frac{1}{\theta},\theta\right), \]
con densidad
\[ g(\omega;\theta) = \frac{ \omega^{1/\theta-1} \exp(-\omega/\theta) }{ \Gamma(1/\theta)\theta^{1/\theta} }. \]
Esta parametrización tiene
\[ E(\omega_i)=1, \qquad \operatorname{Var}(\omega_i)=\theta. \]
Así, \(\theta\) mide el grado de heterogeneidad no observada entre grupos. Si \(\theta=0\), no habría variación de fragilidad entre grupos.
Si las fragilidades fueran observables, la contribución del sujeto \(j\) del grupo \(i\) sería
\[ L_{ij}(\beta,\theta;\omega_i,t_{ij},\delta_{ij},z_{ij}) = g(\omega_i;\theta) \left[ h_0(t_{ij})\omega_i \exp(z_{ij}\beta) \right]^{\delta_{ij}} \exp\left[ -H_0(t_{ij})\omega_i\exp(z_{ij}\beta) \right]. \]
La verosimilitud completa sería
\[ L(\beta,\theta) = \prod_{i=1}^{G} \prod_{j=1}^{n_i} L_{ij}(\beta,\theta;\omega_i,t_{ij},\delta_{ij},z_{ij}). \]
En la práctica, las fragilidades \(\omega_i\) son variables latentes: se supone que existen, pero no se observan directamente. Por eso se requieren métodos de estimación más complejos, como el algoritmo EM.
Algoritmo EM
El algoritmo EM alterna dos pasos:
El procedimiento se repite hasta convergencia.
En modelos paramétricos de supervivencia, plantear el EM puede ser relativamente directo. En el modelo de Cox es más complejo, porque además se debe actualizar el riesgo basal \(h_0(t)\) o el riesgo acumulado basal \(H_0(t)\).
Fragilidad normal
Otra posibilidad es escribir
\[ \omega_i=\exp(\sigma u_i), \]
donde
\[ u_i\sim N(0,1). \]
Entonces
\[ h_{ij}(t_{ij}) = h_0(t_{ij}) \exp(z_{ij}\beta+\sigma u_i). \]
Esta formulación ubica los efectos fijos y aleatorios en la misma escala lineal del predictor del modelo de Cox.
Ejemplo con ashkenazi: conglomerados familiares
Primero ajustamos un modelo de Cox estándar para predecir la edad de inicio del cáncer de mama en función del estado mutacional de la probando (sujeto de prueba).
library(asaur)
library(survival)
data(ashkenazi, package = "asaur")
mod_cox <- coxph(
Surv(age, brcancer) ~ mutant,
data = ashkenazi
)
summary(mod_cox)Call:
coxph(formula = Surv(age, brcancer) ~ mutant, data = ashkenazi)
n= 3920, number of events= 473
coef exp(coef) se(coef) z Pr(>|z|)
mutant 1.1907 3.2895 0.1984 6.002 1.95e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
mutant 3.29 0.304 2.23 4.853
Concordance= 0.523 (se = 0.006 )
Likelihood ratio test= 25.92 on 1 df, p=4e-07
Wald test = 36.03 on 1 df, p=2e-09
Score (logrank) test = 40.49 on 1 df, p=2e-10El modelo es
\[ h_i(t) = h_0(t)\exp(\beta\cdot \texttt{mutant}_i). \]
En el texto original se reporta:
\[ \hat\beta=1.1907, \qquad \exp(\hat\beta)=3.2895, \qquad \operatorname{se}(\hat\beta)=0.1984. \]
Interpretación:
\[ \exp(\hat\beta)=3.29. \]
Las familiares de mujeres portadoras de mutación BRCA tienen un riesgo estimado de desarrollar cáncer de mama aproximadamente 3.29 veces el de las familiares de mujeres no portadoras, bajo el modelo de Cox estándar.
La log-verosimilitud parcial del modelo nulo y del modelo con mutant puede extraerse así:
mod_cox$loglik[1] -3579.707 -3566.745G2 <- 2 * diff(mod_cox$loglik)
p_value <- pchisq(G2, df = 1, lower.tail = FALSE)
G2[1] 25.92229p_value[1] 3.554404e-07En el texto se reporta:
\[ \ell_0=-3579.707, \qquad \ell_1=-3566.745. \]
Por tanto,
\[ G^2 =2(\ell_1-\ell_0)=2(-3566.745+3579.707) =25.924. \]
Comparado con una \(\chi^2_1\), el valor-p es muy pequeño. Esto confirma la importancia de incluir mutant en el modelo.
Modelo marginal con errores estándar robustos por familia
Para corregir los errores estándar por agrupamiento familiar se usa cluster(famID).
mod_cox_cluster <- coxph(
Surv(age, brcancer) ~ mutant + cluster(famID),
data = ashkenazi
)
summary(mod_cox_cluster)Call:
coxph(formula = Surv(age, brcancer) ~ mutant, data = ashkenazi,
cluster = famID)
n= 3920, number of events= 473
coef exp(coef) se(coef) robust se z Pr(>|z|)
mutant 1.1907 3.2895 0.1984 0.2023 5.886 3.96e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
mutant 3.29 0.304 2.213 4.89
Concordance= 0.523 (se = 0.006 )
Likelihood ratio test= 25.92 on 1 df, p=4e-07
Wald test = 34.64 on 1 df, p=4e-09
Score (logrank) test = 40.49 on 1 df, p=2e-10, Robust = 13.09 p=3e-04
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).En el texto se reporta:
\[ \hat\beta=1.1907, \qquad \exp(\hat\beta)=3.2895, \]
con error estándar usual
\[ 0.1984 \]
y error estándar robusto
\[ 0.2023. \]
El estimador puntual no cambia, pero el error estándar robusto es ligeramente mayor. Esto sugiere que el agrupamiento familiar tiene un efecto pequeño en este ejemplo. La asociación sigue siendo altamente significativa.
Modelo de fragilidad gamma
También se puede modelar la dependencia familiar mediante fragilidad gamma:
mod_frail <- coxph(
Surv(age, brcancer) ~ mutant + frailty(famID),
data = ashkenazi
)
summary(mod_frail)Call:
coxph(formula = Surv(age, brcancer) ~ mutant + frailty(famID),
data = ashkenazi)
n= 3920, number of events= 473
coef se(coef) se2 Chisq DF p
mutant 1.272 0.2317 0.2004 30.13 1.0 4.0e-08
frailty(famID) 221.50 211.6 3.1e-01
exp(coef) exp(-coef) lower .95 upper .95
mutant 3.567 0.2804 2.265 5.617
Iterations: 6 outer, 26 Newton-Raphson
Variance of random effect= 0.4937594 I-likelihood = -3564.3
Degrees of freedom for terms= 0.7 211.6
Concordance= 0.938 (se = 0.004 )
Likelihood ratio test= 400.2 on 212.4 df, p=1e-13El modelo es
\[ h_{ij}(t) = h_0(t) \omega_i \exp(\beta\cdot \texttt{mutant}_{ij}), \]
donde \(\omega_i\) es la fragilidad de la familia \(i\).
En el output se reporta aproximadamente:
\[ \hat\beta=1.272. \]
Por tanto,
\[ \exp(1.272)\approx 3.57. \]
El efecto de mutant sigue siendo similar al obtenido con los modelos anteriores.
Por defecto, frailty() usa una distribución gamma. Para usar fragilidad normal se puede especificar:
mod_frail_normal <- coxph(
Surv(age, brcancer) ~ mutant + frailty(famID, distribution = "gaussian"),
data = ashkenazi
)
summary(mod_frail_normal)Call:
coxph(formula = Surv(age, brcancer) ~ mutant + frailty(famID,
distribution = "gaussian"), data = ashkenazi)
n= 3920, number of events= 473
coef se(coef) se2 Chisq DF p
mutant 1.237 0.221 0.1991 31.35 1.0 2.2e-08
frailty(famID, distributi 165.96 152.1 2.1e-01
exp(coef) exp(-coef) lower .95 upper .95
mutant 3.446 0.2902 2.235 5.314
Iterations: 6 outer, 22 Newton-Raphson
Variance of random effect= 0.3572861
Degrees of freedom for terms= 0.8 152.1
Concordance= 0.938 (se = 0.004 )
Likelihood ratio test= 342.5 on 152.9 df, p=<2e-16Modelo mixto de Cox con coxme
El paquete coxme permite ajustar modelos de Cox con efectos aleatorios.
library(coxme)
mod_coxme <- coxme(
Surv(age, brcancer) ~ mutant + (1 | famID),
data = ashkenazi
)
summary(mod_coxme)Mixed effects coxme model
Formula: Surv(age, brcancer) ~ mutant + (1 | famID)
Data: ashkenazi
events, n = 473, 3920
Random effects:
group variable sd variance
1 famID Intercept 0.5912135 0.3495334
Chisq df p AIC BIC
Integrated loglik 30.17 2.0 2.81e-07 26.17 17.85
Penalized loglik 336.37 150.1 2.22e-16 36.16 -588.13
Fixed effects:
coef exp(coef) se(coef) z p
mutant 1.2366 3.4439 0.2205 5.61 2.06e-08La fórmula
Surv(age, brcancer) ~ mutant + (1 | famID)indica que:
mutant es un efecto fijo;(1 | famID) es un intercepto aleatorio por familia.El modelo puede escribirse como
\[ h_{ij}(t) = h_0(t) \exp(\beta \cdot \texttt{mutant}_{ij}+u_i), \]
donde \(u_i\) es el efecto aleatorio de la familia \(i\).
En el output se reporta:
\[ \hat\beta=1.2366, \qquad \exp(\hat\beta)=3.4439, \qquad \operatorname{se}(\hat\beta)=0.2205. \]
También se reporta una varianza del efecto aleatorio familiar de aproximadamente
\[ \widehat{\operatorname{Var}}(u_i)=0.3495, \]
y desviación estándar
\[ \widehat{\operatorname{sd}}(u_i)=0.5912. \]
Interpretación:
mutant, capturada por el efecto aleatorio.Para comparar el modelo mixto con el modelo de Cox ordinario se puede usar la log-verosimilitud integrada. En el texto:
\[ \ell_{\text{cox ordinario}}=-3566.745, \qquad \ell_{\text{integrada}}=-3564.622. \]
Entonces
\[ G^2 = 2(-3564.622+3566.745)= 4.246. \]
Un cálculo aproximado con una \(\chi^2_1\) da:
pchisq(4.246, df = 1, lower.tail = FALSE)El valor-p aproximado es
\[ 0.039. \]
Esto sugiere que el efecto aleatorio familiar mejora el ajuste. Sin embargo, las pruebas para componentes de varianza deben interpretarse con cautela, porque el valor nulo de una varianza está en la frontera del espacio paramétrico.
Ejemplo con diabetes: dos ojos dentro de cada paciente
Ahora se analiza el tiempo hasta la ceguera en pacientes con retinopatía diabética. Cada paciente aporta dos ojos, por lo que las observaciones están agrupadas por paciente.
El objetivo es evaluar:
library(timereg)
library(survival)
library(coxme)
data(diabetes, package = "timereg")
head(diabetes) id time status trteye treat adult agedx
1 5 46.24967 0 2 1 2 28
2 5 46.27553 0 2 0 2 28
3 14 42.50684 0 1 1 1 12
4 14 31.34145 1 1 0 1 12
5 16 42.30098 0 1 1 1 9
6 16 42.27406 0 1 0 1 9str(diabetes)'data.frame': 394 obs. of 7 variables:
$ id : int 5 5 14 14 16 16 25 25 29 29 ...
$ time : num 46.2 46.3 42.5 31.3 42.3 ...
$ status: int 0 0 0 1 0 0 0 0 0 1 ...
$ trteye: int 2 2 1 1 1 1 2 2 2 2 ...
$ treat : int 1 0 1 0 1 0 1 0 1 0 ...
$ adult : int 2 2 1 1 1 1 1 1 1 1 ...
$ agedx : int 28 28 12 12 9 9 9 9 13 13 ...mod_diabetes <- coxme(
Surv(time, status) ~ treat * factor(adult) + (1 | id),
data = diabetes
)
summary(mod_diabetes)Mixed effects coxme model
Formula: Surv(time, status) ~ treat * factor(adult) + (1 | id)
Data: diabetes
events, n = 155, 394
Random effects:
group variable sd variance
1 id Intercept 0.9171924 0.8412418
Chisq df p AIC BIC
Integrated loglik 41.13 4.00 2.521e-08 33.13 20.96
Penalized loglik 213.26 77.99 1.654e-14 57.28 -180.07
Fixed effects:
coef exp(coef) se(coef) z p
treat -0.4998 0.6067 0.2254 -2.22 0.02661
factor(adult)2 0.3995 1.4910 0.2457 1.63 0.10395
treat:factor(adult)2 -0.9681 0.3798 0.3616 -2.68 0.00743El modelo es
\[ h_{ij}(t) = h_0(t) \exp{ \beta_1\texttt{treat}_{ij} + \beta_2\texttt{adult}_{i} + \beta_3(\texttt{treat}_{ij}\times \texttt{adult}_{i}) + u_i }, \]
donde \(u_i\) es un efecto aleatorio del paciente \(i\). Este efecto aleatorio captura la correlación entre los dos ojos de un mismo paciente.
En el output se reportan los siguientes coeficientes:
\[ \hat\beta_{\texttt{treat}}=-0.4998, \qquad \exp(\hat\beta_{\texttt{treat}})=0.6067. \]
Esto indica que, en el grupo de referencia para adult, el tratamiento reduce el riesgo de ceguera. Más específicamente, el riesgo instantáneo del ojo tratado es aproximadamente 0.61 veces el riesgo del ojo no tratado, manteniendo lo demás constante.
También se reporta:
\[ \hat\beta_{\texttt{adult}}=0.3995, \qquad \exp(\hat\beta_{\texttt{adult}})=1.4910. \]
Esto sugiere que los pacientes con diabetes de inicio adulto tienen mayor riesgo de ceguera temprana que los del grupo de referencia, aunque el valor-p reportado es aproximadamente 0.10.
La interacción es
\[ \hat\beta_{\texttt{treat:adult}}=-0.9681, \qquad \exp(\hat\beta_{\texttt{treat:adult}})=0.3798. \]
Esto indica que el efecto protector del tratamiento es mayor en pacientes con diabetes de inicio adulto.
Para obtener las razones de riesgo específicas:
b <- fixef(mod_diabetes)
b treat factor(adult)2 treat:factor(adult)2
-0.4997742 0.3994717 -0.9680873 HR_treat_ref <- exp(b["treat"])
HR_interaction <- exp(b["treat:factor(adult)2"])
HR_treat_adult <- exp(b["treat"] + b["treat:factor(adult)2"])
HR_treat_ref treat
0.6066676 HR_interactiontreat:factor(adult)2
0.3798088 HR_treat_adult treat
0.2304177 Con los valores reportados:
\[ \text{HR}_{\text{tratamiento, referencia}} =\exp(-0.4998)= 0.6067. \]
Para el grupo adulto:
\[ \text{HR}_{\text{tratamiento, adulto}} =\exp(-0.4998-0.9681)= \exp(-1.4679) \approx 0.230. \]
Por tanto, el tratamiento parece reducir mucho más el riesgo de ceguera en el grupo de diabetes de inicio adulto.
El modelo también reporta una varianza del efecto aleatorio por paciente:
\[ \widehat{\operatorname{Var}}(u_i)=0.8412, \]
con desviación estándar
\[ \widehat{\operatorname{sd}}(u_i)=0.9172. \]
Esto indica una heterogeneidad importante entre pacientes y justifica considerar la dependencia entre ojos de un mismo individuo.
Modelo marginal robusto
Un modelo marginal robusto estima el efecto promedio poblacional de las covariables y corrige los errores estándar por agrupamiento.
Ejemplo con diabetes:
mod_diabetes_cluster <- coxph(
Surv(time, status) ~ treat * factor(adult) + cluster(id),
data = diabetes
)
summary(mod_diabetes_cluster)Call:
coxph(formula = Surv(time, status) ~ treat + factor(adult) +
treat:factor(adult), data = diabetes, cluster = id)
n= 394, number of events= 155
coef exp(coef) se(coef) robust se z Pr(>|z|)
treat -0.4248 0.6539 0.2177 0.1850 -2.296 0.0217 *
factor(adult)2 0.3412 1.4067 0.1992 0.1957 1.743 0.0813 .
treat:factor(adult)2 -0.8464 0.4290 0.3509 0.3035 -2.788 0.0053 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
treat 0.6539 1.5293 0.4550 0.9396
factor(adult)2 1.4067 0.7109 0.9585 2.0644
treat:factor(adult)2 0.4290 2.3311 0.2366 0.7777
Concordance= 0.613 (se = 0.019 )
Likelihood ratio test= 28.49 on 3 df, p=3e-06
Wald test = 34.93 on 3 df, p=1e-07
Score (logrank) test = 28.44 on 3 df, p=3e-06, Robust = 30.31 p=1e-06
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).Este modelo mantiene la forma
\[ h_{ij}(t) = h_0(t) \exp{ \beta_1\texttt{treat}_{ij} + \beta_2\texttt{adult}_i + \beta_3(\texttt{treat}_{ij}\times \texttt{adult}_i) }, \]
pero calcula errores estándar robustos por paciente.
Modelo de fragilidad
Un modelo de fragilidad introduce un efecto aleatorio por grupo:
\[ h_{ij}(t) = h_0(t) \exp{ \mathbf{x}_{ij}^{\top}\boldsymbol\beta + u_i }. \]
Este enfoque modela explícitamente la heterogeneidad no observada entre grupos.
Diferencia conceptual
El modelo marginal responde una pregunta de tipo poblacional:
¿Cuál es el efecto promedio de la covariable en la población?
El modelo de fragilidad responde una pregunta condicional al grupo:
¿Cuál es el efecto de la covariable dentro de grupos con fragilidad similar?
Ambos enfoques pueden ser útiles, pero no siempre tienen la misma interpretación.
En análisis de supervivencia clásico se suele estudiar un único desenlace bien definido, por ejemplo muerte, recaída o falla. En muchas aplicaciones, sin embargo, un individuo puede experimentar distintos tipos de eventos, pero solo se observa el primero que ocurre.
Por ejemplo, en pacientes con cáncer de próstata puede interesar el tiempo desde el diagnóstico hasta:
A esta situación se le llama riesgos competitivos. Los eventos compiten porque, una vez que ocurre uno de ellos, los demás ya no pueden observarse para ese individuo.
El problema de tratar eventos competitivos como censura
Si interesa estudiar el tiempo hasta muerte por cáncer de próstata, una estrategia simple sería considerar las muertes por otras causas como censuras. Esto permite usar directamente Kaplan-Meier o modelos de Cox estándar.
Sin embargo, esta estrategia requiere que la censura sea independiente del evento de interés. En riesgos competitivos, esta suposición suele ser dudosa. Por ejemplo, edad, comorbilidades o estado general de salud pueden afectar tanto la probabilidad de morir por cáncer de próstata como la probabilidad de morir por otras causas.
Además, con datos de riesgos competitivos no es posible probar directamente si los tiempos potenciales a las distintas causas son independientes. La evaluación de esta suposición requiere información externa o argumentos científicos sobre las causas de muerte.
Por ello, los análisis que tratan eventos competitivos como censura pueden tener una interpretación ambigua o hipotética.
Supóngase una variable status con la siguiente codificación:
status = 0: censurado;status = 1: muerte por cáncer de próstata;status = 2: muerte por otras causas.Podemos crear indicadores de causa:
library(survival)
library(asaur)
data(prostateSurvival, package = "asaur")
prostateSurvival <- within(prostateSurvival, {
status.prost <- as.numeric(status == 1)
status.other <- as.numeric(status == 2)
})
prostateSurvival.highrisk <- subset(
prostateSurvival,
grade == "poor" & stage == "T2" & ageGroup == "80+"
)
head(prostateSurvival.highrisk) grade stage ageGroup survTime status status.other status.prost
13 poor T2 80+ 21 0 0 0
38 poor T2 80+ 105 0 0 0
41 poor T2 80+ 2 1 0 1
47 poor T2 80+ 67 2 1 0
78 poor T2 80+ 2 0 0 0
93 poor T2 80+ 60 2 1 0Luego se pueden ajustar dos curvas de Kaplan-Meier:
km_prostate <- survfit(
Surv(survTime, status.prost) ~ 1,
data = prostateSurvival.highrisk
)
km_other <- survfit(
Surv(survTime, status.other) ~ 1,
data = prostateSurvival.highrisk
)Para graficar:
plot(
km_prostate$time / 12, 1 - km_prostate$surv,
type = "s", ylim = c(0, 1), lwd = 2, col = "blue",
xlab = "Años desde el diagnóstico",
ylab = "Probabilidad estimada"
)
lines(
km_other$time / 12, km_other$surv,
type = "s", col = "darkgreen", lwd = 2
)
legend(
"right",
legend = c("1 - KM: muerte por próstata",
"KM: supervivencia frente a otras causas"),
col = c("blue", "darkgreen"),
lwd = 2,
bty = "n"
)
El problema es que estas curvas no representan probabilidades reales simultáneamente observables. En el ejemplo descrito, a 10 años la probabilidad estimada de morir por cáncer de próstata es aproximadamente 0.46 y la probabilidad de morir por otras causas es aproximadamente 0.88. La suma excede 1, lo cual muestra que estas estimaciones no pueden interpretarse como probabilidades reales de muerte por una u otra causa para un mismo paciente.
Estas curvas podrían interpretarse hipotéticamente como probabilidades de morir por una causa si la otra causa fuera eliminada, pero eso requiere asumir independencia entre causas.
Supóngase que existen \(K\) causas mutuamente excluyentes de falla. Cada individuo puede experimentar a lo sumo una causa.
La función de riesgo específica por causa para la causa \(j\) se define como
\[ h_j(t) = \lim_{\Delta t \to 0} \frac{ \Pr(t<T<t+\Delta t,\ C=j \mid T>t) }{ \Delta t }. \]
Esta cantidad mide la tasa instantánea de ocurrencia de la causa \(j\), condicionada a que el individuo sigue libre de cualquier evento justo antes de \(t\).
La función de riesgo total es la suma de los riesgos específicos:
\[ h(t)=\sum_{j=1}^{K} h_j(t). \]
Es decir, el riesgo instantáneo de morir por cualquier causa es la suma de los riesgos instantáneos de morir por cada causa específica.
Función de incidencia acumulada
Para cada causa \(j\), se define la función de incidencia acumulada, también llamada función de subdistribución, como
\[ F_j(t)=\Pr(T\leq t,\ C=j). \]
Esta es la probabilidad acumulada de haber experimentado la causa \(j\) antes o en el tiempo \(t\).
Puede escribirse como
\[ F_j(t) = \int_0^t h_j(u)S(u)du, \]
donde \(S(u)\) es la probabilidad de seguir libre de cualquier evento justo antes de \(u\).
A diferencia de una función de distribución ordinaria, \(F_j(t)\) no necesariamente tiende a 1. En cambio,
\[ F_j(\infty)=\Pr(C=j). \]
Es decir, su límite es la probabilidad final de morir por la causa \(j\).
Estimación no paramétrica de la incidencia acumulada
Supóngase que hay tiempos de falla ordenados
\[ t_1,t_2,\ldots,t_D. \]
En el tiempo \(t_i\), sean:
El estimador del riesgo específico por causa \(k\) en \(t_i\) es
\[ \widehat h_k(t_i)=\frac{d_{ik}}{n_i}. \]
El estimador del riesgo total es
\[ \widehat h(t_i)=\frac{d_i}{n_i}. \]
La probabilidad estimada de fallar por la causa \(k\) exactamente en \(t_i\) es
\[ \widehat S(t_{i-1})\widehat h_k(t_i), \]
donde \(\widehat S(t_{i-1})\) es la probabilidad estimada de seguir libre de cualquier evento justo antes de \(t_i\).
Por tanto, la incidencia acumulada para la causa \(k\) se estima mediante
\[ \widehat F_k(t) = \sum_{t_i\leq t} \widehat S(t_{i-1}) \widehat h_k(t_i). \]
Esta fórmula muestra que, para estimar la probabilidad acumulada de morir por una causa específica, no basta con estimar la supervivencia frente a esa causa. Hay que considerar también la probabilidad de seguir libre de cualquier evento hasta cada tiempo.
Ejemplo hipotético con seis individuos
Se consideran seis observaciones y dos causas de falla:
status = 0: censura;status = 1: falla por causa 1;status = 2: falla por causa 2.tt <- c(2, 7, 5, 3, 4, 6)
status <- c(1, 2, 1, 2, 0, 0)
datos <- data.frame(
id = seq_along(tt),
tiempo = tt,
estado = status
)
datos id tiempo estado
1 1 2 1
2 2 7 2
3 3 5 1
4 4 3 2
5 5 4 0
6 6 6 0Primero se estima la supervivencia global, considerando cualquier causa como evento:
library(survival)
status_any <- as.numeric(status >= 1)
km_any <- survfit(Surv(tt, status_any) ~ 1)
summary(km_any)Call: survfit(formula = Surv(tt, status_any) ~ 1)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
2 6 1 0.833 0.152 0.583 1
3 5 1 0.667 0.192 0.379 1
5 3 1 0.444 0.222 0.167 1
7 1 1 0.000 NaN NA NATambién se pueden calcular manualmente las incidencias acumuladas:
event_times <- sort(unique(tt[status > 0]))
S_before <- 1
F1 <- 0
F2 <- 0
tabla_ci <- data.frame()
for (time in event_times) {
n_risk <- sum(tt >= time)
d1 <- sum(tt == time & status == 1)
d2 <- sum(tt == time & status == 2)
d_any <- d1 + d2
h1 <- d1 / n_risk
h2 <- d2 / n_risk
h_any <- d_any / n_risk
F1 <- F1 + S_before * h1
F2 <- F2 + S_before * h2
S_after <- S_before * (1 - h_any)
tabla_ci <- rbind(
tabla_ci,
data.frame(
tiempo = time,
n_risk = n_risk,
eventos_1 = d1,
eventos_2 = d2,
eventos_total = d_any,
S_before = S_before,
h1 = h1,
h2 = h2,
CI_1 = F1,
CI_2 = F2,
S_after = S_after
)
)
S_before <- S_after
}
tabla_ci tiempo n_risk eventos_1 eventos_2 eventos_total S_before h1 h2
1 2 6 1 0 1 1.0000000 0.1666667 0.0
2 3 5 0 1 1 0.8333333 0.0000000 0.2
3 5 3 1 0 1 0.6666667 0.3333333 0.0
4 7 1 0 1 1 0.4444444 0.0000000 1.0
CI_1 CI_2 S_after
1 0.1666667 0.0000000 0.8333333
2 0.1666667 0.1666667 0.6666667
3 0.3888889 0.1666667 0.4444444
4 0.3888889 0.6111111 0.0000000Los valores esperados son:
| Tiempo | En riesgo | Eventos causa 1 | Eventos causa 2 | Supervivencia antes | \(h_1\) | \(h_2\) | \(F_1\) | \(F_2\) |
|---|---|---|---|---|---|---|---|---|
| 2 | 6 | 1 | 0 | 1.000 | \(1/6\) | 0 | 0.167 | 0.000 |
| 3 | 5 | 0 | 1 | 0.833 | 0 | \(1/5\) | 0.167 | 0.167 |
| 5 | 3 | 1 | 0 | 0.667 | \(1/3\) | 0 | 0.389 | 0.167 |
| 7 | 1 | 0 | 1 | 0.444 | 0 | 1 | 0.389 | 0.611 |
Por ejemplo, la probabilidad de un evento de tipo 1 en \(t=2\) es
\[ 1\times \frac{1}{6}=0.167. \]
La probabilidad de un evento de tipo 1 en \(t=5\) es
\[ 0.667\times \frac{1}{3}=0.222. \]
Entonces,
\[ \widehat F_1(5)=0.167+0.222=0.389. \]
Una forma directa de obtener estos resultados es usar Cuminc() del paquete mstate:
library(mstate)
ci <- Cuminc(time = tt, status = status)
ci time Surv CI.1 CI.2 seSurv seCI.1 seCI.2
1 2 0.8333333 0.1666667 0.0000000 0.1521452 0.1521452 0.0000000
2 3 0.6666667 0.1666667 0.1666667 0.1924501 0.1521452 0.1521452
3 5 0.4444444 0.3888889 0.1666667 0.2222222 0.2187224 0.1521452
4 7 0.0000000 0.3888889 0.6111111 0.0000000 0.2187224 0.2187224Incidencia acumulada en los datos de cáncer de próstata
Para el subconjunto de pacientes de alto riesgo se puede estimar la incidencia acumulada de las dos causas:
library(mstate)
ci_prostate <- Cuminc(
time = prostateSurvival.highrisk$survTime,
status = prostateSurvival.highrisk$status
)
head(ci_prostate) time Surv CI.1 CI.2 seSurv seCI.1 seCI.2
1 0 1.0000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
2 1 0.9940758 0.000000000 0.005924171 0.002641510 0.000000000 0.002641510
3 2 0.9880511 0.006024702 0.005924171 0.003756151 0.002686199 0.002641510
4 3 0.9843644 0.008482541 0.007153090 0.004303185 0.003192654 0.002910105
5 4 0.9831168 0.009730151 0.007153090 0.004474935 0.003423736 0.002910105
6 5 0.9780751 0.014771775 0.007153090 0.005112934 0.004233796 0.002910105La salida contiene columnas como:
time: tiempo, en meses;Surv: estimación de Kaplan-Meier para supervivencia frente a cualquier causa;CI.1: incidencia acumulada de muerte por cáncer de próstata;CI.2: incidencia acumulada de muerte por otras causas;seSurv: error estándar de la supervivencia;seCI.1: error estándar de CI.1;seCI.2: error estándar de CI.2.Para graficar ambas incidencias acumuladas:
times <- ci_prostate$time / 12
plot(
times, ci_prostate$CI.1,
type = "s", lwd = 2, col = "blue",
ylim = c(0, 1),
xlab = "Tiempo desde el diagnóstico, en años",
ylab = "Incidencia acumulada"
)
lines(
times, ci_prostate$CI.2,
type = "s", lwd = 2, col = "darkgreen"
)
legend(
"topleft",
legend = c("Muerte por cáncer de próstata", "Muerte por otras causas"),
col = c("blue", "darkgreen"),
lwd = 2,
bty = "n"
)
Estas curvas representan probabilidades reales estimadas:
\[ \widehat F_1(t)=\Pr(T\leq t,\ C=1), \]
\[ \widehat F_2(t)=\Pr(T\leq t,\ C=2). \]
Una representación común es el gráfico apilado:
plot(
times, ci_prostate$CI.1,
type = "s", lwd = 2, col = "blue",
ylim = c(0, 1),
xlab = "Tiempo desde el diagnóstico, en años",
ylab = "Probabilidad acumulada"
)
lines(
times, ci_prostate$CI.1 + ci_prostate$CI.2,
type = "s", lwd = 2, col = "darkgreen"
)
legend(
"topleft",
legend = c(
"Muerte por cáncer de próstata",
"Muerte por cualquier causa"
),
col = c("blue", "darkgreen"),
lwd = 2,
bty = "n"
)
En este gráfico:
\[ \widehat F_1(t)+\widehat F_2(t)=1-\widehat S(t). \]
Cuando hay un único desenlace, el modelo de Cox permite incorporar covariables mediante
\[ h(t\mid z)=h_0(t)\exp(z^\top\beta). \]
Con riesgos competitivos, una estrategia directa consiste en ajustar un modelo de Cox separado para cada causa. Para la causa \(k\), se considera como evento la causa \(k\), y las demás causas se tratan como censura.
Esta estrategia modela el riesgo específico por causa:
\[ h_k(t\mid z) = h_{0k}(t)\exp(z^\top\beta_k). \]
La interpretación de \(\beta_k\) es sobre la tasa instantánea de ocurrencia de la causa \(k\) entre individuos que siguen libres de cualquier evento justo antes de \(t\).
Cox por causa específica: cáncer de próstata
Se restringe el análisis a pacientes con cáncer de próstata en estadio T2. Se estudia el efecto de:
prostateSurvival.T2 <- subset(prostateSurvival, stage == "T2")
cox_prostate <- coxph(
Surv(survTime, status.prost) ~ grade + ageGroup,
data = prostateSurvival.T2
)
summary(cox_prostate)Call:
coxph(formula = Surv(survTime, status.prost) ~ grade + ageGroup,
data = prostateSurvival.T2)
n= 5920, number of events= 410
coef exp(coef) se(coef) z Pr(>|z|)
gradepoor 1.2199 3.3867 0.1004 12.154 < 2e-16 ***
ageGroup70-74 -0.2860 0.7513 0.2595 -1.102 0.2704
ageGroup75-79 0.4027 1.4958 0.2257 1.784 0.0744 .
ageGroup80+ 0.9728 2.6454 0.2148 4.529 5.92e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
gradepoor 3.3867 0.2953 2.7819 4.123
ageGroup70-74 0.7513 1.3310 0.4518 1.249
ageGroup75-79 1.4958 0.6685 0.9611 2.328
ageGroup80+ 2.6454 0.3780 1.7364 4.030
Concordance= 0.74 (se = 0.012 )
Likelihood ratio test= 252.4 on 4 df, p=<2e-16
Wald test = 243.6 on 4 df, p=<2e-16
Score (logrank) test = 278.9 on 4 df, p=<2e-16Para muerte por cáncer de próstata, los resultados reportados son aproximadamente:
| Covariable | Coeficiente | Razón de riesgos |
|---|---|---|
gradepoor |
1.2199 | 3.3867 |
ageGroup70-74 |
-0.2860 | 0.7513 |
ageGroup75-79 |
0.4027 | 1.4958 |
ageGroup80+ |
0.9728 | 2.6454 |
Los pacientes con tumor pobremente diferenciado tienen peor pronóstico que los pacientes con tumor moderadamente diferenciado. La razón de riesgos estimada es
\[ \exp(1.2199)=3.3867. \]
Es decir, entre pacientes libres de evento justo antes de \(t\), aquellos con tumor pobremente diferenciado tienen una tasa instantánea de muerte por cáncer de próstata aproximadamente 3.39 veces la de los pacientes con tumor moderado, ajustando por edad.
Cox por causa específica: muerte por otras causas
Ahora se toma como evento de interés la muerte por otras causas:
cox_other <- coxph(
Surv(survTime, status.other) ~ grade + ageGroup,
data = prostateSurvival.T2
)
summary(cox_other)Call:
coxph(formula = Surv(survTime, status.other) ~ grade + ageGroup,
data = prostateSurvival.T2)
n= 5920, number of events= 1345
coef exp(coef) se(coef) z Pr(>|z|)
gradepoor 0.28104 1.32451 0.05875 4.784 1.72e-06 ***
ageGroup70-74 0.09462 1.09924 0.12492 0.757 0.44879
ageGroup75-79 0.31330 1.36793 0.11709 2.676 0.00746 **
ageGroup80+ 0.79012 2.20367 0.11204 7.052 1.76e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
gradepoor 1.325 0.7550 1.1805 1.486
ageGroup70-74 1.099 0.9097 0.8605 1.404
ageGroup75-79 1.368 0.7310 1.0874 1.721
ageGroup80+ 2.204 0.4538 1.7692 2.745
Concordance= 0.583 (se = 0.008 )
Likelihood ratio test= 159.6 on 4 df, p=<2e-16
Wald test = 159.3 on 4 df, p=<2e-16
Score (logrank) test = 164.6 on 4 df, p=<2e-16Los resultados reportados son aproximadamente:
| Covariable | Coeficiente | Razón de riesgos |
|---|---|---|
gradepoor |
0.2810 | 1.3245 |
ageGroup70-74 |
0.0946 | 1.0992 |
ageGroup75-79 |
0.3133 | 1.3679 |
ageGroup80+ |
0.7901 | 2.2037 |
Tomados literalmente, estos resultados indican que los pacientes con tumor pobremente diferenciado también presentan mayor riesgo de morir por causas distintas al cáncer de próstata. Sin embargo, este efecto es mucho menor que el observado para muerte por cáncer de próstata.
Esta interpretación debe hacerse con cuidado. Los modelos de riesgos específicos por causa modelan tasas instantáneas entre personas que siguen libres de cualquier evento. No modelan directamente probabilidades acumuladas de morir por cada causa.
Fine y Gray propusieron un enfoque alternativo para modelar covariables en presencia de riesgos competitivos. En lugar de modelar el riesgo específico por causa, modelan el riesgo de subdistribución.
Para la causa \(k\), el riesgo de subdistribución se define como
\[ \bar h_k(t) = \lim_{\Delta t\to 0} \frac{ \Pr(t<T_k<t+\Delta t \mid E) }{ \Delta t }, \]
donde el conjunto condicionante es
\[ E= \left\{ T_k>t \right\} \cup \left\{ T_{k'}\leq t,\ k'\neq k \right\}. \]
La diferencia con el riesgo ordinario es que el conjunto de riesgo incluye no solo a quienes siguen libres de evento, sino también a quienes ya tuvieron antes un evento competitivo distinto de \(k\).
Equivalentemente,
\[ \bar h_k(t) = -\frac{d}{dt} \log{(1-F_k(t))}. \]
El modelo de Fine y Gray es
\[ \bar h_k(t\mid z) = \bar h_{0k}(t)\exp(z^\top\beta). \]
Este modelo se relaciona directamente con la incidencia acumulada \(F_k(t)\). Por tanto, sus coeficientes describen efectos sobre la subdistribución de la causa \(k\), no sobre el riesgo instantáneo específico por causa.
Fine y Gray en R
El modelo de Fine y Gray está implementado en la función crr() del paquete cmprsk.
Primero se crea una matriz de covariables. Se elimina la columna de intercepto porque crr() no la requiere.
library(cmprsk)
X <- model.matrix(
~ grade + ageGroup,
data = prostateSurvival.T2
)
X <- X[, -1, drop = FALSE]
head(X) gradepoor ageGroup70-74 ageGroup75-79 ageGroup80+
4 0 1 0 0
6 1 0 1 0
10 1 0 1 0
13 1 0 0 1
15 0 0 1 0
18 0 0 1 0Modelo para muerte por cáncer de próstata:
fg_prostate <- crr(
ftime = prostateSurvival.T2$survTime,
fstatus = prostateSurvival.T2$status,
cov1 = X,
failcode = 1,
cencode = 0
)
summary(fg_prostate)Competing Risks Regression
Call:
crr(ftime = prostateSurvival.T2$survTime, fstatus = prostateSurvival.T2$status,
cov1 = X, failcode = 1, cencode = 0)
coef exp(coef) se(coef) z p-value
gradepoor 1.132 3.102 0.101 11.20 0.00000
ageGroup70-74 -0.272 0.762 0.253 -1.08 0.28000
ageGroup75-79 0.367 1.443 0.219 1.67 0.09400
ageGroup80+ 0.799 2.224 0.208 3.85 0.00012
exp(coef) exp(-coef) 2.5% 97.5%
gradepoor 3.102 0.322 2.544 3.78
ageGroup70-74 0.762 1.312 0.464 1.25
ageGroup75-79 1.443 0.693 0.939 2.22
ageGroup80+ 2.224 0.450 1.481 3.34
Num. cases = 5920
Pseudo Log-likelihood = -3165
Pseudo likelihood ratio test = 204 on 4 df,El argumento failcode = 1 indica que la causa de interés es muerte por cáncer de próstata.
Resultados reportados para Fine y Gray:
| Covariable | Coeficiente | Razón de subriesgos |
|---|---|---|
gradepoor |
1.132 | 3.102 |
ageGroup70-74 |
-0.272 | 0.762 |
ageGroup75-79 |
0.367 | 1.443 |
ageGroup80+ |
0.799 | 2.224 |
Para muerte por otras causas:
fg_other <- crr(
ftime = prostateSurvival.T2$survTime,
fstatus = prostateSurvival.T2$status,
cov1 = X,
failcode = 2,
cencode = 0
)
summary(fg_other)Competing Risks Regression
Call:
crr(ftime = prostateSurvival.T2$survTime, fstatus = prostateSurvival.T2$status,
cov1 = X, failcode = 2, cencode = 0)
coef exp(coef) se(coef) z p-value
gradepoor 0.126 1.13 0.0584 2.154 3.1e-02
ageGroup70-74 0.103 1.11 0.1252 0.824 4.1e-01
ageGroup75-79 0.273 1.31 0.1176 2.323 2.0e-02
ageGroup80+ 0.667 1.95 0.1128 5.917 3.3e-09
exp(coef) exp(-coef) 2.5% 97.5%
gradepoor 1.13 0.882 1.011 1.27
ageGroup70-74 1.11 0.902 0.867 1.42
ageGroup75-79 1.31 0.761 1.044 1.65
ageGroup80+ 1.95 0.513 1.562 2.43
Num. cases = 5920
Pseudo Log-likelihood = -10350
Pseudo likelihood ratio test = 96.8 on 4 df,Resultados reportados:
| Covariable | Coeficiente | Razón de subriesgos |
|---|---|---|
gradepoor |
0.126 | 1.13 |
ageGroup70-74 |
0.103 | 1.11 |
ageGroup75-79 |
0.273 | 1.31 |
ageGroup80+ |
0.667 | 1.95 |
Las razones de riesgos de Cox por causa específica y las razones de subriesgos de Fine-Gray no tienen la misma interpretación.
Una ventaja del enfoque de riesgos específicos por causa es que permite comparar de forma directa si una covariable tiene efectos distintos sobre causas distintas.
Para ello, se transforma el conjunto de datos a formato largo: cada paciente tendrá una fila por cada causa posible.
En el caso de dos causas, cada paciente tendrá dos filas:
Esto puede hacerse con el paquete mstate.
library(mstate)
tmat <- trans.comprisk(
2,
names = c("event-free", "prostate", "other")
)
tmat to
from event-free prostate other
event-free NA 1 2
prostate NA NA NA
other NA NA NALa matriz de transición indica que un paciente puede pasar del estado libre de evento a:
Ahora se prepara el conjunto en formato largo:
prostate_long <- msprep(
data = prostateSurvival.T2,
time = c(NA, "survTime", "survTime"),
status = c(NA, "status.prost", "status.other"),
keep = c("grade", "ageGroup"),
trans = tmat
)
head(prostate_long)An object of class 'msdata'
Data:
id from to trans Tstart Tstop time status grade ageGroup
1 1 1 2 1 0 27 27 0 mode 70-74
2 1 1 3 2 0 27 27 0 mode 70-74
3 2 1 2 1 0 38 38 0 poor 75-79
4 2 1 3 2 0 38 38 1 poor 75-79
5 3 1 2 1 0 13 13 0 poor 75-79
6 3 1 3 2 0 13 13 0 poor 75-79El conjunto resultante contiene el doble de filas que el conjunto original. Sus columnas principales son:
id: identificador del paciente original;
trans: causa o transición;
trans = 1: muerte por cáncer de próstata;trans = 2: muerte por otras causas;time: tiempo de seguimiento;
status: indicador del evento correspondiente;
grade, ageGroup: covariables copiadas del conjunto original.
Para resumir los eventos:
events(prostate_long)$Frequencies
to
from event-free prostate other no event total entering
event-free 0 410 1345 4165 5920
prostate 0 0 0 410 410
other 0 0 0 1345 1345
$Proportions
to
from event-free prostate other no event
event-free 0.00000000 0.06925676 0.22719595 0.70354730
prostate 0.00000000 0.00000000 0.00000000 1.00000000
other 0.00000000 0.00000000 0.00000000 1.00000000En el texto se reporta:
Los modelos separados por causa se pueden reproducir usando subset:
cox_trans1 <- coxph(
Surv(time, status) ~ grade + ageGroup,
data = prostate_long,
subset = trans == 1
)
cox_trans2 <- coxph(
Surv(time, status) ~ grade + ageGroup,
data = prostate_long,
subset = trans == 2
)
summary(cox_trans1)Call:
coxph(formula = Surv(time, status) ~ grade + ageGroup, data = prostate_long,
subset = trans == 1)
n= 5920, number of events= 410
coef exp(coef) se(coef) z Pr(>|z|)
gradepoor 1.2199 3.3867 0.1004 12.154 < 2e-16 ***
ageGroup70-74 -0.2860 0.7513 0.2595 -1.102 0.2704
ageGroup75-79 0.4027 1.4958 0.2257 1.784 0.0744 .
ageGroup80+ 0.9728 2.6454 0.2148 4.529 5.92e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
gradepoor 3.3867 0.2953 2.7819 4.123
ageGroup70-74 0.7513 1.3310 0.4518 1.249
ageGroup75-79 1.4958 0.6685 0.9611 2.328
ageGroup80+ 2.6454 0.3780 1.7364 4.030
Concordance= 0.74 (se = 0.012 )
Likelihood ratio test= 252.4 on 4 df, p=<2e-16
Wald test = 243.6 on 4 df, p=<2e-16
Score (logrank) test = 278.9 on 4 df, p=<2e-16summary(cox_trans2)Call:
coxph(formula = Surv(time, status) ~ grade + ageGroup, data = prostate_long,
subset = trans == 2)
n= 5920, number of events= 1345
coef exp(coef) se(coef) z Pr(>|z|)
gradepoor 0.28104 1.32451 0.05875 4.784 1.72e-06 ***
ageGroup70-74 0.09462 1.09924 0.12492 0.757 0.44879
ageGroup75-79 0.31330 1.36793 0.11709 2.676 0.00746 **
ageGroup80+ 0.79012 2.20367 0.11204 7.052 1.76e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
gradepoor 1.325 0.7550 1.1805 1.486
ageGroup70-74 1.099 0.9097 0.8605 1.404
ageGroup75-79 1.368 0.7310 1.0874 1.721
ageGroup80+ 2.204 0.4538 1.7692 2.745
Concordance= 0.583 (se = 0.008 )
Likelihood ratio test= 159.6 on 4 df, p=<2e-16
Wald test = 159.3 on 4 df, p=<2e-16
Score (logrank) test = 164.6 on 4 df, p=<2e-16Modelo estratificado por causa
Si se estratifica por causa de muerte usando strata(trans), se permite que cada causa tenga su propia función de riesgo basal:
\[ h_k(t\mid z) = h_{0k}(t)\exp(z^\top\beta). \]
Aquí se asume que las covariables tienen el mismo efecto sobre ambas causas.
cox_common <- coxph(
Surv(time, status) ~ grade + ageGroup + strata(trans),
data = prostate_long
)
summary(cox_common)Call:
coxph(formula = Surv(time, status) ~ grade + ageGroup + strata(trans),
data = prostate_long)
n= 11840, number of events= 1755
coef exp(coef) se(coef) z Pr(>|z|)
gradepoor 0.51477 1.67325 0.04963 10.372 < 2e-16 ***
ageGroup70-74 0.02669 1.02705 0.11228 0.238 0.81210
ageGroup75-79 0.33229 1.39416 0.10392 3.198 0.00139 **
ageGroup80+ 0.83347 2.30129 0.09927 8.396 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
gradepoor 1.673 0.5976 1.5182 1.844
ageGroup70-74 1.027 0.9737 0.8242 1.280
ageGroup75-79 1.394 0.7173 1.1373 1.709
ageGroup80+ 2.301 0.4345 1.8944 2.796
Concordance= 0.621 (se = 0.007 )
Likelihood ratio test= 332.2 on 4 df, p=<2e-16
Wald test = 333.7 on 4 df, p=<2e-16
Score (logrank) test = 348.2 on 4 df, p=<2e-16Resultados reportados:
| Covariable | Coeficiente | Razón de riesgos |
|---|---|---|
gradepoor |
0.515 | 1.673 |
ageGroup70-74 |
0.027 | 1.027 |
ageGroup75-79 |
0.332 | 1.394 |
ageGroup80+ |
0.833 | 2.301 |
Este modelo impone que el efecto de grade y ageGroup es común para ambas causas, aunque las funciones basales pueden diferir.
Interacciones entre covariables y causa
Para evaluar si el efecto de una covariable difiere por causa, se incluyen interacciones con trans.
Por ejemplo, para comparar el efecto del grado tumoral entre muerte por próstata y muerte por otras causas:
cox_grade_interaction <- coxph(
Surv(time, status) ~ grade * factor(trans) + ageGroup + strata(trans),
data = prostate_long
)
summary(cox_grade_interaction)Call:
coxph(formula = Surv(time, status) ~ grade * factor(trans) +
ageGroup + strata(trans), data = prostate_long)
n= 11840, number of events= 1755
coef exp(coef) se(coef) z Pr(>|z|)
gradepoor 1.23878 3.45140 0.09997 12.391 < 2e-16 ***
factor(trans)2 NA NA 0.00000 NA NA
ageGroup70-74 0.02637 1.02672 0.11228 0.235 0.81431
ageGroup75-79 0.33266 1.39468 0.10392 3.201 0.00137 **
ageGroup80+ 0.83328 2.30086 0.09927 8.394 < 2e-16 ***
gradepoor:factor(trans)2 -0.96352 0.38155 0.11571 -8.327 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
gradepoor 3.4514 0.2897 2.8373 4.1985
factor(trans)2 NA NA NA NA
ageGroup70-74 1.0267 0.9740 0.8239 1.2795
ageGroup75-79 1.3947 0.7170 1.1377 1.7097
ageGroup80+ 2.3009 0.4346 1.8940 2.7951
gradepoor:factor(trans)2 0.3815 2.6209 0.3041 0.4787
Concordance= 0.62 (se = 0.007 )
Likelihood ratio test= 402.4 on 5 df, p=<2e-16
Wald test = 400.1 on 5 df, p=<2e-16
Score (logrank) test = 436.5 on 5 df, p=<2e-16Los resultados reportados incluyen:
\[ \hat\beta_{\text{gradepoor}}=1.239, \]
que representa el efecto de gradepoor sobre muerte por cáncer de próstata, la categoría de referencia para trans.
Además, la interacción con trans = 2 es aproximadamente
\[ \hat\beta_{\text{gradepoor:trans2}}=-0.963. \]
Entonces, el efecto de gradepoor sobre muerte por otras causas es
\[ 1.239-0.963=0.276. \]
En escala de razón de riesgos:
\[ \exp(0.276)\approx 1.32. \]
La razón entre el efecto sobre muerte por otras causas y el efecto sobre muerte por cáncer de próstata es
\[ \exp(-0.963)=0.382. \]
Esto indica que el efecto del grado pobre sobre muerte por otras causas es mucho menor que su efecto sobre muerte por cáncer de próstata.
Comparación del efecto de edad entre causas
También puede evaluarse si el efecto de la edad difiere entre causas:
cox_age_interaction <- coxph(
Surv(time, status) ~ (grade + ageGroup) * factor(trans) + strata(trans),
data = prostate_long
)
summary(cox_age_interaction)Call:
coxph(formula = Surv(time, status) ~ (grade + ageGroup) * factor(trans) +
strata(trans), data = prostate_long)
n= 11840, number of events= 1755
coef exp(coef) se(coef) z Pr(>|z|)
gradepoor 1.21985 3.38669 0.10037 12.154 < 2e-16 ***
ageGroup70-74 -0.28596 0.75129 0.25947 -1.102 0.2704
ageGroup75-79 0.40266 1.49580 0.22569 1.784 0.0744 .
ageGroup80+ 0.97281 2.64538 0.21479 4.529 5.92e-06 ***
factor(trans)2 NA NA 0.00000 NA NA
gradepoor:factor(trans)2 -0.93881 0.39109 0.11630 -8.072 6.89e-16 ***
ageGroup70-74:factor(trans)2 0.38059 1.46314 0.28797 1.322 0.1863
ageGroup75-79:factor(trans)2 -0.08937 0.91451 0.25426 -0.351 0.7252
ageGroup80+:factor(trans)2 -0.18269 0.83303 0.24225 -0.754 0.4508
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
gradepoor 3.3867 0.2953 2.7819 4.1230
ageGroup70-74 0.7513 1.3310 0.4518 1.2493
ageGroup75-79 1.4958 0.6685 0.9611 2.3280
ageGroup80+ 2.6454 0.3780 1.7364 4.0301
factor(trans)2 NA NA NA NA
gradepoor:factor(trans)2 0.3911 2.5569 0.3114 0.4912
ageGroup70-74:factor(trans)2 1.4631 0.6835 0.8321 2.5728
ageGroup75-79:factor(trans)2 0.9145 1.0935 0.5556 1.5053
ageGroup80+:factor(trans)2 0.8330 1.2004 0.5181 1.3393
Concordance= 0.621 (se = 0.007 )
Likelihood ratio test= 412 on 8 df, p=<2e-16
Wald test = 402.9 on 8 df, p=<2e-16
Score (logrank) test = 443.5 on 8 df, p=<2e-16Resultados reportados para las interacciones de edad con causa:
| Interacción | Coeficiente | Razón de riesgos | Valor-p |
|---|---|---|---|
ageGroup70-74:trans2 |
0.380 | 1.463 | 0.186 |
ageGroup75-79:trans2 |
-0.089 | 0.914 | 0.725 |
ageGroup80+:trans2 |
-0.183 | 0.833 | 0.451 |
Ninguna de estas interacciones es estadísticamente significativa. Por tanto, no hay evidencia clara de que el efecto de la edad difiera entre muerte por cáncer de próstata y muerte por otras causas, después de ajustar por grado tumoral.
En cambio, sí hay evidencia fuerte de que el efecto del grado tumoral difiere entre causas.