1 Introduction

Emergency medical services (EMS) represent the first point of healthcare access for millions of patients worldwide in need of urgent treatments. Despite the heterogeneous organization of EMS agencies globally, they all share the primary goal of minimizing intervention time. To achieve this goal with limited resources, statistical models for forecasting the number of emergency calls and rescue missions play a crucial role in informing EMS agencies about future activity volumes, which are subject to significant daily and seasonal fluctuations. These fluctuations are influenced by seasonal weather variations and extreme weather events (Alessandrini et al. 2011; Bassil et al. 2009; Lin et al. 2013; Attia and Edward 1998; Noble et al. 1971), such as heat waves or cyclones (Noji 2000), social and demographic factors (Schuman et al. 1977; Kamenetzky et al. 1982), and epidemiological phenomena (Al Amiry and Maguire 2021; Diaz et al. 2001), such as influenza epidemics. Additionally, the day of the week and the time of the day impact the volume of EMS activity on different timescales (Cantwell et al. 2015; Batal et al. 2001). In this complex framework, predictive models that account for all these factors can assist EMS organizations in optimal staff management and distributing ambulances and rescue equipment.

The importance of reliable predictive models became evident with the emergence of the COVID-19 pandemic. During that period, EMS organizations faced sudden and significant changes in the volume of emergency calls, and existing predictive models failed to capture this phenomenon. Thus, forecasting the number of calls and rescue missions in the near future is essential to optimize resource allocation and address the rapidly evolving needs of the served populations.

In this study, we describe the development of a predictive model to reliably forecast the number of emergency events in the Lombardy region of Italy. This research effort was driven by a request from AREU (Agenzia Regionale Emergenza Urgenza), the regional agency responsible for the organization and implementation of EMS, for a reliable predictive model to inform organizational decisions during and after the COVID-19 outbreak. Before the COVID-19 pandemic, AREU routinely used a predictive tool based on an unobserved component model for time series data (Harvey 1990). However, this model was only developed to forecast events in the metropolitan area of Milan, the region’s capital and most populous city in the Lombardy region, and it performed poorly when the COVID-19 pandemic began. Therefore, a new family of predictive models covering different areas of Lombardy and accounting for COVID-19-related factors was urgently needed.

The aim was to provide AREU with a predictive model that is (i) computationally efficient and time-effective, allowing for daily updates and new daily predictions, (ii) feasible in terms of accessible and up-to-date covariates, and (iii) flexible to accommodate epidemiological and social changes (Kenett and Gotwalt 2022).

To fulfill this purpose, we propose a generalized additive model (GAM, Hastie and Tibshirani 2017) with a negative binomial family, nonlinear autoregressive terms (Kedem and Fokianos 2002), and exogenous covariates, which can deal over-dispersed count time series (Agresti 2015). The proposed model is flexible, allowing for nonlinear dependence between covariates and outcomes. Although the same family of models (i.e., GAM) has already been used to develop predictive models in a similar framework (e.g., Viglino et al. (2017)), our approach is the first to simultaneously address nonlinearity, auto-regressive effects (Kedem and Fokianos 2002), over-dispersion (Agresti 2015), time (Kedem and Fokianos 2002), weather (Attia and Edward 1998; Noble et al. 1971) and epidemiological (Diaz et al. 2001) factors. For a recent and exhaustive review of the literature on predictive models in the field of EMS, see Al-Azzani et al. (2021); Huang et al. (2021) and the references therein. The developed model was compared, in terms of predictive performance, to alternative forecasting methods, including the predictive tool used by AREU before the COVID-19 pandemic.

