# Robust confidence distributions from proper scoring rules

Erlis Rulia, Laura Venturaa and Monica Musiob a University of Padova, Italy,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
b University of Cagliari, Italy<EMAIL_ADDRESS>
###### Abstract

A confidence distribution is a distribution for a parameter of interest based
on a parametric statistical model. As such, it serves the same purpose for
frequentist statisticians as a posterior distribution for Bayesians, since it
allows to reach point estimates, to assess their precision, to set up tests
along with measures of evidence, to derive confidence intervals, comparing the
parameter of interest with other parameters from other studies, etc. A general
recipe for deriving confidence distributions is based on classical pivotal
quantities and their exact or approximate distributions.

However, in the presence of model misspecifications or outlying values in the
observed data, classical pivotal quantities, and thus confidence
distributions, may be innacurate. The aim of this paper is to discuss the
derivation and application of robust confidence distributions. In particular,
we discuss a general approach based on the Tsallis scoring rule in order to
compute a robust confidence distribution. Examples and simulation results are
discussed for some problems often encountered in practice, such as the two-
sample heteroschedastic comparison, the receiver operating characteristic
curves and regression models.

###### keywords:

AUC; Confidence density; $M$-estimators; Pivotal quantity; Regression model;
Robustness; Scoring rule; $t$-test; Tsallis Score

##  1 Introduction

Suppose data are analysed via some parametric model, and that $\psi$ is a
parameter of interest, such as a location parameter, the difference between
two means, the Area Under the Roc Curve (AUC), a regression coefficient, etc.
For conducting inference on $\psi$, statisticians have many methods in their
toolboxes, such as reaching point estimates, assessing their precision,
setting up tests along with measures of evidence, finding confidence
intervals, comparing $\psi$ with other parameters from other studies, etc. All
these inference topics may be automatically performed using an unique tool
when a frequentist distribution for $\psi$, given the observed data, is
available.

A practical approach on how to reach proper frequentist distributions, without
priors, is based on confidence distributions (CDs) and confidence curves
(CCs); see, among others, Xie and Singh (2013), Schweder and Hjort (2016),
Hjort and Schweder (2018), and references therein. In practice, a confidence
curve analysis is much more informative than providing the prototypical 95%
interval or a $p$-value for an associated hypothesis test. The plot in Figure
1 gives an illustration on making inference using a CC: point estimators
(mode, median and mean), 95% confidence interval and one-sided $p$-value.

Figure 1: Illustration of making inference on the scalar parameter of interest
$\psi$ using a confidence density.

The standard theory for parametric models evolves around the use of likelihood
methods, and this is also partly the case for the theory and applications of
CDs and CCs. Typically, to first-order, CD inference may be based on familiar
large-sample theory for the maximum likelihood estimator (MLE), the Wald
statistic and the likelihood-ratio test. The basic concepts and recipes for
CDs and CCs are however not limited to likelihoods per se, and various
alternatives may be worked with. For instance, it is well-known that for model
misspecifications or in the presence of deviant values in the observed data,
likelihood methods may be innacurate in many applications (see, e.g., Heritier
et al., 2009 and Farcomeni and Ventura, 2012). To deal with model
misspecifications, it may be preferable to base inference on procedures that
are more resistant, that is, which specifically take into account the fact
that the assumed models used by the analysts are only approximate. In order to
produce statistical procedures that are stable with respect to small changes
in the data or to small model departures, robust statistical methods can be
considered.

In this paper, we focus in particular on robust procedures based on proper
scoring rules (SRs). A scoring rule (see, for instance, the overviews by
Machete, 2013, and Dawid and Musio, 2014, and references therein) is a special
kind of loss function designed to measure the quality of a probability
distribution for a random variable, given its observed value. Proper scoring
rules supply unbiased estimating equations for any statistical model, which
can be chosen to increase robustness or for ease of computation. The Brier
score (Brier, 1950), the logarithmic score (Good, 1952), the Tsallis score
(Tsallis, 1988), and the Hyvärinen score (Hyvärinen, 2005) are well-known
instances of scoring rules. In particular, when using the logarithmic score,
the full likelihood and the composite likelihood (Varin et al., 2011) are
obtained as special cases of proper scoring rules (see for instance Dawid and
Musio, 2014). Frequentist scoring rule inference has been widely discussed,
also for robustness (see Ghosh and Basu, 2013, Dawid et al., 2016, and
references therein), while Bayesian inference based on scoring rules has been
considered in Dawid and Musio (2015), Ghosh and Basu (2016) and Giummolé et
al. (2019). To our knowledge, the only application of a particular scoring
rule in the context of CD inference is illustrated in Hjort and Schweder
(2018, Sect. 7). However, the associated scoring rule likelihood-ratio type
test has not a standard $\chi^{2}$ distribution, and thus is not an
appropriate pivotal quantity.

The aim of this paper is to discuss a general approach based on the Tsallis
scoring rule in order to compute a robust CD. In particular, we explore
asymptotic robust SR-CDs obtained by using pivotal quantites from SR inference
(see Dawid et al., 2016). Examples and simulation results are discussed for
three problems which statisticians often encounter in practice, such as the
two-sample heteroschedastic $t$-test, the area under the receiver operating
characteristic (AUC) curve and regression models.

The paper unfolds as follows. Sections 2 and 3 review, respectively, some
background on CDs and scoring rules. Section 4 discusses the construction of
the proposed SR-CDs and the derivation of the corresponding tail area
influence function. Examples and simulations studies are presented in Section
5. Finally, concluding remarks can be found in Section 6.

##  2 Background on confidence distributions

Consider a random sample $y=(y_{1},\ldots,y_{n})$ of size $n$ from a
parametric model with probability density function $f(y;\theta)$, indexed by a
$d-$dimensional parameter $\theta$. Write $\theta=(\psi,\lambda)$, where
$\psi$ is a scalar parameter for which inference is of interest and $\lambda$
represents the remaining $(d-1)$ nuisance parameters.

A modern definition of a confidence curve for $\psi$, say
$cc(\psi)=cc(\psi,y)$, can be found, among others, in Xie and Singh (2013) and
Schweder and Hjort (2016); see also references therein. Write $Y$ for the
random outcome of the data generating mechanism. At the true parameter point
$\theta_{0}=(\psi_{0},\lambda_{0})$, the random variable
$cc(\psi_{0})=cc(\psi_{0},Y)$ should have a uniform distribution on the unit
interval. Then

$\displaystyle P_{\theta_{0}}(cc(\psi_{0},Y)\leq\alpha)=\alpha,\quad\text{for
all}\,\,\alpha.$

Thus confidence intervals can be read off, at each desired level. When
$\alpha$ tends to zero the confidence interval tends to a single point, say
$\tilde{\psi}$, the zero-confidence level estimator of $\psi$.

In regular cases, $cc(\psi)$ is decreasing to the left of $\tilde{\psi}$ and
increasing to the right, in which case the confidence curve $cc(\psi)$ can be
uniquely linked to a full confidence distribution $C(\psi)=C(\psi,y)$, via

