# How do you estimate multiple poisson rates simultaneously?

I have data showing the number of cause-specific deaths in a cohort. I want to use Poisson regression to estimate the cause-specific rate by age. The complication is that I want the predicted cause-specific rates to add-up to the predicted all-cause rate.

Here is some data in R format, showing the total duration of follow-up in the cohort by age (person_years), the number of all-cause deaths (all_cause), and the number of deaths due to three causes: accidents (accidents), non-communicable diseases such as heart disease (ncd), and other (other):

d <- structure(list(age = 18:90, person_years = c(11, 213, 174, 155,
53, 151, 75, 121, 198, 274, 205, 148, 286, 178, 256, 199, 345,
223, 242, 319, 301, 401, 350, 272, 303, 457, 467, 387, 364, 470,
346, 426, 424, 388, 377, 415, 540, 378, 433, 391, 410, 435, 314,
441, 399, 286, 300, 334, 320, 371, 381, 288, 288, 375, 236, 190,
316, 212, 306, 250, 321, 290, 309, 178, 137, 225, 113, 163, 237,
175, 204, 111, 23), all_cause = c(0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 1L, 0L, 2L, 0L, 1L, 0L, 1L, 0L, 1L, 5L, 1L, 1L, 2L,
3L, 2L, 4L, 5L, 7L, 6L, 2L, 7L, 4L, 7L, 9L, 9L, 10L, 17L, 16L,
17L, 10L, 17L, 17L, 17L, 30L, 20L, 17L, 22L, 17L, 25L, 27L, 26L,
29L, 34L, 41L, 29L, 30L, 60L, 36L, 58L, 56L, 94L, 50L, 90L, 58L,
45L, 88L, 36L, 59L, 113L, 75L, 139L, 78L, 21L), accidents = c(0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L,
1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), ncd = c(0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 2L, 0L, 1L, 0L, 1L, 0L, 1L, 5L, 0L,
1L, 1L, 3L, 0L, 3L, 4L, 5L, 6L, 1L, 4L, 1L, 3L, 6L, 4L, 4L, 13L,
15L, 16L, 10L, 12L, 13L, 12L, 26L, 17L, 11L, 20L, 10L, 17L, 19L,
20L, 21L, 28L, 26L, 25L, 23L, 44L, 29L, 45L, 43L, 72L, 31L, 73L,
50L, 36L, 69L, 29L, 36L, 94L, 56L, 118L, 69L, 17L), other = c(0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 3L, 3L, 4L,
2L, 4L, 6L, 4L, 1L, 1L, 0L, 5L, 4L, 4L, 4L, 3L, 5L, 2L, 6L, 8L,
8L, 6L, 8L, 6L, 15L, 4L, 7L, 16L, 7L, 13L, 13L, 22L, 18L, 17L,
8L, 9L, 19L, 7L, 23L, 19L, 19L, 21L, 9L, 4L)), class = "data.frame", row.names = c(NA,
-73L))

We can use Poisson regression to fit a curve on the all-cause rate. In this model, I've used age and age squared:

model1 <- glm(all_cause ~ poly(age, 2) + offset(log(person_years)), data = d, family = 'poisson')

We can predict the age-specific mortality rate and compare this to the observed rates graphically:

model1 <- glm(all_cause ~ poly(age, 2) + offset(log(person_years)), data = d, family = 'poisson')
newdata <- data.frame(age = d$$age, person_years = 100000) d$$predicted_all_cause <- predict(model1, type = 'response', newdata = newdata)
d$$observed_all_cause <- d$$all_cause / d$$person_years * 100000 plot(d$$age, d$$observed_all_cause, ylab = 'mortality rate per 100,000', xlab = 'age') lines(d$$age, d$predicted_all_cause) So that works fine, but if we did something similar with each cause of death, the predicted rates would not add up to the all-cause rate we just estimated. If we fit a separate model for each cause of death ... m_accidents <- glm(accidents ~ poly(age, 2) + offset(log(person_years)), data = d, family = 'poisson') m_ncd <- glm(ncd ~ poly(age, 2) + offset(log(person_years)), data = d, family = 'poisson') m_other <- glm(other ~ poly(age, 2) + offset(log(person_years)), data = d, family = 'poisson') ... Then predict the age-specific rates ... d$$predicted_accidents <- predict(m_accidents, type = 'response', newdata = newdata) d$$predicted_ncd <- predict(m_ncd, type = 'response', newdata = newdata) d$predicted_other <- predict(m_other, type = 'response', newdata = newdata)

... When we add them up, they don't equal the all-cause rate. The head of the table comparing the two is shown below.

d$$sum_predicted_cause_specific <- d$$predicted_accidents + d$$predicted_ncd + d$$predicted_other
head(d[, c('age', 'predicted_all_cause', 'sum_predicted_cause_specific')])

age predicted_all_cause sum_predicted_cause_specific
1  18            102.2936                     108.3501
2  19            111.8447                     118.4173
3  20            122.2923                     129.3918
4  21            133.7212                     141.3514
5  22            146.2239                     154.3805
6  23            159.9019                     168.5714

The fact they are different is not surprising since the curves are estimated separately. But I want them to add up exactly, so that the modelled all-cause rate 'decomposes' into the modelled cause-specific rates. The two approaches I've considered are:

1. a multivariate poisson model, for example using vgam::vglm, but this model will not force the sum of cause-specific rates to add up to the all-cause rate.
2. estimating the all-cause rate using glm and then using multinomial logistic regression to estimate the proportion at each age that are due to each cause of death. This seems like it might work, but feels complicated.

I can't figure out the best approach. I'm happy to use any common statistical software.

• This looks like a case to use multinomial regression to me - you could consider using the functions nnet::multinom or mclogit::mblogit in R. Both are compatible also with emmeans and marginaleffects in case you would like to calculate marginal means and confidence intervals. Commented May 1, 2023 at 22:21
• Do you mean option 2 above? Would you first calculate the proportion of deaths for each cause, then use the proportions as the dependent variable in multinomial regression?
– Dan
Commented May 2, 2023 at 5:18
• Yes - factor cause of death would be your outcome variable and as observation weights you have to pass the count for each. I think you could also pass the log of person years as an offset. mclogit::mblogit also supports estimating overdispersion, or allows adding random effects. You can also use the poisson trick to reparameterize a multinomial model as a poisson one, but it's not so efficient for large datasets. Commented May 2, 2023 at 10:02