The SARS-CoV-2 virus (also known as coronavirus) has emerged between the end of 2019 and the beginning of 2020 as a global threat. The disease caused by the virus has been designated COVID-19. It is hitting most countries in the world at a very quick rate, and with differential timing both across countries and within countries (see, e.g. [7] ). Such differential spreading of the epidemic over different regions provides the opportunity for optimizing the use of the health care resources across a country, thus potentially reducing mortality due to lack of sufficient care.
We have developed a simple software tool to allow for the quick monitoring of province-level data on the cumulative number of reported infected cases in Italy. In Section 2 we describe the methods, and in particular the data sources, the automatic model selection process, and the output provided by the software. In Section 3 we show a few sample plots for some Italian regions, as well as the graphical representation for all provinces. We close with some discussion in Section 4. * I wish to thank Marcello Pagano, Elena Savoia, Rino Bellocco and Alessia Melegaro for useful discussions and comments on this work. I also acknowledge Google for my use of their mapping API, which created the background map in Figure 2 .
Daily data on the reported numbers of infected individuals with the SARS-CoV-2 virus are available at the province level in Italy from the Department of Civil Protection [1] . The data report the daily cumulative number for each of the 107 provinces starting from February 24th, 2020, together with representative coordinates (longitude and latitude) for each province. The resident population of each province on January 1, 2019 was obtained from the Italian Statistical Institute [3] .
Given the resident population N (assumed constant) for a generic province, we let Y t be the number of reported infected at time t in that province. We model the observed rate (Y t /N ) via a Poisson generalized linear model with a first-, second-, or third-degree polynomial regression term in time and the natural logarithm link function [5] . Time is a natural number between one (for February 24, 2020) and the most recent day for which data are available.
Let us focus on the third-degree polynomial model. We thus assume Y t ∼ Poisson(λ t ), with
Note the inclusion of the offset term log(N ) to explicitly extract that (non-time-varying) exposure term from the constant parameter of the regression. The vector parameter β = (β 0 , β 1 , β 2 , β 3 ) T is estimated by maximum liklelihood in R [6] . We do not fit the model for provinces that have experienced a maximum reported infected rate smaller than 0.00005, and just output the "Prevalence too low" message for those cases.
The predicted number of reported infected cases at time t is thus
We can easily compute the first derivative of the rate Y t /N and estimate it by by plugging-in the maximum likelihood estimator (mle) β into the expression, and we obtain ∂ ∂t
Similarly, the second derivative (useful to identify flex points) is estimated by
Note that, by invariance, these estimators are the mles of the quantities that they are estimators for. We did not pursue it here, but their large-sample distribution can be obtained by the delta method from the approximate distribution of β. In particular, they are also approximately normally distributed, and a consistent estimator for the asymptotic variance can easily be constructed, if desired. While we have described the third-degree polynomial model above, the first-degree and the second-degree models follow immediately by dropping either both the quadratic and the cubic term, or just the cubic term from the third-degree model. The predicted rates for these models then follow immediately, and their estimated trends and curvatures can be easily shown to be equal to
for the second-degree model, and
for the first-degree model. Note that the expression of λ t should be adjusted by dropping the appropriate terms from Equation 2 above. Also, one can trivially obtain the expressions for the corresponding derivatives of Y t by dropping the N (constant) terms from the denominators. These formulas allow one to estimate the rate of change in the reported infected rate in particular at the most recent time point. Below we will exploit them to produce a map of the provinces that reflects such trends.
Note that the formulas for the first derivatives also allow for the immediate calculation of the time(s)t at which the local epidemic will reach a maximum (or minimum) for the second-degree model and for the third-degree model. These local maxima (or minima) are, for the two models,
when 4 β 2 2 −12 β 1 β 3 > 0. These predicted local or global maxima and minima are likely to be quite unstable until enough data are collected, and they should be interpreted with caution. As time goes by, one expects such estimated times to have increasing precision. More generally, due to the polynomial shape that we have chosen for the regression component, we expect these models to only be adequate to describe the epidemic during the initial growth phase for the first-degree model, and only up until the maximum prevalence of the epidemic for the second-degree and third-degree models. In particular, any predicted decreases in the cumulative number of cases is clearly meaningless since recovered cases are not accounted for in the data. In general, any forecasting beyond the last observed time point should be performed with caution.
The likelihood maximization algorithm might not converge, and that needs to be monitored and accounted for in the model selection process described below. Specifically, we initially try to fit the three nested models (first-, second-, and third-degree) to the data, and recover the convergence status for all three models. We then analyze the 8 possible cases separately. If only one of the models converges, we select that model. If two models converge, we consider the significance of the larger model through either the appropriate Wald test statistic (if the two models differ by just one term) or by the likelihood ratio test statistic if the two models that converge are the first-degree and the third-degree models. We then select the best model between the two. If all three models converge, we perform a backward elimination process starting from the third-degree polynomial model: if that is significantly better than the second-degree polynomial we select the former; if not, we drop the third-degree model and compare the second-degree model to the first-degree model, and select the final model between those. For all model selection tests we use a p-value cutoff equal to 0.05.
We first produce a separate plot for each province. In each plot we show the reported cases over time together with a line showing the fitted values from the selected model for that province. The plot also contains the indication of the degree of the polynomial that was selected, and the estimated rate of increase at the last time point, and the possible messages signaling that no modeling was performed due to the low prevalence in that province, or for lack of convergence of all models.
We then plot the estimated slopes on a map of Italy as obtained from google maps using the ggmap package ( [2] ; [4] ). For each province we plot an empty circle centered at the representative coordinates for that province, with the area of the circle proportional to the most recent fraction of reported positives over the resident population. Importantly, we use color-coded circles to indicate the value of the estimated growth rates as described in Subsection 2.1. We use the coding green=≤ 0.1 or N/A; yellow=(0.1, 0.2]; blue=(0.2, 0.3]; orange=(0.3, 0.4]; red=> 0.4. Lastly, the text output contains the observed rate of cases over the population; all the details about the models that were fit to each province; the estimated maxima / minima; and a final tables containing the names of all provinces and the estimated slope at the last time point, together with the class used to color the circle in the country map.
We ran our software on March 16, 2020, and data until March 16, 2020 were available. We show a sample of the plots that are generated. The complete output contains one plot for each province, one figure for the whole country, and a detailed text output (not shown here). Figure 1 shows a representative example of some province-specific plots, with the annotated date of last data point in the title, the degree of the selected model, and any warnings. Figure 2 shows the whole country, with colored circles for the provinces as described in Section 2.2. This figure clearly shows that the current stage of the epidemic consists of a highly prevalent region in the Padana valley (with large circles) but with lower rates of increase in the number of reported infected cases. The surrounding ares have lower prevalence, but faster growth. The center-South and the islands are still quite silent, although some provinces are starting to show the start of local epidemics. 
It should be noted that the data that is being analyzed here is just the number of reported positives, which clearly underestimates the number of positives in the population. As such, the use of the terms prevalence and incidence is not quite appropriate when analyzing these data. Also, we are not using the date of onset of the infection (which in most cases is unknown), so that any inference will suffer from a delay relative to the trend of the infections in the population (see, e.g. [8] ). On the other hand, we believe that these analyses still have the potential for capturing trends that may be useful for the monitoring and management of the disease burden on the health care system. Indeed, the software can be easily modified to be used on other counts which may be available, such as the number of the deceased persons or the number of patients in intensive care units.
One could envision more sophisticated models that take into account the time autocorrelation aspects, or the spatial correlation. While the former approach would be a reasonable choice, here we have focused on the behavior of the growth curves that the last available time point rather than on forecasting, and simplicity thus seemed to be a valuable criterion. In As far as the inclusion of spatial features, the spread of the epidemic likely follows transportation routes rather than pure distance among provinces, so that it would not be trivial to implement such an extended, global model over a country. In addition, the starting times of the local epidemics appear to differ. Indeed, as shown in Figure 1 , the choice of fitting separate third-degree polynomial regression models on the log scale estimated allows one to capture local epidemics that may (and do in our case) start at different times in different locations. The approach that we have followed here thus uses a common starting time for all provinces and automatically adjusts to different starting times for the local epidemics, while providing raw estimates of growth within each province on the same day. Note that any models used in such ongoing monitoring will need to be assessed closely as the days go by, so that they may continue to provide a satisfactory fit for all provinces.
Similarly to what was done for the local and global maxima, for the third-degree model one could estimate the time when the flex point was reached, or when it will be reached. From equation 4 this requires solving a fourth-degree polynomial, which for simplicity one could solve numerically. We did not pursue this here.
Perhaps most importantly, the software runs in a few seconds. Repeating the analyses daily allows for the production of a new map for each day, and the comparison of such maps over time allows for the appreciation of the spatial dynamics of the epidemic across the provinces.
The R code and the file that contains the resident population of the Italian provinces are available from https://sites.google.com/view/epilocal/home, or from the author. Comments and contributions are welcome.