$\displaystyle
cc(\psi)=|1-2C(\psi,y)|=\left\\{\begin{array}[]{ll}1-2C(\psi,y),&\text{if}\,\,\psi\leq\tilde{\psi}\\\
2C(\psi,y)-1,&\text{if}\,\,\psi\geq\tilde{\psi}.\end{array}\right.$

With $C(\psi)$ a CD, $[C^{-1}(0.05),C^{-1}(0.95)]$ becomes an equi-tailed 90%
confidence interval, etc. Also, solving $cc(\psi)=0.90$ yields two cut-off
points for $\psi$, precisely those of a 90% confidence interval.
Correspondingly one may start with a given set of nested confidence intervals,
for all levels $\alpha$, and convert these into, precisely, a CD.

A general recipe to derive a CD is based on pivotal quantities. Suppose
$q(\psi;y)$ is a function monotone increasing in $\psi$, with a distribution
not depending on the underlying parameter, i.e. $q(\psi;y)$ is a pivotal
quantity. Thus $Q(x)=P_{\theta}(q(\psi;Y)\leq x)$ does not depend on $\theta$,
or on $\psi$, which implies that

$\displaystyle C(\psi)=Q(q(\psi;y))$

is a CD. The corresponding confidence curve for $\psi$ is

$\displaystyle cc(\psi)=\frac{\partial Q(q(\psi;y))}{\partial
q(\psi;y)}\,\frac{\partial q(\psi;y)}{\partial\psi}.$

###  2.1 Likelihood-based CDs

In various classical setups for parametric models, there are well-working
large-sample approximations for the behaviour of estimators, etc., and these
lead to constructions of CDs and CCs. For instance, if an estimator
$\hat{\psi}$ is such that $(\hat{\psi}-\psi)\,\,\dot{\sim}\,\,N(0,\tau^{2})$
to first-order, and $\hat{\tau}$ is a consistent estimator for $\tau$, then
$(\hat{\psi}-\psi)/\hat{\tau}\,\,\dot{\sim}\,\,N(0,1)$. Writing

$\displaystyle
C(\psi)\,\dot{=}\,\Phi\left(\frac{\psi-\hat{\psi}}{\hat{\tau}}\right),$

we have $C(\psi)\,\,\dot{\sim}\,\,U(0,1)$. Hence such $C(\psi)$ is an
asymptotically first-order valid CD, allowing us to write

$\displaystyle\psi|y\,\,\dot{\sim}\,\,N(\hat{\psi},\hat{\tau}^{2}),$

in the CD sense. In particular, if $\hat{\psi}$ is the MLE of $\psi$, then the
CD is derived from the profile Wald statistic

$\displaystyle
w_{p}(\psi)=\frac{\hat{\psi}-\psi}{\sqrt{j_{p}(\hat{\psi})^{-1}}},$ (2)

with $j_{p}(\psi)$ profile observed information, and it coincides with the
asymptotic first-order Bayesian posterior distribution for $\psi$.

A recipe that typically works better than (2) is the following. Let
$\ell(\theta)$ be the log-likelihood function for $\theta$, and let
$\ell_{p}(\psi)=\ell(\psi,\hat{\lambda}_{\psi})$ be the profile log-likelihood
for $\psi$, where $\hat{\lambda}_{\psi}$ is the MLE for $\lambda$ given
$\psi$. The profile log-likelihood ratio test
$W_{p}(\psi)=2(\ell_{p}(\hat{\psi})-\ell_{p}(\psi))$, under mild regularity
conditions, has an asymptotic null $\chi^{2}_{1}$ distribution. Hence
$\Gamma_{1}(W_{p}(\psi))\,\,\dot{\sim}\,\,U(0,1)$, with $\Gamma_{1}(\cdot)$
denoting the $\chi^{2}_{1}$ distribution function, and

$\displaystyle C(\psi)\,\dot{=}\,\Gamma_{1}(W_{p}(\psi))$

is an asymptotic CD. It can reflect asymmetry and also likelihood
multimodality in the underlying distributions, unlike the simpler Wald-type
confidence distribution. Similarly, the profile likelihood root

$\displaystyle
r_{p}(\psi)=\text{sign}(\hat{\psi}-\psi)\sqrt{2(\ell_{p}(\hat{\psi})-\ell_{p}(\psi))}$

can be used to derive a first-order CD, since it has a first-order standard
normal null distribution. Improved CD inference based on higher-order
asymptotics is discussed in Ruli and Ventura (2020).

##  3 Background on scoring rules

It is well-known that for complex models and/or model misspecification,
likelihood methods may be innaccurate in many applications. To deal with
complex models or model misspecifications, useful surrogate likelihoods are
given by proper scoring rules.

A scoring rule is a loss function which is used to measure the quality of a
given probability distribution $Q$ for a random variable $Y$, in view of the
result $y$ of $Y$; see Dawid (1986). The function $S(y;Q)$ takes values in
$\mathbb{R}$ and its expected value under $P$ will be denoted by $S(P;Q)$. The
scoring rule $S$ is called proper relative to the class of distributions
$\mathcal{P}$ if the following inequality is satisfied for all
$P,\,Q\in\mathcal{P}$:

$S(P;Q)\geq S(P;P).$

It is strictly proper relative to $\mathcal{P}$ if equation (3) is satisfied
with equality if and only if $Q=P$. Note that in the following we identify a
distribution $Q$ by its probability density $q$ with respect some measure
$\mu$; so the two notations $S(y;q)$ and $S(y;Q)$ are indistinguishable.

When working with a parametric model with probability density function
$f(y;\theta)$, an important example of proper scoring rules is the log-score,
which is defined as $S(y;\theta)=-\log{f(y;\theta)}$ (Good, 1952) and which
corresponds to minus the log-likelihood function.

In this paper we focus on the Tsallis score (Tsallis, 1988), given by

$S(y;\theta)=(\gamma-1)\int\\!f(y;\theta)^{\gamma}\,dy-\gamma
f(y;\theta)^{\gamma-1},\quad\gamma>1.$

The Tsallis score gives in general robust procedures (Ghosh and Basu, 2013,
Dawid et al., 2016), and the parameter $\gamma$ is a trade-off between
efficiency and robustness. Applications of the Tsallis score for robust
inference has been discussed in, among others, Ghosh and Basu (2013), Pak
(2014), Basu et al. (2016), Ghosh et al. (2019), and references therein.

###  3.1 Inference based on scoring rules

The validity of inference about $\theta$ 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), the maximum
composite likelihood estimator (see e.g. Varin et al., 2011), and robust
estimators (see e.g. Huber and Ronchetti, 2009 and references therein) among
others.

Given a proper scoring rule $S(y;\theta)$, let us denote by
$S(\theta)=\sum_{i=1}^{n}S(y_{i};\theta)$ the total empirical score. Moreover,
let $s(y;\theta)$ be the gradient vector of $S(y;\theta)$ with respect to
$\theta$, i.e. $s(y;\theta)=\partial S(y;\theta)/\partial\theta$. Under broad
regularity conditions (see Mameli and Ventura, 2015, and references therein),
the scoring rule estimator $\tilde{\theta}$ is the solution of the unbiased
estimating equation

$s(\theta)=\sum_{i=1}^{n}s(y_{i};\theta)=0$

and it is asymptotically normal, with mean $\theta$ and covariance matrix

$\displaystyle
V(\theta)=K(\theta)^{-1}J(\theta)(K(\theta)^{-1})^{{\footnotesize{{}^{{}_{\sf
T}}}}},$

where $K(\theta)=E_{\theta}(\partial
s(\theta)/\partial\theta^{{\footnotesize{{}^{{}_{\sf T}}}}})$ and
$J(\theta)=E_{\theta}(s(\theta)s(\theta)^{{\footnotesize{{}^{{}_{\sf T}}}}})$
are the sensitivity and the variability matrices, respectively. The matrix
$G(\theta)=V(\theta)^{-1}$ is known as the Godambe information and its form is
due to the failure of the information identity since, in general,
$K(\theta)\neq J(\theta)$.

Asymptotic inference on the parameter $\theta$ can be based on the Wald-type
statistic

$\displaystyle
w_{S}(\theta)=(\tilde{\theta}-\theta)^{{\footnotesize{{}^{{}_{\sf
T}}}}}V(\tilde{\theta})^{-1}(\tilde{\theta}-\theta),$

which has an asymptotic chi-square distribution with $d$ degrees of freedom.
In contrast, the asymptotic distribution of the scoring rule ratio statistic

$\displaystyle W_{S}(\theta)=2\left\\{S(\theta)-S(\tilde{\theta})\right\\}$

is a linear combination of independent chi-square random variables with
coefficients related to the eigenvalues of the matrix
$J(\theta)K(\theta)^{-1}$ (Dawid et al., 2016). More formally,

$W_{S}(\theta)\,\,\dot{\sim}\,\,\sum_{j=1}^{d}\mu_{j}Z_{j}^{2},$

with $\mu_{1},\ldots,\mu_{d}$ eigenvalues of $J(\theta)K(\theta)^{-1}$ and
$Z_{1},\ldots,Z_{d}$ independent standard normal variables. Adjustments of the
scoring rule ratio statistic have received consideration in Dawid et al.
(2016). In particular, using the rescaling factor
$A(\theta)=(s(\theta)^{{\footnotesize{{}^{{}_{\sf
T}}}}}J(\theta)s(\theta))/(s(\theta)^{{\footnotesize{{}^{{}_{\sf
T}}}}}K(\theta)s(\theta))$, we have

$\displaystyle
W_{S}^{adj}(\theta)=A(\theta)W_{S}(\theta)\,\dot{\sim}\,\chi^{2}_{d}.$

Analogous limiting results can be shown to hold for inference on the scalar
parameter $\psi$. With the partition $(\psi,\lambda)$, the scoring rule
estimating function is similarly partitioned as
$s(y;\theta)=(s_{\psi}(y;\theta),s_{\lambda}(y;\theta))$, where
$s_{\psi}(y;\theta)=(\partial/\partial\psi)S(y;\theta)$ and
$s_{\lambda}(y;\theta)=(\partial/\partial\lambda)S(y;\theta)$. Moreover,
consider the further partitions

$K=\left[\begin{array}[]{cc}K_{\psi\psi}&K_{\psi\lambda}\\\
K_{\lambda\psi}&K_{\lambda\lambda}\end{array}\right]\ ,\quad
K^{-1}=\left[\begin{array}[]{cc}K^{\psi\psi}&K^{\psi\lambda}\\\
K^{\lambda\psi}&K^{\lambda\lambda}\end{array}\right]\ ,$

and similarly for $G$ and $G^{-1}$. Finally, let $\tilde{\lambda}_{\psi}$ be
the constrained scoring rule estimate of $\lambda$, let
$\tilde{\theta}_{\psi}=(\psi,\tilde{\lambda}_{\psi})$, and let $\tilde{\psi}$
be the $\psi$ component of $\tilde{\theta}$.

