Calculer les taux d'incidence à l'aide du modèle de poisson: relation avec le rapport de risque du modèle Cox PH

Je veux calculer les taux d'incidence à présenter le long des ratios de risque afin de présenter des mesures de risque relatives et absolues. J'ai vu dans d'autres études que de tels taux d'incidence peuvent être calculés en utilisant des modèles de poisson avec un temps de suivi dans le modèle comme décalage. J'ai donc essayé cela dans R comme suit:

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

Je m'attendrais à ce que les ratios de taux d'incidence soient similaires aux ratios de risque du modèle Cox PH avec les mêmes termes, mais ils diffèrent en quelque sorte. Suis-je en train d'utiliser la bonne approche pour calculer les taux d'incidence?

Toute aide serait grandement appréciée!

Pour autant que je puisse voir, il n'y a rien de mal avec votre code ou vos calculs. Vous pouvez cependant sauter quelques lignes de code en obtenant les ratios de taux d'incidence par . Les deux modèles font des hypothèses différentes, ce qui conduit potentiellement à des résultats différents.${\tt exp(coef(mod))}$

La régression de Poisson suppose des risques constants. Le modèle de Cox suppose seulement que les dangers sont proportionnels. Si l'hypothèse de dangers constants est remplie, cette question

La régression de Cox a-t-elle une distribution de Poisson sous-jacente?

explique le lien entre Cox et la régression de Poisson.

Nous pouvons utiliser la simulation pour étudier deux situations: les dangers constants et les dangers non constants (mais proportionnels). Simulons d'abord les données d'une population à risque constant. Le rapport de risque a la forme

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

où est un vecteur de paramètres, est un vecteur de covariables et est un nombre positif fixe. Ainsi, l'hypothèse du risque constant de la régression de Poisson est remplie. Maintenant, nous simulons à partir de ce modèle en utilisant le fait (trouvé dans de nombreux livres, par exemple Modeling Survival Data de Therneau, p.13) que la fonction de distribution, , du temps de survie avec comme aléa peut être trouvée comme$\beta$$x$$\lambda_0$$F$$\lambda$

$F(t) = 1 - \exp\left(\int_0^t \lambda(s)\text{ d}s\right)$ .

Avec cela, nous pouvons également trouver l'inverse de , . Avec cette fonction, nous simulons les temps de survie avec l'aléa correct en dessinant des variables uniformes sur et en les transformant en utilisant . Faisons le.$F$$F^{-1}$$(0,1)$$F^{-1}$

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)

Pour le modèle de Cox, les moyennes des estimations des paramètres sont

        sex       rxLev   rxLev+5FU         age 
-0.03826301  0.04167353  0.09069553  0.10025534 

et pour le modèle de Poisson

(Intercept)         sex       rxLev   rxLev+5FU         age 
-1.23651275 -0.03822161  0.03678366  0.08606452  0.09812454 

Pour les deux modèles, nous voyons que cela est proche des vraies valeurs, en se souvenant que la différence entre les hommes et les femmes était de -0,04, par exemple, et elle est estimée à -0,038 pour les deux modèles. Nous pouvons maintenant faire de même avec la fonction de hasard non constante

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

Nous simulons maintenant comme avant.

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)

Pour le modèle Cox, nous obtenons maintenant

        sex       rxLev   rxLev+5FU         age 
-0.04220381  0.04497241  0.09163522  0.10029121  

et pour le modèle de Poisson

(Intercept)         sex       rxLev   rxLev+5FU         age 
-0.12001361 -0.01937333  0.02028097  0.04318946  0.04908300

Dans cette simulation, les moyennes du modèle de Poisson sont nettement plus éloignées des vraies valeurs que celles du modèle de Cox. Cela n'est pas surprenant car nous avons violé l'hypothèse de dangers constants.

Lorsque le danger est constant, la fonction de survivant, , est de la forme$S$

$S(t) = \exp(-\alpha * t)$ ,

pour certains positifs dépendant du sujet spécifique, donc est convexe. Si nous utilisons l'estimateur de Kaplan-Meier pour obtenir une estimation de pour les données d'origine, nous voyons ce qui suit.$\alpha$$S$$S$

Cette fonction semble concave. Cela ne prouve rien, mais cela pourrait donner à penser que l'hypothèse de dangers constants n'est pas remplie pour cet ensemble de données, ce qui pourrait à son tour expliquer les écarts entre les deux modèles.

Une dernière remarque sur les données, pour autant que je sache contient des données sur le temps de récidive du cancer et le temps de mort (il y a deux observations pour chaque valeur de ). Dans ce qui précède, nous l'avons modélisé comme si c'était la même chose. Ce n'est probablement pas une bonne idée.${\tt colon}$${\tt id}$