Statistical comparison of counts between two different time periods

I have data on the drugs that people are taking at a given time point. I would like to observe the differences in which drugs people are taking at two different time periods e.g. comparing summer/winter. What I have is individual prescriptions for each patient, so I can count the number of times each particular drug was prescribed.

What I've currently done is calculate the distribution over drugs for each time period and gotten two probability distributions. Then, to identify which drugs are prescribed differently between the time periods, I've used a two-sample z-test of proportions. Because I'm looking at a lot of different drugs, I used FDR correction for multiple comparisons. To find interesting drugs, I look at the log-fold change and the p-value.

My first questions is: Does this sound reasonable?

My other concern is that, for instance, since the COVID pandemic started, I get the felling that a lot more prescriptions were given out of a particular medication e.g. vitamin D supplement. This could end up skewing the overall distribution. In this case, I'd guess that looking at the raw count data would be better than comparing distributions. What might be a reasonable way to approach this?

For using a z-test for proportion here, is the idea to look at "Patient was prescribed drug at least once in time period" vs. "Patient was not prescribed drug at least once in time period"? If so, that's a reasonable thing to look at, as long as the time periods are equally long and as long as intervals between prescriptions don't cause trouble (worst case: yearly renewal and all doctors try to do it in January, because budget is available, or other weird effects like that; best case: pretty much random prescription times renewed at substantially shorter intervals than your time periods of interest). However, comparing that using a two-sample test seems wrong, because you'll presumably have a lot of patients occuring in both time periods. I.e. these could be looked at as paired proportions, which you can analyze in a number of ways. For example, you can look at this as a contingency table of "prescibed in summer" vs. "prescribed in winter" and methods for contingency tables - there are extensions if you only have data for some patients in one season. You can also use a random effects logistic regression with a random patient effect, e.g. in R this would look like this:

library(lme4)
glmer(y ~ (1|patient) + season, family = binomial)


The nice thing about the logistic regression is that it extends sensibly to the setting with multiple types of prescriptions e.g. by adding a prescription type by season interaction and allowing the patient random effect to vary by prescription type, e.g. like this:

glmer(y ~ 0 + (prescription_type|patient) + season + prescription_type + season*prescription_type,
family = binomial)


You could then look at the season by prescription type interaction (and its p-value) and apply multiplicity corrections as needed. You can also adjust for other variables, if that makes sense. Perhaps taking a Bayesian perspective and embedding the possibility for diverging effects inside a hierarchical model (see e.g. these comments by Frank Harrell or e.g. the horseshoe prior, which let's you encode your a-priori belief about likely number of meaningful interactions could be an approach).

If there are additional things going on beyond seasonality (e.g., as you mention, COVID-19, or a general long-term trend for more prescriptions, or the average person in your database having gotten a bit older etc.), these will confound any of the simple ways of looking at this, no matter whether you look at aggregated numbers or individual counts. Especially if you look across many years, you can try to adjust for such effects (e.g. by having a "this time period is during the COVID-19 pandemic"-effect with interactions to specific prescriptions in the model), but this is why it's hard to causally attribute differences in prescriptions to a causal effect of it being a particular season (as opposed to something else that happened at the same time). In that sense the models I describe above are perhaps more descriptive.

You may also have the issue that within the same person, many chronic drugs may be taken for a long time once initiated, so having them in the first time interval for some patients likely means that more people (the original ones plus some more) will take them in the second time interval. If your conclusion changes when you take summer first or winter first, then this is a sign this might be going on (another argument for looking at multiple years).

If I understand you correctly you have two measurement times of drug prescription, at summer and at winter for each patient. There is d possible types of drugs.

Your question is if there is any difference in which drugs are taken at the specific time.

Assuming that each drug is measured as yes/no:

One method is to consider a mixed binominal regression model with summer/winter as outcome, drug as independent variable, and id as a mixed effect. You can then test for no association between drug and time a year.

A second method, is to consider each drug individually. Form a 2*2 table and use Macnemars test for paired observations, to se if there is any changes in prescription of that specific drug.

Assuming the drug is a continuous measurement:

You can create a linear mixed model, with drug amount as outcome, drug type and season as fixed effects, and id as random effect, and test for no effect of season.

Or you can use a paired t-test to compare the changes of drugs prescribed in different seasons.

About the method you suggests If you use the proportion test, you will ignore that individuals may have a measurement at both seasons. If you can argue to ignore this dependence, you are fine. However that means that you ignore that a person having drug a prescribed in summer might be more likely to have the same drug prescribed in winter.

Correcting for multiple testing I am not familiar with the FDR correction, however I would recommend not to do this unless you are doing a lot if test - say 1000. You can read the following paper to see why. In short, when you decrease the chance of type 1 errors doing multiple testing, you simultaneously increase the chance of type 2 errors. At the very least I would strongly recommend not to present adjusted p-values without also presenting the unadjusted p-values.

About you concerns with COVID If you have an expected effect of, lets say COVID on the prescription of drugs, you can correct for this (I will leave the how for another question, to keep to subject). If you have no prior knowledge about the effect, I would recommend to leave this part of data out of the analysis if possible. If not, try to do sensitivity analysis where you leave that part of the data set out, again if possible.

I hope this helps you with you question. If I have misunderstood anything feel free to add correction comments.

• The viewpoint that no multiplicity adjustments are needed is a topic of ongoing debate. If one cares about the validity of p-values and controlling either the familywise type I error rate (or the false discovery rate), then this is of course required in a frequentist paradigm (otherwise, go Bayesian, but with proper joint priors on effects). Jun 28, 2021 at 11:46
• An interesting view-point, if you have any articles which clearly states the gain of adjusting is higher than the dis-advantages of adjusting the p-values in simple tests/ or regressions, I would very much like to read them. As far as I know most of the gain by adjusting p-values can be obtained simply by taking care in interpretation of the result, when remembering that a significant result is not a final solution but a mere limit made to simplify the test-result to a binary conclusion option. Jun 28, 2021 at 12:08