Modélisation des taux de mortalité à l'aide de la régression de Poisson

J'examine les tendances (entre 1998 et 2011) des taux de mortalité chez les patients atteints de la maladie de Crohn. Chaque patient (cas) a été inclus entre 1998 et 2011. À l'inclusion, chaque patient a été apparié à un contrôle sain de même âge et de même sexe. J'analyse les tendances des taux de mortalité. En faisant cela directement, sans aucun ajustement, j'obtiens des taux de mortalité fluctuants dans le temps, ce qui est probablement dû au fait que les individus inclus une année donnée ne seront pas comparables à ceux inclus une autre année. Je vise donc à ajuster les taux de mortalité. Je m'attends à ce que les taux de mortalité dans les deux groupes (cas et témoins) diminuent avec le temps et l'écart entre les cas et les témoins se réduise successivement.

Mon idée est de faire l'ajustement par régression de Poisson. Mes données sont au niveau individuel. Je souhaite obtenir une estimation du taux d'incidence (pour 1 000 années-personnes) des cas et des témoins chaque année de 1998 à 2011. Le temps de survie serait inclus comme décalage dans le modèle. Quelque chose de similaire a été fait ici .

J'ai joint les 200 premières lignes de mon ensemble de données, qui comprend 1500 personnes. Voici les données . Explication variable:

• mort = si le patient est décédé ou non pendant le suivi
• surv = durée de survie en jours
• groupe d'âge = groupe d'âge catégorisé (4 groupes)
• sexe = homme / femme
• diagnostic = 0 pour un contrôle sain, 1 pour la maladie de Crohn
• âge = âge en années
• inclusion_year = année d'inclusion dans l'étude

Qu'est-ce que j'ai essayé jusqu'à présent? J'ai essayé d'adapter les modèles de Poisson avec la fonction glm () dans R, en utilisant des observations individuelles (log (surv) comme décalage), mais j'ai reçu une erreur ou je n'ai pas pu comprendre comment utiliser les ajustements. J'ai également agrégé les données en groupes, puis analysé le nombre de décès dans glm (); lorsque j'ai utilisé l'ajustement pour obtenir des taux d'incidence, je ne pouvais obtenir que des taux pour un âge / groupe d'âge et un sexe spécifiques (comme cela devait être spécifié dans la fonction prédire ()).

J'apprécierais vraiment quelques conseils statistiques et des exemples de codage, qui peuvent être faits sur l'ensemble de données joint.

Sans voir l'ensemble de données (non disponible), il semble généralement correct. Ce qui est bien avec les régressions de Poisson, c'est qu'elles peuvent fournir des taux lorsqu'elles sont utilisées comme suggéré. Une chose qui peut être utile de garder à l'esprit est qu'il peut y avoir une surdispersion où vous devez passer à une régression binomiale négative (voir le package MASS).

La régression de Poisson ne se soucie pas de savoir si les données sont agrégées ou non, mais en pratique, les données non agrégées sont fragiles et peuvent provoquer des erreurs inattendues. Notez que vous ne pouvez pas en avoir surv == 0pour aucun des cas. Quand j'ai testé les estimations sont les mêmes:

set.seed(1)
n <- 1500
data <- 
  data.frame(
    dead = sample(0:1, n, replace = TRUE, prob = c(.9, .1)),
    surv = ceiling(exp(runif(100))*365),
    gender = sample(c("Male", "Female"), n, replace = TRUE),
    diagnosis = sample(0:1, n, replace = TRUE),
    age = sample(60:80, n, replace = TRUE),
    inclusion_year = sample(1998:2011, n, replace = TRUE)
  )

library(dplyr)
model <- 
  data %>% 
  group_by(gender, 
           diagnosis,
           age,
           inclusion_year) %>% 
  summarise(Deaths = sum(dead),
            Person_time = sum(surv)) %>%
  glm(Deaths ~ gender + diagnosis + I(age - 70) + I(inclusion_year - 1998) + offset(log(Person_time/10^3/365.25)), 
      data = . , family = poisson)

alt_model <- glm(dead ~ gender + diagnosis + I(age - 70) + I(inclusion_year - 1998) + offset(log(surv/10^3/365.25)), 
    data = data , family = poisson)
sum(coef(alt_model) - coef(model))
# > 1.779132e-14
sum(abs(confint(alt_model) - confint(model)))
# > 6.013114e-11

Lorsque vous obtenez un taux, il est important de centrer les variables afin que l'interception soit interprétable, par exemple:

> exp(coef(model)["(Intercept)"])
(Intercept) 
    51.3771

Peut être interprété comme le taux de base, puis les covariables sont des ratios de taux. Si nous voulons le taux de base après 10 ans:

> exp(coef(model)["(Intercept)"] + coef(model)["I(inclusion_year - 1998)"]*10)
(Intercept) 
     47.427 

J'ai actuellement modélisé l'année d'inclusion comme une variable de tendance, mais vous devriez probablement vérifier les non-linéarités et parfois il est utile de faire une catégorisation des points temporels. J'ai utilisé cette approche dans cet article:

D. Gordon, P. Gillgren, S. Eloranta, H. Olsson, M. Gordon, J. Hansson et KE Smedby, «Tendances temporelles de l'incidence du mélanome cutané selon l'emplacement anatomique détaillé et les modèles d'exposition aux rayonnements ultraviolets: une population rétrospective -based study », Melanoma Res., vol. 25, non. 4, p. 348–356, août 2015.