A profile scoring rule Wald-type statistic for the $\psi$ component may be
defined as

$w_{Sp}(\psi)=(\tilde{\psi}-\psi)(\tilde{G}^{\psi\psi})^{-1/2},$

and it has an asymptotic $N(0,1)$ null distribution. Similarly, the profile
scoring rule score-type statistic
$s_{\psi}(\tilde{\theta}_{\psi})^{T}K^{\psi\psi}(G^{\psi\psi})^{-1}K^{\psi\psi}s_{\psi}(\tilde{\theta}_{\psi})$
has an asymptotic $\chi^{2}_{1}$ null distribution. Finally, we have that the
asymptotic distribution of the profile scoring rule ratio statistic for
$\psi$, given by
$W_{Sp}(\psi)=2\left(S(\tilde{\theta}_{\psi})-S(\tilde{\theta})\right)$, is
$\nu\chi^{2}_{1}$, where
$\nu=(\tilde{K}^{\psi\psi})^{-1}\tilde{G}^{\psi\psi}$. In view of this, an
adjusted profile scoring rule ratio statistic can be computed as

$\displaystyle
W_{Sp}^{adj}(\psi)=\frac{W_{Sp}(\psi)}{\nu}\,\,\dot{\sim}\,\,\chi^{2}_{1}.$

The adjusted profile scoring rule root, analogous to (2.1), can be defined as

$\displaystyle
r_{Sp}(\psi)=\text{sign}(\tilde{\psi}-\psi)\sqrt{W_{Sp}^{adj}(\psi)},$

which has an asymptotic standard normal distribution.

###  3.2 Examples

In this section, the Tsallis scoring rule is illustrated for two well-known
models: the regression model and the exponential family. All the quantities
necessary to compute $w_{Sp}(\psi)$ and $r_{Sp}(\psi)$ are derived.

Tsallis score for regression. Consider the general regression model

$\displaystyle
y_{i}=\mu(x_{i},\beta)+\varepsilon_{i},\quad\text{with}\,\,i=1,\ldots,n,$ (3)

with $x_{i}$ vector of fixed covariates, $\beta$ an unknown $p$-dimensional
parameter, and $\varepsilon_{i}$ independent and identically distributed
$N(0,\sigma^{2})$ random variables. The classical linear model is obtained
with $\mu(x_{i},\beta)=x_{i}^{{\footnotesize{{}^{{}_{\sf T}}}}}\beta$,
$i=1,\ldots,n$.

For this model, the total Tsallis score for $\theta=(\beta,\sigma^{2})$ is

$\displaystyle
S(\theta)=\sum_{i=1}^{n}\left[-\gamma\left(\frac{1}{\sqrt{2\pi\sigma^{2}}}\right)^{(\gamma-1)}\exp{\left(-\frac{(\gamma-1)}{2\sigma^{2}}(y_{i}-\mu_{i}(x_{i},\beta))^{2}\right)}+c_{\gamma}\right],$
(4)

with
$c_{\gamma}=(\gamma-1)/\sqrt{\gamma}(2\pi\sigma^{2})^{\frac{(\gamma-1)}{2}}$.
It is possibile to show that (Girardi et al., 2020)

$K(\theta)=\begin{pmatrix}\frac{\xi_{\alpha}}{n}\frac{\partial\mu}{\partial\beta}^{T}\frac{\partial\mu}{\partial\beta}&0\\\
0&\varsigma_{\alpha}\\\ \end{pmatrix},$

where
$\frac{\partial\mu}{\partial\beta}^{T}=(\frac{\partial\mu_{1}}{\partial\beta},\cdots,\frac{\partial\mu_{n}}{\partial\beta})$
is a $p\times n$ matrix, with $\mu_{i}=\mu(x_{i},\beta)$, $i=1,\ldots,n$, and
$\xi_{\alpha}$ and $\varsigma_{\alpha}$ are the same as given in Gosh and Basu
(2013) for the linear regression model (see Sect.6), namely
$\xi_{\alpha}=(2\pi)^{-\alpha/2}\sigma^{-(\alpha+2)/2}(1+\alpha)^{-3/2}$ and
$\varsigma_{\alpha}=\frac{1}{4}(2\pi)^{-\alpha/2}\sigma^{-(\alpha+4)/2}\frac{2+\alpha^{2}}{(1+\alpha)^{5/2}}$.
Moreover

$J(\theta)=\begin{pmatrix}\frac{\xi_{2\alpha}}{n}\frac{\partial\mu}{\partial\beta}^{T}\frac{\partial\mu}{\partial\beta}&0\\\
0&\varsigma_{2\alpha}-\frac{\alpha^{2}}{4}\xi_{\alpha}\\\ \end{pmatrix}.$

For the Tsallis score (4), it can be shown that the IF is bounded in $y$ for
all $\gamma>1$ (Girardi et al., 2020).

Tsallis score for exponential family. Let $Y$ have a distribution belonging to
the canonical exponential family

$f(y;\theta)=e^{\sum_{i=1}^{s}\theta_{i}t_{i}(y)-c(\theta)+d(y)},$ (5)

where $c(\theta)$ is a strictly convex $C^{\infty}$ function, $\theta$ is the
$s$-dimensional natural parameter, and $t(y)=(t_{1}(y),\ldots,t_{s}(y))$ is a
vector of sufficient statistics. It can be shown that

$\int\\!f(y;\theta)^{\gamma}\,dy=e^{c(\gamma\theta)-\gamma
c(\theta)}E_{\theta}\Big{[}e^{(\gamma-1)d(Y)}\Big{]}$

(see Nielsen and Nock, 2012).

The Tsallis score can be written as

$S(y;\theta)=(\gamma-1)e^{c(\gamma\theta)-\gamma
c(\theta)}E_{\theta}\Big{[}e^{(\gamma-1)d(Y)}\Big{]}-\gamma
f(y;\theta)^{(\gamma-1)}.$ (6)

In particular, when $d(y)$ is equal to $0$ (as in the Normal, Gamma and Beta
models for instance), (6) becomes a closed-form formula since
$E_{\theta}\Big{[}e^{(\gamma-1)d(y)}\Big{]}=E_{\theta}[1]=1$ and reduces to

$S(y;\theta)=(\gamma-1)e^{c(\gamma\theta)-\gamma c(\theta)}-\gamma
f(y;\theta)^{(\gamma-1)}.$ (7)

For instance, if $Y\sim\Gamma(\alpha,\beta)$, we have

$f(y;\alpha,\beta)=e^{-y\beta+(\alpha-1)\ln
y+\alpha\ln\beta-\ln\Gamma(\alpha)},$

with $\theta=(\alpha-1,-\beta)$ and

$c(\theta)=\ln\Gamma(\alpha)-\alpha\ln\beta=\ln\Gamma(\theta_{1}+1)-(\theta_{1}+1)\ln|-\theta_{2}|.$

Substituting (5) in (7), we find the following expression for the Tsallis
score

$S(y;\theta)=(\gamma-1)e^{c(\gamma\theta)-\gamma c(\theta)}-\gamma
e^{(\gamma-1)(\theta t(y)^{T}-c(\theta))}.$ (8)

Define $c_{i}(\theta)=\partial c(\theta)/\partial(\theta_{i})$ (note that
$c_{i}(\theta)=E_{\theta}(t_{i}(y))$), we have

$\displaystyle s(y;\theta_{i})$ $\displaystyle=$ $\displaystyle\frac{\partial
S(y;\theta)}{\partial\theta_{i}}=\gamma(\gamma-1)\Big{[}e^{c(\gamma\theta)-\gamma
c(\theta)}(c_{i}(\gamma\theta)-c_{i}(\theta)))+$ (9) $\displaystyle-$
$\displaystyle
e^{(\gamma-1)(\sum_{i=1}^{s}\theta_{i}t_{i}(y)-c(\theta))}\Big{(}t_{i}(y)-c_{i}(\theta)\Big{)}\Big{]}.$
(10)

The Tsallis estimator $\tilde{\theta}$ is the solution of the system

$\sum_{j=1}^{n}s(y_{j};\theta_{i})=0\quad i=1,\cdots,s.$

The estimator is robust if and only if $s(y;\theta_{i}),i=1,\ldots,s,$ is
bounded in $y$ for all $\theta$’s. In our case this condition requires that

$e^{(\gamma-1)(\sum_{i=1}^{s}\theta_{i}t_{i}(y)-c(\theta))}\Big{(}t_{i}(y)-E_{\theta}(t_{i}(y))\Big{)}$

is a bounded function of $y$ for each $i=1,\cdots,s$ and for each value of
$\theta=(\theta_{1},\ldots,\theta_{s})$. If $\gamma>1$, these conditions are
satisfied in the normal model, in the Gamma model if $\alpha>1$ and in the
Beta model if $\alpha$ and $\beta$ are both $>1$.

##  4 Confidence distributions from scoring rules

In this section we discuss how to derive CDs from proper scoring rules. In
particular, we discuss asymptotic CDs based on first-order approximations of
SR pivotal quantities.

