affil0affil0affiliationtext: Department of Statistical Science, Duke
University

# Optimal Conformal Prediction for Small Areas

Elizabeth Bersson  Peter D. Hoff

###### Abstract

Existing inferential methods for small area data involve a trade-off between
maintaining area-level frequentist coverage rates and improving inferential
precision via the incorporation of indirect information. In this article, we
propose a method to obtain an area-level prediction region for a future
observation which mitigates this trade-off. The proposed method takes a
conformal prediction approach in which the conformity measure is the posterior
predictive density of a working model that incorporates indirect information.
The resulting prediction region has guaranteed frequentist coverage regardless
of the working model, and, if the working model assumptions are accurate, the
region has minimum expected volume compared to other regions with the same
coverage rate. When constructed under a Normal working model, we prove such a
prediction region is an interval and construct an efficient algorithm to
obtain the exact interval. We illustrate the performance of our method through
simulation studies and an application to EPA radon survey data.

Keywords: Exchangeability; Hierarchical model; Kriging; Nonparametric;
Prediction region; Frequentist coverage.

##  1 Introduction

Precise and accurate inference on a sample obtained from non-overlapping
subpopulations, referred to as areas or domains, is an important goal in a
wide range of fields including economics (Berg2014), ecology (Sinha2009), and
others (Chattopadhyay1999) where localized inference for various socio-
demographic groups or refined geographic regions is of interest. In such
applications, it is common to have small area-specific sample sizes, which
presents challenges in making precise area-specific inferences via direct
methods which only make use of within-area samples. Direct methods can be
unbiased and achieve target frequentist coverage rates for all areas, but
don’t take advantage of all available information, and so may be inefficient.
As a result, researchers often turn to indirect or model-based methods that
allow information to be shared across areas. Borrowing information across
areas may decrease variability of point estimates and volume of confidence and
prediction regions, but doing so can introduce bias and thus alter area-level
frequentist coverage rates from their nominal level (Carlin1990). For more on
small area inference, see Rao2015, or, for information on multilevel modeling
more broadly, see Gelman2006a. Recently, some ‘frequentist and Bayes’ (FAB)
methods have been developed for confidence intervals that maintain area-level
frequentist error rate control and allow for information sharing (Yu2018;
Burris2020).

We focus on the task of predicting a future response in each small area.
Although it is not the focus of much of the literature, response prediction is
particularly useful for small areas. It is more general than quantifying area-
level effects and allows for a more natural comparison across different types
or levels of groups. Additionally, it may allow for more straightforward
interpretation of differences across areas than estimates of effects,
particularly in non-normal populations. More broadly, the underlying goal of
most statistical inference can be framed as being motivated by prediction
(Shmueli2010). To quote Dr. Hirotugu Akaike, “The predictive point of view is
a prototypical point of view to explain the basic activity of statistical
analysis” (Findley1995).

Consider an exchangeable sample
$\boldsymbol{Y}_{j}=\\{Y_{1,j},\ldots,Y_{n_{j},j}\\}$ from
$j\in\\{1,\ldots,J\\}$ independent subpopulations. We wish to obtain a
prediction region $A_{j}$ for a future response from area $j$,
$Y_{n_{j}+1,j}$, that is accurate, in the sense that it maintains
$(1-\alpha)100\%$ frequentist coverage,

$P_{\theta}\left(Y_{n_{j}+1,j}\in A_{j}\right)=\alpha,$

and precise, in that the expected volume is comparatively small. For data that
appear approximately Normally distributed, a commonly used direct method is
the classical pivot prediction interval. For area $j$, that is

$\bar{y}_{j}\pm z_{1-\alpha/2}\left({\sigma}^{2}(1+1/n_{j})\right)^{1/2},$

where $\bar{y}_{j}$ is the sample mean of area $j$, ${\sigma}_{j}^{2}$ is
known population variance, and $z_{q}$ is the $q$th quantile of the standard
Normal distribution. If the parametric assumptions hold true, this interval
will have the desired frequentist coverage rate. As this interval may be
prohibitively wide due to a small sample size $n_{j}$, a researcher may turn
to a Bayesian interval:

$\tilde{\theta}_{j}\pm
z_{1-\alpha/2}\left(\left(1/\tau^{2}+n/{\sigma}_{j}^{2}\right)^{-1}+\hat{\sigma}^{2}\right)^{1/2}$
(1)

where
$\tilde{\theta}_{j}:=\frac{\mu/\tau^{2}+\bar{y}n/{\sigma}^{2}}{1/\tau^{2}+n_{j}/{\sigma}^{2}}$
for prior parameters $\mu,\tau^{2}$ that are typically estimated in an
empirical Bayesian manner. This interval is narrower than the pivot by a
factor of A, and will have the nominal coverage rate on average, but, as a
result of a biased center and constant interval width, the corresponding
frequentist coverage rate declines as $|\theta_{j}-\mu|$ increases (see the
solid lines in Figure 1). Furthermore, if the parametric assumptions are not
accurate, the frequentist coverage rate of both methods suffers. Consider, for
example, a sample from a single population arising from a distribution with
mean $\theta$ and variance $1$ with a point mass equidistant on either side of
$\theta$:

$f_{y}(y)=\lambda\delta_{\theta-1}(y)+(1-\lambda)\delta_{\theta+1}(y)$ (2)

for $\lambda\sim Bernoulli(1/2)$ where $\delta_{c}$ is the Dirac delta
distribution at point $c$. Under such a sampling model, the frequentist
coverage rate attained declines drastically from the desired rate for the
pivot method, and has strange step-wise behavior dependent on
$|\theta_{j}-\mu|$ )see dashed lines in Figure 1).

Figure 1:  Frequentist coverage percent of classical pivot prediction (blue)
and Bayesian prediction (red) for parameters $n=3,\variance=1$ and Bayesian
prior parameters $\mu=0,\tau^{2}=1/2$. Results when parametric assumptions are
true (solid lines) and false (dashed lines).

When area-specific inference is of interest, such as in, for example, policy
or planning applications, maintaining the desired frequentist coverage rate
within each area may be particularly important. In this article, we develop a
non-parametric FAB prediction method that has area-level frequentist coverage
rate control while incorporating indirect information to improve prediction
region precision.

There is an extensive body of literature on estimating (or, often,
“predicting”) random effects for small areas (see, for a review, Skrondal2009
or Pfeffermann2013), but less work has been done on response prediction.
Afshartous2005 offer a review of parametric point prediction methods. The
accuracy and theoretical guarantees of these methods rely on modeling
assumptions. Other prediction methods such as those presented in Vidoni2006,
as well as empirical or fully Bayesian prediction methods produce precise
prediction intervals (Gelman2006), but do not maintain the desired coverage
level at each area.

Conformal prediction, introduced in Gammerman1998 and further developed in
Vovk2005, is a non-parametric method which relys solely on the assumption of
exchangeability to produce prediction regions guaranteed to have the desired
coverage. Candidate predicted values are included in a conformal prediction
region if they are similar to the observed data, as judged by a chosen
conformity measure. A brief review of the conformal prediction method is
included in Section 2.1. Much of the work on the conformal approach up to now
has been focused on methods for a single population (Lei2013;
Papadopoulos2011; Vovk2019). These developments may be used to construct
“direct” conformal prediction regions for each area separately, but doing so
could be inefficient, as information is not shared across groups. We construct
an information sharing conformity measure such that the corresponding expected
prediction region volume is minimized for a given error rate, as shown in
Hoff2021a. Such a prediction region, a FAB prediction region, maintains
desired frequentist coverage and may result in narrow prediction regions in
cases where auxiliary information is shared appropriately. When the proposed
conformity measure is constructed under a Normal working model, we prove the
resulting prediction region is an interval that contains a standard Bayesian
point prediction estimator. This implies a coherent method of classically
Bayesian point prediction while providing uncertainty quantification which
maintains frequentist coverage. Additionally, we develop a computationally
efficient procedure that makes full use of the data to obtain the exact FAB
conformal prediction region for each area.

This article proceeds as follows. In Section 2, we briefly review the generic
conformal prediction method and detail the motivation, computation, and
properties of the Bayes-optimal, or FAB, conformity measure for a single
population when indirect information is available. A comparison of various
prediction procedures is presented via a simulation study in Section 3. In
Sections 5 and 6, we extend the information sharing conformal algorithm to a
small area regime and illustrate our method’s performance through an
application to an EPA radon dataset. We conclude with a discussion in Section
7.

##  2 Bayes Optimal Conformal Prediction

###  2.1 Review of Conformal Prediction

Conformal prediction is a method of obtaining a prediction region for a new
observation $Y_{n+1}$ based on an exchangeable sample
$\Ybf=\\{Y_{1},...,Y_{n}\\}$ from a real-valued population. Having observed
$\Ybf=\boldsymbol{y}$, a candidate value $y_{n+1}$ of $Y_{n+1}$ is included in
the conformal prediction region if it sufficiently “conforms” to the sample,
as measured by a _conformity measure_
$C:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}$. The reader is
referred to Vovk2005 §2.1-2.2 for more details. The conformal prediction
region can be constructed to have the desired frequentist coverage rate by
including only those $y_{n+1}$-values with corresponding conformity score
$c_{n+1}$ greater than or equal to that of some fraction of the conformity
scores of the observed elements of the sample, $\\{c_{1},...,c_{n}\\}$.
Specifically, a $100(1-\alpha)\%$ prediction region for $Y_{n+1}$ can be
constructed as follows:

To determine if a candidate value $y_{n+1}$ is included in the prediction
region,

  1. 1.

compute
$c_{i}(y_{n+1}):=C\left(\\{y_{1},..,y_{n},y_{n+1}\\}\backslash\\{y_{i}\\},y_{i}\right)$
for $i=1,\ldots,n+1$;

  2. 2.

set $p_{y}:=\frac{\text{$\\#\big{\\{}i=1,...,n+1:c_{i}(y_{n+1})\leq
c_{n+1}(y_{n+1})\big{\\}}$}}{n+1}$.

The value $y_{n+1}$ is included in the region if $p_{y}>\alpha$. Note that
each conformity score $c_{i}(y_{n+1})$ is a function of the candidate
$y_{n+1}$. More compactly, the conformal prediction region may be expressed as

$A^{c}\left(\Ybf\right)=\left\\{y_{n+1}\in\mathcal{Y}:\frac{\\#\left(i=1,...,n+1:c_{i}(y_{n+1})\leq
c_{n+1}(y_{n+1})\right)}{n+1}>\frac{k}{n+1}\right\\},$ (3)

where $k=\lfloor\alpha(n+1)\rfloor$. The resulting prediction region,
$A^{C}(\Ybf)$, has greater than or equal to the desired coverage, referred to
as conservative coverage, and has exact coverage if $\epsilon=l/(n+1)$ for
some integer $l\in\\{0,1,2,...,n+1\\}$. The frequentist coverage guarantee
follows from the exchangeability assumption as all permutations of the
collection of random variables $\left\\{Y_{1},...,Y_{n+1}\right\\}$ are
equiprobable, and thus, all permutations of conformity scores
$\left\\{c_{1},...,c_{n+1}\right\\}$ are equiprobable (Balasubramanian2014
§1.3).

In summary, the conformal prediction method quantifies how well a candidate
prediction conforms with previously observed data. If the candidate is more
similar than some specified number of samples, it is included in the
prediction interval. Formally, the candidate is accepted if it’s corresponding
conformity score is greater than $k$ of the entire bag of conformity scores. A
key advantage of the conformal algorithm is that the frequentist coverage
guarantee of the conformal method holds regardless of both the true
distribution of the random variables and the choice of conformity measure.
While the coverage guarantee does not depend on $C$, the choice of $C$ does
affect the prediction region’s size. Thus, under a thoughtfully chosen
conformity measure, we can obtain a narrow prediction region that maintains
the desired coverage rate.

###  2.2 Bayes-Optimal Conformal Prediction via a Normal Working Model

