key: cord-0573153-pan6p4e9
authors: Girardi, Paolo; Greco, Luca; Mameli, Valentina; Musio, Monica; Racugno, Walter; Ruli, Erlis; Ventura, Laura
title: Robust inference for nonlinear regression models from the Tsallis score: application to Covid-19 contagion in Italy
date: 2020-04-07
journal: nan
DOI: nan
sha: e7090f75b94ec838f6387245f521c35a7002dd00
doc_id: 573153
cord_uid: pan6p4e9

We discuss an approach for fitting robust nonlinear regression models, which can be employed to model and predict the contagion dynamics of the 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 nonlinear regression models.

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.

Robust and efficient estimation by minimising a density power divergence

Nonlinear regression analysis and its applications

Minimum scoring rule inference

Robust estimation for independent non-homogeneous observations using density power divergence with applications to linear regression

Robust Bayes estimation using the density power divergence

Objective Bayesian inference with proper scoring rules

Robust Statistics

Bootstrap adjustments of signed scoring rule root statistics

Higher-order asymptotics for scoring rules

Dose-Response Analysis using R

Possible generalization of Boltzmann-Gibbs statistics

This research work was partially supported by University of Padova (BIRD197903) and by PRIN 2015 (Grant 2015EASZFS003).