Robust inference for non‐linear regression models from the Tsallis score: Application to coronavirus disease 2019 contagion in Italy

We discuss an approach of robust fitting on non‐linear regression models, in both frequentist and Bayesian approaches, which can be employed to model and predict the contagion dynamics of the coronavirus disease 2019 (COVID‐19) in Italy. The focus is on the analysis of epidemic data using robust dose–response curves, but the functionality is applicable to arbitrary non‐linear regression models.


| INTRODUCTION
We aim to discuss a robust approach to model and predict the spread of the coronavirus disease 2019  in Italy, due to the worldwide epidemic outbreak of the new coronavirus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). In particular, we focus on deaths and intensive care unit (ICU) hospitalizations data that are expected to aid the detection of the time when the peaks and the upper asymptotes of contagion, in both daily new cases and total cases, are 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 with model misspecifications and data reliability simultaneously.
Non-linear regression is an extension of classical linear regression, in which data are modelled by a function, which is a non-linear combination of unknown parameters and depends on an independent variable. A relevant application of non-linear regression models concerns the modelling of so-called dose-response relation, useful in toxicology, pharmacology, and 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.
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 & Watts, 2007) 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 non-linear models. The log-likelihood function for θ = (β, σ 2 ) is and all likelihood quantities (maximum likelihood estimates, tests, confidence intervals, prediction, etc.) can be easily derived. Using the statistical environment R, the package drc (Ritz, Baty, Streibig, & Gerhard, 2015) provides a user-friendly interface to specify the model assumptions about the non-linear relationship and comes with a number of extractors for summarizing fitted models and carrying out inference on derived parameters.

| TSALLIS SCORE
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 a probability density function f(y;θ), with θ 2 Θ ⊆ R d , an important example of proper scoring rules is the log-score, which corresponds to minus the log-likelihood function. In this paper, to deal with robustness, we focus on the Tsallis score (Tsallis, 1988), given by Sðy; θÞ = ðγ −1Þ ð fðy; θÞ γ dy −γfðy; θÞ γ −1 , γ > 1: The density power divergence d α of Basu et al. (1998) is just (4), with γ = α + 1, multiplied by 1/α. The Tsallis score gives in general robust procedures (Ghosh & Basu, 2013), and the parameter γ is a trade-off between efficiency and robustness.

| Inference
The validity of inference about θ = (β, σ 2 ) using scoring rules can be justified by invoking the general theory of unbiased M-estimating functions.
Indeed, inference based on proper scoring rules is a special kind of M-estimation (see, e.g. Dawid et al., 2016, and references therein 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 & Ventura, 2015, and references therein), the scoring rule estimatorθ is the solution of the unbiased estimating equation sðθÞ = P n i = 1 sðy i ; θÞ = 0, and it is asymptotically normal, with mean θ and covariance matrix V(θ)/n, where 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. In contrast, the asymptotic distribution of the scoring rule ratio statistic W S ðθÞ = 2 SðθÞ−SðθÞ È É is a linear combination of independent chi-square random variables with coefficients related to the eigenvalues of the matrix J(θ)K(θ) −1 (Dawid et al., 2016). More formally, can be obtained by using a parametric bootstrap; see Varin, Reid, and Firth (2011) for a detailed discussion of the issues related to the estimation of J(θ) and K(θ). However, for Tsallis score (5), the matrices K(θ) and J(θ) can be derived analytically. Indeed, under the same assumptions of Theorem 3.1 in Ghosh and Basu (2013), it is possible to show that , n, and ξ α and ς α are the same as given in Ghosh 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 (6) to derive the asymptotic distribution ofθ. Note thatβ andσ 2 are asymptotically independent.

| Influence function
From the general theory of M-estimators, the influence function (IF) of the estimatorθ is given by and it measures the effect on the estimatorθ of an infinitesimal contamination at the point y, standardized by the mass of the contamination. The estimatorθ is B-robust if and only if s(y;θ) is bounded in y (see Huber & Ronchetti, 2009). Note that the IF of the MLE is proportional to the score function; therefore, in general, MLE has an unbounded IF; i.e. it is not B-robust. Sufficient conditions for the robustness of the Tsallis score are discussed in Basu et al. (1998) and Dawid et al. (2016).
For the Tsallis score (5), straightforward calculation shows 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 sexpf− s 2 g and s 2 expf−s 2 g are bounded in s 2 R, then both the influence functions IFðy;βÞ and IFðy;σ 2 Þ are bounded in y for all γ > 1.