Paralleling results in Section 2.1 for likelihood based CDs, a recipe to
derive an asymptotic CD from scoring rules is based on pivotal quantites. Let
$q_{S}(\psi;y)$ a scoring rule pivotal quantity, such as the profile Wald-type
statistic $w_{Sp}(\psi)$ or the adjusted profile scoring rule root
$r_{Sp}(\psi)$.

Thus,

$\displaystyle
C_{S}^{w}(\psi)\,\dot{=}\,\Phi\left((\psi-\tilde{\psi})(\tilde{G}^{\psi\psi})^{-1/2}\right)$
(11)

and

$\displaystyle
C_{S}^{r}(\psi)\,\dot{=}\,\Phi\left(\text{sign}(\psi-\tilde{\psi})\sqrt{W_{Sp}^{adj}(\psi)}\right)$
(12)

are first-order asymptotic CDs. The corresponding CCs are, respectively,

$cc_{S}^{w}(\psi)\,\dot{=}\,\frac{\phi\left((\psi-\tilde{\psi})(\tilde{G}^{\psi\psi})^{-1/2}\right)}{\sqrt{G^{\psi\psi}}}$

and

$cc_{S}^{r}(\psi)\,\dot{=}\,\phi\left(\text{sign}(\psi-\tilde{\psi})\sqrt{W_{Sp}^{adj}(\psi)}\right)\,\left|\frac{\partial
W_{Sp}^{adj}(\psi)^{1/2}}{\partial\psi}\right|,$

where $\phi(\cdot)$ is the density function of the standard normal
distribution. As for likelihood based CDs, (12) can reflect asymmetry, unlike
the simpler Wald-type confidence distribution (11).

For instance, using (12), the confidence median is $\tilde{\psi}$ and an
$(1-\alpha)$ equi-tailed confidence intervals can be obtained as
$\\{\psi:|r_{Sp}(\psi)|\leq z_{1-\alpha/2}\\}$. When testing, for instance,
$H_{0}:\psi=\psi_{0}$ against $H_{1}:\psi<\psi_{0}$, the $p$-value is
$p=C_{S}^{r}(\psi_{0})$, while when testing $H_{0}:\psi=\psi_{0}$ against
$H_{1}:\psi\neq\psi_{0}$ the $p$-value is $p=2(1-\Phi(|r_{Sp}(\psi_{0})|))$. A
measure of evidence for $\psi_{1}<\psi<\psi_{2}$ can be computed as
$C_{S}^{r}(\psi_{2})-C_{S}^{r}(\psi_{1})$.

Finally, note that (11) coincides with the asymptotic first-order Bayesian
posterior for $\psi$ discussed in Giummolé et al. (2019) and in Pauli et al.
(2011), and that, when in particular $S(\theta)$ is the logarithmic score,
(11) reduces to (2).

###  4.1 Robustness of the tail area

From the general theory of $M$-estimators, the influence function ($IF$) of
the estimator $\tilde{\theta}$ is given by

$\displaystyle IF(y;\tilde{\theta})=K(\theta)^{-1}s(y;\theta),$ (13)

and it measures the effect on the estimator $\tilde{\theta}$ of an
infinitesimal contamination at the point $y$, standardised by the mass of the
contamination. The estimator $\tilde{\theta}$ is B-robust if and only if
$s(y;\theta)$ is bounded in $y$. 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. The general theory of robust tests has been discussed in
Heritier and Ronchetti (1994).
Sufficient conditions for the robustness of the Tsallis score are discussed,
for instance, in Basu et al. (1998) and Dawid et al. (2016). In this section
we investigate the effects of model deviations on CDs. In this respect, let us
write the scoring rule pivotal quantity more generally as
$q_{S}(\psi;T(\hat{F}_{n}))$, where $\hat{F}_{n}$ is the empirical
distribution function and $T(F)$ is the functional defined by the scoring rule
estimating equation $\int s(y;T(F))\,dF(y)=0$, where $F=F(y;\theta)$ is the
assumed parametric model. Indeed, both the scoring rule pivotal quantities
$w_{Sp}(\psi)$ and $r_{Sp}(\psi)$ are functions of the scoring rule estimator.
For CD inference it is of interest the tail area. For a fixed value of $\psi$
the CD tail area is given by

