Objective Priors for Invariant e-Values in the Presence of Nuisance Parameters

1. Introduction

In this article we discuss an objective matching prior, that maintains the invariance of the posterior mode when testing specific hypotheses for parametric models, especially in the presence of nuisance parameters. These parameters are often introduced to establish flexible and realistic models, although the primary focus of inference is typically limited to a parameter of interest. The proposed approach offers the advantage of eliminating the need for eliciting information on the nuisance components and for conducting multidimensional integration, and it produces invariant $e-$values in the presence of nuisance parameters.

The parametric framework that we consider can be described as follows. Consider a random sample $y=(y_1,\dots ,y_n)$ of size n from a random variable Y with parametric model $f(y;\theta )$, indexed by a paramater $\theta$, with $\theta \in \mathsf{\Theta }\subseteq \mathrm{I}{\mathrm{R}}^d$. Given a prior $\pi \left(\theta \right)$ on $\theta$, Bayesian inference for $\theta$ is based on the posterior density

$$\begin{array}{c}\hfill \pi \left(\theta \right|y)\propto \pi (\theta \left)L\right(\theta ),\end{array}$$

(1)

where $L\left(\theta \right)$ represents the likelihood function based on $f(y;\theta )$. Interest is in particular in the situation in which $\theta =(\psi ,\lambda )$, where $\psi$ is a scalar parameter for which inference is required and $\lambda$ represents the remaining $(d-1)$ nuisance parameters. In such case, Bayesian inference for $\psi$ is based on the marginal posterior density

$$\begin{array}{c}\hfill {\pi }_m\left(\psi \right|y)=\int \pi (\psi ,\lambda |y)d\lambda \propto \int \pi (\psi ,\lambda )L(\psi ,\lambda )d\lambda ,\end{array}$$

(2)

which for its computation requires both elicitation on the nuisance parameter $\lambda$ and multidimensional integration.

Asymptotic arguments are widely used in Bayesian inference through (1) and (2), based on developments of so-called higher-order asymptotics (see, e.g., [1,2,3,4,5]). Indeed, the theory of asymptotic expansions provides very accurate approximations to posterior distributions, and to various summary quantities of interest, including tail areas, credible regions and for the Full Bayesian Significance Test (see, e.g., [6,7]). Moreover, they are particularly useful for sensitivity analyses (see [8,9,10]) and also for the derivation of matching priors (see [11], and references therein). For instance, focusing on the presence of nuisance parameters, the Laplace approximation to (2) provides

$$\begin{array}{c}\hfill {\pi }_m\left(\psi \right|y)\ddot{=}\frac{1}{\sqrt{2\pi }}{\left|j_p\left(\widehat{\psi }\right)\right|}^{1/2}exp\{{\ell }_p\left(\psi \right)-{\ell }_p\left(\widehat{\psi }\right)\}\frac{|j_{\lambda \lambda }(\widehat{\psi },\widehat{\lambda }){|}^{1/2}}{|j_{\lambda \lambda }(\psi ,{\widehat{\lambda }}_{\psi }){|}^{1/2}}\frac{\pi (\psi ,{\widehat{\lambda }}_{\psi })}{\pi (\widehat{\psi },\widehat{\lambda })},\end{array}$$

(3)

where ${\ell }_p\left(\psi \right)=log{L}_p\left(\psi \right)=logL(\psi ,{\widehat{\lambda }}_{\psi })$ is the profile log-likelihood for $\psi$, with ${\widehat{\lambda }}_{\psi }$ the constrained maximum likelihood estimate (MLE) of $\lambda$ given $\psi$, $(\widehat{\psi },\widehat{\lambda })$ is the full MLE, and $j_p\left(\psi \right)=-{\partial }^2{\ell }_p\left(\psi \right)/\partial {\psi }^2$ is the profile observed information. Moreover, $j_{\lambda \lambda }(\psi ,\lambda )$ is the $(\lambda ,\lambda )-$block of the observed information from the full log-likelihood $\ell (\psi ,\lambda )=logL(\psi ,\lambda )$, and the notation $\ddot{=}$ indicates that the approximation is accurate to order $O\left(n^{-3/2}\right)$ in moderate deviation regions (see, e.g., [12], Chapter 2). One appealing feature of higher-order approximations like (3) is that they may routinely be applied in practical Bayesian inference, since they require little more than standard likelihood quantities for their implementation, and hence they may be available at little additional computational cost over simple first-order approximations.

In the presence of nuisance parameters, starting from approximation (3) it is possible to define a general posterior distribution for $\psi$ of the form

$$\begin{array}{c}\hfill {\pi }^{*}\left(\psi \right|y)\propto {\pi }^{*}\left(\psi \right)L_p\left(\psi \right),\end{array}$$

(4)

where ${\pi }^{*}\left(\psi \right)$ is now a prior distribution on $\psi$ only. Bayesian inference based on pseudo-likelihood functions – i.e. functions of $\psi$ only and of the data y with properties similar to those of a genuine likelihood function, such as the profile likelihood – has been widely used and discussed in the recent statistical literature. Moreover, it has been theoretically motivated in several papers (see, for instance, [13,14,15,16], and references therein), also focusing on the derivation of suitable objective priors. Especially when the dimension of $\lambda$ is large, there are two advantages in using (4) instead of the marginal posterior distribution (2). First, the elicitation over $\lambda$ is not necessary and, second, the computation of the integrals in (2) is circumvented.

Focusing on (4), in this paper it is of interest to test the precise (or sharp) null hypothesis