We focused our modeling efforts on predicting the volume of emergency events, defined as the dispatch of the most appropriate transport until the rescue is completed and/or the patient is referred to the most suitable hospital facilities. Planning emergency transportation dispatches is the most crucial organizational effort from a managerial and planning perspective. While most of the literature on forecasting models in EMS has focused on predicting emergency call volumes, trends in emergency call volumes do not always align with changes in actual emergency events, which are the ultimate organizational targets. In Lombardy, approximately half of the emergency calls are not forwarded to the second-level public safety access point because they are deemed non-urgent upon initial assessment (AREU 2022). Various factors, such as media communications and social alerts, influence the proportion of calls that do not require emergency services. For example, immediately after the first COVID-19 case was confirmed in Italy, an unprecedented increase in emergency calls was observed in Lombardy. Many of these calls were requests for information and recommendations about the new disease (Marrazzo et al. 2020; Dattner et al. 2022; Molenberghs 2023).

The fact that our study aims to predict the number of future emergency events is a strength. Specifically, having a reliable model to forecast such events allows the central direction unit to activate additional ambulances and increase the number of emergency teams on shift. As reported by AREU, the organization in charge of EMS in Lombardy, implementing such actions requires at least 24 hours. This designated time frame allows sufficient flexibility to adjust the number of ambulances. It ensures effective planning and response to potential emergencies, maintaining an optimal balance of resources and providing timely medical assistance to those in need (Berchi et al. 2010; Gilardi et al. 2021). Therefore, we considered predicting events one day in advance as the shortest forecasting time frame. In addition, since that data is updated hourly, the proposed model allows for considering shorter prediction intervals if deemed necessary. The model presented in this paper has been implemented into user-friendly software. Following the initial testing phase reported here, the software has been successfully used by AREU for several months.

The paper is structured as follows. Subsec. 2.1 describes the Lombardy EMS, which is the target of our research. The data sources utilized in the model development are outlined in Subsec. 2.2. The model development process is described in Subsec. 2.3, while Sect. 3 presents the estimates and predictive performance of the model. In Sect. 4, we discuss the strengths and limitations of our study.

2 Methods

2.1 EMS organization in Lombardy

Lombardy has a population of 10 million, the most populous region In Italy, accounting for approximately \(17\%\) of the Italian people. Before the COVID-19 pandemic, the regional EMS received around 1 million emergency calls annually, corresponding to over 900,000 emergency events. AREU, the public agency responsible for coordinating, guiding, managing, and monitoring the out-of-hospital emergency network and emergency number service in Lombardy, oversees these operations. The agency is structured into peripheral units known as Territorial Business Units (AAT, Articolazioni Aziendali Territoriali) as well as Regional Emergency Operations Rooms (SOREU, Sale Operative Regionali dell’Emergenza Urgenza). The twelve AATs are distributed throughout the region, approximately following the provincial division, while the four SOREUs coordinate rescue operations in super-provincial competence areas (i.e., AAT aggregations). The geographical position and organization of the SOREUs and the distribution of the AATs are depicted in Fig. 1.

Fig. 1
figure 1

Description of the AREU organizational structures: a four SOREUs: AREU emergency management; b twelve ATTs: AREU management of territorial resources

Full size image

2.2 Data

A comprehensive set of variables pertaining to emergency events in Lombardy from January 1st, 2015, until May 9th, 2021 were considered potential predictors for model development. We focused on predictors that are continuously available to AREU, enabling decision-makers to effectively utilize the predictive models for short-term forecasting (Kenett and Gotwalt 2022). A total of \(134\) variables were collected from various data sources: the AREU databases for emergency response data; the Regional Agency for Environmental Protection (ARPA) registry for weather data; the Department of Civil Protection repository for COVID-19 data; the Higher Institute of Health (ISS) estimates for other epidemiological information.

The AREU database contains information about all received calls: the SOREU receiving the call, the exact time of the call (i.e., date and hour/minutes), the call classification (e.g., first aid), the AREU administrative area (i.e., province, zone, AAT), call initiation location (e.g., home, street) with geographic coordinates, the reason for the call, the severity code (i.e., triage), and whether the call triggered an aid response, resulting in an event. The data are aggregated at an hourly level.

Weather data from ARPA data are obtained through the Open Data project (ARPA - Regione Lombardia 2014), which collects temperature, rainfall, and snowfall data from sensors distributed across the Lombardy region. For each AAT, central sensors in the city center are selected, and the weighted average is computed at the SOREU level. The weights equal the proportion of calls (from January 1st, 2015 until May 9th, 2021) received from each province in the SOREU. The weather data are available at an hourly level but were aggregated at the daily level for our analyses, using indicators such as the daily average or the difference between the maximum and minimum temperature of the day.