$\displaystyle C_{S}(\psi)=\Phi(q_{S}(\psi;T(\hat{F}_{n})).$

To study the effects of deviations from the assumed parametric model $F$, we
use the tail area influece function (see, e.g., Field and Ronchetti,1990, and
Ronchetti and Ventura, 2001)

$\displaystyle
TAIF(y;T)=\left.\frac{\partial}{\partial\varepsilon}\Phi(q_{S}(\psi;T(F_{\varepsilon})))\right|_{\varepsilon=0},$
(14)

where $F_{\varepsilon}=(1-\varepsilon)F+\epsilon\Delta_{y}$ and $\Delta_{y}$
is the probability measure which puts mass 1 at the point $y$. The $TAIF(y;T)$
thus describes the normalized influence on the CD tail area of an
infinitesimal observation at $y$ and, by considering its supremum, it can be
used to evaluate the maximum bias of the tail area on the
$\varepsilon$-neighborhood of $F$.

By computing (14) we can identify the functional, defined by $s(y;T(F))$, for
which (4.1) is robust in the sense that the tail area influence function is
bounded and therefore the maximum bias of the corresponding $p$-value is
bounded in the $\varepsilon$-neighborhood of the model. After some
calculations, we obtain

$\displaystyle TAIF(y;T)$ $\displaystyle=$
$\displaystyle\phi(q_{S}(\psi;T(F)))\,\,\frac{\partial
q_{S}(\psi;T(F))}{\partial T(F)}\,\,\left.\frac{\partial
T(F_{\varepsilon})}{\partial\varepsilon}\right|_{\varepsilon=0},$ (15)

where the last term in (15) is the IF (13) of the scoring rule estimator.
Thus, the tail area influence function for the CD tail area at the statistical
model $F$ is proportional to the scoring rule estimating function and this
gives an immediate handle on robustness. Furthermore, it is bounded with
respect to $y$ when the scoring rule estimating function is bounded.

When considering robust scoring rules, i.e. scoring rules that lead to
estimators with bounded IF, such as the Tsallis scoring rule (see, e.g., Dawid
et al., 2016), then both (11) and (12) are robust CDs. On the contrary, when
considering the logarithmic score, typically (15) is not bounded since the MLE
is not B-robust.

##  5 Classical case studies in clinical research and simulation studies

In this section the practical usage of robust CDs will be presented in some
classical contexts, both by simulated and real-life data.

In particular, we consider the following frameworks: a regression analysis,
the comparison of two means in the presence of heteroschedasticity and
inference on the AUC.

###  5.1 Two sample comparison

Many experimental measurements are reported as realizations from a normal
distribution, and the simplest comparison we can make is between two groups.
The independent samples $t$-test is used when two separate sets of independent
and identically distributed samples are obtained, one from each of the two
populations being compared.

Let us assume that data are heteroscedastic, i.e. that standard deviations are
different from each other. In particular, let us assume that
$x=(x_{1},\ldots,x_{n_{x}})$ and that $y=(y_{1},\ldots,y_{n_{y}})$ are two
random samples from, respectively, $X\sim N(\mu_{x},\sigma^{2}_{x})$ and
$Y\sim N(\mu_{y},\sigma^{2}_{y})$. Writing
$\theta=(\mu_{x},\mu_{y},\sigma^{2}_{x},\sigma^{2}_{y})$, the Tsallis
empirical score is

$S(\theta)=(\gamma-1)c_{x}+\gamma\sum_{i=1}^{n_{x}}f_{x}(x_{i};\mu_{x},\sigma^{2}_{x})^{\gamma-1}+(\gamma-1)c_{y}+\gamma\sum_{i=1}^{n_{y}}f_{y}(y_{i};\mu_{y},\sigma^{2}_{y})^{\gamma-1},$

where $c_{x}$ and $c_{y}$ are suitable constants, depending on $\sigma^{2}x$
and $\sigma^{2}_{y}$, respectively. Typically, the parameter of interest is
$\psi=\mu_{x}-\mu_{y}$.

Simulation results. In order to assess the quality of CD inference for $\psi$
based on the Tsallis scoring rule in comparison to likelihood based CD
inference, we ran a simulation experiment in which the robustness constant
$\gamma$ is fixed in such a way that the resulting estimator is 10% less
efficient than the MLE, under the true model. We generated $10^{5}$ datasets
with sizes of the two samples $(n_{1}=10,n_{2}=20)$ form the true model with
parameter values $\mu_{x}=2,\mu_{y}=0,\sigma_{x}=\sigma_{x}=1$, i.e. $\psi=2$,
with and without contamination. Contaminated data were generated shiftwise, by
adding -7 to the last observation of the first sample, i.e. the sample with
size $n_{1}$. For each dataset we computed the coverage of Wald- and
$r_{p}$-type CD intervals at various confidence levels. Furthermore, the
uniformity of the $p$-values when testing $H_{0}:\psi=2$ against
$H_{1}:\psi<2$ is also checked.

From the simulation results shown in Figure 2 and 3 we note that, under the
central model, the two likelihood-based CDs show a reasonably good
performance. Apparently, the CD based on the Tsallis $r_{p}$ is slightly
better than its corresponding likelihood-based quantity. This sub-optimal
behaviour of $r_{p}$ under the true model could be explained by the fact that
MLE of $\sigma_{1},\sigma_{2}$ is biased in finite samples. Thus, we can
conclude that Tsallis, i.e. robust, CDs based on the $r_{p}$ statistic is a
valid alternative to its likelihood-based counterpart. Nevertheless, with
contaminated data, only robust $r_{p}$ CD gives improved inference, i.e. CIs
coverage closer to the nominal value, less biased estimators and $p$-values
closer to uniform.

Figure 2: Two normal random variables. Empirical coverage of CIs and
distribution of the CD estimates of $\psi$ with data generated under the true
model without contamination (first row) and with contamination (second row).
Dashed lines represent 10 $\times$ Monte Carlo standard error from the
theoretical confidence level.

Figure 3: Two normal random variables. Empirical distribution of $p$-values
for $H_{0}:\psi=2$ against $H_{1}:\psi<2$ using either non contaminated data
(first row) or contaminated (second row).

Case study. The debate whether statins can have adverse effects on cognitive
decline in elderly been raging since their introduction in 1987. The dataset
considered here (see Mandas et al., 2014) contains measurements on the Mini
Mental Score (MMSE) on 329 subjects aged 65 or older, living in a little
Sardinia village, collected in 2014. Two groups of subjects are considered: 59
cases treated with statins (average and sample standard deviation MMSE 24.76
and 3.47, respectively) and 270 controls (average and sample standard
deviation 26.05 and 3.41, respectively). The left plot of Figure 7 illustrates
the boxplots of the MMSE in the two groups. Due to the presence of several
outliers, the normal assumption for the control may be questionable. To
highlight the impact of the outlying observations, we consider a ”cleaned”
version of the dataset by eliminating three control patients with smallest
MMSE.

Figure 4: Left: boxplot of MMSE for case and control patients. Right:
confidence curves based the pivot $r_{p}(\psi)$ obtained from the likelihood
and the Tsallis scoring rule, using either complete (all) data or ”cleaned”
data.

In this application interest is on the group mean difference, i.e. for cases
minus controls and the aim is to test for the efficacy of the treatment by
testing $H_{0}:\psi=0$ against $H_{1}:\psi\neq 0$, or by computing a measure
of evidence for $\psi<0$. From the right plot of Figure 7 we notice that the
three outlying control observations have a substantial impact on the
likelihood-based CC. On the other hand, the CC based on the Tsallis scoring
rule with or without the outlying observations remains essentially unchanged.

Table 1 gives the $p$-value for testing $H_{0}:\psi=0$ against $H_{0}:\psi\neq
0$, the measure of evidence for $\psi<0$ and the median of the CCs, with the
complete and the cleaned datasets. Note that the CDs based on the complete
data and on the Tsallis scoring rule have similar summaries. On the other
hand, the CD based on the log-likelihood with the elimination of the outliers
gives quite different evidences.

Data | $r_{p}(\psi)$ |  $p$-value | Evidence for $\psi<0$ | Median
---|---|---|---|---
complete | log-likelihood | $5.7\times 10^{-3}$ | $2.9\times 10^{-3}$ | -1.30
cleaned | $1.5\times 10^{-3}$ | $7.7\times 10^{-4}$ | -1.47
complete | Tsallis | $5.5\times 10^{-4}$ | $2.8\times 10^{-4}$ | -1.51
cleaned | $4.9\times 10^{-4}$ | $2.5\times 10^{-4}$ | -1.52
Table 1: Summaries of CCs based on the MMSE data obtained from the
$r_{p}(\psi)$ pivot.

###  5.2 Area under the ROC curve

Let $X_{1}$ and $X_{2}$ be independent random variables with distributions
$F_{X_{1}}(x_{1};\theta_{1})$ and $F_{X_{2}}(x_{2};\theta_{2})$, respectively.
A stress-strength model is concerned with the problem of evaluating
$P(X_{1}<X_{2})$. For instance, in a clinical study, $X_{1}$ may be the
response of a control group, $X_{2}$ the response of a treatment group and the
reliability parameter $P(X_{1}<X_{2})$ measures the effectiveness of the
treatment as given by the area under the ROC curve.

By the definition of reliability, $P(X_{1}<X_{2})$ can be evaluated as a
function of the parameter $\theta=(\theta_{1},\theta_{2})$, through the
relation

$\psi=\psi(\theta)=P(X_{1}<X_{2})=\int
F_{X_{1}}(t;\theta_{1})\,dF_{X_{2}}(t;\theta_{2}).$

Theoretical expressions for $\psi$ are available under several distributional
assumptions both for $X_{1}$ and $X_{2}$ (see Kotz et al., 2003). For
instance, if $X_{1}$ and $X_{2}$ are independent normal random variables, i.e.
$X_{1}\sim N(\mu_{1},\sigma^{2}_{1})$ and $X_{2}\sim
N(\mu_{2},\sigma^{2}_{2})$, the reliability parameter is
$\psi=\Phi\left(\frac{\mu_{2}-\mu_{1}}{\sqrt{\sigma^{2}_{1}+\sigma^{2}_{2}}}\right)$.

Simulation results. Let $(x_{11},\ldots,x_{1n_{1}})$ and
$(x_{21},\ldots,x_{2n_{2}})$ be independent samples from two exponential
distributions with parameters $\lambda_{1}$ and $\lambda_{2}$. Since
$E(X_{1})=1/\lambda_{1}$ and $E(X_{2})=1/\lambda_{2}$, then
$\psi=\frac{\lambda_{1}}{\lambda_{1}+\lambda_{2}}$, which is the parameter of
interest.

In this situation the total Tsallis score is

$S(\lambda_{1},\lambda_{2})=\left(1-\frac{1}{\gamma}\right)\left(n_{1}\lambda_{1}^{\gamma-1}+n_{2}\lambda_{2}^{\gamma-1}\right)-\gamma\left[\lambda_{1}^{\gamma-1}\sum_{j=1}^{n_{1}}e^{-(\gamma-1)\lambda_{1}x_{1j}}+\lambda_{2}^{\gamma-1}\sum_{j=1}^{n_{2}}e^{-(\gamma-1)\lambda_{2}x_{2j}}\right].$

In the simulations, we set the true AUC to $\psi=0.85$ and the nuisance
parameter $\lambda_{2}=2/3$; thus the implied true rates for the exponential
distributions are $\lambda_{1}=3.778$ and $\lambda_{2}=2/3$. The robustness
tuning parameter $\gamma$ is fixed in such a way that the resulting estimator
is 10% less efficient than the MLE, under the true model. We generated
$10^{5}$ datasets with sizes of the two samples $(n_{1}=20,n_{2}=40)$ and
computed the coverage of Wald- and $r_{p}(\psi)$-type confidence intervals at
various confidence levels. Furthermore, the uniformity of the $p$-values when
testing $H_{0}:\psi=\psi_{0}$ against $H_{1}:\psi<\psi_{0}$ is also checked.
Data are generated from the true model, with and without contamination. To
generate a contaminated dataset, the last observation of the first sample,
i.e. the one with $n_{1}$ observations, is shifted by adding 3. The Wald-type
confidence intervals and the associated $p$-values were computed with $\psi$
reparametrized on the logit scale; the confidence intervals were then re-
transformed back to the original $[0,1]$ scale.

Figure 5: AUC for exponential random variables. Empirical coverage of CIs and
distribution of the CD estimates with data generated under the true model
without contamination (first row) and with contamination (second row). Dashed
lines represent 10 $\times$ Monte Carlo standard error from the theoretical
confidence level.

Figure 6: AUC for exponential random variables. Empirical distribution of
p-values for $H_{0}:\psi=0.85$ against $H_{1}:\psi<0.85$. Figures on the first
row are based on non contaminated data, those of the second row are based on
contaminated data.

From the simulation study, illustrated by Figures 5 and 6, we can conclude the
following. With non contaminated data, the performance of the various CD
considered is quite similar. Also the distribution of the $p$-values based on
$r_{p}(\psi)$-type pivots and Wald-type pivots (see Fig. 6) are very similar,
with the former being closer to uniform than those based on Wald-type pivots.
Under contaminated data, CDs based on $r_{p}(\psi)$ and Wald pivots perform
poorly, whereas their robust counterparts perform substantially better, with
the Tsallis $r_{p}(\psi)$ being preferred over the Tsallis Wald-type pivot.
The formers also lead to $p$-values being closer to uniform than the
$p$-values obtained with the likelihood-based pivots.

Case study. The data come from a study which aimed to assess the role of the
HSP70 (Heat Shock Protein 70 kilodaltons) protein on the presence of
Anaplastic Large Cell Lymphoma (ALCL, see Ventura and Racugno, 2011). Diseased
patients seem to have higher HSP70 levels than healthy subjects. Thus, HSP70
protein levels can be studied as a biomarker for detecting early ALCL lymphoma
and, therefore, its effectiveness in diagnosing the disease can be evaluated
by $\psi=P(X<Y)$. The data at hand consist of a small sample: 10 patients with
ALCL (sample mean and sample standard deviation equal 1.437 and 1.549) and 4
healthy subjects (sample mean and sample standard deviation equal 0.235 and
0.151). Two independent exponential random variables were assumed for the
protein level in both groups of patients. There appears (see the boxplots in
Figure 7) to be a diseased patient with HSP70 level very different from the
rest of the sample, thus in the subsequent analyses we will compute CDs for
$\psi$ using both the original dataset and the dataset obtained by deleting
the outlying diseased patient.

Figure 7: Left: boxplot of HSP70 for case and control patients. Right:
confidence curves based the pivot $r_{p}(\psi)$ obtained from the likelihood
and the Tsallis scoring rule, using either complete (all) data or with a
”cleaned” version.

From the CCs (shown on the right plot of Figure 7) we note that the CC based
on the Tsallis scoring rule is in between the two log-likelihood CCs, which
are obtained from the complete and cleaned data, respectively. Furthermore,
there also appear a substantial difference between the two Tsallis-based CCs,
which is presumably due to the small sample sizes. In particular, the 0.95 CI
for the AUC with the Tsallis CCs with complete and cleaned data are (0.596,
0.944) and (0.531 0.929), respectively, whereas those for log-likelihood based
are (0.627, 0.948) with the complete data and (0.54, 0.93) with the cleaned
data.

###  5.3 Linear regression model

Let us consider a linear regression model of the form (3) with
$\mu(x_{i},\beta)=x_{i}^{{\footnotesize{{}^{{}_{\sf T}}}}}\beta$,
$i=1,\ldots,n$. Usually, the parameter of interest is one of the regression
coefficients $\beta$.

The Tsallis score is given by (4). Th Tsallis score estimator is $B$-robust
since the influence function is bounded (see Ghosh and Basu, 2013).

Simulation results. For $p=3$, let $\psi=\beta_{2}$ be the scalar parameter of
interest and let $\lambda=(\beta_{1},\beta_{3},\sigma)$ be the nuisance
parameter. As in the previous examples, we ran a simulation experiment with
$n=50$ and the robustness constant $\gamma$ is fixed to 1.22.

Figure 8 (first column) reports the empirical coverages of bilateral
confidence intervals based on the four considered CDs, both under the central
model and under a contaminated model. The contamination is obtained with shift
contamination. As in the previous examples, under the central model, the
Tsallis CDs prove to be good competitors to the likelihood based CDs, while,
under the contaminated model, the robust CDs present a better and robust
performance. Figure 8 (second column) reports the boxplots of the median of
the likelihood and Tsallis based CDs, both under the central model and under a
contaminated model. Also in this example, under the central model, the two
estimators present a similar behaviour, while only the Tsallis estimator
present a robust performance with respect to contamination.

Figure 8: Regression model. Empirical coverage of CIs and distribution of the
CD estimates with data generated under the true model without contamination
(first row) and with contamination (second row). Dashed lines represent
10$\times$Monte Carlo standard error from the theoretical confidence level.

Finally, Figure 9report the uniform quantile-quantile plots of the p-values
from the four CDs when testing $H_{0}:\psi=\psi_{0}$ against
$H_{1}:\psi<\psi_{0}$, under the central model and under a contaminated model.
We note that, under the central model all the CDs present a reasonable
performance, while, under the contaminated model, only the robust CD presents
a good performance.

Figure 9: Regression model. Empirical distribution of $p$-values for
$H_{0}:\psi=0$ against $H_{1}:\psi<0$. Figures on the first row are based on
non contaminated data, those of the second row are based on contaminated data.

Case study. The GFR dataset contains measurements of the glomerular filtration
rate ($GFR$) and serum creatinine ($CR$) on $n=30$ subjects. The $GFR$ is the
volume of fluid filtered from the renal glomerular capillaries into the
Bowmans capsule per unit of time (typically in millilitres per minute) and,
clinically, it is often used to determine renal function. Its estimation is of
clinical importance and several techniques are used for that purpose. One of
them is based on $CR$, an endogenous molecule, synthesized in the body, which
is freely filtered by the glomerulus (but also secreted by the renal tubules
in very small amounts). Several models have been proposed in the literature to
explain $GFR$ as a function of $CR$. Here, following Heritier et al. (2009),
we consider a model for $GFR$ based on $CR^{-1}$ and $AGE$, i.e.
$GFR=\beta_{1}+\beta_{2}CR^{-1}+\beta_{3}AGE+\varepsilon$. The data are
illustrated in Fig. 10: note that there are some observations which look like
outliers.

Figure 10: Scatterplot diagram of GFR data.

Figure 11 gives the CCs based on $r_{p}(\psi)$ and $w_{p}(\psi)$ and on
$r_{Sp}(\psi)$ and $w_{Sp}(\psi)$ for $\psi=\beta_{3}$, i.e. the parameter of
$AGE$. It can be noted that the robust CCs are quite different, giving
different inferential conclusions about the effect on $AGE$ on $GFR$.

Figure 11: CDs based on $r_{p}(\psi)$ and $w_{p}(\psi)$ and on $r_{Sp}(\psi)$
and $w_{Sp}(\psi)$ for $\psi=\beta_{3}$.

Table gives the p-values for $H_{0}:\psi=0$ against $H_{0}:\psi\neq 0$, the
measure of evidence for $\psi<0$ and the median of the CCs, based on
$r_{p}(\psi)$ and on $r_{Sp}(\psi)$. Note that the CDs based on $r_{p}(\psi)$
and on $r_{Sp}(\psi)$ give quite different conclusions about the effect of
$AGE$ on $GFR$.

CD | $CC(0)$ | p-value $H_{0}:\psi=0$ | median
---|---|---|---
$\Phi(r_{p}(\psi))$ | 0.85 | 0.298 | -0.33
$\Phi(r_{Sp}(\psi))$ | 0.98 | 0.039 | -0.50
Table 2: CDs summaries based on $r_{p}(\psi)$ and on $r_{Sp}(\psi)$ for the
GFR data and $\psi=\beta_{3}$.

##  6 Discussion

In practical applications, CDs are more informative than a simpler confidence
interval or a $p$-value, since they describe the complete distribution
estimator for the parameter of interest, as the posterior distribution for
bayesians. We would like to stress that, under appropriate smoothness
conditions, the Tsallis scoring rule can be applied to any statistical model
$f(y;\theta)$ and delivers an associated $M$-estimator. While this may lead to
a loss of efficiency in comparison with full likelihood methods, it can
exhibit improved robustness or computational advantages. Moreover, under
smoothness conditions, any proper scoring rule can be used to derive a CD,
using the first-order approximations of SR pivotal quantities.

Higher-order asymptotic expansions for scoring rules have been recently
discussed by Mameli and Ventura (2015) and Mameli et al. (2017). These higher-
order asymptotic expansions to the distribution of the scoring rule estimator,
of the scoring rule ratio test statistic and of the signed scoring rule root
statistic for a scalar parameter allow to derive higher-order pivotal
quantities, which improves the first-order approximations. The use of these
higher-order expansions could be investigated to derive CDs when dealing with
small sample sizes.

Finally, in this paper only the median of the CDs has been considered. We are
planning to investigate also other point estimators derived from the CDs, such
as the mode, in particular in situation in which the CD exhibits a strong
asymmetry.

#### Founding

This research work was partially supported by the University of Padova
(BIRD197903) and by MIUR (PRIN 2015, grant 2015EASZFS_003).

## References

  * [1] Basu A, Harris IR, Hjort NL, Jones MC. Robust and efficient estimation by minimising a density power divergence. Biometrika, 1998; 85: 549–559.
  * [2] Basu A, Mandal A, Martin N, Pardo L. Generalized Wald-type tests based on minimum density power divergence estimators. Statistics 2016; 50: 1–26.
  * [3] Brier GW (1950). Verification of forecasts expressed in terms of probability. Mon. Weather Rev., 1950; 78: 1–3.
  * [4] Dawid AP. Probability forecasting. In: Encyclopedia of Statistical Sciences (S. Kotz, N. L. Johnson, and C. B. Read eds.) 1986; 210–218.
  * [5] Dawid AP. Musio M. Theory and Applications of Proper Scoring Rules. Metron, 2014; 72: 169–183.
  * [6] Dawid AP, Musio M. Bayesian model selection based on proper scoring rules (with discussion). Bayesian Analysis, 2015; 10: 479–521.
  * [7] Dawid AP, Musio M, Ventura L. Minimum scoring rule inference. Scand. J. Statist., 2016; 43: 123–138.
  * [8] Farcomeni A, Ventura L (2012). An overview of robust methods in medical research. Stat. Meth. Med. Res., 2012; 21: 111–133.
  * [9] Field CA, Ronchetti E. Small Sample Asymptotics. IMS Monograph Series, Hayward (CA); 1991.
  * [10] Ghosh M, Basu A. Robust estimation for independent non-homogeneous observations using density power divergence with applications to linear regression. Electr. J. Statist., 2013; 7: 2420–2456.
  * [11] Ghosh M, Basu A. Robust Bayes estimation using the density power divergence. Ann. Inst. Stat. Math., 2016; 68: 413–437.
  * [12] Ghosh A, Martin N, Basu A, Pardo L. A new class of robust two-sample Wald-type tests. Int. J. Biostat., 2019; 20170023.
  * [13] Girardi P, Greco L, Mameli V, Musio M, Racugno W, Ruli E, Ventura L. Robust inference for nonlinear regression models from the Tsallis score: application to Covid-19 contagion in Italy. Stat , 2020; 9:e309.
  * [14] Giummolé F, Mameli V, Ruli E, Ventura L. Objective Bayesian inference with proper scoring rules. Test, 2019; 28: 728–755.
  * [15] Good IJ. Rational decisions. J. Roy. Statist. Soc. B, 1952; 14: 107–114.
  * [16] Heritier S, Cantoni E, Copt S, Victoria-Feser MP ). Robust Methods in Biostatistics. Wiley; 2009.
  * [17] Heritier S, Ronchetti EM. Robust bounded-influence tests in general parametric models. J. Americ. Statist. Assoc., 1994; 89: 897–904.
  * [18] Huber PJ, Ronchetti EM. Robust Statistics. Wiley, New York; 2009.
  * [19] Hyvärinen A. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 2005; 6: 695–709.
  * [20] Hjort NL, Schweder T. Confidence distributions and related themes.J. Statist. Plan. Infer., 2018; 195: 1–13.
  * [21] Machete R. Contrasting probabilistic scoring rules. J. Statist. Plann. Inf., 2013; 143: 1781–1790.
  * [22] Mameli V, Musio M, Ventura L. Bootstrap adjustments of signed scoring rule root statistics. Comm. Statist. - Simul. Comput., 2018; 47: 4, 1204–1215.
  * [23] Mameli V, Ventura L. Higher-order asymptotics for scoring rules. J. Statist. Plann. Inf., 2015; 165: 13–26.
  * [24] Mandas A. Congiu MG, Abete C, Dessì S, Manconi PE, Musio M, Columbu S., Racugno. Cognitive decline and depressive symptoms in late-life are associated with statin use: evidence from a population-based study of Sardinian old people living in their own home. Neur. Res., 2014; 36: 3, 247–254.
  * [25] Nielsen F, Nock R. A closed-form expression for the Sharma-Mittal entropy of exponential families. J. Phys. A: Math. Theor., 2012; 45: 032003.
  * [26] Pak RJ. The minimum density power divergence estimation for the lognormal density. Comm. Stat. - Theory and Methods, 2014; 43: 4582–4588.
  * [27] Pauli F, Racugno W, Ventura L. Bayesian composite marginal likelihoods. Statistica Sinica, 2011; 21: 149–164.
  * [28] Ronchetti E, Ventura L. Between stability and higher-order asymptotics. Stat. and Comput., 2001; 11: 67–73.
  * [29] Ruli E, Ventura L. Can Bayesian, confidence distribution and frequentist inference agree?. Statistical Methods & Applications, 2020; DOI s10260-020-00520-y.
  * [30] Schweder T, Hjort NL . Confidence, Likelihood, Probability: Statistical Inference with Confidence Distributions. Cambridge University Press; 2016.
  * [31] Tsallis C. Possible generalization of Boltzmann-Gibbs statistics. J. Statist. Physics, 1988; 52: 479–487.
  * [32] Varin C, Reid N, Firth D. An overview of composite likelihood methods. Statist. Sinica, 2011; 21: 5–42.
  * [33] Ventura L, Racugno W. Recent advances on Bayesian inference for $P(X<Y)$, Bayesian Analysis, 2011; 6: 411–428.
  * [34] Xie M, Singh K. Confidence distribution, the frequentist distribution estimator of a parameter: a review. Int. Statist. Rev., 2013; 81: 3–39.