$$\begin{array}{c}\hfill H_0:\psi ={\psi }_0\mathrm{against}{H}_1:\psi \ne {\psi }_0\end{array}$$

(5)

using the measure of evidence for the Full Bayesian Significance Test (see, e.g., [6,7]). The Full Bayesian Significance Test (FBST) quantifies evidence by considering the posterior probability associated with the least probable points in the parameter space under $H_0$. Higher-order asymptotic computation of the FBST for precise null hypotheses in the presence of nuisance parameters has been discussed in [2].

The original measure of evidence for the FBST is not invariant under suitable transformations of the parameter, a property which has been reached however in the more recent definition of the $e-$value (see [17,18], and references therein). Neverthless, when working on a scalar parameter of interest, in the presence of nuisance parameters, the $e-$value is not invariant with respect to marginalisations of the nuisance parameter and it must be used in the full dimensionality of the parameter space. This requires elicitation on the complete parameters, numerical optimization and numerical integration, that can be computationally heavy especially when the dimension of $\lambda$ is large.

The aim of this paper is to consider the $e-$value in the context of the pseudo-posterior distribution ${\pi }^{*}\left(\psi \right|y)$, suggesting in this respect a suitable objective prior ${\pi }^{*}\left(\psi \right)$ to be used in (4). More precisely, focus is on a particular matching prior, which ensure invariance of the posterior mode of the pseudo-posterior distribution. As a consequence also Highest Probability Density credible (HPD) sets are invariant, as well as the $e-$value.

The paper is organized as follows. Section 2 provides a short review on the FBST for testing precise null hypotheses and also illustrates asymptotic approximations for the $e-$value, extending results of [2]. Section 3 discusses the derivation of the objective matching prior for the parameter of interest only, called median matching prior, that produces invariant $e-$values. Also several targeted and challenging examples are discussed. Finally, Section 4 closes with some concluding remarks.

2. The FBST measure of evidence

Suppose that we need to decide between two hypotheses: the null $H_0$ and the alternative $H_1$. The usual Bayesian testing procedure is based on the well-known Bayes factor (BF), defined as the ratio of the posterior odds to the prior odds in favor of the null hypothesis. A high BF or its logarithm suggest evidence in favor of $H_0$. However, it is well-known that, when improper priors are used, the BF can be undetermined and, when the null hypothesis is precise (as specified in (5)), the BF can lead to the so-called Jeffreys-Lindley’s paradox (see, e.g. [19]). Moreover, the BF is not calibrated, i.e. its finite sampling distribution is unknown and it may depend on the nuisance parameter.

To avoid these drawbacks, in recent years an alternative Bayesian procedure, called FBST, has been introduced by [7] in case of sharp hypothesis $H_0$ identified by the null set ${\mathsf{\Theta }}_0$, a submanifold of $\mathsf{\Theta }$ of lower dimension. The FBST quantifies evidence by considering the posterior probability associated with the least probable points in the parameter space ${\mathsf{\Theta }}_0$. When this probability is high, it favors the null hypothesis, providing a clear and interpretable measure of support for $H_0$ (see, e.g. [6,18,20], and references therein). The FBST is based on a specific loss function, and thus the decision made under this procedure is the action that minimizes the corresponding posterior risk.

The FBST operates by determining the $e-$value, a representation of Bayesian evidence associated to $H_0$. To construct the $e-$value, the authors introduced the posterior surprise function and its supremum given, respectively, by

$${\pi }_s\left(\theta \right|y)=\frac{\pi \left(\theta \right|y)}{r\left(\theta \right)}\mathrm{and}{s}^{*}={\pi }_s\left({\theta }^{*}\right|y)=\underset{\theta \in {\mathsf{\Theta }}_0}{sup}{\pi }_s\left(\theta \right|y),$$

where $r\left(\theta \right)$ is a suitable reference function to be chosen. Then, they introduce the tangential set $T_y\left({\theta }^{*}\right)$ defined as the set of parameter values for which the posterior surprise function exceeds the supremum $s^{*}$, that is

$$T_y\left({\theta }^{*}\right)=\{\theta \in \mathsf{\Theta }:{\pi }_s\left(\theta \right|y)>s^{*}\}.$$

This set, often referred to as the Highest Relative Surprise Set, includes parameter values with higher surprise than those within the null set ${\mathsf{\Theta }}_0$. The $e-$value is then computed as

$$ev=1-{\int }_{T_y\left({\theta }^{*}\right)}{\pi }_s\left(\theta \right|y)d\theta ,$$

and $H_0$ is rejected for small values of $ev$.

The original FBST, as proposed by [7,21], relies on a flat reference function $r\left(\theta \right)\propto 1$, so that this first version involved the determination of the tangential set $T_y\left(\theta \right)$ starting only from the posterior distribution $\pi \left(\theta \right|y)$. However, this initial version lacked invariance under reparameterizations. Subsequent refinements of the FBST introduced the importance of reference density functions, making the $e-$value explicitly invariant under appropriate transformations of the parameter. Common choices for the reference function include uninformative priors, like the uniform distribution, maximum entropy densities, or Jeffreys’ invariant prior. In [22], the use of the Jeffreys’ prior, $\pi \left(\theta \right)\propto {\left|i\left(\theta \right)\right|}^{1/2}$, where $i\left(\theta \right)$ is the Fisher information derived from $L\left(\theta \right)$, is discussed as the reference function to derive invariant HPD sets and Maximum A Posteriori (MAP) estimators that are invariant under reparameterizations. Note that the $ev$ uses the full dimensionality of the parameter space. Moreover, this measure is not invariant with respect to transformations of the nuisance parameters and the use of high posterior densities to construct credible sets may produce inconsistencies.

Concerning the asymptotic behavior of the $ev$ it can be proven that, under suitable regularity conditions as the sample size increases, with ${\theta }_0$ representing the true parameter value (see [18]), it holds:

• If $H_0$ is false, i.e. ${\theta }_0\notin H_0$ , then $ev$ converges in probability to 1.
• If $H_0$ is true, i.e. ${\theta }_0\in H_0$ , then, denoting by $V\left(c\right)=Pr(ev\le c)$ the cumulative distribution function of $ev$, we have that $V\left(c\right)\approx Q(d-h,Q^{-1}(d,c))$, with $d=dim\left(\mathsf{\Theta }\right)$, $h=dim\left({\mathsf{\Theta }}_0\right)$ and $Q(k,x)$ the cumulative chi−square distribution with k degrees of freedom.

In practice, the computation of $ev$ is performed in two steps: (a) a numerical optimization and (b) a numerical integration. The numerical optimization step consists in finding the maximizer ${\theta }^{*}$ of ${\pi }_s\left(\theta \right|y)$ under the null hypothesis. The numerical integration step consists of integrating the posterior surprise function over the region where it is greater than ${\pi }_s\left({\theta }^{*}\right|y)$, to obtain the $e-$value. These computational steps make the FBST a computationally intensive procedure. Despite efficient computational algorithms for local and global optimization, as well as numerical integration, obtaining precise results for hypotheses like (5) is highly demanding, especially with large nuisance parameter dimensions. Numerical integration can be tackled by resorting to higher-order tail area approximations, as reviewed in the Bayesian framework in [4]. An application of asymptotic approximation to the FBST in its first formulation, i.e. with reference function $r\left(\theta \right)\propto 1$, has been discussed in [2].

2.1. Asymptotic approximations for the $e-value$

A first-order approximation for the $e-$value, when testing (5), is simply given by (see, e.g., [23,24])

$$\begin{array}{c}\hfill ev\dot{=}2\left(1-\mathsf{\Phi }\left(\left|\frac{{\psi }_0-\widehat{\psi }}{\sqrt{j_p{\left(\widehat{\psi }\right)}^{-1}}}\right|\right)\right),\end{array}$$

(6)

where the symbol "$\dot{=}$" indicates that the approximation is accurate to $O\left(n^{-1/2}\right)$ and $\mathsf{\Phi }(·)$ is the standard normal distribution function. Thus, to first-order, $ev$ agrees with the $p-$value based on the profile Wald statistic

$$\begin{array}{c}\hfill w_p\left(\psi \right)=\frac{(\widehat{\psi }-{\psi }_0)}{\sqrt{j_p{\left(\widehat{\psi }\right)}^{-1}}}.\end{array}$$

(7)

In practice, the approximation (6) of $ev$ may be inaccurate, in particular when the dimension of $\lambda$ is large with respect to the sample size, because it forces the marginal posterior distribution to be symmetric.

The practical computation of $ev$ requires the evaluation of integrals of the marginal posterior distribution. In order to have more accurate evaluations of $ev$, it may be useful to resort to higher-order asymptotics based on tail area approximations (see, e.g., [3,4], and references therein). Indeed, the measure of evidence involves integrals of the marginal surprise posterior density ${\pi }_{ms}\left(\psi \right|y)$. In particular, extending the application of the tail area argument to the marginal surprise posterior density, we can derive a $O\left(n^{-3/2}\right)$ approximation to the marginal surprise posterior tail area probability, given by

$$\begin{array}{c}\hfill {\int }_{{\psi }_0}^{\infty }{\pi }_{ms}\left(\psi \right|y)d\psi \ddot{=}\mathsf{\Phi }\left(r_B^{*}\left({\psi }_0\right)\right),\end{array}$$

(8)

where

$$r_B^{*}\left(\psi \right)=r_p\left(\psi \right)+\frac{1}{r_p\left(\psi \right)}log\frac{q_B\left(\psi \right)}{r_p\left(\psi \right)},$$

with

$$r_p\left(\psi \right)=\mathrm{sign}(\widehat{\psi }-\psi ){\left[2({\ell }_p\left(\widehat{\psi }\right)-{\ell }_p\left(\psi \right))\right]}^{1/2}$$

profile likelihood root and

$$q_B\left(\psi \right)={\ell }_p^{\prime }\left(\psi \right){\left|j_p\left(\widehat{\psi }\right)\right|}^{-1/2}\frac{|j_{\lambda \lambda }(\psi ,{\widehat{\lambda }}_{\psi }){|}^{1/2}}{|j_{\lambda \lambda }(\widehat{\psi },\widehat{\lambda }){|}^{1/2}}\frac{\pi (\widehat{\psi },\widehat{\lambda })}{\pi (\psi ,{\widehat{\lambda }}_{\psi })}\frac{r(\psi ,{\widehat{\lambda }}_{\psi })}{r(\widehat{\psi },\widehat{\lambda })}.$$

In the expression of $q_B\left(\psi \right)$, ${\ell }_p^{\prime }\left(\psi \right)=\partial {\ell }_p\left(\psi \right)/\partial \psi$ is the profile score function.

Using the tail area approximation (8), a third-order approximation of the measure of evidence $ev$ can be derived. The approximation, assuming without loss of generality that ${\psi }_0$ is smaller than the MAP of ${\pi }_{ms}\left(\psi \right|y)$, is given by

