# On approximate robust confidence distributions

Bortolato, E., Ventura, L.
Department of Statistical Sciences
University of Padova
Via C. Battisti 241, 35121 Padova, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A confidence distribution is a complete tool for making frequentist inference
for a parameter of interest $\psi$ based on an assumed parametric model.
Indeed, it allows to reach point estimates, to assess their precision, to set
up tests along with measures of evidence for statements of the type
”$\psi>\psi_{0}$” or ”$\psi_{1}\leq\psi\leq\psi_{2}$”, to derive confidence
intervals, comparing the parameter of interest with other parameters from
other studies, etc.

The aim of this contribution is to discuss robust confidence distributions
derived from unbiased $M-$estimating functions, which provide robust inference
for $\psi$ when the assumed distribution is just an approximate parametric
model or in the presence of deviant values in the observed data. Paralleling
likelihood-based results and extending results available for robust scoring
rules, we first illustrate how robust confidence distributions can be derived
from the asymptotic theory of robust pivotal quantities. Then, we discuss the
derivation of robust confidence distributions via simulation methods. An
application and a simulation study are illustrated in the context of non-
inferiority testing, in which null hypotheses of the form
$H_{0}:\psi\leq\psi_{0}$ are of interest.

###### keywords:

Bayesian simulation , Confidence density , Discrepancies , $M$-estimators ,
Pivotal quantity , Robustness , Non-inferiority testing.

††journal: EcoSta journal

##  1 Introduction

In recent years there has been considerable interest in frequentist inference
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
distribution analysis is much more informative than providing a $(1-\alpha)$%
interval or a $p$-value for an associated hypothesis test.

Let the scalar parameter of interest be $\psi$. With inference on $\psi$ we
shall understand statements of the type ”$\psi>\psi_{0}$” or
”$\psi_{1}\leq\psi\leq\psi_{2}$”, etc., where $\psi_{0}$, $\psi_{1}$ and
$\psi_{2}$ are given fixed values. To each statement, a CD allows us to
associate how much confidence the data have in the statement. The plot in Fig.
1 gives an illustration on making inference using a confidence density for a
scalar parameter of interest $\psi$: point estimators, $(1-\alpha)$%
confidence interval, one-sided $p$-value and measure of evidence for
”$\psi_{1}\leq\psi\leq\psi_{2}$”.

Figure 1: Illustration of making inference on the scalar parameter of interest
$\psi$ using a confidence density: point estimators (mode, median, mean),
$(1-\alpha)\%$ quantile-type confidence intervals, one-sided $p$-value and
measure of evidence for $\psi_{1}<\psi<\psi_{2}$.

The standard theory for CDs evolves around the use of likelihood methods for a
scalar parameter of interest $\psi$ of a parametric model. 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. However, it is well-known that likelihood-based methods are not
robust when the assumed distribution is just an approximate parametric model
or in the presence of deviant values in the observed data. In this case, 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 (see,
e.g., Hampel et al., 1986, Huber and Ronchetti, 2009, Heritier et al., 2009,
and Farcomeni and Ventura, 2012, and references therein).

The aim of this paper is to discuss the derivation of robust confidence
distributions, together with their associated densities. Paralleling
likelihood-based results, we first illustrate that asymptotic robust CDs can
be obtained by using pivotal quantities derived from unbiased $M$-estimating
functions, extending the theory of robust scoring rules discussed in Hjort and
Schweder (2018) and Ruli et al. (2022) to the more general setting of
$M$-estimating functions. Secondly, we explore two simulation-based approaches
to derive robust CDs. The first one relies on a frequentist reinterpretation
of Approximate Bayesian Computation (ABC) techniques (see, e.g., Rubio and
Johansen, 2013, and Bee et al., 2017, Ruli et al., 2020, Thornton et al.,
2022). The second approach leverages a Montecarlo rejection algorithm for
obtaining a one-sided $p$-value function. An application and a simulation
study are illustrated in the context of non-inferiority testing (see, e.g.,
Rothman et al., 2012), in which null hypotheses of the form
$H_{0}:\psi\leq\psi_{0}$ are of interest in order to establish if a new
product is not unacceptably worse than a product already in use. Finally, in
conclusion, we mention the possibility to resort to a non-parametric
derivation of CDs based on integral probability semimetrics (Muller, 1997) or
pseudo-metrics (Huber and Ronchetti, 2009, Chap. 2). This approach makes use
of discrepancies among cumulatives density functions as nonparametric pivots
for performing inference on $\psi$, and an example is illustrated in a
misspecification context, as in Legramanti et al. (2022).

The paper is organized as follows. Section 2 reviews some background on CDs.
Section 3 discusses the derivation of first-order robust CDs from unbiased
$M$-estimating functions and of simulation-based CDs. An application and a
simulation study are presented in Section 4 in the context of non-inferiority
testing. Finally, Section 5 mentions the non-parametric derivation of CDs
based on integral probability semimetrics, and Section 6 discusses some
concluding remarks.

##  2 Background on confidence distributions

###  2.1 Approximate likelihood-based confidence distributions

Consider a sample $y=(y_{1},\ldots,y_{n})$ of size $n$ from a random variable
$Y$ with assumed parametric model $f(y;\theta)$, indexed by a $d-$dimensional
parameter $\theta$. Let $\theta=(\psi,\lambda)$, where $\psi$ is a scalar
parameter of primary interest and $\lambda$ represents the remaining $(d-1)$
nuisance parameters.

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

$\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
$\psi^{*}$, the zero-confidence level estimator of $\psi$. In regular cases,
$cc(\psi)$ is decreasing to the left of $\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\psi^{*}\\\
2C(\psi,y)-1,&\text{if}\,\,\psi\geq\psi^{*}.\end{array}\right.$

With $C(\psi)$ a CD, $[C^{-1}(0.05),C^{-1}(0.95)]$ becomes an equi-tailed 90%
confidence interval. Also, solving $cc(\psi)=0.90$ yields two cut-off points
for $\psi$, precisely those of a 90% confidence interval.

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 pivot
(Barndorff-Nielsen and Cox, 1994). Thus $Q(x)=P_{\theta}(q(\psi;Y)\leq x)$
does not depend on $\psi$, which implies that

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

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

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

If the natural pivot is decreasing in $\psi$, then $C(\psi)=1-Q(q(\psi;y))$.

In the likelihood framework, there are well-working large-sample
approximations for the behaviour of pivotal quantities and these lead to
constructions of CDs. For instance, 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}}},$ (3)

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

A pivotal quantity that typically works better than (3) 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))$ (4)

is a first-order asymptotic CD, which 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 (see, among others, Severini, 2000, Reid, 2003, Brazzale et al.,
2007, and references therein) is discussed, for instance, in Schweder and
Hjort (2016, Chap. 7); see also Ruli and Ventura (2021). One key formula is
the modified profile likelihood root

$\displaystyle
r_{p}^{*}(\psi)=r_{p}(\psi)+\frac{1}{r_{p}(\psi)}\log\frac{q_{p}(\psi)}{r_{p}(\psi)},$
(5)

which has a third-order standard normal null distribution. In (5), the
quantity $q_{p}(\psi)$ is a suitably defined correction term (see, e.g.,
Severini, 2000, Chapter 9). In practice, $r^{*}_{p}(\psi)$ is a higher-order
pivotal quantity obtained as a refinement of the likelihood root
$r_{p}(\psi)$, which allow us to obtain an asymptotically third-order accurate
CD, i.e. with error of order $O(n^{-3/2})$.

###  2.2 Simulation-based confidence distributions

In the framework of CDs obtained from pivotal quantities, when the
distribution of the pivot (2) is not known, approximate confidence
distributions can be obtained by bootstrapping (Efron, 1979). Indeed, this
method provides an estimate of the sampling distribution of a statistic, and
this empirical sampling distribution can be turned into an approximate
confidence distribution in several ways (see, e.g., Schweder and Hjort, 2016,
Chap. 7).