# Immutability and Design Patterns in Ruby

Seamus Brady<EMAIL_ADDRESS>
(28/03/2013)

###### Abstract

Functional Programming has seen a resurgence in interest in the last few years
and is often mentioned in opposition to Object-Orientated Programming.
However, Object-Orientated Programming can actually absorb some of the lessons
of Functional Programming to improve performance and reliability. This paper
examines three ways in which Immutability, a common Functional Programming
technique, can be used to develop immutable objects in Ruby.

##  1 Introduction

Functional Programming has left it’s mathematical, academic environment and
become much more visible in recent years. Languages such as Erlang, Haskell
and Clojure are being used in many corporate environments.[1]

Immutable data is a standard feature of Functional Programming [2].
Immutability, when applied to objects, creates objects that cannot be modified
after they have been created [6]. Immutable Objects have certain advantages
over standard Mutable Objects. This paper considers how a programmer using
Ruby, a popular Object-Orientated programming language, can create such
Immutable Objects.

Section 2 outlines the Functional Programming approach and how some of it’s
advantages may be appropriated by Object-Orientated programmers.

There are recommended strategies for creating immutable objects in strongly
typed languages such as Java [5]. In Section 3 three of these strategies have
been reworked into “Proto Design Patterns” to fit Ruby’s dynamically typed
nature.