$$\begin{array}{c}\hfill ev\left(\psi \right)\ddot{=}1-\mathsf{\Phi }\left(r_B^{*}\left({\psi }_0\right)\right)+\mathsf{\Phi }\left(r_B^{*}\left({\psi }_0^{*}\right)\right),\end{array}$$

(9)

with ${\psi }_0^{*}$ the value of the parameter such that ${\pi }_{ms}\left({\psi }_0^{*}\right|y)={\pi }_{ms}\left({\psi }_0\right|y)$. Note that

$$\mathsf{\Phi }\left(r_B^{*}\left({\psi }_0\right)\right)-\mathsf{\Phi }\left(r_B^{*}\left({\psi }_0^{*}\right)\right)\ddot{=}{\int }_{{\psi }_0^{*}}^{{\psi }_0}{\pi }_{ms}\left(\psi \right|y)d\psi =1-ev$$

in (9) gives the posterior probability of the HPD credible interval $({\psi }_0,{\psi }_0^{*})$. Note also that the higher-order approximation (9) does not call for any condition on the prior $\pi (\psi ,\lambda )$, i.e. it can be also improper. Finally, when ${\pi }_{ms}\left(\psi \right|y)$ is symmetric, Equation (9) reduces to $ev\ddot{=}2(1-\mathsf{\Phi }\left(r_B^{*}\left({\psi }_0\right)\right))$.

While tail area approximations require little more than standard likelihood quantities for their implementation and, in this respect, they are available at little additional computational cost over the first-order approximation, they require elicitation on the complete parameter $\theta$ and to choose the reference function $r\left(\theta \right)$.

3. An invariant objective prior

The aim of this section is to derive a default prior ${\pi }^{*}\left(\psi \right)$ to be used in (4). To this end, following [11] we use the shrinkage argument, which is a crucial procedure in the development of matching priors, i.e. priors that ensure, up to the desired order of asymptotics, an agreement between Bayesian and frequentist procedures. Examples of matching priors are (see [11]) for posterior quantiles, for credible regions and for prediction. Here, we focus on a specific matching prior, that ensures the invariance of the posterior mode in the posterior distribution (4). As a consequence, the invariance extends to HPDs, as well as the $e-$value, achieved incorporating the reference function within the prior.

The proposed choice of the prior ${\pi }^{*}\left(\psi \right)$, that makes the MAP and thus also HPDs and the $e-$value invariant under 1-1 reparameterization, will depend on the log-likelihood $\ell \left(\theta \right)$ and on its derivatives. In regular parametric estimation problems, both the MLE and the score estimating function exhibit an asymptotically symmetric distribution centered at the true parameter value and at zero, respectively. However, these asymptotic behaviours may poorly reflect exact sampling distributions particularly in cases with small or moderate sample information, sparse data, or complex models. Several proposals have been developed to correct the estimate or the estimating function. Most available methods are aimed at approximate bias adjustment, either of the MLE or of the profile score function, also when nuisance parameters are present (see [25] for a review of bias reduction for the MLE and [26] and subsequent literature for bias correction of the profile score). Lack of equivariance impacts the so-called implicit bias reduction methods, which achieve first-order bias correction by modifying the score equation (see [25,27]). To avoid this drawback, in this paper we focus on the median modification of the score, or profile score equation, whose solution respects equivariance under monotone reparameterizations ([28]). Similar to Firth’s implicit method ([27]), the median modification of the score, or profile score, does not rely on finiteness of the MLE, thereby effectively preventing infinite estimates.

In practice, to derive the median matching prior prior ${\pi }^{*}\left(\psi \right)$, we impose that the MAP of ${\pi }^{*}\left(\psi \right|y)$ coincides with a refined version of the MLE, obtained as the solution of the median modified score function ([28]). To introduce this new invariant prior, we initially explore the scenario without nuisance parameters and then the situation in which nuisance parameters are present.

3.1. No nuisance parameters

Let’s explore first the scenario where $\theta$ is scalar. In order to obtain median bias reduction of the MLE, it is possible to resort to a modified version of the score function of the form

$$\begin{array}{c}\hfill t\left(\theta \right)={\ell }_{\theta }\left(\theta \right)+m\left(\theta \right),\end{array}$$

(10)

where ${\ell }_{\theta }\left(\theta \right)={\ell }_{\theta }(\theta ;y)=\partial \ell (\theta ;y)/\partial \theta$ is the score function and $m\left(\theta \right)$ a suitable correction term of order $O\left(1\right)$. In particular, the median modified score function assumes for $m\left(\theta \right)$ the expression

$$m\left(\theta \right)=-\frac{E\left({\ell }_{\theta }{\left(\theta \right)}^3\right)}{6i\left(\theta \right)}.$$

The solution $\tilde{\theta }$ to the equation $t\left(\theta \right)=0$ not only upholds equivariance under componentwise monotone reparameterizations but also approximates median unbiasedness ([28]). Note that likelihood inference based on (10) does not depend explicitely on the MLE. Indeed, the modified score function has been found to overcome infinite estimate problems. Likewise the MLE, also $\tilde{\theta }$ is asymptotically $N(\theta ,i{\left(\theta \right)}^{-1})$, so that the Wald-type statistics only differ in location.

