En muchas aplicaciones de análisis de supervivencia, además de comparar dos curvas marginales de supervivencia, interesa estudiar cómo influyen variables exógenas (covariables) sobre el riesgo de falla. En el capítulo anterior vimos el test de rango logarítmico para comparar dos grupos sin asumir explícitamente una forma paramétrica para cada función de supervivencia. Para ello introdujimos la alternativa de Lehmann, definida mediante
\[ S_{1}(t)=\bigl[S_{0}(t)\bigr]^{\psi}, \]
y observamos que esta relación entre supervivencias se traduce, a través de la ecuación
\[ S(t)=\exp\bigl(-\!\!\int_{0}^{t}h(u)\,du\bigr), \]
en el siguiente supuesto sobre tasas de riesgo (hazard):
\[ \begin{equation*} h_{1}(t)=\psi\,h_{0}(t)\,. \tag{5.1.1} \end{equation*} \]
Esta última ecuación se conoce como la suposición de riesgos proporcionales entre los dos grupos. Bajo ella, la razón de riesgos \(\psi\) cuantifica de manera uniforme la diferencia entre sus funciones de riesgo a lo largo del tiempo: si \(\psi>1\), el grupo 1 tiene un riesgo mayor que el grupo 0 en todo instante, y si \(\psi<1\) ocurre lo contrario.
Para generalizar este planteamiento e incorporar información de covariables, se escribe la constante de proporcionalidad \(\psi\) en función de un vector de covariables \(z=(z_{1},\dots,z_{p})\). La forma más utilizada es la parametrización exponencial:
\[ \begin{equation*} \psi = e^{\,z\beta}, \tag{5.1.2} \end{equation*} \]
donde \(\beta=(\beta_{1},\dots,\beta_{p})^{T}\) es el vector de coeficientes que mide el efecto de cada covariable sobre el riesgo. Bajo esta formulación, la covariable \(j\) con valor \(z_{j}\) multiplica el riesgo basal \(h_{0}(t)\) por el factor \(e^{z_{j}\beta_{j}}\), manteniendo la proporción constante a lo largo del tiempo.
En conjunto, el modelo de riesgos proporcionales se expresa como
\[ h(t\mid z) \;=\; h_{0}(t)\,\exp\bigl(z\beta\bigr), \]
donde:
Este modelo, propuesto originalmente por Cox, permite estimar los coeficientes \(\beta\) sin especificar una forma paramétrica para \(h_{0}(t)\), utilizando la denominada verosimilitud parcial. De este modo, se pueden incorporar variables explicativas al análisis de supervivencia con datos censurados, extendiendo el alcance de las regresiones clásicas (lineal o logística) a escenarios temporales donde no se asume distribución paramétrica para la supervivencia basal.
Para ajustar un modelo de riesgos proporcionales sin especificar la forma de la función de riesgo basal, se recurre a la verosimilitud parcial, que sólo involucra los tiempos de falla (censuras no aportan factores). A cada tiempo de falla \(t_{j}\) (ordenados de menor a mayor) le corresponde un conjunto de individuos “en riesgo” justo antes de \(t_{j}\), denotado por \(R_{j}\). Bajo el modelo de riesgos proporcionales, la tasa de riesgo del sujeto \(i\) en \(t_{j}\) se escribe como
\[ h_{i}(t_{j}) \;=\; h_{0}(t_{j})\,\psi_{i}, \quad\text{siendo}\quad \psi_{i} = e^{z_i\beta}. \]
En el caso de comparar dos grupos (control versus tratamiento), basta con tomar \(z_{i}=0\) para controles y \(z_{i}=1\) para tratados, de modo que \(\psi_{i}=1\) en el grupo control y \(\psi_{i}=\psi=e^{\beta}\) en el grupo experimental.
Supongamos que, en el primer tiempo de falla \(t_{1}\), el paciente \(i\in R_{1}\) es el que falla. La probabilidad condicional de que justo ese individuo falle antes que los otros en riesgo es
\[ p_{1} \;=\; \frac{h_{i}(t_{1})}{\sum_{k\in R_{1}} h_{k}(t_{1})} \;=\; \frac{h_{0}(t_{1})\,\psi_{i}}{\sum_{k\in R_{1}} h_{0}(t_{1})\,\psi_{k}} \;=\; \frac{\psi_{i}}{\sum_{k\in R_{1}} \psi_{k}}, \]
donde el factor \(h_{0}(t_{1})\) se anula en numerador y denominador. Tras ese evento, todos los sujetos que fallaron o fueron censurados antes del siguiente fallo forman el nuevo riesgo residual \(R_{2}\), y se repite el mismo razonamiento para obtener
\[ p_{2} \;=\;\frac{\psi_{i_2}}{\sum_{k\in R_{2}}\psi_{k}}, \quad \dots,\quad p_{D} = \frac{\psi_{i_D}}{\sum_{k\in R_{D}}\psi_{k}}, \]
donde \(D\) es el número total de tiempos de falla.
La verosimilitud parcial se define como el producto de estos factores condicionales (uno por cada falla):
\[ L(\psi)\;=\; \prod_{j=1}^{D} p_{j} \;=\; \prod_{j=1}^{D} \frac{\psi_{i_j}}{\sum_{k\in R_{j}} \psi_{k}}. \]
Obsérvese que los tiempos exactos no intervienen, sino sólo el orden de las fallas y la pertenencia de cada sujeto al grupo control ( \(\psi_{i}=1\) ) o tratamiento ( \(\psi_{i}=\psi\) ). Además, la función de riesgo basal \(h_{0}(t)\) queda completamente fuera de esta expresión.
En un ejemplo del capítulo anterior teníamos seis pacientes: tres en control (\(\psi_{i}=1\)) y tres en tratamiento (\(\psi_{i}=\psi\)). Los tiempos de falla ordenados fueron 6, 10, 15 y 25, sin empates. El riesgo inicial es \(R_{1}={1,2,3,4,5,6}\), con \(\psi_{1}=\psi_{2}=\psi_{4}=1\) y \(\psi_{3}=\psi_{5}=\psi_{6}=\psi\). En cada etapa:
En \(t_{1}=6\) falla el paciente 1 (control), por lo que
\[ p_{1} \;=\; \frac{\psi_{1}}{\sum_{k\in R_{1}}\psi_{k}} \;=\;\frac{1}{3\,\psi + 3}. \]
Tras eso, el paciente 1 sale de riesgo y, además, el paciente 2 (control) queda censurado en \(t=7\). Queda entonces \(R_{2}={3,4,5,6}\), con dos controles (\(\psi=1\)) y dos tratados (\(\psi\)). En \(t_{2}=10\) falla el paciente 3 (tratado), así que
\[ p_{2} \;=\; \frac{\psi_{3}}{\sum_{k\in R_{2}}\psi_{k}} \;=\;\frac{\psi}{3\,\psi + 1}. \]
Quedan en riesgo \(R_{3}={4,5,6}\); en \(t_{3}=15\) falla paciente 4 (control), entonces
\[ p_{3} \;=\; \frac{1}{2\,\psi + 1}. \]
Finalmente, en \(t_{4}=25\) sólo queda en riesgo el paciente 6 (tratado), que falla con probabilidad 1, de modo que \(p_{4}=1\).
En consecuencia, la verosimilitud parcial (en función de \(\psi\)) es
\[ L(\psi) \;=\; \frac{1}{3\psi + 3}\;\times\;\frac{\psi}{3\psi+1}\;\times\;\frac{1}{2\psi + 1}\;\times\;1 \;=\; \frac{\psi}{(3\psi+3)\,(3\psi+1)\,(2\psi+1)}. \]
Para facilitar la estimación, se suele parametrizar \(\psi=e^{\beta}\) y trabajar con la log–verosimilitud parcial:
\[ \ell(\beta) \;=\;\log L(e^{\beta}) \;=\;\beta \;-\; \log(3e^{\beta}+3)\;-\;\log(3e^{\beta}+1)\;-\;\log(2e^{\beta}+1). \]
El estimador de máxima verosimilitud parcial ( MPLE ) \(\hat\beta\) es el valor que maximiza \(\ell(\beta)\). Nótese que, aunque ésta es una “verosimilitud”, no es una función de probabilidad en el sentido clásico, pues no contiene factores para tiempos censurados; sin embargo, todas las derivaciones estadísticas (puntuales, errores estándar, tests) se obtienen de forma análoga a una verosimilitud ordinaria.
Para hallar \(\hat\beta\) usamos optim(), indicando que maximizemos la función. Por ejemplo:
plsimple <- function(beta){
psi <- exp(beta)
result <- (beta-log(3*psi+3)-log(3*psi +1)-log(2*psi +1))
return(result)
}
# Se busca el máximo de plsimple en el intervalo [-3, 1]
res <- optim(par=0, fn=plsimple, method="L-BFGS-B",
lower=-3, upper=1,
control=list(fnscale=-1))
res$par # valor de beta que maximiza la log-verosimilitud parcial[1] -1.326129En este ejemplo el resultado fue \(\hat\beta\approx -1.326\), de donde \(\hat\psi = e^{-1.326}\approx 0.265\) mide la razón de riesgos estimada entre el grupo experimental y el control.
Al derivar \(\ell(\beta)\) respecto a \(\beta\) y evaluarlo en \(\beta=0\), se obtiene el estadístico score correspondiente. En efecto, esa derivada coincide con la suma de \((d_{0j}-e_{0j})\) a lo largo de los tiempos de falla y, por tanto, equivale al estadístico clásico del test de rango logarítmico. Esto explica por qué el test de rango log—rank surge como prueba score en este modelo de Cox: para \(H_{0}:\beta=0\), el score \(U(\beta)|_{\beta=0}\) es asintóticamente normal con varianza dada por la información de Fisher. Así, sin conocer \(h_{0}(t)\), podemos contrastar \(H_{0}:\beta=0\) frente a \(H_{A}:\beta\neq0\) usando el estadístico score o la razón de verosimilitudes parciales.
En el contexto del modelo de riesgos proporcionales, se pueden derivar analogías de las tres pruebas habituales en teoría de la verosimilitud—Wald, score y razón de verosimilitudes—empleando la verosimilitud parcial. Aunque la teoría subyacente es más compleja que en el caso de verosimilitud completa, en la práctica estas pruebas suelen arrojar resultados muy similares.
Denotemos por
\[ l(\beta) \]
la log–verosimilitud parcial, por
\[ S(\beta)=l'(\beta) \]
la función score, y por
\[ I(\beta) = -l''(\beta) = -S'(\beta) \]
la información observada, que mide la curvatura de la log–verosimilitud en cada valor de \(\beta\). El valor \(\hat\beta\) que maximiza \(l(\beta)\) es el estimador de máxima verosimilitud parcial (MPLE).
La prueba de Wald se basa en la aproximación normal de \(\hat\beta\). Su estadístico es
\[ Z_{W} \;=\;\frac{\hat\beta}{\mathrm{s.e.}(\hat\beta)}, \]
donde la varianza de \(\hat\beta\) se estima como
\[ \operatorname{Var}(\hat\beta)\approx \bigl[I(\hat\beta)\bigr]^{-1}, \qquad \mathrm{s.e.}(\hat\beta) = \frac{1}{\sqrt{I(\hat\beta)}}. \]
Así, para contrastar \(H_{0}:\beta=0\) frente a una alternativa bilateral, se rechaza \(H_{0}\) si \(\lvert Z_{W}\rvert > z_{\alpha/2}\). Equivalentemente,
\[ Z_{W}^{2}=\frac{\hat\beta^{2}}{1/I(\hat\beta)} \sim \chi^{2}_{1}\quad (\text{bajo } H_{0}). \]
Un intervalo de confianza \(1-\alpha\) para \(\beta\) es
\[ \hat\beta \pm z_{\alpha/2}\,\mathrm{s.e.}(\hat\beta). \]
La prueba score emplea la función score evaluada bajo \(H_{0}\) (normalmente \(\beta=0\)). Denotemos
\[ U(\beta) = S(\beta) = l'(\beta), \quad I(\beta) = -l''(\beta). \]
Bajo \(H_{0}:\beta=0\), el estadístico es
\[ Z_{S} \;=\; \frac{S(0)}{\sqrt{\,I(0)\,}}, \]
y se rechaza \(H_{0}\) si \(\lvert Z_{S}\rvert > z_{\alpha/2}\), o equivalentemente
\[ Z_{S}^{2} = \frac{[S(0)]^{2}}{I(0)} \sim \chi^{2}_{1}. \]
En particular, en este modelo de Cox, la prueba score coincide con el test de rango logarítmico.
La prueba de razón de verosimilitudes compara los valores de la log–verosimilitud parcial bajo la hipótesis nula \(\beta=0\) y bajo la estimación \(\hat\beta\). El estadístico es
\[ \Lambda \;=\; 2\bigl[l(\hat\beta)\;-\;l(0)\bigr] \;\sim\;\chi^{2}_{1}\quad(\text{aprox. bajo }H_{0}). \]
Una ventaja de este test es su invarianza ante transformaciones monótonas de \(\beta\) (por ejemplo, usar \(\psi=e^{\beta}\)).
En la sección anterior se obtuvo la log–verosimilitud parcial
\[ \ell(\beta) \;=\; \beta - \log(3e^{\beta}+3)\;-\;\log(3e^{\beta}+1)\;-\;\log(2e^{\beta}+1), \]
correspondiente a los datos sintéticos de seis pacientes (tres controles y tres tratados). Usando la función coxph() de R se obtienen los siguientes resultados:
library(survival)
tt <- c(6, 7, 10, 15, 19, 25)
delta <- c(1, 0, 1, 1, 0, 1)
trt <- c(0, 0, 1, 0, 1, 1)
result.cox <- coxph(Surv(tt, delta) ~ trt)
summary(result.cox)Call:
coxph(formula = Surv(tt, delta) ~ trt)
n= 6, number of events= 4
coef exp(coef) se(coef) z Pr(>|z|)
trt -1.3261 0.2655 1.2509 -1.06 0.289
exp(coef) exp(-coef) lower .95 upper .95
trt 0.2655 3.766 0.02287 3.082
Concordance= 0.7 (se = 0.116 )
Likelihood ratio test= 1.21 on 1 df, p=0.3
Wald test = 1.12 on 1 df, p=0.3
Score (logrank) test = 1.27 on 1 df, p=0.3A continuación se muestran los cálculos numéricos en R de cada prueba, usando funciones para derivar la log–verosimilitud parcial definida como plsimple:
library(numDeriv)
## 1) Prueba score (Score) en beta = 0
score0 <- grad(func=plsimple, x=0) # S(0)
info0 <- -hessian(func=plsimple, x=0) # I(0), la información
chisq_S <- (score0^2) / info0
pval_S <- pchisq(chisq_S, df=1, lower.tail=FALSE)
## 2) Prueba de Wald
betahat <- result.cox$coefficients # ̂β obtenido por optim o coxph
info_b <- -hessian(func=plsimple, x=betahat)
se_b <- 1 / sqrt(info_b) # desviación estándar de ̂β
Z_wald <- betahat / se_b
pval_W <- 2 * pnorm(abs(Z_wald), lower.tail=FALSE)
## 3) Prueba de razón de verosimilitudes
loglik_hat <- plsimple(betahat)
loglik_0 <- plsimple(0)
chisq_LR <- 2 * (loglik_hat - loglik_0)
pval_LR <- pchisq(chisq_LR, df=1, lower.tail=FALSE)Cuando deseamos incluir varios predictores en el modelo de riesgos proporcionales, extendemos lo visto en la sección anterior definiendo para cada individuo \(i\):
La razón de riesgos de \(i\) respecto al riesgo “base” se escribe \[ \psi_{i} \;=\; \exp\bigl(z_{i}'\,\beta\bigr), \] de modo que el riesgo instantáneo del sujeto \(i\) en tiempo \(t\) es \[ h_{i}(t) \;=\; h_{0}(t)\,\psi_{i}, \] donde \(h_{0}(t)\) es el riesgo basal no parametrizado.
Sea \(\{t_{1} < t_{2} < \cdots < t_{D}\}\) el conjunto de los \(D\) tiempos de falla observados. Para cada tiempo de falla \(t_{j}\), definimos el conjunto en riesgo \[ R_{j} \;=\; \{\,\text{todos los individuos aún “vivos” justo antes de }t_{j}\}. \] De ese conjunto, precisamente uno falla en \(t_{j}\).
La verosimilitud parcial se construye producto de razones condicionales: en cada tiempo de falla \(t_{j}\), la probabilidad de que el sujeto \(i\in R_{j}\) sea quien falle es \[ \frac{h_{i}(t_{j})}{\sum_{k\in R_{j}} h_{k}(t_{j})} \;=\; \frac{h_{0}(t_{j})\,\psi_{i}}{\sum_{k\in R_{j}} h_{0}(t_{j})\,\psi_{k}} \;=\; \frac{\psi_{i}}{\sum_{k\in R_{j}} \psi_{k}}, \] pues el factor \(h_{0}(t_{j})\) se cancela en numerador y denominador.
Si etiquetamos con \(\delta_{j}=1\) al individuo que efectivamente falla en \(t_{j}\), la verosimilitud parcial total resulta \[ L(\beta) \;=\; \prod_{j=1}^{D} \frac{\psi_{j}}{\sum_{k\in R_{j}} \psi_{k}}, \] donde \(\psi_{j} = \exp\bigl(z_{j}'\beta\bigr)\) corresponde al individuo que falla en \(t_{j}\).
Su logaritmo es \[ \ell(\beta) \;=\; \sum_{j=1}^{D} \biggl[ \ln\bigl(\psi_{j}\bigr) \;-\; \ln\bigl(\!\sum_{k\in R_{j}} \psi_{k}\bigr) \biggr] \;=\; \sum_{j=1}^{D} z_{j}'\beta \;-\; \sum_{j=1}^{D} \ln\Bigl(\sum_{k\in R_{j}} e^{\,z_{k}'\beta}\Bigr). \]
Observación.
Para derivar estimadores y pruebas, necesitamos la función score
\[
S(\beta)
\;=\;
\frac{\partial \ell(\beta)}{\partial \beta}
\;\in\;\mathbb{R}^{p},
\] cuyas componentes son, para \(l=1,\dots,p\), \[
S_{l}(\beta)
\;=\;
\frac{\partial \ell(\beta)}{\partial \beta_{l}}
\;=\;
\sum_{j=1}^{D}
\Biggl[
z_{j\,l}
\;-\;
\frac{\sum_{k\in R_{j}}
z_{k\,l}\,\exp(z_{k}'\beta)}
{\sum_{k\in R_{j}} \exp(z_{k}'\beta)}
\Biggr].
\] Intuitivamente, en cada tiempo de falla \(t_{j}\) se resta el valor observado \(z_{j\,l}\) menos el “promedio ponderado” de \(z_{\,\cdot\,l}\) en el riesgo \(R_{j}\). Cuando sólo hay una covariable binaria, \(S(\beta=0)\) coincide con el estadístico de rango logarítmico.
La matriz de información observada (o Hessiano negativo) es
\[
I(\beta)
\;=\;
-\,\frac{\partial^{2}\ell(\beta)}{\partial \beta\,\partial \beta'}
\;=\;
-\,\frac{\partial S(\beta)}{\partial \beta'}
\;\in\;\mathbb{R}^{p\times p}.
\] Cada entrada \([I(\beta)]_{l,m} = -\,\frac{\partial^{2}\ell}{\partial\beta_{l}\,\partial\beta_{m}}\). En la práctica, \(I(\beta)\) se evalúa numéricamente o a través del software, y su inversa aproxima la matriz de varianzas y covarianzas de \(\hat\beta\).
Sea \(\hat\beta\) el valor que maximiza \(\ell(\beta)\). Para contrastar
\[
H_{0}: \beta = 0
\quad\text{versus}\quad
H_{A}: \beta \neq 0,
\] podemos usar las tres formas clásicas:
Prueba de Wald
\[
X_{W}^{2}
\;=\;
\hat{\beta}' \,I(\hat{\beta})\,\hat{\beta}
\;\sim\;\chi^{2}_{p}\quad(\text{aprox. bajo }H_{0}).
\] Equivalente a formar \(Z_{l} = \hat\beta_{l}/\mathrm{s.e.}(\hat\beta_{l})\) y comparar contra \(\chi^{2}\).
Prueba score (Score)
\[
X_{S}^{2}
\;=\;
S(0)'\;I(0)^{-1}\;S(0)
\;\sim\;\chi^{2}_{p}\quad(\text{aprox. bajo }H_{0}).
\] Evalúa la pendiente de \(\ell(\beta)\) en \(\beta=0\) versus su varianza.
Prueba de razón de verosimilitudes
\[
X_{LR}^{2}
\;=\;
2\,\bigl[\ell(\hat\beta)\;-\;\ell(0)\bigr]
\;\sim\;\chi^{2}_{p}\quad(\text{aprox. bajo }H_{0}).
\] Es invariante a reparametrizaciones (por ejemplo, si se usa \(\psi=\exp(\beta)\)).
En todos los casos, si \(p>1\), el estadístico se compara contra \(\chi^{2}_{p}\).
Nota práctica.
La mayoría de paquetes (por ejemplo, coxph() en R) reportan directamente \(\hat\beta\), sus errores estándar, y los valores de \(\chi^{2}\) para las tres pruebas, junto con el valor–\(p\).
Tras ajustar un modelo de riesgos proporcionales, queremos reconstruir la supervivencia “base” (cuando todos los covariables valen cero). Primero, la tasa de riesgo basal en cada tiempo de falla \(t_{i}\) se estima como
\[
\hat{h}_{0}(t_{i})
\;=\;
\frac{d_{i}}
{\sum_{j\in R_{i}} \exp\bigl(z_{j}'\,\hat\beta\bigr)},
\]
donde:
Si sólo hay un grupo (\(\beta=0\)), esto coincide con el estimador de Nelson–Aalen:
\[
\hat{h}_{0}(t_{i}) = \frac{d_{i}}{|R_{i}|}.
\]
Para obtener la supervivencia basal acumulada, definimos
\[
\hat{H}_{0}(t)
\;=\;
\sum_{t_{j}\le t} \hat{h}_{0}(t_{j}),
\qquad
\hat{S}_{0}(t)
\;=\;
\exp\bigl[-\hat{H}_{0}(t)\bigr].
\]
Esta \(\hat{S}_{0}(t)\) corresponde a la supervivencia de un individuo con \(z=0\). Si algún covariable no puede significar “cero” (ej., edad), su interpretación es limitada.
Para un individuo con vector de covariables \(z\), la supervivencia condicional es
\[
\hat{S}(t \mid z)
\;=\;
\bigl[\hat{S}_{0}(t)\bigr]^{\exp\bigl(z'\,\hat\beta\bigr)}.
\]
En R, la función basehaz() extrae la hazard acumulada basal; conviene usar centered=FALSE para que calcule \(H_{0}(t)\) en \(\beta=0\) (el valor por defecto centra en la media de los covariados, lo cual no siempre es apropiado para variables categóricas).
En la práctica, pueden aparecer “ties” (empates) en los tiempos de falla por dos razones distintas:
Tiempo continuo registrado con redondeo.
Si los tiempos verdaderos son continuos pero se registran en días, semanas o meses enteros, varios individuos pueden compartir el mismo valor registrado aunque en realidad sus tiempos difieran ligeramente.
Tiempo verdaderamente discreto.
Cuando el evento sólo puede ocurrir en intervalos discretos (p. ej., número de ciclos o de compresiones hasta fallo).
En ambos casos, el tratamiento de los empates en el ajuste del modelo de riesgos proporcionales difiere ligeramente. A continuación se describen los dos enfoques principales, y luego se mencionan dos aproximaciones prácticas para datos con muchos empates.
Si suponemos que los tiempos son contínuos pero están redondeados, podemos aplicar el modelo usual
\[
h(t;z) \;=\; e^{z\,\beta}\,h_{0}(t),
\]
y escribir la verosimilitud parcial como producto de factores “marginales” para cada tiempo de falla distinto. Cuando \(r\) sujetos fallan simultáneamente en \(t_i\),
Por ejemplo, si en \(t=1\) dos pacientes del grupo tratado fallan de entre 10 en riesgo (4 tratados y 6 control):
Este método se llama marginal o Breslow exacto, y coincide con la solución “exacta” para el caso continuo redondeado.
Si los tiempos son verdaderamente discretos, conviene formular un modelo de riesgos discretos tipo logístico:
\[
\frac{h(t;z)}{1 - h(t;z)} \;=\; e^{z\,\beta}\,\frac{h_{0}(t)}{1 - h_{0}(t)}.
\]
Cuando \(r\) sujetos fallan simultáneamente en \(t_i\) y hay \(\binom{n_i}{r}\) formas de elegirlos:
Por ejemplo, si en \(t=1\) dos tratados fallan de entre 10 en riesgo:
Este enfoque se denomina exacto discreto (o “Breslow exacto discreto”). Con pocas observaciones se puede enumerar, pero si hay muchos empates se vuelve computacionalmente costoso.
Cuando hay demasiados empates para enumerar, conviene usar aproximaciones que eviten la suma sobre permutaciones completas:
Aproximación de Breslow.
Se simplifica el denominador para que sólo dependa del total en riesgo \(n_i\), elevándolo a la potencia \(r\). En el ejemplo con dos empates en \(t=1\):
\[
L_{1}^{\text{Breslow}}(\beta)
\;=\;
\frac{\binom{2}{2}\,e^{2\beta}}{\bigl(4\,e^{\beta} + 6\bigr)^{2}}
\;=\;
\frac{2\,e^{2\beta}}{\bigl(4\,e^{\beta} + 6\bigr)^{2}},
\]
y análogamente para otros tiempos con \(r=2\).
Aproximación de Efron.
Refina la anterior corrigiendo parcialmente el denominador para reflejar cuántos sujetos “ya han fallado” dentro del mismo tiempo. Por ejemplo, para dos empates en \(t=1\):
\[
L_{1}^{\text{Efron}}(\beta)
\;=\;
\frac{e^{\beta}}{4\,e^{\beta}+6}\;\times\;
\frac{e^{\beta}}{\bigl(4\,e^{\beta}+6\bigr) \;-\; \tfrac{1}{2}\,(2\,e^{\beta})}
\;=\;
\frac{e^{\beta}}{4\,e^{\beta} + 6}\;\times\;
\frac{e^{\beta}}{3\,e^{\beta} + 6},
\]
donde en el segundo factor se restan “a mitad” los dos sujetos que fallarán. De modo similar para otros tiempos con empates.
Tanto Breslow como Efron dan valores muy cercanos al método exacto discreto, pero con notable ahorro computacional cuando hay muchos empates.
En muchos estudios clínicos es relevante medir el tiempo desde el diagnóstico hasta el evento (muerte o censura) en lugar de desde el ingreso al ensayo. Cuando los pacientes se inscriben en el estudio un tiempo después del diagnóstico, los casos con supervivencia muy corta (que fallecerían antes de inscribirse) no se observan. Este sesgo se conoce como truncamiento a la izquierda.
Para ilustrar, consideremos un ensayo hipotético con seis pacientes, donde:
tt es el tiempo desde el ingreso hasta el evento,status indica muerte (1) o censura (0),grp es el indicador de tratamiento (0=control, 1=experimental),backTime es el tiempo retrospectivo hasta el diagnóstico (negativo):tt <- c(6, 7, 10, 15, 19, 25)
status <- c(1, 0, 1, 1, 0, 1)
grp <- c(0, 0, 1, 0, 1, 1)
backTime <- c(-3, -11, -3, -7, -10, -5)Si ignoramos el truncamiento, ajustamos un modelo de Cox estándar:
coxph(Surv(tt, status) ~ grp)Call:
coxph(formula = Surv(tt, status) ~ grp)
coef exp(coef) se(coef) z p
grp -1.3261 0.2655 1.2509 -1.06 0.289
Likelihood ratio test=1.21 on 1 df, p=0.2715
n= 6, number of events= 4 que suele dar un hazard ratio estimado \(\exp(\hat\beta)<1\) (beneficio del tratamiento) pero no significativo (\(p\approx0.29\)).
Para analizar desde el diagnóstico, realineamos los tiempos:
tm.enter <- -backTime
tm.exit <- tt - backTimey volvemos a ajustar contabilizando el truncamiento:
coxph(Surv(tm.enter, tm.exit, status, type="counting") ~ grp)Call:
coxph(formula = Surv(tm.enter, tm.exit, status, type = "counting") ~
grp)
coef exp(coef) se(coef) z p
grp -1.073 0.342 1.235 -0.869 0.385
Likelihood ratio test=0.81 on 1 df, p=0.3677
n= 6, number of events= 4 El resultado confirma de nuevo un efecto no significativo (\(p\approx0.37\)), pero ahora referenciado a la supervivencia desde diagnóstico.
De manera similar, en los datos de Channing House, al comparar supervivencia de hombres y mujeres condicionando a alcanzar cierta edad de ingreso, es imprescindible usar el enfoque de start time en Surv() para evitar estimaciones sesgadas cuando el riesgo inicial varía con la edad:
library(asaur)
data(ChanningHouse)
ChanningHouse <- within(ChanningHouse, {
entryYears <- entry/12
exitYears <- exit/12
})
channing68 <- ChanningHouse[ChanningHouse$exitYears >= 68,]
coxph(Surv(entryYears, exitYears, cens, type="counting") ~ sex, data=channing68)Call:
coxph(formula = Surv(entryYears, exitYears, cens, type = "counting") ~
sex, data = channing68)
coef exp(coef) se(coef) z p
sexMale 0.2733 1.3143 0.1762 1.552 0.121
Likelihood ratio test=2.3 on 1 df, p=0.1292
n= 451, number of events= 172 En análisis de supervivencia, además del tratamiento o intervención principal, a menudo disponemos de múltiples variables clínicas, demográficas o biomarcadores que pueden influir en el tiempo hasta el evento. En ensayos aleatorizados, la aleatorización debe equilibrar los factores de confusión, pero suele ser conveniente incluirlas en el modelo para corregir posibles desbalances y comprender el efecto de estas covariables. En estudios observacionales, el ajuste por covariables es esencial para estimar sin sesgos el efecto de la intervención y explorar el papel de otros factores.
Consideremos los datos simulados donde el genotipo (“mutant”/“wt”) confunde la estimación del efecto del tratamiento. Primero ajustamos un modelo de Cox sin tener en cuenta el genotipo:
set.seed(4321)
lambda.mutant.0 <- 0.03
lambda.mutant.1 <- 0.03 * 0.55
lambda.wt.0 <- 0.03 * 0.2
lambda.wt.1 <- 0.03 * 0.2 * 0.55
tt.control.mutant <- rexp(25, rate = lambda.mutant.0)
tt.treat.mutant <- rexp(125, rate = lambda.mutant.1)
tt.control.wt <- rexp(125, rate = lambda.wt.0)
tt.treat.wt <- rexp( 25, rate = lambda.wt.1)
ttAll <- c(tt.control.mutant, tt.treat.mutant,
tt.control.wt, tt.treat.wt)
status <- rep(1, length(ttAll)) # sin censura
genotype <- factor(c(rep("mutant", 150), rep("wt", 150)))
trt <- factor(c(rep("control", 25), rep("treatment", 125),
rep("control", 125), rep("treatment", 25)))
coxph(Surv(ttAll, status) ~ trt)Call:
coxph(formula = Surv(ttAll, status) ~ trt)
coef exp(coef) se(coef) z p
trttreatment 0.4639 1.5902 0.1173 3.955 7.64e-05
Likelihood ratio test=15.51 on 1 df, p=8.198e-05
n= 300, number of events= 300 Este resultado sugiere (incorrectamente) que el tratamiento aumenta el riesgo (\(e^{0.464}=1.59\)). Al estratificar por genotipo, corregimos esta interpretación:
coxph(Surv(ttAll, status) ~ trt + strata(genotype))Call:
coxph(formula = Surv(ttAll, status) ~ trt + strata(genotype))
coef exp(coef) se(coef) z p
trttreatment -0.4530 0.6357 0.1643 -2.757 0.00584
Likelihood ratio test=7.66 on 1 df, p=0.005655
n= 300, number of events= 300 Ahora \(\hat\beta<0\), indicando que el tratamiento reduce la tasa de fallo dentro de cada genotipo. Finalmente, podemos incluir el genotipo como covariable:
coxph(Surv(ttAll, status) ~ trt + genotype)Call:
coxph(formula = Surv(ttAll, status) ~ trt + genotype)
coef exp(coef) se(coef) z p
trttreatment -0.4523 0.6361 0.1634 -2.768 0.00563
genotypewt -1.5676 0.2085 0.1825 -8.589 < 2e-16
Likelihood ratio test=93.39 on 2 df, p=< 2.2e-16
n= 300, number of events= 300 Aquí, además del efecto corregido del tratamiento, cuantificamos que el genotipo “wt” tiene un riesgo menor que el “mutant”.
En modelos de riesgos proporcionales, además de variables indicadoras para comparar dos grupos, es habitual incorporar covariables adicionales que reflejen características del paciente, como edad, género o raza. Una covariable continua (por ejemplo, edad) cuantifica cambios en el log–hazard por unidad de la variable, mientras que una covariable categórica de \(L\) niveles se introduce mediante \(L-1\) variables dummy comparadas con un nivel de referencia.
Sea \(z_{i}=(z_{i1},\dots,z_{ik})^\top\) el vector de \(k\) covariables del sujeto \(i\), donde algunos \(z_{ij}\) son indicadores (0/1) y otros numéricos. El modelo de riesgos proporcionales se escribe
\[ \log\bigl(\psi_{i}\bigr)=z_{i}^{\prime}\beta, \]
de modo que el hazard ratio relativo a la línea base es \(\psi_{i}=e^{z_{i}^{\prime}\beta}\). Entonces:
Es posible mejorar el modelo incluyendo:
En este ejemplo generamos un conjunto de datos de supervivencia con dos covariables —raza (tres niveles: “white”, “black” y “other”) y edad (continua)— para ilustrar cómo se ajusta un modelo de riesgos proporcionales de Cox y cómo se recuperan los valores verdaderos de los efectos.
Primero creamos 60 edades uniformes entre 40 y 80 años, y asignamos por igual 20 observaciones a cada nivel de raza, tomando “white” como categoría de referencia. A continuación, definimos el tasa de fallo de cada sujeto mediante un modelo lineal en la escala logarítmica:
\[ \log(\lambda_{i}) \;=\; -4.5 \;+\; \begin{cases} 0 & \text{si raza = white},\\ 1 & \text{si raza = black},\\ 2 & \text{si raza = other}, \end{cases} \;+\;0.05 \times \text{edad}_{i}, \]
de modo que la edad incrementa la tasa de fallo en un 5 % por año y las razas “black” y “other” elevan la tasa en un factor \(e^1\) y \(e^2\), respectivamente. Simulamos entonces tiempos de supervivencia exponenciales sin censura y ajustamos el modelo de Cox:
set.seed(123)
age <- runif(n = 60, min = 40, max = 80)
race <- factor(rep(c("white","black","other"), each = 20))
race <- relevel(race, ref = "white")
X <- model.matrix(~ race + age)[,-1] #Ilustración de la matroz de diseño
log.rate <- -4.5 + c(rep(0,20), rep(1,20), rep(2,20)) + 0.05 * age
tt <- rexp(n = 60, rate = exp(log.rate))
status <- rep(1, 60)
library(survival)
result.cox <- coxph(Surv(tt, status) ~ race + age)
summary(result.cox)Call:
coxph(formula = Surv(tt, status) ~ race + age)
n= 60, number of events= 60
coef exp(coef) se(coef) z Pr(>|z|)
raceblack 1.37357 3.94941 0.37288 3.684 0.00023 ***
raceother 2.55178 12.82995 0.45654 5.589 2.28e-08 ***
age 0.07686 1.07989 0.01555 4.943 7.70e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
raceblack 3.949 0.25320 1.902 8.202
raceother 12.830 0.07794 5.243 31.393
age 1.080 0.92602 1.047 1.113
Concordance= 0.759 (se = 0.031 )
Likelihood ratio test= 44.54 on 3 df, p=1e-09
Wald test = 36.16 on 3 df, p=7e-08
Score (logrank) test = 38.87 on 3 df, p=2e-08El ajuste recupera coeficientes muy cercanos a los valores de simulación. Estos resultados muestran que, ajustando por edad, la raza “black” multiplica el riesgo de muerte por ≈ 3.94 y “other” por ≈ 12.82 respecto a “white”, mientras que cada año adicional aumenta el hazard en un 8% aproximadamente. Los valores p confirman la significación estadística de todos los efectos.
Cuando ajustamos varios modelos de Cox que difieren en las covariables incluidas, sólo podemos comparar aquellos cuyos conjuntos de predictores sean anidados: es decir, que las covariables de un modelo reducido estén contenidas en el modelo completo.
En el ejemplo de los datos pharmacoSmoking consideramos tres modelos:
ageGroup4employmentageGroup4 + employmentModelos A y B quedan ambos anidados en C, pero A y B entre sí no lo están, por lo que no pueden compararse directamente.
Al ajustar cada modelo con coxph obtenemos estas log-verosimilitudes parciales (logLik):
data("pharmacoSmoking")
modelA.coxph <- coxph(Surv(ttr, relapse) ~ ageGroup4,data = pharmacoSmoking)
modelB.coxph <- coxph(Surv(ttr, relapse) ~ employment,data = pharmacoSmoking)
modelC.coxph <- coxph(Surv(ttr, relapse) ~ ageGroup4+employment,data = pharmacoSmoking)
logLik(modelA.coxph) 'log Lik.' -380.043 (df=3)logLik(modelB.coxph) 'log Lik.' -385.1232 (df=2)logLik(modelC.coxph) 'log Lik.' -377.7597 (df=5)Para evaluar si, por ejemplo, ageGroup4 aporta información adicional a un modelo que ya contiene employment, comparamos A (reducido) contra C (completo). El estadístico de razón de verosimilitudes es
\[ \Lambda \;=\; 2\,\bigl[\ell(\hat\beta_{\text{C}})\;-\;\ell(\hat\beta_{\text{A}})\bigr] \;=\; 2\,\bigl[-377.7597 \;-\;(-380.043)\bigr] \;=\;4.567 \]
Bajo \(H_{0}\) (coeficientes de ageGroup4 = 0) este estadístico sigue una \(\chi^{2}\) con \(5-3=2\) grados de libertad, de donde
pchisq(4.567, df = 2, lower.tail = FALSE) [1] 0.1019268anova(modelA.coxph,modelC.coxph)Analysis of Deviance Table
Cox model: response is Surv(ttr, relapse)
Model 1: ~ ageGroup4
Model 2: ~ ageGroup4 + employment
loglik Chisq Df Pr(>|Chi|)
1 -380.04
2 -377.76 4.5666 2 0.1019Como \(p>0.05\), no rechazamos \(H_{0}\): ageGroup4 no mejora significativamente el modelo que ya incluye employment.
De igual manera, para probar si employment añade información más allá de ageGroup4, comparamos B contra C:
\[ \Lambda \;=\;2\,\bigl[-377.7597 \;-\;(-385.1232)\bigr]\;=\;14.727, \]
con \(5-2=3\) grados de libertad:
pchisq(14.727, df = 3, lower.tail = FALSE) [1] 0.002065452anova(modelB.coxph,modelC.coxph)Analysis of Deviance Table
Cox model: response is Surv(ttr, relapse)
Model 1: ~ employment
Model 2: ~ ageGroup4 + employment
loglik Chisq Df Pr(>|Chi|)
1 -385.12
2 -377.76 14.727 3 0.002065 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Aquí \(p<0.01\), por lo que sí rechazamos \(H_{0}\): employment aporta significativamente al modelo que incluye ageGroup4.
Finalmente, para saber si ageGroup4 debe entrar “al menos” en el modelo, lo comparamos con el modelo nulo (sin covariables). Si
modelNull.coxph <- coxph(Surv(ttr, relapse) ~ 1,data = pharmacoSmoking)
logLik(modelNull.coxph) 'log Lik.' -386.1533 (df=0)entonces
\[ \Lambda \;=\;2\,\bigl[-380.043\;-\;(-386.1533)\bigr]\;=\;12.221 \]
con \(3-0=3\) grados de libertad,
pchisq(12.221, df = 3, lower.tail = FALSE)[1] 0.006663207anova(modelNull.coxph,modelA.coxph)Analysis of Deviance Table
Cox model: response is Surv(ttr, relapse)
Model 1: ~ 1
Model 2: ~ ageGroup4
loglik Chisq Df Pr(>|Chi|)
1 -386.15
2 -380.04 12.221 3 0.006664 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1y concluimos que ageGroup4 por sí sola es altamente significativa.
Cuando disponemos de muchos predictores potenciales, el procedimiento paso a paso (stepwise) basado en p-valores presenta varios problemas:
Una alternativa más robusta para buscar el modelo “óptimo” es el Akaike Information Criterion (AIC), definido como
\[ \mathrm{AIC} = -2\,\ell(\hat\beta)\;+\;2\,k \]
donde
- \(\ell(\hat\beta)\) es la log-verosimilitud parcial evaluada en el M.P.L.E.,
- \(k\) el número de parámetros.
El primer término mide el ajuste al dato (menor \(-2\,\ell\) = mejor ajuste) y el segundo penaliza la complejidad (más parámetros = mayor penalización). Por tanto, un modelo con AIC más baja es preferible.
pharmacoSmokingHabiendo ajustado los tres modelos anidados de la sección anterior:
AIC(modelA.coxph) [1] 766.086AIC(modelB.coxph) [1] 774.2464AIC(modelC.coxph) [1] 765.5194El mejor según AIC es el Modelo C (ageGroup4 + employment), seguido muy de cerca por A.
En lugar de usar p-valores, podemos usar la función step() con criterio AIC:
modelAll <- coxph(Surv(ttr, relapse) ~ grp + gender + race +
employment + yearsSmoking + levelSmoking +
ageGroup4 + priorAttempts + longestNoSmoke,
data = pharmacoSmoking)
result.step <- step(modelAll,
scope = list(upper = ~ grp + gender + race + employment +
yearsSmoking + levelSmoking + ageGroup4 +
priorAttempts + longestNoSmoke,
lower = ~ grp))Start: AIC=770.2
Surv(ttr, relapse) ~ grp + gender + race + employment + yearsSmoking +
levelSmoking + ageGroup4 + priorAttempts + longestNoSmoke
Df AIC
- race 3 766.98
- yearsSmoking 1 768.20
- gender 1 768.20
- priorAttempts 1 768.24
- levelSmoking 1 768.47
- longestNoSmoke 1 769.04
<none> 770.20
- employment 2 772.45
- ageGroup4 3 774.11
Step: AIC=766.98
Surv(ttr, relapse) ~ grp + gender + employment + yearsSmoking +
levelSmoking + ageGroup4 + priorAttempts + longestNoSmoke
Df AIC
- levelSmoking 1 764.98
- gender 1 765.00
- priorAttempts 1 765.01
- yearsSmoking 1 765.04
- longestNoSmoke 1 766.29
<none> 766.98
- employment 2 768.37
- ageGroup4 3 770.16
+ race 3 770.20
Step: AIC=764.98
Surv(ttr, relapse) ~ grp + gender + employment + yearsSmoking +
ageGroup4 + priorAttempts + longestNoSmoke
Df AIC
- gender 1 763.00
- priorAttempts 1 763.01
- yearsSmoking 1 763.06
- longestNoSmoke 1 764.29
<none> 764.98
- employment 2 766.37
+ levelSmoking 1 766.98
- ageGroup4 3 768.18
+ race 3 768.47
Step: AIC=763
Surv(ttr, relapse) ~ grp + employment + yearsSmoking + ageGroup4 +
priorAttempts + longestNoSmoke
Df AIC
- priorAttempts 1 761.02
- yearsSmoking 1 761.08
- longestNoSmoke 1 762.31
<none> 763.00
- employment 2 764.42
+ gender 1 764.98
+ levelSmoking 1 765.00
- ageGroup4 3 766.32
+ race 3 766.48
Step: AIC=761.02
Surv(ttr, relapse) ~ grp + employment + yearsSmoking + ageGroup4 +
longestNoSmoke
Df AIC
- yearsSmoking 1 759.10
- longestNoSmoke 1 760.34
<none> 761.02
- employment 2 762.42
+ priorAttempts 1 763.00
+ gender 1 763.01
+ levelSmoking 1 763.02
- ageGroup4 3 764.50
+ race 3 764.52
Step: AIC=759.1
Surv(ttr, relapse) ~ grp + employment + ageGroup4 + longestNoSmoke
Df AIC
- longestNoSmoke 1 758.42
<none> 759.10
- employment 2 760.42
+ yearsSmoking 1 761.02
+ gender 1 761.08
+ levelSmoking 1 761.08
+ priorAttempts 1 761.08
+ race 3 762.52
- ageGroup4 3 766.90
Step: AIC=758.42
Surv(ttr, relapse) ~ grp + employment + ageGroup4
Df AIC
<none> 758.42
+ longestNoSmoke 1 759.10
- employment 2 760.31
+ yearsSmoking 1 760.34
+ gender 1 760.39
+ priorAttempts 1 760.40
+ levelSmoking 1 760.41
+ race 3 761.53
- ageGroup4 3 767.24La salida muestra en cada paso qué término, al eliminarse, más reduce la AIC, hasta llegar al modelo final:
Surv(...) ~ grp + employment + ageGroup4
AIC = 758.42Cuya tabla de coeficientes (result.step) confirma la significancia de los predictores retenidos.
Nota: