Finite sample inference for empirical Bayesian methods

TL;DR

This paper introduces a universal finite-sample inference method for empirical Bayesian models using holdout likelihood ratios, enabling valid confidence sets and hypothesis tests in complex, high-dimensional settings.

Contribution

It proposes a new, general approach based on holdout likelihood ratios for constructing finite-sample valid confidence sets and hypothesis tests in empirical Bayesian inference.

Findings

Method provides finite-sample valid inference in complex models.

Numerical studies show effective and meaningful inferential statements.

Real data applications demonstrate practical utility.

Abstract

In recent years, empirical Bayesian (EB) inference has become an attractive approach for estimation in parametric models arising in a variety of real-life problems, especially in complex and high-dimensional scientific applications. However, compared to the relative abundance of available general methods for computing point estimators in the EB framework, the construction of confidence sets and hypothesis tests with good theoretical properties remains difficult and problem specific. Motivated by the universal inference framework of Wasserman et al. (2020), we propose a general and universal method, based on holdout likelihood ratios, and utilizing the hierarchical structure of the specified Bayesian model for constructing confidence sets and hypothesis tests that are finite sample valid. We illustrate our method through a range of numerical studies and real data applications, which demonstrate that the approach is able to generate useful and meaningful inferential statements in the relevant contexts.

Tables

Table 1. Table 1: Stein’s problem simulation results reported as average performances over 1000 1000 1000 replications.

$n$	${\psi }^2$	$\alpha$	Coverage of (9)	Coverage of (10)	Relative Width
10	$1^2$	0.05	0.948∗	1.000∗	1.979∗
		0.005	0.988∗	1.000∗	1.738∗
		0.0005	0.993∗	1.000∗	1.641∗
	$5^2$	0.05	0.943	1.000	1.902
		0.005	0.994	1.000	1.543
		0.0005	0.999	1.000	1.388
	${10}^2$	0.05	0.947	1.000	2.058
		0.005	0.994	1.000	1.633
		0.0005	0.999	1.000	1.455
100	$1^2$	0.05	0.937	0.999	2.068
		0.005	0.997	1.000	1.806
		0.0005	1.000	1.000	1.697
	$5^2$	0.05	0.949	1.000	1.912
		0.005	0.995	1.000	1.540
		0.0005	1.000	1.000	1.395
	${10}^2$	0.05	0.947	1.000	2.068
		0.005	0.995	1.000	1.635
		0.0005	0.999	1.000	1.455
1000	$1^2$	0.05	0.949	0.999	2.087
		0.005	0.991	1.000	1.815
		0.0005	1.000	1.000	1.705
	$5^2$	0.05	0.963	1.000	1.910
		0.005	0.997	1.000	1.544
		0.0005	1.000	1.000	1.399
	${10}^2$	0.05	0.942	1.000	2.066
		0.005	0.995	1.000	1.632
		0.0005	0.999	1.000	1.455

Table 2. Table 2: Experimental results for CIs constructed for Poisson–gamma count models. The Coverage and Length columns report the coverage proportion and average lengths in each scenario, as computed from 1000 replications.

$n$	$𝝍$	$\alpha$	Coverage	Length
10	$(2,2)$	0.05	0.998	3.632
		0.005	1.000	5.484
		0.0005	1.000	6.919
	$(2,5)$	0.05	0.999	2.976
		0.005	0.999	3.910
		0.0005	1.000	5.481
	$(5,2)$	0.05	0.997∗	5.468∗
		0.005	0.999∗	7.118∗
		0.0005	1.000∗	8.349∗
100	$(2,2)$	0.05	0.998	3.898
		0.005	0.999	5.277
		0.0005	1.000	6.883
	$(2,5)$	0.05	0.999	2.958
		0.005	1.000	3.914
		0.0005	1.000	5.374
	$(5,2)$	0.05	1.000	5.628
		0.005	1.000	7.124
		0.0005	1.000	8.529
1000	$(2,2)$	0.05	1.000	4.070
		0.005	1.000	5.424
		0.0005	1.000	6.344
	$(2,5)$	0.05	0.999	3.049
		0.005	1.000	3.960
		0.0005	1.000	5.479
	$(5,2)$	0.05	0.998	5.297
		0.005	1.000	7.205
		0.0005	1.000	8.714

Table 3. Table 3: Experimental results for testing the hypothesis H 0 subscript H 0 \text{H}_{0} : θ n − 1 ∗ = θ n ∗ superscript subscript 𝜃 𝑛 1 superscript subscript 𝜃 𝑛 \theta_{n-1}^{*}=\theta_{n}^{*} for Poisson–gamma count models. The Rejection Proportion columns report the average number of rejections, from 1000 tests, at levels of significance α ∈ { 0.05 , 0.005 , 0.0005 } 𝛼 0.05 0.005 0.0005 \alpha\in\left\{0.05,0.005,0.0005\right\} .

			Rejection Proportion at level $\alpha$
$n$	$𝝍$	$\Delta$	$0.05$	$0.005$	$0.0005$
10	$(2,2)$	0	0.000	0.000	0.000
		1	0.004	0.000	0.000
		5	0.280	0.193	0.128
		10	0.413	0.363	0.317
	$(2,5)$	0	0.000	0.000	0.000
		1	0.007	0.002	0.000
		5	0.143	0.096	0.064
		10	0.222	0.192	0.170
	$(5,2)$	0	0.001	0.000	0.000
		1	0.001	0.000	0.000
		5	0.177	0.107	0.052
		10	0.389	0.320	0.254
100	$(2,2)$	0	0.000	0.000	0.000
		1	0.014	0.003	0.000
		5	0.401	0.289	0.194
		10	0.562	0.489	0.427
	$(2,5)$	0	0.000	0.000	0.000
		1	0.015	0.000	0.000
		5	0.208	0.127	0.074
		10	0.296	0.235	0.179
	$(5,2)$	0	0.000	0.000	0.000
		1	0.004	0.000	0.000
		5	0.264	0.150	0.090
		10	0.500	0.425	0.344
1000	$(2,2)$	0	0.001	0.000	0.000
		1	0.021	0.001	0.000
		5	0.423	0.300	0.216
		10	0.576	0.513	0.450
	$(2,5)$	0	0.000	0.000	0.000
		1	0.012	0.000	0.000
		5	0.185	0.108	0.061
		10	0.321	0.254	0.197
	$(5,2)$	0	0.000	0.000	0.000
		1	0.003	0.001	0.000
		5	0.276	0.168	0.088
		10	0.507	0.428	0.354

Table 4. Table 4: Experimental results for testing the hypothesis H 0 subscript H 0 \text{H}_{0} : θ i ​ 1 ∗ = θ i ​ 2 ∗ superscript subscript 𝜃 𝑖 1 superscript subscript 𝜃 𝑖 2 \theta_{i1}^{*}=\theta_{i2}^{*} for Beta–binomial count series models. The Rejection proportion columns report the average number of rejections, from 100 test replicates, at levels of significance α ∈ { 0.05 , 0.005 , 0.0005 } 𝛼 0.05 0.005 0.0005 \alpha\in\left\{0.05,0.005,0.0005\right\} .

			Rejection proportion at level $\alpha$
$n$	$𝝍$	$\Delta$	$0.0005$	$0.005$	$0.05$
10	$(2,2)$	0	0.000	0.000	0.000
		0.2	0.000	0.004	0.039
		0.5	0.305	0.471	0.709
		0.9	0.980	1.000	1.000
	$(2,5)$	0	0.000	0.000	0.000
		0.2	0.000	0.001	0.025
		0.5	0.249	0.464	0.692
		0.9	0.995	1.000	1.000
	$(5,2)$	0	0.000	0.000	0.000
		0.2	0.000	0.006	0.052
		0.5	0.281	0.459	0.690
		0.9	0.993	0.993	1.000
100	$(2,2)$	0	0.000	0.000	0.000
		0.2	0.000	0.004	0.037
		0.5	0.272	0.459	0.700
		0.9	0.996	0.998	1.000
	$(2,5)$	0	0.000	0.000	0.000
		0.2	0.000	0.003	0.032
		0.5	0.267	0.459	0.693
		0.9	0.994	0.999	1.000
	$(5,2)$	0	0.000	0.000	0.000
		0.2	0.000	0.004	0.047
		0.5	0.269	0.459	0.697
		0.9	0.987	0.998	0.999
1000	$(2,2)$	0	0.000	0.000	0.000
		0.2	0.000	0.003	0.031
		0.5	0.280	0.476	0.707
		0.9	0.982	0.992	0.998
	$(2,5)$	0	0.000	0.000	0.000
		0.2	0.000	0.003	0.030
		0.5	0.264	0.459	0.693
		0.9	0.989	0.996	1.000
	$(5,2)$	0	0.000	0.000	0.000
		0.2	0.000	0.005	0.047
		0.5	0.279	0.474	0.706
		0.9	0.986	0.995	0.999

Table 5. Table 5: Comparison of differential methylation calling results between different methods: (i) FSEB (ii) Fisher tests with FDR-adjusted p 𝑝 p -values (FF) (iii) Fisher tests, unadjusted (F) (iv) score tests with FDR-adjusted p 𝑝 p -values (SF) and (v) score tests, unadjusted (S). The upper table gives the proportions of sites called to be differentially expressed under the tests of sizes α ∈ { 0.0005 , 0.05 } 𝛼 0.0005 0.05 \alpha\in\{0.0005,0.05\} . The lower table gives the proportion of overlaps between differential methylation calls from each pair of methods at a fixed level α ∈ { 0.0005 , 0.05 } 𝛼 0.0005 0.05 \alpha\in\{0.0005,0.05\} .

	Proportion of rejections at level
	$\alpha =0.0005$	$\alpha =0.05$
FSEB	0.0012	0.0154
FF	0.0003	0.0097
F	0.0098	0.1102
SF	0.1333	0.1528
S	0.1457	0.2926

Equations

$$f\left({\mathbf{x}}_i|{\mathbf{\Theta }}_i={\mathbf{\theta }}_i\right)=f\left({\mathbf{x}}_i|{\mathbf{\theta }}_i\right)\text{,}$$

$${\text{Pr}}_{{\mathbf{\theta }}_i^{\ast }}\left[{\mathbf{\theta }}_i^{\ast }\in {\mathcal{C}}_i^{\alpha }\left({\mathbf{D}}_n\right)\right]\ge 1-\alpha \text{,}$$

$$\underset{{\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }\in {\mathbb{T}}_{\mathbb{I},0}}{\sup }{\text{Pr}}_{{\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }}\left[P_{\mathbb{I}}\left({\mathbf{D}}_n\right)\le \alpha \right]\le \alpha \text{,}$$

$${\text{Pr}}_{{\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }}\left(\mathcal{A}\right)$$

$$\displaystyle{\color[rgb]{0,0,0}L_{\mathbb{I}}\left(\bm{\psi}\right)=\prod_{i\in\mathbb{I}}\int_{\mathbb{T}}f\left(\bm{X}_{i}|\bm{\theta}_{i}\right)\pi\left(\bm{\theta}_{i};\bm{\psi}\right)\text{d}\mathfrak{n}(\bm{\theta}_{i})}$$

$$l_{\mathbb{I}}\left({\mathbf{\vartheta }}_{\mathbb{I}}\right)=\prod _{i\in \mathbb{I}}f\left({\mathbf{X}}_i|{\mathbf{\theta }}_i\right)\text{,}$$

$$R_{\mathbb{I},n}\left({\mathbf{\vartheta }}_{\mathbb{I}}\right)=L_{\mathbb{I}}\left({\hat{\mathbf{\psi }}}_{\mathbb{I},n}\right)/l_{\mathbb{I}}\left({\mathbf{\vartheta }}_{\mathbb{I}}\right)\text{,}$$

$${\mathcal{C}}_i^{\alpha }\left({\mathbf{D}}_n\right)=\left\{\mathbf{\theta }\in \mathbb{T}:R_{i,n}\left(\mathbf{\theta }\right)\le 1/\alpha \right\}\text{.}$$

$${\text{E}}_{{\mathbf{\theta }}_{\mathbb{I}}^{\ast }}\left[R_{\mathbb{I},n}\left({\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }\right)\right]$$

$$\underset{\text{(i)}}{=}{\int }_{{\mathbb{X}}^{n-|\mathbb{I}|}}{\int }_{{\mathbb{X}}^{|\mathbb{I}|}}\frac{L_{\mathbb{I}}\left({\hat{\mathbf{\psi }}}_{\mathbb{I},n}\right)}{l_{\mathbb{I}}\left({\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }\right)}\prod _{i\in \mathbb{I}}f\left({\mathbf{x}}_i|{\mathbf{\theta }}_i^{\ast }\right)\text{d}\mathfrak{m}\left({\mathbf{d}}_{\mathbb{I}}\right)\prod _{j\notin \mathbb{I}}f\left({\mathbf{x}}_j|{\mathbf{\theta }}_j\right)\text{d}\mathfrak{m}\left({\bar{\mathbf{d}}}_{\mathbb{I}}\right)$$

$$\underset{\text{(ii)}}{=}{\int }_{{\mathbb{X}}^{n-|\mathbb{I}|}}{\int }_{{\mathbb{X}}^{|\mathbb{I}|}}{L}_{\mathbb{I}}\left({\hat{\mathbf{\psi }}}_{\mathbb{I},n}\right)\text{d}\mathfrak{m}\left({\mathbf{d}}_{\mathbb{I}}\right)\prod _{j\notin \mathbb{I}}f\left({\mathbf{x}}_j|{\mathbf{\theta }}_j\right)\text{d}\mathfrak{m}\left({\bar{\mathbf{d}}}_{\mathbb{I}}\right)$$

$$\underset{\text{(iii)}}{=}{\int }_{{\mathbb{X}}^{n-|\mathbb{I}|}}\prod _{j\notin \mathbb{I}}f\left({\mathbf{x}}_j|{\mathbf{\theta }}_j\right)\text{d}\mathfrak{m}\left({\bar{\mathbf{d}}}_{\mathbb{I}}\right)$$

$$\underset{\text{(iv)}}{=}1\text{.}$$

$${\mathrm{Pr}}_{{\mathbf{\theta }}_i^{\ast }}\left[{\mathbf{\theta }}_i^{\ast }\in {\mathcal{C}}_i^{\alpha }\left({\mathbf{D}}_n\right)\right]\ge 1-\alpha \text{.}$$

$${\mathrm{Pr}}_{{\mathbf{\theta }}_i^{\ast }}\left[R_{i,n}\left({\mathbf{\theta }}_i^{\ast }\right)\ge 1/\alpha \right]\le \alpha {\text{E}}_{{\mathbf{\theta }}_i^{\ast }}\left[R_{i,n}\left({\mathbf{\theta }}_i^{\ast }\right)\right]=\alpha \text{,}$$

$${\mathrm{Pr}}_{{\mathbf{\theta }}_i^{\ast }}\left[{\mathbf{\theta }}_i^{\ast }\in {\mathcal{C}}_i^{\alpha }\left({\mathbf{D}}_n\right)\right]={\mathrm{Pr}}_{{\mathbf{\theta }}_i^{\ast }}\left[R_{i,n}\left({\mathbf{\theta }}_i^{\ast }\right)\le 1/\alpha \right]\ge 1-\alpha$$

$${\tilde{\mathbf{\vartheta }}}_{\mathbb{I}}\in \left\{{\tilde{\mathbf{\vartheta }}}_{\mathbb{I}}\in {\mathbb{T}}_{\mathbb{I},0}:l_{\mathbb{I}}\left({\tilde{\mathbf{\vartheta }}}_{\mathbb{I}}\right)=\underset{{\mathbf{\vartheta }}_{\mathbb{I}}\in {\mathbb{T}}_{\mathbb{I},0}}{\sup }{l}_{\mathbb{I}}\left({\mathbf{\vartheta }}_{\mathbb{I}}\right)\right\}\text{.}$$

$$T_{\mathbb{I}}\left({\mathbf{D}}_n\right)=L_{\mathbb{I}}\left({\hat{\mathbf{\psi }}}_{\mathbb{I},n}\right)/l_{\mathbb{I}}\left({\tilde{\mathbf{\vartheta }}}_{\mathbb{I}}\right)\text{.}$$

$${\mathrm{Pr}}_{{\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }}\left[T_{\mathbb{I}}\left({\mathbf{D}}_n\right)\ge 1/\alpha \right]$$

$$=\alpha {\text{E}}_{{\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }}\left[\frac{L_{\mathbb{I}}\left({\hat{\mathbf{\psi }}}_{\mathbb{I},n}\right)}{l_{\mathbb{I}}\left({\tilde{\mathbf{\vartheta }}}_{\mathbb{I}}\right)}\right]\underset{\text{(i)}}{\le }\alpha {\text{E}}_{{\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }}\left[\frac{L_{\mathbb{I}}\left({\hat{\mathbf{\psi }}}_{\mathbb{I},n}\right)}{l_{\mathbb{I}}\left({\mathbf{\vartheta }}_{\mathbb{I}}^{\ast }\right)}\right]\underset{\text{(ii)}}{=}\alpha \text{,}$$

$$\text{N}\left({\bar{g}}_n{X}_n,{\bar{g}}_n+2\left[X_n{\left(1-{\bar{g}}_n\right)}^2\right]/\left[n-2\right]\right).$$

$$\text{E}\left[{\Theta }_n|{\mathbf{D}}_n\right]={\bar{g}}_n{X}_n\text{ and }\text{var}\left[{\Theta }_n|{\mathbf{D}}_n\right]={\bar{g}}_n+\frac{2X_n{\left(1-{\bar{g}}_n\right)}^2}{n-2}\text{,}$$

$${\bar{g}}_n{X}_n\pm {\zeta }_{1-\alpha /2}\sqrt{{\bar{g}}_n+\frac{2\left[X_n{\left(1-{\bar{g}}_n\right)}^2\right]}{n-2}}\text{,}$$

$$f\left(x_n|{\theta }_n\right)=\varphi \left(x_n;{\theta }_n,1\right)\text{ and }\pi \left({\theta }_n;\psi \right)=\varphi \left({\theta }_n;0,{\psi }^2\right)\text{,}$$

$$L_{\mathbb{I}}\left(\psi \right)={\int }_{\mathbb{R}}\varphi \left(X_n;\theta ,1\right)\varphi \left(\theta ;0,{\psi }^2\right)\text{d}\theta =\varphi \left(X_n;0,1+{\psi }^2\right)$$

$$R_{\mathbb{I},n}\left({\theta }_n\right)$$

$$=\varphi \left(X_n;0,1+{\hat{\psi }}_{-n}^2\right)/\varphi \left(X_n;{\theta }_n,1\right)\text{,}$$

$${\left(X_n-\theta \right)}^2\le 2\log \left(1/\alpha \right)+2\log \left(1+{\hat{\psi }}_{-n}^2\right)+\frac{X_n^2}{\left(1+{\hat{\psi }}_{-n}^2\right)}\text{,}$$

$$X_n\pm \sqrt{2\log \left(1/\alpha \right)+2\log \left(1+{\hat{\psi }}_{-n}^2\right)+\frac{X_n^2}{\left(1+{\hat{\psi }}_{-n}^2\right)}}\text{.}$$

$${\hat{\psi }}_{-n}^2=\max \left\{0,s_{-n}^2-1\right\}\text{,}$$

$$f\left(x_i,{\theta }_i;\mathbf{\psi }\right)=\frac{b^a}{\Gamma \left(a\right)}{\theta }_i^{a-1}\exp \left(b{\theta }_i\right)\frac{{\left({\theta }_i{w}_i\right)}^{x_i}\exp \left(-{\theta }_i{w}_i\right)}{x_i}\text{,}$$

$$f\left(x_i;\mathbf{\psi }\right)=\left(\genfrac{}{}{0ex}{}{x_i+a+1}{x_i}\right){\left(\frac{b}{w_i+b}\right)}^a{\left(\frac{w_i}{w_i+b}\right)}^{x_i}\text{,}$$

$${\hat{\theta }}_n=\frac{X_i+{\hat{a}}_n}{w_i+{\hat{b}}_n}\text{.}$$

Full text

1 Introduction

Let $\mathbf{D}_{n}=\left(\bm{X}_{i}\right)_{i\in\left[n\right]}$ be our data, presented as a sequence of $n\in\mathbb{N}=\left\{1,2,\dots\right\}$ random variables $\bm{X}_{i}\in\mathbb{X}$ ($i\in\left[n\right]=\left\{1,\dots,n\right\}$). For each $i\in\left[n\right]$, let $\bm{\Theta}_{i}\in\mathbb{T}$ be a random variable with probability density function (PDF) $\pi\left(\bm{\theta}_{i};\bm{\psi}\right)$, where $\bm{\psi}\in\mathbb{P}$ is a hyperparameter. Furthermore, suppose that $\left[\bm{X}_{i}|\bm{\Theta}_{i}=\bm{\theta}_{i}\right]$ arises from a family of data generating processes (DGPs) with conditional PDFs

[TABLE]

and that the sequence $\left(\left(\bm{X}_{i},\bm{\Theta}_{i}\right)\right)_{i\in\left[n\right]}$ is independent.

Suppose that $\left(\bm{\Theta}_{i}\right)_{i\in\left[n\right]}$ is realized at $\bm{\vartheta}_{n}^{*}=\left(\bm{\theta}_{i}^{*}\right)_{i\in\left[n\right]}$, where each realization $\bm{\theta}_{i}^{*}$ ($i\in\left[n\right]$) is unknown, and where $\bm{\psi}$ is also unknown. Let $\mathbb{I}\subset\left[n\right]$, and write $\bm{\vartheta}_{\mathbb{I}}^{*}=\left(\bm{\theta}_{i}^{*}\right)_{i\in\mathbb{I}}$. When $\mathbb{I}=\left\{i\right\}$, we shall use the shorthand $\mathbb{I}=i$, where it causes no confusion.

Under this setup, for significance level $\alpha\in\left(0,1\right)$, we wish to draw inference regarding the realized sequence $\bm{\vartheta}_{n}^{*}$ by way of constructing $100\left(1-\alpha\right)\%$ confidence sets $\mathcal{C}_{i}^{\alpha}\left(\mathbf{D}_{n}\right)$ that satisfy:

[TABLE]

and $p$-values $P_{\mathbb{I}}\left(\mathbf{D}_{n}\right)$ for testing null hypotheses $\text{H}_{0}:\bm{\vartheta}_{\mathbb{I}}^{*}\in\mathbb{T}_{\mathbb{I},0}\subset\mathbb{T}^{\left|\mathbb{I}\right|}$ that satisfy:

[TABLE]

where $\text{Pr}_{\bm{\theta}_{i}^{*}}$ and $\text{Pr}_{\bm{\vartheta}_{\mathbb{I}}^{*}}$ denote probability measures consistent with the PDF $f\left(\bm{x}_{i}|\bm{\theta}_{i}^{*}\right)$, for each $i\in\left[n\right]$, and for all $i\in\mathbb{I}$, respectively. That is, for a measurable set $\mathcal{A}\subset\mathbb{X}^{n}$, and assuming absolute continuity of $\text{Pr}_{\bm{\vartheta}_{\mathbb{I}}^{*}}$ with respect to some measure $\mathfrak{m}$ (typically the Lebesgue or counting measure), we can write

[TABLE]

where $\bm{\theta}_{j}$ is an arbitrary element of $\mathbb{T}$, for each $j\notin\mathbb{I}$.

The setup above falls within the framework of empirical Bayesian (EB) inference, as exposited in the volumes of Maritz and Lwin (1989), Ahmed and Reid (2001), Serdobolskii (2008), Efron (2010), and Bickel (2020). Over the years, there has been a sustained interest in the construction and computation of EB point estimators for $\bm{\vartheta}_{n}^{*}$, in various contexts, with many convenient and general computational tools now made available, for instance, via the software of Johnstone and Silverman (2005), Leng et al. (2013), Koenker and Gu (2017), and Narasimhan and Efron (2020). Unfortunately, the probabilistic properties of $\bm{\vartheta}_{n}^{*}$ tend to be difficult to characterize, making the construction of confidence sets and hypothesis tests with good theoretical properties relatively less routine than the construction of point estimators. When restricted to certain classes of models, such constructions are nevertheless possible, as exemplified by the works of Casella and Hwang (1983), Morris (1983a), Laird and Louis (1987), Datta et al. (2002), Tai and Speed (2006), Hwang et al. (2009), Hwang and Zhao (2013), and Yoshimori and Lahiri (2014), among others.

In this work, we adapt the universal inference framework of Wasserman et al. (2020) to produce valid confidence sets and $p$-values with properties (1) and (2), respectively, for arbitrary estimators of $\bm{\vartheta}_{n}^{*}$. As with the constructions of Wasserman et al. (2020), the produced inferential methods are all valid for finite sample size $n$ and require no assumptions beyond correctness of model specification. The confidence sets and $p$-values arise by construction of holdout likelihood ratios that can be demonstrated to have the $e$-value property, as described in Vovk and Wang (2021) (see also the $s$-values of Grunwald et al., 2020 and the betting values of Shafer, 2021). Here, we are able to take into account the hierarchical structure of the Bayesian specified model by using the fact that parameterized $e$-values are closed when averaged with respect to an appropriate probability measure (cf. Vovk, 2007 and Kaufmann and Koolen, 2018). Due to the finite sample correctness of our constructions, we shall refer to our methods as finite sample EB (FSEB) techniques.

Along with our methodological developments, we also demonstrate the application of our FSEB techniques in numerical studies and real data applications. These applications include the use of FSEB methods for constructing confidence intervals (CIs) for the classic mean estimator of Stein (1956), and testing and CI construction in Poisson–gamma models and Beta–binomial models, as per Koenker and Gu (2017) and Hardcastle and Kelly (2013), respectively. Real data applications are demonstrated via the analysis of insurance data from Haastrup (2000) and differential methylation data from Cruickshanks et al. (2013). In these real and synthetic applications, we show that FSEB methods, satisfying conditions (1) and (2), are able to generate useful and meaningful inferential statements.

We proceed as follows. In Section 2, we introduce the confidence set and $p$-value constructions for drawing inference regarding EB models. In Section 3, numerical studies of simulated data are used to demonstrate the applicability and effectiveness of FSEB constructions. In Section 4, FSEB methods are applied to real data to further show the practicality of the techniques. Lastly, in Section 5, we provide discussions and conclusions regarding our results.

2 Confidence sets and hypothesis tests

We retain the notation and setup from Section 1. For each subset $\mathbb{I}\subset\left[n\right]$, let us write $\mathbf{D}_{\mathbb{I}}=\left(\bm{X}_{i}\right)_{i\in\mathbb{I}}$ and $\overline{\mathbf{D}}_{\mathbb{I}}=\left(\bm{X}_{i}\right)_{i\in\left[n\right]\backslash\mathbb{I}}$.

Suppose that we have available some estimator of $\bm{\psi}$ that only depends on $\overline{\mathbf{D}}_{\mathbb{I}}$ (and not $\mathbf{D}_{\mathbb{I}}$), which we shall denote by $\hat{\bm{\psi}}_{\mathbb{I},n}$. Furthermore, for fixed $\bm{\psi}$, write the integrated and unintegrated likelihood of the data $\mathbf{D}_{\mathbb{I}}$, as

[TABLE]

and

[TABLE]

respectively, where $\bm{\vartheta}_{\mathbb{I}}=\left(\bm{\theta}_{i}\right)_{i\in\mathbb{I}}$ (here, $\bm{\vartheta}_{\left\{i\right\}}=\bm{\theta}_{i}$). We note that in (4), we have assumed that $\pi(\cdot;\bm{\psi})$ is a density function with respect to some measure on $\mathbb{T}$, $\mathfrak{n}$.

Define the ratio statistic:

[TABLE]

and consider sets of the form

[TABLE]

The following Lemma is an adaptation of the main idea of Wasserman et al. (2020) for the context of empierical Bayes estimators, and allows us to show that $\mathcal{C}_{i}^{\alpha}\left(\mathbf{D}_{n}\right)$ satisfies property (1).

Lemma 1**.**

For each $\mathbb{I}\subset\left[n\right]$ and fixed sequence $\bm{\vartheta}_{n}^{*}\in\mathbb{T}^{n}$, $\mathrm{E}_{\bm{\vartheta}_{\mathbb{I}}^{*}}\left[R_{\mathbb{I},n}\left(\bm{\vartheta}_{\mathbb{I}}^{*}\right)\right]=1$.

Proof.

Let $\mathbf{d}_{\mathbb{I}}$ and $\bar{\mathbf{d}}_{\mathbb{I}}$ be realizations of $\mathbf{D}_{\mathbb{I}}$ and $\overline{\mathbf{D}}_{\mathbb{I}}$, respectively. Then, using (3), write

[TABLE]

Here, (i) is true by definition of (6), (ii) is true by definition of (5), (iii) is true by the fact that (4) is a probability density function on $\mathbb{X}^{\left|\mathbb{I}\right|}$, with respect to $\mathfrak{m}$, and (iv) is true by the fact that $\prod_{j\notin\mathbb{I}}f\left(\bm{x}_{j}|\bm{\theta}_{j}\right)$ is a probability density function on $\mathbb{X}^{n-\left|\mathbb{I}\right|}$, with respect to $\mathfrak{m}$.

∎

Proposition 1**.**

For each $i\in\left[n\right]$, $\mathcal{C}_{i}^{\alpha}\left(\mathbf{D}_{n}\right)$ is a $100\left(1-\alpha\right)\%$ confidence set, in the sense that

[TABLE]

Proof.

For each $i$, Markov’s inequality states that

[TABLE]

which implies that

[TABLE]

by Lemma 1. ∎

Next, we consider the testing of null hypotheses $\text{H}_{0}\text{: }\bm{\vartheta}_{\mathbb{I}}^{*}\in\mathbb{T}_{\mathbb{I},0}$ against an arbitrary alternative $\text{H}_{1}\text{: }\bm{\vartheta}_{\mathbb{I}}^{*}\in\mathbb{T}_{\mathbb{I},1}\subseteq\mathbb{T}^{\left|\mathbb{I}\right|}$. To this end, we define the maximum unintegrated likelihood estimator of $\bm{\vartheta}_{\mathbb{I}}^{*}$, under $\text{H}_{0}$ as

[TABLE]

Using (7), and again letting $\hat{\bm{\psi}}_{\mathbb{I},n}$ be an arbitrary estimator of $\bm{\psi}$, depending only on $\overline{\mathbf{D}}_{\mathbb{I}}$, we define the ratio test statistic

[TABLE]

The following result establishes the fact that the $p$-value $P_{\mathbb{I}}\left(\mathbf{D}_{n}\right)=1/T_{\mathbb{I}}\left(\mathbf{D}_{n}\right)$ has the correct size, under $\text{H}_{0}$.

Proposition 2**.**

For any $\alpha\in\left(0,1\right)$ and $\bm{\vartheta}_{\mathbb{I}}^{*}\in\mathbb{T}_{\mathbb{I},0}$, $\Pr_{\bm{\vartheta}_{\mathbb{I}}^{*}}\left[P_{\mathbb{I}}\left(\mathbf{D}_{n}\right)\leq\alpha\right]\leq\alpha$.

Proof.

Assume that $\bm{\vartheta}_{\mathbb{I}}^{*}\in\mathbb{T}_{\mathbb{I},0}$. By Markov’s inequality, we have

[TABLE]

where the (i) is true due to the fact that $l_{\mathbb{I}}\left(\tilde{\bm{\vartheta}}_{\mathbb{I}}\right)\geq l_{\mathbb{I}}\left(\bm{\vartheta}_{\mathbb{I}}^{*}\right)$, by the definition of (7), and the (ii) is true due to Lemma 1. ∎

We note that Propositions 1 and 2 are empirical Bayes analogues of Theorems 1 and 2 from Wasserman et al. (2020), which provide guarantees for universal inference confidence set and hypothesis test constructions, respectively. Furthermore, the use of Lemma 1 in the proofs also imply that the CIs constructed via Proposition 1 are $e$-CIs, as defined by Xu et al. (2022), and the $p$-values obtained via Proposition 2 can be said to be $e$-value calibrated, as per the definitions of Wang and Ramdas (2022).

3 FSEB examples and some numerical results

To demonstrate the usefulness of the FSEB results from Section 2, we shall present a number of synthetic and real world applications of the confidence and testing constructions. All of the computation is conducted in the R programming environment (R Core Team, 2020) and replicable scripts are made available at https://github.com/hiendn/Universal_EB. Where unspecified, numerical optimization is conducted using the optim() or optimize() functions in the case of multivariate and univariate optimization, respectively.

3.1 Stein’s problem

We begin by studying the estimation of normal means, as originally considered in Stein (1956). Here, we largely follow the exposition of Efron (2010, Ch. 1) and note that the estimator falls within the shrinkage paradigm exposited in Serdobolskii (2008). We consider this setting due to its simplicity and the availability of a simple EB-based method to compare our methodology against.

Let $\left(\left(X_{i},\Theta_{i}\right)\right)_{i\in\left[n\right]}$ be IID and for each $i\in\left[n\right]$, $\Theta_{i}\sim\text{N}\left(0,\psi^{2}\right)$ ($\psi^{2}>0$) and $\left[X_{i}|\Theta_{i}=\theta_{i}\right]\sim\text{N}\left(\theta_{i},1\right)$, where $\text{N}\left(\mu,\sigma^{2}\right)$ is the normal law with mean $\mu\in\mathbb{R}$ and variance $\sigma^{2}>0$. We assume that $\psi^{2}$ is unknown and that we observe data $\mathbf{D}_{n}$ and wish to construct CIs for the realizations $\theta_{n}^{*}$, which characterize the DGP of the observations $X_{n}$.

Following Efron (2010, Sec. 1.5), when $\psi^{2}$ is known, the posterior distribution of $\left[\Theta_{n}|X_{n}=x_{n}\right]$ is $\text{N}\left(g\left(\psi^{2}\right)x_{n},g\left(\psi^{2}\right)\right)$, where $g\left(\psi^{2}\right)=\psi^{2}/\left(1+\psi^{2}\right)$. Using the data $\mathbf{D}_{n}$, we have the fact that $\sum_{i=1}^{n-1}X_{i}^{2}\sim\left(\psi^{2}+1\right)\chi_{n-1}^{2}$, where $\chi_{\nu}^{2}$ is the chi-squared distribution with $\nu$ degrees of freedom. This implies a method-of-moment estimator for $g$ of the form: $\bar{g}_{n}=1-\left(n-2\right)/\sum_{i=1}^{n}X_{i}^{2}$, in the case of unknown $\psi^{2}$.

We can simply approximate the distribution of $\left[\Theta_{n}|\mathbf{D}_{n}\right]$ as $\text{N}\left(\bar{g}_{n}X_{n},\bar{g}_{n}\right)$, although this approximation ignores the variability of $\bar{g}_{n}$. As noted by Efron (2010, Sec. 1.5), via a hierarchical Bayesian interpretation using an objective Bayesian prior, we may instead deduce the more accurate approximate distribution:

[TABLE]

Specifically, Efron (2010) considers the hyperparameter $\psi^{2}$ as being a random variable, say $\Psi^{2}$, and places a so-called objective (or non-informative) prior on $\Psi^{2}$. In particular, the improper prior assumption that $\Psi^{2}+1\sim\text{Uniform}\left(0,\infty\right)$ is made. Then, it follows from careful derivation that

[TABLE]

and thus we obtain (8) via a normal approximation for the distribution of $\left[\Theta_{n}|\mathbf{D}_{n}\right]$ (cf. Morris 1983b, Sec. 4).

The approximation then provides $100\left(1-\alpha\right)\%$ posterior credible intervals for $\Theta_{n}$ of the form

[TABLE]

where $\zeta_{1-\alpha/2}$ is the $\left(1-\alpha/2\right)$ quantile of the standard normal distribution. This posterior result can then be taken as an approximate $100\left(1-\alpha\right)\%$ confidence interval for $\theta_{n}^{*}$.

Now, we wish to apply the FSEB results from Section 2. Here, $\mathbb{I}=\left\{n\right\}$, and from the setup of the problem, we have

[TABLE]

where $\phi\left(x;\mu,\sigma^{2}\right)$ is the normal PDF with mean $\mu$ and variance $\sigma^{2}$. Thus,

[TABLE]

and $l_{\mathbb{I}}\left(\theta_{n}\right)=\phi\left(x_{n};\theta_{n},1\right)$, which yields a ratio statistic of the form

[TABLE]

when combined with an appropriate estimator $\hat{\psi}_{-n}^{2}$ for $\psi^{2}$, using only $\bar{\mathbf{D}}_{\mathbb{I},n}=\mathbf{D}_{n-1}$. We can obtain the region $\mathcal{C}_{\mathbb{I}}^{\alpha}\left(\mathbf{D}_{n}\right)$ by solving $R_{\mathbb{I},n}\left(\theta_{n}\right)\leq 1/\alpha$ to obtain:

[TABLE]

which, by Proposition 1, yields the $100\left(1-\alpha\right)\%$ CI for $\theta_{n}^{*}$:

[TABLE]

We shall consider implementations of the CI of form (10) using the estimator

[TABLE]

where $s_{-n}^{2}$ is the sample variance of the $\bar{\mathbf{D}}_{\mathbb{I},n}$, and $s_{-n}^{2}-1$ is the method of moment estimator of $\psi^{2}$. The maximum operator stops the estimator from becoming negative and causes no problems in the computation of (10).

We now compare the performances of the CIs of forms (9) and (10). To do so, we shall consider data sets of sizes $n\in\left\{10,100,1000\right\}$, $\psi^{2}\in\left\{1^{2},5^{2},10^{2}\right\}$, and $\alpha\in\left\{0.05,0.005,0.0005\right\}$. For each triplet $\left(n,\psi^{2},\alpha\right)$, we repeat the computation of (9) and (10) $1000$ times and record the coverage probability and average relative widths of the intervals (computed as the width of (10) divided by that of (9)). The results of our experiment are presented in Table 1.

From Table 1, we observe that the CIs of form (9) tended to produce intervals with the desired levels of coverage, whereas the FSEB CIs of form (10) tended to be conservative and contained the parameter of interest in almost all replications. The price that is paid for this conservativeness is obvious when viewing the relative widths, which implies that for $95\%$ CIs, the EB CIs of form (10) are twice as wide, on average, when compared to the CIs of form (9). However, the relative widths decrease as $\alpha$ gets smaller, implying that the intervals perform relatively similarly when a high level of confidence is required. We further observe that $n$ and $\psi^{2}$ had little effect on the performances of the intervals except in the case when $n=10$ and $\psi^{2}=1$, whereupon it was possible for the intervals of form (9) to not be computable in some cases.

From these results we can make a number of conclusions. Firstly, if one is willing to make the necessary hierarchical and objective Bayesian assumptions, as stated in Efron (2010, Sec. 1.5), then the intervals of form (9) provide very good performance. However, without those assumptions, we can still obtain reasonable CIs that have correct coverage via the FSEB methods from Section 2. Furthermore, these intervals become more efficient compared to (9) when higher levels of confidence are desired. Lastly, when $n$ is small and $\psi^{2}$ is also small, the intervals of form (9) can become uncomputable and thus one may consider the use of (10) as an alternative.

3.2 Poisson–gamma count model

The following example is taken from Koenker and Gu (2017) and was originally studied in Norberg (1989) and then subsequently in Haastrup (2000). In this example, we firstly consider IID parameters $\left(\Theta_{i}\right)_{i\in\left[n\right]}$ generated with gamma DGP: $\Theta_{i}\sim\text{Gamma}\left(a,b\right)$, for each $i\in\left[n\right]$, where $a>0$ and $b>0$ are the shape and rate hyperparameters, respectively, which we put into $\bm{\psi}$. Then, for each $i$, we suppose that the data $\mathbf{D}_{n}=\left(X_{i}\right)_{i\in\left[n\right]}$, depending on the covariate sequence $\mathbf{w}_{n}=\left(w_{i}\right)_{i\in\left[n\right]}$, has the Poisson DGP: $\left[X_{i}|\Theta_{i}=\theta_{i}\right]\sim\text{Poisson}\left(\theta_{i}w_{i}\right)$, where $w_{i}>0$. We again wish to use the data $\mathbf{D}_{n}$ to estimate the realization of $\Theta_{n}$: $\theta_{n}^{*}$, which characterizes the DGP of $X_{n}$.

Under the specification above, for each $i$, we have the fact that $\left(X_{i},\Theta_{i}\right)$ has the joint PDF:

[TABLE]

which we can marginalize to obtain

[TABLE]

and which can be seen as a Poisson–gamma mixture model. We can then construct the likelihood of $\mathbf{D}_{n}$ using expression (12), from which we may compute maximum likelihood estimates $\hat{\bm{\psi}}_{n}=\left(\hat{a}_{n},\hat{b}_{n}\right)$ of $\bm{\psi}$. Upon noting that (11) implies the conditional expectation $\text{E}\left[\Theta_{i}|X_{i}=x_{i}\right]=\left(x_{i}+a\right)/\left(w_{i}+b\right)$, we obtain the estimator for $\theta_{n}^{*}$:

[TABLE]

3.2.1 Confidence intervals

We again wish to apply the general result from Section 2 to construct CIs. Firstly, we have $\mathbb{I}=\left\{n\right\}$ and

[TABLE]

As per (12), we can write

[TABLE]

Then, since $l_{\mathbb{I}}\left(\theta_{n}\right)=f\left(X_{n}|\theta_{n}\right)$, we have

[TABLE]

when combined with an estimator $\hat{\bm{\psi}}_{-n}=\left(\hat{a}_{-n},\hat{b}_{-n}\right)$ of $\bm{\psi}$, using only $\bar{\mathbf{D}}_{\mathbb{I},n}=\mathbf{D}_{n-1}$.

For any $\alpha\in\left(0,1\right)$, we then obtain a $100\left(1-\alpha\right)\%$ CI for $\theta_{n}$ by solving $R_{\mathbb{I},n}\left(\theta_{n}\right)\leq 1/\alpha$, which can be done numerically. We shall use the MLE of $\bm{\psi}$, computed with the data $\bar{\mathbf{D}}_{\mathbb{I},n}$ and marginal PDF (12), as the estimator $\hat{\bm{\psi}}_{-n}$.

To demonstrate the performance of the CI construction, above, we conduct the following numerical experiment. We generate data sets consisting of $n\in\left\{10,100,1000\right\}$ observations characterized by hyperparameters $\bm{\psi}=\left(a,b\right)=\left\{\left(2,2\right),\left(2,5\right),\left(5,2\right)\right\}$, and we compute intervals using significance levels $\alpha\in\left\{0.05,0.005,0.0005\right\}$. Here, we shall generate $\mathbf{w}_{n}$ IID uniformly between 0 and 10. For each triplet $\left(n,\bm{\psi},\alpha\right)$, we repeat the construction of our CIs 1000 times and record the coverage probability and average width for each case. The results of the experiment are reported in Table 2.

From Table 2, we observe that the empirical coverage of the CIs are higher than the nominal value and are thus behaving as per the conclusions of Proposition 1. As expected, we also find that increasing the nominal confidence level also increases the coverage proportion, but at a cost of increasing the lengths of the CIs. From the usual asymptotic theory of maximum likelihood estimators, we anticipate that increasing $n$ will decrease the variance of the estimator $\hat{\bm{\psi}}_{-n}$. However, as in Section 3.1, this does not appear to have any observable effect on either the coverage proportion nor lengths of the CIs.

3.2.2 Hypothesis tests

Next, we consider testing the null hypothesis $\text{H}_{0}$: $\theta_{n-1}^{*}=\theta_{n}^{*}$. To this end, we use the hypothesis testing framework from Section 2. That is, we let $\mathbb{I}=\left\{n-1,n\right\}$ and estimate $\bm{\psi}$ via the maximum likelihood estimator $\hat{\bm{\psi}}_{\mathbb{I},n}=\left(a_{\mathbb{I},n},b_{\mathbb{I},n}\right)$, computed from the data $\bar{\mathbf{D}}_{\mathbb{I},n}=\mathbf{D}_{n-2}$.

We can write

[TABLE]

[TABLE]

and $\bm{\vartheta}_{\mathbb{I}}^{*}=\left(\theta_{n-1}^{*},\theta_{n}^{*}\right)$. We are also required to compute the maximum likelihood estimator of $\bm{\vartheta}_{\mathbb{I}}^{*}$, under $\text{H}_{0}$, as per (7), which can be written as

[TABLE]

Using the components above, we define the test statistic $T_{\mathbb{I}}\left(\mathbf{D}_{n}\right)=L_{\mathbb{I}}\left(\hat{\bm{\psi}}_{\mathbb{I},n}\right)/l_{\mathbb{I}}\left(\tilde{\bm{\vartheta}}_{\mathbb{I}}\right)$, from which we can derive the $p$-value $P_{\mathbb{I}}\left(\mathbf{D}_{n}\right)=1/T_{\mathbb{I}}\left(\mathbf{D}_{n}\right)$ for testing $\text{H}_{0}$.

To demonstrate the application of this test, we conduct another numerical experiment. As in Section 3.2.1, we generate data sets of sizes $n\in\left\{10,100,1000\right\}$, where the data $\mathbf{D}_{n-1}$ are generated with parameters $\left(\Theta_{i}\right)_{i\in\left[n-1\right]}$ arising from gamma distributions with hyperparameters $\bm{\psi}=\left(a,b\right)=\left\{\left(2,2\right),\left(2,5\right),\left(5,2\right)\right\}$. The final observation $X_{n}$, making up $\mathbf{D}_{n}$, is then generated with parameter $\Theta_{n}=\Theta_{n-1}+\Delta$, where $\Delta\in\left\{0,1,5,10\right\}$. As before, we generate the covariate sequence $\mathbf{w}_{n}$ IID uniformly between 0 and 10. For each triplet $\left(n,\bm{\psi},\Delta\right)$, we test $\text{H}_{0}$: $\theta_{n-1}^{*}=\theta_{n}^{*}$ 1000 times and record the average number of rejections under at the levels of significance $\alpha\in\left\{0.05,0.005,0.0005\right\}$. The results are then reported in Table 3.

The results for the $\Delta=0$ cases in Table 3 show that the tests reject true null hypotheses at below the nominal sizes $\alpha$, in accordance with Proposition 2. For each combination of $n$ and $\bm{\psi}$, as $\Delta$ increases, the proportion of rejections increase, demonstrating that the tests become more powerful when detecting larger differences between $\theta_{n-1}^{*}$ and $\theta_{n}^{*}$, as expected. There also appears to be an increase in power due to larger sample sizes. This is an interesting outcome, since we can only be sure that sample size affects the variability of the estimator $\bm{\psi}_{\mathbb{I},n}$. Overall, we can be confident that the tests are behaving as required, albeit they may be somewhat underpowered as they are not achieving the nominal sizes.

3.3 Beta–binomial data series

Data from genome-level biological studies, using modern high-throughput sequencing technologies (Krueger et al., 2012), often take the form of a series of counts, which may be modelled through sets of non-identical (possibly correlated) binomial distributions, with beta priors, in a Bayesian framework. The question of interest may vary, for example, from assessing the range of likely values for the binomial parameter in a particular region of the data, to comparing whether two sections of one or more data series are generated from identical distributions. For purposes of demonstrating the performance of the FSEB method in these scenario, we will make the simplifying assumption that all data points are independently distributed, within, as well as across, any of $G$ data series that may be observed.

3.3.1 Confidence Sets

First, let us assume that we only have a single series, i.e. $G=1$. Then, we can assume $X_{i}\sim\mbox{Bin}(m_{i},\theta_{i})$, and propose a common prior distribution for $\Theta_{i}$ $(i=1,\ldots,n)$: Beta($\gamma,\beta$). Using the techniques described in Section 2, we can find confidence sets for $\theta^{*}_{i}$, $(i=1,\ldots,n)$. For each $i$, we define, as previously, a subset $\mathbb{I}=\{i\}$, so that $\mathbf{D}_{\mathbb{I}}=X_{i}$ and $\overline{\mathbf{D}}_{\mathbb{I}}=\left(X_{i}\right)_{i\in\left[n\right]\backslash\{i\}}$. We then have,

[TABLE]

where

[TABLE]

and

[TABLE]

which gives the ratio

[TABLE]

Here, $\hat{\gamma}_{-n}$ and $\hat{\beta}_{-n}$ are the empirical Bayes estimates of $\gamma$ and $\beta$, given by

[TABLE]

and

[TABLE]

where

[TABLE]

$\bar{m}=\frac{1}{n-1}\sum_{j\in\left[n\right]\backslash i}m_{j}$, and $\hat{V}_{x}=\frac{1}{n-1}\sum_{j\in\left[n\right]\backslash i}(\frac{x_{j}}{m_{j}}-\hat{\mu}_{\text{EB}})^{2}$. Further, $B\left(a,b\right)=\int_{0}^{1}t^{a-1}\left(1-t\right)^{b-1}\text{d}t$ is the Beta function, taking inputs $a>0$ and $b>0$.

We simulated data from the binomial model under two cases: (a) setting beta hyperparameters $(\alpha,\beta)=(10,10)$, and hierarchically simulating $\theta_{i}^{*}$, $i\in[n]$, and then $x_{i}$ from a binomial distribution; and (b) setting a range of $\theta_{i}^{*}$ ($i\in[n]$) values equidistantly spanning the interval $(0.1,0.9)$ for $n=10,100$. Here, $m_{i}$ ($i\in[n]$) were given integer values uniformly generated in the range $[15,40]$. In all cases, it was seen that the CIs had perfect coverage, always containing the true value of $\theta_{i}^{*}$. An example of the $n=10$ case is shown in Figure 1.

3.3.2 Hypothesis testing

Aiming to detect genomic regions that may have differing characteristics between two series, a pertinent question of interest may be considered by testing the hypotheses: $H_{0}$: $\theta_{i1}^{*}=\theta_{i2}^{*}$ vs. $H_{1}$: $\theta_{i1}^{*}\neq\theta_{i2}^{*}$, for every $i\in[n]$ (with $G=2$ series). Then, $\mathbf{D}_{n}=\left(\bm{X}_{i}\right)_{i\in\left[n\right]}$, where $\bm{X}_{i}=(X_{i1},X_{i2})$. From Section 2, the ratio test statistic takes the form

[TABLE]

where $\hat{\gamma}_{\mathbb{I},n}$ and $\hat{\beta}_{\mathbb{I},n}$ are EB estimators of $\gamma$ and $\beta$, depending only on $\bar{\mathbf{D}}_{\mathbb{I},n}=\mathbf{D}_{n}\backslash\{X_{i1},X_{i2}\}$. With $\tilde{{\vartheta}}_{\mathbb{I}}=\frac{x_{i1}+x_{i2}}{m_{i1}+m_{i2}}=\tilde{\theta}_{i}$, write $l_{\mathbb{I}}\left(\tilde{{\vartheta}}_{\mathbb{I}}\right)=f(x_{i1},x_{i2}|\tilde{\theta}_{i})$, and

[TABLE]

which gives

[TABLE]

where $\hat{\gamma}_{\mathbb{I},n}$ and $\hat{\beta}_{\mathbb{I},n}$ are calculated in a similar fashion to Section 3.3.1 except that data from both sequences should be used to estimate $\hat{\mu}_{\text{EB}}$ and $\hat{\phi}_{\text{EB}}$, in the sense that

[TABLE]

where

[TABLE]

In our first simulation, we assessed the performance of the test statistic in terms of the Type I error. Assuming a window size of $n=20$, realized data $(x_{i1},x_{i2})$ ($i\in[n]$), were simulated from independent binomial distributions with $\theta_{i1}^{*}=\theta_{i2}^{*}=\theta_{i}^{*}$ ($i=1,\ldots,n$), with $\theta_{i}^{*}$ ranging between $0.1$ and $0.9$, and $m_{i1},m_{i2}\in\mathbb{N}$ uniformly and independently sampled from the range $[15,40]$. The first panel of Figure 2 shows the calculated test statistic values $T_{\mathbb{I}}\left(\mathbf{D}_{n}\right)$ for the $20$ genomic indices on the logarithmic scale, over 100 independently replicated datasets, with horizontal lines displaying values $\log(1/\alpha)$, for significance levels $\alpha\in\{0.01,0.02,0.05\}$. No points were observed above the line corresponding to $\alpha=0.01$, indicating that the Type I error of the test statistic does not exceed the nominal level. Next, we assessed the power of the test statistic at three levels of significance ($\alpha\in\{0.01,0.02,0.05\}$) and differing effect sizes. For each $i$ ($i\in[n]$), $\theta_{i1}^{*}$ was set to be a value between $0.05$ and $0.95$, and $\theta_{i2}^{*}=\theta_{i1}^{*}+\Delta$, where $0.1<\Delta<0.9$ (with $\theta_{i2}^{*}<1$). A sample of $20$ replicates were simulated under each possible set of values of $(\theta_{1}^{*},\theta_{2}^{*})$. The second panel of Figure 2 shows that the power functions increased rapidly to 1 as the difference $\Delta$ was increased.

In our next numerical experiment, we generated data sets of sizes $n\in\{10,100,1000\}$, where realized observations $x_{i1}$, and $x_{i2}$ are simulated from independent binomial distributions with parameters $\theta_{i1}^{*}$ and $\theta_{i2}^{*}$, respectively ($i\in[n]$). For each $i$, $\theta_{i1^{*}}$ was generated from a beta distribution, in turn, with hyperparameters $\bm{\psi}=(\gamma,\beta)\in\{(2,2),(2,5),(5,2)\}$; and $\theta_{i2}^{*}=\theta_{i1}^{*}+\Delta$, where $\Delta\in\{0,0.2,0.5,0.9\}$. We generated 100 instances of data under each setting and assessed the power of the FSEB test statistic through the number of rejections at levels $\alpha\in\{0.0005,0.005,0.05\}$. The results are shown in Table 4.

Similarly to the Poisson–gamma example, it can be seen that the tests reject true null hypotheses at below the nominal sizes $\alpha$, in each case. For each combination of $n$ and $\bm{\psi}$, as $\Delta$ increases, the rejection rate increases, making the tests more powerful as expected, when detecting larger differences between $\theta_{i1}^{*}$ and $\theta_{i2}^{*}$, frequently reaching a power of 1 even when the difference was not maximal. There did not appear to be a clear increase in power with the sample size, within the settings considered. Overall, we may conclude, as previously, that the tests are behaving as expected, although both this example and the Poisson–gamma case show that the tests may be underpowered as they do not achieve the nominal size for any value of $\alpha$.

As an additional assessment of how FSEB performs in comparison to other tests in a similar setting, we carried out a number of additional simulation studies, in which FSEB was compared with Fisher’s exact test and a score test, over various settings of $n$, $\bm{\psi}$ and $\Delta$, as well as for different ranges of $m_{i}$ ($i=1\in[n]$). Comparisons were made using the $p$-values as well as false discovery rate (FDR) corrected $p$-values arising from FDR control methods (Wang and Ramdas, 2022), and are presented in the online Supplementary Materials (Tables S1–S8 and Figures S1–S8). It is evident in almost all cases (and especially in case C, which most closely resembles the real life application scenario) that (i) the power levels are very similar across methods, especially as values of $n$, $m_{i}$ ($i\in[n]$) and effect sizes increase, and (ii) in every case, there are some settings in which Fisher’s test and the score test are anti-conservative (even after FDR correction), with their Type I error greatly exceeding the nominal levels of significance, while this never occurs for FSEB, even without FDR correction.

4 Real-data applications

4.1 The Norberg data

We now wish to apply the FSEB CI construction from Section 3.2.1 to produce CIs in a real data application. We shall investigate the Norberg data set from the REBayes package of Koenker and Gu (2017), obtained from Haastrup (2000). These data pertain to group life insurance claims from Norwegian workmen. Here, we have $n=72$ observations $\mathbf{D}_{n}$, containing total number of death claims $X_{i}$, along with covariates $\mathbf{w}_{n}$, where $w_{i}$ is the number of years of exposure, normalized by a factor of 344, for $i\in\left[n\right]$. Here each $i$ is an individual occupation group.

To analyze the data, we use the Poisson–gamma model and estimate the generative parameters $\bm{\vartheta}_{n}^{*}$ using estimates of form (13). Here, each $\theta_{i}^{*}$ can be interpreted as an unobserved multiplicative occupation specific risk factor that influences the number of claims made within occupation group $i$. To obtain individually-valid $95\%$ CIs for each of the $n$ estimates, we then apply the method from Section 3.2.1. We present both the estimated risk factors and their CIs in Figure 3.

From Figure 3, we notice that most of the estimates of $\bm{\vartheta}_{n}^{*}$ are between zero and two, with the exception of occupation group $i=22$, which has an estimated risk factor of $\theta_{22}^{*}=2.59$. Although the risk factors are all quite small, the associated CIs can become very large, as can be seen in the top plot. This is due to the conservative nature of the CI constructions that we have already observed from Section 3.1.

We observe that wider CIs were associated with observations where $X_{i}=0$, with $w_{i}$ being small. In particular, the largest CI, occurring for $i=55$, has response $X_{55}=0$ and the smallest covariate value in the data set: $w_{55}=4.45$. The next largest CI occurs for $i=5$ and also corresponds to a response $X_{5}=0$ and the second smallest covariate value $w_{5}=11.30$.

However, upon observation of the bottom plot, we see that although some of the CIs are too wide to be meaningful, there are still numerous meaningful CIs that provide confidence regarding the lower limits as well as upper limits of the underlying risk factors. In particular, we observe that the CIs for occupation groups $i=26$ and $i=54$ are remarkably narrow and precise. Of course, the preceding inferential observations are only valid when considering each of the $n$ CIs, individually, and under the assumption that we had chosen to draw inference regarding the corresponding parameter of the CI, before any data are observed.

If we wish to draw inference regarding all $n$ elements of $\bm{\vartheta}_{n}^{*}$, simultaneously, then we should instead construct a $100\left(1-\alpha\right)\%$ simultaneous confidence set $\bar{\mathcal{C}}^{\alpha}\left(\mathbf{D}_{n}\right)$, with the property that

[TABLE]

Using Bonferroni’s inequality, we can take $\bar{\mathcal{C}}^{\alpha}\left(\mathbf{D}_{n}\right)$ to be the Cartesian product of the individual $100\left(1-\alpha/n\right)\%$ (adjusted) CI for each parameter $\theta_{i}^{*}$:

[TABLE]

Using the $\alpha=0.05$, we obtain the $95\%$ simultaneous confidence set that appears in Figure 4. We observe that the simultaneous confidence set now permits us to draw useful inference regarding multiple parameters, at the same time. For example, inspecting the $n$ adjusted CIs, we observe that the occupations corresponding to indices $i\in\left\{8,22,50\right\}$ all have lower bounds above $0.5$. Thus, interpreting these indices specifically, we can say that each of the three adjusted confidence intervals, which yield the inference that the risk factors $\theta_{i}^{*}>0.5$ for $i\in\left\{8,22,50\right\}$, contains the parameter $\theta_{i}^{*}$ with probability $0.95$, under repeated sampling.

Since our individual CI and adjusted CI constructions are $e$-CIs, one can alternatively approach the problem of drawing simultaneously valid inference via the false coverage rate (FCR) controlling techniques of Xu et al. (2022). Using again the parameters $\theta_{i}^{*}$ corresponding to $i\in\left\{8,22,50\right\}$, as an example, we can use Theorem 2 of Xu et al. (2022) to make the statement that the three adjusted CIs $\mathcal{C}_{i}^{3\alpha/n}\left(\mathbf{D}_{n}\right)$, for $i\in\left\{8,22,50\right\}$, can be interpreted at the FCR controlled level $\alpha\in(0,1)$, in the sense that

