Independence in Poisson regression when used for rates estimation

Question

I'm interested in using Poisson regression for estimating mortality rates. I've observations for a rather long period, and have to face the fact that covariates often vary over time. This is especially the case when the age of the people comes into play and the observation period exceeds a few months.

In this case, here is my way to deal with this issue : the observation period is divided into intervals where the covariates are constant. Each observation then becomes a set of observations and the exposures are calculated for each interval. However, the log-likelihood in Poisson regression models assume independence between the observations. But, as I have, say for a person named Jack, a period where he is forty and another one where he is forty one, I will have for Jack two observations in my model, and I can't really say that those two observations are independent.

Can I just assume that the death rates are independent anyway, or are there specific methods to deal with this issue in a better way?

EDIT

In this lesson, at the paragrah 7.4.3 The Equivalent Poisson Model there is something that could let think my approach is not invalid. Nevertheless, as far as I'am aggregating data in each bunch of pseudo-observation (and not as in the text , one for each combination of individual and interval) I wonder if it still stand.

EDIT About grouping data : I have a set of thousands of observations. Each one relate to an individual, and for each one I have the date of entry in the study, the date of exit, the date of death, an indicator for smoking or not, the birth date, the gender (and a few other one but that don't matter here). I cut the whole interval of time covered form the first entry to the last at each time something happen (someone entre or leave the study, or a covariate changes). For each interval of time, and for each set of similar age and other covaraiates, I compute the total exposition, and the number of deaths. And I this way I get "pseudo-observations" with a number of death, an exposition, and covariates. Then I fit a Poisson Regression (generalised linear model, with expostion as offset).

Answer 1

The notes you posted explain why creating these pseudo-observations and fitting a Poisson log-linear regression (Poisson GLM with log-link) is equivalent to a piece-wise exponential survival model. They also note that it is valid to aggregate the pseudo-observations with equal covariates. This is mentioned at the very end of the "The Equivalent Poisson Model" section. So I will stick to the question about independence as it relates to these intervals.

Why doesn't creating many small interval observations out of one big interval observation in Poisson regression violate some independence assumption?

The key thing is that by assuming your response is Poisson distributed, you already assume something about the independence of these smaller intervals.

Suppose your data set were simpler, and all observations represented a single unit interval of time, say January - December 2020, and each row represented a different subject. And suppose your response variable were a discrete count of something for each observation within this interval of time. You decide to do some form of Poisson regression, having made the following assumptions about the data:

1. Your response variable is a Poisson distributed random variable. That is, you are counting events which occur independently of one another at a fixed rate in the interval.
2. Your observations are independent of one another. That is, you have pairs $(X_i, y_i)$, and each pair is independent.

You then decide to split Jack's row into two sub-intervals. You partition it into a January - June sub-interval, and a July - December sub-interval, and distribute the count of something between the sub-intervals based on whether each something happened in Jan-Jun, or Jul-Dec. Maybe you are doing this for fun, or maybe it is because there is some covariate you think is important which changes half-way through the year (which, note, implies that you think the rate changes midway through the original interval, and also that it wouldn't be possible to include this covariate in a regression using the original data).

You started with one observation $(X_j, y_j)$ for Jack, now you have two observations $(X_{j1}, y_{j1})$ and $(X_{j2}, y_{j2})$. It sure "feels like" these aren't independent. Indeed, $X_{j1}$ and $X_{j2}$ are definitely not independent, they will be necessarily identical on all the covariates which you could have included before splitting into sub-intervals. But crucially, $y_{j1}$ and $y_{j2}$ are independent, by assumption. It is a property of the Poisson distribution that the occurrences of something we are counting in the original interval are occurring independently from one another. So necessarily, the count of something for Jack in the Jan-Jun sub-interval is independent of the count of something for him in the Jul-Dec sub-interval. Again, you got this for free once you assumed the response was Poisson.

What about the fact that $X_{j1}$ and $X_{j2}$ are not independent? Aren't we pretending we have two observations, when really it is only one? Well, this is part of what the exposure offset is for. When we maximize likelihood in a Poisson regression without exposure offsets (constant duration intervals), the log-likelihood function we're maximizing looks like this:

$$\sum_i{y_i \log{\mu_i} - \mu_i}.$$

But when we have differing interval lengths and exposure offsets, it looks like this:

$$\sum_i{y_i \log{(\mu_ie_i)} - \mu_ie_i} = \sum_i{y_i\log{(\mu_i)} + y_i\log(e_i)} - \mu_ie_i.$$

Which rearranges as: $$\sum_i{e_i\left(\frac{y_i}{e_i}\log{(\mu_i)} - \mu_i\right) + y_i\log(e_i)}$$

The term on the right doesn't depend on $\mu$ and so is irrelevant. So, our likelihood function is equivalent to the likelihood function of a Poisson regression with $y_i/e_i$ as the response, and $e_i$ as the weight. Conceptually we are standardizing the observed rate of occurrence so that all observations have the same measurement unit (observations per unit time), and then we are weighting each observation by the amount of time it represents.

All that is to say, if we partition Jack's interval into two equal length sub-intervals, then each interval is only being treated as half an observation.

Answer 2

If you want to divide the observation period for each subject into epochs you can use generalized estimating equations (GEE) to model this as a quasi-Poisson process while accounting for the correlation between subject measurements across epochs. I don't think this is standard practice, though. It might be more appropriate to adjust a cross-sectional Poisson model for baseline age or other factors to investigate if the mortality rate differs across groups. In either model I recommend including an additional scale parameter or using the empirical sandwich covariance estimator when estimating standard errors.

Answer 3

The answer by @JonnyLomond answers your question as posed (+1), but also raises some questions about how you will do the modeling. You say:

For each interval of time, and for each set of similar age and other covariates, I compute the total exposition, and the number of deaths,

which suggests that you are binning continuous predictors like age, which is not usually a good idea. You generally are better off modeling individuals directly with their actual covariate values, then displaying results as desired for particular combinations of covariate values.

Also, it sounds like you might not be using time-interval-specific piecewise-constant hazards as in those notes, but rather a single constant baseline hazard for all your data (related to the intercept in the Poisson modeling). As the course notes you cite say (Section 7.2.2):

treating the deaths as Poisson conditional on exposure time leads to exactly the same estimates (and standard errors) as treating the exposure times as censored observations from an exponential distribution.

If you have a single intercept shared by all time intervals, you are thus assuming an exponential survival model, which doesn't always work well. Unless there's a compelling reason to force an exponential model onto the data, you might do better with a Cox model, which allows the data to tell you the form of the baseline survival function as well as the hazards associated with the covariates.

Your data seem already to be arranged in a way close to the counting-process format used to fit Cox models with time-varying covariates. From your description, each individual is evidently represented by multiple rows covering non-overlapping time intervals. During each time interval, the individual's covariate values, included in that row, are constant. Each row contains an indication of an event time. The only necessary change would seem to be to set the end time for the row with the event to the event time itself, with an indicator variable for the event.

The independence argument for the Poisson model does require that there be at most 1 event possible per individual* and that covariate values only change at time-interval boundaries. Under those assumptions, the argument holds more generally for proportional-hazards models. Under proportional hazards, an individual contributes to the (partial) likelihood at all event times in the data set while at risk. At each event time, however, only the instantaneous covariate values are used to estimate the relative hazards; covariate values at prior times are irrelevant, unlike for other approaches like accelerated failure time models. As the R time-dependent survival vignette explains (page 4):

One common question with this data setup is whether we need to worry about correlated data, since a given subject has multiple observations. The answer is no, we do not. The reason is that this representation is simply a programming trick. The likelihood equations at any time point use only one copy of any subject, the program picks out the correct row of data at each time.

If you need a fully parametric model for some reason, you can use tools provided by the R flexsurv or eha packages with counting-process data. The latter can take individual IDs into account in model types that require it.

---

*If an individual can experience more than 1 event, then you do have to take intra-individual correlations into account with methods like cluster variances, frailty terms, or a random-effects model.