Since Bayes’ theorem is a statement of adittivity on the log scale $log\pi \left(\theta \right|y)=log\pi (\theta )+logL(\theta )+$constant, we observe that in the Bayesian framework $m\left(\theta \right)$ can be interpreted as the derivative of the logarithm of a prior, that is $m\left(\theta \right)=\partial log\pi \left(\theta \right)/\partial \theta$. We are thus looking for a matching prior ${\pi }^{*}\left(\theta \right)$ such that

$$\frac{\partial log{\pi }^{*}\left(\theta \right)}{\partial \theta }=-\frac{E\left({\ell }_{\theta }{\left(\theta \right)}^3\right)}{6i\left(\theta \right)}.$$

In the scalar parameter case, it is straightforward to show that the proposed median matching prior takes the form

$$\begin{array}{ccc}\hfill {\pi }^{*}\left(\theta \right)& \propto & exp\left(-\frac{1}{6}\int i{\left(\theta \right)}^{-1}E\left({\ell }_{\theta }{\left(\theta \right)}^3\right)d\theta \right)\hfill \\ & \propto & exp\left(\frac{1}{6}\int i{\left(\theta \right)}^{-1}(3E\left({\ell }_{\theta \theta }\left(\theta \right){\ell }_{\theta }\left(\theta \right)\right)+E\left({\ell }_{\theta \theta \theta }\left(\theta \right)\right))d\theta \right),\hfill \end{array}$$

where ${\ell }_{\theta \theta }\left(\theta \right)=\partial {\ell }_{\theta }\left(\theta \right)/\partial \theta$ and ${\ell }_{\theta \theta \theta }\left(\theta \right)=\partial {\ell }_{\theta \theta }\left(\theta \right)/\partial \theta$. The posterior based on the median matching prior is thus

$${\pi }^{*}\left(\theta \right|y)\propto exp\left(\ell \left(\theta \right)-\frac{1}{6}\int i{\left(\theta \right)}^{-1}E\left({\ell }_{\theta }{\left(\theta \right)}^3\right)d\theta \right).$$

A first-order approximation for the $e-$value, when testing $H_0:\theta ={\theta }_0$, is simply given by

$$\begin{array}{c}\hfill ev\dot{=}2\left(1-\mathsf{\Phi }\left(\left|\frac{{\theta }_0-\tilde{\theta }}{\sqrt{i{\left({\theta }_0\right)}^{-1}}}\right|\right)\right),\end{array}$$

(11)

which differs in location with respect to the classical first-order approximation for the $e-$value based on the MLE. A second approximation for the $e-$value, when testing $H_0:\theta ={\theta }_0$, can be obtained from the asymptotic distribution of the modified score function (10), that is

$$\begin{array}{c}\hfill ev\dot{=}2\left(1-\mathsf{\Phi }\left(\left|\frac{t\left({\theta }_0\right)}{\sqrt{i\left({\theta }_0\right)}}\right|\right)\right).\end{array}$$

(12)

Although the first-order equivalence between (11) and (12), note that (11) is based on an easily understandable comparison between estimated value and hypothetical value, taking estimation error into account, and is widely used in applications, but does not satisfy the principle of parameterization invariance. On the other hand, $t\left(\theta \right)/\sqrt{i\left(\theta \right)}$ is parameterization invariant.

Note that, when using a predictive matching prior, i.e. a prior ensuring asymptotic equivalence of higher-order frequentist and Bayesian predictive densities (see, e.g., [11]), the term $m\left(\theta \right)$ in (10) corresponds to the Firth’s adjustment ([27])

$$m_F\left(\theta \right)=-\frac{(E\left({\ell }_{\theta }{\left(\theta \right)}^3\right)+E\left({\ell }_{\theta \theta }\left(\theta \right){\ell }_{\theta }\left(\theta \right)\right))}{2i\left(\theta \right)}.$$

In view of this, for general regular models, Firth’s estimate coincides with the mode of the posterior distribution obtained using the default predictive matching prior. However, lack of invariance affects this kind of adjustment ([25]), unless dealing with linear transformations.

Example 1: 

One parameter exponential family. For a one parameter exponential family with canonical parameter $\theta$, i.e. with density

$$f(y;\theta )=exp\left\{\theta a\right(y)-K(\theta \left)\right\}b\left(y\right),$$

the median modified score function has the form

$$t\left(\theta \right)={\ell }_{\theta }\left(\theta \right)+\frac{K_{\theta \theta \theta }}{6K_{\theta \theta }},$$

where $K_{\theta \theta \theta }={\partial }^{3}K\left(\theta \right)/\partial {\theta }^3$ and $K_{\theta \theta }={\partial }^{2}K\left(\theta \right)/\partial {\theta }^2=i\left(\theta \right)$. In this parameterization, $t\left(\theta \right)$ can be seen as the first derivative of the log-posterior

$$log\pi \left(\theta \right|y)=\ell (\theta )+logi(\theta )/6.$$

On the other hand, Firth’s modified score takes the form $t_F\left(\theta \right)={\ell }_{\theta }\left(\theta \right)+K_{\theta \theta \theta }/\left(2K_{\theta \theta }\right)$. The effect of the median modification is to consider the median matching prior ${\pi }^{*}\left(\theta \right)\propto i{\left(\theta \right)}^{1/6}$, while $t_F\left(\theta \right)$ implies a Jeffreys’ prior ${\pi }_J\left(\theta \right)\propto i{\left(\theta \right)}^{1/2}$. Note that, for a one parameter exponential family with canonical parameter $\theta$, both ${\pi }^{*}\left(\theta \right)$and ${\pi }_J\left(\theta \right)$ belong to the family of invariant priors discussed in [29,30].