Two main criteria of the usefulness of a prediction region are validity and
precision. As frequentist validity is guaranteed by the conformal method, we
focus on constructing an optimally precise prediction region through a
conformity measure that takes advantage of auxiliary or prior information.
This information will enter the conformity measure through a working model. In
practice, we often wish to make a prediction based on a random sample arising
from some unknown distribution. There may be reason to believe the population
arises from a working model, that is, some specific distribution that possibly
incorporates prior information on model parameters. In this section, we derive
a Bayes-optimal conformity measure for a single population using a Normal
working model. The conformity measure is Bayes-optimal in the sense that it
results in a prediction region with minimum expected volume for the specified
coverage rate if the working model is true, as shown in Hoff2021a. Regardless
of the accuracy of the working model, the resulting prediction region will
maintain the desired coverage rate.

Hoff2021a shows that the Bayes-optimal conformity measure is the posterior
predictive density of the model:

$C_{B}\left({\ybf},y_{n+1}\right)=p\left(y_{n+1}|\ybf\right)=\int_{\Theta}p(y_{n+1}|\theta)p(\theta|\ybf)d\theta.$
(4)

A FAB conformal prediction region, denoted $A^{fab}(\Ybf)$, can be constructed
for a population via the Bayes-optimal conformity measure obtained for a
working model that incorporates auxiliary information. Recall that the choice
of conformity measure, and, therefore, the choice of working model, does not
affect the coverage guarantee of the conformal algorithm. A well chosen
working model can, however, result in narrower prediction regions.

In this article, we focus on properties and analysis of the FAB prediction
region under a Normal working model. Suppose we think the population might be
close to Normally distributed, and prior information regarding the mean and
variance of the distribution is available. This information can be exploited
to improve the precision of $A^{fab}$ by utilizing the following Normal
working model in the construction of $C_{B}$,

$\displaystyle Y_{1},...,Y_{n}\sim{}$ $\displaystyle N(\theta,\variance)$ (5)
$\displaystyle\theta\sim{}$ $\displaystyle N(\mu,\tau^{2}\variance)$
$\displaystyle\variance\sim{}$ $\displaystyle IG(a/2,b/2).$

For now, we assume all hyperparameters of the working model,
$\\{\mu,\tau^{2},a,b\\}$, are known. In practice, they may be estimated or
otherwise obtained from auxiliary information. This will be more thoroughly
discussed in the Section 5.2. A standard calculation shows that the Bayes-
optimal conformity measure under this working model is

$C_{B}\left({\ybf},y_{n+1}\right)=p(y_{n+1}|\ybf)=dt\left(y_{n+1}|a_{\sigma},\mu_{\theta},\frac{b_{\sigma}}{a_{\sigma}}(1+\tau^{2}_{\theta})\right),$
(6)

where $dt$ is the non-standard, non-central $t$ density111 Explicitly, for
$\variance_{t}=\frac{b_{\sigma}}{a_{\sigma}}(1+\tau^{2}_{\theta})$,
$\displaystyle p(y_{i}|a_{\sigma},\mu_{\theta},\variance_{t})={}$
$\displaystyle\frac{\Gamma\left(\frac{2a_{\sigma}+1}{2}\right)}{\sqrt{2a_{\sigma}\pi}\Gamma\left(\frac{2a_{\sigma}}{2}\right)}\left(\frac{1}{\sqrt{\variance_{t}}}\left(1+\frac{1}{2a_{\sigma}}\frac{(y_{i}-\mu_{\theta})^{2}}{\variance_{t}}\right)^{-(2a_{\sigma}+1)/2}\right)$
and

$\displaystyle\tau^{2}_{\theta}=(1/\tau^{2}+n)^{-1}{}$
$\displaystyle\mu_{\theta}=(\mu/\tau^{2}+\onent\ybf)\tau^{2}_{\theta}$
$\displaystyle a_{\sigma}=a+n{}$ $\displaystyle
b_{\sigma}=b+\ybf^{T}\ybf+\mu^{2}/\tau^{2}-\left(\tau^{2}_{\theta}\right)^{-1}\mu_{\theta}^{2}.$

For the methodology that follows, it will be useful to formalize a definition
of equivalent conformity measures (ECM). As we will show, computations may be
facilitated by considering alternative representations of a conformity
measure, or ECMs. Formally,

######  Definition 1 (equivalent conformity measure).

Two conformity measures are called equivalent conformity measures (ECM) if the
resulting conformal prediction regions are equivalent.

The idea of an ECM and its computational benefits have been discussed before
in the conformal literature. For example, standard computation of the
prediction region resulting from the popular non-conformity measure distance
to the average, $C_{avg}(\ybf,y_{n+1})=\left|y_{n+1}-\overline{\ybf}\right|,$
typically requires re-computing the mean of a set $n+1$ times during the
execution of the conformal algorithm. As discussed in Shafer2008, this can be
avoided upon consideration of an ECM,
$C_{avg}\left(\\{\ybf,y_{n+1}\\},y_{n+1}\right)$. Considering this
representation in place of the classically defined measure allows each
conformity score to be defined in terms of the sample mean and the unknown
candidate $y_{n+1}$ which can simplify computation. In more complex settings
such as with the Bayes-optimal measure we consider, the computational gain
obtained from constructing an algorithm under an ECM may be substantial. Our
conformity measure has the same property results under the hierarchical Normal
working model. That is,

######  Theorem 1.

Under Model 5, $C_{B}\left(\ybf,y_{n+1}\right)$ and
$C_{B}\left(\\{\ybf,y_{n+1}\\},y_{n+1}\right)$ are ECM.

All proofs are contained in Appendix 1.1. Given Theorem 1, the structure of
the fundamental inequality in the conformal algorithm,

$C_{B}\left(\\{\boldsymbol{y},y_{n+1}\\},y_{i}\right)\leq
C_{B}\left(\\{\boldsymbol{y},y_{n+1}\\},y_{n+1}\right),$

is greatly simplified as each side of the inequality is a $t$ density with the
same parameters. Specifically, this inequality simplifies to an expression
quadratic in the candidate $y_{n+1}$ which aids in simplifying computation and
proofs of properties of the prediction region.

###  2.3 Computation and Properties

In order to obtain the conformal prediction region, in principle, the
conformal algorithm must be evaluated for each candidate in the sample space
$\mathcal{Y}$. As such, unless $\mathcal{Y}$ is a finite set, the conformal
prediction method could, in general, be prohibitively computationally
expensive. However, by making use of properties of the form of $C_{B}$ under
the Normal working model (Equations 5), it turns out the exact conformal
prediction region can be obtained by a procedure that involves evaluating a
simple function of the sample. Additionally, we prove the FAB conformal
prediction region under the Normal working model will be an interval that
contains the posterior mean estimator of the population mean,
$\tilde{\theta}:=(\mu/\tau^{2}+\sum_{k=1}^{n}y_{k})(1/\tau^{2}+n)^{-1}$.

The FAB conformal prediction region can be obtained via a two step process.
First, for each $i=1,...,n+1$, find the sub-region of acceptance,

$S_{i}:=\left\\{y_{n+1}\in\mathcal{Y}:c_{B,i}(y_{n+1})\leq
c_{B,n+1}(y_{n+1})\right\\}.$ (7)

Then, translate the information in the set $\\{S_{1},...,S_{n+1}\\}$ to obtain
the number of $i=1,...,n+1$ such that $c_{i}\leq c_{n+1}$ at each point in the
domain $\mathcal{Y}$. As made clear by the representation of a generic
conformal prediction region given in Equation 3, this information fully
classifies the conformal prediction region for a given error rate $\epsilon$.

This process is visualized for a dummy sample of size $n=4$ in Figure 2. For
clarity, the conformity scores for each value in the sample and the candidate
prediction are plotted as a function of the candidate prediction in panel (a).
The corresponding sub-regions of acceptance are the regions in the sample
space where each conformity score is less than or equal to the conformity
score of the candidate. Under the Normal working model, each sub-region of
acceptance, plotted in panel (b), is an interval that contains
$\tilde{\theta}$. Finally, this information can be directly translated to the
number of conformity scores less than or equal to the candidate conformity
score over the sample space. Dividing this value by $(n+1)$ yields the
conformal p-value, $p_{y}$ (shown in panel (c)). From Figure 2(c), it is easy
to see for a prediction error rate of, for example, $\alpha=0.2$, the
resulting conformal prediction region is the region where $\\#(i:c_{i}\leq
c_{n+1})>1$, or, $[-3.1,2.4]$.

Figure 2: Visualization of the process to obtain the FAB conformal region: (a)
conformity scores for each value in the sample (red dashed curves) and the
candidate (thick black curve) over the sample space; (b) sub-regions of
acceptance corresponding to the conformity scores; (c) number of conformity
scores less than or equal to the candidate conformity score over the sample
space; vertical black lines are drawn at $\tilde{\theta}$.

Given the standard form of the Bayes-optimal conformity measure (Equation 6),
the sub-regions of acceptance are difficult to obtain analytically as the
candidate $y_{n+1}$ appears non-linearly in each
$c_{B,1}(y_{n+1}),...,c_{B,n+1}(y_{n+1})$. Upon consideration of the
equivalent representation of $C_{B}$ given in Theorem 1, the regions
$S_{1},...,S_{n}$ can be expressed in closed form. This allows for efficient
and exact computation of the prediction region and, in turn, can be used to
prove the FAB prediction region is an interval under the Normal working model.
These results are formalized below. All proofs are contained in Appendix 1.1.

We first present two lemmas (1 and 2) used to prove a conformal prediction
region is an interval. If the conformal p-value $p_{y}$ is a step function
over the domain $\mathcal{Y}$ with a symmetric number of unit steps to and
from $1/(n+1)$ and $1$, as in Figure 2(c), then the prediction region will be
an interval. The following two lemmas may be used to prove this is the case.

######  Lemma 1.

Under the conformal algorithm, if for each $i=1,...,n+1$, the region

$\\{y_{n+1}\in\mathcal{Y}:c_{i}(y_{n+1})\leq c_{n+1}(y_{n+1})\\}$

is an interval and contains some common value, then
$f(y_{n+1})=\\#\\{i\in\\{1,...,n+1\\}:c_{i}\leq c_{n+1}\\}$ is a step function
over the domain $\mathcal{Y}$ with ordered range
$\\{1,2,...,n,n+1,n,...,2,1\\}$.

######  Lemma 2.

Under the conformal algorithm, if $f(y_{n+1})$ is a step function over the
domain $\mathcal{Y}$ with ordered range $\\{1,2,...,n,n+1,n,...,2,1\\}$, then
the resulting prediction region is an interval.

It turns out the hypothesis of Lemma 1 holds for the Normal working model case
under $C_{B}$. Specifically, by utilizing the equivalent form of $C_{B}$ given
by Theorem 1, we are able to conclude that each $S_{i}$ is an interval for
$i\in\\{1,...,n\\}$ and obtain the closed form expression of the bounds. This
result is presented in Lemma 3.

######  Lemma 3.

Under Model 5, for each $i=1,...,n$, the region

$S_{i}=\\{y_{n+1}\in\mathcal{Y}:C_{B}\left(\\{\boldsymbol{y},y_{n+1}\\}\backslash
y_{i},y_{i}\right)\leq C_{B}\left(\\{\boldsymbol{y},y_{n+1}\\}\backslash
y_{n+1},y_{n+1}\right)\\},$

is an interval
$\left[\min\\{y_{i},g(y_{i})\\},\max\\{y_{i},g(y_{i})\\}\right]$ where

$g(y_{i}):=\frac{2\left(\mu/\tau^{2}+\sum_{k\in\\{1:n\\}}y_{k}\right)(1/\tau^{2}+n+1)^{-1}-y_{i}}{1-2(1/\tau^{2}+n+1)^{-1}}.$

