Por lo que veo no hay nada malo en tu código o en tus cálculos. Sin embargo, podrías saltarte algunas líneas de código obteniendo los ratios de incidencia mediante ${\tt exp(coef(mod))}$ . Los dos modelos parten de supuestos diferentes, lo que puede llevar a resultados distintos.
La regresión de Poisson asume riesgos constantes. El modelo de Cox sólo asume que los riesgos son proporcionales. Si se cumple el supuesto de riesgos constantes esta pregunta
¿Tiene la regresión de Cox una distribución de Poisson subyacente?
explica la conexión entre la regresión de Cox y la de Poisson.
Podemos utilizar la simulación para estudiar dos situaciones: riesgos constantes y riesgos no constantes (pero proporcionales). Primero vamos a simular los datos de una población con un riesgo constante. La razón de riesgo tiene la forma
$\lambda(t) = \lambda_0\exp(\beta'x)$ ,
donde $\beta$ es un vector de parámetros, $x$ es un vector de covariables y $\lambda_0$ es un número fijo positivo. Por tanto, se cumple el supuesto de riesgo constante de la regresión de Poisson. Ahora simulamos a partir de este modelo utilizando el hecho (que se encuentra en muchos libros, por ejemplo, Modeling Survival Data de Therneau, p.13) de que la función de distribución $F$ del tiempo de supervivencia con $\lambda$ como peligro se puede encontrar como
$F(t) = 1 - \exp\left(\int_0^t \lambda(s)\text{ d}s\right)$ .
Con esto también podemos encontrar la inversa de $F$ , $F^{-1}$ . Con esta función simulamos los tiempos de supervivencia con el peligro correcto dibujando variables que son uniformes en $(0,1)$ y transformarlos utilizando $F^{-1}$ . Hagámoslo.
library(survival)
data(colon)
data <- with(colon, data.frame(sex = sex, rx = rx, age = age))
n <- dim(data)[1]
# defining linP, the linear predictor, beta*x in the above notation
linP <- with(colon, log(0.05) + c(0.05, 0.01)[as.factor(sex)] + c(0.01,0.05,0.1)[rx] + 0.1*age)
h <- exp(linP)
simFuncC <- function() {
cens <- runif(n) # simulating censoring times
toe <- -log(runif(n))/h # simulating times of events
event <- ifelse(toe <= cens, 1, 0) # deciding if time of event or censoring is the smallest
data$time <- pmin(toe, cens)
data$event <- event
mCox <- coxph(Surv(time, event) ~ sex + rx + age, data = data)
mPois <- glm(event ~ sex + rx + age, data = data, offset = log(time))
c(coef(mCox), coef(mPois))
}
sim <- t(replicate(1000, simFuncC()))
colMeans(sim)
Para el modelo de Cox, las medias de las estimaciones de los parámetros son
sex rxLev rxLev+5FU age
-0.03826301 0.04167353 0.09069553 0.10025534
y para el modelo de Poisson
(Intercept) sex rxLev rxLev+5FU age
-1.23651275 -0.03822161 0.03678366 0.08606452 0.09812454
Para ambos modelos, vemos que se acerca a los valores reales, recordando que la diferencia entre hombres y mujeres era de -0,04, por ejemplo, y se estima en -0,038 para ambos modelos. Ahora podemos hacer lo mismo con la función de riesgo no constante
$\lambda(t) = \lambda_0 t \exp(\beta'x)$ .
Ahora simulamos como antes.
simFuncN <- function() {
cens <- runif(n)
toe <- sqrt(-log(runif(n))/h)
event <- ifelse(toe <= cens, 1, 0)
data$time <- pmin(toe, cens)
data$event <- event
mCox <- coxph(Surv(time, event) ~ sex + rx + age, data = data)
mPois <- glm(event ~ sex + rx + age, data = data, offset = log(time))
c(coef(mCox), coef(mPois))
}
sim <- t(replicate(1000, simFuncN()))
colMeans(sim)
Para el modelo de Cox obtenemos ahora
sex rxLev rxLev+5FU age
-0.04220381 0.04497241 0.09163522 0.10029121
y para el modelo de Poisson
(Intercept) sex rxLev rxLev+5FU age
-0.12001361 -0.01937333 0.02028097 0.04318946 0.04908300
En esta simulación, las medias del modelo de Poisson están claramente más alejadas de los valores reales que las del modelo de Cox. Esto no es sorprendente, ya que hemos violado el supuesto de riesgos constantes.
Cuando el peligro es constante, la función de supervivencia, $S$ es de la forma
$S(t) = \exp(-\alpha * t) $ ,
para algún positivo $\alpha$ depende del tema específico, por lo que $S$ es convexo. Si utilizamos el estimador de Kaplan-Meier para obtener una estimación de $S$ para los datos originales, vemos lo siguiente.
![enter image description here]()
Esta función parece cóncava. Esto no demuestra nada, pero podría ser un indicio de que el supuesto de riesgos constantes no se cumple para este conjunto de datos, lo que a su vez podría explicar las discrepancias entre los dos modelos.
Una última observación sobre los datos, por lo que sé ${\tt colon}$ contiene datos sobre el tiempo hasta la reaparición del cáncer y el tiempo hasta la muerte (hay dos observaciones para cada valor de ${\tt id}$ ). En lo anterior, lo hemos modelado como si fuera lo mismo. Probablemente no sea una buena idea.