COVID-19 data from the Department of Civil Protection are obtained from Dipartimento della Protezione Civile (2020). The data cover the Lombardy region, except for the total number of positive cases reported at the province level. The data are available daily. In addition to the basic epidemiological information, we compute the effective reproduction number (Rodpothong and Auewarakul 2012; Dabbaghian and Mago 2014) using the method employed by the ISS, available at Istituto Superiore di Sanità (2020).

Finally, we consider the weekly flu incidence at the national level from 2015 to 2021 as communicated by ISS (Istituto Superiore di Sanità 2020).

2.3 Model development

Let \(Y\) be a random vector of dimension \(T\times 1\), with its components \(Y_{t}\) representing the counts of events at time \(t\in {\mathscr {T}}\), \({\mathscr {T}}:= \left\{ 1,\dots , T\right\}\). In our study, the time index \(t\) is on the hourly scale, and the spatial aggregation unit was the SOREU. Therefore, we targeted the prediction of the number of emergency events for each hour of the day in each of the four SOREUs. Let \({\textbf{X}}=\left( X_1, \dots , X_K \right)\) be a set of \(T\)-dimensional covariates associated with the response variable \(Y\), indexed by \(k \in \mathscr{K}, \mathscr{K} := \{1, \dots ,K\}\). Covariates can be year-, month- or day-specific and can be lags of \(Y\), such as the total number of events in the previous days.

Given the important territorial and urban differences of the SOREUs, characterized by extremely different orographic characteristics and distributions of metropolitan and rural areas, we developed separate models for each area.

We opted for a GAM model to capture nonlinear relationships between \(X\) and \(Y\) with smoothing splines (Wood 2004, 2011, 2017). Since the response was a count variable and preliminary exploratory data analyses showed evidence of over-dispersion (Agresti 2015), we assumed a negative binomial distribution (Casella and Berger 1990) for the response. Specifically, we assumed \(Y_t\sim \mathscr{N}\mathscr{B}\left( r_t, \pi \right)\), where \(t \in \mathscr{T}\), with \(E\left( Y_t\right) = \mu _t=\frac{r_t \pi }{1-\pi }\) and \(V\left( Y_t\right) = \mu _t + \frac{\mu _t^2}{\pi }\). The expected value \(\mu _t\) of \(Y_t\) is assumed to depend on the covariates as follows:

$$\begin{aligned} \ln (\mu _t) = \alpha + \sum _{k = 1}^{K} f_k(X_{tk}), \quad \text {where} \quad f_{k}(X_{tk}) = \sum _{d = 1}^{D_k} \beta _{kd} b_{kd}(X_{tk}). \end{aligned}$$

Here, the coefficients \(\beta _{kd}\) are the unknown parameters to be estimated, while \(b_{kd}(\cdot )\) are known basis functions, and \(D_{k}\) is the number of bases for the covariate \(X_{tk}\).

We had to perform variable selection to find the optimal model to predict the number of events in each SOREU. The variable selection process was based on the data from the Plain SOREU, the second largest SOREU. We did not use the Metropolitan SOREU, which is the largest among the four, because AREU already had a predictive model for this area, and the development of a model for the remaining SOREU was their priority.

Given the large number of available explanatory variables (i.e., \(134\)), we relied on assumptions a priori, based on the existing literature, about their effects on the count of events and their possible interactions. We assumed the presence of a time effect at the annual, weekly, and daily levels. In addition, we assumed an effect of weather (Attia and Edward 1998; Noble et al. 1971), epidemiological (Diaz et al. 2001) and auto-regressive factors (Kedem and Fokianos 2002). We constructed several sets of covariates to explain seasonal, weather, epidemiological, and auto-regressive effects with the available variables. For example, we explored whether the seasonal effect was better modeled at the month or quarter levels. In the first case, the variable “month” was entered into the model, while in the second case, we entered the variable “quarter”. Similarly, the meteorological impact was analyzed by considering the average daily temperature and the difference between the daily minimum and maximum temperature. Again, the impact of the COVID-19 spread was analyzed through the number of hospitalized patients with symptoms and the number of positive swabs. These are just a few examples of the combinations analyzed in this variable selection process.

