Modelling rare outcome in treatment evaluation This question is related to my previous question.
I am conducting a treatment evaluation in a retrospective cohort study. My dataset has 2000+ cases, each with 48 monthly observations (24 pre- and 24 post-treatment). Treatment was delivered individually and the study has been conceptualised using a multiple baseline design.
I would like to assess the change (rates) in the outcome pre- / post-treatment and compare it to the change in the outcome among non-treated patients over the same period. The outcome occurs rarely in both treated and non-treated cases (the vast majority of cases do not experience the outcome at any point).
I am seeking guidance on modelling techniques that may be appropriate. I have considered multilevel modelling (time within patient) and interrupted time series; however, given the outcome occurs rarely, I am unclear whether the data should be modelled at the group or individual level. I have also considered using conditional poisson regression, but I am unsure how to implement this in R.
Assistance required
Would anyone be able to provide some guidance on modelling techniques that are appropriate for use when the outcome is rare?
Also, would anyone be able to suggest how to implement conditional poisson regression in R?
I have attached a mock (case-level) dataset below. I am aware it is not in the correct format for longitudinal analyses, but I am hoping my question can be answered without needing access to data.
structure(list(id = 1:20, tx = structure(c(1L, 2L, 2L, 1L, 2L, 
1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L), .Label = c("No", 
"Yes"), class = "factor"), preNbr = c(0L, 0L, 4L, 2L, 10L, 2L, 
0L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 9L, 0L, 1L, 4L, 5L, 0L), postNbr = c(2L, 
0L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 5L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L)), .Names = c("id", "tx", "preNbr", "postNbr"), row.names = c(4L, 
31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 40L, 41L, 42L, 43L, 46L, 
47L, 48L, 49L, 51L, 52L, 56L), class = "data.frame")

 A: Caveat: I am no medical nor statistical expert, this answer reflects my best judgement, but when medical decisions are involved you should not depend on random Internet advice. You should only use analysis methods that you understand and that you can verify to make sense in your case. Please don't use my answer blindly. Also it would be much better had you determined analysis methods prior to data gathering, but I understand the reality of medical research is difficult and optimal choices are not always possible. 

The actual answer:
When you are in doubt whether to model at individual or group level, use some kind of partial pooling - preferably a fully Bayesian hierarchical model and you get the best of both worlds. Also, fully Bayesian approach handles rare outcomes well. 
Lets assume there is a population-level risk $\alpha$, and each patient has an additional fixed individual risk term $\nu_i$ (where $E_i(\nu) = 0$). This patient risk is then modified by treatment and the relative treatment effect $\beta$ is the same for all patients. We will use a Poisson link (because you assume it), but I am not sure it is the right fit for your case - really depends on your data. If the outcomes are correlated a negative binomial or even something completely different might be better. The resulting model could look something like:
$Y_i^{post} \sim Poisson(exp(\mu_i^{post}))$
$\mu_i^{post} = \alpha + \beta X_i + \nu_i$
$Y_i^{pre} \sim Poisson(exp(\mu_i^{pre}))$
$\mu_i^{pre} = \alpha + \nu_i$
$\nu_i \sim N(0,\tau)$
$\tau \sim HalfNormal(0,\sigma_{\tau})$
Here $X_i$ is zero or one, for the control and treatment group respectively. Note that treating $Y_i^{pre}$ as a random variable lets us handle measurement error in the pre-treatment data. When fitting such a model, $\tau$ binds together all the $\nu_i$ and lets them share information - if patients differ a lot, it will reflect in larger $\tau$ and higher uncertainties in $\nu_i$, if the patients are more similar, $\tau$ will be smaller and you will be more certain about the values of $\nu_i$. Here $\sigma_{\tau}$ is a hyperparameter that you need to set, using your judgement on how much between-patient variability is reasonable (If we don't expect the log-rate to differ by more than 10, we may put $\sigma_{\tau} = 10/2 = 5$ ). Since you have quite a lot of data, the choice of $\sigma_{\tau}$ should not matter that much.
Note that such model cannot be fit by classical/frequentist methods (the ML estimate may be highly unreliable), and you have to go Bayesian. But you get the added value of being able to say things as "The probability that the treatment effect is > X is Y%" and having a measure of between-patient variability which may be of interest on its own.
Now how to actually fit such a model? If you want full control and don't mind the learning curve, the above probabilistic model can be translated almost line by line to Stan (which integrates seamlessly with R).  The Stan manual has a multiple examples of hierarchical models (e.g.  "Hierarchical Logistic Regression", and "Hierarchical Models and the Non-Centered Parameterization"). There is also a nice case study on the topic with Python code and I have coded a very similar model for my blogpost. If you go down this rabbit hole, don't be afraid to ask at Stan's Discourse, the community is helpful. Also in this case, it is easy to modify the model, if your assumptions are different than the ones I have stated.
A simpler, but somehow limited approach is to use the INLA R package where you can use the familiar formula syntax (with a few extensions). The main limitation is that we now have to treat $Y_i^{pre}$ as a predictor and not consider the measurement error. The maximum likelihood estimate for $Y_i^{pre}$ is $log(Y_i^{pre})$ but since there are zeroes in the data, we need to use something like $log(1 + Y_i^{pre})$ instead. The following R code works with the sample data you provided:
data <- structure(list(id = 1:20, tx = structure(c(1L, 2L, 2L, 1L, 2L, 
1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L), .Label = c("No", 
"Yes"), class = "factor"), preNbr = c(0L, 0L, 4L, 2L, 10L, 2L, 
0L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 9L, 0L, 1L, 4L, 5L, 0L), postNbr = c(2L, 
0L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 5L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L)), .Names = c("id", "tx", "preNbr", "postNbr"), row.names = c(4L, 
31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 40L, 41L, 42L, 43L, 46L, 
47L, 48L, 49L, 51L, 52L, 56L), class = "data.frame")

library(INLA)

fit_inla <- inla(postNbr ~ tx + f(id,model="iid") + log1p(preNbr), family = c("poisson"), 
           data = data, control.predictor=list(link=1)) 

#Show the results
fit_inla$summary.fixed
fit_inla$summary.hyperpar

INLA also treates $\tau$ slightly differently than in the above probabilistic model, but this should not matter since you have a lot of data. Read the docs on priors in INLA for more details.
Good luck with your modelling.