Additionally, the posterior mean estimator of the population mean,
$\tilde{\theta}$, is contained in each sub-region of acceptance
$S_{1},...,S_{n}$. Not only is this result useful for proving the prediction
region is an interval, but it also suggests that using $\tilde{\theta}$ as the
estimator for a new prediction and taking $A^{fab}(\Ybf)$ as the prediction
interval is a coherent method to predict the unknown value in a classically
Bayesian manner while providing uncertainty quantification that maintains the
desired frequentist coverage. Formally,

######  Lemma 4.

For each $i=1,...,n+1$, the interval

$S_{i}=\\{y_{n+1}\in\mathcal{Y}:C_{B}\left(\\{\boldsymbol{y},y_{n+1}\\}\backslash
y_{i},y_{i}\right)\leq C_{B}\left(\\{\boldsymbol{y},y_{n+1}\\}\backslash
y_{n+1},y_{n+1}\right)\\},$

contains the posterior mean estimator of the population mean,
$\tilde{\theta}:=(\mu/\tau^{2}+\sum_{k=1}^{n}y_{k})(1/\tau^{2}+n)^{-1}$.

In total, these results prove the Bayes-optimal conformal prediction region
under the Normal working model is be an interval that contains
$\tilde{\theta}$. Formally,

######  Theorem 2.

The prediction region obtained from the conformal algorithm will be an
interval under Model 5 and conformity measure $C_{B}$.

The nature of $C_{B}$ under the Normal working model suggests a simple,
efficient algorithm to obtain $A^{fab}(\Ybf)$. In particular, as $p_{y}$ is an
incremental step function over the sample space characterized by a symmetry in
the number of steps on either side of $p_{y}=1$ (e.g. Figure 2(c)), the
corresponding prediction region can be obtained by taking the $k$th ordered
step location from either end of the collection of step locations, where
$k=\lfloor\alpha(n+1)\rfloor$. Specifically, collect the bounds of each sub-
region of acceptance $S_{1},...,S_{n}$ in a vector $\wbf$. Recall that the
desired prediction error $\epsilon$ enters in the conformal algorithm through
$k$. The region in $\mathcal{Y}$ where $\\#\\{c_{i}\leq c_{n+1}\\}\geq k+1$ is
the $k$th ordered entry of $\wbf$ and the reversed entries of $\wbf$. To
summarise, to obtain $A^{fab}(\Ybf)$,

  1. 1.

for each $i=1,...,n$, obtain $y_{i},g(y_{i})$, the two critical values of
$S_{i}$,

  2. 2.

set
$\wbf=\begin{bmatrix}y_{1}&\cdots&y_{n}&g(y_{1})&\cdots&g(y_{n})\end{bmatrix}^{T}$
and $k=\lfloor\alpha(n+1)\rfloor$,

  3. 3.

acquire the bounds of the prediction region, the $k$th and $(2n-k+1)$th order
statistics of $\wbf$.

Then, the Bayes-optimal conformal prediction region with $(1-\alpha)100\%$
coverage is

$A^{fab}(\Ybf)=\left(\wbf_{(k)},\wbf_{(2n-k+1)}\right).$

A brief note that if $y_{i}=g(y_{i})$ for at least one $i\in\\{1,...,n\\}$,
the resulting prediction region will be a singular point for certain error
rates. In the case of continuous random variables, this happens with
probability 0. double check.

##  3 Numerical Comparisons

To demonstrate properties of the FAB prediction method and compare to
alternative prediction methods, we numerically evaluate prediction region
properties resulting from various populations.horrible. The benefit of
incorporating prior or indirect information in the construction of prediction
intervals as in the Bayes-optimal conformal method is most clearly seen
through a comparison with a non-parametric method with the same frequentist
coverage guarantee, such as the distance to average conformal prediction
method. The difference between these two methods is the ability to utilize
prior information in the construction of the prediction interval. Our
numerical results compare the risk and Bayes risk of the FAB conformal and
distance to average conformal prediction intervals, where the risk function is
taken to be expected interval width $E_{\theta}[|A(\Ybf)|]$.

The prior parameters used in the FAB conformal method are $\\{\mu,\tau^{2}\\}$
and respectively represent the prior expected population mean and confidence
around this expectation. To assess the effect of these prior parameters, we
allow them to vary in this study. Additionally, we compare simulation results
for various sample sizes $n\in\\{3,7,11,15,19\\}$ and consider a prediction
error rate of $0.25$, chosen such that the conformal methods will result in
regions with exact coverage. In general, we expect the FAB conformal method to
outperform competing methods when there is concentrated and accurate prior
information regarding the mean of the population, but a limited amount of
information in the sample itself. More specifically, the FAB conformal
interval can be expected to produce narrower intervals than default methods
when $|\theta-\mu|$ is small, $\tau^{2}$ is small, and $n$ is small, or a
combination of these properties.

Figure 3: Expected width ratio of FAB conformal interval to distance to
average conformal interval for increasing $n\in\\{3,7,11,15,19\\}$ in
decreasing darkness: (a) expectation taken with respect to $Y$ and $\theta$;
(b) expectation taken over $Y$ conditional on $\theta$ for
$\mu=0,\tau^{2}=1/2$; (c) same as (b) for $\tau^{2}=2$.maybe add bounds image

The ratio of expected interval widths of the FAB method relative to the
distance to average method are displayed in Figure 3. Recall that the FAB
method incorporates prior information while the distance to average method
does not. We compute the expected interval widths via Monte Carlo
approximation using 25,000 independently generated replications for each
combination of values of $\theta,\tau^{2},$ and $n$. The effects of sample
size and prior variance of the population mean are the focus of Figure 3(a).
This figure plots the ratio of Bayes expected interval widths of the FAB
conformal to distance to average conformal intervals, where the expectations
are taken with respect to $Y$ and $\theta$ under the prior $\theta\sim
N(\mu,\tau^{2})$. The Bayes risks of the conformal intervals do not depend on
$\mu$. As informed by the theoretical results of FAB conformal prediction
(Hoff2021a), the Bayes risk of the FAB interval is smaller than that of the
distance to average interval, with the overall deviation between the methods’
expected interval widths decreasing as the sample size increases. The FAB
interval is substantially narrower than the distance to average interval for a
wide range of $\tau^{2}$ values under very small sample sizes. Intuitively,
for small sample sizes, even a low level of confidence around the prior value
of the population mean $\mu$ is useful information and will translate to
narrower prediction intervals if utilized in the construction of $C_{B}$. More
confidence, as conveyed through a smaller $\tau^{2}$ value, translates to a
more substantive increase in precision. Even for larger sample sizes, a
nontrivial gain in precision occurs under small (less than about 0.5)
$\tau^{2}$ values, representing very informative prior information about the
population mean through a prior with tight concentration around $\mu$.

For a given concentration level $\tau^{2}$, we may wish to understand how the
accuracy of our mean estimate $\mu$ affects the resulting FAB prediction
interval width. Figures 3(b)-(c) display the expected interval width ratio of
the FAB conformal to the distance to average conformal method for varying
population means and sample sizes when $\mu=0$ and $\tau^{2}=1/2,2$. Under
this set-up, by the Bayes-optimal property of the conformal measure, the FAB
interval will outperform alternatives when $|\theta-\mu|\approx 0$ and
$\tau^{2}$ is very small. The numerical results match this conclusion, and, as
the distance between the prior mean and the population mean increases in
absolute value, the FAB intervals become wider, and thus the benefit of
utilizing this type of prior information declines. As seen in panel (b), for a
moderately small $\tau^{2}$ value, and for the smallest sample size
considered, the FAB conformal method results in an interval width that is
17.6% narrower than the competitor when $\theta=\mu$ exactly. For larger
$\tau^{2}$, as seen in the panel (c), there is a less substantial benefit to
utilizing this auxiliary information in $C_{B}$, but some benefit is seen
nonetheless for a wider range of $\theta$ divergences from $\mu$.

Overall, when auxiliary information is available to accurately inform
conformal parameter values, FAB prediction intervals outperform the commonly
utilized prediction intervals considered as judged by precision and coverage
guarantees. For a multiple population regime, when sharing information
correctly informs about the heterogeneity of the various groups, there will be
a benefit to exploiting this information in obtaining prediction intervals.
The benefit is particularly stark when obtaining prediction intervals based on
small sample sizes.return

##  4 DELETE

Figure 4: Comparison of parametric and nonparametric prediction regions
obtained for parameters $n=3,\mu=0,\tau^{2}=1/2,\theta\in[-2.5,2.5]$: (a)
expected interval width; (b) prediction interval endpoints. Black solid lines
represent the FAB conformal method, blue dotted lines represent the distance
to average conformal method, and thick grey lines represent the classical
pivot method. Diagonal black dashed line represents the 45°line.

The classical frequentist prediction interval based off of a pivotal quantity
results in prediction regions with the desired frequentist coverage rate if
the sampling model is correctly specified. In the case of assumed Normality,
the corresponding pivot prediction interval is

$\bar{y}\pm t_{1-\alpha/2,n-1}\left(\hat{\sigma}^{2}(1+1/n)\right)^{1/2},$

where $\bar{y}$ is the sample mean, $\hat{\sigma}^{2}$ is the sample variance,
and $t_{q,d}$ is the $q$th quantile of the standard $t$ distribution with $d$
degrees of freedom. The expected pivot interval is symmetric about the
population mean $\theta$ and maintains a constant interval width regardless of
$\theta$. The resulting pivot prediction interval will be slightly narrower
than the distance to average interval as a result of taking into account the
underlying parametric structure. The FAB conformal method produces narrower
prediction intervals than both when the prior information utilized in $C_{B}$
is accurate, as characterized by a small $|\theta-\mu|$ (Figure 4(a)).

The FAB prediction interval exhibits standard ’Bayesian’ behavior by shifting
the interval towards the prior information, in this case, $\mu=0$. As this
prior information declines in accuracy, the FAB interval must become wider to
maintain the desired frequentist coverage rate. As such, the benefit of
utilizing prior information in the construction of a prediction interval
decreases as the accuracy of that information declines (Figure 4). If the
underlying assumption of Normality holds true for the population considered,
as in these experiments, the pivot method has the same Frequentist coverage
rate as the non-parametric methods, so any gain in interval precision from the
FAB method is desirable.all over place

Oftentimes, particularly in small sample size situations, it is difficult to
draw accurate parametric conclusions regarding a population.return Consider,
for example, a mixture model with mean $\theta$ and variance 1 with a point
mass equidistant on either side of $\theta$:

$f_{x}(x)=\lambda\delta_{\theta-1}(x)+(1-\lambda)\delta_{\theta+1}(x)$ (8)

for $\lambda\sim Bernoulli(1/2)$ where $\delta_{c}$ is the Dirac delta
distribution at point $c$. For a very small sample size, it may be difficult
to distinguish this density from a Normal density with mean $\theta$ and
variance 1. The corresponding oracle pivot interval for this mean and variance
combination can, however, be within the two possible point masses resulting in
a frequentist coverage rate of 0. Consider, for example, $n=2,\alpha=1/2$. The
resulting pivot interval is

$\left(\theta-
z_{1-\alpha/2}\left(1+1/n\right)^{1/2},\theta+z_{1-\alpha/2}\left(1+1/n\right)^{1/2}\right)\approx\left(\theta-0.826,\theta+0.826\right).$

Very small sample sizes like this are quite common in practice, and in small
area applications in particular (cite), so the non-parametric coverage
guarantee of the FAB conformal method is especially appealing when applied in
a similar scenario.

###  4.1 Comparison with Indirect Parametric Prediction

Figure 5: Comparison of Bayes and FAB prediction intervals obtained for
parameters $n=3,\mu=0,\tau^{2}=1/2,\theta\in[-2.5,2.5]$: (a) prediction
interval endpoints. Black solid lines represent the FAB conformal method and
red dashed lines represent the empirical Bayesian method. Diagonal black
dashed line represents the 45°line. (c) empirical Bayesian coverage rate.
intervals

As a final comparison, we evaluate differences in properties between the FAB
conformal interval and Bayesian prediction, both of which utilize prior
information, but have differing frequentist coverage guarentees. The Bayesian
prediction interval based on a Normal population with prior $\theta\sim
N(\mu,\tau^{2})$ is

$\tilde{\theta}\pm
t_{1-\alpha/2,n-1}\left(\left(1/\tau^{2}+n/\hat{\sigma}^{2}\right)^{-1}+\hat{\sigma}^{2}\right)^{1/2}$
(9)

where
$\tilde{\theta}:=\frac{\mu/\tau^{2}+\bar{y}n/\hat{\sigma}^{2}}{1/\tau^{2}+n/\hat{\sigma}^{2}}$.
This interval is narrower than the pivot by a factor of A and has constant
interval width for a given prior concentration $\tau^{2}$. As a result of a
biased center and constant interval width, the corresponding frequentist
coverage rate varies with $|\theta-\mu|$ (5).

Panel (b) displays the corresponding interval bounds for each prediction
method. In this small $\tau^{2}$, small sample size case, for small
$|\theta-\mu|$, the FAB prediction interval has been shown to be notably more
precise on average than the distance to average conformal interval.

As $|\theta-\mu|$ grows, the FAB interval width blows up in a manner that the
Bayesian interval doesn’t. If maintaining the desired coverage regardless of
accuracy of prior information is a priority, this may be more desirable than
the behavior of the empirical Bayesian interval which maintains the same
interval width regardless of the inaccuracy of the prior information
available, as evidenced by a decline in coverage (Figure
LABEL:single_pop_allmethods(c)). As seen in panel (b), the FAB conformal
bounds follow a similar pattern as the empirical Bayesian bounds, albeit more
extreme. Both methods are shifted away from the population mean $\theta$ and
towards the prior hypermean $\mu=0$.

##  5 FAB Small Area Prediction

###  5.1 Information Sharing Working Model

In inference on small areas, utilizing indirect methods that share information
across areas has been shown to improve precision compared with direct methods,
particularly for areas with small sample sizes (Gelman2006). With this
motivation, we extend the FAB conformal prediction method explored thus far to
a small area regime. In the construction of the Bayes-optimal conformity
measure, information is shared across groups via a multilevel working model in
order to increase prediction region precision while maintaining area-level
frequentist coverage guarantees.

For group $j\in\\{1,...,J\\}$, we observe an exchangeable sample
$(Y_{1,j},...,Y_{n_{j},j})=(y_{1,j},...,y_{n_{j},j})$ of length $n_{j}$ such
that the samples are independent across groups. Suppose a reasonable working
model for the populations is a spatial Fay-Herriot model (Fay1979).
Specifically, for each area $j\in\\{1,...,J\\}$, the working sampling model is
a Normal distribution with area-specific mean and variance,
$Y_{1,j},...,Y_{n_{j},j}\sim{}N\left(\theta_{j},\variance_{j}\right)$, and the
working linking model for the area-level parameters is

$\displaystyle\thetabf\sim{}$ $\displaystyle
N\left(\boldsymbol{X}\boldsymbol{\beta},\eta^{2}\boldsymbol{G}\right)$ (10)
$\displaystyle\variance_{1},...,\variance_{J}\sim{}$ $\displaystyle
IG(a/2,b/2),$

where $\boldsymbol{G}$ is a spatial covariance matrix such as, e.g., that
which results from the popular simultaneous (SAR) autoregressive model
$\boldsymbol{G}=[(\boldsymbol{I}-\rho\boldsymbol{W})(\boldsymbol{I}-\rho\boldsymbol{W}^{T})]^{-1}$
(Singh2005). $\boldsymbol{W}$ is a distance matrix among areas that is
typically row-standardized to sum to 1, and $\rho\in(-1,1)$ is a spatial
correlation parameter. For more on spatial modeling, see Banerjee2014. This
flexible set-up allows for inclusion of an array of auxiliary information
including area-level covariates $\boldsymbol{X}$ and spatial relationships in
the linking model for the population means which can be exploited to improve
precision of prediction regions.

For population $j$, we may estimate unknown parameters in the working model
(Eqns 10) given all available data outside of area $j$, and use this
information to inform unknown parameters in the Bayes optimal conformity
measure. A FAB conformal predictive interval may then be constructed given
this information. By the assumed independence across populations, the form of
the Bayes optimal conformity measure for population $j$ is as given in
Equation 6. The resulting prediction region is constructed from a measure that
shares information and will maintain the desired Frequentist coverage for each
area.

###  5.2 FAB Conformal Parameter Estimation Procedure

All that remains to proceed with the FAB conformal method is to acquire values
for the unknown prior parameters needed for the Bayes-optimal conformity
score. For area $j$, as a result of Lemma 3, values are only required for
$\mu$, the prior mean of area $j$’s population mean $\theta_{j}$, and
$\tau^{2}$, the ratio of the prior variance of $\theta_{j}$ to population
$j$’s variance $\variance_{j}$. If the parameters of the working model are
known, we can take $\\{\mu,\tau^{2}\\}$ to be the conditional mean and
conditional variance proportion of $\theta_{j}$ and proceed with evaluation of
the Bayes-optimal conformal algorithm. In practice, these values are not
known, but they can be estimated via standard techniques from the spatial Fay-
Herriot literature. The theoretical guarantees of the conformal algorithm do
not depend on the accuracy of the conformity measure, and thus are upheld
regardless of the values of the parameters $\\{\mu,\tau^{2}\\}$. However, as
conveyed in the simulation study and proven theoretically in Hoff2021a, if the
values of the estimates are an accurate representation of heterogeneity across
groups, then the FAB conformal prediction intervals achieve improved precision
over alternative prediction regions with the same coverage. As a result, we
suggest considering computational efficiency in determining an estimation
method.

We propose an empirical Bayesian approach whereby values of the prior
parameters $\\{\mu,\tau^{2}\\}$ are obtained for each group $j$ using samples
from all populations except for population $j$. As an overview, for area $j$,
our estimation procedure proceeds with first computing estimates of unknown
parameters in the working model using all external data. Then, given these
estimates, obtaining the conditional mean of $\theta_{j}$ and the proportion
of the conditional variance of $\theta_{j}$ to an estimate of area $j$s
population variance, which are taken as $\mu$ and $\tau^{2}$, respectively.

In more detail, first, obtain the maximum likelihood estimates (MLEs) of $a,b$
by maximizing the marginal density of $\\{S_{k}^{2}\\}_{k\in K}$ for
$K=\\{1,...,J\\}\backslash j$ where
$S^{2}_{k}=\sum_{i=1}^{n_{k}}(\Ybf_{ik}-\bar{\Ybf}_{k})$ . Under the
assumptions of the working sampling model, the marginal density for the entire
population can be shown to be

$p(s_{1}^{2},...,s_{J}^{2}|a,b)=\prod_{j=1}^{J}f(s_{j}^{2})\frac{\Gamma\left(\frac{a+n_{j}-1}{2}\right)\left(\frac{b}{2}\right)^{a/2}}{\Gamma\left(\frac{a}{2}\right)\left(\frac{b+s^{2}_{j}}{2}\right)^{(a+n_{j}-1)/2}}$

for a function $f$ that does not depend on the hyperparameters $a,b$. Use the
resulting MLEs to obtain empirical Bayes estimates of each area’s population
variance. That is, take
$\hat{\sigma}^{2}_{k}=(\hat{b}+s_{k}^{2})(\hat{a}+(n_{k}-1)+1)$ for
$k\in\\{1,....K\\}$ and $\hat{\sigma}_{j}^{2}=\hat{b}/(\hat{a}+1)$. Taking
$\\{\hat{\sigma}_{k}^{2}\\}$ as plug-in values of the population variances,
obtain the marginal maximum likelihood estimates of the mean prior
hyperparameters $\boldsymbol{\beta},\rho$, and $\eta^{2}$ through standard
REML or ML procedures (Pratesi2008). Maximum likelihood estimates
$\\{\hat{\boldsymbol{\beta}},\hat{\eta}^{2},\hat{\rho}\\}$ of
$\\{\boldsymbol{\beta},\eta^{2},\rho\\}$ may then be used as to obtain an
empirical Bayes estimate $\hat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta}$.

For area $j$, obtain estimates
$\\{\hat{\boldsymbol{\beta}},\hat{\eta}^{2},\hat{\rho},\hat{\boldsymbol{\theta}},\\{\hat{\sigma}_{k}^{2}\\}_{k\in
K},\hat{\sigma}_{j}^{2}\\}$ as specified above using all data outside of area
$j$. Then, take the conformal parameters to be the conditional mean and the
proportion of conditional variance of $\theta_{j}$ obtained under these
estimates:

$\displaystyle{\mu}_{j}={}$ $\displaystyle
E[\theta_{j}|\boldsymbol{\theta}_{-j}]$ (11) $\displaystyle={}$
$\displaystyle\boldsymbol{x}_{j}^{T}\hat{\boldsymbol{\beta}}+\boldsymbol{V}_{[j,-j]}\boldsymbol{V}_{[-j,-j]}^{-1}\left(\hat{\boldsymbol{\theta}}_{-j}-\boldsymbol{X}_{[-j]}\hat{\boldsymbol{\beta}}\right)$
$\displaystyle{\tau}_{j}^{2}={}$ $\displaystyle
Var[{\theta}_{j}|{\boldsymbol{\theta}}_{-j}]/\hat{\sigma}^{2}_{j}$ (12)
$\displaystyle={}$
$\displaystyle\left({\boldsymbol{V}}_{[j,j]}-{\boldsymbol{V}}_{[j,-j]}{\boldsymbol{V}}_{[-j,-j]}^{-1}{\boldsymbol{V}}_{[-j,j]}\right)/\hat{\sigma}^{2}_{j}$

where
$\boldsymbol{V}=\hat{\eta}^{2}[(\boldsymbol{I}-\hat{\rho}\boldsymbol{W})(\boldsymbol{I}-\hat{\rho}\boldsymbol{W}^{T})]^{-1}$.
Finally, given these values for the conformal parameters obtained from
information independent of area $j$, the conformal algorithm proceeds as
described in Section 2.3. For each area, the algorithm yields an interval that
may have improved precision over other methods as a result of information
sharing and maintains the specified frequentist coverage rate.

##  6 Radon Data Example

In 1985, indoor radon levels were identified as a significant risk factor to
public health. Shortly thereafter, the U.S. Environmental Protection Agency
(EPA) conducted a national survey in an effort to better understand the risk
(USEnvironmentalProtectionAgency1992). Indoor radon levels were measured from
a stratified random sample of households across 9 states of interest. We limit
our scope to a subset of the available data and explore the sample collected
in Minnesota. This sample consists of 919 observations in total throughout the
state’s 85 counties. Within-county sample sizes range from 1 to 116, 46% of
which are 5 or less.

Price1996 explored modeling the household radon levels in Minnesota with a
goal of improving estimated county-level means, and accurate county-specific
predictions are cited as being of particular interest. Due to the frequency of
small within-county sample sizes, these are difficult tasks. Given the
abundance of auxiliary information including county-level covariates and
apparent spatial relationships among radon levels across counties, sharing
information in the construction of confidence or prediction intervals is a
natural tool to improve inferential precision. In what follows, we compare
properties of prediction intervals resulting from direct, indirect, and FAB
methods. As county-specific predictive inference is of primary interest, an
ideal prediction interval will have the desired $(1-\alpha)100\%$ coverage
while maintaining an interval width that is practically informative.