Example 2: 

Scale model. Consider the scale model $f(y;\theta )=(1/\theta )p_0(y/\theta )$, where $p_0(·)$ is a given function. Let $g(·)=-log{p}_0(·)$. We have $E\left({\ell }_{\theta }^3\right)=c_1/{\theta }^3$, $E\left({\ell }_{\theta }^3{\ell }_{\theta }\right)=c_2/{\theta }^3$ and $i\left(\theta \right)=c_3/{\theta }^2$, with $c_1=\int (y^3{g}^{\prime \prime \prime }\left(y\right)+6y^2{g}^{\prime \prime }\left(y\right)+6y{g}^{\prime }\left(y\right)-2)p_0\left(y\right)dy$, $c_2=\int (3y{g}^{\prime }\left(y\right)+y^2{g}^{\prime \prime }\left(y\right)-2y^2{g}^{\prime }{\left(y\right)}^2-y^3{g}^{\prime }\left(y\right)g^{\prime \prime }\left(y\right)-1)p_0\left(y\right)dy$ and $c_3=\int (2y{g}^{\prime }\left(y\right)+y^2{g}^{\prime \prime }\left(y\right)-1)p_0\left(y\right)dy$. The median matching prior is thus ${\pi }^{*}\left(\theta \right)\propto {\theta }^{-c_1/6c_3}$, while the Jeffreys’ prior for a one-parameter scale model is ${\pi }_J\left(\theta \right)\propto {\theta }^{-1}$ and the prior associated to the Firth’s adjustment is ${\pi }_F\left(\theta \right)\propto {\theta }^{-(c_1+c_2)/2c_3}$.

Example 3: 

Skew-normal distribution. Consider a skew-normal distribution with shape parameter $\theta \in \mathrm{I}\mathrm{R}$, with density $f(y;\theta )=2\varphi \left(y\right)\mathsf{\Phi }\left(y\theta \right)$, where $\varphi (·)$ is the standard normal density function. The median correction term for the score function associated to the median matching prior is (see [28,31])

$$m\left(\theta \right)=\frac{E(y^3\varphi {\left(y\theta \right)}^3/\mathsf{\Phi }{\left(y\theta \right)}^3)}{6E(y^2\varphi {\left(y\theta \right)}^2/\mathsf{\Phi }{\left(y\theta \right)}^2)}.$$

Numerical integration must be performed to obtain the expected values involved in $m\left(\theta \right)$.

In order to illustrate the proposed prior, we consider draws from the skew-normal model with true parameter ${\theta }_0=3$ and increasing sample sizes $n=20,30,50,200$ (Figure 1). The posterior distributions are obtained with the method by [32], i.e. drawing ${10}^5$ values and accepting the best $5\%$. The $e-$values associated to the null (true) hypothesis $H_0:\theta =3$ and the (false) hypothesis $H_0:\theta =4$ are reported in Table 1. For comparison, also the Jeffreys’ prior ([33]), the predictive matching prior ([31]) and the flat prior, with uniform reference function, are considered. Progressive agreement among evidence values obtained with the proposed median matching prior and the other priors for larger sample size is shown. Also, as expected, progressively increasing n, the evidence values indicate agreement with the true hypothesis $H_0:{\theta }_0=3$ and disagreement with $H_0:{\theta }_0=4$ for all the priors used. Anyway, note that the posterior distribution obtained with a flat prior, and a uniform reference function, is proportional to the likelihood function that can be monotone. In view of this, while the MAPs of the posterior based on the default priors are always finite, in some samples the MAP of the posterior with the non-informative prior may be infinite. An example of this effect is illustrated in Figure 2.

The properties of first-order approximations of the $e-$values have been investigated by a simulation study, with sample sizes $n=20,30,50,200$. Results are displayed in Figure 3. Distributions of the $e-$value from the posterior based on the median matching prior are better, both for small and moderate sample sizes, in terms of convergence to the Uniform distribution. Moreover, score-type $e-$values (12) are also preferable than Wald-type $e-$values (11). For the results with the posterior distribution obtained with a flat prior we found $4.3\%,4.2\%,0.9\%,0\%$ of infinite estimates for the sample sizes considered in the simulation study, and in these cases the $e-$value was considered as 0.

3.2. Presence of nuisance parameters

In the presence of nuisance parameters, in order to obtain median bias reduction of the MLE, it is possible to resort to a modified version of the profile score function of the form

$$\begin{array}{c}\hfill t_p\left(\psi \right)={\ell }_p^{\prime }\left(\psi \right)+m(\psi ,{\widehat{\lambda }}_{\psi }),\end{array}$$

(13)

where $m(\psi ,\lambda )$ a suitable correction term of order $O\left(1\right)$. In particular, for the median modified profile score function, the adjustment $m(\psi ,\lambda )$ assumes the expression

$$m(\psi ,\lambda )=-{\kappa }_{1\psi }+\frac{{\kappa }_{3\psi }}{6{\kappa }_{2\psi }},$$

where ${\kappa }_{1\psi }$, ${\kappa }_{2\psi }$ and ${\kappa }_{3\psi }$ are the first three cumulants of ${\ell }_p^{\prime }\left(\psi \right)$ (see [28], Section 2.2, for their expression). For the estimator ${\tilde{\psi }}_p$, defined as the solution of $t_p\left(\psi \right)=0$, parameterization equivariance holds under interest respecting reparameterizations ([28]).