After that, we analyzed this set of covariates through exploratory analysis to understand the underlying relationships between the dependent variable (i.e., count of events) and covariates and among the covariates themselves. We finally selected five models with different sets of covariates \({\textbf{X}}\) cross-validated across different time periods. Table 1 briefly describes the set of variables inserted in these five final models. Some of them are considered categorical variables in some models, and others as continuous variables inside a spline function. The \(\star\) symbol in Table 1 expresses if the ordinal variable is considered a continuous variable inside a spline function.

Table 1 Description of the five models considered in the cross-validation step
Full size table

The observations in time–series data are temporally dependent, which means that traditional cross-validation procedures do not apply directly. To tackle this issue, a rolling cross-validation approach suitable for temporal data was preferred. Particularly, the candidate models were fitted to data \(Y_1,\dots ,Y_t\), forecasts \({\hat{Y}}_{t+1}\) were obtained, and the relative prediction errors \(E_t=\frac{{\hat{Y}}_{t+1}-{Y_{t+1}}}{{Y_{t+1}}}\) were calculated for each model. This was repeated for \(t=m,\ldots ,m+N-1\), where \(m\) corresponds to the minimum size of the dataset used to fit the models and \(N\) is the number of forecasts. The mean absolute percentage error (MAPE) for each model, computed across \(N\) forecasts, was used as a performance metric in the cross-validation step:

$$\begin{aligned} MAPE = \frac{\sum _{i = 1}^{N}|E_i|}{N} 100, \end{aligned}$$
(1)

where \(|\cdot |\) is the absolute value, and \(E_i\) is the relative prediction error from forecast \(i\), \(i=1,\dots ,N\). The choice of the relative prediction error for assessing forecasting performance was based on matching the performance measure used by AREU. However, the user-friendly software where the final model is implemented also reports other performance metrics, such as the minimum absolute percentage error, the maximum absolute percentage error, and the median absolute percentage error.

Table 2 shows the MAPE considering the five different models described in Table 1 cross-validated on a rolling basis across different time periods for the Plain SOREU. The MAPE is computed across 30 days considering one, two, five, and seven days ahead predictions in 2016, 2018, and 2020. Therefore, for each year, 30 days were randomly selected. Then, in the first step, the training set is composed of the period before the first day of the 30 days randomly selected. The test set comprises the first, second, fifth, and seventh days of the 30 days selected to have predictions one, two, five, and seven days ahead. In the second step, the first day of the 30 days randomly selected is added to the training set, the remaining 29 days are used as a test set, and so on. Model 1 has the lowest cross-validated MAPE overall compared to the other four models.

Table 2 MAPE from five different models described in Table 1 calculated across different time periods (30 days) considering respectively one, two, five, and seven days ahead predictions for the Plain SOREU
Full size table

The predictive performance of the developed models was compared to that of other alternative approaches. The first comparator was the model used by AREU for the Metropolitan SOREU, which was based on an unobserved component time series model (Harvey 1990). The model included several covariates such as a dichotomous variable for increased calls (e.g., snowfall), a dichotomous variable for decreased calls (e.g., midweek holiday), predicted average daily temperature, flu incidence value (estimated), and the square of the difference between the average daily temperature and a reference temperature (\(16^\circ\)C). The model also considered stochastic seasonality at the weekly level. The developed model was also compared to other benchmarking models estimated on the same training data. These included a naive, deterministic model where predictions were simply set to the latest available observation, an ARIMA model based on the time series of previous events (Makridakis and Hibon 1997), and a generalized linear model (GLM) for count time series (Christou and Fokianos 2014).

