Modellierung der Sterblichkeitsraten mithilfe der Poisson-Regression

Ich untersuche Trends (zwischen 1998 und 2011) der Sterblichkeitsraten bei Patienten mit Morbus Crohn. Jeder Patient (Fall) wurde in den Jahren 1998 bis 2011 eingeschlossen. Bei der Aufnahme wurde jeder Patient einer gesunden Kontrolle mit demselben Alter und Geschlecht zugeordnet. Ich analysiere Trends bei den Sterblichkeitsraten. Wenn ich dies direkt und ohne Anpassungen mache, erhalte ich im Laufe der Zeit schwankende Sterblichkeitsraten, was vermutlich darauf zurückzuführen ist, dass Personen, die ein bestimmtes Jahr eingeschlossen haben, nicht mit denen vergleichbar sind, die ein anderes Jahr enthalten. Daher möchte ich die Sterblichkeitsraten anpassen. Ich gehe davon aus, dass die Sterblichkeitsraten in beiden Gruppen (Fälle und Kontrollen) im Laufe der Zeit sinken und sich die Kluft zwischen Fällen und Kontrollen sukzessive verringern wird.

Meine Idee ist es, die Anpassung mittels Poisson-Regression vorzunehmen. Meine Daten sind auf individueller Ebene. Ich möchte eine Schätzung der Inzidenzrate (pro 1000 Personenjahre) für Fälle und Kontrollen jedes Jahr von 1998 bis 2011 erhalten. Die Überlebenszeit würde als Offset in das Modell aufgenommen. Etwas Ähnliches ist geschehen hier .

Ich habe die 200 ersten Zeilen aus meinem Datensatz angehängt, der aus 1500 Personen besteht. Hier sind die Daten . Variable Erklärung:

• tot = ob der Patient während der Nachsorge gestorben ist oder nicht
• Überleben = Überlebenszeit in Tagen
• Altersgruppe = kategorisierte Altersgruppe (4 Gruppen)
• Geschlecht = männlich / weiblich
• Diagnose = 0 für gesunde Kontrolle, 1 für Morbus Crohn
• Alter = Alter in Jahren
• Einschlussjahr = Jahr der Aufnahme in die Studie

Was habe ich bisher versucht? Ich habe versucht, Poisson-Modelle mit der Funktion glm () in R unter Verwendung einzelner Beobachtungen (log (Surv) als Offset) anzupassen, aber ich habe entweder einen Fehler erhalten oder konnte nicht herausfinden, wie die Anpassungen verwendet werden sollen. Ich habe die Daten auch in Gruppen zusammengefasst und dann die Todeszahlen in glm () analysiert. Wenn ich die Anpassung verwendete, um Inzidenzraten zu erhalten, konnte ich nur Raten für ein bestimmtes Alter / eine bestimmte Altersgruppe und ein bestimmtes Geschlecht erhalten (wie in der Funktion "dict ()" angegeben).

Ich würde mich sehr über statistische Ratschläge und Codierungsbeispiele freuen, die für den angehängten Datensatz durchgeführt werden können.

Ohne den Datensatz zu sehen (nicht verfügbar), scheint er größtenteils korrekt zu sein. Das Schöne an Poisson-Regressionen ist, dass sie bei Verwendung wie vorgeschlagen Raten liefern können. Beachten Sie möglicherweise, dass es zu einer Überdispersion kommen kann, bei der Sie zu einer negativen binomialen Regression wechseln sollten (siehe MASS-Paket).

Der Poisson-Regression ist es egal, ob die Daten aggregiert sind oder nicht, aber in der Praxis sind nicht aggregierte Daten schwach und können einige unerwartete Fehler verursachen. Beachten Sie, dass Sie surv == 0für keinen der Fälle haben können. Wenn ich getestet habe, sind die Schätzungen gleich:

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

Wenn Sie eine Rate erhalten, ist es wichtig, die Variablen so zu zentrieren, dass der Achsenabschnitt interpretierbar ist, z.

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

Kann als Basisrate interpretiert werden und dann sind die Kovariaten Ratenverhältnisse. Wenn wir den Basiszinssatz nach 10 Jahren wollen:

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

Ich habe derzeit das Einschlussjahr als Trendvariable modelliert, aber Sie sollten wahrscheinlich nach Nichtlinearitäten suchen, und manchmal ist es nützlich, eine Kategorisierung der Zeitpunkte vorzunehmen. Ich habe diesen Ansatz in diesem Artikel verwendet:

D. Gordon, P. Gillgren, S. Eloranta, H. Olsson, M. Gordon, J. Hansson und KE Smedby -basierte Studie ”Melanoma Res., vol. 25, nein. 4, S. 348–356, August 2015.