Note that, also in the context of nuisance parameters, we are in the situation in which the proposed prior ${\pi }^{*}\left(\psi \right)$ is known through its first derivative; this is typically the situation with matching priors (see, e.g., [11]). Since the parameter of interest is scalar, the posterior based on the median matching prior can be written as

$$\begin{array}{c}\hfill {\pi }^{*}\left(\psi \right|y)\propto exp\left({\ell }_p\left(\psi \right)+\int m(\psi ,{\widehat{\lambda }}_{\psi })d\psi \right).\end{array}$$

(14)

A simple analytical way of approximating to first-order the posterior distribution (14) based on the median matching prior is to resort to a quadratic form of $t_p\left(\psi \right)$. In particular, the approximate posterior distribution for $\psi$ takes the form

$$\begin{array}{c}\hfill {\pi }^{*}\left(\psi \right|y)\dot{\propto }exp\left(-\frac{1}{2}{s}_p(\psi ;y)\right),\end{array}$$

(15)

where $s_p\left(\psi \right)=t_p{\left(\psi \right)}^2{j}_p{\left(\psi \right)}^{-1}$ is a Rao score-type statistic based on (13) and the symbol "$\dot{\propto }$" means asymptotic proportionality to first-order. In this case, a first-order approximation of the $e-$value, when testing $H_0:\psi ={\psi }_0$, is given by

$$\begin{array}{c}\hfill ev\dot{=}2\left(1-\mathsf{\Phi }\left(\left|\frac{t_p\left({\psi }_0\right)}{\sqrt{j_p\left({\psi }_0\right)}}\right|\right)\right).\end{array}$$

(16)

In this case, an higher order approximation via (3) would be impractical since a closed form prior is not available. As an alternative, simulation-based approaches may be used to derive the implied posterior distribution (14) based on the median matching prior. The first one relies on Approximate Bayesian Computation (ABC) techniques, using ${\tilde{\psi }}_p$ or the modified profile score function $t_p\left(\psi \right)$ as summary statistics; see [34] for the modification of the algorithm of [32] by using a profile estimating equation. This first method introduces an approximation at the level of the posterior estimation. The second one still relies on (13) but considers use of Manifold MCMC methods (see ,e.g., [35]) to conditioning exactly on the profile equation and not up to a tolerance level, as in ABC (see [36,37]). The algorithm moves on the constrained space $\{(y,\psi )\in \mathcal{Y}\times \Psi |t_p\left({\tilde{\psi }}_p\right)=0\}$, where ${\tilde{\psi }}_p$ is the solution of the estimating equation on the original data. For the latter method, we need minimal regularity assumptions on $m(\psi ,\lambda )$, which is assumed to be continuous, differentiable and available in closed form expression. Note, for instance, that in the skew-normal example in Section 3.1 these conditions are not met.

Example 4:

Exponential family. If $f(y;\theta )$ is an exponential family of order d with canonical parameter $(\psi ,\lambda )$, i.e. $f(y;\psi ,\lambda )=exp\{\psi t\left(y\right)+{\lambda }^{T}s\left(y\right)-K(\psi ,\lambda )\}h\left(y\right)$, quantities involved in $m(\psi ,\lambda )$ are simply obtained from derivatives of $K(\psi ,\lambda )$ ([28]). Note that, in this framework, ${\ell }_p^{\prime }\left(\psi \right)-{\kappa }_{1\psi }$ is an approximation with error of order $O\left(n^{-1}\right)$ of the score for $\psi$ in the conditional model given $s\left(y\right)$ (see e.g. [38],Section 10.2). Then, in the continuous case, the MAP ${\tilde{\psi }}_p$ is an approximation of the optimal conditional median unbiased estimator, and ${\pi }^{*}\left(\psi \right|y)$ is related to the conditional likelihood for $\psi$ given by $L_c\left(\psi \right)=exp(\psi t\left(y\right)-K_s\left(\psi \right))$; see [39] for a Bayesian interpretation of such pseudo-likelihoods.

Example 5:

Multivariate regression model. Consider a regression model of the form

$$Y_{ij}={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+{ϵ}_{ij},i=1,\dots ,n,j=1,2,$$

where it is assumed that ${ϵ}_i\sim N_2(0,\mathsf{\Sigma })$, with $\mathsf{\Sigma }={\sigma }^2\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right)$ positive definite matrix, and $({\beta }_0,{\beta }_1,{\beta }_2,{\sigma }^2,\rho )$ are unknown parameters. This model is widely used for instance in time series analysis, in which as regression covariates the past of the observables y are used. We focus on the problem of testing hypothesis on the correlation coefficient $\rho$.

Consider a draw with true parameter ${\rho }_0=0.95$ and $n=20$. For obtaining the proposed posterior (14) for $\rho$ we first compute the MAPs with the median matching prior and also, for comparison, with the predictive matching prior, which are respectively 0.953 and 0.92. Note that the expression of the predictive matching prior for $(\psi ,\lambda )$ is given in [40] and it still corresponds to the Firth’s adjustment to the score function. The expressions of the modified profile estimating functions $t_p\left(\psi \right)$ and $t_F\left(\psi \right)$ are obtained from [41] and are available in closed form expressions. Hence, the Manifold MCMC method can be used to obtain the implied posteriors, whose approximation is comparable to that of any MCMC sampler. In particular we used 20000 iterations.