Bootstrap procedures can be distinguished into two major categories:
parametric and non parametric. With parametric bootstrap, the seeked
distribution is taken to be that of the statistics given pseudo-data,
generated by a consistent estimate of the parametric model. With the second,
instead, the pseudo-data generating scheme is a multinomial resampling of data
points. The error in bootstrap confidence intervals is in most of the cases of
order $O_{p}(n^{-1})$, becoming of order $O_{p}(n^{-3/2})$ when correction
procedures can be implemented, and for $t-$bootstrap (see DiCiccio and Romano,
1988, DiCiccio and Efron, 1996, Schweder and Hjort, 2016, Chap. 7).

For instance, consider a monotone transformation $h(\psi)$ of a scalar
parameter of interest, and let $q(\psi,y)=(h(\psi)-h(\hat{\psi}))/\hat{\tau}$
an approximate studentized-pivot, where $\hat{\psi}$ is the MLE of $\psi$ and
$\hat{\tau}$ is a suitable estimate of the pivot standard deviation. Let
$Q(\cdot)$ be the distribution function of $q(\psi,y)$. Then, a confidence
distribution for the parameter of interest is
$C(h(\psi))=Q\left((h(\psi)-h(\hat{\psi}))/\hat{\tau}\right)$, with
appropriate confidence density $cd(h(\psi))=\partial C(h(\psi))/\partial\psi$.
When $Q(\cdot)$ is unknown, it can be estimated via bootstrapping. Let
$h(\psi^{*})$ and $\hat{\tau}^{*}$ be the result of bootstrapping, then the
$Q(\cdot)$ distribution can be estimated as $\hat{Q}$, via bootstrapped values
of $q^{*}=q(\psi^{*},y^{*})=(h(\psi)-h(\hat{\psi}^{*}))/\hat{\tau}^{*}$. The
approximate CD is then

$C_{t-boot}(\psi)=\hat{Q}\left(\frac{h(\psi)-h(\hat{\psi})}{\hat{\tau}}\right).$

This bootstrap method applies even when $q(\psi,y)$ is not a perfect pivot,
but is especially successful when it is, because $q^{*}$ then has exactly the
same distribution $Q(\cdot)$ as $q(\psi,y)$. Note that the method
automatically takes care of bias and asymmetry in $Q(\cdot)$.

##  3 Robust confidence distributions

###  3.1 Asymptotic derivation from $M$-estimating functions

The class of $M$-estimators is broad and includes a variety of well-known
estimators. For example it includes the MLE, the maximum composite likelihood
estimator (see e.g. Varin et al., 2011), estimators based on proper scoring
rules (see, e.g., Dawid et al., 2016, and references therein), and classical
robust estimators (see e.g. Huber and Ronchetti, 2009, and references
therein).

Under broad regularity conditions, an $M$-estimator $\tilde{\theta}$ is the
solution of the unbiased estimating equation

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