[TABLE]

where $\mathbb{I}\left(\mathbf{D}_{n}\right)$ is a data-dependent subset of parameter indices. In particular, we observe the realization $\left\{8,22,50\right\}$ of $\mathbb{I}\left(\mathbf{D}_{n}\right)$, corresponding to the data-dependent rule of selecting indices with adjusted CIs $\mathcal{C}_{i}^{\alpha/n}\left(\mathbf{D}_{n}\right)$ with lower bounds greater than $0.5$. Here, \left\llbracket\mathsf{A}\right\rrbracket=1 if statement $\mathsf{A}$ is true and [math], otherwise.

Clearly, controlling the FCR at level $\alpha$ yields narrower CIs for each of our the three assessed parameters than does the more blunt simultaneous confidence set approach. In particular, the $95\%$ simultaneous adjusted CIs obtained via Bonferroni’s inequality are $\left(0.775,4.485\right)$, $\left(1.375,5.520\right)$, and $\left(0.505,3.565\right)$, and the $0.05$ level FCR controlled adjusted CIs are $\left(0.810,4.300\right)$, $\left(1.430,5.390\right)$, and $\left(0.555,3.390\right)$, for the parameters $\theta_{i}^{*}$ corresponding to the respective parameters $i\in\left\{8,22,50\right\}$. Overall, these are positive results as we do not know of another general method for generating CIs in this EB setting, whether individually or jointly.

4.2 Differential

methylation detection in bisulphite sequencing data

DNA methylation is a chemical modification of DNA caused by the addition of a methyl ($CH_{3}$-) group to a DNA nucleotide – usually a C that is followed by a G – called a CpG site, which is an important factor in controlling gene expression over the human genome. Detecting differences in the methylation patterns between normal and ageing cells can shed light on the complex biological processes underlying human ageing, and hence has been an important scientific problem over the last decade (Smith and Meissner, 2013). Methylation patterns can be detected using high-throughput bisulphite sequencing experiments (Krueger et al., 2012), in which data are generated in the form of sequences of numbers of methylated cytosines, $x_{ig}$, among the total counts of cytosines, $m_{ig}$, for $n$ CpG sites on a genome $(i\in[n])$, for $G$ groups of cell types $g\in[G]$. Often, there are $G=2$ groups, as in our example that follows, for which the question of interest is to detect regions of differential methylation in the DNA of normal and ageing cells. Based on the setup above, a set of bisulphite sequencing data from an experiment with $G$ groups might be considered as $G$ series of (possibly correlated) observations from non-identical binomial distributions. The degree of dependence between adjacent CpG sites typically depends on the genomic distance between these loci, but since these are often separated by hundreds of bases, for the moment it is assumed that this correlation is negligible and is not incorporated into our model.

4.2.1 Application to Methylation data from Human chromosome 21

We evaluated the test statistic $T_{\mathbb{I}}\left(\mathbf{D}_{n}\right)$ over a paired segment of methylation data from normal and ageing cells, from $100,000$ CpG sites on human chromosome 21 (Cruickshanks et al., 2013). After data cleaning and filtering (to remove sites with too low or too high degrees of experimental coverage, that can introduce errors), $58,361$ sites remained for analysis. Figure 5 shows the predicted demarcation of the data into differentially and non-differentially methylated sites over the entire region, at three cutoff levels of significance, overlaid with a moving average using a window size of 10 sites. It was observed that large values of the test statistic were often found in grouped clusters, which would be biologically meaningful, as loss of methylation in ageing cells is more likely to be highly region-specific, rather than randomly scattered over the genome. The overall rejection rates for the FSEB procedure corresponding to significance levels of $\alpha=0.0005,0.05,0.02$ and $0.01$ were found to be $0.0012$, $0.0154$, $0.0092$, and $0.0064$, respectively.

As a comparison to other methods for detecting differential methylation, we also applied site-by-site Fisher tests and score tests as implemented for bisulphite sequencing data in the R Bioconductor package DMRcaller (Catoni et al., 2018). For purposes of comparison, we used two significance level cutoffs of $0.05$ and $0.0005$ for our FSEB test statistic, along with the same cutoffs subject to a Benjamini–Hochberg FDR correction for the other two testing methods. Figure 6 shows the comparison between the calculated site-specific $p$-values of the Fisher and score tests with the calculated FSEB test statistic (all on the logarithmic scale) over the entire genomic segment, which indicates a remarkable degree of overlap in the regions of differential methylation. There are, however, significant differences as well, in both the numbers of differential methylation calls and their location. In particular, the FSEB test statistic appeared to have stronger evidence for differential methylation in two regions, one on the left side of the figure, and one towards the centre. The Fisher test, being the most conservative, almost missed this central region (gave a very weak signal), while the score test gave a very high proportion of differential methylation calls compared to both other methods – however, the results from the score test may not be as reliable as many cells contained small numbers of counts which may render the test assumptions invalid. Table 5 gives a summary of the overlap and differences of the results from the different methods at two levels of significance, indicating that with FDR corrections, the Fisher test appears to be the most conservative, the score test the least conservative, and the FSEB procedure in-between the two. We also calculated, for each pair of methods, the proportion of matching calls, defined as the ratio of the number of sites predicted by both methods as either differentially methylated, or non-differentially methylated, to the total number of sites. These proportions indicated a high degree of concordance, especially between FSEB and Fisher tests, with the score test showing the least degree of concordance at both levels of significance. As expected, the degree of concordance decreased with an increase in $\alpha$, but only slightly so, between the FDR-corrected Fisher test and FSEB.

5 Conclusion

EB is a powerful and popular paradigm for conducting parametric inference in situations where the DGP can be assumed to possess a hierarchical structure. Over the years, general frameworks for point estimation have been developed for EB, such as via the shrinkage estimators of Serdobolskii (2008) or the various method of moments and likelihood-based methods described in Maritz and Lwin (1989, Sec. 3). Contrastingly, the construction of interval estimators and hypothesis tests for EB parameters rely primarily on bespoke derivations and analysis of the specific models under investigation.

In this paper, we have adapted the general universal inference framework for finite sample valid interval estimation and hypothesis testing of Wasserman et al. (2020) to construct a general framework within the EB setting, which we refer to as the FSEB technique. In Section 2, we proved that these FSEB techniques generate valid confidence sets and hypothesis tests of the correct size. In Section 3, we demonstrated via numerical simulations, that the FSEB methods can be used in well-studied synthetic scenarios. There, we highlight that the methods can generate meaningful inference for realistic DGPs. This point is further elaborated in Section 4, where we also showed that our FSEB approach can be usefully applied to draw inference from real world data, in the contexts of insurance risk and the bioinformatics study of DNA methylation.

We note that although our framework is general, due to it being Markov inequality-based, it shares the same general criticism that may be laid upon other universal inference methods, which is that the confidence sets and hypothesis tests can often be conservative, in the sense that the nominal confidence level or size is not achieved. The lack of power due to the looseness of Markov’s inequality was first mentioned and discussed in Wasserman et al. (2020), where it is also pointed out that, in the universal inference setting, the logarithm of the analogous ratio statistics to (6) have tail probabilities that scale, in $\alpha$, like those of $\chi^{2}$ statistics. The conservativeness of universal inference constructions is further discussed in the works of Dunn et al. (2021), Tse and Davison (2022), and Strieder and Drton (2022), where the topic is thoroughly explored via simulations and theoretical results regarding some classes of sufficiently regular problems. We observe this phenomenon in the comparisons in Sections 3.1 (and further expanded in the Supplementary Materials). We also explored subsampling-based tests within the FSEB framework, along the lines proposed by Dunn et al. (2021), which led to very minor increases in power in some cases with small sample sizes without affecting the Type I error. With such an outcome not entirely discernible from sampling error, and with the substantial increase to computational cost, it does not seem worthwhile to employ the subsampling-based approach here. A possible reason for the lack improvement in power observed, despite subsampling, can be attributed to the fact that the sets $\mathbb{I}$, and their complements, are not exchangeable; since the indices fundamentally define the hypotheses and parameters of interest.

However, we note that since the methodology falls within the $e$-value framework, it also inherits desirable properties, such as the ability to combine test statistics by averaging (Vovk and Wang, 2021), and the ability to more-powerfully conduct false discovery rate control when tests are arbitrarily dependent (Wang and Ramdas, 2022).

Overall, we believe that FSEB techniques can be usefully incorporated into any EB-based inference setting, especially when no other interval estimators or tests are already available, and are a useful addition to the statistical tool set. Although a method that is based on the careful analysis of the particular setting is always preferable in terms of exploiting the problem specific properties in order to generate powerful tests and tight intervals, FSEB methods can always be used in cases where such careful analyses may be mathematically difficult or overly time consuming.