We compare the posterior distribution based on the proposed median matching prior, with those obtained with the predictive matching prior and with an inverse-Wishart prior for the covariance matrix $\mathsf{\Sigma }$ with one degree of freedom and identity position, and uniform prior on the regression parameters. The posterior distributions are displayed in Figure 4. The hypothesis of interest is $H_0:\rho =0.9$, and a smaller e-value indicating disagreement with the hypothesis should be preferable. The $e-$values are 0.25 with the median matching prior, 0.36 for the predictive matching prior and 0.60 with the inverse-Wishart prior. Note that the $e-$value based on the inverse-Wishart prior involves the constrained maximization and multidimensional integration thus is not directly readable in Figure 4. Indeed, one crucial difference is that the original $e-$value formulation links the evidence of the null hypothesis to the evidence of a more refined hypothesis, choosing the MAP under the null hypothesis for all the nuisance parameters, while in the alternative (tangential) set all values are used, and integration is performed on the full dimensionality of the space. On the contrary, in the proposed posterior based on the median matching prior, the maximizer of nuisance parameters are taken both in the null and non null sets.

Finally, for the posterior based on the inverse-Wishart prior, we also computed the $e-$value based on high-order tail area approximation (9) of the marginal surprise posterior, which is equal to 0.27. This procedure still avoids the multidimensional integration but the result is not invariant to changes of parametrization.

Example 6:

Logistic regression model. Let $y_i$, $i=1,\dots ,n$, be independent realizations of binary random variables with probability ${\pi }_i$, where $log({\pi }_i/(1-{\pi }_i))={\eta }_i=x_i\beta$ and $x_i=(x_{i1},\dots ,x_{ip})$ is a row vector of covariates. We assume that a generic scalar component of $\beta$ is of interest and we treat the remaining components as nuisance parameters.

As an example, we consider the endometrial cancer grade dataset analyzed, among others, in [42]. The aim of the clinical study was to evaluate the relationship between the histology of the endometrium (HG), the binary response variable, of $n=79$ patients and three risk factors: 1. Neovascularization (NV), that indicates the presence or extent of new blood vessel formation; 2. Pulsatility Index (PI), that measures blood flow resistance in the endometrium; 3. Endometrium Height (EH), that indicates the thickness or height of the endometrium. A logistic model for HG, including an intercept and using all the covariates (NV, PI, EH), has been fitted, but maximum likelihood leads to infinite MLE of the coefficient ${\beta }_2$ related to NV, due to quasi-complete separation. This phenomenon prohibits the use of diffuse priors for ${\beta }_2$, since the corresponding posterior wouldn’t concentrate. Moreover, the $e-$value with non-informative priors cannot be obtained also for any hypothesis concerning parameters different from ${\beta }_2$.

If we consider ${\beta }_2$ as the parameter of interest, while the remaining regression coefficients are treated as nuisance parameters, the analysis with the median matching prior allows to obtain a global proper posterior, with MAP equal to 3.86, open to interpretation both in the original scale and in terms of odds ratios. Similarly, the posterior based on the predictive matching prior, which in this model coincides with Jeffreys’ prior $\pi \left(\beta \right)\propto {\left|i\left(\beta \right)\right|}^{1/2}$ is proper, with the MAP set at 2.92. The latter suffers from lack of interpretability on different scales, since a different parametrization in estimation phase would affect the results.

If we consider ${\beta }_3$ as the parameter of interest, related to the risk factor PI, the MAPs are -0.038 when using the median matching prior and -0.035 when using the predictive matching prior. The $e-$values for the hypothesis ${\mathcal{H}}_0:{\beta }_3=0$ are 0.60 and 0.55, respectively (see Figure 5). Likewise, the interpretation of e-values remains consistent and independent of parametrization solely in the first case.

4. Conclusions

Although (14) cannot always be considered as orthodox in a Bayesian setting, the use of alternative likelihoods is nowadays widely shared, and several papers focus on the Bayesian application of some well-known pseudo-likelihoods. In particular, the proposed posterior ${\pi }^{*}\left(\psi \right|y)$ has the advantages of avoiding the elicitation on the nuisance parameter $\lambda$ and of the computation of multidimensional integrations. Moreover, it provides invariant MAPs, HPDs and $e-$values, without the adoption of a reference function . The use of matching priors, i.e. priors that ensure, up to the desired order of asymptotics, an agreement between Bayesian and frequentist procedures, can be extended to compute other Bayesian measures of evidence in the presence of sharp hypotheses with nuisance parameters (see, for example, [43]). Finally, we remark that frequentist properties of the MAP of the posterior based of the proposed median matching prior in comparison with the MAP of the posterior based of the predictive matching prior have been investigated in [28,41] for some of the examples discussed in this paper.

For inference on a full vector parameter $\theta$, with $d>1$ components, a direct extension of the rationale leading to (10) does not seem to be practicable due to lack of a manageable definition of multivariate median. Actually, in [28,44] it is shown how the method can be extended to a vector parameter of interest, in the presence of nuisance parameters, by simultaneously solving median bias corrected score equations for all parameter components. This leads to componentwise third-order median unbiasedness and parameterization equivariance. Moreover, the use of default priors involving all parameter components, also the nuisance, becomes necessary to regularize likelihoods in case of monotonicity.

As a final remark, we note that among the possible objective priors that ensures invariance of the posterior, we did not focus on the Jeffreys’ in the multidimensional case, since it often exhibits poor convergence properties. Conversely, the default matching priors considered in this paper are easily generalizable to the multidimensional case [44] preserving good convergence properties. Data and code used in the examples are available in https://github.com/elenabortolato/invariant_evalues.