###  1.1 Resources and Methods

This problem was approached by researching existing advice for creating
Immutable Objects in Java and redefining the approaches suggested to work with
Ruby.

##  2 Background

###  2.1 The Problem of Mutable Objects

In Object-Orientated programming, objects contain both program state (data)
and the means to operate on that state (methods). This simple encapsulation
allows an Object-Orientated language to model complex systems very
efficiently.

However, this approach has two difficulties:

> Mutable objects are not thread-safe by default and it’s easy for clients to
> change their state outside our control. [7]

In an era when more and more systems have multiple processors, thread safe
code is essential for certain areas of a system that demand concurrent
programming [8]. Also, as systems get ever more complex, subtle bugs caused by
mutable objects in complex data structures are very difficult to fix [4].

###  2.2 Monkey Patching in Ruby

Ruby is a general purpose Object-Orientated programming language [12]. As it
is a dynamically typed “scripting” language it is often used in preference to
Java. The challenge in designing Mutable Objects in Ruby is that, unlike Java,
classes can be changed (“monkey patched”) as any point in the program [13].

This means that Ruby is endlessly flexible, but monkey patching can create
very subtle and difficult to find bugs [14].

If an object is being used in part of the system that demands thread safety,
or where unwanted changes in state could cause problems, using an Immutable
Object could help avoid bugs.

###  2.3 Functional Programming Approach

Functional Programming is a distinct programming approach that treats a
computer program as a series of functions to be evaluated:

> Functional programming is so called because a program consists entirely of
> functions. The main program itself is written as a function which receives
> the program’s input as its argument and delivers the program’s output as its
> result.[1]

One of the main reasons that Functional Programming is recommended by
developers is that changes in the program’s state are avoided as part of the
Functional approach:

> …functional programs contain no side-effects at all. A function call can
> have no effect other than to compute its result. This eliminates a major
> source of bugs, and also makes the order of execution irrelevant \- since no
> side-effect can change the value of an expression, it can be evaluated at
> any time. [1]

To a commercial programmer who spends too much time debugging problems that
arise from inadvertent changes to program state, this seems like a dream come
true [4]. However, Functional Programming is often dismissed as an viable
alternative for two reasons (amongst others):

  * •

Functional Programming has a reputation for being difficult [3].

  * •

Most developers have spent time and money investing in alternative approaches
to programming and may be reluctant to change to something new.

###  2.4 Immutable Objects By Default

> Classes should be immutable unless there’s a very good reason to make them
> mutable….If a class cannot be made immutable, limit its mutability as much
> as possible. [9]

However, Object-Orientated programmers do not have to take up Functional
Programming completely to avoid the issues identified above [11]. The “no side
effects” advantage of Functional Programming can be achieved very simply in
Object-Orientated by creating all objects as Immutable by default. This
addresses the problems identified above, as well as providing other advantages
[10]:

  * •

Immutable Objects are thread safe by design.

  * •

They are simple to create and copy.

  * •

They do not need to be validated after creating as they do not change.

  * •

Objects can be made mutable only when needed.

##  3 Immutable Design Patterns

###  3.1 Proto-Patterns Versus Design Pattern

The three Immutable Design Patterns defined below are called “design patterns”
as part of this paper, but they may be more correctly called “proto-patterns”
rather than full Design Patterns.

The standard definition of a Design Pattern, according to the Portland Pattern
Repository [15], is that they occur in at “least three systems”. The patterns
described here may be better described as “proto-patterns”, defined as
“something documented in a pattern form, but lacks supporting known uses”
[16]. I leave it to the reader to decide whether they have seen these patterns
in use. The term Immutable Design Patterns is used with this caveat.

