Statistical analysis:
An interrupted time series analysis using ARIMA forecasts, a specific form of a “before” and “after” analysis, was used in order to quantify the pandemic (as population-level health “interventions”) and its outcomes by relying on the natural expirement that results from the intervention (Schaffer et. al, 2021). However, results relying on simple pre- and post intervention designs are prone to bias because they do not account for underlying short and long-term trends. It is therefore important to determine whether the “effects” of the intervention are different from what would have happened anyway (Soumerai et. al, 2015). Interrupted time series (ITS) analysis is more robust to this bias as it controls for issues such as underlying trends by longtitudinally tracking outcomes before and after the intervention (Schaffer et. al, 2021). ITS designs are highly regarded for their rigor and can arguably be considered the best possible approach for evaluating the impact of new policies implemented nationwide, with an increasing use in health research (Ewusi et al., 2020). By including another time series as a control, the internal validity of ITS designs can be further increased. Lopez Bernal et. al (2018) suggest that controlled interrupted time series (CITS) analysis can achieve similar results to randomized controlled trial studies (RCT) benchmarks, offering an option when RCTs are neither feasible nor ethical (Schaffer et. al, 2021). An ITS design using Autoregressive Integrated Moving Average (ARIMA) models regresses the outcome Yt only on the outcome measured at previous time points, with the advantage of being able to handle autocorrelation, trend and seasonality in a time series (Schaffer et. al, 2021).
For the CITS design we define the implementation of Covid-19 measurements in Austria, such as restrictions and lockdowns, as a population-level health intervention that can be analysed as a natural expirement. In line with our hypothesis, we expect that this intervention is affecting the mental health of the Austrian population, with a much more pronounced effect among adolescents and especially female adolescents, based on the international literature about the pandemic impact on the mental heath of different age groups. We therefore expect to see this effect in the amount of antidepressant, antipsychotic and benzodiazepane patient rates.
The population group of interest, Austrian adolescents, have seen a long period of distance schooling from the beginning of November 2020 to February 8th 2021, but the first lockdown took place already in March 2020. Since the first lockdown was restricted to a few weeks, we hypothesized that an expected increase of prescriptions of antidepressants, antipsychotics and benzodiazepines will be seen starting in 2020 Q3. Further, since the mental health strains continue with the duration of Covid-19 measurements, we expect the effect to be cumulative. So we define the start of the intervention as 2020 Q2 and expect patient rates to rise from 2020 Q3 onwards until the end of our observation period. School closings affected a large part of the adolescent population in Austria (including apprentices at their training sites), while the whole population was subdued to lockdowns. Since both groups were affected by Covid-19 measurements, we used the patient rates in the whole population as a control series to see whether adolescents were affected stronger by the measurements.
The dataset contained quarterly data on patient rates and total patient numbers from 2013 to 2021 for antidepressants (ATC Code N06A), antipsychotics (ATC Code N05A) and benzodiazepines (ATC Code N05BA) as well as individual medications belonging to these groups. Each pharmaceutical group (and each medication) was split into three age groups and two genders: 10–14 years, 15–19 years and the whole population (including all age groups) as well as male and female gender. For SSRIs (ATC code N06AB), we decided to analyze the group with regard to specific substances, which have been recommended and included in different systematic reviews and meta-analyses on antidepressants in adolescence (Cipriani et al., 2016; Hetrick et al., 2021). We included antipsychotic agents with a described off-label use in children and adolescents (e.g. van Schalwyk et al., 2017; Shafiq & Pringsheim, 2018) as well as selecting Diazepam, which has a license for use in acute aggressive behavior in children and adolescents. To strenghten the forecast, we defined a minimum sample size of 100 patients per quarter and eliminated groups in which less than 100 patients per quarter were observed on average. This led to the exclusion of the groups in Table 2 in the appendix.
We used a variation of the ITS method demonstrated by Niederkrotentaler et al. (2019) and a more detailed description by Schaffer et al. (2021) to forecast patient rates. First, the data was split into 72 time series, one for each respective age group, gender and ATC combination. Then the time series development as well as ACF and PACF plots were inspected between Q1 2013 and Q2 2020 to identify trends, seasonal or cyclical patterns and autocorrelations in the data for subsequent ARIMA modelling. Where heteroscedascity was visible, log transformations were applied. Second, a seasonal ARIMA model was fitted to each time series from 2013 Q1 up to the intervention in 2020 Q2. The analysis was performed using R (R Core Team, 2022), the packages tidyverse, lubridate and tsibble for transformation of the dataset (Wickham et al., 2019; Grolemund & Wickham, 2011; Wang et al., 2020) and the package fable to fit and select ARIMA models (O'Hara-Wild et al., 2021 A). Each model was fit by choosing the optimal model according to the fable packages auto model selection (with transformations applied according to prespecification), which uses a variation of the Hyndman-Khandakar algorithm selecting for the smallest AICc value (Hyndman & Athanasopoulos, 2021; Hyndman & Khandakar, 2008). The stepwise procedure and approximation options of the model selection function were disabled to ensure the consideration of a larger modelspace in the search and a ljung-box test was applied to the selected model using the feasts package (O'Hara-Wild et al., 2021 B) to see whether the residuals resembled white noise. Third, 97.5% confidence bands were forecasted using the model for the six quarters after the intervention date. These forecasts can be interpreted as the expectable development path as they show how patient rates would have likely developed if they followed their historical developments in the absence of an intervention. We then compared the observed data from 2020 Q3 to 2021 Q4 to the confidence intervals obtained from the models. Each upwards deviation of the observed values outside of the confidence bands was classified as a significant deviation from the expectable path at the 5% level, as demonstrated for example by Niederkrotentaler et al. (2019).
We further illustrate the developments with absolute numbers of patients in the respective age groups before the onset of the Covid-19 Pandemic (2019 Q4) and at the end of our observation period (2021 Q4) as well as the index developments between these dates with the reference point of the index set to the 2019 Q4 (see Table 1). All absolute numbers behind the developments for the categories antidepressants, antipsychotics and benzodiazepanes are shown in table 3 in the appendix, while absolute numbers for individual medications could not be published due to data restrictions.
As the dataset provided by the Federation of Austrian Social Insurances consisted of accumulated data without possibility for identifying individualized data, a waiver was received from the ethical review committee of the Medical University of Vienna.