Log radon levels, $Y_{ij}$ for household $i$ in county $j$, appear
approximately Normally distributed, so we assume a Normal working sampling
model $Y_{ij}\sim N(\theta_{j},\variance_{j})$, independently across both
household $i\in\\{1,...,n_{j}\\}$ and county $j\in\\{1,...,J\\}$. $\theta_{j}$
and $\variance_{j}$ represent the subpopulation mean and variance,
respectively, of log-radon levels in homes within county $j$. We are
interested in obtaining a prediction interval for radon at one new randomly
sampled household in each county, $Y_{n_{j}+1,j}$. To illustrate the effect of
sharing information across counties in constructing prediction regions, we
compute FAB prediction intervals for each county $j$ by utilizing the spatial
Fay-Herriot model as a working linking model (Equations 10). Estimates of the
conformal parameters are obtained using data from all counties outside of $j$,
as described in Section 5.2. We include as a covariate a county-level soil
uranium measurement and allow for spatial effects under the (row-standardized)
squared exponential distance matrix between county centroids. That is, before
row-standardization, the matrix entries are
$\\{w_{ik}\\}=e^{-||x_{i}-x_{k}||^{2}}$ for counties $i,k$. The resulting
Bayes-optimal conformal prediction interval that incorporates covariates and
spatial information is labeled FAB-S in figures and tables. Additionally, we
obtain FAB conformal prediction intervals under a standard multilevel working
model that does not incorporate covariate or spatial auxiliary information.
That is, each population mean is linked through a shared prior mean and
variance ($\theta_{j}\sim N(\mu,\variance_{j}\tau^{2})$, for
$j\in\\{1,...,J\\}$).

| % Wider than FAB | % Wider than FAB-S
---|---|---
Piv | 58.5% | 62.2%
EB | 51.2% | 57.0%
C- Avg | 63.4% | 68.3%
Table 1: Comparisons of interval width of common prediction methods with
various Bayes-optimal conformal intervals from the radon data study: FAB
intervals constructed under standard multilevel model and FAB-S intervals
incorporate a county level covariate and spatial information.

We compare the Bayes-optimal conformal intervals to a few commonly used
prediction methods. For each county $j$, with prediction error rate
$\epsilon_{j}$, we obtain the classical direct prediction method, the pivot
interval (Piv) of the form $\bar{y}_{j}\pm
t_{1-\alpha_{j}/2,n_{j}-1}\left(\hat{\sigma}^{2}_{j}(1+1/n_{j})\right)^{1/2}$
where $\bar{y}_{j}$ and $\hat{\sigma_{j}}^{2}$ are the $j$th county sample
mean and variance, respectively. Additionally, we obtain direct conformal
prediction regions that only make use of within-county information via the
distance to average non-conformity measure (C-Avg). We may also obtain an
empirical Bayesian (EB) prediction interval, $\tilde{\theta}_{j}\pm
t_{1-\alpha_{j}/2,n_{j}-1}\left(\left(1/\hat{\tau}_{-j}^{2}+n_{j}/\hat{\sigma}_{j}^{2}\right)^{-1}+\hat{\sigma}_{j}^{2}\right)^{1/2}$,
where
$\tilde{\theta}_{j}:=\frac{\hat{\mu}_{-j}/\hat{\tau}_{-j}^{2}+\bar{y}_{j}n_{j}/\hat{\sigma}_{j}^{2}}{1/\hat{\tau}_{-j}^{2}+n_{j}/\hat{\sigma}_{j}^{2}}$,
and $\hat{\mu}_{-j},\hat{\tau}_{-j}^{2}$ are taken to be MLEs of obtained with
all data outside of county $j$ under the classical mean multilevel model with
a linking prior of $\theta_{1},...,\theta_{J}\sim N(\mu,\tau^{2})$. The
intervals arising from all methods are obtained for counties with sample sizes
greater than 1 under a county-specific error rate
$\alpha_{j}=\lfloor\frac{1}{3}(n_{j}+1)\rfloor/(n_{j}+1)$ to allow for within-
county consistency across the various prediction intervals considered.

A comparison of interval widths is presented in Table 1. In summary, including
relevant auxiliary information in the construction of prediction intervals via
the Bayes-optimal conformal procedure results in increased interval precision
overall. All prediction methods considered are wider than the FAB and FAB-S
prediction intervals for over half of the counties. As more relevant auxiliary
information is included, the percent wider only increases. In fact, the
intervals constructed without spatial or covariate information are wider than
those with this additional information for 69.5% of counties.

Figure 6: County-level radon prediction intervals: black solid lines represent
the FAB conformal intervals and thick grey lines represent the distance to
average conformal intervals. Dotted lines drawn at the state-wide sample mean
$\sum\bar{y}_{j}/J$ and 45° line.

In Figure 6, the prediction intervals obtained from the FAB-S conformal
method, the pivot method, and the empirical Bayesian method are plotted.
Previous analysis informs that Bayes-optimal prediction intervals will be
narrower than alternative intervals for counties where the working linking
model assumptions are closest to truth. That is, the FAB-S intervals perform
best for counties where the heterogeneity across county-specific mean radon
values are well described by the provided spatial and covariate information.
In comparison to the pivot intervals, FAB prediction intervals are frequently
narrower even though they theoretically maintain the same coverage level.
Empirical Bayesian intervals do not maintain the desired prediction error rate
within each county, and thus, in some cases, the intervals may be narrower
than necessary to achieve the desired coverage. This is particularly the case
for counties with a county-specific sample mean far from the state-wide sample
mean as frequentist coverage rates of empirical Bayesian intervals wane as
$|\bar{y}_{j}-\sum_{j}\bar{y}_{j}/J|$ grows (as shown in Figure 5(c)). Even
so, the various FAB conformal intervals perform similarly or better in terms
of interval width overall while maintaining the desired frequentist coverage
rate. Furthermore, the pivot and empirical Bayesian prediction intervals will
not maintain their respective notions of coverage if the parametric
assumptions underlying the methods do not hold true. As previously discussed,
the frequentist coverage guarantee of the conformal prediction intervals does
not rely on parametric assumptions.

In this data example, FAB prediction intervals outperform the commonly
utilized prediction intervals considered as judged by overall precision and
coverage guarantees. FAB prediction intervals attain the desired county-level
frequentist coverage guarantee, and are narrower than the other considered
prediction intervals for the majority of counties, particularly when
constructed with relevant auxiliary information. Sharing information in a
useful way can result in narrower prediction intervals than otherwise possible
while still achieving the desired coverage rate within a FAB prediction
interval. There are, however, some counties where the FAB intervals are
notably wider than the intervals constructed via parametric methods. In these
counties, we sacrifice precision for the non-parametric frequentist coverage
guarantee. We see clearly the trade-off between achieving narrow prediction
interval width and maintaining desired frequentist coverage.

##  7 Discussion

##  1 Appendix

One straightforward method for proving two measures are ECM is to show, for
each $i=1,...,n+1$, the sub-region of acceptance, $a_{i}$,
$a_{i}=\\{y_{n+1}\in\mathcal{Y}:c_{i}(y_{n+1})\leq c_{n+1}(y_{n+1})\\}$, is
the same for both measures. Formally,

######  Lemma 5.

For conformity measures $C,D$, if

$\\{y_{n+1}:c_{i}(y_{n+1})\leq
c_{n+1}(y_{n+1})\\}=\\{y_{n+1}:d_{i}(y_{n+1})\leq
d_{n+1}(y_{n+1})\\}\;\forall\;i=1,...,n+1$

then $C$ and $D$ are ECM.

###### Proof of Lemma 5.

$\displaystyle\\{y_{n+1}:{}$ $\displaystyle c_{i}\leq
c_{n+1}\\}=\\{y_{n+1}:d_{i}\leq d_{n+1}\\}\;\forall i=1,...,n+1$
$\displaystyle\Rightarrow{}$ $\displaystyle\\#\\{i:c_{i}\leq
c_{n+1}\\}=\\#\\{i:d_{i}\leq d_{n+1}\\}$ $\displaystyle\Rightarrow{}$
$\displaystyle p_{y,c}=p_{y,d}$

Thus each candidate $y$ will be treated the same under both conformity
measures $C,D$. ∎

###  1.1 Proofs

###### Proof of Lemma 1.

$f(y_{n+1}:y_{n+1}<=\gamma)$ increases stepwise from 1 to $n+1$ over the range
$(-\infty,\gamma]$. Likewise, $f(y_{n+1}:y_{n+1}>=\gamma)$ decreases stepwise
from $n+1$ to 1 over the range $[\gamma,\infty)$. Therefore, $f(y_{n+1}$ is a
step function over the domain $\mathcal{Y}$ with ordered range
$\\{1,2,...,n,n+1,n,...,2,1\\}$. ∎

###### Proof of Lemma 2.

The region in $\mathcal{Y}$ where $f(y_{n+1})\geq a$ for $a\in\\{1,...,n+1\\}$
is an interval. Therefore,

$A(\Ybf)=\\{y_{n+1}\in\mathcal{Y}:f(y_{n+1})/(n+1)>\alpha\\}$

is an interval. ∎

###### Proof of Thm 1.

WLOG, by Lemma 6, we consider the conformity between two values $y_{1},y_{2}$.
In this case,

$C_{B}\left(y_{1},y_{2}\right):=p\left(y_{2}|y_{1}\right)=\frac{\Gamma\left(\frac{a_{12}+1}{2}\right)}{\sqrt{a_{12}\pi}\Gamma\left(\frac{a_{12}}{2}\right)}\left(\left(\frac{b_{2|1}}{a_{12}}(1+\tau^{2}_{12})\right)^{-1/2}\left(1+\frac{1}{a_{12}}\frac{(y_{2}-\mu_{2|1})^{2}}{\frac{b_{2|1}}{a_{12}}(1+\tau^{2}_{12})}\right)^{-(a_{12}+1)/2}\right)$

where

$\displaystyle\tau^{2}_{12}=(1/\tau^{2}+1)^{-1}$
$\displaystyle\mu_{2|1}=(\mu/\tau^{2}+y_{1})\tau^{2}_{12}$ $\displaystyle
a_{12}=a+1$ $\displaystyle
b_{2|1}=b+y_{1}^{2}+\mu^{2}/\tau^{2}-(\mu/\tau^{2}+y_{1})^{2}\tau^{2}_{12}$

Now, suppose the conformal algorithm requires that we identify the region of
$y_{1}$ s.t. $C\left(y_{1},y_{2}\right)\leq C\left(y_{2},y_{1}\right)$. This
region can be shown to be the same as that obtained from
$C\left(\\{y_{1},y_{2}\\},y_{2}\right)\leq
C\left(\\{y_{1},y_{2}\\},y_{1}\right)$, which is the definition of equivalent
conformity measures:

$\displaystyle\hskip 19.91692ptC\left(y_{1},y_{2}\right)\leq
C\left(y_{2},y_{1}\right)$ $\displaystyle\Leftrightarrow
C\left(y_{1},y_{2}\right)/C\left(y_{2},y_{1}\right)\leq 1$
$\displaystyle\Leftrightarrow\frac{\left(b_{1|2}\right)\left(1+\frac{(y_{1}-\mu_{1|2})^{2}}{b_{1|2}(1+\tau^{2}_{12})}\right)^{a_{12}+1}}{\left(b_{2|1}\right)\left(1+\frac{(y_{2}-\mu_{2|1})^{2}}{b_{2|1}(1+\tau^{2}_{12})}\right)^{a_{12}+1}}\leq
1$
$\displaystyle\Leftrightarrow\frac{\left(b_{1|2}\right)\left(\frac{b_{1|2}(1+\tau^{2}_{12})+(y_{1}-\mu_{1|2})^{2}}{b_{1|2}}\right)^{a_{12}+1}}{\left(b_{2|1}\right)\left(\frac{b_{2|1}(1+\tau^{2}_{12})+(y_{2}-\mu_{2|1})^{2}}{b_{2|1}}\right)^{a_{12}+1}}\leq
1$
$\displaystyle\Leftrightarrow\left(\frac{b_{2|1}}{b_{1|2}}\right)^{a_{12}}\left(\frac{b_{1|2}(1+\tau^{2}_{12})+(y_{1}-\mu_{1|2})^{2}}{b_{2|1}(1+\tau^{2}_{12})+(y_{2}-\mu_{2|1})^{2}}\right)^{a_{12}+1}\leq
1$ (13)
$\displaystyle\Leftrightarrow\left(\frac{b_{2|1}}{b_{1|2}}\right)^{a_{12}}\leq
1$ $\displaystyle\Leftrightarrow b_{2|1}-b_{1|2}\leq 0{}$
$\displaystyle:=\left(\beta+y_{1}^{2}+\mu^{2}/\tau^{2}-(\mu/\tau^{2}+y_{1})^{2}\tau^{2}_{12}\right)$
$\displaystyle\hskip
28.45274pt-\left(\beta+y_{2}^{2}+\mu^{2}/\tau^{2}-(\mu/\tau^{2}+y_{2})^{2}\tau^{2}_{12}\right)\leq
0$
$\displaystyle\Leftrightarrow\left(y_{1}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}-\left(y_{2}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}\leq
0$ (14)

Since the second term in Equation 13 simplifies to 1:

$\displaystyle\text{numerator}\left(\frac{b_{1|2}(1+\tau^{2}_{12})+(y_{1}-\mu_{1|2})^{2}}{b_{2|1}(1+\tau^{2}_{12})+(y_{2}-\mu_{2|1})^{2}}\right)$
$\displaystyle:={}$ $\displaystyle
b_{1|2}(1+\tau^{2}_{12})+(y_{1}-\mu_{1|2})^{2}$ $\displaystyle:={}$
$\displaystyle\left(\beta+{y_{2}}^{2}+\mu^{2}/\tau^{2}-(\mu/\tau^{2}+{y_{2}})^{2}\tau_{12}^{2}\right)(1+\tau_{12}^{2})$
$\displaystyle\hskip
28.45274pt+({y_{1}}-\left(\mu/\tau^{2}+{y_{2}}\right)\tau_{12}^{2})^{2}$
$\displaystyle\hskip
28.45274pt+y_{2}\left[-2(\mu/\tau^{2})\tau_{12}^{2}(1+\tau_{12}^{2})+2(\mu/\tau^{2})(\tau_{12}^{2})^{2}\right]$
$\displaystyle\hskip 28.45274pt+y_{1}^{2}$ $\displaystyle\hskip
28.45274pt+y_{1}\left[-2(\mu/\tau^{2})\tau_{12}^{2}\right]$
$\displaystyle\hskip 28.45274pt+y_{1}y_{2}\left[-2\tau_{12}^{2}\right]$
$\displaystyle={}$ $\displaystyle D+y_{2}^{2}$ (15) $\displaystyle\hskip
28.45274pt+y_{2}\left[-2(\mu/\tau^{2})\tau_{12}^{2})\right]$
$\displaystyle\hskip 28.45274pt+y_{1}^{2}$ $\displaystyle\hskip
28.45274pt+y_{1}\left[-2(\mu/\tau^{2})\tau_{12}^{2}\right]$
$\displaystyle\hskip 28.45274pt+y_{1}y_{2}\left[-2\tau_{12}^{2}\right]$