###  3.2 Immutable Object Pattern

####  3.2.1 Description

This design pattern is the simplest of the three as it just describes the
creation of an Immutable Object. Listing 1 illustrates this pattern.

The class ImmutablePerson has the following properties:

  * •

The constructor takes all the necessary values to construct the object
completely.

  * •

There are no mutator methods on the class.

  * •

All values returned by accessor method return a copy of the instance variable
value, rather than the original value. This is to avoid changes to attribute
values as per Listing 2.

  * •

The constructor calls the Ruby freeze method which is somewhat similar to the
final class attribute in Java. Calling freeze on a Ruby object makes it
Immutable. This means that after the constructor has completed processing the
class cannot be changed.

Listing 1: Immutable Object Example 1

⬇

class ImmutablePerson

attr_reader :name

def initialize name, address

@name = name

@address = address

self.freeze

end

def address

return @address.clone

end

end

class Address

attr_accessor :line1, :line2, :line3

def initialize line1, line2, line3

@line1 = line1

@line2 = line2

@line3 = line3

end

end

####  3.2.2 Advantages and Disadvantages

  * •

This Design Pattern creates a completely Immutable Object in a simple way.

  * •

However, in a situation where you are working with existing Mutable Objects
this Design Pattern will be of no benefit. The next two patterns will be of
more use here.

Listing 2: Immutable Object Example 2

⬇

address = Address.new("Foo Street", "Bar Town", "Test County")

person = ImmutablePerson.new("Foo McBar", address)

puts person.address.line1 # "Foo Street"

person.address.line1 = "Bar Street"

puts person.address.line1 # still "Foo Street"

###  3.3 Immutable Subclass Pattern

####  3.3.1 Description

Listing 3 illustrates this pattern. This Design Pattern illustrates a
situation where a Mutable class has a Mutable and an Immutable subclass.

  * •

The class AbstractPerson is an abstract superclass.

  * •

AbstractPerson has two subclasses - MutablePerson, a standard Mutable Object.
ImmutablePerson is a Immutable Object.

  * •

ImmutablePerson calls the constructor on the superclass and then freezes the
class.

  * •

Rather than passing a copy of the attribute (as in Listing 1) to avoid changes
to attribute values, this example freezes the instance variable object. This
achieves the same outcome as illustrated in Listing 4. This may be more useful
in some situations where an instance variable value is difficult to clone
(memory resources etc.). This will also throw an error when changes to
instance variables are attempted. This error may be useful in debugging
unwanted changes to state.

Listing 3: Immutable Subclass Pattern Example 1

⬇

class AbstractPerson

attr_accessor :name, :address

def initialize name, address

@name = name

@address = address

end

end

class Address

attr_accessor :line1, :line2, :line3

def initialize line1, line2, line3

@line1 = line1

@line2 = line2

@line3 = line3

end

end

class MutablePerson < AbstractPerson

# additional behaviour here...

end

class ImmutablePerson < AbstractPerson

def initialize name, address

super name, address

self.freeze

@address.freeze # freeze inst vars

end

def address

return @address

end

end

####  3.3.2 Advantages and Disadvantages

  * •

This Design Pattern allows a system to differentiate between Mutable and
Immutable objects using inheritance. For instance this may be useful where
part of the object hierarchy are needed only for reporting and can be read
only.

  * •

In some situations inheritance may not be the best way to differentiate
between Mutable and Immutable objects. In these situations the next Design
Pattern would be more beneficial.

Listing 4: Immutable Subclass Pattern Example 2

⬇

address = Address.new("Foo Street", "Bar Town", "Test County")

person1 = MutablePerson.new("Foo McBar", address)

person1.address.line1 = "Bar Street"

puts person1.address.line1 # "Bar Street"

person2 = ImmutablePerson.new("Bar McFoo", address)

person2.address.line1 = "Bar Street" # error: can’t modify frozen

###  3.4 Immutable Adapter Pattern

####  3.4.1 Description

Listing 5 illustrates this pattern. This Design Pattern illustrates a
situation where a Mutable class is wrapped in an Immutable class.

  * •

The class Person is a standard Mutable class.

  * •

ImmutablePerson contains an instance of Person in an instance variable.

  * •

The constructor of ImmutablePerson constructs a new Person class and then
freezes the objects and it’s instance variable values.

  * •

The accessors of ImmutablePerson query the wrapped Person class. As per
Listing 6, this example freezes the instance variable objects. This will also
throw an error when changes to instance variables are attempted.

Listing 5: Immutable Adapter Pattern Example 1

⬇

class Person

attr_accessor :name, :address

def initialize name, address

@name = name

@address = address

end

end

class Address

attr_accessor :line1, :line2, :line3

def initialize line1, line2, line3

@line1 = line1

@line2 = line2

@line3 = line3

end

end

class ImmutablePerson

def initialize name, address

@mutable_person = Person.new name, address

@mutable_person.freeze

@mutable_person.address.freeze

end

def name

return @mutable_person.name

end

def address

return @mutable_person.address

end

end

####  3.4.2 Advantages and Disadvantages

  * •

This Design Pattern allows any Mutable Object to be wrapped and passed around
as an Immutable Object.

  * •

This Design Pattern, in keeping with the Adapter Design Pattern, will allow
the Mutable Object to be extended if necessary.

  * •

This is the most useful of the three Immutable Design Pattern in existing
systems but does involve indirection. For simplicity of understanding the
other two Immutable Design Patterns should be considered first.

Listing 6: Immutable Adapter Pattern Example 2

⬇

address = Address.new("Foo Street", "Bar Town", "Test County")

person1 = Person.new("Foo McBar", address)

person1.address.line1 = "Bar Street"

puts person1.address.line1 # "Bar Street"

person2 = ImmutablePerson.new("Bar McFoo", address)

person2.address.line1 = "Bar Street" # error: can’t modify frozen

##  4 Related Work

There has been no formal definition of Immutable Design Patterns in Ruby as
far as the author is aware but there is plenty information on Immutability
across the web:

  * •

The Portland Pattern Repository [17] contains much discussion on Mutable and
Immutable Objects across many programming languages.

  * •

Hamster is a set of “Efficient, Immutable, Thread-Safe Collection classes for
Ruby”. All Hamster collections are Immutable making them thread safe [18].

  * •

The Functional Ruby blog offers a wider perspective on Functional Programming
in Ruby [19].

##  5 Summary

In this paper, the three Immutable Design Patterns in Table 1 were discussed.
Implementing these Design Patterns will allow a programmer to take advantage
of the “no side-effects” properties of Functional Programming within an
Object-Orientated programming language such as Ruby.

Name | Description
---|---
Immutable Object Pattern | Object itself is immutable
Immutable Subclass Pattern | Immutable subclass of mutable parent class
Immutable Adapter Pattern | Immutable wrapper around mutable class
Table 1: Summary of Ruby Immutable Design Patterns

## References

  * [1] http://cufp.galois.com/2009/abstracts.html#ChristopherPiroEugeneLetuchy
  * [2] Why Functional Programming Matters, http://www.cse.chalmers.se/~rjmh/Papers/whyfp.pdf
  * [3] Functional Programming Is Hard, That’s Why It’s Good, http://dave.fayr.am/posts/2011-08-19-lets-go-shopping.html
  * [4] Consequences when using Mutable Fields in hashCode(), http://blog.mgm-tp.com/2012/03/hashset-java-puzzler/
  * [5] Use Inheritance or Delegation to Define Immutable Classes, http://www.informit.com/articles/article.aspx?p=20522
  * [6] Immutable Classes, http://user.cs.tu-berlin.de/~stephan/OOPL/Sather-Chap4-12.html#MARKER-9-135
  * [7] Imperative, Mutable Code, https://www.inkling.com/read/functional-programming-java-dean-wampler-1st/chapter-5/imperative-mutable-code
  * [8] Thread-safe design, http://www.javaworld.com/javaworld/javaqa/1999-04/01-threadsafe.html
  * [9] Bloch, Joshua. (2008). Effective Java (2nd Edition) . ISBN-10: 0321356683.
  * [10] Immutable objects, http://www.javapractices.com/topic/TopicAction.do?Id=29
  * [11] Functional programming in object oriented languages, http://www.harukizaemon.com/blog/2010/03/01/functional-programming-in-object-oriented-languages/
  * [12] http://www.ruby-lang.org/en/
  * [13] http://stackoverflow.com/questions/394144/what-does-monkey-patching-exactly-mean-in-ruby
  * [14] Monkeypatching is destroying Ruby, http://devblog.avdi.org/2008/02/23/why-monkeypatching-is-destroying-ruby/
  * [15] http://c2.com/cgi/wiki?RuleOfThree
  * [16] http://c2.com/cgi/wiki?ProtoPattern
  * [17] http://c2.com/cgi/wiki?WelcomeVisitors
  * [18] https://github.com/harukizaemon/hamster#readme
  * [19] http://functionalruby.com/