I have 400 discrete geographic areas of vastly different sizes, where I am recording the number of specific events. Every event is timestamped so I can summerise from daily counts to weekly, monthly, etc.

The number of events varies significantly between areas, mostly as a product of their size. Bigger area = more events. Not all areas will see these events, even after several months, however, the larger areas may see several tens per week.

These events, when all the areas are combined, have a very apparent seasonal pattern. Some of the larger areas will display this individually due to the number of events.

I have two main questions;

I changed something in one of these areas at a given date, where the intensity of the change will increase over time. Given the seasonality, which is the most appropriate test to demonstrate the change had an effect? How could I demonstrate that the rate of change is more significant in this area? Would simply tracking the monthly % change across all areas be sufficient?

Something changed in all of these areas at a given date, given the seasonality, which is the most appropriate test to demonstrate the change had an effect?

I am looking at Poisson regression, repeated measures ANOVA and ANCOVA.

The time series element, the seasonality and the very biased treatment vs non-treatment groups are throwing me off a bit.

It may be worth noting these changes all occurred in the real world, not in a controlled way.

In the following example, df would create the example data, where treatment1 refers to the data for the first question and treatment2 for the data for the second question. The actual treatment months do not align for each treatment, it serves only to create an example.

Please also note I have not created a seasonal variable for event_count.

df = data.frame(
  year = rep(2011:2020,
             each = 12),
  month = rep(month, 10),
  area = paste("Area", rep(1:6, each = 120), sep = "_"),
  treatment1 = c(rep("not treated", length.out = 100), rep("treated", times = 20), rep("not treated", times = 600)),
  treatment2 = c(rep("not treated", length.out = 100), rep("treated", times = 20)),
  event_count = round(rnorm(120, mean = 50, sd = 10), digits = 0)

Thank you in advance


2 Answers 2


Firstly, if you had made these changes in multiple places randomly, then a Poisson random effects model (random area effect that accounts for correlation within regions) for the number of events by time interval with a log total person-time at risk in the time interval (time * population at risk; or in your case perhaps time * area) offset would make sense. This essentially ends up being a model for events per area unit per time unit. If it is not likely/plausible that the event number is proportional to the size of an area, then you might have to model the area as a covariate (possibly using splines to allow for non-linear effects - there are methods for enforcing a monotonic effect, which would presumably be plausible).

The model would include an intervention effect (assuming the intervention effect is immediately there at full strength after intervention) and a fixed effect for every time interval. You could either include the logarithm of the total time * area at risk (or events per time, if the proportionality with area is implausible) before change (for the same time interval for all regions) into the model, or alternatively model those preceding time periods as observations. The second option is in contrast to treating that as fixed covariates; this could be quite appropriate, if there is no selection of areas you look at based on some previous outcomes but you rather took e.g. all US states or all counties in some states without selecting on previous outcomes but rather e.g. based on whether these locations have data available. The idea would be that other areas, where the change has not happened would account for seasonal trends. You can think of this as the count data version of the difference in difference model.

Secondly, this becomes a lot trickier, if intervention was not assinged by randomization. Then you end up having to consider the selection process that let to what units were selected for intervention. There's a bunch of options for doing that. These include propensity scores or covariate adjustment.


Under the assumption that you don't have readings for every day or that the data is "sparse" using monthly sums might be appropriate. In either case (daily or monthly) post one of your series showing starting date and country of origin and a brief description and I will try and help further. Note that weekly sums are not usually useful as activity in week x of year 1 is not usually duplicated in year x+1 except for very systematic processes AND of course how to deal with the leap year conundrum AND the impact of when exactly the holidays and or significant events have occurred in the past and will occur in the future.


120 monthly values for AREA1 were analyzed the first 100 were without any treatment ..the last 20 months reflected two treatments (BOTH) . Here is the data in graph form enter image description here and partially in list form enter image description here

To test the hypothesis that the treatment had an effect one needs to be concerned about possible anomalies clouding/obfuscating that test. A number of possible anomalies were identified and the model is here.

enter image description here with the treatment variable (called BOTH) not-significant. The Actual/Fit graph is hereenter image description here with the Actual/Cleansed here enter image description here

If there was seasonal activity be it seasonal deterministic effects or seasonal arima effects these would have to be identified and included in the model in order to assess the effect of the de jure intervention point. Software exists to handle this sometimes onerous task and would be mandatory for a large-scale problem as you have described.

If there was some untreated but important seasonal structure it would be in evidence in the acf of the model residuals. enter image description here

Often the de facto date of the intervention is different from the de jure date due to delays in the response to the treatment or a gradual response to the treatment. his is handled via INTERVENTION DETECTION and is discussed here and nearly everywhere http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html except in introductory/basic courses in time series analysis,.

  • $\begingroup$ Yes, data is sparse. I will add an example - thank you $\endgroup$
    – Will
    Commented Apr 21, 2020 at 12:12
  • $\begingroup$ actually post the data . $\endgroup$
    – IrishStat
    Commented Apr 21, 2020 at 12:47
  • $\begingroup$ What's the best way to format? Looks awful! $\endgroup$
    – Will
    Commented Apr 21, 2020 at 12:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.