for some function $D$ that does not depend on $y_{1}$ or $y_{2}$. Since the
coefficients on the $y_{1},y_{2}$ terms are the same, by symmetry, the
denominator of the second term in Eqn 13 will also be equal to Equation 15.

To complete the proof, we now show that the inequality
$C(\\{y_{1},y_{2}\\},y_{2})\leq C(\\{y_{1},y_{2}\\},y_{1})$ simplifies to
Equation 14. Notice that all parameters of the $t$ distribution that define
this chosen conformity measure will be the same regardless of which variable
we are obtaining the conformity score for. The conformity score is:

$C_{B}\left(\\{y_{1},y_{2}\\},X\right)=p(X|y_{1},y_{2})=dt_{\nu_{12}}\left(X|\mu_{12},{\color[rgb]{1,0,0}\variance_{12}}\right)$

where $\mu_{12}=\left(\mu/\tau^{2}+y_{1}+y_{2}\right)\tau^{2}_{12^{\prime}}$
and $\tau^{2}_{12^{\prime}}=\left(1/\tau^{2}+2\right)^{-1}$ for any
$X\in\mathcal{Y}$.

Then,

$\displaystyle\hskip 19.91692ptC_{B}\left(\\{y_{1},y_{2}\\},y_{2}\right)\leq
C_{B}\left(\\{y_{1},y_{2}\\},y_{1}\right)$
$\displaystyle\Leftrightarrow\left(y_{1}-\left(\mu/\tau^{2}+y_{1}+y_{2}\right)\tau^{2}_{12^{\prime}}\right)^{2}\leq\left(y_{2}-\left(\mu/\tau^{2}+y_{1}+y_{2}\right)\tau^{2}_{12^{\prime}}\right)^{2}$
$\displaystyle\Leftrightarrow\left(y_{1}-\mu/\tau^{2}\frac{\tau^{2}_{12^{\prime}}}{1-2\tau^{2}_{12^{\prime}}}\right)^{2}-\left(y_{2}-\mu/\tau^{2}\frac{\tau^{2}_{12^{\prime}}}{1-2\tau^{2}_{12^{\prime}}}\right)^{2}\leq
0$
$\displaystyle\Leftrightarrow\left(y_{1}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}-\left(y_{2}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}\leq
0$ (16)

since

$\displaystyle\frac{\tau^{2}_{12^{\prime}}}{1-2\tau^{2}_{12^{\prime}}}:={}$
$\displaystyle\frac{\left(1/\tau^{2}+2\right)^{-1}}{1-2\left(1/\tau^{2}+2\right)^{-1}}$
$\displaystyle={}$
$\displaystyle\frac{1}{\left(1/\tau^{2}+2\right)-2\left(1/\tau^{2}+2\right)\left(1/\tau^{2}+2\right)^{-1}}$
$\displaystyle={}$ $\displaystyle\frac{1}{1/\tau^{2}}$ $\displaystyle={}$
$\displaystyle\frac{1}{1/\tau^{2}+1-\left(\tau^{2}_{12}\right)^{-1}\tau^{2}_{12}}$
$\displaystyle=:{}$
$\displaystyle\frac{1}{\left(\tau^{2}_{12}\right)^{-1}-\left(\tau^{2}_{12}\right)^{-1}\tau^{2}_{12}}$
$\displaystyle={}$ $\displaystyle\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}.$

In summary, by Equations 14 and 16 we have that

$\displaystyle\hskip 19.91692ptC_{B}\left(y_{1},y_{2}\right)\leq
C_{B}\left(y_{2},y_{1}\right)$
$\displaystyle\Leftrightarrow\left(y_{1}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}-\left(y_{2}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}\leq
0$ $\displaystyle\Leftrightarrow C\left(\\{y_{1},y_{2}\\},y_{2}\right)\leq
C\left(\\{y_{1},y_{2}\\},y_{1}\right)$

This result holds for each $i\in\\{1,...,n+1\\}$ and so by Lemma 5,
$C_{B}\left(\\{Y_{1},...,Y_{n+1}\\}\backslash Y_{i},Y_{i}\right)$ and
$C_{B}\left(\\{Y_{1},...,Y_{n+1}\\},Y_{i}\right)$ are ECD under the Normal
working model. ∎

###### Proof of Lemma 3.

As shown in the proof of Theorem 1,

$\displaystyle\hskip 19.91692ptC\left(y_{1},y_{2}\right)\leq
C\left(y_{2},y_{1}\right)$ $\displaystyle\Leftrightarrow
C\left(\\{y_{1},y_{2}\\},y_{2}\right)\leq
C\left(\\{y_{1},y_{2}\\},y_{1}\right)$
$\displaystyle\Leftrightarrow\left(y_{1}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}-\left(y_{2}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}\leq
0$

Therefore, the inequality reduces to a quadratic function of the unknown
$y_{a}$. Label this function $h$:

$h(y_{1}):=\left(y_{1}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}-\left(y_{2}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}\leq
0$

By standard quadratic formula theory, we can draw a few conclusions:

  1. 1.

Notice that $h(y_{1})$ is in vertex form. In $h(y_{1})$, the coefficient on
the first term is positive (1), so the parabola will be upward facing.

  2. 2.

The discriminant obtained via the quadratic formula is
$4\left(y_{2}-\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)^{2}$.
For $P_{\thetabf}$ continuous, $P_{\theta}(y_{2}=c)=0$, so the discriminant
will be positive with probability 1.

$Pr(\text{discriminant
positive})=1-Pr\left(y_{2}=\mu/\tau^{2}\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}\right)=1$
(17)

Therefore, with probability 1, there will be two solutions to the equation
$h(y_{1})=0$.

By items (1) and (2), the solution to the inequality $h(y_{1})\leq 0$ will be
an interval with probability 1. ∎

###### Proof of Lemma 4.

Following notation used thus far in proofs, we aim to show (WLOG)

$y_{2}<\tilde{\mu}_{\theta}<f(y_{2})$

where $\tilde{\mu}_{\theta}:=(\mu/\tau^{2}+y_{2})\tau^{2}_{12}$ and
$f(y_{2}):=(\mu/\tau^{2})\tau^{2}_{12}(1-\tau^{2}_{12})^{-1}-y_{2}$. By the
quadratic formula, it is easy to show that the vertex in between the two
solutions ( $y_{2},f(y_{2})$) is
$v_{b}:=(\mu/\tau^{2})\tau^{2}_{12}(1-\tau^{2}_{12})^{-1}$.

Now, we are considering the case where $y_{2}<f(y_{2})$. Then, we know
$y_{2}<v_{b}$.

  1. 1.

First, WTS $y_{2}<\tilde{\mu}_{\theta}$.

$\displaystyle\hskip
19.91692pty_{2}<(\mu/\tau^{2})\tau^{2}_{12}(1-\tau^{2}_{12})^{-1}{}$ known
$\displaystyle\Leftrightarrow
y_{2}(1-\tau^{2}_{12})<(\mu/\tau^{2})\tau^{2}_{12}$
$\displaystyle\Leftrightarrow
y_{2}<(\mu/\tau^{2})\tau^{2}_{12}+y_{2}\tau^{2}_{12}$
$\displaystyle=:y_{2}<\tilde{\mu}_{\theta}$

  2. 2.

Now, WTS $f(y_{2})>\tilde{\mu}_{\theta}$.

$\displaystyle\hskip 19.91692pty_{2}<\tilde{\mu}_{\theta}{}$ shown above
$\displaystyle:=y_{2}<(\mu/\tau^{2})\tau^{2}_{12}+y_{2}\tau^{2}_{12}$
$\displaystyle\Leftrightarrow
y_{2}\left(1-\tau^{2}_{12}\right)<(\mu/\tau^{2})\tau^{2}_{12}$
$\displaystyle\Leftrightarrow
y_{2}<(\mu/\tau^{2})\frac{\tau^{2}_{12}}{\left(1-\tau^{2}_{12}\right)}$
$\displaystyle\Leftrightarrow
y_{2}<(\mu/\tau^{2})\left((1-\tau^{2}_{12})^{-1}-1\right)$
$\displaystyle\Leftrightarrow
y_{2}+\mu/\tau^{2}<(\mu/\tau^{2})(1-\tau^{2}_{12})^{-1}$
$\displaystyle\Leftrightarrow\left(y_{2}+\mu/\tau^{2}\right)\tau^{2}_{12}<(\mu/\tau^{2})(1-\tau^{2}_{12})^{-1}\tau^{2}_{12}$
$\displaystyle=:\tilde{\mu}_{\theta}<v_{b}$

since

$(1-\tau^{2}_{12})^{-1}-1=\frac{1-\left(1-\tau^{2}_{12}\right)}{1-\tau^{2}_{12}}=\frac{\tau^{2}_{12}}{1-\tau^{2}_{12}}$

and by quadratic equation theory, we know $v_{b}<f(y_{2})$! To summarize, by
rearranging terms, we can show that the following ordering of terms exists,
given that $y_{2}<f(y_{2})$:

$y_{2}<\tilde{\mu}_{\theta}<v_{b}<f(y_{2})$

∎

###### Proof of Lemma 2.

Thus, the region of $\mathcal{Y}$ where $\\#\\{c_{i}\leq c_{n+1}\\}\geq k$ for
integers $k=1,...,n+1$ is an interval.COMPLETE ∎