| Bayesian inference
The use of surrogate likelihoods in the Bayes formula has received considerable attention in the last decade (see the review by Ventura & Racugno, 2016, and references therein).
Paralleling the derivation of posterior distributions based on composite likelihoods, a Tsallis posterior distribution can be obtained by using the scoring rule instead of the full likelihood in the Bayes formula. Let π(θ) be a prior distribution for the parameter θ. The proposed SR-posterior distribution is defined as with The choice of a prior distribution π(θ) to be used in (9) involves the same problems typical of the standard Bayesian perspective. For instance, for objective Bayesian inference, the expected α-divergence to the Tsallis posterior distribution can be used (Giummolé et al., 2019), and it is given by π G (θ) / jG(θ)j 1/2 . We aim to build a data-driven model that can provide support to policymakers engaged in contrasting the spread of the COVID-19. In particular, we have applied robust inference for Model (1) to the available data, which covers the period from February 24 to April 24, 2020. The data sources are the daily report of the Protezione Civile (https://github.com/pcm-dpc/COVID-19/). We consider two independent applications to:

| APPLICATION TO COVID-19 CONTAGION DYNAMICS IN ITALY
• Daily deaths (DD) for COVID-19, deaths confirmed by the Istituto Superiore di Sanità (ISS); • ICU hospitalizations with positive COVID-19 swab; these data can be interpreted as a 'department use index' .
Moreover, our analyses are limited to two different geographical extensions: Italy and Lombardia. However, the proposed methodology has been applied to all the Italian regions.
We illustrate statistical modelling of cumulative DD and cumulative intensive ICU in Italy and in the Italian region Lombardia using the Tsallis scoring rule. Following Dawid et al. (2016), for DD and ICU data, the Tsallis score (5) represents a composite scoring rule, based on only marginals.
In particular, it can be interpreted as an independent scoring rule (see also Varin et al., 2011, for composite likelihoods).
Figures 1 and 2 display the robust fitted models for Italy and Lombardia, respectively, for both DD and ICU data. Here, a three-parameter log-logistic curve has been used, i.e. obtained by setting c = 0, f = 1 in (3). In this setting, the parameter e represents the median time. In each plot, the blue points denote the observed data, and the red curves are the estimates/forecasts with our robust non-linear models. Furthermore, the plots on the left show the cumulated data, whereas those on the right give the daily data (new deaths and new hospitalizations, every day). By looking at the plots about the cumulative data, we appreciate the characteristic S-shape of the log-logistic curve describing the evolution of total cases. The upper asymptote of the curves represents the expected number of total deaths or ICU hospitalizations due to the COVID-19 outbreak.
On the other hand, the inspection of the panels devoted to the evolution of daily data is informative about the peak of daily new cases. The peak in the right panels corresponds to the mode of the distribution of cumulative cases, and it is always dominated by the parameter e due to the right 4 of 9 GIRARDI ET AL.
asymmetry. We remark that, in all the figures, the models fit well the observed data. In particular, the plots on the right show that the peak in new cases has already been reached for both DD and ICU data, hence highlighting the effectiveness of the restrictions. It is also evident that such counts grew faster at the beginning of the outbreak than they are decreasing after the peak. The plots on the left reveal that the end of the outbreak is still to come, and total numbers are expected to increase. When the cumulative data attain the upper asymptote, then the daily data decrease to zero.
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  counts, the total is about 76k occupancy-days, the inflection point is on April 10 (Day 47), and the peak on March 31. By the end of May, the death tolls will go below 15 units, whereas by the end of June, the number of ICU will be well below 100.
In order to assess the accuracy of the fitted models, some numerical studies have been carried out to investigate the actual sampling distribution of the proposed estimator. To this end, 10 000 samples have been drawn from the fitted model on Italy death data. The sampling distributions of the Tsallis estimates for (b, d, e, σ) are displayed in Figure 3. They all exhibit reasonable accuracy and precision, as confirmed by the comparison with the normal approximation to the distribution of the Tsallis estimator based on the fitted model.

| Estimative density and simulation study
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 F I G U R E 3 Sampling distribution of the Tsallis estimator for (b, d, e, σ). The solid curve is the normal approximation; the dashed curve is a kernel density estimate. The vertical dashed line gives the original estimate p e ðzÞ = p Z ðz;θÞ, obtained by substituting the unknown θ with a consistent estimator, such as the Tsallis estimator or the MLE. Figure 4 reports the estimative predictive densities based on both the estimators for DD and the cumulative ICU. Note the Tsallis estimative predictive density is shifted on the right and exhibits 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 = 10 000 Monte Carlo replications. Figure 5 reports the boxplots of point estimators of the mean μ(z, β) for the future observation z based on the predictive density estimated with the MLE and of the predictive density estimated with the Tsallis scoring rule, both under the central model (left) and under a contaminated model (right).
We note that, under the central model, the two estimators present a similar behaviour. On the contrary, under the contaminated model, only the Tsallis estimator presents a robust performance. The values of the bias (and sd) are quite similar under the central model, contrary to the contaminated model.

| Bayesian analysis
Mixing the robust procedure with the Bayesian approach allows us to include prior information (objective or subjective) on the parameters of the model. Moreover, the plots of the posterior distributions for the parameters of the model may be quite useful in practice since they are more informative than a simpler point or interval estimator. For instance, the marginal posterior distributions allow us to assign probabilities to intervals on the parameters. Figure 6 gives the violin plots of the marginal posterior distributions for parameters of the model and the expected mean, both for DD and for ICU data, using the reference prior (Giummolé et al., 2019). The posterior medians (the white dots) for the upper asymptotes and the inflection points are consistent with the frequentist robust estimates. Note that the classical marginal posterior distribution shows too narrow tails with respect to robust posterior distributions, which, on the contrary, can take into account the actual large uncertainty in the available data.

| FINAL REMARKS
To conclude, we believe that our procedure can constitute a useful statistical tool in modelling the Italian COVID-19 contagion data. Indeed, the Tsallis robust procedures allow us to take into account the inevitable inaccuracy of the Italian COVID-19 data, which are often underestimated.
Moreover, an example is those deaths of patients who died with symptoms compatible with COVID-19 but who have not had a symptom, 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 they allow us to deal at the same time with model misspecifications (with respect to the normal assumption, independence or homoscedasticity) and data reliability.
F I G U R E 4 Estimative predictive densities based on the Tsallis and the ML estimators for DD (left) and cumulated ICU (right) GIRARDI ET AL.

of 9
The estimation of the model and, as a consequence, the calculation of the expected asymptote and the inflection point are based on the assumptions that the adopted restrictions will not be subject to change. For these reasons, these fitted models cannot be used for predictive purposes, since it is not possible to predict how the data will change when the restrictions are modified at the end of the lockdown, scheduled for May 3. The day-by-day monitoring of the model fit stability will allow us to evaluate deviations from the actual lockdown situation with respect to the next reopening.
As a final remark, since the variables are daily counts, we will investigate the use of the Tsallis scoring rule in the context of non-linear Poisson regression models.