We use the open-source statistical software R (R Core Team 2015) for preprocessing and model fitting. Among the available R packages for fitting GAMs, we utilized the mgcv (Wood 2011) package, which enables the fitting of GAMs on large datasets through the bam function. The ARIMA and GLM benchmarking models were fitted using the auto.arima R function and the tsglm function from tscount package (Liboschik et al. 2017), respectively. The full code used in this paper is available at https://github.com/angeella/Tsunami_project.

3 Results

3.1 Model estimates

This section will focus on presenting the model developed specifically for the Plain SOREU, which was used for variable selection. The covariates included in the optimal model are described in the first model of Table 1.

We employed cubic regression splines (Green and Silverman 2019) to capture the seasonal effects of hours and quarters. The cubic splines were chosen because they match up to the second derivative with their start, which is suitable for modeling cyclical effects. Specifically, we used \(24\) basis functions for modeling the seasonal effects of hours and \(4\) basis functions for modeling the seasonal effects of quarters.

For the day variable, we opted for P-splines with \(7\) basis functions (Eilers and Marx 1996) to capture the day’s seasonal effects. Finally, the tensor product smooth (Ramsay et al. 1997) (two-dimensional smooth, where the shape of one dimension varies smoothly over the other dimension) between day and hour covariates is applied to analyze the interaction between these two variables. Linear associations were assumed for the remaining covariates, which were included in the model without any smoothing terms. P-splines performed similarly to cubic regression splines but are numerically more stable and easier to fit. In addition, this permits a fast estimate of the smooth tensor product. The time complexity of the model fitting is essential in this case since AREU needs a computationally fast model to forecast the number of emergency events (Kenett and Gotwalt 2022).

For each covariate, we analyze the estimated association with the outcome. Figure 2a illustrates the estimated seasonal effects of hours within a day. We observe a positive effect during the late morning and afternoon when people are typically awake and active. In contrast, the effect becomes negative during the night and early morning when people are usually asleep. Figure 2b displays the effects of the day of the week. Mondays show the highest incidence of events compared to the other days of the week. Figure 2c presents the seasonal effects of quarters. The number of events is generally lower in the summer when the schools are closed, and people are on vacation. This observation may be associated with rising temperatures during that period. Finally, looking at the interaction between hours and days in Fig. 2d, we can see a positive effect in the morning (9 a.m.–12 p.m.) of weekdays and a strong negative effect in the morning of the weekend. Furthermore, we can see a positive effect during the night on the weekend and a negative effect around 5 a.m. on Monday, Tuesday, and Wednesday.

Fig. 2
figure 2

Plain model results: coefficients plots during the COVID-19 pandemic, a Effect of the hour of the day; b Effect of the day of the week, where \(1\) stands for Monday; c Effect of the quarters of the year, where \(1\) stands for the first quarter of the year; d Effect of the interaction between days and hours, where \(1\) stands for Monday. The black dotted lines represent two standard error bounds

Full size image

We estimate the model considering data before the COVID-19 era to see if the temporal dynamics, particularly the seasonal effects, have also changed due to the pandemic. Figure 3 shows the same plots presented in Fig. 2. We can see that the effect of the hour of the day (Fig. 3a) and of the quarter (Fig. 3d) do not change markedly, unlike the effect of the day of the week (Fig. 3b). One possible interpretation could be that before the COVID-19 pandemic, many emergency events were concentrated during the weekend. Therefore, with the arrival of the COVID-19 pandemic, the emergencies referred to AREU have shifted by a few days. People become infected during the week and weekend and seek care on Monday. Before the COVID-19 pandemic, on the other hand, the primary demand, as we have said, was at weekends, when people went out to parties/hiking or when general practitioners were simply not available. Furthermore, Fig. 3d shows a positive effect on early mornings during all weeks and a night effect during the weekend. Again, before the COVID-19 pandemic, people most needed AREU interventions during weeknights and weekends.

Fig. 3
figure 3

Plain model results: coefficients plots using data before the COVID-19 pandemic, a Effect of the hour of the day; b Effect of the day of the week, where \(1\) stands for Monday; c Effect of the quarters of the year, where \(1\) indicates the first quarter of the year; d Effect of the interaction between days and hours, where \(1\) stands for Monday. The black dotted lines represent two standard error bounds

