We aim to discuss a robust model which can be useful to model and predict the spread of the Coronavirus disease SARS in Italy. In particular, we focus on a statistical model for daily deaths and cumulated intensive care unit hospitalizations and that can help to understand when the peak and the upper asymptote of contagion is reached, so that preventive measures (such as mobility restrictions) can be applied and/or relaxed. For these data robust procedures are particularly useful since they allow us to deal at the same time with model misspecifications and data reliability.
Nonlinear regression is an extension of classical linear regression, in which data are modeled by a function, which is a nonlinear combination of unknown parameters and depends on an independent variable. A relevant application of non-linear regression models concerns the modeling of so called dose-response relation, useful in toxicology, pharmacology and in the analysis of epidemic data. In these frameworks, the parameters of the model have a relevant interpretation, such as the upper limit and the inflection point.
Robust inference for nonlinear regression models 3 A normal non-linear regression model is obtained by replacing the linear predictor x T β by a known non-linear function µ(x, β), called the mean function. The model (see, e.g., Bates and Watts, 2007) y i = µ(x i , β) + ε i , with i = 1, . . . , n,
(1.1)
is called a non-linear regression model, where x i is a scalar covariate, β is an unknown p-dimensional parameter, and ε i are independent and identically distributed N (0, σ 2 ) random variables.
Likelihood inference is the usual approach to deal with nonlinear models. The log- Table 1 ). These models are parameterized using a unified structure with a coefficient b denoting the steepness of the curve, c and d the lower and upper asymptotes or limits of the response, and, for some models, e the inflection point. For instance, the five parameter log-logistic model assumes
However, in the presence of model misspecifications or deviations in the observed data, classical likelihood inference may be innacurate (see, e.g., Huber and Ronchetti, 2009, and references therein) . The aim of this paper is to discuss the use of robust inference on nonlinear regression models. In particular, we discuss a general approach based on the Tsallis score (Basu et al., 1988 , Ghosh and Basu, 2013 , Dawid et al., 2016 , Mameli et al., 2018 .
To deal with model misspecifications, useful surrogate likelihoods are given by proper scoring rules (see Dawid et al., 2016, and references therein) . A scoring rule is a loss function which is used to measure the quality of a given probability distribution for a random variable Y , in view of the result y of Y .
When working with a parametric model with probability density function f (y; θ), with θ ∈ Θ ⊆ IR d , an important example of proper scoring rules is the log-score, which corresponds to minus the log-likelihood function (Good, 1952) . In this paper, to deal with robustness, we focus on the Tsallis score (Tsallis, 1988) , given by
The density power divergence d α of Basu et al. (1998) is just (2.1), with γ = α + 1, multiplied by 1/α. The Tsallis score gives in general robust procedures Basu, 2013, Dawid et al., 2016) , and the parameter γ is a trade-off between efficiency and robustness.
Robust inference for nonlinear regression models
For the nonlinear regression model (1.1), the total Tsallis score for θ = (β, σ 2 ) is
(2.2)
The validity of inference about θ = (β, σ 2 ) using scoring rules can be justified by Let s(y; θ) be the gradient vector of S(y; θ) with respect to θ, i.e. s(y; θ) = ∂S(y; θ)/∂θ.
Under broad regularity conditions (see Mameli and Ventura, 2015 , and references therein), the scoring rule estimatorθ is the solution of the unbiased estimating equation s(θ) = n i=1 s(y i ; θ) = 0 (see Dawid et al., 2016) and it is asymptotically normal, with mean θ and covariance matrix V (θ)/n, where
are the sensitivity and the variability matrices, respectively. The matrix G(θ) = V (θ) −1 is known as the Godambe information and its form is due to the failure of the information identity since, in general, K(θ) = J(θ).
Asymptotic inference on the parameter θ can be based on the Wald-type statistic
which has an asymptotic chi-square distribution with d degrees of freedom; see Dawid et al. (2016) . In contrast, the asymptotic distribution of the scoring rule ratio statistic 
A consistent estimate of V (θ) can be obtained by using a parametric bootstrap; see
Varin et al. (2011) for a detailed discussion of the issues related to the estimation of J(θ) and K(θ). However, for Tsallis score (2.2) the matrices K(θ) and J(θ) can be derived analitically. Indeed, under the same assumptions of Theorem 3.1 in Gosh and Basu (2013), it is possibile to show that
and ξ α and ς α are the same as given in Gosh and Basu (2013) for the linear regression model (see Sect.6), namely ξ α = (2π) −α/2 σ −(α+2)/2 (1 + α) −3/2 and ς α = 1 4 (2π) −α/2 σ −(α+4)/2 2+α 2 (1+α) 5/2 . Moreover, the computation of J(θ) leads to
These matrices can then be used in (2.3) to derive the asymptotic distribution ofθ.
Note thatβ andσ 2 are asymptotically independent. For the Tsallis score (2.2), straightforward calculation show that the IF for the Tsallis estimator of the regression coefficients becomes
and that the IF for the Tsallis estimator of the error variance becomes
Since the functions s exp{−s 2 } and s 2 exp{−s 2 } are bounded in s ∈ IR, than both the influence functions IF (y;β) and IF (y;σ 2 ) are bounded in y for all γ > 1.
Paolo Girardi et al.
We illustrate statistical modelling of daily deaths ( We aim to build a data-driven model which can provide support to policy makers engaged in contrasting the spread of the Covid-19. In particular, we have applied robust inference for the model (1.1) to the available data, which cover the period from January 20 to April 4, 2020. The data sources are the daily report of the Protezione Civile (opendatadpc.maps.arcgis.com/ apps/opsdashboard/index.html).
The reader is also pointed to the work of the StatGroup-19 research group that models the mean function by using the five parameters Richards curve (Divino et al., 2020) .
Here, we consider two independent applications to:
• Daily Deaths (DD) for Covid-19, deaths confirmed by the Istituto Superiore di Robust inference for nonlinear regression models 9 Sanitá (ISS);
• Cumulative Intensive Care Unit (ICU) hospitalizations with positive Covid-19 swab; this data can be interpreted as a "department use index".
Here we model the considered data in two different geographical extensions: Italy and Lombardia. However, the proposed methodology has been applied to all the italian regions. We remark that in all the figures the models fit well the observed data. Moreover, it
can be also noted that with respect to the DD the peaks have been reached few day ago, while for the expected number of daily ICU we are on to the peaks. The peak is not a single day; it is a stabilization period that can last even several days, before a stable decrease.
The robust fits (Tsallis estimates and 95% confidence intervals) of the parameters e (inflection point) and d (upper asymptote) for the models are summarized in Tables   1 and 2 for DD and ICU, respectively. With respect to DD data the estimated 
In the simplest instance of prediction, the object of inference is a future or a yet unobserved random variable Z. Let p Z (z; θ) be the density of Z. The basic frequentist approach to prediction of Z, on the basis of the observed y from Y , consists in using the estimative predictive density function p e (z) = p Z (z;θ), obtained by substituting the unknown θ with a consistent estimator of θ, such as the Tsallis estimator or the MLE. Figure 5 reports the estimative predictive densities based on both the estimators for DD and cumulated ICU. Note the the Tsallis estimative predictive density is shifted on the right and exhibits a larger variability.
To compare the predictive performances of the Tsallis method with respect to the MLE, a simulation study has been performed, based on N = 10000 Monte Carlo replications. Figure 6 reports the boxplots of point estimators of the mean µ(z, β) 
To conclude, we believe that our procedure can constitute a useful statistical tool in modelling Italian Covid-19 contagion data. Indeed, the Tsallis robust procedures allow us to take into account for the inevitable inaccuracy of the Italian Covid-19 data, which are often underestimated. And example are those deaths of patients who died with symptoms compatible with Covid-19 but who have not had a tampon, or what has been described by many media regarding the growing number of elderly people who remain in their homes while needing to be hospitalized in intensive care. Thus, for these data robust procedures are particularly useful since the allow us to deal at the same time with model misspecifications (with respect to the normal assumption, independence or with respect to homoschedasticicy) and data reliability.
With respect to the italian Covid-19, we are now exploring the Bayesian approach (see Giummolé et al., 2019) , since it allow us to include prior information on the parameters of the model. Moreover, the plots of the posterior distributions for the parameters of the model and of the predictive distibutions may be quite useful in practice.
As a final remarks, since the variables are daily counts we will investigate the use of the Tsallis scoring rule in the context of nonlinear Poisson regression models.
Updates on the results and on the Italian regions may be found in the web page of the Robbayes-C19 research group:
https://homes.stat.unipd.it/lauraventura/content/ricerca.