Calcule las tasas de incidencia utilizando el modelo de Poisson: relación con la razón de riesgo del modelo Cox PH

Quiero calcular las tasas de incidencia para presentar a lo largo de las razones de riesgo para presentar medidas de riesgo relativas y absolutas. Vi en otros estudios que tales tasas de incidencia pueden calcularse utilizando modelos de Poisson con tiempo de seguimiento en el modelo como compensación. Así que intenté eso en R de la siguiente manera:

library(survival)

# Get example data
data(colon)
colon$status <- ifelse(colon$etype==1,0,1) # set to 0/1 (needed for poisson later on)

# Fit cox model for rx (age + sex adjusted)
coxph(Surv(time,status)~rx+sex+age, data=colon)
# HR (rxLev): 0.92  
# HR (rxLev+5FU): 0.74

# Get incidence rates using poisson models with same terms and log(time) as offset
mod <- glm(status~offset(log(time))+rx+sex+age, data=colon, family=poisson)

# Get rates using predict-function
Obs <- predict(mod, data.frame(time=1, rx="Obs", age=mean(colon$age),
                                   sex=mean(colon$sex)),  type="response")
Lev <- predict(mod, data.frame(time=1, rx="Lev", age=mean(colon$age), 
                                   sex=mean(colon$sex)),  type="response")
Lev5FU <- predict(mod, data.frame(time=1, rx="Lev+5FU", age=mean(colon$age), 
                                      sex=mean(colon$sex)),  type="response")

# Calculate incidence rate ratio's:
Lev/Obs # 0.98
Lev5FU/Obs # 0.84

Esperaría que las tasas de incidencia sean similares a las razones de riesgo del modelo Cox PH con los mismos términos, pero de alguna manera difieren. ¿Estoy usando el enfoque correcto para calcular las tasas de incidencia?

Cualquier ayuda sería muy apreciada!

Por lo que puedo ver, no hay nada malo con su código o sus cálculos. Sin embargo, puede omitir algunas líneas de código obteniendo los índices de incidencia por${\tt exp(coef(mod))}$. Los dos modelos hacen suposiciones diferentes, y esto potencialmente conduce a resultados diferentes.

La regresión de Poisson supone riesgos constantes. El modelo de Cox solo supone que los riesgos son proporcionales. Si se cumple el supuesto de riesgos constantes, esta pregunta

¿La regresión de Cox tiene una distribución de Poisson subyacente?

explica la conexión entre la regresión de Cox y Poisson.

Podemos utilizar la simulación para estudiar dos situaciones: peligros constantes y peligros no constantes (pero proporcionales). Primero simulemos datos de una población con un peligro constante. La razón de riesgo tiene la forma

$\lambda(t) = \lambda_0\exp(\beta'x)$,

dónde $\beta$ es un vector de parámetros, $x$ es un vector de covariables y $\lambda_0$Es un número positivo fijo. Por lo tanto, se cumple el supuesto de riesgo constante de la regresión de Poisson. Ahora simulamos a partir de este modelo utilizando el hecho (encontrado en muchos libros, por ejemplo, Modelando datos de supervivencia 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 el inverso de $F$, $F^{-1}$. Con esta función simulamos tiempos de supervivencia con el peligro correcto al dibujar variables que son uniformes en$(0,1)$ y transformándolos usando $F^{-1}$. Vamos a hacerlo.

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, los promedios 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 esto está cerca de los valores verdaderos, recordando que la diferencia entre hombres y mujeres era -0.04, por ejemplo, y se estima en -0.038 para ambos modelos. Ahora podemos hacer lo mismo con la función de peligro 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 Cox ahora tenemos

        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, los promedios del modelo de Poisson están claramente más lejos de los valores reales que los 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)$,

por algo positivo $\alpha$ depende del tema específico, por lo tanto $S$es convexo Si usamos el estimador de Kaplan-Meier para obtener una estimación de$S$ para los datos originales, vemos lo siguiente.

Esta función se ve cóncava. Esto no prueba nada, pero podría ser una pista de que no se cumple el supuesto de riesgos constantes para este conjunto de datos, lo que a su vez podría explicar las discrepancias entre los dos modelos.

Una observación final sobre los datos, que yo sepa ${\tt colon}$ contiene datos sobre el tiempo hasta la recurrencia 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 la misma cosa. Probablemente no sea una buena idea.