Interrupted time series when most observations are zero

I am evaluating the impact of a medical intervention on an adverse outcome using administrative data collected over several years. The treatment was delivered individually. Therefore, each patient has a unique study period. There are 48 monthly observations (24 pre- and 24 post-intervention) for each patient. I have data on both treated and non-treated individuals. The groups are not equivalent - patients at greater risk of the outcome were more likely to be treated and spent less time at risk during any given month. Some patients died before the end of their 24-month post-intervention period. The outcome occurred rarely in both groups; in any given month, no more than five (< 1%) patients experienced the outcome. The majority of patients did not experience the outcome at any observation (i.e., all observations = 0).

A mock dataset is available below. Since the number of outcome events (out) is so rare, the data have been aggregated by treatment condition (tx).

structure(list(tx = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"),
ptNbr = c(1468L, 1468L, 1468L, 1468L, 1468L, 1468L, 1468L,
1468L, 1468L, 1468L, 1468L, 1468L, 1468L, 1468L, 1468L, 1468L,
1468L, 1468L, 1468L, 1469L, 1469L, 1469L, 1469L, 1469L, 1469L,
1469L, 1469L, 1469L, 1469L, 1469L, 1469L, 1469L, 1469L, 1469L,
1469L, 1469L, 1466L, 1466L, 1466L, 1465L, 1465L, 1465L, 1464L,
1388L, 1290L, 1222L, 1196L, 1176L, 674L, 674L, 674L, 674L,
674L, 674L, 674L, 674L, 674L, 674L, 674L, 674L, 674L, 674L,
674L, 674L, 674L, 674L, 674L, 674L, 674L, 674L, 674L, 674L,
674L, 670L, 670L, 669L, 669L, 668L, 665L, 664L, 650L, 635L,
627L, 613L, 602L, 589L, 575L, 567L, 560L, 556L, 549L, 544L,
534L, 527L, 512L, 502L), year = c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L
), month = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L,
12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L,
4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L,
8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L,
10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L,
12L), time = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L,
12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L,
24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L,
36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L,
48L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L,
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L,
26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L,
38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L), int = c(0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L), post = c(0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L,
13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L,
19L, 20L, 21L, 22L, 23L, 24L), out = c(10L, 16L, 2L, 14L,
11L, 9L, 11L, 1L, 11L, 4L, 6L, 2L, 4L, 6L, 5L, 0L, 8L, 15L,
8L, 2L, 2L, 8L, 2L, 3L, 18L, 2L, 5L, 2L, 12L, 13L, 3L, 11L,
5L, 5L, 9L, 11L, 9L, 9L, 3L, 2L, 20L, 5L, 0L, 1L, 3L, 13L,
6L, 8L, 0L, 7L, 5L, 7L, 0L, 20L, 3L, 2L, 21L, 8L, 17L, 4L,
1L, 16L, 0L, 2L, 29L, 27L, 8L, 8L, 5L, 18L, 8L, 12L, 5L,
0L, 0L, 7L, 0L, 5L, 0L, 2L, 32L, 6L, 2L, 1L, 3L, 3L, 0L,
3L, 2L, 0L, 0L, 0L, 0L, 2L, 0L, 23L), risk = c(44561L, 44609L,
44579L, 44520L, 44520L, 44561L, 44608L, 44523L, 44490L, 44407L,
44565L, 44541L, 44589L, 44645L, 44615L, 44529L, 44585L, 44614L,
44629L, 44651L, 44598L, 44427L, 44511L, 44539L, 44572L, 44678L,
44612L, 44520L, 44539L, 44558L, 44642L, 44590L, 44491L, 44371L,
44492L, 44294L, 44096L, 44419L, 44393L, 44254L, 44301L, 44169L,
43041L, 40778L, 37804L, 36706L, 36038L, 35313L, 19829L, 19986L,
19939L, 19936L, 19960L, 19939L, 19902L, 19844L, 19801L, 19869L,
19679L, 19720L, 19675L, 19785L, 19705L, 19802L, 19817L, 19790L,
19671L, 19591L, 19316L, 18273L, 16592L, 8660L, 20190L, 20067L,
20185L, 20252L, 20207L, 20133L, 20152L, 19852L, 19516L, 18978L,
18858L, 18414L, 17959L, 17505L, 17261L, 17189L, 16916L, 16745L,
16574L, 16279L, 15975L, 15748L, 15389L, 15090L)), row.names = c(NA,
-96L), class = "data.frame", .Names = c("tx", "ptNbr", "year",
"month", "time", "int", "post", "out", "risk"))


Questions

Since the treatment occurred at different times for each patient, aggregating the data would appear to prevent an analysis of any seasonal behaviour present in the outcome. Actually, it is my understanding that one advantage of using a multiple-baseline design is that there is no need to check for seasonality. Is this correct or an overly simplistic interpretation? If seasonality should be checked and modelled, must this be done at the patient-level?

Although the data have been aggregated to maximise the number of events that occurred in any given month, it seems more appropriate to model the zero and non-zero outcomes separately (e.g., through a multilevel [time within patient] approach with random intercepts and slopes for each patient), given so few patients experienced the outcome. If this is a more legitimate approach, should I model separately a time series for each patient and pool the results (as seems to occur in meta-regression)? If so, how would I model the majority of patients for whom the outcome at all observations is zero?