Full size image

Finally, we report the quantiles at levels \(0.8\) and \(0.99\) of the predictive distribution considering hourly predictions one day ahead. Figure 4 shows these quantiles for each day from May, 9th, 2020 until May, 9th, 2021. Knowing the quantiles of predicted values can help AREU with resource planning and making informed emergency management decisions as a risk measure. For example, AREU knows that there is \(80\%\) probability that on March, 29th 2021, the hourly events will be at most 23 (and 29 at most with probability \(99\%\)), represented by the red dotted line in Fig. 4. In other words, \(23\) (\(29\)) represents the cut-off point such that only \(20\%\) (\(1\%\)) of the predicted emergency events exceed this value.

Fig. 4
figure 4

Plain model results: Quantiles at level \(0.8\) (black) and \(0.99\) (light green) considering the predicted events at hourly level one day ahead. The red dotted line shows the maximum quantile at level \(0.8\) (\(0.99\)) corresponding to \(23\) (\(29\)) events expected on March 29, 2021

Full size image

3.2 Model performance

We first focus on the model’s performance developed for the Plain SOREU. For each date from May, 9th, 2020 until May, 9th, 2021, the model described in the previous section was fitted on one year of retrospective data. The resulting model was used to forecast the number of events one, two, five, and seven days ahead. Figure 5a compares the predictions one, two, five, and seven days ahead with the true number of emergency events, while Fig. 5b shows the forecast errors for the one day ahead prediction. The periods where the error appears to be slightly larger are during the second pandemic wave (November-February 2021) and around the holiday season. In general, the mean absolute percentage error over the whole year was \(4.53\%\). This is a good result, considering that AREU required the mean absolute percentage error to be at most \(5\%\) for the predictions to be helpful in efficiently planning future activities. Appendix 2 shows the MAPE values grouped by month and year to have a better view of Fig. 5.

Fig. 5
figure 5

Plain model results: a Predictions one, two, five, and seven days ahead across one year and relative true values; b Forecast errors one day ahead, the predictions inside the two red dotted horizontal lines have absolute errors below \(5\%\). The black dotted line represents an absolute error equals \(0\)

Full size image

The model selected for the Plain SOREU data was then applied to the data of the other three SOREUs. The one day ahead mean absolute percentage error was still acceptable for the Alps (\(6.241\%\)), Lakes (\(5.815\%\)), and Metropolitan (\(4.309\%\)) SOREUs. Figure 6a shows the predictions one, two, five, and seven days ahead for the Metropolitan SOREU, while Fig. 6b describes the prediction errors for the one day ahead forecast. Similarly, Fig. 7a and b illustrate the results for the Lakes SOREU, and Fig. 8a and b for the Alps SOREU. As the Plain case, Appendix 2 reports the MAPE values considering all SOREUs grouped by month and year. The changes in the effects after fitting the model to pre and post-COVID data for the Metropolitan, Lakes and Alps SOREUs are reported in Appendix 1.

Fig. 6
figure 6

Metropolitan model results: a Predictions one, two, five, and seven days ahead across one year and relative true values; b Forecast errors one day ahead, the predictions inside the two red dotted horizontal lines have absolute errors below \(5\%\). The black dotted line represents an absolute error equals \(0\)

Full size image
Fig. 7
figure 7

Lakes model results: a Predictions one, two, five, and seven days ahead across one year and relative true values; b Forecast errors one day ahead, the predictions inside the two black dotted horizontal lines have absolute errors below \(5\%\). The black dotted line represents an absolute error equals \(0\)

Full size image
Fig. 8
figure 8

Alps model results: a Predictions one, two, five, and seven days ahead across one year and relative true values; b Forecast errors one day ahead, the predictions inside the two black dotted horizontal lines have absolute errors below \(5\%\). The black dotted line represents an absolute error equals \(0\)

Full size image

