At the end of January 2020, two Chinese tourists were hospitalized in Rome and tested positive to SARS-CoV-2. At the beginning of February, a group of Italian citizens was repatriated from Wuhan -among them, one tested positive. As the news media reported these headlines, neither the Italian public nor the Italian authorities appeared to perceive an imminent threat, though retrospective analyses now suggest that the virus may have been circulating in the north of the country as far back as December 2019 (e.g., detection of SARS-CoV-2 in the wastewater of Milan and Turin 1 ). The first recorded non-travel related COVID-19 case occurred in Codogno (Lombardia) -where a 38 years old male visited the hospital first on February 17, and then again on February 19 with worsening respiratory symptoms; in this date, he was tested and diagnosed. On February 20, two individuals tested positive in Vo' Euganeo (Veneto). Notably, the outbreaks in Lombardia and Veneto took two very different paths, something many observers attributed to the early response and aggressive testing strategy adopted by the regional authorities in Veneto 2, 3 . After some initial, much debated inconsistencies (e.g., hesitations in implementing local lock-downs in areas hosting major industrial production hubs, contested decisions to move patients between hospitals and nursing homes and to keep major sports events open to the public in Lombardia), starting in early March, local and central authorities took progressively more stringent measures to limit mobility and social gatherings -culminating with a general nationwide lock-down on March 9 and the suspension of all nonessential production activities on March 23 (starting in early May, activities restarted and mobility and gathering restrictions were gradually loosened).
Lock-down notwithstanding, based on official records, Italy saw a total of ≈ 35, 200 COVID-19 deaths as of the beginning of August. While other countries (e.g., the U.S. and Brazil) have reached much higher death counts, Italy's relative death toll remains rather stark at 58.25 per 100, 000 inhabitants. This may be partially attributable to the fact that Italy's population is very old (nationally, the median age is almost 46 years and the percentage of individuals over 65 almost 22%), and that age itself correlates with conditions such as type II diabetes, hypertension and chronic respiratory ailments, which substantially worsen illness and increase the likelihood of death for individuals affected by the virus. But perhaps the most striking aspect of the COVID-19 epidemic in Italy has been its heterogeneity. Some parts of Lombardia and of other regions in the industrialized north were hit early and especially hard, yet other demographically and socio-economically similar areas fared better. Moreover, most of the central and southern regions of the country experienced a much milder epidemic -notwithstanding waves of relocations from employment-related domiciles in the north back to family homes in the center and south around the time of the . ISTAT curves "take off" earlier and in some regions are as much as twice as high at their peak -possibly due to many COVID-19 deaths happening at home and/or not being recorded as such in hospitals, especially in the early stages of the epidemic. (b) MAX mortality curves (per 100,000 inhabitants) in the 20 Italian regions, before (top) and after (bottom) the shifts produced by probKMA run with K = 2. In the bottom panel, time is marked as a day number (as opposed to a date); this represents the region-specific time of the epidemic unfolding, and corresponds to actual time (starting on February 16 and ending on April 30) only for regions with no shifts, e.g., Lombardia). Curves are again smoothed with splines, with degree of smoothing selected by generalized cross-validation. Lombardia, Veneto, Emilia Romagna and Campania, also shown in (a), are highlighted in color. In all panels, vertical lines mark the dates of the national lock-down (March 9) and of the suspension of all nonessential production activities (March 23). In the bottom panel of (b) vertical lines still show these dates without shifts; stars on the curves mark the lock-down after the region specific shifts.
(the "exponential" pattern) followed by a slower descent from the peak; it includes many northern regions. Based on the shifts, Lombardia was first, followed by Emilia Romagna, Marche, Liguria, Piemonte, Trento/Bolzano, and last Valle d'Aosta. Lombardia and Valle d'Aosta presented the most extreme peaks -but Valle d'Aosta's descent was steeper (with a second late ascent likely due to data recording imprecisions; Valle d'Aosta is a very small region with only ≈ 125, 000 inhabitants). Group 2 follows a "flat(tened)" pattern; it includes all regions in southern and central Italy and, remarkably, Veneto -where the curve was successfully curbed. The shifts produced for this group are less stable and less meaningful in terms of interpretation, as flatter profiles leave more leeway in aligning curves against each other. All results (except for the shifts in Group 2) are rather consistent when using DPC and ISTAT curves (see Fig. S3a and Fig. S3c ), and when using distances defined on derivatives instead of curve levels. The solution with K = 3 places Lombardia (ISTAT curves) or Lombardia and Valle d'Aosta (MAX and DPC curves) in a cluster of their own (see Fig. S4 ). We also validated our results using a modification of funBI 16 , a functional biclustering technique, and IWTomics 17 , a functional testing technique which contrasts two sets of aligned curves pinpointing the locations and scales at which they differ (see Methods). Figure 2 (b) shows how, starting a little over two weeks from the beginning of their motif (wherever that was in each curve), Group 1 and Group 2 differ significantly at all temporal scales (see also Fig. S3b and Fig. S3d ). Why the two epidemics? The pattern of deaths characterizing Group 1 may be due, in large part, to the fact that the virus had circulated silently in the north of Italy for a long period of time before any kind of behavioral changes by the general public, medical protocols, or mitigation policies by local and central authorities were put in place. Mounting evidence suggests that a large share of COVID-19 cases are asymptomatic and yet contagious 3 Gray areas in the middle and bottom panels mark days when the difference between the two groups is significant (adjusted p-value < 5%). Starting a little over two weeks from the beginning of their epidemic, curves in the two groups differ significantly at all temporal scales.
reservoir of virus found its way to vulnerable individuals (some researchers also hypothesize Antibody-Dependent-Enhancement of SARS-CoV-2 18 , and thus a role for re-infections). But a variety of additional factors may have contributed to shaping the two epidemics; we explore some below.
Non-parametric FDA methods can be used to rank curves based on the notion of depth -from the innermost to the most extreme, and to identify outliers 19, 20 . Figure 3 shows a functional box plot of the MAX mortality curves and a depth ranking of the curves in the DPC, ISTAT and MAX data sets -shifted based on probKMA run with K = 2 and restricted to their aligned 65-day portions. The ranking is directional; we attributed signs to the depth measurements, so that curves far over or under the median curve are at the top or bottom of the ranking, respectively (see Methods). The top portion of the ranking comprises regions with "exponential" epidemics (Group 1) and is rather stable across data sets; Lombardia and Valle d'Aosta are consistently among the most extreme curves (they are also identified as outliers in the MAX and DPC data sets). The midand bottom portions of the ranking comprise regions with "flat(tened)" epidemics (Group 2) and are less stable across data sets, as the flatter profiles can more easily switch in their depth ranks. However, Toscana (which is the median in the MAX and ISTAT data sets) and Veneto are consistently among the deepest, most central curves. This analysis highlights again the tragic epidemic unfolding in Lombardia, and, by contrast, confirms how Veneto managed to "flatten" its curve back into the bulk.
Local mobility and positivity as statistical predictors of mortality Next, we focus on two key variables. The first is one of the most discussed policy-actionable variables, mobility, which has been curtailed to various degrees through lock-down measures in most of the countries affected by COVID-19. The second is Segments joining the regions across the three rankings show how the top portion remains rather stable, while the mid-and bottom portions contain several crossings. Regions at the top are those characterized by "exponential" epidemics (Group 1), while regions in the middle and at the bottom are those with "flat(tened)" epidemics (Group2), whose curves can more easily switch in their depth ranks.
one of the most discussed sentinel indicators, positivity, i.e. the fraction of performed tests returning positive results. For both these variables daily values for the period February 16 -April 30 can be obtained from data in the public domain at regional resolution. We considered differential mobility curves provided by Google for the category "Grocery & pharmacy". These express the fractional reduction with respect to January 2020 levels, and refer to mobility linked to first necessities -such as buying food, medicine, etc. For Italy, they are provided at the resolution of regions. Even though individuals were allowed to leave their homes for these necessities also during the most restrictive phase of the lock-down, the reduction captured by Google's "Grocery & pharmacy" was substantial. Mobility in weekdays fell by roughly 0.30, i.e. 30%, in the week after the lock-down Curves are smoothed with splines, with degree of smoothing selected by generalized cross-validation, and shifted based on probKMA run on the MAX mortality curves with K = 2; time is marked as a day number representing the region-specific time of the epidemic unfolding, and corresponds to actual time (starting on February 16 and ending on April 30) only for regions with no shifts, e.g., Lombardia. Vertical lines show the days corresponding to the nationwide lock-down (March 9) and the suspension of all nonessential production activities (March 23) without shifts, stars on the curves mark the lock-down after the region specific shifts. The example regions of Fig. 1 (a) are highlighted in color. (b) Estimated effect surfaces from the joint function-on-function regression of MAX mortality on local mobility and positivity shown in 3D and as contour plots (March 9, without shift, is again marked on both). Early and mid-period local mobility levels are strong positive predictors of mortality at its peak. Positivity has similar but much weaker predictive signals, likely because the effects are subsumed by mobility. Late local mobility has a negative association with mortality at its peak (mobility resumed faster in regions with milder epidemics), and late positivity a strong positive one (positivity remained elevated in regions with worse epidemics). The regression captures a large share of the variability in mortality curves (in-sample R 2 = 0.90, LOO-CV R 2 = 0.52), with substantial and comparable contributions of the two predictors (partial R 2 s = 0.62, 0.53).
(March 9), and further decreased in following weeks -reaching the lowest levels (between approximately −0.60 and −0.40 depending on the region) in the week after the suspension of nonessential production activities (March 23). It then slowly increased, getting back in a range between approximately −0.40 and −0.20 at the very end of April (see Fig. S5 ). In Lombardia, the peak MAX mortality was between March 20 and 25 -i.e., roughly, simultaneous to the lowest mobility and two weeks after its first substantial drop. Notably, in most Italian regions mobility during lock-down weekends reached −1.00, i.e. −100%. For comparison, in the state of New York, which had among the strongest restriction measures in the U.S., Google's "Grocery & pharmacy" never fell below −0.40. We refer to Google's "Grocery & pharmacy" curves as local mobility because they measure how much individuals move around where they live, as opposed to how much individuals move from place to place -e.g., to go from Wuhan to Milan, or from Milan to Palermo, or New York City. Obviously both types of mobility are relevant for the spread of a virus, and definitions depend on scale/resolution, but the first one is the one we analyzed.
To construct positivity curves, we combined daily public records on number of tests performed and number of new cases, which are also provided by the Italian Civil Protection agency. Taking daily ratios of new cases on tests performed is clearly imperfect, because of (variable and unreported) delays in test results. But regularizing and smoothing these ratios (see Methods)
produced a reasonable proxy. Smoothed positivity surpassed 0.1, i.e. 10%, as early as February 20 in some hard hit regions, peaked in a staggered fashion throughout March, and fell below 0.10 for all regions by around April 22 (see Fig. S5 ). Lombardia surpassed 0.10 around February 22 and peaked around ; that is, roughly, about a month and about a week prior to the peak of MAX mortality, respectively. Though we cannot draw exact parallels (our positivity curves are approximate and smoothed), this is consistent with what was observed, e.g., in New York City -where positivity was above 0.10 approximately from March 6-7 to May 12-13 and peaked at about 0.70 around March 28, with deaths peaking between April 5 and 13.
To anchor local mobility and positivity curves to the epidemic unfolding in each region, we shifted them congruently with the mortality curves. Figure 4a displays shifted curves based on probKMA run on MAX data with K = 2 ( Fig. S6 displays shifted curves based on probKMA run on DPC and ISTAT data). The horizontal axis now indicates again days in the region-specific epidemic unfolding, restricted to the 65-day portions where mortality curves align forming the two probKMA motifs.
We then used function-on-function regressions 8, 21 to model the statistical dependence of mortality on local mobility and positivity; in symbols, we fit the joint model y(t) = α(t) + β mob (s,t)x mob (s)ds + β pos (s,t)x pos (s)ds + ε(t), where y(t) is the response curve, i.e. mortality, α(t) is the intercept, ε(t) is the model error, and x mob (s) and x pos (s) are the predictor curvesmobility and positivity, respectively. These predictors are integrated over time, with "effects" represented by surfaces; β mob (t, s) is the association of mortality at time t with local mobility at time s, and similarly β pos (t, s) for positivity (see Methods). Figure 4b shows the effect surfaces for local mobility and positivity estimated using the MAX curves as response.β mob (t, s) suggests that local mobility levels early on and mid-way through the epidemic (e.g., around the March 9 lock-down date for Lombardia) are strong positive predictors of mortality at its peak, with the early predictive signal stronger than the mid-way one. In contrast, the local mobility level late in the epidemic has a negative association with mortality at its peak, likely reflecting a faster resumption of mobility in regions with milder epidemics.β pos (t, s) suggests that positivity levels early on and mid-way through the epidemic are also positive predictors of mortality at its peak -though the predictive signals are substantially weaker than those of mobility, likely because they are confounded with the latter. However, the positivity level late in the epidemic has a marked positive association with mortality at its peak. Here the signal is "detangled" from that of mobility, and one finds a sort of retrospective signature; regions which fared worse still had hightened positivity in the late stages of their epidemics. Estimated effect surfaces are remarkably similar across the three data sets (MAX, DPC and ISTAT), and the joint models all have in-sample R 2 s above 90% and leave-one-out cross-validated (LOO-CV) R 2 s above 50% (see Table S2 ), with strong and comparable contributions of local mobility and positivity (e.g., for the MAX curves, the partial R 2 s are 62% and 53%, respectively). Also, while this is not the case for all regions, residuals are rather consistent across data sets for Veneto, whose mortality is well predicted, and for Lombardia, whose mortality is always and sizably underestimated (see Fig. S7a and Fig. S10 ).
In order to further assess the roles of local mobility and positivity, we also considered marginal function-on-function regressions for mortality on each, separately; in symbols, y(t) = α(t) + β mob (s,t)x mob (s)ds + ε(t) and y(t) = α(t) + β pos (s,t)x pos (s)ds + ε(t). Effect surface estimates for local mobility are very similar to those in the joint models for all three data sets (see Fig. S8 ). Those for positivity confirm a strong association with mortality at its peak, but are less defined in terms of time profile (see Fig. S9 ). In summary, we find substantial evidence that local mobility and positivity are associated with COVID-19 mortality, and can predict it with some lag-time. Though the data at our disposal does not allow us to pinpoint lag lengths with accuracy, our analysis does support their roles as policy-actionable and monitoring variables, respectively. We also find that, even when considered jointly, these variables are not enough to fully account for the massive numbers of COVID-19 deaths recorded in Lombardia, the worst hit region in the country.
We considered several scalar (non-longitudinal) covariates proxying for socio-demographic, infrastructural and environmental factors debated by scientists and policy-makers during the epidemic (see Table S3 ). Many of these are suboptimal proxies retrieved from public sources. They refer to the closest times we could find data for (usually 2017 or 2018) and are, too, at the coarse resolution of regions. After an initial screen to guarantee data quality and limit the amount of collinearities, we focused on 12 covariates capturing aging of the population; prevalence of pre-existing conditions believed to affect disease severity; quality of distributed primary health care vs. centralized hospital-based health care; the potential of hospitals and nursing homes, but also schools, work places, households and public transport to act as contagion hubs; and pollution levels (see Table 1 ; Fig. S14 provides marginal densities, pair-wise scatter plots and correlations).
Even the restricted set of 12 covariates presents a distinct interdependence structure (see covariates dendrogram in Fig. 5 Table S4 ). For instance, our contagion hubs proxies for hospitals, schools and work places, and our (inverse) proxy for quality of distributed, primary health care (number of adults per family doctor), tend to vary closely together across regions. Also, our contagion hub proxy for public transport and pollution levels tend to vary together (this is not counter-intuitive, as both increase in more industrialized regions with large metropolitan areas), as do the percentages of individuals affected by diabetes and allergies, and our proxy for quality of centralized, hospital-based health care (ICU beds Ability of households to act as contagion hubs [4] 2017, ASR Lombardia Public transport rides per capita Ability of public transport to act as contagion hub 2017, ISTAT PM10
Pollution levels (particulates) 2018, ISTAT Table 1 . Scalar covariates potentially affecting COVID-19 mortality.
[1] The percentages of over 65, 70, 75 and 80 are highly correlated at the resolution of regions; we took over 65 as representative. [2] The prevalence of diabetes, hypertension and chronic bronchitis are highly correlated at the resolution of regions; we took diabetes as representative (allergies are not as highly collinear and were retained as separate). [3] Availability of ICU beds is also directly relevant for withstanding the impact of COVID-19 surges. [4] Average members per household is not a direct proxi of inter-generational contacts, but it may capture some of its effects.
per 100, 000 inhabitants) and the percentage of individuals over 65. Conversely, some regions show similar profiles across covariates (see regions dendrogram in Fig. 5(a) ). For instance, Lombardia, Veneto, Emilia Romagna and Piemonte have strong similarities, as do groups of southern regions (e.g., Sicilia, Campania, Puglia and Calabria; Basilicata, Abruzzo and Molise). An interesting characterization is produced using the Cheng and Church's biclustering algorithm 22 , which we implement with an adjusted mean squared residue, or H-score 23 . A bicluster is a subset of regions which exhibit similar behavior across a subset of covariates. Figure 5 (b) shows two biclusters with similar adjusted H-score values, obtained through the same run of the algorithm. The first bicluster comprises central and southern regions, all with "flat(tened)" epidemics (Group 2). Its regions have low ratios of adults to family doctors, limited concentrations in hospitals, nursing homes, work places and public transport, and low pollution levels. They also have high percentages of diabetic individuals and limited availability of ICU beds. The second bicluster comprises northern regions with "exponential" epidemics (Group 1), such as Lombardia, Emilia-Romagna and Piemonte, but also northern and central regions with "flat(tened)" epidemics (Group 2), such as Veneto, Friuli Venezia Giulia and Toscana. Its regions have high ratios of adults to family doctors, high concentrations in hospitals, work places and classrooms, and tend to have large percentages of individuals over 65. They also have low percentages of diabetic individuals and medium or small-sized households. Next, we used functional regressions with a two-fold aim: pursue a more direct, systematic assessment of the associations between the scalar covariates and COVID-19 mortality; and use the scalar covariates as controls in models comprising mobility and positivity to re-assess these key predictors. We stress again that the coarse resolution of the data poses serious limitations for these analyses, because it may dilute some predictive signals and because it bounds us to a small sample size. With only n = 20 observational units (the regions), fitting functional regression models comprising many terms (e.g., several scalar covariates and possibly their interactions; mobility and positivity curves along with more than one scalar covariate) produces unstable, overfit outcomes. Thus, we evaluate only the marginal effects of the scalar predictors, and the effects of mobility and stability with one scalar control at a time. The marginal function-on-scalar regressions of mortality curves on each of the 12 covariates have in-sample R 2 s ranging between ≈ 20 and 65%. Here the "effects" are curves; β x (t) represents the association of mortality at time t with the covariate x. For 8 of the covariates theβ x (t)'s show the expected signs throughout the peak period of the epidemic. In particular, the (inverse) proxy for quality of distributed, primary health care is the strongest marginal predictor; adults per family doctor shows a very large positive association with mortality. Also hospital, school and work place contagion hub proxies show strong positive associations with mortality. Nursing homes and public transport contagion hub proxies, pollution and the percentage of individuals over 65 are positive but comparatively weaker marginal predictors. For 4 of the covariates theβ x (t)'s show unexpected signs. The percentages of diabetics and allergic individuals show negative associations with mortality, likely due to the fact that their prevalence is high(er) in areas which were spared the brunt of the epidemic. In fact, estimated effect curves become positive when a differential intercept is included in the model to account for different overall mortality levels in Group 1 and Group 2 regions (see Fig. S12 ). Also the average number of members per household shows a negative association with mortality. Its small range of variation across regions (≈ 2.0 − 2. There are also marked associations among groups of covariates. The contagion hubs proxies for hospitals, schools and work places, and number of adults per family doctor, vary closely together. So do the contagion hub proxy for public transport and pollution levels; the percentages of individuals affected by diabetes and allergies; and ICU beds and the percentage of individuals over 65. (b): restricted heat-maps further illustrating interdependencies through two biclusters of regions and covariates. Color-coding within cells corresponds to that in (a), and each bicluster is identified by a border color and its adjusted H-score (an inverse measure of bicluster strength; see Methods). The first bicluster (adjusted H-score = 0.0902) comprises central and southern regions with "flat(tened)" epidemics (Group 2). The second bicluster (adjusted H-score = 0.0942) comprises northern regions with "exponential" epidemics (Group 1) but also northern and central regions from Group 2.
not allow it to properly proxy the effect of household contagions. At the same time, a strong negative correlation with the percentage of individuals over 65 may not allow it to properly proxy inter-generational contacts; regions with more elderly people are in fact those with smaller households. The negative association of average number of members per household with mortality, which persists even when including a differential intercept for Group 1 and Group 2 in the model (see Fig. S12 ), may simply be a "shadow" of this its negative correlation with the percentage of individuals over 65. Finally, ICU beds per 100, 000 inhabitants shows a positive association with mortality which, too, persists when including a differential intercept for Group 1 and Group 2 in the model (see Fig. S12 ), and may be in part a "shadow" of positive correlations with percentage of individuals over 65 and average number of beds per hospital. However, this proxy for quality of centralized, hospital-based health care, so prominent to the public debate during the epidemic, is not a negative predictor of mortality in our analysis. In conclusion, better proxies and finer resolution may reveal stronger aggravating roles for age, nursing homes, public transport and pollution 4, 5 and better dissect the roles of chronic conditions, households and inter-generational contacts, and ICU availability 24, 25 . But our analysis, notwithstanding limitations in the data, suggests important roles of primary care in mitigating mortality, and of 9/29 q 0 20 40 60 β(t) for marginal regressions
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Table 1 . The top plot displays the signs of the effect curves estimated on the MAX data. Time, marked as the 65 days of the region specific epidemic unfoldings, is on the vertical axis (the nationwide lock down on March 9, without shift, is marked by a horizontal line. Red, blue and green indicate, respectively, positive, negative and non-significant portions (i.e., where 95% confidence bands around the estimated effect curve are entirely above, entirely below, or contain 0; see Methods). The bottom plot displays in-sample R 2 s for the regressions fitted on MAX, ISTAT and DPC data; these are remarkably consistent. The names in red on the horizontal axes indicate the top 5 covariates selected by SnNAL-EN on all three data sets (see Methods); these are also the ones with the largest R 2 s. (b) Results from the joint function-on-function regression of MAX mortality on local mobility, positivity, and the first principal component (pc1) of the top 5 covariates, used as a "summary" control. This control does not modify the shapes of the estimated effect surfaces for mobility and positivity (shown on top) -which are very similar to the ones in Fig. 4(b) . The estimated effect curve for pc1 shows a positive and significant association with mortality at its peak (bottom right; 95% confidence band in dashes, gray corresponds to non significant portions, vertical dashed line corresponds to March 9, without shift). The sign of this effect is consistent with marginal findings, based on the loadings of the first principal component (bottom left; positive for adults per family doctor, average beds per hospital, average students per classroom and average employees per firm, and negative for average members per household). With the addition of pc1, the regression reaches an in-sample R 2 = 0.94 and a LOO-CV R 2 = 0.7. The contributions of local mobility and positivity remain high (partial R 2 = 0.66 and 0.61, respectively). That of our "summary" covariate is also substantial (partial R 2 = 0.39).
contacts in hospitals, schools and work places in aggravating it.
The results of our marginal function-on-scalar regressions, which are summarized in Fig. 6 (a) for MAX mortality curves, are also consistent across data sets (see Fig. S12 ) -which lends them support, at least at the resolution of regions. To further validate their stability we ran a functional generalization of SsNAL-EN 26 -an Elastic Net-like algorithm that performs feature selection for regressions with many predictors, producing reasonably stable outcomes even with small sample sizes and collinear features. Reassuringly, the output of SsNAL-EN is consistent with the marginal analysis, and again consistent across data sets (see Table S1 ): the top feature is always adults per family doctor, and the top 5 always include, in addition to it, average beds per hospital, average students per classroom, average employees per firm, and average members per household.
Finally, we ran again the function-on-function regression of mortality on local mobility and positivity, and re-evaluated the effects of these predictors introducing in the model one of the top 5 scalar covariates at a time (see Fig. S13 for results on 10/29 DPC, ISTAT and MAX data), as well as their first principal component, which explains ≈ 68% of their variability and can act as a "summary" control (see Fig. 6 (b) for MAX curves and Fig. S11 for DPC and ISTAT curves). Remarkably, while the control covariate "subsumes" some of the predictive power in each model, the estimated effect surfaces of local mobility and positivity retain the same shapes, and they remain very strong and comparable contributors (e.g., for MAX curves in Fig. 6(b) , the overall in-sample R 2 reaches 94%, the LOO-CV R 2 is 70%, and the partial R 2 s are 66, 61 and 39%, respectively, for local mobility, positivity and the first principal component; see also Table S2 ). Thus, with all the limitations of the data at our disposal, controlling for relevant covariates does not modify how the epidemic unfolding is associated to local mobility and positivity over time. Introducing socio-demographic, infrastructural and environmental factors in the modeling also does not change what we observed concerning residuals: mortality in Veneto is well predicted, and mortality in Lombardia remains sizably underestimated (see Fig. S7b ) for MAX and Fig. S11 for DPC and ISTAT).
Notwithstanding the limitations of the data employed in this study, using FDA techniques we were able to characterize heterogeneous and staggered epidemics in different areas of Italy --recapitulating and quantitating what scientists, policy makers and the public saw unfolding during the months of February, March and April 2020. In addition, we were able to document strong associations of COVID-19 mortality with local mobility and positivity, which persist in models that control for other relevant covariates. Investigating local mobility and positivity as, respectively, an actionable effector and a sentinel indicator of epidemic strength and progression, possibly to be used to adapt mitigation and containment efforts in real time, will require more and better data. In particular, accurate data on cases and hospitalizations in addition to deaths, and at a resolution much finer than that of Italian regions. Such data would also be critical to better capture predictive signals in a number of covariates -which may weaken and/or become confounded when aggregating data over broad, internally heterogeneous areas. But our results, along with those of other recent studies 27 , do support a role for mobility as a key modulator of COVID-19 spread and for positivity as a monitoring variable. Moreover, they support a role for distributed, primary health care in mitigating mortality, and for hospitals, schools and work places as contagion hubs that may aggravate the epidemic. If confirmed and fine-tuned on higher resolution data, also these findings could inform decision making -e.g., on short-and medium-term investments to boost distributed health care, or "pod" patients, students or employees. Finally, an extension of the temporal span of the data would also be of great interest to properly characterize different phases of the Italian epidemic -including its evolution after the gradual weakening of lock-down measures in May 2020. We believe that our work demonstrates the potential of FDA techniques for analyzing epidemiological data. Our pipelines and the mix of FDA tools used in this study could be applied to COVID-19 data from other parts of the world.
Daily cumulative COVID-19 death counts per region were retrieved from the Italian Civil Protection agency (Dipartimento della Protezione Civile; DPC A ). DPC mortality curves from February 24 to April 30 were computed for each region as the daily increments in COVID-19 death counts, divided by the population of the region as of January 1, 2019 (as recorded by ISTAT B ). DPC mortality curves were set to zero for the period February 16-23, before the Civil Protection agency started releasing data. Daily death counts from all causes in 7270 Italian municipalities (about 93.5% of the Italian population) for the years 2015-20 were downloaded from the Italian National Institute of Statistics (ISTAT C ) on June 4, 2020. Data were aggregated by region, and ISTAT differential mortality curves from February 16 to April 30 were computed for each region as the daily difference between 2020 deaths and the average daily deaths in 2015-19, divided by the total population of the municipalities included in the death counts as of January 1, 2019 D ). MAX mortality curves were created taking, for each region and each day, the maximum between DPC mortality and ISTAT differential mortality. Daily measurements concerning "Grocery & pharmacy" mobility from February 16 to April 30 were downloaded for each region from the Google Mobility Report E (local mobility curves). These measurements express percent changes with respect to the corresponding daily mobility levels in the first five weeks of 2020 (January 3 to February 6). Positivity curves were constructed using raw data from the Italian Civil Protection agency A . For each day from February 24 to April 30 and each region, we took the ratio between the number of new positive cases and the number of new tests performed. The ratios were truncated at 0 and 1 to account for irregularities in the row data A https://github.com/pcm-dpc/COVID-19/tree/master/dati-regioni B http://asti.istat.it/asti C https://www.istat.it/it/files/2020/03/Dataset-decessi-comunali-giornalieri-e-tracciato-record-4giugno.zip (e.g., positive cases = −1, or positive cases exceeding tests performed, presumably due to delays in test results). Like DPC mortality, positivity curves were set to zero for the period before the Civil Protection agency started releasing data (February [16] [17] [18] [19] [20] [21] [22] [23] . For all functional data sets, the two self-governing provinces of Trento and Bolzano were considered together as the Trento/Bolzano region, since not all data were available for both provinces separately. The 20 curves in each functional data set were smoothed using cubic smoothing B-splines with knots at each day and roughness penalty on the curve second derivative 8 . For each functional data set the smoothing parameter was selected minimizing the average generalized cross-validation error (GCV 28 ) across the 20 curves. All computations were performed using the R package fda 29 .
We considered a large number of scalar covariates of potential interest (see Table S3 ), and focused on the 12 listed in Table  1 and below. In retrieving and computing various measurements, as was done for the functional variables, the provinces of Trento and Bolzano were aggregated into the Trento/Bolzano region. % Over 65 was retrieved from ISTAT B at the regional level for the year 2018. % Diabetics and % Allergics were retrieved from ISTAT F at the regional level for the year 2018. Adults per family doctor was retrieved from the Ministry of Health G at the regional level for the year 2017. To compute ICU beds per 100, 000 inhabitants, we collected the total number of ICU beds in each region in 2018 from the Ministry of Health H , multiplied by 100,000 and divided by the population of the region as of January 1, 2019 B . To compute Ave. beds per hospital (whole) we used data from the Ministry of Health I , which provides the number of beds per ward in each hospital in 2018. We first aggregated them over wards belonging to the same hospital, and then averaged over hospitals in each region. Ave. beds per nursing home (ward) was also obtained based on data for the year 2018 from the Ministry of Health J -here we considered regional averages at the level of wards, without aggregating over wards inside the same nursing home (the ward-level covariate had a slightly higher association with mortality outcomes). To compute Ave. students per classroom we used data from the Ministry of Education K LM N , which provides the number of students in each classroom of each school in the country (public or private, at every level of education), for the year 2018. We averaged them over schools in each region. Data for Trento/Bolzano and Valle d'Aosta were missing, and were imputed through random forest imputation 30 , with default parameters maxiter=10 (maximum number of iterations to be performed given the stopping criterion is not met beforehand) and ntree=100 (number of trees to grow in each forest). To compute Ave. employees per firm we used data from ISTAT O , which provides number of employees per firm at the level of municipalities. We averaged them over firms in each region. Data for Valle d'Aosta were missing, and were again imputed through random forest imputation with default parameters. Ave. members per household was retrieved from ASR Lombardia P at the regional level for the year 2017. To compute Public transport rides per capita we used data from ISTAT Q , which provides the number of rides per capita for each Italian province in 2017. We multiplied these by the provinces' population as of January 1, 2019 B , summed up over provinces in the same region, and divided by the region population as of January 1, 2019 B . To compute PM10 we used data from ISTAT Q , which provides the average annual concentrations of PM10 (in µg/m3) detected by air quality meters distributed over the Italian territory. We averaged them over meters located in each region.
We used a number of standard multivariate techniques to analyze data on the 12 scalar covariates -including the extraction of Principal Components 31 , the calculation of Variance Inflation Factors 32 to evaluate multicollinearities, and clustering based on hierarchical agglomeration 31 . The latter was used both to agglomerate covariates with similar behavior across regions and to agglomerate regions with similar behavior across covariates. Agglomerative hierarchical clustering groups elements in a set with a bottom-up procedure that results in a dendrogram. Each element starts in its own cluster, and pairs of clusters are merged iteratively with a chosen distance for elements and linkage criterion for clusters. We employed the correlation distance, defined as d(x 1 , x 2 ) = 1 − corr(x 1 , x 2 ) for two generic elements x 1 and x 2 , and the complete linkage, defined as D(X 1 , X 2 ) = max x 1 ∈X 1 ,x 2 ∈X 2 d(x 1 , x 2 ) for two generic clusters X 1 and X 2 (thus, the distance between two clusters is defined as the furthest distance between their elements). We also used biclustering on the 20 (regions) by 12 (covariates) data matrix, to identify subsets of regions exhibiting similar behaviors across subsets of covariates. Following standard literature, we sought sub-matrices of the data whose entries are consistent with the "ideal" additive model x i, j = µ + α i + τ j , where µ is the typical value within the bicluster, and α i and τ j are additive adjustments for row i and column j, but we set all α i s to 0 in order to find constant column biclusters, i.e., sub-matrices with constant columns (covariates). We employed the Cheng and Church Biclustering Algorithm 22 , a greedy algorithm which finds the largest sub-matrices whose departure from the additive model is below a user-defined threshold. The departure is computed using the H-score (or mean squared residue score); in symbols, H(I, J) = 1 |I||J| ∑ i∈I, j∈J (x i, j − x I, j ) 2 , where I and J index the sets of rows and columns composing the bicluster, x i, j is a generic cell in the bicluster and x I, j is the mean of column j. We implemented this algorithm with a recently proposed adjustment to the H-score 23 that corrects a bias towards smaller biclusters in the original formulation. The adjusted H-score is defined as H ad j (I, J) = (∏ I−1 r=2 r 2
We performed local clustering of smoothed mortality curves (DPC, ISTAT and MAX, separately) using probabilistic K-mean with local alignment (probKMA 15 ). ProbKMA is a K-mean-like algorithm for functional data that finds K groups in a set of curves based on a local similarity among portions of the curves themselves. This allows the discovery of functional motifs, i.e. of typical local shapes that recur within and across the curves. In symbols, the algorithm finds K motifs v 1 , . . . , v K , membership probabilities p k,i and shifts s k,i (i.e. the starting points of the motif instances) for each cluster-curve pair that minimize the generalized least-squares functional J
, wherex i is the portion of the curve i corresponding to the shift s k,i , and d is the distance used to capture local similarity. For each data set, we considered K = 2 and K = 3. ProbKMA is probabilistic; it returns as output a membership probability p k,i for each cluster-curve pair. However, such an output can be turned into a hard partition by assigning each curve to the group with highest membership probability -which is what we did here. Notably, for K = 2, membership probabilities showed that Lombardia's and Valle d'Aosta's extreme mortality patterns were not well accommodated even in the "exponential" group 15 . The algorithm can employ different definitions of similarity d and thus capture different aspects of curve shapes. We used Euclidean (L 2 ) distance between curve levels for our main analysis -in symbols, d = 1 c c 0 (x(t) − v(t)) 2 dt for two generic curves x and v -though using Euclidean distance between curve derivatives produced similar results (not shown). ProbKMA allows the length of the motifs to be extended endogenously starting from a minimal one fixed in input. However, to identify epidemic patterns we ran it with a fixed motif length of 65 days -hence allowing a maximum shift of 10 days between curves (the mortality curves are 75 days long). The same clusters and very similar shifts were obtained with a fixed motif length of 50 days, which allows a maximum shift of 25 days (results not shown). The shifts produced by probKMA with K = 2 on the three mortality data sets (DPC, ISTAT and MAX) were employed to align, in addition to the mortality curves themselves, local mobility and positivity curves. All subsequent analyses employing shifted curves (tests contrasting groups of curves, functional boxplots and depth analyses, and functional regression models) were therefore restricted to the 65-day portions where mortality curves aligned following the two probKMA motifs. We also validated the groups produced by probKMA with a modified version of funBI 16 , an algorithm tipically used for finding functional biclusters. We used the modified funBI to identify groups of curves characterized by group-specific fixed length motifs, considering all possible sub-curves of a fixed length and clustering them with a divisive hierarchical algorithm (results not shown).
We employed an Interval-Wise Testing algorithm developed for omics data (IWTomics 17 ) to test for differences between the two groups of shifted mortality curves produced by probKMA with K = 2 (again, separately for DPC, ISTAT and MAX). IWTomics is a non-parametric, permutation-based functional hypothesis test. It contrasts two sets of curves aligned on a common domain to detect locations where the two sets differ significantly, and scales at which such significant differences are displayed (scales correspond to varying degrees of adjustment for multiple testing on intervals of varying lengths). Here locations are represented by the 65 days where the shifted mortality curves are defined, while scales vary from 1 day to the whole 65 days. The test was performed with the R package IWTomics 17 . The package allows the user to select among various possible test statistics; we employed the mean.
The functional boxplot 19 is an exploratory tool used to visualize functional data. It is constructed after ordering a set of curves based on a depth measure, such as the modified band depth 20 . The statistics employed to construct a functional boxplots are: the 50% central region envelope, the median curve, and the maximum non-outlying envelope. The 50% central region envelope corresponds to the box in a classical boxplot; it contains the 50% deepest, most centrally located curves. The median, i.e. the deepest curve, is inside this box and represents a robust "center" of the functional data set. The maximum non-outlying envelope is obtained by inflating the 50% central region envelope by 1.5 times its range. All curves extending outside of this envelope are flagged as outliers (the fact that the ISTAT data set in Fig. 3(b) lacks outlying curves based on this definition is due to the
width of its 50% central region envelope). We ranked the curves based on their depth measurements, after attributing a sign to such measurements with an ad hoc procedure. We subtract the median from each curve, and consider the share of the domain on which the difference is positive. If this is larger than 50%, we attribute a positive sign to the curve's depth -otherwise, we attribute a negative sign. Curves can thus be ranked from the most outlying above the median (labeled as positive), down to those close to the median, down to the most outlying below the median (labeled as negative) -see Fig. 3(b) . While this is not a fully general procedure, it works well on the DPC, ISTAT and MAX mortality curves we considered, which are rather unambiguously above/below the median (the share of the domain where the difference from the median is positive is ≥ 70 or ≤ 30% for all curves in all three data sets). Note also that the median curve of a data set, defined as the deepest, does not necessarily have half of the curves above it and half of the curves below it in the signed ranking we created (e.g., Toscana is the median curve in both ISTAT and MAX data sets, but the number of curves above/below it differs).
We consider models where a functional response variable is regressed against functional predictors and/or scalar covariates 8, 9 . All are special cases of the general equation 21 Fig. S12 we fit marginal regressions of this type allowing the estimation of two different intercepts: α 1 (t) for curves in Group 1 and α 2 (t) for curves in Group 2. Finally, for the joint regression of mortality on local mobility, positivity and one scalar control variable, we have L = 2 and J = 1. To fit all these functional regressions we used the R package refund 33 , which estimates the functional coefficients as well as their standard errors. We used these standard errors to construct confidence bands around the estimated functional coefficients. To gauge the explanatory power of each model, we computed the in-sample R 2 as well as the Leave-One-Out Cross-Validation (LOO-CV) R 2 . The former is a functional generalization of the classical coefficient of determination defined as SS reg /(SS reg + SS res ), where SS reg and SS res are the regression and the residual sum of squares, respectively. To compute the latter, for each observation i, one replaces the fitted response curveŷ i (t) (from the model fitted on all observations) with the predicted response curveŷ pred,i (t) obtained for i from the model fitted withholding i itself. Finally, for models with multiple terms (predictor and/or covariate), the partial R 2 of each term is computed as (R 2 − R 2 red )/(1 − R 2 red ), where R 2 is the coefficient of determination of the complete model, and R 2 red that of the model comprising all terms but the one being evaluated.
SsNAL-EN 26 is an efficient algorithm to perform Elastic Net 34 feature selection in a standard regression framework, i.e. when both response and features are scalars. The Elastic Net is a hybrid between LASSO and Ridge, which penalizes both the L 1 and the L 2 (Euclidean) norm of the regression coefficients. The L 1 penalty induces sparsity selecting only the most predictive among the features. The L 2 penalty regularizes coefficient estimates mitigating variance inflation due to collinearity. To perform feature selection in the functional regression setting, we applied a generalization of SsNAL-EN which incorporates a group structure in the Elastic Net objective function and uses the Functional Principal Components basis expansion to represent a functional response. In particular, we performed feature selection for the regression of mortality against all 12 scalar covariates in Table 1 . Notably, we selected the same top 5 features across all three data sets (DPC, ISTAT and MAX) (see Table S1 )lending strong support to their association with mortality. In the bottom panels vertical lines still show these dates without shifts, stars on the curves mark the lock-down after the region-specific shifts. Table S2 . Functional regression models (in-sample) R R R 2 2 2 , LOO-CV R R R 2 2 2 and partial R R R 2 2 2 s. For each functional linear model which regresses mortality on the covariates listed in the first column, the table reports the (in-sample) R 2 , the LOO-CV R 2 and the partial R 2 s. Figure S8 . Associating mortality to mobility. Results from the function-on-function regression of mortality on local mobility. The top row displays the estimated effect surface (the March 9 date is marked) with respective in-sample R 2 (for LOO-CV R 2 see Table S2 ). The bottom row displays the regression residuals (for barplots interpretation see Fig. S7 ). Figure S9 . Associating mortality to positivity. Results from the function-on-function regression of mortality on positivity. The top row displays the estimated effect surface (the March 9 date is marked) with respective in-sample R 2 (for LOO-CV R 2 see Table S2 ). The bottom row displays the regression residuals (for barplots interpretation see Fig. S7 ). Figure S10 . Associating mortality to mobility and positivity -ISTAT and DPC. Results from the joint function-on-function regression of ISTAT and DPC mortality on mobility and positivity. The top row shows the estimated effect surfaces (the March 9 date is marked) with respective partial R 2 (for in-sample R 2 and LOO-CV R 2 see Table S2 ). The bottom row shows the regression residuals (for barplots interpretation see Figure S11 . Associating mortality to mobility, positivity and first principal component -ISTAT and DPC. Results from the joint function-on-function regression of ISTAT and DPC mortality on mobility, positivity, and the first principal component (pc1) of the top 5 covariates, used as a "summary" control. Each panels shows: the estimated effect surfaces for mobility and positivity and the estimated effect curve for pc1 with respective partial R 2 (for in-sample R 2 and LOO-CV R 2 see Table S2 ). For interpretating the regression residuals see Fig. S7 .
MAX ISTAT DPC covariates R 2 LOO-CV R 2 partial R 2 s R 2 LOO-CV R 2 partial R 2 s R 2 LOO-CV R 2 partial R
marginal β β β (t q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 ISTAT % q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 DPC % q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q marginal β β β (t q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 ISTAT % q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 DPC % q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure S12 . Marginal function-on-scalar regressions. Results for marginal function-on-scalar regressions. Mortality curves are regressed against each of the scalar covariates in Table 1 . The top-row displays the signs of the effect curves estimated when just one intercept is included in the model. The bottom-row displays the signs of the effect curves estimated when we consider two different intercepts for curves in Group 1 and curves in Group 2. Time is on the vertical axis (the national lock down on March 9, without shift, is marked by a horizontal line). Red, blue and green indicate, respectively, positive, negative and non-significant portions (i.e., where 95% confidence bands around the estimated effect curve are entirely above, entirely below, or contain 0).
MAX ISTAT DPĈ β β β ( ( (t t t, , , s s s) ) ) mobilityβ β β ( ( (t t t, , , s s s) ) ) positivityβ β β ( ( (t t t, , , s s s) ) ) mobilityβ β β ( ( (t t t, , , s s s) ) ) positivityβ β β ( ( (t t t, , , s s s) ) ) mobilityβ β β ( ( (t t t 