###### Proof of Thm 2.

By Lemmas 3 and 4, the region where $\mathcal{Y}$ where

$C\left(\\{Y_{1},...,Y_{n+1}\\}\backslash Y_{i},Y_{i}\right)\leq
C\left(\\{Y_{1},...,Y_{n+1}\\}\backslash Y_{n+1},Y_{n+1}\right)$

holds will be an interval that contains the shrinkage estimator of the sample
$Y_{1},...,Y_{n}$. As a result, the distribution of $\\#\\{c_{i}\leq
c_{n+1}\\}$ step-wise continuously increases from 1 to $n+1$ and then step-
wise continuously decreases from $n+1$ to $1$. Then, by Lemma 2, the
prediction region is an interval. ∎

###  1.2 Notation simplification to be used in proofs

We begin by introducing a lemma that will greatly simplify notation used in
proofs to follow.The intuition here is that, in our problem set-up, for two
random variables of interest, $Y_{j},Y_{k}$, conditioning on extra data
$\Ybf_{-j,-k}$ simply requires an update of the posterior of model parameters
$\thetabf$. Lemma 6 allows us to consider a simplified regime where we have a
single data point $y_{j}$ and a single candidate prediction $y_{n+1}$ All
other data can be considered to be part of the prior and is accounted for in
prior hyperparameters.

######  Lemma 6.

The posterior predictive distribution
$p(Y_{j}|\\{Y_{1},...,Y_{n+1}\\}\backslash Y_{j})$ under an IID hierarchical
model, i.e.:

$\displaystyle Y_{1},...,Y_{n+1}\distiid$ $\displaystyle P_{\thetabf}$ (18)
$\displaystyle\thetabf\sim{}$ $\displaystyle q(\thetabf).$

is equivalent to the posterior predictive distribution $p(Y_{j}|Y_{k})$ under
the prior

$\tilde{q}(\thetabf)=p_{q}(\thetabf|\Ybf_{-j,-k})=\frac{p(\Ybf_{-j,-k}|\thetabf)q(\thetabf)}{\int_{\Theta}p(\Ybf_{-j,-k}|\thetabf)q(\thetabf)d\thetabf}$

where $\Ybf_{-j,-k}:=\\{Y_{1},...,Y_{n+1}\\}\backslash\\{Y_{j},Y_{k}\\}$ for
any $k\in\\{1,...,n+1\\},\;k\neq j$.

###### Proof of Lemma 6.

We aim to show the posterior predictive distribution of $Y_{j}$ under Model 18
is equivalent to that under following model:

$\displaystyle Y_{j},Y_{k}\distiid{}$ $\displaystyle P_{\thetabf}$ (19)
$\displaystyle\tilde{q}(\thetabf)={}$ $\displaystyle
p_{q}(\thetabf|\Ybf_{-j,-k})=\frac{p(\Ybf_{-j,-k}|\thetabf)q(\thetabf)}{\int_{\Theta}p(\Ybf_{-j,-k}|\thetabf)q(\thetabf)d\thetabf}.$

where $p_{q}(\thetabf|\cdot)$ refers to the posterior distribution of
$\thetabf$ conditional on observed data $(\cdot)$ under the prior $q$. First,

$\displaystyle p_{\tilde{q}}(\thetabf|Y_{k})\propto{}$ $\displaystyle
p(Y_{k}|\thetabf)\tilde{q}(\thetabf)$ $\displaystyle:={}$ $\displaystyle
p(Y_{k}|\thetabf)p_{q}(\thetabf|\Ybf_{-j,-k})$ $\displaystyle\propto{}$
$\displaystyle p(Y_{k}|\thetabf)p(\Ybf_{-j,-k}|\thetabf)q(\thetabf)$
$\displaystyle\equiv{}$ $\displaystyle p(\Ybf_{-j}|\thetabf)q(\thetabf)$
$\displaystyle\propto{}$ $\displaystyle p_{q}(\thetabf|\Ybf_{-j})$

Then,

$\displaystyle p_{q}(Y_{j}|\Ybf_{-j})={}$
$\displaystyle\int_{\Theta}p_{q}(Y_{j},\thetabf|\Ybf_{-j})d\thetabf$
$\displaystyle={}$
$\displaystyle\int_{\Theta}p(Y_{j}|\thetabf)p_{q}(\thetabf|\Ybf_{-j})d\thetabf$
$\displaystyle\equiv{}$
$\displaystyle\int_{\Theta}p(Y_{j}|\thetabf)p_{\tilde{q}}(\thetabf|Y_{k})d\thetabf$
$\displaystyle={}$ $\displaystyle p_{\tilde{q}}(Y_{j}|Y_{k})$

∎

The results from Lemma 6 will be used in most proofs of the following results.


# Is one brick enough to break the wall of spoken dialogue state tracking?

###### Abstract

In Task-Oriented Dialogue (TOD) systems, correctly updating the system’s
understanding of the user’s needs (a.k.a dialogue state tracking) is key to a
smooth interaction. Traditionally, TOD systems perform this update in three
steps: transcription of the user’s utterance, semantic extraction of the key
concepts, and contextualization with the previously identified concepts. Such
cascade approaches suffer from cascading errors and separate optimization.
End-to-End approaches have been proved helpful up to the semantic extraction
step. This paper goes one step further paving the path towards completely
neural spoken dialogue state tracking by comparing three approaches: (1) a
state of the art cascade approach, (2) a locally E2E approach with rule-based
contextualization and (3) a completely neural approach. Our study highlights
that although they all outperform the recent DSTC11 best model, especially
with a filtering post-processing step, (1) remains the most accurate approach.
Indeed, both (2) and (3) have trouble propagating context as dialogues unfold
showing that context propagation in completely neural approaches is an open
challenge.

Index Terms—  spoken dialogue systems, context adaptation, end-to-end,
dialogue state tracking

##  1 Introduction

Digitization enables many tasks to be automated, nevertheless users sometimes
require assistance to perform a specific task such as making a reservation at
a restaurant or booking a hotel room. Task-Oriented Dialogue (TOD) systems are
designed to assist such users. A common approach to implement them is to break
the problem down to three iterative steps [1]: updating the system’s
understanding of the users’ needs, reasoning over a database and domain
knowledge to choose the next action and providing the user an answer. This
paper focuses on the first step.

Fig. 1: Spoken Dialogue State Tracking alternatives. Red characters indicate
potential cascading errors.

Traditionally the user’s needs update consists of three components,
respectively performing the transcription of the user’s utterance, semantic
extraction and contextualization of the extracted concepts [2]. Unfortunately,
this method presents the inconvenience of propagating errors of a component on
to the next one(s) (i.e. cascading errors) and of not optimizing all
components on the final objective (i.e. separate optimization) [3]. End-to-End
(E2E) approaches may address these issues by designing models in which the
gradient (i.e. error signal) can back-propagate from the output all the way to
the input [4].

On the one hand, with the advent of deep-learning and textual embeddings,
state of the art Dialogue State Tracking (DST) models work directly on
automatic transcriptions [5]. However, such approaches require careful,
dataset specific, mechanisms to catch and correct transcription errors
together with data augmentation to increase the model’s robustness to specific
upstream errors.

On the other hand, Spoken Language Understanding (SLU) working directly on the
speech signal, has successfully been applied to tasks which process utterances
individually such as voice command slot filling [6] and dialogue act
classification [7, 8]. Such systems often leverage transfer learning of
previously trained models which is challenging because it requires a trade-off
between learning new knowledge (e.g. domain’s vocabulary, domain’s ontology
structure) and keeping previous capabilities (e.g. transcription of open
vocabulary concepts) on a small amount of data [3].

In TOD such models must also adapt to the dialogue’s current context. While
E2E SLU models have been efficiently designed for independent single user
utterance processing [6, 9, 10], contextually dependant E2E SLU remains
unexplored to the best of our knowledge. Indeed, dialogue history integration
to guide the current turn’s prediction (e.g. better spelling of technical
vocabulary) has been already implemented [7, 8, 11, 12]. Yet, for the tasks
described in these studies, the context was not mandatory for solving the task
(e.g. the task does not need to process cross-turn references resolution).

Producing contextual semantic annotations is expensive because of the
cognitive load required to analyse the context and adapt the annotation. Such
datasets exist for chat based dialogue understanding [13, 14] but lack for
spoken dialogues, explaining the gap between E2E SLU and DST. To the best of
our knowledge, the vocalized version of Multi-Woz [15] is the only dataset
providing speech, transcriptions and contextual semantic annotations (i.e.
representing the user’s needs from the beginning of the dialog up to the
current turn).

This paper lies at the intersection of both directions. We focus on
contextually dependant semantic extraction, such as DST, in which the previous
dialogue turns are mandatory to correctly process the current one. In fact,
the spoken DST models presented in this paper output a summary of the user’s
needs since the beginning of the dialogue. State of the art DST systems use
cascade approaches which create a textual bottleneck both in terms of data and
model inference. E2E approaches do not require ground-truth transcriptions and
can be jointly optimized. We pave the path towards E2E spoken DST by comparing
(1) a state of the art cascade approach, (2) a locally E2E approach with rule-
based contextualization and (3) a completely neural approach.

Fig. 2: Three approaches for context propagation in spoken DST: (1) SOTA
cascade, (2) local E2E model with rule-based contextualization and (3)
completely neural model. (1) and (2) present cascading errors in red
characters. Hatched components are speech-related while solid ones are text-
related.

##  2 Method

###  2.1 Task-Oriented Dialogues

In TOD users require assistance from an agent to complete a task such as
making a reservation at a restaurant or booking a hotel room. More formally,
let us define a TOD as a sequence of $t$ dialogue turns
$U_{1},A_{2},\dots,A_{t-1},U_{t}$ where $A_{t-1}$ and $U_{t}$ respectively
correspond to textual agent’s turn $t-1$ and spoken user’s turn $t$. The goal
of DST is to keep up to date a condensed representation of the user’s needs.
In this paper, users needs are represented as Dialogue States (DS) and
correspond to a list of $n$ slot-value pairs linearized as
slot1=value1;...;slotn=valuen. At a given turn $t$, a TOD system is thus
inputted the previous context111For $t=1$, the context is empty. and the
current user turn from which it should output the updated user needs $DS_{t}$.

###  2.2 Context Propagation

As the dialogue unfolds, the user might refer to previously mentioned
entities. In order to design a contextually dependant SLU model, we need to
propagate the context of the previous turns $DS_{t-2}+A_{t-1}$ to inform the
prediction $\hat{DS}_{t}$ of the current user turn $U_{t}$. This paper
compares three alternatives of contextually dependant SLU models shown in Fig.
2. For each approach we train the model(s) over 10 epochs on a single 24Gb GPU
and use the last checkpoint at inference222Code will be made available upon
publication..

####  2.2.1 Cascade DST Approach

The cascade approach consists of an Automatic Speech Recognition (ASR) model
which transcribes the user’s turn $U_{t}$ and concatenates it to the previous’
turns context $DS_{t-2}+A_{t-1}$ as the input of a Natural Language
Understanding (NLU) model which then predicts the next Dialogue State
$DS_{t}$.

Both components are trained separately: WavLM [16], with two additional linear
layers outputting tokens’ probabilities, to transcribe the user turns (fine-
tuned with CTC loss) and T5 Encoder-Decoder [17] to output dialogue states.
Note that the NLU model is trained with user turns transcriptions of the ASR
model in order to be as close as possible to its inference regime.

####  2.2.2 Local E2E Approach

DS are updated through three operations: addition of a new slot-value pair,
modification and suppression of a previously mentioned slot-value pair. In
order to update $DS_{t-2}$ and obtain $DS_{t}$ with a rule based system we
need to encode these operations into local DS. While additions and
modifications of slots can remain the same in local DS, we mark suppressed
slots by assigning them the value <unk>.

Finally, for references to previously mentioned slots, as shown in Fig. 2, we
explicit the reference through the name of the referred slot’s value present
in $DS_{t-2}$.