Figure 9 compares the one day ahead predictions obtained with the proposed model to the predictions of the model currently used by AREU on the Metropolitan SOREU, which was the only SOREU for which the AREU model was calibrated. The predictions are compared across \(145\) days (i.e., from 2020–10-14 until 2021–03-07). We observed a significant improvement: our model achieved a MAPE of \(4.44\%\) whereas the AREU model achieved a MAPE of \(5.11\%\). Therefore, the developed GAM showed a significantly lower error than the unobserved component time series model used by AREU.

Fig. 9
figure 9

Metropolitan model results: Predictions one day ahead across \(145\) days and relative true values

Full size image

Regarding the other benchmarking models, when looking at the one day ahead predictions over one year, the naive deterministic, ARIMA and model and tsglm models achieved a MAPE of \(7.359\%\), \(11.177\%\) and \(10.461\%\) respectively. Thus, the proposed model outperformed all of the considered benchmark models.

To test the flexibility of the proposed model, we trained it using data from the pre-COVID-19 period only, spanning from December 31, 2017, to December 31, 2018. Considering the Plain SOREU, we obtained MAPE values (calculated from December 31, 2017, to December 31, 2018): \(4.093\) for one day ahead, \(4.341\) for two days ahead, \(4.391\) for five days ahead and \(4.465\) for seven days ahead. The results for the other SOREUs are provided in Appendix 2. The proposed model demonstrates its ability to capture the dynamic changes in epidemiological and societal circumstances thanks to the incorporation of selected covariates such as the effective reproduction number (Rodpothong and Auewarakul 2012; Dabbaghian and Mago 2014)), as well as the adaptability of the splines in handling variations in the number of emergency events across different periods.

4 Discussion

We have described the development of a statistical model developed with state-of-the-art techniques to predict emergency events in the most populous Italian region. Because the forecasting errors were compatible with the organizational needs of AREU, the model is a valuable tool for efficiently planning emergency activities. A user-friendly dashboard providing the model estimates was implemented to make this model available to decision-makers in AREU, enabling the use of the model in everyday EMS coordination. Such an application is represented in Fig. 10.

Fig. 10
figure 10

Screenshot of one of the many user-friendly interfaces that make the proposed model fully operational on a daily basis to support the decisions by AREU

Full size image

Due to the different organization schemes of EMS around the globe, generalizing models predicting emergency events or calls to other contexts is often complex and, in some cases, impossible. However, the general methodology applied in this study and the overall model development process can also be used in other settings. Additionally, we found that our model performed very well on all of the SOREUs in Lombardy, even though the variables were selected based on the Plain SOREU data only. This suggests that the form of the model may generalize to slightly different contexts, such as other Italian regions or other European countries. This possibility is already under consideration by AREU decision-makers.

Our study focused on the GAM methodology for model development. Further research should evaluate the application of generalized additive mixed models (GAMMs), i.e., extensions of GAM incorporating random effects. GAMMs are better suited to deal with autocorrelation structures at a price, however, of higher computational costs (Lin and Zhang 1999). In addition, further analyses should explore the possible addition of interaction terms. Bayesian extensions could also be explored. Such models could be used to model trend and seasonal components with appropriate Markov random field priors with different forms and degrees of smoothness (Fahrmeir and Lang 2001). Finally, considering additional epidemiological variables, especially those related to COVID-19 or other emerging health crises, could be valuable and appropriate.

The proposed model was able to capture the dramatic daily and seasonal variations that emerged during the COVID-19 pandemic, attaining good performances both during one of the peaks of the COVID-19 epidemic in Lombardy and in moments when the volume of cases was limited. Nonetheless, we emphasize how predictive models like the one we presented need periodic, rigorous, quality evaluations and, eventually, substantial updates. This is particularly true in the COVID-19 era, when the rapidly increasing proportion of vaccinated individuals and the sudden spread of new viral variants may change the effect of COVID-19-related predictors in the model. Even beyond COVID-19 factors, structural and organizational modifications in the EMS apparatus may heavily impact the role of the predictors in the model. For these reasons, strict, continuous collaboration between public health stakeholders and data analysts is deemed essential.