and it is asymptotically normal, with mean $\theta$ and covariance matrix
$V(\theta)=K(\theta)^{-1}J(\theta)(K(\theta)^{-1})^{{\footnotesize{{}^{{}_{\sf
T}}}}}$, where $K(\theta)=E_{\theta}(\partial
g(\theta)/\partial\theta^{{\footnotesize{{}^{{}_{\sf T}}}}})$ and
$J(\theta)=E_{\theta}(g(\theta)g(\theta)^{{\footnotesize{{}^{{}_{\sf T}}}}})$
are the sensitivity and the variability matrices, respectively. The matrix
$V_{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)$. Let us denote with
$G(\theta)=\sum_{i=1}^{n}G(y_{i};\theta)$ the function such that $g(\theta)$
is the gradient vector, i.e. $g(y;\theta)=\partial
G(y;\theta)/\partial\theta$.

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}g(y;\theta),$ (6)

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
$g(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.

Paralleling likelihood-based results, asymptotic robust inference on the
scalar parameter of interest $\psi$ can be based on first-order pivots. With
the partition $\theta=(\psi,\lambda)$, the $M$-estimating function is
similarly partitioned as
$g(y;\theta)=(g_{\psi}(y;\theta),g_{\lambda}(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 $V_{g}$ and $V_{g}^{-1}$. Finally, let
$\tilde{\lambda}_{\psi}$ be the constrained $M$-estimate of $\lambda$, let
$\tilde{\theta}_{\psi}=(\psi,\tilde{\lambda}_{\psi})$, and let $\tilde{\psi}$
be the $\psi$ component of $\tilde{\theta}$. Then, a profile Wald-type
statistic for the $\psi$ may be defined as

$w_{R}(\psi)=(\tilde{\psi}-\psi)(\tilde{V}_{g}^{\psi\psi})^{-1/2},$

and it has an asymptotic $N(0,1)$ null distribution. Similarly, the profile
score-type statistic

$w_{sR}(\psi)=g_{\psi}(\tilde{\theta}_{\psi})^{T}K^{\psi\psi}(V_{g}^{\psi\psi})^{-1}K^{\psi\psi}g_{\psi}(\tilde{\theta}_{\psi})$

has an asymptotic $\chi^{2}_{1}$ null distribution, while the asymptotic
distribution of the profile ratio-type statistic for $\psi$, given by
$W_{R}(\psi)=2\left(G(\tilde{\theta}_{\psi})-G(\tilde{\theta})\right)$, is
$\nu\chi^{2}_{1}$, where
$\nu=(\tilde{K}^{\psi\psi})^{-1}\tilde{V}_{g}^{\psi\psi}$. In view of this,
for the adjusted profile ratio-type statistic to first-order it holds

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

Finally, the adjusted profile root, analogous to (2.1), can be defined as

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

which has an asymptotic standard normal distribution. For the general theory
of robust tests see Heritier and Ronchetti (1994).

Paralleling results in Section 2.1 for likelihood based CDs, a recipe to
derive an asymptotic CD from robust $M$-estimating functions is based on
pivotal quantites, extending the theory illustrated for robust scoring rules
in Hjort and Schweder (2018) and Ruli et al. (2022). To this end, let us
denote with $q_{R}(\psi;y)$ a robust pivotal quantity, such as the profile
Wald-type statistic $w_{R}(\psi)$ or the adjusted profile scoring rule root
$r_{R}(\psi)$. Then,

$\displaystyle
C_{R}^{w}(\psi)\,\dot{=}\,\Phi\left((\psi-\tilde{\psi})(\tilde{V}_{g}^{\psi\psi})^{-1/2}\right)$
(7)

and

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

are first-order asymptotic CDs, and the corresponding confidence densities
are, respectively,

$cd_{R}^{w}(\psi)\,\dot{=}\,\frac{\phi\left((\psi-\tilde{\psi})(\tilde{V}_{g}^{\psi\psi})^{-1/2}\right)}{\sqrt{\tilde{V}_{g}^{\psi\psi}}}$

and

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

where $\phi(\cdot)$ is the density function of the standard normal
distribution. Note that the Wald-type based confidence density
$cd_{R}^{w}(\psi)$ coincides with the asymptotic first-order robust Bayesian
posterior distribution for $\psi$ (see, e.g., Greco et al., 2008, and Ventura
and Racugno, 2016).

In practice, using for instance (8), the confidence median is $\tilde{\psi}$
and an $(1-\alpha)$ equi-tailed confidence interval can be obtained as
$\\{\psi:|r_{R}(\psi)|\leq z_{1-\alpha/2}\\}$, where $z_{1-\alpha/2}\\}$ is
the $(1-\alpha/2)$-quantile of the standard normal density. When testing, for
instance, $H_{0}:\psi=\psi_{0}$ against $H_{1}:\psi<\psi_{0}$, the $p$-value
is $p=C_{R}^{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_{R}(\psi_{0})|))$.
Furthermore, a measure of evidence for a statement of the form
”$\psi_{1}<\psi<\psi_{2}$” can be computed as
$C_{R}^{r}(\psi_{2})-C_{R}^{r}(\psi_{1})$.

To study the stability of robust CDs, let us write the robust pivotal quantity
more generally as $q_{R}(\psi;T(\hat{F}_{n}))$, where $\hat{F}_{n}$ is the
empirical distribution function and $T(F)$ is the functional defined by the
unbiased $M$-estimating equation $\int g(y;T(F))\,dF(y)=0$, where
$F=F(y;\theta)$ is the assumed parametric model. In CD inference the tail
area, given by $C_{R}(\psi)=\Phi(q_{R}(\psi;T(\hat{F}_{n}))$, plays a central
role and thus we can consider the tail area influence function (see, e.g.,
Field and Ronchetti, 1990, and Ronchetti and Ventura, 2001), given by

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

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$. It can be shown that

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

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

The application of (7) and (8) in the particular context of a robust scoring
rule has been discussed in Ruli et al. (2022). In particular, the Tsallis
score (Tsallis, 1988) is considered, which is given by

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

with corresponding unbiased $M$-estimating function $g(\theta)=\partial
G(y;\theta)/\partial\theta$ (Ghosh and Basu, 2013, Dawid et al., 2016), and
with the parameter $\gamma$ which gives a trade-off between efficiency and
robustness.

###  3.2 Derivation of robust confidence distributions via simulation

In the context of robust procedures, suitable modifications of the bootstrap
have been explored for improving the stability of inference. In fact, both
families of parametric and non parametric bootstrap presents some drawbacks.
For instance, with the non parametric bootstrap, despite the direct
specification of the data generating process is avoided and this robustifies
the analysis, the distribution of estimators might be unstable since the
amount of outliers of pseudo-samples can be higher than that in the original
dataset. For remedying this drawback, some authors proposed the use of weigths
associated to observations before performing repeatedly likelihood
maximization under the assumed model, as in weighted Likelihood Bootstrap
(Newton and Raftery, 1994). Lyddon et al. (2019) showed that the corresponding
estimator is asymptotically normal and with covariance matrix with the
classical structure of sandwich form. Also, Chen and Zhou (2020) introduced a
similar approach based on estimating equations and provided non-asymptotic
guarantees for the resulting errors. Moreover, a further problem when using a
bootstrap approach together with robust procedures as $M$-estimating equations
is the need of repeatedly solving numerically non-convex or complex
optimization problems, which may be computationally expensive.

In this section we inspect two alternative methods for computing CDs based on
robust $M$-estimating functions, that go beyond some limitations of the
bootstrap. Broadly speaking, the first method is based on a frequentist
reinterpretation of the ABC machinery (see, e.g., Bee et al., 2017, Ruli et
al., 2020, Thornton et al., 2022), whose properties have been derived by Rubio
and Johansen (2013) in a general setup. The idea consists in generating
candidate parameter values from an uniform distribution, computing a robust
suitable summary statistic using the simulated data and then accepting only
the parameter values such that the corresponding summary statistic is ”close”
to its observed counterpart (see Algorithm 1).

Algorithm 1 Accept-reject robust ABC

Input: proposal $p(\psi,\lambda)$, number of iterations $R$, robust summary
statistic $t(\cdot)$, $t^{obs}=t(y^{obs})$, where $y^{obs}$ is the observed
sample, tolerance $\varepsilon$, distance $\rho(\cdot;\cdot)$

for $j\in 1,\ldots,R$ do

Sample $(\psi_{j}^{*},\lambda_{j}^{*})\sim p(\psi,\lambda)$ and $y_{j}^{*}\sim
f(y;\psi_{j}^{*},\lambda_{j}^{*})$

Compute $t_{j}^{*}=t(y_{j}^{*})$

Accept $\psi_{j}^{*}$ if $\rho(t_{j}^{*};t^{obs})\leq\varepsilon$ else reject

end for

resample the accepted $(\psi^{*},\lambda^{*})$ with probability $\propto
1/p(\psi^{*},\lambda^{*})$

return robust approximate normalized pseudo-likelihood $\propto$ confidence
density $\hat{cd}_{R}^{abc}(\psi)$

In Algorithm 1, the summary statistics of Soubeyrand and Haon-Lasportes (2015)
or of Ruli et al. (2016, 2020) can be used. In particular, the first one is
based directly on the $M$-estimator $\tilde{\psi}$ as the summary statistic
$t(y)$ and a, possibly rescaled, distance among the observed and the simulated
value of the statistic. In the second one, a rescaled version of the
$M$-estimating function $g(\theta)$, evaluated at a fixed value of the
parameter, is used as a summary statistic $t(y)$; this avoids repeated
evaluations of the consistency correction involved in the $M$-estimating
function. Note also that when using the $M$-estimator as a summary statistic,
the algorithm for solving the estimating equation might not converge after a
prefixed number of iterations, thus causing additional noise in the results.
For a single parameter of interest, with the partition
$\theta=(\psi,\lambda)$, we propose to modify the algorithm of Ruli et al.
(2020) by using the profile estimating equation $g_{\psi}(y;\psi,\lambda)$ and
plugging in the value of proposals $\lambda^{*}$ for nuisance parameters used
to generate pseudo-data. The treatment of the nuisance parameters resembles as
a generalized profile likelihood computation.

Note that, assuming the regularity assumptions of Soubeyrand and Haon-
Lasportes (2015) and the usual regularity conditions on $M$-estimators (Huber
and Ronchetti, 2009, Chap. 4), then for $n\to\infty$ the robust confidence
densities derived via simulation are asymptotically equivalent to the Wald-
type confidence density $cd_{R}^{w}(\psi)$. Moreover, following Ruli et al.
(2020), if $g(y;\theta)$ is bounded in $y$, i.e. if the $M$-estimator is
B-robust, then asymptotically the posterior mode, as well as other posterior
summaries of the robust confidence density $\hat{cd}_{R}^{abc}(\psi)$ have
bounded $IF$.

The second method for computing CDs based on robust $M$-estimating functions
is similar but aims at computing directly a Montecarlo $p$-value or a
significance function using a different acceptance rule (see Bortolato and
Ventura, 2022): again the parameter values are sampled from a uniform
distribution, a robust summary statistic $t(y^{*})$ using the simulated data
is obtained, and then the proposed parameter is accepted if the robust summary
statistic $t^{*}$ is greater of its observed counterpart $t^{obs}$ (see
Algorithm 2). For obtaining the confidence density, if the obtained CD is
indeed monotone increasing, hence it is a proper CD, Algorithm 3 can be used.

Algorithm 2 Accept-reject confidence curve and confidence distribution
computing

Input: proposal $p(\psi,\lambda)$, number of iterations $R$, robust summary
statistic $t(\cdot)$, $t^{obs}=t(y^{obs})$, where $y^{obs}$ is the observed
sample.

for $j\in 1,\ldots,R$ do

Sample $(\psi_{j}^{*},\lambda_{j}^{*})\sim p(\psi,\lambda)$ and $y_{j}^{*}\sim
f(y;\psi_{j}^{*},\lambda_{j}^{*})$

Compute $t_{j}^{*}=t(y_{j}^{*})$

Accept $\psi_{j}^{*}$ if $t_{j}^{*}\geq t^{obs}$ else reject

end forreturn $\psi^{*}\text{ with density }\propto cc_{R}(\psi)$ a confidence
curve/distribution

Algorithm 3 Confidence density

Input: robust confidence distribution $C_{R}(\psi),$ desired size $R$, grid of
values $G=\\{\frac{1}{R},\frac{2}{R},\ldots,\frac{R-1}{R},1\\}$

for $j\in 1,\ldots,R$ do

Compute the empirical quantiles of the CD: $\psi_{j}^{*}=C_{R}^{-1}(G_{j})$

end for

return $\psi^{*}\sim\hat{cd}_{R}(\psi)$.

As a final remark, we observe that, for obtaining stable results with the
accept-reject schemes, it is suggestable to increase the number of proposals
as the dimesion of the parameter space increases, in order to have resonable
acceptance rate with ABC-type algorithm, and in general for obtaining more
precise estimation of the confidence distributions.

##  4 Applications to non-inferiority tests

The aim of this section is to introduce and apply CDs inference in the context
of non-inferiority testing, in which interest is in establishing if a new
product is not unacceptably worse than a product already in use. Applications
of non-inferiority testing has revealed an attractive problem in medical
statistics, biostatistics, statistical quality control and engineering
statistics, among others. Here we focus in non-inferiority clinical trials
where the aim is to show that an experimental treatment is not (much) worse
than a standard treatment. Clinical practice, however, is not the only field
of application of these tests: in comparing the performance of sensors in
industrial environment, for instance, the margin may be linked to some
difference in costs due to sensor functioning. Other applications can be found
in machine learning literature, where instead the meaningful margin is related
to the accuracy or to the speed in classification tasks.

In the process of evaluating the efficacy of an experimental treatment, it is
common to develop studies in which the two arms are the new and the standard
therapy, respectively, rather than the new and the placebo. This is because it
is considered unethical to deprive patients from a therapy that has already
been proven to be beneficial. The underlying research hypothesis to be
verified is that new therapies have equivalent or non-inferior efficacies to
the ones currently in use. Both non-inferiority and superiority tests are
examples of directional (one-sided) tests (see, e.g., D’Agostino et al., 2003,
Rothmann et al., 2012, and references therein). In particular, the non-
inferiority test wants to test that the treatment mean $\mu_{N}$ is not worse
than the reference mean $\mu_{S}$ by more than a given equivalence margin
$\delta$. The actual direction of the hypothesis depends on the response
variable being studied. This question can be formulated into a test procedure
for which the null hypothesis is

$H_{0}:\mu_{S}-\mu_{N}\geq\delta,{}$

where $\delta>0$ is the equivalence margin, when higher values of the response
variable mean better results, versus

$H_{1}:\mu_{S}-\mu_{N}<\delta.$

The scalar parameter of interest in this context is thus
$\psi=\mu_{S}-\mu_{N}$, and non-inferiority is claimed when the null
hypothesis is rejected.

Figure 2: Testing procedure with traditional comparative studies and non-
inferiority studies using confidence densities: vertical lines represent the
equivalence/non-inferiority margin ($\delta=-5$).

The equivalence margin $\delta$ corresponds to the practical acceptable
difference and should be pre-specified before the data is recorded (see e.g.
Garret, 2003). An overly conservative margin might result in a high risk of
not being able to claim non-inferiority when it actually is non-inferior.
Conversely, overly liberal margins could result in a high risk of claiming
non-inferiority when it actually is not non-inferior. A reasonable margin
would be best derived from a combination of factors: the expected event rate,
the duration of follow-up, and the number and nature of the events. However,
arbitrary clinical judgment and the sponsor budget are of a great influence,
resulting in a somewhat subjective non-inferiority margin. It is not clear in
some situations how to perform the choice, and multiple thresholds could be
plausible; in this respect, CDs are particularly useful to perform sensitivity
analyses. Indeed, in this situation a confidence distribution on the
difference $\psi=\mu_{S}-\mu_{N}$ will simultaneously show the evidence of the
$p$-value against the null for a series of values $\delta$, and decide for a
reasonable $\delta$ with the nominal control of the rejection level and
possible alternatives.

Here we consider an example of trial where higher levels of the response
variable mean that the new treatment is effective. The aim is verifying that
the new treatment ($N$) is not unacceptably worse to the standard ($S$). Let
us assume that $n=80$ patients are randomized into two groups, and the model
for the data is assumed to be

$Y_{S}=\mu_{N}+\psi+u,\\\ \quad Y_{N}=\mu_{N}+u,\quad u\sim
N(0,\,\sigma^{2}).$ (11)

The normal distribution on the error term is often the basis of statistical
analyses in medicine, genetics and in related sciences. Under this assumption,
parametric inferential procedures based on the sample means, standard
deviations, two-samples $t$-test, and so on, are the most efficient. However,
it is well known that they are not robust when the normal distribution is just
an approximate parametric model or in the presence of deviant values in the
observed data (see, e.g., Farcomeni and Ventura, 2012). In the framework
described by (11), we inspect the effect of adding some contamination in the
data of the new treatment group. In particular, in the contaminated scenario,
$10\%$ of the error terms in the new treatment group are half-Cauchy
distributed (see Fig. 3).

Figure 3: Boxplots of the new treatment group and of the standard group under
the two scenarios. Red dots indicate group means, green dots group medians.

It is of interest to compare CDs inference for $\psi$ based on the following
approaches (abbreviations are also used in Fig. 4 and in the following) used
to derive confidence densities:

  1. 1.

exact classical Wald-type confidence density based on $w_{p}(\psi)$, which is
related to the classical two sample $t$-test (Wald/Mean)

  2. 2.

robust asymptotic Wald-type confidence density $cd_{R}^{w}(\psi)$ based on the
Huber’s estimator (Wald/M-test)

  3. 3.

approximate confidence density based on ABC (Algorithm 1) with robust Huber’s
estimator as summary statistics (ABC/M-est)

  4. 4.

approximate confidence density based on ABC (Algorithm 1) with the robust
Huber’s estimating equation as summary statistic (ABC/M-EE)

  5. 5.

simulated confidence density (Algorithm 2) based on the robust Huber’s
estimator (CDensity/M-est)

  6. 6.

simulated confidence density (Algorithm 2) based on the robust Huber’s
estimating equation (CDensity/M-EE)

  7. 7.

approximate confidence density based on ABC (Algorithm 1) with the difference
of medians as summary statistics (ABC/Median)

  8. 8.

simulated confidence density (Algorithm 2) based on the difference of medians
(CDensity/Median)

  9. 9.

parametric bootstrap confidence density (Boot/Basic)

  10. 10.

parametric bootstrap with normal intervals confidence density (Boot/Norm)

  11. 11.

parametric bootstrap with percentiles confidence density (Boot/Perc)

The nominal value of the mean difference between the treatment effects is
$\psi_{0}$ is fixed to 2.6, and for simulation-based confidence distribution
as well as for those obtained by the ABC-type algorithm we used $10^{5}$
proposals and a tolerance level of $0.1$. In the Huber’s estimator we fix the
tuning constant which controls the desired degree of robustness to 1.345,
which imply that the estimator is 5% less efficient than the corresponding MLE
under the assumed model.

From the resulting confidence densities illustrated in Fig. 4 we note that,
when the data come from the central model (left column) all the confidence
densities are in reasonable agreement, even if the confidence densities based
on the medians behave slightly worse, with a greater variability. When the
data are contaminated (right column), the non robust confidence density
Wald/Mean is less trustworthy as it drifts away from the true parameter value
(green dotted line). This is not the case however for the robust confidence
densities which remains centred around the true parameter value. We further
note that in the contaminated case, the robust confidence densities based on
the $M$-estimating equation (ABC/M-EE and CDensity/M-EE) are the ones with the
smallest variability. For all these confidence densities, Table 1 gives the
measures of evidence for the equivalence margin $\delta$ taken equal to 4
(black dotted line in Fig. 4), that is for the statement ”$\psi>\delta$”. As a
reference, consider the result derived by the exact $t-$distribution of the
exact classical Wald-type confidence density in the non contaminated case,
which is $0.08$. The results, without and with the contamination, confirm the
behaviour of the confidence densities in Fig. 4, in particular the non
robustness of the likelihood-based confidence density (Wald/Mean). The most
stable values under contamination seem to be those obtained with M-EE
approaches (0.09 with ABC/M-EE and 0.05 with CD/M-EE). The same analysis could
be done in principle, given the CD, for any margin $\delta$.

Figure 4: Confidence densities for $\psi$ based on $10^{5}$ proposals, without (left column) and with contamination (right column). Results according to a different types of statistic are reported for each row, inferential techniques are represented with different colors. Black vertical dotted line is $\delta$, green vertical dotted line is $\psi_{0}$. Method |  $0\%$ cont. |  $10\%$ cont.
---|---|---
Wald/Mean | 0.08 | 0.42
Wald/M-test | 0.07 | 0.03
ABC/Median | 0.24 | 0.20
ABC/M-EE | 0.11 | 0.09
ABC/M-est | 0.12 | 0.09
CDensity/Median | 0.20 | 0.21
CDensity/M-EE | 0.08 | 0.05
CDensity/M-est | 0.14 | 0.11
Table 1: Confidence measures of evidence for the null hypothesis
$H_{0}:\psi>\delta$, with $\delta=4$ associated to Fig. 4, without and with
contamination.

###  4.1 Simulation study

For investigating the behaviour of the several confidence densities, we
perform a simulation study under two sample sizes $n=40,80$ (20, 40 per group)
and for each of them we investigate two scenarios: one in which the
assumptions of the model in (11) are met by the true data generating model,
and the second one where $10\%$ of the error terms in the new treatment group
are half-Cauchy distributed (as in Fig. 3). The families of methods to derive
the confidence densities considered are the same as in the example above, and
confidence distributions construction are again based on exact and asymptotic
pivotal quantities or simulation-based. For the Rejection-ABC-type confidence
distributions the tolerance for the discrepancy was still set to $0.1$, the
true value of parameter of interest to $\psi_{0}=2.6$, and the Huber’s tuning
constant to 1.345. The proposals for $\psi$ were taken Uniform in $[-3,9]$,
for the parameter $\mu_{N}$ we sample from a Uniform $[110,130]$, while for
$\sigma$ we generated values from a Uniform $[1,8]$. For the simulations, 4000
values were generated from the proposals and a total of 2000 simulations were
performed.

We compare the empirical coverages of $90\%$ and $95\%$ equi-tailed confidence
intervals. Results are synthetized in Tables 2 and 3. We also report in Tables
4 and 5 the error associated to confidence median point estimators, in terms
of bias ($b=\sum_{r=1}^{R}\tilde{\theta}_{r}-\theta_{0}$), probability of
overestimation ($PO=\sum_{r=1}^{R}1_{\\{\tilde{\theta}_{r}>\theta_{0}\\}}$)
and I type error with $\alpha=0.05$. We note that, under the central model,
the Wald/Mean CD shows a good performance, as well as some robust CDs
(Wald/M-test, CDensity/M-EE and ABC/M-EE). With contaminated data, the
Wald/Mean CD tends to be affected by contamination, whereas the robust CDs
perform substantially better, with the CDs based on $M$-estimating equations
being preferred over those based on $M$-estimators. Asymptotic symmetric
confidence densities based on Wald-type robust CDs and ABC-type confidence
densities seem to be affected more by bias than the simulated CDs (see Tables
4 and 5). Note finally that ABC-type results, even if behaving well, depend on
a tolerance choice, hence the results may degradate when the latter is not
well calibrated.

As a final remark, note that an interesting aspect of this simulation study
was the difference among the approach of using robust $M$-estimating functions
instead of robust $M$-estimates, especially in the treatment of nuisance
parameters. For the M-est CD based on $\tilde{\psi}$, for each fixed values
$(\psi^{*},\lambda^{*})$, the corresponding acceptance probability is
$Pr_{\psi^{*},\lambda^{*}}\left(\tilde{\psi}^{*}>\tilde{\psi}\right)$. This
might be recognized to be similar to a bootstrap $p$-value, with the exception
of what is the model considered for the simulation. In contrast, for the CD
with the profile $M$-estimating function as summary statistic, the acceptance
probability for same values $(\psi^{*},\lambda^{*})$ is
$Pr_{\psi^{*},\lambda^{*}}\left(g_{\psi}(y;\tilde{\psi},\lambda^{*})\geq
0\right)$, as if $\lambda^{*}$ was the oracle estimate. Hence, the CD for each
$\psi$ is associated to an average $p$-value, i.e. $CD(\psi)=\int
Pr_{\psi,\lambda}(g_{\psi}(y;\tilde{\psi},\lambda)\geq 0)d\lambda.$ Note that
under the true model, with $(\psi_{0},\lambda_{0})$ true parameter value, the
corresponding $p$-value would be equal to the one without the nuisance
parameter $Pr_{\psi_{0},\lambda_{0}}(g_{\psi}(y;\tilde{\psi},\lambda_{0})\geq
0).$ Note that, instead, keeping fixed the nuisance parameters in the
simulations would correspond to consider them as known.

Contamination | 0% | 10%
---|---|---
$n=40$ | 95% CI | 90% CI | 95% CI | 90% CI
Wald/Mean | 93.9 | 89.1 | 97.1 | 94.0
Wald/M-test | 93.7 | 88.4 | 94.1 | 88.3
ABC/Median | 97.1 | 93.4 | 97.7 | 93.7
ABC/M-EE | 92.7 | 87.2 | 93.7 | 88.9
ABC/M-est | 97.0 | 93.1 | 97.6 | 93.9
CDensity/Median | 99.5 | 97.6 | 99.2 | 98.0
CDensity/M-EE | 95.8 | 90.5 | 96.7 | 92.1
CDensity/M-est | 99.4 | 97.3 | 99.2 | 98.0
Boot/basic | 93.4 | 88.0 | 92.3 | 86.1
Boot/Norm | 93.5 | 88.2 | 92.4 | 86.0
Boot/Perc | 93.4 | 87.9 | 92.3 | 86.1
Table 2: Empirical coverages in a simulation study without and with 10% of contamination and $n=40$.  Contamination | 0% | 10%
---|---|---
$n=80$ | 95% CI | 90% CI | 95% CI | 90% CI
Wald/Mean | 95.5 | 90.0 | 95.9 | 92.2
Wald/M-test | 95.1 | 89.6 | 93.9 | 87.9
ABC/Median | 97.2 | 93.5 | 97.3 | 93.5
ABC/M-EE | 93.3 | 86.9 | 92.7 | 87.3
ABC/M-est | 97.5 | 93.6 | 97.2 | 93.9
CDensity/Median | 99.1 | 97.6 | 99.2 | 97.5
CDensity/M-EE | 95.9 | 89.4 | 96.4 | 92.5
CDensity/M-est | 99.1 | 97.3 | 99.2 | 97.7
Boot/Basic | 94.1 | 89.5 | 92.3 | 87.5
Boot/Norm | 94.2 | 89.5 | 92.4 | 87.4
Boot/Perc | 94.3 | 89.6 | 92.3 | 87.5
Table 3: Empirical coverages in a simulation study without and with 10% of contamination and $n=80$.  Contamination | 0% | 10%
---|---|---
$n=40$ | $|b|$ | $PU$ | I type err. | $|b|$ | $PU$ | I type err.
Wald/Mean | 0.01 | 0.51 | 0.06 | 5.57 | 0.65 | 0.03
Wald/M-test | 0.00 | 0.51 | 0.06 | 0.23 | 0.42 | 0.08
ABC/Median | 0.03 | 0.52 | 0.01 | 0.09 | 0.46 | 0.01
ABC /M-EE | 0.00 | 0.51 | 0.03 | 0.23 | 0.42 | 0.03
ABC/M-est | 0.01 | 0.51 | 0.01 | 0.23 | 0.42 | 0.01
CDensity/Median | 0.15 | 0.55 | 0.01 | 0.03 | 0.51 | 0.01
CDensity/M-EE | 0.11 | 0.55 | 0.03 | 0.11 | 0.46 | 0.03
CDensity/M-est | 0.13 | 0.56 | 0.01 | 0.09 | 0.46 | 0.01
Boot/Basic | 0.84 | 0.75 | 0.06 | 1.07 | 0.79 | 0.10
Boot/Norm | 0.84 | 0.75 | 0.06 | 1.07 | 0.79 | 0.10
Boot/Perc | 0.84 | 0.75 | 0.06 | 1.07 | 0.79 | 0.09
Table 4: Measures of stability of CDs: absolute bias ($|b|$), probability of underestimation ($PU$) and I type error $(\alpha=0.05)$ of confidence estimators (medians) in the simulation study with $n=40$. Contamination | 0% | 10%
---|---|---
$n=80$ | $|b|$ | $PU$ | I type err. | $|b|$ | $PU$ | I type err.
Wald/Mean | 0.02 | 0.49 | 0.05 | 1.76 | 0.58 | 0.05
Wald/M-test | 0.01 | 0.49 | 0.05 | 0.19 | 0.42 | 0.08
ABC/Median | 0.02 | 0.50 | 0.02 | 0.05 | 0.49 | 0.02
ABC/M-EE | 0.02 | 0.48 | 0.03 | 0.19 | 0.42 | 0.03
ABC/M-est | 0.01 | 0.49 | 0.01 | 0.19 | 0.42 | 0.02
CDensity/Median | 0.07 | 0.53 | 0.01 | 0.02 | 0.51 | 0.02
CDensity/M-EE | 0.08 | 0.53 | 0.03 | 0.11 | 0.45 | 0.03
CDensity/M-est | 0.08 | 0.54 | 0.01 | 0.11 | 0.46 | 0.02
Boot/Basic | 0.03 | 51.1 | 0.05 | 0.39 | 0.33 | 0.09
Boot/Norm | 0.03 | 51.1 | 0.05 | 0.39 | 0.33 | 0.09
Boot/Perc | 0.03 | 51.1 | 0.05 | 0.39 | 0.33 | 0.09
Table 5: Measures of stability of CDs: absolute bias ($|b|$), probability of
underestimation ($PU$) and I type error $(\alpha=0.05)$ of confidence
estimators (medians) in the simulation study with $n=80$.

###  4.2 Real data application

A class of problems requiring similar considerations to those of non-
inferiority tests, i.e. sensitivity analysis with respect to the reference
margin $\delta$, is that of superiority studies (see Fig. 5).

Figure 5: Example of making inference with confidence densities in superiority
test, with margin $\delta=-3.5$.

Here we analyze the data collected in a randomized controlled trial (see
Carhart-Harris et al., 2021, Nayak et al., 2022) with the aim of assessing the
superiority of a new therapy with psilocybin (P) versus that with escitalopram
(E), in treating major depressive disorder. The dataset contains the scores
obtained by $n=57$ patients on a questionnaire, before and after a 6-week
period of therapy.

Figure 6: Real data example: boxplots of pre-post differences of scores in the
group using psilocybin (P) versus escitalopram (E) (left), confidence
densitites for the parameter $\beta_{2}$, measuring the difference in efficacy
of the two therapies (right).

The model considered for the scores at the time of follow-up (FU) is the
following

$y_{FU}=\beta_{0}+\beta_{1}y_{BL}+\beta_{2}\text{P}+u,\quad u\sim
N(0,\sigma^{2}),$

where $y_{BL}$ represents the value at the baseline and P is a dummy variable
that equals 1 if the subject belongs to the group treated with the new therapy
(psilocybin), and thus the coefficient relates to the additional change with
respect to the control group (escitalopram) after the therapy. A reduction of
the score indicates a clinical improvement; thus superiority is claimed if the
estimate of the coefficient $\beta_{2}$ is sufficiently lower than 0. In
particular, in order to conclude in favour of meaningful superiority, the
clinicians considered as reference a margin $\delta=-5.3$. It is of interest
to provide stable measures of evidence for the statement ”$\beta_{2}>\delta$”,
with $\delta=-5.3$ ($H_{0}$).

The MLE for the parameter $\beta_{2}$ and its standard error are,
respectively, $-5.32$ and $1.44$, while the robust counterparts are $-6.18$
and $1.33$. Note that after removing two outliers the MLE become -6.23, with
standard error 1.35. We resume the whole confidence densities based on Wald-
type methods together with simulated confidence densities based on Huber’s
estimators and Huber’s estimating equations in Fig. 6. As it can be noted the
classical confidence density (Wald/Mean) is shifted to the right, because of
the presence of outliers. Evidence measures for different margins are reported
in Table 6. Using the margin chosen by the clinicians (-5.3) there is no
evidence of superiority at level $\alpha=0.1$; however note that the measure
of evidence computed with the Wald-type confidence density (Wald/Mean) is the
double of the ones computed with the robust confidence densities. With a
margin of $\delta=-3.5$ all the robust procedure would agree in claiming
superiority with $\alpha=0.1$, while according to classical Wald-type
confidence density (Wald/Mean) there would not be enough evidence to conclude
superiority.

$-\delta$ | -3.5 | -4 | -4.5 | -5 | -5.3
---|---|---|---|---|---
Wald/Mean | 0.11 | 0.18 | 0.29 | 0.41 | 0.49
Wald/M-test | 0.02 | 0.05 | 0.10 | 0.19 | 0.26
ABC/M-est | 0.00 | 0.00 | 0.14 | 0.14 | 0.29
ABC/M-EE | 0.08 | 0.12 | 0.18 | 0.24 | 0.30
CDensity/M-est | 0.03 | 0.09 | 0.14 | 0.21 | 0.28
CDensity/M-EE | 0.07 | 0.14 | 0.17 | 0.22 | 0.28
Table 6: Measures of evidence for the hypothesis ”$\beta_{2}>\delta$” for
several margins.

##  5 Confidence distributions based on integral probability semimetrics

For mitigating the effects of small departures from model assumptions, and
possible dramatic changes of inferential conclusions due to unconvenient
choices of pivotal quantities and summaries, the use of non parametric
procedures may be an alternative.

Here, we mention the construction of CDs based on discrepancy measures defined
on the space of distributions, belonging to the class of integral probability
semimetrics (Muller, 1997) or pseudo-metrics (Huber and Ronchetti, 2009,
Chapt. 2). These divergences are classically associated to the concept of
stability, used as global tests and studied in the context of misspecified
models where the meaningfulness of model features is uncertain, wherealse
directly comparing the distributions happens to be more natural (Bernton et
al., 2019, Frazier et al., 2020 and Legramanti et al., 2022). In particular,
we focus on the Kolmogorov-Smirnov distance ($d_{KS}$) and the Wasserstein
distance ($d_{W}$) for one-dimensional distributions, defined respectively as

$d_{KS}(P,Q)=\underset{y}{\sup}|P(y)-Q(y)|,$
$d_{W}(P,Q)=\int_{Y}|P(y)-Q(y)|dy.$

The discrepancy furnishes global indication of potential agreement between the
two distributions $P$ and $Q$, analogously to a likelihood ratio test as in
the likelihood-based inference for correctly specified models. For estimating
the distances, empirical cumulative distribution functions ($\hat{P}_{n}(y)$
and $\hat{Q}_{n}(y)$, respectively) are used.

In general models, the non-asymptotic distributions of the statistics
$d_{KS}(\cdot,\cdot)$ and $d_{W}(\cdot,\cdot)$ are complex, and numerical
methods are employed to compute exact $p$-values. Here we suggest to rely on
Algorithm 1 to derive CDs, once identifyed as summary statistics suitable
discrepancies with distribution stochastically monotone in a scalar parameter
of interest $\theta$. In particular, let us consider the observed sample
$y^{\text{obs}}$ and a fixed reference sample $y^{\text{ref}}$, drawn from a
completely known model $f(y;\theta^{\text{ref}})$. A sequence of unilateral
tests can be built by using as observed summary statistic in Algorithm 2 the
quantity
$d(y^{\text{obs}},y^{\text{ref}})=d(\hat{P}(y^{\text{obs}}),\hat{Q}(y^{\text{ref}}))$,
where $d(\cdot,\cdot)$ may be the Kolmogorov-Smirnov distance ($d_{KS}$) or
the Wasserstein distance ($d_{W}$). Then, the CD is obtained with the Accept-
Reject scheme of Algorithm 2, evaluating

$Pr_{\theta^{*}}(d(y^{*},y^{\text{ref}})>d(y^{\text{obs}},y^{\text{ref}})),$

where $y^{*}$ is simulated from the central model $y^{*}\sim f(y;\theta^{*})$.
Also, by Algorithm 3 a confidence density can be retrieved. To obtain a proper
confidence distribution, the distribution of the summary statistic should be
stochastically ordered in the parameter of interest. Hence it is convenient to
draw $y^{\text{ref}}$ from the model $f(y;\theta^{\prime})$, with
$\theta^{\prime}$ being the supremum of the proposal distribution support in
Algorithm 2.

Otherwise, a serie of bilateral tests, direclty comparing
$d(y^{*},y^{\text{obs}})$ to zero, without a reference sample, can also be
performed, for obtaining a confidence curve instead of a proper confidence
distribution (Legramanti et al., 2022).

###  5.1 Example

As in Legramanti et al. (2022) we consider a contamination study. The data
$y^{\text{obs}}=(y_{1},\ldots,y_{n})$, with sample size $n=100$, are
realizations of a Gaussian random variable $N(\theta,1)$, with nominal value
$\theta_{0}=1$. Within this setting, some scenarios of contamination are
investigated: for each one a percentage of observations is substituded with
the most extreme positive realization of a Cauchy of the same size. The amount
of contamination here is ${(5\%,10\%,15\%)}$, respectively. In particular,
using Algorithm 2 we simulated uniformly $\theta$ in $[-3,3]$ and used as a
pivot the distance $d(y^{\text{obs}},y^{\text{ref}})$, where $y^{\text{ref}}$
is drawn from $N(3,1)$.

As shown in Fig. 7, although the sample mean is dragged, the confidence
distributions remain close to and concentrated around the nominal value of 1.
In particular, the test based on the Wasserstein distance is higly stable up
to the $15\%$ of contamination. Compared to the approximate posteriors in
Legramanti et al. (2022), the CD based on Wasserstein distance seems even more
stable.

Let us denote with $\tilde{\theta}^{m}$ the confidence median and let us focus
on the Wasserstein distance. Under the non contaminated sample
($y_{\theta_{0}}$) the confidence median satisfies

$Pr(d_{W}(y_{\tilde{\theta}^{m}},y_{\theta^{ref}})>d_{W}(y_{\theta_{0}},y_{\theta^{ref}}))=0.5.$

When considering a $\epsilon$-contaminated sample
($y_{\theta_{0}^{c_{\epsilon}}}$, with $\epsilon<1$% of the data are not
generated from the assumed model), we look for $\theta^{*}$ that satisfies

$Pr(d_{W}(y_{\theta^{*}},y_{\theta^{ref}})>d_{W}(y_{\theta_{0}^{c_{\epsilon}}},y_{\theta^{ref}}))=0.5.$
(12)

The difference $\theta^{*}-\tilde{\theta}^{m}$ is the shift due to the
contamination. Writing $d_{W}(y_{\theta_{0}^{c_{\epsilon}}},y_{\theta^{ref}})$
as

$d_{W}(y_{\theta_{0}^{c_{\epsilon}}},y_{\theta^{ref}})=(1-\epsilon)d_{W}(y_{\theta_{0}},y_{\theta^{ref}})+\epsilon\cdot
d_{W}(c,y_{\theta^{ref}}),$

we can rewrite (12) as

$Pr\left(d_{W}(y_{\theta^{*}},y_{\theta^{ref}})>d_{W}(y_{\theta_{0}},y_{\theta^{ref}})+\underbrace{\epsilon[d_{W}(c,y_{\theta^{ref}})-d_{W}(y_{\theta_{0}},y_{\theta^{ref}})]}_{\Delta}\right)=0.5.$

As the term $\Delta\rightarrow 0$, the confidence median is recovered. In
particular this happens in the trivial case, when $\epsilon\rightarrow 0$ or
if $\theta^{ref}$ minimizes
$d_{W}(c,y_{\theta^{ref}})-d_{W}(y_{\theta_{0}},y_{\theta^{ref}})$, that means
it parametrizes the model which corresponds to the barycenter between the
central one and the model that generates the contamination. The optimal value
cannot be known in advance, but as an initial guess a nonrobust estimate could
be considered.

For analysing the behaviour of resulting confidence densities under the
extreme case in which the contamination amount is $\epsilon=0.2$, the data
$y^{\text{obs}}$ are still realizations of a Gaussian random variable $N(1,1)$
and for the contamination a percentage of observations is substituded with
realizations from a Cauchy. For the derivation of the CDs, we consider
different choices for the reference parameter
$\theta^{ref}=(3,4,5,6,10,20,40,100)$ (see boxplot of the data and confidence
densities in Fig. 8).

Figure 7: Confidence densities and ABC posterior based on Kolmogorv-Smirnov
(KS) and Wassserstein (W) distances; vertical lines represent the sample means
for increasing level of contamination.

Figure 8: Boxplot of the data (left) and confidence densities based on the
Wasserstein distance with $20\%$ of the data ($\epsilon=0.2$) from a
contamination model and different choices for the reference parameter. The
empirical mean (dotted line) is 2.81, while the true value is 1, and the
empirical mean of the uncontaminated sample (1-$\epsilon\;\times 100$ % of the
data) is 1.16.

##  6 Discussion

Robust statistics is an extension of classical parametric statistics that
specifically takes into account the fact that the assumed parametric models
used by the researchers are only approximate. The contribution of this paper
is to fill the gap between robust inference and confidence distributions
analyses. Indeed, 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. In particular, CDs allow to compute measures of
evidence for statements of the type ”$\psi>\psi_{0}$” or
”$\psi_{1}\leq\psi\leq\psi_{2}$”, which are of particular interest in many
real data applications, as for instance the one considered in the paper on
non-inferiority testing. Remark that the application to non-inferiority trials
discussed here can be easily extended to the comparison of other parameters
than means, such as odds ratios, hazard ratios, etc., where stability of
inferential conclusions with respect to small changes in the data or to small
model departures is essential.

The derivation of robust CDs discussed in Section 3.2, based on a frequentist
reinterpretation of ABC techniques to obtain a normalized pseudo-likelihood
function for the parameter of interest, represents a practical alternative to
other robust pseudo-likelihood functions, such as the empirical likelihood or
the quasi-likelihood (see e.g. Greco et al., 2008, and references therein).
Indeed, the aforementioned functions may present some drawbacks: the empirical
likelihood is not computable for small sample sizes and quasi-likelihoods can
be easily obtained only for scalar parameters. On the contrary, the approach
that directly incorporates robust estimating functions into ABC techniques,
with respect to available approaches based on pseudo-likelihoods, can be
computationally faster when the evaluation of the estimating function is
expensive, can be computed even for small sample sizes and for
multidimensional parameters of interest, and the derived normalized pseudo-
likelihood has the right curvature (see Ruli et al., 2020). Obviously, the
proposed method can be applied to any unbiased robust estimating equations,
such as $S$-estimating equations.

Related to the above comment, we highlight that for a scalar parameter of
interest in the presence of nuisance parameters we propose to modify the
algorithm of Ruli et al. (2020) by using the a profile estimating equation and
plugging in the value of proposals for nuisance parameters used to generate
pseudo-data. The treatment of nuisance parameters is a theme of ever-renewed
attention in frequentist inference. An interesting extension of our proposal,
not only limited to robust inference but for general estimating functions, is
the further study of ways to handling nuisance parameters with CD-based
inference, in particular with the options available with the bootstrap (see
DiCiccio and Romano, 1988) and with ABC techniques based on general profile-
type estimating functions.

As a final remark, we note that in this paper we mention the possibility to
adopt non parametric criteria and statistics, other than just centrality
measures, for deriving confidence distributions for a scalar parameter of
interest in presence of contamination. A central parametric model is assumed
but the observed data are evaluated in terms of non parametric pseudo-
distrances from a reference model, directly based on the empirical cdfs. At
present our preliminary study is limited to models with a scalar parameter of
interest, since adapting these kind of test procedures to more complex models
to deal with nuisance parameters may require information about the context,
situation-dependent considerations and also give rise to confidence curves
instead of confidence distributions.

## References

Barndorff-Nielsen, O.E., Cox, D.R. (1994). Inference and Asymptotics, CRC
Press.

Bee, M., Benedetti, R., Espa, G. (2017). Approximate maximum likelihood
estimation of the Bingham distribution. Computational Statistics & Data
Analysis, 108, 84–96.

Bortolato, E., Ventura, L. (2022). Confidence distributions and fusion
inference for intractable likelihoods. Book of Short Papers SIS 2022, Caserta
2022, to appear.

Brazzale, A.R., Davison, A.C., Reid, N. (2007). Applied asymptotics. Case-
studies in small sample statistics. Cambridge University Press.

Bernton, E., Jacob, P. E., Gerber, M., Robert, C. P. (2019). On parameter
estimation with the Wasserstein distance. Information and Inference: A Journal
of the IMA, 8, 657–676.

Carhart-Harris, R., Giribaldi, B., Watts, R., Baker-Jones, M., Murphy-Beiner,
A., Murphy, R., Nutt, D. J. (2021). Trial of psilocybin versus escitalopram
for depression. New England Journal of Medicine, 384, 1402–1411.

Chen, X., Zhou, W. (2020). Robust inference via multiplier bootstrap. The
Annals of Statistics, 48, 1665–1691.

D’Agostino, R.Bb, Massaro, J.M., Sullivan, L.M. (2003). Non- inferiority
trials: design concepts and issues the encounters of academic consultants in
statistics. Statistics in Medicine, 22, 169–186.

Dawid, A.P., Musio, M., Ventura, L. (2016). Minimum scoring rule inference.
Scand. J. Statist., 43, 123–138.

DiCiccio, T. J., Romano, J. P. (1988). A review of bootstrap confidence
intervals. Journal of the Royal Statistical Society: Series B
(Methodological), 50, 338–354.

DiCiccio, T. J., Efron, B. (1996). Bootstrap confidence intervals. Statistical
science, 11, 189-228.

Efron, B. (1979). Bootstrap methods: another look at the jackknife. The Annals
of Statistics, 7, 1–26.

Farcomeni, A., Ventura, L. (2012). An overview of robust methods in medical
research. Statistical Methods in Medical Research, 21, 111–133.

Field, C.A., Ronchetti, E. (1991). Small Sample Asymptotics. IMS Monograph
Series, Hayward (CA).

Frazier, D. T., Robert, C. P., Rousseau, J. (2020). Model misspecification in
approximate Bayesian computation: consequences and diagnostics. Journal of the
Royal Statistical Society: Series B (Statistical Methodology), 82, 421–444.

Garret, A.D. (2003). Therapeutic equivalence: fallacies and falsification.
Statistics in Medicine, 22, 741–762.

Ghosh, M., Basu, A. (2013). Robust estimation for independent non-homogeneous
observations using density power divergence with applications to linear
regression. Electronic Journal of Statististics, 7, 2420–2456.

Greco, L., Racugno, W., Ventura L. (2008). Robust likelihood functions in
Bayesian inference. Journal of Statistical Planning and Inference, 138,
1258–1270.

Hampel, F.R., Ronchetti, E.M., Rousseeuw, P.J., Stahel, W.A. (1986). Robust
Statistics. The Approach Based on Influence Functions. Wiley.

Heritier, S., Cantoni, E., Copt, S., Victoria-Feser, M.P. (2009). Robust
Methods in Biostatistics. Wiley.

Heritier, S. and Ronchetti, E.M.(1994). Robust bounded-influence tests in
general parametric models. Journal of the American Statististical Association,
89, 897–904.

Hjort, N.L., Schweder, T. (2018). Confidence distributions and related themes.
Journal of Statististical Planning and Inference, 195, 1–13.

Huber, P.J., Ronchetti, E.M. (2009). Robust Statistics. Wiley, New York.

Legramanti, S., Durante, D., Alquier, P. (2022). Concentration and robustness
of discrepancy-based ABC via Rademacher complexity, arXiv:2206.06991.

Lyddon, S. P., Holmes, C. C., Walker, S. G. (2019). General Bayesian updating
and the loss-likelihood bootstrap. Biometrika, 106, 465–478.

Muller, A. (1997). Integral probability metrics and their generating classes
of functions. Advances in Applied Probability, 29, 429–443.

Nayak, S., Bari, B. A., Yaden, D. B., Spriggs, M. J., Rosas, F., Peill, J. M.,
… , Carhart-Harris, R. (2022). A Bayesian Reanalysis of a Trial of Psilocybin
versus Escitalopram for Depression. PsyArXiv preprint

Newton, M. A., Raftery, A. E. (1994). Approximate Bayesian inference with the
weighted likelihood bootstrap. Journal of the Royal Statistical Society:
Series B (Methodological), 56, 3–26.

Reid, N. (2003). The 2000 Wald memorial lectures: asymptotics and the theory
of inference. Annals of Statistics, 31, 1695–1731.

Ronchetti, E., Ventura, L. (2001). Between stability and higher-order
asymptotics. Statistics and Computing, 11, 67–73.

Rothmann, M.D., Wiens, B.L., Chan, I.S.F (2012). Design and Analysis of Non-
Inferiority Trials. CRC Press.

Rubio, F.J., Johansen, A.M. (2013). A simple approach to maximum intractable
likelihood estimation. Electronic Journal of Statistics, 7, 1632–54.

Ruli, E., Sartori, N., Ventura, L. (2016)- Approximate Bayesian Computation
with composite score functions. Statistics and Computing, 26, 679–692.

Ruli, E., Sartori, N., Ventura, L. (2020). Robust approximate Bayesian
inference. Journal of Statistical Planning and Inference, 205, 10–22.

Ruli, E., Ventura, L. (2021). Can Bayesian, confidence distribution and
frequentist inference agree?. Statistical Methods & Applications, 30, 359–373.

Ruli, E., Ventura, L., Musio., M. (2022). Robust confidence distributions from
proper scoring rules. Statistics, 56, 455-478.

Schweder, T., Hjort, N.L. (2016). Confidence, Likelihood, Probability:
Statistical Inference with Confidence Distributions. Cambridge University
Press.

Severini, T.A. (2000). Likelihood methods in statistics. Oxford University
Press.

Soubeyrand, S., Haon-Lasportes, E., (2015). Weak convergence of posteriors
conditional on maximum pseudo-likelihood estimates and implications in ABC.
Statistics and Probability Letters, 107, 84–92.

Thornton, S., Li, W., Xie, M. (2022), Approximate confidence distribution
computing, arXiv:2206.01707.

Tsallis, C. (1988). Possible generalization of Boltzmann-Gibbs statistics.
Journal of Statistical Physics, 52, 479–487.

Varin, C., Reid, N., Firth, D. (2011). An overview of composite likelihood
methods. Statistica Sinica, 21, 5–42.

Ventura, L., Racugno, W. (2016). Pseudo-likelihoods for Bayesian inference.
In: Topics on Methodological and Applied Statistical Inference, Series Studies
in Theoretical and Applied Statistics, Springer-Verlag, 205–220.

Ventura, L., Sartori, N., Racugno, W. (2013). Objective Bayesian higher-order
asymptotics in models with nuisance parameters. Computational Statistics and
Data Analysis, 60, 90–96.

Xie, M., Singh, K. (2013). Confidence distribution, the frequentist
distribution estimator of a parameter: a review. International Statistical
Review, 81, 3–39.