This model relies on a Whisper [18] backbone model with a fine-tuned decoder
on the user turns. At inference time, the context $DS_{t-2}+A_{t-1}$ is added
before the decoder input tokens to condition the decoding. Note that erroneous
references and suppression (e.g. non-existing slots) are discarded by post-
processing the outputs.

####  2.2.3 Completely Neural Approach

The completely neural approach leverages both previous approaches by fusing
audio and semantic encoder’s outputs and feeding them to a semantic decoder.
The goal of this approach is to enable joint optimization of all components.
More formally, we have:

$\displaystyle h_{1}$ $\displaystyle=\text{Sem-Enc}(DS_{t-2}+A_{t-1})$
$\displaystyle h_{2}$ $\displaystyle=\text{Audio-Enc}(U_{t})$ $\displaystyle
h$ $\displaystyle=\text{Fusion}(h_{1}+\text{Conv}(h_{2}))$
$\displaystyle\hat{DS}_{t}$ $\displaystyle=\text{Sem-Dec}(h)$

WavLM [16] and T5’s encoder [17] respectively encode the current turn (audio)
and the context (text). Given that both models do not have the same processing
windows, two convolution layers (stride $3$, kernel size $9$) are added to
down-sample the audio encoder’s outputs. The fusion layer is a MLP which
enables the model to select and mix the information from both encoders.
Finally, a T5 decoder outputs $DS_{t}$ conditioned on the fusion of both
encoders’ outputs.

##  3 Results

###  3.1 Spoken Multi-Woz dataset

Multi-Woz is a human-human chat-based English Task-Oriented Dialogue (TOD)
dataset commonly used for training and evaluating dialogue systems [14]. A
spoken version of Multi-Woz with vocalized user turns was published in the
context of the Speech Aware Dialogue Systems track of the 11 edition of the
Dialogue System Technology Challenge333https://dstc11.dstc.community/ (DSTC11)
[15]. The user utterances in the training set are available as synthetic
speech, whereas the dev and test sets (Dev$|$Test) include both synthetic and
human speech versions (TTS$|$Human). The dataset contains close to 10,000
dialogues with a 80/10/10 train-dev-test split and an average of 13.3 turns
per dialogue. Among the pre-defined slots, we can distinguish 3 groups:
categorical slots with a closed set of values ($\sim$60%), non-categorical
slots with an open set of values ($\sim$30%) and time slots ($\sim$10%). Note
that, in order to reduce the value overlap across sets, non-categorical slots
were replaced and time slots offset in the Dev and Test sets.

###  3.2 Evaluation

We evaluate all approaches with a turn-level exact match metric known as
Joint-Goal Accuracy (JGA$\uparrow$) [19]. This metric requires to post-process
the coma separated slot-value output format to convert it into a valid
dictionary which does not take into account the order of the slot-value pairs.
Given that generative models, such as T5, are prone to hallucinations, a
filtering step discards all slots which are not part of the predefined slots.
We present here the results of all three approaches in two scenarios: with
ground truth context $DS_{t-2}$, in Table 1 and with the previous prediction
$\hat{DS}_{t-2}$ in Fig. 4.

| Dev | Test
---|---|---
Text | 81.2 | 80.3
w/o filtering | 64.8 | 63.8
| TTS | Human | TTS | Human
DSTC11 baseline [15] | n/a | n/a
w/o filtering | 38.4 | 31.8 | n/a
DSTC11 best [5] | n/a | n/a
w/o filtering | 47.2 | 43.2 | 44.0 | 39.5
(1) Cascade | 75.2 | 71.9 | 75.4 | 71.8
w/o filtering | 53.5 | 48.6 | 53.2 | 47.8
(2) Local E2E | 60.5 | 61.8 | 62.3 | 62.6
w/o filtering | 43.6 | 43.8 | 42.8 | 42.7
(3) E2E | 56.8 | 54.2 | 57.1 | 53.9
w/o filtering | 14.9 | 12.6 | 14.3 | 12.3

Table 1: JGA of spoken DST with ground-truth previous state $DS_{t-2}$. Text
line shows upper-bound performance. Note that [15] and [5] use dialogue
history as input to the DST model.

We find that the cascade approach remains a tough competitor only 6 points
behind the text oracle model. Although both the local and global E2E
approaches achieve higher accuracies than the recent DSTC11 best model,
especially with the filtering post-processing step, they are not competitive
with a carefully designed cascade approach. It is noteworthy that filtering
out the undefined slots has a significant impact for all approaches. The
completely E2E approach seems particularly prone to hallucinations which
indicates that it has more trouble selecting the relevant information from the
audio encoder’s hidden states.

In order to get a more precise understanding of the differences between these
approaches, we further evaluate each slot’s F1-measure and present each slot
groups average F1-measure on the Test-Human set in Fig. 3. Categorical slots
present no difficulty while non-categorical and time slots are more
challenging. Interestingly the backbone model seems to play a role in the
output format given that the local E2E and completely E2E models’ F1-measures
are reversed for the Time slot group which requires careful formatting.

Fig. 3: Test-Human slot group average F1 measure. Fig. 4: Test-Human turn
accuracy with and without ground-truth previous state for each approach. Note
that there are fewer and fewer dialogues as the number of turns increases.

In a more realistic scenario where we base our next prediction on the previous
one444Note that $A_{t-1}$ and $U_{t}$ remain unchanged which might lead to
some incoherences. we observe that, as illustrated in Fig. 4, both the cascade
and the completely E2E approaches perform as good with their predictions as
with the ground-truth whereas the local E2E approach collapses. In fact it
achieves a higher accuracy on the first turn and collapses on the next ones
which highlights that using decoder prefixing to propagate the dialogue’s
context might not be the best method. Also note that the global E2E approach
has more trouble handling long dialogues than its alternatives.

##  4 Conclusion

In order to pave the path towards E2E spoken DST, we compare a state of the
art cascade approach with local and global E2E approaches. Our study
highlights that although they all outperform the recent DSTC11 winner model,
especially with a filtering post-processing step, the cascade approach remains
the most accurate approach and context propagation in completely neural
approaches an open challenge. Local E2E and global E2E approaches behave quite
differently: the former is very accurate on the first dialogue turn but
collapses when contextualizing its predictions and the latter hallucinates
much more and has trouble with long dialogues.

Our results remain to be confirmed on other contextually dependant SLU
datasets to come. Given the chat-based origin of Multi-Woz, the turns are
assumed to be perfectly separable which is not the case in general. A more
fine-grained evaluation to assess which errors are low-impact errors (e.g.
rectified with the help of a database, with no impact on the dialogue
trajectory) and improving the post processing with a more careful filtering is
left as future work.

## References

  * [1] Gokhan Tur and Renato De Mori,  SLU in Commercial and Research Spoken Dialogue Systems,  Wiley Telecom, 2011.
  * [2] Jason D. Williams, Antoine Raux, and Matthew Henderson,  “The Dialog State Tracking Challenge Series: A Review,”  Dialogue & Discourse, 2016.
  * [3] Salima Mdhaffar, Valentin Pelloin, Antoine Caubrière, Gaëlle Laperrière, Sahar Ghannay, Bassam Jabaian, Nathalie Camelin, and Yannick Estève,  “Impact Analysis of the Use of Speech and Language Models Pretrained by Self-Supersivion for Spoken Language Understanding,”  in Conference on Language Resources and Evaluation (LREC), 2022\.
  * [4] Dmitriy Serdyuk, Yongqiang Wang, Christian Fuegen, Anuj Kumar, Baiyang Liu, and Yoshua Bengio,  “Towards end-to-end spoken language understanding,”  in ICASSP, 2018.
  * [5] Léo Jacqmin, Lucas Druart, Valentin Vielzeuf, Lina Maria Rojas-Barahona, Yannick Estève, and Benoît Favre,  “OLISIA: a Cascade System for Spoken Dialogue State Tracking,”  in Special Interest Group on Discourse and Dialogue (SIGDIAL), 2023\.
  * [6] Siddhant Arora, Hayato Futami, Shih-Lun Wu, Jessica Huynh, Yifan Peng, Yosuke Kashiwagi, Emiru Tsunoo, Brian Yan, and Shinji Watanabe,  “A study on the integration of pipeline and e2e slu systems for spoken semantic parsing toward stop quality challenge,”  in ICASSP, 2023.
  * [7] Vishal Sunder, Samuel Thomas, Hong-Kwang J. Kuo, Jatin Ganhotra, Brian Kingsbury, and Eric Fosler-Lussier,  “Towards end-to-end integration of dialog history for improved spoken language understanding,”  in ICASSP, 2022.
  * [8] Jatin Ganhotra, Samuel Thomas, Hong-Kwang J. Kuo, Sachindra Joshi, George Saon, Zoltán Tüske, and Brian Kingsbury,  “Integrating Dialog History into End-to-End Spoken Language Understanding Systems,”  in Interspeech, 2021.
  * [9] Guangzhi Sun, Chao Zhang, and Philip C. Woodland,  “End-to-end spoken language understanding with tree-constrained pointer generator,”  in ICASSP, 2023.
  * [10] Valentin Pelloin, Nathalie Camelin, Antoine Laurent, Renato De Mori, Antoine Caubrière, Yannick Estève, and Sylvain Meignier,  “End2end acoustic to semantic transduction,”  in ICASSP, 2021.
  * [11] Natalia Tomashenko, Christian Raymond, Antoine Caubrière, Renato De Mori, and Yannick Estève,  “Dialogue history integration into end-to-end signal-to-concept spoken language understanding systems,”  in ICASSP, 2020.
  * [12] Vishal Sunder, Eric Fosler-Lussier, Samuel Thomas, Hong-Kwang J Kuo, and Brian Kingsbury,  “ConvKT: Conversation-Level Knowledge Transfer for Context Aware End-to-End Spoken Language Understanding,”  in Interspeech, 2023.
  * [13] Xiangkun Hu, Junqi Dai, Hang Yan, Yi Zhang, Qipeng Guo, Xipeng Qiu, and Zheng Zhang,  “Dialogue meaning representation for task-oriented dialogue systems,”  in Conference on Empirical Methods in Natural Language Processing (EMNLP), 2022.
  * [14] Paweł Budzianowski, Tsung-Hsien Wen, Bo-Hsiang Tseng, Iñigo Casanueva, Stefan Ultes, Osman Ramadan, and Milica Gasic,  “Multiwoz - a large-scale multi-domain wizard-of-oz dataset for task-oriented dialogue modelling,”  in Conference on Empirical Methods in Natural Language Processing (EMNLP), 2018.
  * [15] Hagen Soltau, Izhak Shafran, Mingqiu Wang, Abhinav Rastogi, Jeffrey Zhao, Ye Jia, Wei Han, Yuan Cao, and Aramys Miranda,  “Speech aware dialog system technology challenge (dstc11),”  in Special Interest Group on Discourse and Dialogue (SIGDIAL), 2023\.
  * [16] Sanyuan Chen, Chengyi Wang, Zhengyang Chen, Yu Wu, Shujie Liu, Zhuo Chen, Jinyu Li, Naoyuki Kanda, Takuya Yoshioka, Xiong Xiao, Jian Wu, Long Zhou, Shuo Ren, Yanmin Qian, Yao Qian, Micheal Zeng, and Furu Wei,  “Wavlm: Large-scale self-supervised pre-training for full stack speech processing,”  IEEE Journal of Selected Topics in Signal Processing (JSTSP), 2021\.
  * [17] Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu,  “Exploring the limits of transfer learning with a unified text-to-text transformer,”  Journal Machine Learning Research (JMLR), 2020.
  * [18] Alec Radford, Jong Wook Kim, Tao Xu, Greg Brockman, Christine McLeavey, and Ilya Sutskever,  “Robust speech recognition via large-scale weak supervision,”  ArXiv, vol. abs/2212.04356, 2022.
  * [19] Victor Zhong, Caiming Xiong, and Richard Socher,  “Global-locally self-attentive encoder for dialogue state tracking,”  in Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics, 2018.