key: cord-350240-bmppif8g
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-08-12
journal: Stat (Int Stat Inst)
DOI: 10.1002/sta4.309
sha: 
doc_id: 350240
cord_uid: bmppif8g

We discuss an approach of robust fitting on nonlinear regression models, both in a frequentist and a Bayesian approach, which can be employed to model and predict the contagion dynamics of 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 approach to model and predict the spread of the COronaVIrus Disease 2019 in Italy, due to the worldwide epidemic outbreak of the new coronavirus SARS-CoV-2. In particular, we focus on deaths and intensive care unit hospitalizations data, that are expected to aid the detection of the time when the peaks and the upper asymptotes of contagion, both in 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.

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.

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) 

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-likelihood function for θ = (β, σ 2 ) is

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1002/sta4.309 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 et al., 2015) provides a user-friendly interface to specify the model assumptions about the nonlinear relationship and comes with a number of extractors for summarizing fitted models and carrying out inference on derived parameters.

A large number of more or less well-known mean functions are available (log-logistic, Weibull, Gamma, etc, see for instance Ritz et al., 2015 , Table 1 ). These functions 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 curve assumes

However, in the presence of model misspecifications or deviations in the observed data from the assumed model, classical likelihood inference may be inaccurate (see, e.g., Huber and Ronchetti, 2009 , Farcomeni and Ventura, 2012 , Farcomeni and Greco, 2015 , and references therein). This paper aims 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 , Giummolé et al., 2019 , both in a frequentist and Bayesian frameworks.

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 (4), with γ = α + 1, multiplied by 1/α. The Tsallis score gives in general robust procedures (Ghosh and Basu, 2013) , and the parameter γ is a trade-off between efficiency and robustness.

For the nonlinear regression model (1), the total Tsallis score for θ = (β, σ 2 ) is

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) . The class of M-estimators is broad and includes a variety of well-known estimators. For example it includes the maximum likelihood estimator (MLE) and robust estimators (see e.g. Huber and Ronchetti, 2009 ) among others.

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 and it is asymptotically normal, with mean θ and covariance matrix V(θ)/n, where

where K(θ) = E θ (∂s(Y; θ)/∂θ T ) and J(θ) = E θ (s(Y; θ)s(Y; θ) T ) 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 statis-

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, W S (θ)∼ d j=1 µ j Z 2 j , with µ 1 , . . . , µ d eigenvalues of J(θ)K(θ) −1 and Z 1 , . . . , Z d independent standard normal variables. Adjustments of the scoring rule ratio statistic have received consideration in Dawid et al. (2016) , extending results of Pace et al. (2011) for composite likelihoods. In particular, using the rescaling factor A(θ) = (s(θ) T J(θ)s(θ))/(s(θ) T K(θ)s(θ)), we have

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 (5) the matrices K(θ) and J(θ) can be derived analytically. Indeed, under the same assumptions of Theorem 3.1 in Gosh and Basu (2013) , it is possibile to show that

. . , n, 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 (6) to derive the asymptotic distribution ofθ. Note thatβ andσ 2 are asymptotically independent.

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, standardised by the mass of the contamination.

The estimatorθ is B-robust if and only if s(y; θ) is bounded in y (see Hampel et al., 1986) . Note that the IF of the MLE is proportional to the score function; therefore, in general, MLE has 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 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.

The use of surrogate likelihoods in the Bayes formula has received considerable attention in the last decade (see the review by Ventura and 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 Bayes formula. Let π(θ) be a prior distribution for the parameter θ. The proposed SR-posterior distribution is defined as

A possible choice of the matrix C is given by Giummolé et al. (2019) and references therein.

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 (θ) ∝ |G(θ)| 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 the model (1) to the available data, which covers the period from February 24th to April 24th, 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:

• Daily Deaths (DD) for COVID-19, deaths confirmed by the Istituto Superiore di Sanitá (ISS);

• Intensive Care Unit (ICU) hospitalizations with positive COVID-19 swab; this 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-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 nonlinear 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 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 intensive care unit 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 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 both for 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 for DD and ICU, respectively. For DD Italy data, the model predicts an impressive total of more than 31k deceased at the end of the outbreak. According to the fitted model, the expected number of new deaths will be below 15 counts by the end of June. The inflection point is on April 5th (day 42), whereas the fitted peak is on March 30th. For what concerns ICU, the total expected number of ICU occupancy-days is about 186k, whereas the inflection point is on April 8th (day 45) and the peak on April 1st. The number of ICU will decline at the same level as the end of February, at the very beginning of the outbreak in Italy, during the first days of July. Moving to Lombardia, the model estimates a total of about 148k deaths at the end of the epidemic, an inflection point on April 1st (day 38), and a peak on March 27th. For ICU counts, the total is about 76k occupancy-days, the inflection point is on April 10th (day 47), the peak on March 31st. 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, 10000 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 and precision, as confirmed by the comparison with the normal approximation to the distribution of the Tsallis estimator based on the fitted model. 

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 pe(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 = 10000 Monte Carlo replications. Figure 

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.

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 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 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 they allow us to deal at the same time with model misspecifications (with respect to the normal assumption, independence, or with respect to homoscedasticity) and data reliability.

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 lock-down, scheduled for May 3rd. The day by day monitoring of the model fit stability will allow us to evaluate deviations from the actual lock-down 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 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

In this web page also the R codes are available.

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) .

The data that support the findings of this study are available in GitHub repository at the Italian Protezione Civile account (pcm-dpc/ and described in Morettini et al. (2020) . These data were derived from the following resources available in the public domain:

https://github.com/pcm-dpc/COVID-19/.

Robust and efficient estimation by minimising a density power divergence

Nonlinear regression analysis and its applications

Minimum scoring rule inference

An overview of robust methods in medical research

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

COVID-19 in Italy: Dataset of the Italian Civil Protection Department

Dose-Response Analysis using R

Possible generalization of Boltzmann-Gibbs statistics

An overview of composite likelihood methods

Pseudo-likelihoods for Bayesian inference