# Cox model for recurrent events (with estimation of residuals / of an individual effect)

Consider some individuals that are followed during a period $T = 1$ (the same for all individuals). The individual indexed by $i$ have $n_i \ge 0$ events, at times $t_{i1}, t_{i2}, \dots, t_{in_i}$. There is a vector $X_{i}$ of covariates that influence the probability of event occurrence.

A simple model is a Poisson regression: $n_i \sim \mathcal P(\lambda_i)$ with $E(n_i) = \lambda_i = \lambda_0 e^{X_i \beta}$.

If after estimating the coefficients, I want some kind of individual effect, I could consider for example response residuals $n_i - \widehat{\lambda_i} = n_i - \widehat{\lambda_0} e^{X_i \widehat{\beta}}$.

Now, what's happening is that $X_i$ depends on $t$, and moreover the baseline risk $\lambda_0$ might as well depend on $t$. I can consider a counting process with intensity $\lambda_i(t) = \lambda_0(t) e^{X_i(t) \beta}$.

Do you think it’s the way to go? Can I fit this model (and, if possible, obtain response residuals or anything that can serve as a proxy of an individual effect) with a classical R package like survival?

• Shouldn't the expected number of events be the cumulative intensity instead of just the intensity? i.e. E(ni) over (a,b) = lambda_i (b-a)? Jun 16, 2016 at 12:57
• The period is the same for all individuals, so you can use at as unit time. This is irrelevant. Jun 16, 2016 at 20:39

What you mean by individual effect are actually the martingale residuals, that is $n_i - E(n_i)$. There is a whole theory on this and they are commonly used to assess model fit in Cox-type models. I am not sure about the interpretation of individual effect. They do have nice asymptotic properties, and you can find a chapter about them in most serious survival analysis books (I can give some references if you want to).

If you want to fit a semi-parametric model time-dependent model, that is possible with coxph(), as long as you have the data in the correct format. That is, one row should be (tstart, tstop, status, x) where status == 1 at event time points, and 0 for the end of follow-up, i.e. tstop == 1, or time points where x changes value. A data set that is in this format that I can think of is readmission in the frailtypack package.

To think about an individual effect, I would suggest to use an individual indicator as a factor and look at the interaction between this and $X$. However, I think that you need a lot of data to get something useful out of something like this.

• Thank you Theodor, this is helpful (the format of readmission is very clear). Concerning your final comment about introducing individual indicators to retrieve individual effects, would it be possible to model this effect as random and retrieve the equivalent of a BLUP? (I am more familiar with the linear mixed model...) Jun 17, 2016 at 9:16
• Yes these martingale residuals (thank you for pointing out the correct term) are not really individual effects but I think they still can be a useful proxy... does that sound stupid? Jun 17, 2016 at 9:21
• You can use random effects with the option +frailty(id), but the random effects are assumed to be independent of X. These are multiplicative effects on the hazard scale, i.e. lambda(t|z) = z lambda(t). You can get posterior estimates for z directly from the fitted object Jun 17, 2016 at 9:48
• As for martingale residuals, I do not think that they can be interpreted as such, but then again I am not an expert in that (so maybe it's actually possible) Jun 17, 2016 at 9:49
• Thank you for this nice answer and the reactivity in comments! Jun 17, 2016 at 10:19