# Conformal Multi-Target Hyperrectangles

Max Sampson<EMAIL_ADDRESS>
Department of Statistics & Actuarial Science
University of Iowa
Iowa City, Iowa 52242-1409, USA Kung-Sik Chan<EMAIL_ADDRESS>
Department of Statistics & Actuarial Science
University of Iowa
Iowa City, Iowa 52242-1409, USA

###### Abstract

We propose conformal hyperrectangular prediction regions for multi-target
regression. We propose split conformal prediction algorithms for both point
and quantile regression to form hyperrectangular prediction regions, which
allow for easy marginal interpretation and do not require covariance
estimation. In practice, it is preferable that a prediction region is
balanced, that is, having identical marginal prediction coverage, since
prediction accuracy is generally equally important across components of the
response vector. The proposed algorithms possess two desirable properties,
namely, tight asymptotic overall nominal coverage as well as asymptotic
balance, that is, identical asymptotic marginal coverage, under mild
conditions. We then compare our methods to some existing methods on both
simulated and real data sets. Our simulation results and real data analysis
show that our methods outperform existing methods while achieving the desired
nominal coverage and good balance between dimensions.

Keywords: Balance, Conformal Prediction, Multi-Output Prediction, Multi-target
regression, Uncertainty Quantification,

Max Sampson & Kung-Sik Chan

##  1 Introduction

We consider the problem of forming semi-parametric prediction regions with a
$p$-dimensional response for $p\geq 2$. Assume $n$ i.i.d. (more generally,
exchangeable) copies of
$(({\boldsymbol{Y}}_{1},{\boldsymbol{X}}_{1}),({\boldsymbol{Y}}_{2},{\boldsymbol{X}}_{2})\ldots,({\boldsymbol{Y}}_{n},{\boldsymbol{X}}_{n}))$,
are observed. Then, we want a set,
${\boldsymbol{C}}({\boldsymbol{x}})={\boldsymbol{C}}_{n}(({\boldsymbol{Y}}_{1},{\boldsymbol{X}}_{1}),$
$({\boldsymbol{Y}}_{2},{\boldsymbol{X}}_{2}),\ldots,$
$({\boldsymbol{Y}}_{n},{\boldsymbol{X}}_{n}),{\boldsymbol{x}})$ such that for
a new pair $({\boldsymbol{Y}}_{n+1},{\boldsymbol{X}}_{n+1})$,

$\mathbb{P}({\boldsymbol{Y}}_{n+1}\in{\boldsymbol{C}}({\boldsymbol{X}}_{n+1}))\geq
1-\alpha$ (1)

If a suitable non-conformity score is used, this coverage can be guaranteed in
finite samples using conformal prediction (Vovk et al., 2005; Lei and
Wasserman, 2013). As with other conformal prediction methods, we aim to make
no distributional assumptions with our prediction sets. Because of that, we
focus on a hyperrectangular shape for the prediction sets. This is a similar
idea to hyperrectangular confidence intervals and balanced bootstrap intervals
from Robertson et al. (2021); Beran (1988); Romano and Wolf (2010). One
advantage that hyperrectangular regions have compared to ellipsoidal regions
is that they are easier to visualize in more than two dimensions. Along with
being easier to visualize, a marginal interval can also easily be found. For
example, in a two dimensional response case, the marginal interval for the
first response is the same regardless of the second response’s value. This
would not be the case with an ellipsoidal prediction region. Though, this
marginal interval may have a coverage rate higher than $1-\alpha$. Another
advantage is that no estimate of a covariance structure between dimensions
needs to be made for rectangular regions.

Borrowing from conformal prediction for multivariate functional data, one
option would be to use a non-conformity score of
$V=\max\limits_{i=1,\ldots,q}|Y_{i}-\hat{Y}_{i}|$ (Diquigiovanni et al.,
2022). This guarantees the coverage in (1), but fails to take into account the
fact that one (or more) of the $p$-dimensions may have a higher variability
than the others, causing overcoverage and poor balance. There are also methods
proposed in Kuleshov et al. (2018); Johnstone and Cox (2021) that use distance
functions, like the L2 norm or Mahalanobis distance, for their non-conformity
scores. While in theory, this is a valid approach, it requires a
$p-$dimensional grid search which limits the practical use to a few
dimensions. The Mahalanobis distance also requires a covariance matrix to be
estimated. Other methods use a Bonferroni or similar correction for
independently modeled responses (Neeven and Smirnov, 2018; Messoudi et al.,
2020). This can cause overcoverage and wider than necessary regions. We
outline a new non-conformity score that attempts to take the differences in
variability into account, while guaranteeing the coverage in (1). Moreover,
the new non-conformity score is shown to achieve an asymptotically balanced
conformal prediction hyperrectangle under mild conditions.

For multiple comparisons, balance is defined as having the same type 1 error
rate for each comparison (Beran, 1988; Tukey, 1953; Scheffe, 1953; Roy and
Bose, 1953). We extend this idea to hyperrectangular prediction regions for a
$p$-dimensional response. Balance for hyperrectangular prediction regions
means that the marginal prediction intervals enjoy identical coverage rates.
For multiple comparisons, balance has a natural interpretation, fairness. For
example, the importance of one follow-up test from an ANOVA model is not more
important than another. One wouldn’t conduct follow-up tests where the first
comparison has a marginal type 1 error of $\alpha$ and all others have a
marginal type 1 error of $0$. The same is true for a $p$-dimensional
prediction region, coverage in one dimension is not less important than
coverage in another.

We can imagine a scenario where we construct a joint prediction region such
that the overall coverage rate is $1-\alpha$. The first $k-1$ dimensions have
marginal prediction intervals of $\mathbb{R}$, and the last has a marginal
prediction set that has a coverage rate equal to $1-\alpha$. This would be a
practically useless prediction region. But, we can imagine a more realistic
scenario where we disregard balance and have $k-1$ prediction sets that have
marginal coverage slightly less than $1$, and one prediction set that has
marginal coverage slightly larger than $1-\alpha$. Clearly, this would be
undesirable from a practical standpoint (Bruder and Wolf, 2017).

We consider point regression, including mean regression, in Section 2, and
detail our proposed algorithm in Algorithm 1. A drawback of the point
regression approach is that the conformal hyperrectangles with a fixed
coverage rate are of identical (hyper)volume for any set of covariates. For
conditionally heteroscedastic data, we would want the volume of the
hyperrectangle to vary based on the observed covariates. We consider the use
of quantile regression for handling conditional heteroscedascity in Section 3,
and detail the algorithm in Algorithm 2. For ease of exposition, we assume
that $p$ = 3 in both cases, but it is easy to see how it could be extended to
dimensions higher than 3. We discuss the validity and the balance of the
proposed methods in Sections 4 and 5. We illustrate the empirical performance
of the proposed methods via simulations in Section 6, and a real application
in Section 7. We conclude in Section 8.

##  2 Conditional Point Regression

We now describe conformal hyperrectangular regression (CHR). As with split
conformal prediction, we begin by splitting our data into sets used for
training and calibration (Lei et al., 2018; Vovk et al., 2005). However,
unlike classical split conformal prediction, we split using two calibration
sets instead of one. The training set is indexed by $\mathcal{I}_{tr}$, the
first calibration set by $\mathcal{I}_{cal1}$, and the second calibration set
by $\mathcal{I}_{cal2}$. Given any point regression algorithm,
${\boldsymbol{f}}$, we fit a $p-dimensional$ function $\hat{{\boldsymbol{f}}}$
on the training set:

$\\{\hat{{\boldsymbol{f}}}\\}\leftarrow{\boldsymbol{f}}(\\{({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i}):i\in\mathcal{I}_{tr}\\}).$

In the next step, for each of the $p$-dimensions of the response, we compute
non-conformity scores. These scores will depend on the non-conformity measure
one wishes to use; we’ll use the absolute difference to demonstrate how to
form CHR prediction regions. The scores are computed on the first calibration
set:

$V_{i,j}=|\hat{f}_{j}({\boldsymbol{X}}_{i,j})-Y_{i,j}|,\text{ }\forall
i\in\mathcal{I}_{cal1}.$

Define the vectors
${\boldsymbol{V_{l}}}=(V_{1l},V_{2l},\ldots,V_{n{cal_{1}}l}),\text{
}l=1,\ldots,p$. Then, we form initial intervals based on the second
calibration set for each $j=1,\ldots,p$. The $p$th dimension’s intervals take
the form

$C_{p}({\boldsymbol{X}}_{k})=[\hat{f}_{p}({\boldsymbol{X}}_{k})-Q_{1-\alpha}({\boldsymbol{V}}_{p};Z_{cal1}),\hat{f}_{p}({\boldsymbol{X}}_{k})+Q_{1-\alpha}({\boldsymbol{V}}_{p};Z_{cal1})],\text{
}\forall k\in\mathcal{I}_{cal2}$ (2)

where,

$Q_{\delta}({\boldsymbol{V}};Z_{cal}):=(\delta)(1+\frac{1}{|Z_{cal}|})-\text{th
empirical quantile of }\\{V_{i}\\}$

The other $p-1$ prediction intervals take the same form with their
corresponding non-conformity scores. Define the lower-bound of the prediction
interval in (2) as $\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{k})$ and the
upper-bound as $\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{k})$.

At this stage, in an effort to have some balance between the dimensions, we
record the length of the $p$ sides. Let $A$ denote the length of the first
interval, $B$ denote the length of the second interval, and $C$ the length of
the third interval. If we use the absolute non-conformity score or the signed
error non-conformity score, these lengths will be constant for all of the
$ncal_{1}$ intervals. How to deal with non-conformity scores that give
different sized prediction intervals for difference covariates is covered in
Section 3.

Now, using the intervals found in (2), we compute a non-conformity score to
make sure our rectangular prediction regions will have the coverage described
in (1). For each of the $p$-dimensions of the response, compute the following
non-conformity scores:

$E_{k,j}=\max\\{\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{k})-Y_{k,j},Y_{k,j}-\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{k})\\},\
\forall k\in\mathcal{I}_{cal2}\text{ and }j=1,\ldots,p.$

We could just use these as the final non-conformity scores, but to help with
balance and reducing the size of the overall prediction regions, we now
convert $E_{k,j}$ to the length of the first dimension using the side lengths
we found earlier. Denote these converted scores as $A_{ik}$. For clarity, some
examples of conversions would be $A_{11}=E_{11}$,
$A_{12}=E_{12}\times\frac{A}{B}$, $A_{13}=E_{13}\times\frac{A}{C}$, etc. Now
define $W_{k}=\max\limits_{j=1,\ldots,p}A_{k,j}$.

Our final steps are to find the adjustments for each of the dimensions, and
construct the final prediction region. Define $Adj_{1}$ to be the
$(1-\alpha)(1+\frac{1}{|Z_{cal2}|})-$th empirical quantile of $\\{W_{k}\\}$.
This is the final adjustment for the first (more generally, the reference)
dimension. We can find the adjustments for the other dimensions by undoing the
conversion to the reference side. For clarity, some examples of this with the
first dimension as the reference would be, $Adj_{2}=Adj_{1}\times\frac{B}{A}$,
$Adj_{3}=Adj_{1}\times\frac{C}{A}$, etc. Finally, given new input data,
${\boldsymbol{X}}_{n+1}$, the prediction region for ${\boldsymbol{Y}}_{n+1}$
is constructed as

$\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1})=[\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{X}}_{n+1})-{\boldsymbol{Adj}},\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{X}}_{n+1})+{\boldsymbol{Adj}})].$

For reference, the procedure is summarized in Algorithm 1.

Algorithm 1 P-dimensional Point Conformal Regression

Input: level $\alpha$, data =
$\mathcal{Z}=({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})_{i\in\mathcal{I}}$,
score function $V$, test point ${\boldsymbol{x}}$, and regression algorithm
${\boldsymbol{f}}({\boldsymbol{X}};\mathcal{D})$ for the point estimate using
$\mathcal{D}$ as data
Procedure:

1:Split $\mathcal{Z}$ into a training fold
$\mathcal{Z}_{tr}\overset{\Delta}{=}({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})_{i\in\mathcal{I}_{tr}}$
and two calibration folds
$\mathcal{Z}_{cal1}\overset{\Delta}{=}({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})_{i\in\mathcal{I}_{cal1}}$
and
$\mathcal{Z}_{cal2}\overset{\Delta}{=}({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})_{i\in\mathcal{I}_{cal2}}$

2:Train the model ${\boldsymbol{\hat{f}}}({\boldsymbol{x}};Z_{tr})$ using the
training set

3:For each $i\in\mathcal{I}_{cal1}$ and each $j\in 1,\ldots,p$, compute the
score $V_{i,j}$

4:Find $p\times|\mathcal{Z}_{cal2}|$ $1-\alpha$ conformal prediction intervals
using the second calibration set’s covariates. Define these intervals as
$(\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{k}),\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{k}))$
for each $k\in\mathcal{I}_{cal2}$

5:Find the $p-1$ ratios of these intervals and define them as
$\frac{A}{B}=\frac{\hat{q}_{11,\text{high}}({\boldsymbol{X}}_{1})-\hat{q}_{11,\text{low}}({\boldsymbol{X}}_{1})}{\hat{q}_{12,\text{high}}({\boldsymbol{X}}_{1})-\hat{q}_{12,\text{low}}({\boldsymbol{X}}_{1})}$,
$\frac{A}{C}=\frac{\hat{q}_{11,\text{high}}({\boldsymbol{X}}_{1})-\hat{q}_{11,\text{low}}({\boldsymbol{X}}_{1})}{\hat{q}_{13,\text{high}}({\boldsymbol{X}}_{1})-\hat{q}_{13,\text{low}}({\boldsymbol{X}}_{1})}$,
…

6:For each $k\in\mathcal{I}_{cal2}$ and each $j\in 1,\ldots,p$ compute the
score
$E_{k,j}=\max\\{\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{k})-Y_{k,j},Y_{k,j}-\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{k})\\}$

7:For each $k\in\mathcal{I}_{cal2}$ convert the $E_{k,j}$ to the first
dimension using the ratios found in step 5. For example, $A_{11}=E_{11}$ and
$A_{12}=E_{12}*\frac{A}{B}$

8:Define $W_{k}=\max_{j}A_{k,j}$

9:Find $Adj_{1}=(1-\alpha)(1+\frac{1}{|Z_{cal2}|})-\text{th empirical quantile
of }\\{W_{k}\\}$

10:Using the ratios from step 5, find $Adj_{2}=Adj_{1}*\frac{B}{A}$,
$Adj_{3}=Adj_{1}*\frac{C}{A}$, …

Output:
$\hat{{\boldsymbol{C}}}({\boldsymbol{x}})=(\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{x}})-{\boldsymbol{Adj}},\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{x}})+{\boldsymbol{Adj}})$,
where $\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{x}})$ and
$\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{x}})$ are the p-dimensional
prediction sets for ${\boldsymbol{x}}$ defined in step 4 and
${\boldsymbol{Adj}}$ is the p-dimensional vector of adjustments found in steps
9 and 10

###  2.1 Initial Hyperrectangle Forming Methods

Here, we describe some existing methods for forming conformal prediction
intervals with conditional mean estimates (more generally, point estimates).
One of the most common non-conformity scores for a point estimate is the
absolute difference, $V_{i}=|Y_{i}-\hat{f}({\boldsymbol{X}}_{i})|$
(Papadopoulos et al., 2002; Vovk et al., 2005). The prediction interval in the
case of the absolute difference the interval is

$C({\boldsymbol{X}}_{n+1})=[\hat{f}({\boldsymbol{X}}_{n+1})-Q_{1-\alpha}({\boldsymbol{V}};Z_{cal1}),\hat{f}({\boldsymbol{X}}_{n+1})+Q_{1-\alpha}({\boldsymbol{V}};Z_{cal1})]$

where

$Q_{1-\alpha}({\boldsymbol{V}};Z_{cal1}):=(1-\alpha)(1+\frac{1}{|Z_{cal1}|})-\text{th
empirical quantile of }\\{V_{i}\\}.$

and $|Z_{cal1}|$ is the size of the calibration set. We refer to this method
as Mean Absolute method in Section 6.

A second option that allows one to control the tail error rates is signed-
conformal regression. Define the signed error non-conformity scores,
$R_{i}=Y_{i}-\hat{f}({\boldsymbol{X}}_{i})$ and $V_{i}=-1\times R_{i}$. Then,
the signed error conformal prediction region (SECPR) is given by

$C({\boldsymbol{X}}_{n+1})=[\hat{f}({\boldsymbol{X}}_{n+1})-Q_{1-\alpha_{1}}({\boldsymbol{V}};Z_{cal1}),\hat{f}({\boldsymbol{X}}_{n+1})+Q_{1-\alpha_{2}}({\boldsymbol{R}};Z_{cal1})]$

where $\alpha=\alpha_{1}+\alpha_{2}.$ For example, if one wanted a conformal
prediction interval with equal tailed errors, they could take
$\alpha_{1}=\alpha_{2}=\alpha/2$ (Linusson et al., 2014; Romano et al., 2019).
We refer to this method as Mean Signed in Section 6.

##  3 Quantile Regression

We now describe an extension of conformal hyperrectangle for regression to
non-constant prediction interval widths, conformal quantile hyperrectangular
regression. As with other split conformal prediction methods, we begin by
splitting our data into sets used for training and calibration. The training
set is indexed by $\mathcal{I}_{tr}$ and the calibration set by
$\mathcal{I}_{cal}$. Given any quantile regression algorithm,
$\mathcal{{\boldsymbol{A}}}$, we fit 2 $\times$ $p$ conditional quantile
functions $\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{x}})$ and
$\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{x}})$ on the training set.
Any quantiles will work as long as the quantile for
$\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{x}})$ is smaller than that
of $\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{x}})$. In practice,
splitting $\alpha$ in two and giving the tails equal probability tends to work
well.

$\\{\hat{{\boldsymbol{q}}}_{\text{low}},\hat{{\boldsymbol{q}}}_{\text{high}}\\}\leftarrow\mathcal{{\boldsymbol{A}}}(\\{({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i}):i\in\mathcal{I}_{tr}\\})$

We now form initial intervals for the calibration set based on the trained
model. The intervals take the form

${\boldsymbol{C}}({\boldsymbol{X}}_{k})=[\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{X}}_{k}),\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{X}}_{k})],\text{
}\forall k\in\mathcal{I}_{cal}.$ (3)

To keep some balance, we record the length of the $p\times|\mathcal{I}_{cal}|$
sides. Let $A_{i}$ ($B_{i}$, $C_{i}$) denote the length of the first (second,
third) interval for the $i$-th data pair in the calibration set. Note that the
$A$’s are generally distinct real numbers and so are the $B$’s and $C$’s. That
is, the lengths of the sides for different covariates will generally not be
the same, unlike in Section 2.

Now, using the intervals found in (3), we compute a non-conformity score to
make sure our rectangular prediction regions will have the coverage described
in (1). For each of the $p$-dimensions of the response, compute the following
non-conformity scores:

$E_{k,j}=\max\\{\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{k})-Y_{k,j},Y_{k,j}-\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{k})\\},\
\forall k\in\mathcal{I}_{cal}\text{ and }j=1,\ldots,p.$

We now convert $E_{k,j}$ to the length of the first dimension using the side
lengths we found earlier. Denote these converted scores as $W_{k,j}$. For
clarity, some examples of conversions would be $W_{11}=E_{11}$,
$W_{12}=E_{12}\times\frac{A_{1}}{B_{1}}$,
$W_{21}=E_{21},W_{22}=E_{22}\times\frac{A_{2}}{B_{2}}$, etc. Now, define
$W_{k}=\max\limits_{j=1,\ldots,p}W_{k,j}$.

Our final steps are to find the adjustments for each of the dimensions, and
construct the final prediction region. First, we form the intervals for the
data with unseen responses,

${\boldsymbol{C}}({\boldsymbol{X}}_{n+1})=[\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{X}}_{n+1}),\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{X}}_{n+1})].$
(4)

Now, define $Adj_{1}$ to be the $(1-\alpha)(1+\frac{1}{|Z_{cal}|})-\text{th
empirical quantile of }\\{W_{k}\\}$. This is the final adjustment for the
first dimension. We can find the adjustments for the other dimensions by
undoing the conversion to the first side. To do this, we’ll need to find the
lengths of the 1st, 2nd, 3rd, etc. sides in (4). Denote these lengths as
$A_{n+1}$, $B_{n+1}$, $C_{n+1}$, etc. Then, the final adjustments for the
other sides would be $Adj_{2}=Adj_{1}\times\frac{B_{n+1}}{A_{n+1}}$,
$Adj_{3}=Adj_{1}\times\frac{C_{n+1}}{A_{n+1}}$, etc. Finally, given new input
data, ${\boldsymbol{X}}_{n+1}$, the prediction region for
${\boldsymbol{Y}}_{n+1}$ is constructed as

$\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1})=[\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{X}}_{n+1})-{\boldsymbol{Adj}},\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{X}}_{n+1})+{\boldsymbol{Adj}})].$

For reference, the procedure is summarized in Algorithm 2.

Algorithm 2 P-dimensional Quantile Conformal Regression

Input: level $\alpha$, data =
$\mathcal{Z}=({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})_{i\in\mathcal{I}}$,
score function $V$, test point ${\boldsymbol{x}}$, and regression algorithm
$\mathcal{{\boldsymbol{A}}}({\boldsymbol{X}};\mathcal{D})$ for the conditional
quantiles using $\mathcal{D}$ as data
Procedure:

1:Split $\mathcal{Z}$ into a training fold
$\mathcal{Z}_{tr}\overset{\Delta}{=}({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})_{i\in\mathcal{I}_{tr}}$
and a calibration fold
$\mathcal{Z}_{cal}\overset{\Delta}{=}({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})_{i\in\mathcal{I}_{cal}}$

2:Train the models
$\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{x}};Z_{tr})$ and
$\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{x}};Z_{tr})$ using the
training set

3:Find and record the length of the $p\times|\mathcal{Z}_{cal}|$ sides for the
hyperrectangles based on the trained models from step 2. Denote these as
$A_{k}$ for the first side length on the $k$th data set point, $B_{k}$ for the
second side length on the $k$th data set point, etc.

4:For each $k\in\mathcal{I}_{cal}$ and each $j\in 1,\ldots,p$, compute the
score
$E_{k,j}=\max\\{\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{k})-Y_{k,j},Y_{k,j}-\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{k})\\}$

5:Convert the $E_{k,j}$ to the first dimension using the side lengths found in
step 3. $W_{11}=E_{11}$, $W_{12}=E_{12}\times\frac{A_{1}}{B_{1}}$,
$W_{21}=E_{21},W_{22}=E_{22}\times\frac{A_{2}}{B_{2}}$, etc.

6:Define $W_{k}=\max\limits_{j=1,\ldots,p}W_{k,j}$

7:Find $Adj_{1}=(1-\alpha)(1+\frac{1}{|Z_{cal}|})-\text{th empirical quantile
of }\\{W_{k}\\}$. This is the adjustment for the first dimension

8:Find the lengths of the $p$ test point sides using the trained models from
step 2 for the test point, ${\boldsymbol{x}}$. Denote the lengths of these
sides as $A_{n+1}$, $B_{n+1}$, $C_{n+1}$, etc.

9:The final adjustments for the other sides become
$Adj_{2}=Adj_{1}\times\frac{B_{n+1}}{A_{n+1}}$,
$Adj_{3}=Adj_{1}\times\frac{C_{n+1}}{A_{n+1}}$, etc.

Output:
$\hat{{\boldsymbol{C}}}({\boldsymbol{x}})=(\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{x}})-{\boldsymbol{Adj}},\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{x}})+{\boldsymbol{Adj}})$,
where $\hat{{\boldsymbol{q}}}_{\text{low}}({\boldsymbol{x}})$ and
$\hat{{\boldsymbol{q}}}_{\text{high}}({\boldsymbol{x}})$ are the p-dimensional
prediction sets for ${\boldsymbol{x}}$ defined in step 2 and
${\boldsymbol{Adj}}$ is the p-dimensional vector of adjustments found in steps
7 and 9

###  3.1 Remarks

One unanswered question is, what coverage rate should the initial
hyperrectangle be formed with for both point and quantile regression? We
recommended $1-\alpha$ in both Algorithm 1 and Algorithm 2. We recommended
$1-\alpha$ and not $1-\alpha/\binom{p}{2}$ because it is easier for a maximum
non-conformity score (which ours is) to expand rather than contract. If a
smaller volume is preferred, and at least one dimension has a skewed error,
one can do a grid search starting at a value less than $1-\alpha$ and
increasing up to $1-\alpha/\binom{p}{2}$. As a general note, this search can
also be done for quantile regression in the one-dimensional case to reduce the
interval length.

A note on which dimension to use as the reference. First, we will assume that
the ratio is constant, as it is in Section 2. Denote the non-conformity scores
for the $k$-th observation from the second calibration set in the $j$-th
dimension as $E_{k,j}$. The overall $k$-th non-conformity score with the first
dimension as the reference is then,

$W^{1}_{k}=\max\Big{\\{}E_{k,1},E_{k,2}\times\frac{A}{B},\ldots,E_{k,p}\times\frac{A}{P}\Big{\\}}.$

With the second dimension as the reference, the score would then be,

$W^{2}_{k}=\max\Big{\\{}E_{k,1}\times\frac{B}{A},E_{k,2},\ldots,E_{k,p}\times\frac{B}{P}\Big{\\}}.$

With the last dimension as the reference, the score would then be,

$W^{p}_{k}=\max\Big{\\{}E_{k,1}\times\frac{P}{A},E_{k,2}\times\frac{P}{B},\ldots,E_{k,p}\Big{\\}}.$

It is clear that, though $W^{j}_{k}$ will be different depending on which of
the $j$ dimensions is the reference, the (potentially) scaled $E_{k,j}$ that
is selected will be the same. That is, if $E_{k,2}\times\frac{A}{B}$ was the
maximum with the first dimension as the reference, $E_{k,2}$ will be the
maximum with the second dimension as the reference. So, the final score that
is chosen when the first dimension is the reference will always be the same as
the final score selected (scaled by a constant) when the $j$-th dimension is
the reference. Now, let $W_{q}$ represent the chosen final non-conformity
score. If the reference dimension is the first dimension, we have that,

$W^{1}_{q}=\max\Big{\\{}E_{q,1},E_{q,2}\times\frac{A}{B},\ldots,E_{q,p}\times\frac{A}{P}\Big{\\}}.$

With no loss of generality, assume that $E_{q,1}$ is the maximum of the two.
Then, the final adjustments are,

$Adj_{1}=E_{q,1},$ $Adj_{2}=E_{q,1}\times\frac{B}{A},$

and

$Adj_{p}=E_{q,1}\times\frac{P}{A}.$

Next, assume that the second dimension was the reference. We then have that,

$W^{2}_{q}=\max\Big{\\{}E_{q,1}\times\frac{B}{A},E_{q,2},\ldots,E_{q,p}\times\frac{B}{P}\Big{\\}}.$

Because $E_{q,1}$ was the maximum before, $E_{q,1}\times\frac{B}{A}$ will be
the maximum for the above and our final adjustments are,

$Adj_{1}=E_{q,1}\times\frac{B}{A}\times\frac{A}{B}=E_{q,1},$
$Adj_{2}=E_{q,1}\times\frac{B}{A},$

and

$Adj_{p}=E_{q,1}\times\frac{B}{A}\times\frac{P}{B}=E_{q,1}\times\frac{P}{A}.$

Because the second dimension is general, the reference dimension does not
change the results when the length ratios are a fixed quantity.

If our ratios are not the same, as in Section 3, then the reference dimension
can make a minor difference. Though, we’ve found in simulation studies that
the difference in volume between different reference dimensions is small. It
ultimately depends on the ratio output by the quantile regression models that
are being used. Because the reference dimension has a fixed length adjustment,
regardless of how much it varies with the covariates, it does make some
intuitive sense to set the “most homoskedastic” dimension to be the reference.
One potential approach to finding the most homoskedastic dimension is to look
at the variability of the side lengths output in the calibration set. The less
variable the side lengths are, the more homoskedastic the dimension should be,
assuming a reasonably good model.

##  4 Validity of Conformal Hyperrectangles

###### Theorem 1

For both conformal hyperrectangular regression and conformal quantile
hyperrectangular regression,

$\mathbb{P}({\boldsymbol{Y}}_{n+1}\in\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1}))\geq
1-\alpha$

If the $V_{i}$’s are almost surely distinct, then

$\mathbb{P}({\boldsymbol{Y}}_{n+1}\in\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1}))\leq
1-\alpha+\frac{1}{n_{cal}+1}.$

The proof is provided in Section 9.

##  5 Balance

Here, we introduce two assumptions under which the proposed methods can
achieve balance in some sense. Throughout this section, we assume that
$({\boldsymbol{X}}_{i},{\boldsymbol{Y}}_{i})\overset{i.i.d.}{\sim}P$.
Assumption ICED: There exist $\mu_{i},c_{i},i=1,\ldots,p$, where the $c$’s are
positive constants and the $\mu$’s are functions of ${\boldsymbol{X}}$ such
that $c_{1}\\{Y_{1}-\mu_{1}({\boldsymbol{X}})\\}|{\boldsymbol{X}}\sim
c_{2}\\{Y_{2}-\mu_{2}({\boldsymbol{X}})\\}|{\boldsymbol{X}}\sim\cdots\sim
c_{p}\\{Y_{p}-\mu_{p}({\boldsymbol{X}})\\}|{\boldsymbol{X}}$, that is, the
conditional distributions of $c_{i}\\{Y_{i}-\mu_{i}({\boldsymbol{X}})\\}$
given ${\boldsymbol{X}}$ are identical for $i=1,\ldots,p$, almost surely.
Assumption CQ: Assume the ICED assumption holds and denote $G$ as the common
distribution of $c_{i}\\{Y_{i}-\mu_{i}({\boldsymbol{X}})\\}$. Assume that the
quantile function of $G$ is a continuous function.

The first assumption will be simply referred to as the identical conditional
error distribution (ICED) assumption, while the second assumption is referred
to as the continuous quantile (CQ) assumption.

We now introduce the concept of asymptotic balance for conformal
hyperrectanges.

###### Definition 2

Asymptotic Balance

The hyperrectangle $\prod_{j}\hat{C}_{j}$ is said to achieve asymptotic
balance as $n_{cal}$ and $n_{tr}\to\infty$ if

$\max\limits_{j}\\{|\eta_{j}-\frac{1}{p}\sum\limits_{i=1}^{p}\eta_{i}|\\}\to
0\text{ in probability},$

where $\eta_{j}:=P(Y_{j}\notin\hat{C}_{j}),\ \forall j=1,\ldots,p$, where
$\hat{C}_{j}$ is the marginal conformal prediction interval for the $j$-th
dimension.

We now state the main results whose proofs are given in Lemma 9.

###### Theorem 3

For ${\boldsymbol{X}}={\boldsymbol{x}}\in\mathcal{C}$, where $\mathcal{C}$ is
a compact set, under the ICED assumption when the initial unconformalized
hyperrectangle is formed with a consistent estimator for the conditional
quantiles (that is, any conditional $\alpha\times 100\%$-quantile estimate
with fixed $\alpha\in(0,1)$ is uniformly consistent for
${\boldsymbol{x}}\in\mathcal{C}$), conformal quantile hyperrectangular
regression achieves asymptotic balance between dimensions for $\alpha\in(0,1)$
as $n_{tr}\to\infty$.

###### Theorem 4

For ${\boldsymbol{X}}={\boldsymbol{x}}\in\mathcal{C}$, where $\mathcal{C}$ is
a compact set, under the CQ assumption when the initial unconformalized
hyperrectangle is formed with an estimator for $\mu_{j}$ such that
$\sup_{x\in\mathcal{C}}|\hat{\mu}_{j}({\boldsymbol{x}})-\mu_{j}({\boldsymbol{x}})|\overset{a.s.}{\to}0\text{
as }n_{tr}\to\infty\ \forall j$ and an initial non-conformity score of
$V_{i,j}=|Y_{i,j}-\mu_{j}({\boldsymbol{X}})|$, conformal hyperrectangular
regression achieves asymptotic balance between dimensions for
$\alpha\in(0,1)$, as $n_{tr}\to\infty$ and $n_{cal1}\to\infty$.

The preceding results show that, under the ICED assumption or the CQ
assumption and with a large training and calibration sample size and a
reasonably good algorithm, it is expected that the hyperrectangles’ coverage
rates are expected to be reasonably balanced. Under the same assumptions, we
could use a Bonferroni correction to achieve balance, but our methods achieve
balance as well as the desired nominal coverage thanks to Theorem 1. A
Bonferroni correction would overcover both marginally and overall. In
otherwords, a Bonferroni correction is marginally (and overall) conservative,
while our methods allow us to find the right marginal coverage rates that lead
to tight overall coverage. This is corroborated by the numerical results given
in Section 6.

Next, we look at a favorable property of balance besides fairness. If instead
of the ICED assumption, we make the much stronger assumption that the
responses between dimensions are conditionally independent and identically
distributed, that is, $Y_{j}|{\boldsymbol{X}}$ are i.i.d.

Consider a hyperrectangle of the form $\prod_{i=1}^{p}R_{i}$. Suppose the
average marginal conditional coverage rate of the the marginal prediction sets
$R_{1},\ldots,R_{p}$ is equal to a fixed constant, say, $\tau>0$, that is,
$\sum_{j=1}^{p}P(Y_{j}\in R_{j}|{\boldsymbol{X}})\equiv p\tau$. We claim that
the maximum joint coverage probability of such hyperrectangles equals
$\tau^{p}$, which is attained if the hyperrectangle is furthermore balanced,
that is, $\forall j,\ P(Y_{j}\in R_{j}|{\boldsymbol{X}})\equiv\tau$. To see
this, by independence we have that,

$1-\alpha=P(Y_{1}\in R_{1},\ldots,Y_{p}\in R_{p}|{\boldsymbol{X}})$
$=\prod_{j=1}^{p}P(Y_{j}\in R_{j}|{\boldsymbol{X}}).$

Because the geometric mean is not greater than the arithmetic mean, it is
furthermore

$\leq\Bigg{[}\sum_{j=1}^{p}\frac{P(Y_{j}\in
R_{j}|{\boldsymbol{X}})}{p}\Bigg{]}^{p},$

with equality if and only if $P(Y_{j}\in R_{j}|{\boldsymbol{X}})$ are the same
$\forall j$, which must then be $\tau$. This completes the proof of the claim.
By achieving balance under the independence between dimensions assumption, we
are minimizing the sum of marginal coverage rates while maintaining the
overall coverage of the prediction region. This allows our prediction
intervals to be sharper than other methods with the same overall coverage.

The inequality that the product of coverage probabilities is less than the
joint probability was proven in Sidak (1967) for Normal random variables.
Dykstra et al. (1973) extended it to responses that are conditionally
independent and identically distributed between dimensions, the conditioning
variable need not be observed, so the case where the response variables are
(conditionally) exchangeable between dimensions applies to the bound:

$P(Y_{1}\in R_{1},\ldots,Y_{p}\in
R_{p}|{\boldsymbol{X}})\geq\prod_{j=1}^{p}P(Y_{j}\in R_{j}|{\boldsymbol{X}}).$
(5)

Indeed, even when the independence assumption is broken, our simulation
results show that when the responses are exchangeable between dimensions,
marginal coverage rates are just slightly below $(1-\alpha)^{1/p}$.

As a note of interest, Dykstra et al. (1973) first observed that as the
probabilities become large, the product of the marginal coverage probabilities
approaches the joint coverage probability when the marginal probabilities are
pairwise negatively dependent. Formally that is,

$|P(Y_{1}\in R_{1},\ldots,Y_{p}\in
R_{p}|{\boldsymbol{X}})-\prod_{j=1}^{p}P(Y_{j}\in
R_{j}|{\boldsymbol{X}})|\leq\frac{p-1}{2p}\Big{[}\log\Big{\\{}\prod_{j=1}^{p}P(Y_{j}\in
R_{j}|{\boldsymbol{X}})\Big{\\}}\Big{]}^{2}.$

Marginal coverage probabilities that make up a joint prediction set are not
necessarily pairwise negatively dependent, though, they can be in the
exchangeable response case.

##  6 Simulation Studies

###  6.1 Comparison of Four Methods

We looked at three dimensional responses with a varying range of correlation
between the error terms. In each simulation, the model was correctly
specified. We compared five approaches to forming the conformal prediction
sets. The first two used the conditional mean approach found in Section 2.
Data were split as follows: 500 for training, 250 for the first calibration,
and 250 for the second calibration. To form the initial hyperrectangles, we
used the absolute difference non-conformity score (Mean Absolute) and the
signed error non-conformity score (Mean Signed) with equal tails. The third
approach was our proposed quantile regression approach found in Section 3
(Quantile), with equal tails for the initial quantile model. Data were split;
500 for training and 500 for calibration. The final was the maximum absolute
difference,
$V=\max\limits_{j=1,\ldots,q}|Y_{j}-\hat{f}_{j}({\boldsymbol{X}})|$ (Absolute
Max), where $\hat{f}_{j}({\boldsymbol{X}})$ was an estimate of the conditional
mean. Data were split; 500 for training and 500 for calibration. All results
were rounded to three decimal places, with simulation standard errors in
parentheses if they were larger than 0.0005.

The data were generated as follows, $X_{1}\sim\text{exp}(\lambda=0.2)$,
$X_{2}\sim\text{Unif}(-5,5)$, $E(Y_{1}|{\boldsymbol{X}})=5+2X_{1}$,
$E(Y_{2}|{\boldsymbol{X}})=3X_{1}+X_{1}X_{2}$, and
$E(Y_{3}|{\boldsymbol{X}})=X_{2}^{2}$. The errors for the responses were
generated as $\epsilon_{1}\sim\text{Gamma}(\alpha=2,\lambda=0.2)$,
$\epsilon_{2}\sim\text{Gamma}(\alpha=3,\lambda=0.5)$, and
$\epsilon_{3}\sim\mathcal{N}(5X_{2},1)$ where $\lambda$ represents the rate
parameter and $\alpha$ represents the shape parameter. We started off with the
uncorrelated errors for all of the responses that are listed above,
${\boldsymbol{\epsilon}}=(\epsilon_{1},\epsilon_{2},\epsilon_{3})^{T}$, and a
desired correlation matrix between them, $R$. To generate the correlated
errors, we multiplied ${\boldsymbol{\epsilon}}$ by the upper-triangular
Cholesky decomposition of $R$. For each simulation, the mean or quantile
structure of the model was correctly specified. All parameters in the
specified model were estimated. The miscoverage rate, $\alpha$, was set to be
0.10 in all simulations. The simulation size was $1,000,000$ in all
simulations. In the four simulations the error correlation matrix was

$R_{1}=\begin{bmatrix}1&0.3&0.6\\\ 0.3&1&0.5\\\
0.6&0.5&1\end{bmatrix}\,\,R_{2}=\begin{bmatrix}1&0.8&0.8\\\ 0.8&1&0.8\\\
0.8&0.8&1\end{bmatrix}\,\,R_{3}=\begin{bmatrix}1&0.6&0.3\\\ 0.6&1&0.5\\\
0.3&0.5&1\end{bmatrix}\,\,R_{4}=\begin{bmatrix}1&0.2&0.2\\\ 0.2&1&0.2\\\
0.2&0.2&1\end{bmatrix}.$

The results of the simulations are given in Tables 1 \- 4. Coverage is the
average coverage over all simulations, volume is the average volume. Len i and
Marg i represent the length of the interval in the i-th dimension and the
marginal coverage in the i-th dimension, respectively.

Approach | Coverage | Volume | Len 1 | Len 2 | Len3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---|---
Mean Absolute | 0.899 | 23,499.007 (4.656) | 25.144 | 15.360 | 60.307 | 0.939 | 0.954 | 0.983
Mean Signed | 0.900 | 21,008.494 (2.377) | 25.711 | 14.776 | 55.242 | 0.959 | 0.950 | 0.968
Quantile | 0.900 | 15,274.272 (2.224) | 25.199 | 14.409 | 42.165 | 0.957 | 0.943 | 0.968
Absolute Max | 0.900 | 110,634.273 (8.948) | 47.972 | 47.972 | 47.972 | 0.991 | 1.000 | 0.906
Table 1: Simulation Result: Setup Number 1 Approach | Coverage | Volume | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---|---
Mean Absolute | 0.899 | 28,927.063 (5.099) | 23.457 | 21.044 | 58.178 | 0.930 | 0.937 | 0.965
Mean Signed | 0.900 | 29,657.267 (3.652) | 24.841 | 21.556 | 55.271 | 0.953 | 0.944 | 0.952
Quantile | 0.900 | 22,086.655 (3.165) | 24.120 | 20.885 | 43.764 | 0.946 | 0.934 | 0.951
Absolute Max | 0.900 | 121,806.996 (10.909) | 49.527 | 49.527 | 49.527 | 0.992 | 0.997 | 0.904
Table 2: Simulation Result: Setup Number 2 Approach | Coverage | Volume | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---|---
Mean Absolute | 0.899 | 24,423.676 (4.319) | 23.804 | 18.501 | 55.077 | 0.932 | 0.943 | 0.980
Mean Signed | 0.899 | 23,606.321 (2.528) | 24.883 | 18.357 | 51.620 | 0.954 | 0.942 | 0.968
Quantile | 0.900 | 16,050.704 (2.470) | 24.372 | 17.941 | 36.675 | 0.949 | 0.935 | 0.969
Absolute Max | 0.900 | 100,115.258 (7.246) | 46.408 | 46.408 | 46.408 | 0.990 | 0.999 | 0.908
Table 3: Simulation Result: Setup Number 3 Approach | Coverage | Volume | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---|---
Mean Absolute | 0.899 | 21,768.531 (4.820) | 25.522 | 14.445 | 58.368 | 0.941 | 0.954 | 0.993
Mean Signed | 0.899 | 18,133.270 (2.036) | 25.441 | 13.365 | 52.243 | 0.958 | 0.948 | 0.983
Quantile | 0.900 | 12,000.470 (1.938) | 25.068 | 13.410 | 35.678 | 0.956 | 0.943 | 0.985
Absolute Max | 0.900 | 96,679.263 (6.951) | 45.874 | 45.874 | 45.874 | 0.990 | 1.000 | 0.909
Table 4: Simulation Result: Setup Number 4

We can see from the simulation results that all four methods tended to achieve
the desired $90\%$ overall coverage. The absolute maximum approach has poor
length and volume compared to the other methods, and the balance (spreading
out the miscoverage between dimensions) is poor when the variance of the error
term for at least one dimension is much larger than other dimensions. The
signed and quantile approaches tend to have the best balance, but because of
the heteroskedastic error for the third response, the quantile rectangular
approach does a significantly better job with respect to the volume and the
length of the intervals.

###  6.2 Ten-dimensional Response

A second simulation was done for a ten-dimensional response. The simulation
setup is given in Table 5 and Table 6. The results are given in Table 7. Only
the quantile method of forming the rectangular prediction region was used.
Data were split; 500 for training and 500 for calibration. All of the
$\epsilon$ terms were generated from a multivariate Normal distribution with a
mean of 0, a variance of 1, and a correlation of 0.5. That is,
$\epsilon_{i}\sim\mathcal{N}(0,1)\text{ }i=1,\ldots,10$ and
$corr(\epsilon_{i},\epsilon_{j})=0.5$ for $i\neq j$ and $i,j=1,\ldots,10$. All
quantile models were fit with the correctly specified conditional quantiles,
which can be found in Table 6. We can also see that some of the error terms
had a heteroskedastic term added to them. The miscoverage rate, $\alpha$, was
set to be 0.10. The simulation size was $1,000,000$. The results of the
simulation can be found in Table 7.

Variable | Distribution | Parameters
---|---|---
$X_{1}$ | Uniform | min = -2; max = 5
$X_{2}$ | Uniform | min = -5; max = -1
$X_{3}$ | Uniform | min = -6; max = 10
$X_{4}$ | Uniform | min = 0; max = 4
$X_{5}$ | Uniform | min = $X_{2}$; max = $X_{4}$
Table 5: Simulation Setup: Ten Dimensions Variable | Conditional Mean | Error
---|---|---
$Y_{1}$ | $2X_{1}$ | $\epsilon_{1}$
$Y_{2}$ | $X_{1}+X_{1}X_{2}$ | $\epsilon_{2}$
$Y_{3}$ | $X_{2}^{2}$ | $\epsilon_{3}+5X_{2}$
$Y_{4}$ | $X_{2}X_{5}$ | $\epsilon_{4}+X_{5}^{2}$
$Y_{5}$ | $X_{5}^{2}$ | $\epsilon_{5}$
$Y_{6}$ | $X_{1}^{2}$ | $\epsilon_{6}-2X_{1}$
$Y_{7}$ | $X_{4}^{2}$ | $\epsilon_{7}-X_{4}$
$Y_{8}$ | $X_{3}^{2}$ | $\epsilon_{8}$
$Y_{9}$ | $X_{4}^{2}$ | $\epsilon_{9}$
$Y_{10}$ | $X_{1}X_{2}$ | $\epsilon_{10}$
Table 6: Simulation Setup: Ten Dimensions Dimension | Coverage | Length
---|---|---
Overall | 0.900 | NA
First | 0.987 | 5.043
Second | 0.986 | 5.040
Third | 0.990 | 7.048
Fourth | 0.981 | 17.506 (0.001)
Fifth | 0.987 | 5.045
Sixth | 0.995 | 11.457 (0.001)
Seventh | 0.987 | 5.251
Eighth | 0.987 | 5.044
Ninth | 0.987 | 5.045
Tenth | 0.987 | 5.045
Table 7: Simulation Result: Ten Dimensions

We can see, when the error terms are similar (in this case, the same for
dimensions 1, 2, 8, 9, 10), the balance is nearly perfect. Marginal coverage
tends to be higher for the dimensions with a heteroskedastic error, though the
length is not unreasonable. It is also clear that the overall desired coverage
of $90\%$ is achieved.

###  6.3 Checking Balance Between Dimensions

Two simulation scenarios were run with a 10 dimensional response. Only the
quantile method of forming the rectangular prediction region was used, one
with a homoskedastic error and one with a heteroskedastic error. Data were
split; 500 for training and 500 for calibration. In the homoskedastic
scenario, all of the $\epsilon$ terms were generated from a multivariate
Normal distribution with a mean of 0, a variance of 1, and a correlation of
0.9. That is, $\epsilon_{i}\sim\mathcal{N}(0,1)\text{ }i=1,\ldots,10$ and
$corr(\epsilon_{i},\epsilon_{j})=0.9$ for $i\neq j$ and $i,j=1,\ldots,10$. In
the heteroskedastic scenario, all of the $\epsilon|{\boldsymbol{X}}$ terms
were generated from a multivariate Normal distribution with a mean of 0, a
variance of $|X_{1}|$, and a correlation of 0.9. That is,
$\epsilon_{i}|{\boldsymbol{X}}\sim\mathcal{N}(0,|X_{1}|)\text{ }i=1,\ldots,10$
and $corr(\epsilon_{i}|{\boldsymbol{X}},\epsilon_{j}|{\boldsymbol{X}})=0.9$
for $i\neq j$ and $i,j=1,\ldots,10$. All quantile models were fit with the
correctly specified conditional quantiles, which can be found in Table 8 and
Table 9. Covariates were generated according to Table 5. The results of the
simulation scenario for the homoskedastic setup can be found in Table 10. For
the heteroskedastic scenario, they can be found in Table 11. The simulation
size was 1,000,000. We also compared our conformal quantile hyperrectangular
regression with marginal conformalized quantile regression and a Bonferroni
correction. Those results can be found in Table 12 and Table 13.

Variable | Conditional Mean | Error
---|---|---
$Y_{1}$ | $2X_{1}$ | $\epsilon_{1}$
$Y_{2}$ | $X_{1}+X_{1}X_{2}$ | $2\epsilon_{2}$
$Y_{3}$ | $X_{2}^{2}$ | $3\epsilon_{3}$
$Y_{4}$ | $X_{2}X_{5}$ | $4\epsilon_{4}$
$Y_{5}$ | $X_{5}^{2}$ | $5\epsilon_{5}$
$Y_{6}$ | $X_{1}^{2}$ | $\epsilon_{6}$
$Y_{7}$ | $X_{4}^{2}$ | $2\epsilon_{7}$
$Y_{8}$ | $X_{3}^{2}$ | $3\epsilon_{8}$
$Y_{9}$ | $X_{4}^{2}$ | $4\epsilon_{9}$
$Y_{10}$ | $X_{1}X_{2}$ | $5\epsilon_{10}$
Table 8: Simulation Setup: Homoskedastic Ten Dimensions Variable | Conditional Mean | Error
---|---|---
$Y_{1}$ | $2X_{1}$ | $\epsilon_{1}$
$Y_{2}$ | $X_{1}+X_{1}X_{2}$ | $2\epsilon_{2}$
$Y_{3}$ | $X_{2}^{2}$ | $3\epsilon_{3}$
$Y_{4}$ | $X_{2}X_{5}$ | $4\epsilon_{4}$
$Y_{5}$ | $X_{5}^{2}$ | $5\epsilon_{5}$
$Y_{6}$ | $X_{1}^{2}$ | $\epsilon_{6}$
$Y_{7}$ | $X_{4}^{2}$ | $2\epsilon_{7}$
$Y_{8}$ | $X_{3}^{2}$ | $3\epsilon_{8}$
$Y_{9}$ | $X_{4}^{2}$ | $4\epsilon_{9}$
$Y_{10}$ | $X_{1}X_{2}$ | $5\epsilon_{10}$
Table 9: Simulation Setup: Heteroskedastic Ten Dimensions Dimension | Coverage | Length
---|---|---
Overall | 0.900 | NA
First | 0.965 | 4.257
Second | 0.964 | 8.508
Third | 0.965 | 12.777 (0.001)
Fourth | 0.965 | 17.037 (0.001)
Fifth | 0.965 | 21.295 (0.001)
Sixth | 0.965 | 4.258
Seventh | 0.965 | 8.517
Eighth | 0.965 | 12.776 (0.001)
Ninth | 0.965 | 17.036 (0.001)
Tenth | 0.965 | 21.285 (0.001)
Table 10: Hyperrectangle Balance Simulation Result: Homoskedastic Ten Dimensions Dimension | Coverage | Length
---|---|---
Overall | 0.901 | NA
First | 0.958 | 6.396
Second | 0.960 | 12.703 (0.001)
Third | 0.960 | 19.054 (0.001)
Fourth | 0.960 | 25.394 (0.001)
Fifth | 0.960 | 31.747 (0.002)
Sixth | 0.960 | 6.350
Seventh | 0.960 | 12.699 (0.001)
Eighth | 0.960 | 19.047 (0.001)
Ninth | 0.960 | 25.395 (0.001)
Tenth | 0.960 | 31.738 (0.002)
Table 11: Hyperrectangle Balance Simulation Result: Heteroskedastic Ten Dimensions Dimension | Coverage | Length
---|---|---
Overall | 0.963 | NA
First | 0.990 | 5.381
Second | 0.990 | 10.979 (0.001)
Third | 0.990 | 16.152 (0.001)
Fourth | 0.990 | 21.552 (0.001)
Fifth | 0.990 | 26.945 (0.002)
Sixth | 0.990 | 5.384
Seventh | 0.990 | 10.768 (0.001)
Eighth | 0.990 | 16.155 (0.001)
Ninth | 0.990 | 21.529 (0.001)
Tenth | 0.990 | 26.933 (0.002)
Table 12: Bonferroni Balance Simulation Result: Homoskedastic Ten Dimensions Dimension | Coverage | Length
---|---|---
Overall | 0.967 | NA
First | 0.990 | 9.346 (0.001)
Second | 0.990 | 18.357 (0.002)
Third | 0.990 | 26.997 (0.002)
Fourth | 0.990 | 36.033 (0.003)
Fifth | 0.990 | 45.052 (0.004)
Sixth | 0.990 | 8.995 (0.001)
Seventh | 0.990 | 17.996 (0.002)
Eighth | 0.990 | 26.999 (0.002)
Ninth | 0.990 | 35.996 (0.003)
Tenth | 0.990 | 44.993 (0.004)
Table 13: Bonferroni Balance Simulation Result: Heteroskedastic Ten Dimensions

These simulation results provide corroboration with Section 5, that if the
conditional error terms follow the same distribution, up to a scale change,
our method will achieve good balance and the desired overall coverage with a
consistent model.

A second set of simulations were conducted to check the balance when the
responses are not conditionally exchangeable between $p$-dimensions, but still
follow the ICED assumption. Only the quantile method of forming the
hyperrectangular prediction region was used. As in the other simulations, the
models were fit with the correctly specified conditional quantiles. The data
were generated as follows, $X_{1}\sim\text{exp}(\lambda=0.2)$,
$X_{2}\sim\text{Unif}(-5,5)$, $E(Y_{1}|{\boldsymbol{X}})=5+2X_{1}$,
$E(Y_{2}|{\boldsymbol{X}})=3X_{1}+X_{1}X_{2}$, and
$E(Y_{3}|{\boldsymbol{X}})=5X_{2}+X_{2}^{2}$. In the homoskedastic scenario,
all of the error terms were generated from a standard Normal distribution. The
correlation between the first and second dimension was -0.8. All other
correlations between dimensions were 0.8. In the heteroskedastic scenario, all
of the error terms were generated from a Normal distribution with a mean of 0
mean and a variance of $|X_{1}|$. The correlation between the first and second
dimension was -0.8. All other correlations between dimensions were 0.8. The
simulation sizes were $1,000,000$. The results of the simulation are given in
Table 14.

We looked at a third set of simulations to see if balance is close to holding
when the ICED assumption is slightly violated. In each simulation,
$X_{1}\sim\text{exp}(\lambda=0.2)$, $X_{2}\sim\text{Unif}(-5,5)$,
$E(Y_{1}|{\boldsymbol{X}})=X_{1}$, $E(Y_{2}|{\boldsymbol{X}})=X_{1}$, and
$E(Y_{3}|{\boldsymbol{X}})=X_{1}$. In the homoskedastic case, the errors were
generated as
$\epsilon_{1}|{\boldsymbol{X}}\sim\text{Gamma}(\alpha=2,\lambda=0.2)$,
$\epsilon_{2}|{\boldsymbol{X}}\sim\text{Gamma}(\alpha=2,\lambda=0.2)$, and
$\epsilon_{3}|{\boldsymbol{X}}\sim\text{Gamma}(\alpha=2,\lambda=0.2)$ where
$\lambda$ represents the rate parameter and $\alpha$ represents the shape
parameter. In the heteroskedastic case, the errors were
$\epsilon_{1}|{\boldsymbol{X}}\sim\text{Gamma}(\alpha=2|X_{2}|,\lambda=0.2)$,
$\epsilon_{2}|{\boldsymbol{X}}\sim\text{Gamma}(\alpha=2|X_{2}|,\lambda=0.2)$,
and
$\epsilon_{3}|{\boldsymbol{X}}\sim\text{Gamma}(\alpha=2|X_{2}|,\lambda=0.2)$.
We started off with the uncorrelated errors for all of the responses that are
listed above,
${\boldsymbol{\epsilon}}=(\epsilon_{1},\epsilon_{2},\epsilon_{3})^{T}$, and a
desired correlation matrix between them, $R$. To generate the correlated
errors, we multiplied ${\boldsymbol{\epsilon}}$ by the upper-triangular
Cholesky decomposition of $R$. Because linear combinations of gamma
distributions do not generally follow a gamma distribution, the ICED
assumption is violated. For each simulation, the correctly specified quantile
structure of the model was estimated. Data were split; 2,000 for training and
200 for calibration. The miscoverage rate, $\alpha$, was set to be 0.10 in all
simulations. The simulation size was $1,000,000$ in all simulations. In the
four simulations the error correlation maxtrix was

$R_{1}=\begin{bmatrix}1&0.3&0.6\\\ 0.3&1&0.5\\\
0.6&0.5&1\end{bmatrix}\,\,R_{2}=\begin{bmatrix}1&0.8&0.8\\\ 0.8&1&0.8\\\
0.8&0.8&1\end{bmatrix}\,\,R_{3}=\begin{bmatrix}1&0.6&0.3\\\ 0.6&1&0.5\\\
0.3&0.5&1\end{bmatrix}\,\,R_{4}=\begin{bmatrix}1&0.2&0.2\\\ 0.2&1&0.2\\\
0.2&0.2&1\end{bmatrix}.$

The results of the simulations are given in Tables 15 \- 18. Coverage is the
average coverage over all simulations, volume is the average volume. Len i and
Marg i represent the length of the interval in the i-th dimension and the
marginal coverage in the i-th dimension, respectively.

Error | Coverage | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---
Homoskedastic | 0.900 | 3.981 | 3.977 | 3.978 | 0.951 | 0.950 | 0.950
Heteroskedastic | 0.900 | 7.970 | 7.961 | 7.968 | 0.951 | 0.949 | 0.949
Table 14: Simulation Result: Normal Error Three Dimensions Error | Coverage | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---
Homoskedastic | 0.896 | 26.155 | 26.499 | 27.073 | 0.963 | 0.959 | 0.954
Heteroskedastic | 0.897 | 41.718 | 41.949 | 42.403 | 0.960 | 0.959 | 0.958
Table 15: Simulation Result: Three Dimension Balance Comparison 1 Error | Coverage | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---
Homoskedastic | 0.897 | 24.733 | 25.413 | 25.504 | 0.953 | 0.944 | 0.941
Heteroskedastic | 0.897 | 42.063 | 42.326 | 42.765 | 0.961 | 0.959 | 0.958
Table 16: Simulation Result: Three Dimension Balance Comparison 2 Error | Coverage | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---
Homoskedastic | 0.896 | 26.195 | 26.913 | 26.904 | 0.963 | 0.958 | 0.955
Heteroskedastic | 0.897 | 41.751 | 42.286 | 42.267 | 0.960 | 0.958 | 0.958
Table 17: Simulation Result: Three Dimension Balance Comparison 3 Error | Coverage | Len 1 | Len 2 | Len 3 | Marg 1 | Marg 2 | Marg 3
---|---|---|---|---|---|---|---
Homoskedastic | 0.897 | 26.771 | 26.966 | 27.081 | 0.965 | 0.963 | 0.961
Heteroskedastic | 0.896 | 42.661 | 42.794 | 42.879 | 0.964 | 0.963 | 0.962
Table 18: Simulation Result: Three Dimension Balance Comparison 4

It is clear from these simulation results that balance is nearly achieved in
each simulation scenario. A note from these studies is that it appears the
correlation between the conditional error terms may impact how good the
balance is in finite samples. In Table 18 the marginal coverage rates are much
closer to each other, compared to Table 16, which are the setups with the
lowest correlation and highest correlation between dimensions, respectively.
We can also see that, although close, the bound in (5) holds as a strict
inequality for both the ten dimensional and three dimensional simulations.
Balance also appears to nearly hold with minor violations to the ICED
assumption.

One point that should be noted from all of the simulations is that the
rectangular approach to conformal prediction intervals maintains reasonable
balance without sacrificing the volume of the prediction regions. As we can
see from the simulations, conformalized hyperrectangles do a much better job
than marginal conformal prediction intervals with a Bonferroni correction. The
Bonferroni correction suffers from severe overcoverage when the marginal
distributions have correlated errors, which is likely to occur with a
multivariate response.

##  7 Real Data Analysis

All analyses were performed using R Statistical Software (R Core Team, 2023).
We looked at how well our conformal quantile hyperrectangle method worked in
predicting systolic and diastolic blood pressure for children. The blood
pressure data was obtained on children aged 9 and 11 from Bristol, UK via the
bp R package (Schwenck and Gaynanova, 2022). Covariates included in the
analysis were the child’s gender, age, height, weight, and BMI. The responses
of interest were the average systolic and diastolic blood pressure during the
subject’s first visit.

The number of subjects in the data set was 1,289. The data were randomly
permuted 200 times. For each permutation, 900 observations were used to train
the quantile neural network models, 200 were used to compute non-conformity
scores, and 189 were used to test the out of sample prediction interval
coverage for a 90% prediction interval. All models were computed with the qrnn
R package with one hidden layer with three nodes (Cannon, 2023). We also
compared the conformal quantile hyperrectangle results to the naïve Bonferroni
quantile model that used the trained Bonferroni adjusted quantile models
without a conformal adjustment as well as marginal conformalized quantile
regression that used a Bonferroni correction with a univariate conformal
adjustment. That is, the marginal coverage rates were set to 95% and we looked
at how often they both correctly captured the blood pressure. One without a
conformal adjustment and one with a conformal adjustment. Those comparisons
can be found in Table 19. Coverage is the average coverage for the out of
sample data over all permutations, volume is the average volume. Len i and
Marg i represent the length of the interval in the i-th dimension and the
marginal coverage in the i-th dimension, respectively.

Approach | Coverage | Marg 1 | Marg 2 | Len 1 | Len 2
---|---|---|---|---|---
Quantile Conformal Hyperrectangles | 0.899 | 0.941 | 0.935 | 45.104 | 42.573
Naïve Bonferroni | 0.859 | 0.914 | 0.917 | 42.685 | 40.085
Quantile Conformal Bonferroni | 0.915 | 0.947 | 0.949 | 48.820 | 46.240
Table 19: Blood Pressure Data Results

It is clear from the data analysis that, although the conformal quantile
hyperrectangle approach gives longer sides and a larger volume than the naïve
quantile neural network approach with a Bonferroni adjustment, the coverage is
much closer to the desired 90%. In fact, the naïve method undercovers while
the conformal Bonferroni correction overcovers. Looking at a boxplot of the
overall coverage comparisons Figure 1, it is easy to see that the naïve
approach frequently undercovers out of sample predictions while the conformal
Bonferroni correction overcovers. The conformal quantile hyperrectangular
approach is the only one that achieves the desired 90% nominal coverage on
average without being overly conservative.

Figure 1: Comparison of overall coverage. QCH represents the quantile
conformal hyperrectangle approach, NB represents the naïve Bonferroni
approach, and QCB represents the quantile conformal bonferroni approach.

##  8 Conclusion

The conformal rectangular approach introduced here is a new way of
constructing multivariate prediction intervals that do not require a
computationally expensive grid search. It controls the miscoverage rate in
finite samples, while allowing for heteroskedastic or homoskedastic errors. It
also helps maintain balance between the dimensions of the response. Though, we
may lose some precision in the model fitting process because we require the
data to be split, it’s clear from the real data analysis that slight loss in
model precision is worth the coverage guarantees.

More work can be done exploring other ways to compare the dimensions.
Conditional variance or standard deviation also seem like reasonable measures
of comparison. Another direction is to explore the use of Cartesian product of
prediction sets, for example, in the case of bimodal marginal distributions.
Also, based on the simulations in Section 6, there may also be weaker
assumptions to guarantee balance that can be further explored.

Acknowledgments and Disclosure of Funding

This project was partially funded by National Institutes of Health Predoctoral
Training Grant T32 HL 144461.

## References

  * Beran (1988) Rudolf Beran.  Balanced simultaneous confidence sets.  _Journal of the American Statistical Association_ , 83(403):679–686, 1988.  ISSN 01621459.  URL http://www.jstor.org/stable/2289291.
  * Bruder and Wolf (2017) Stefan Bruder and Michael Wolf.  Balanced bootstrap joint confidence bands for structural impulse response functions.  Working Paper 246, University of Zurich, Zurich, 2017.  URL https://hdl.handle.net/10419/162448.
  * Cannon (2023) Alex J. Cannon.  _qrnn: Quantile Regression Neural Network_ , 2023.  URL https://CRAN.R-project.org/package=qrnn.  R package version 2.1.
  * Diquigiovanni et al. (2022) Jacopo Diquigiovanni, Matteo Fontana, and Simone Vantini.  Conformal prediction bands for multivariate functional data.  _Journal of Multivariate Analysis_ , 189:104879, 2022.  ISSN 0047-259X.  doi: https://doi.org/10.1016/j.jmva.2021.104879.  URL https://www.sciencedirect.com/science/article/pii/S0047259X21001573.
  * Dykstra et al. (1973) R. L. Dykstra, John E. Hewett, and W. A. Thompson Jr.  Events which are Almost Independent.  _The Annals of Statistics_ , 1(4):674 – 681, 1973.  doi: 10.1214/aos/1176342462.  URL https://doi.org/10.1214/aos/1176342462.
  * Johnstone and Cox (2021) Chancellor Johnstone and Bruce Cox.  Conformal uncertainty sets for robust optimization.  In _Proceedings of the Tenth Symposium on Conformal and Probabilistic Prediction and Applications_ , volume 152 of _Proceedings of Machine Learning Research_ , pages 72–90. PMLR, 08–10 Sep 2021.  URL https://proceedings.mlr.press/v152/johnstone21a.html.
  * Kuleshov et al. (2018) Alexander Kuleshov, Alexander Bernstein, and Evgeny Burnaev.  Conformal prediction in manifold learning.  In _Proceedings of the Seventh Workshop on Conformal and Probabilistic Prediction and Applications_ , volume 91 of _Proceedings of Machine Learning Research_ , pages 234–253. PMLR, 11–13 Jun 2018.  URL https://proceedings.mlr.press/v91/kuleshov18a.html.
  * Lei and Wasserman (2013) Jing Lei and Larry Wasserman.  Distribution-free Prediction Bands for Non-parametric Regression.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 76(1):71–96, 07 2013.  ISSN 1369-7412.  doi: 10.1111/rssb.12021.  URL https://doi.org/10.1111/rssb.12021.
  * Lei et al. (2018) Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J. Tibshirani, and Larry Wasserman.  Distribution-free predictive inference for regression.  _Journal of the American Statistical Association_ , 113(523):1094–1111, 2018.  doi: 10.1080/01621459.2017.1307116.  URL https://doi.org/10.1080/01621459.2017.1307116.
  * Linusson et al. (2014) Henrik Linusson, Ulf Johansson, and Tuve Löfström.  Signed-error conformal regression.  In _Pacific-Asia Conference on Knowledge Discovery and Data Mining_ , pages 224–236. Springer, 2014.
  * Messoudi et al. (2020) Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau.  Conformal multi-target regression using neural networks.  In _Proceedings of the Ninth Symposium on Conformal and Probabilistic Prediction and Applications_ , volume 128 of _Proceedings of Machine Learning Research_ , pages 65–83. PMLR, 09–11 Sep 2020.  URL https://proceedings.mlr.press/v128/messoudi20a.html.
  * Neeven and Smirnov (2018) Jelmer Neeven and Evgueni Smirnov.  Conformal stacked weather forecasting.  In _Proceedings of the Seventh Workshop on Conformal and Probabilistic Prediction and Applications_ , volume 91 of _Proceedings of Machine Learning Research_ , pages 220–233. PMLR, 11–13 Jun 2018.  URL https://proceedings.mlr.press/v91/neeven18a.html.
  * Papadopoulos et al. (2002) Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman.  Inductive confidence machines for regression.  In _European Conference on Machine Learning_ , pages 345–356, 2002.
  * R Core Team (2023) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2023.  URL https://www.R-project.org/.
  * Robertson et al. (2021) Nathan Robertson, James M. Flegal, Dootika Vats, and Galin L. Jones.  Assessing and visualizing simultaneous simulation error.  _Journal of Computational and Graphical Statistics_ , 30(2):324–334, 2021.  doi: 10.1080/10618600.2020.1824871.  URL https://doi.org/10.1080/10618600.2020.1824871.
  * Romano and Wolf (2010) Joseph P. Romano and Michael Wolf.  Balanced control of generalized error rates.  _The Annals of Statistics_ , 38(1):598 – 633, 2010.  doi: 10.1214/09-AOS734.  URL https://doi.org/10.1214/09-AOS734.
  * Romano et al. (2019) Yaniv Romano, Evan Patterson, and Emmanuel Candes.  Conformalized quantile regression.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper_files/paper/2019/file/5103c3584b063c431bd1268e9b5e76fb-Paper.pdf.
  * Roy and Bose (1953) S. N. Roy and R. C. Bose.  Simultaneous Confidence Interval Estimation.  _The Annals of Mathematical Statistics_ , 24(4):513 – 536, 1953.  doi: 10.1214/aoms/1177728912.  URL https://doi.org/10.1214/aoms/1177728912.
  * Scheffe (1953) Henry Scheffe.  A method for judging all contrasts in the analysis of variance.  _Biometrika_ , 40(1/2):87–104, 1953.  ISSN 00063444.  URL http://www.jstor.org/stable/2333100.
  * Schwenck and Gaynanova (2022) John Schwenck and Irina Gaynanova.  _bp: Blood Pressure Analysis in R_ , 2022.  URL https://CRAN.R-project.org/package=bp.  R package version 2.1.0.
  * Sidak (1967) Zbynek Sidak.  Rectangular confidence regions for the means of multivariate normal distributions.  _Journal of the American Statistical Association_ , 62(318):626–633, 1967.  ISSN 01621459.  URL http://www.jstor.org/stable/2283989.
  * Tibshirani et al. (2019) Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas.  Conformal prediction under covariate shift.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper_files/paper/2019/file/8fb21ee7a2207526da55a679f0332de2-Paper.pdf.
  * Tukey (1953) John W. Tukey.  The problem of multiple comparisons.  Unpublished Manuscript, 1953.
  * Van Der Vaart (1998) A.W. Van Der Vaart.  _Asymptotic Statistics_.  Cambridge Series in Statistical and Probabilistic Mathematics, 3. Cambridge University Press, 1998.  ISBN 9780521496032.  URL https://books.google.com/books?id=udhfQgAACAAJ.
  * Vovk et al. (2005) Vladimir Vovk, Alex Gammerman, and Glenn Shafer.  _Algorithmic Learning in a Random World_.  Springer-Verlag, Berlin, Heidelberg, 2005.  ISBN 0387001522.

##  9 Appendix

The next two lemmas come from Romano et al. (2019) Similar results can be
found in Tibshirani et al. (2019); Lei et al. (2018); Vovk et al. (2005).

###### Lemma 5

(Quantiles and exchangeability). Suppose $Z_{1},\ldots Z_{n}$ are exchangeable
random variables.

For any $\alpha\in(0,1)$,

$\mathbb{P}(Z_{n}\leq\hat{Q}_{n}(\alpha))\geq\alpha,$

where $\hat{Q}_{n}$ is the empirical quantile function,
$\hat{Q}_{n}(\alpha)=Z_{(\lceil\alpha n\rceil)}$.

Moreover, if the random variables $Z_{1},\ldots,Z_{n}$ are almost surely
distinct, then

$\mathbb{P}(Z_{n}\leq\hat{Q}_{n}(\alpha))\leq\alpha+\frac{1}{n}.$

###### Lemma 6

(Inflation of Quantiles). Suppose $Z_{1},\ldots Z_{n}$ are exchangeable random
variables.

For any $\alpha\in(0,1)$,

$\mathbb{P}(Z_{n+1}\leq\hat{Q}_{n}((1+\frac{1}{n})\alpha))\geq\alpha.$

Moreover, if the random variables $Z_{1},\ldots,Z_{n}$ are almost surely
distinct, then

$\mathbb{P}(Z_{n+1}\leq\hat{Q}_{n}((1+\frac{1}{n})\alpha))\leq\alpha+\frac{1}{n+1}.$

The proof of Theorem 1

Proof  The validity of the proposed conformal hyperrectangles can be seen as
follows. For CHR, let $n_{cal}$ denote $n_{cal2}$ and $\mathcal{I}_{cal}$
denote $\mathcal{I}_{cal2}$. Call the ratios for converting between
dimensions, $R_{i}$, $i=1,\ldots,p$. If they depend on the covariates, call
them $R_{i}({\boldsymbol{X}}_{n+1})$, $i=1,\ldots,p$.

Denote $k$ as the smallest $\lceil(n_{cal}+1)(1-\alpha)\rceil$ value in
$\\{W_{i};i\in\mathcal{I}_{cal}\\}$.

Note that $W_{n+1}\leq k$ is equivalent to

$\iff\max\\{\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{n+1})-Y_{j,n+1},Y_{j,n+1}-\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{n+1})\\}\leq
k\times R_{j}^{-1}({\boldsymbol{X}}_{n+1}),\ \forall j\in\\{1,\ldots,p\\}$
$\iff Y_{j,n+1}\in[\hat{q}_{j,\text{low}}({\boldsymbol{X}}_{n+1})-k\times
R_{j}^{-1}({\boldsymbol{X}}_{n+1}),\hat{q}_{j,\text{high}}({\boldsymbol{X}}_{n+1})+k\times
R_{j}^{-1}({\boldsymbol{X}}_{n+1})],\ \forall j\in\\{1,\ldots,p\\}$
$\iff{\boldsymbol{Y}}_{n+1}\in\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1}).$

So we have that,

${\boldsymbol{Y}}_{n+1}\in\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1})\iff
W_{n+1}\leq k.$

Which tells us,

$\mathbb{P}({\boldsymbol{Y}}_{n+1}\in\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1}))=\mathbb{P}(W_{n+1}\leq
k).$

Because the original pairs of $({\boldsymbol{Y}}_{i},{\boldsymbol{X}}_{i})$
are exchangeable, so are the scores $W_{i}$ for $i\in\mathcal{I}_{cal}$ and
$i=n+1$. So, by Lemma 6,

$\mathbb{P}({\boldsymbol{Y}}_{n+1}\in\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1}))\geq
1-\alpha,$

and, under the assumption that the $W_{i}$’s are almost surely distinct,

$\mathbb{P}({\boldsymbol{Y}}_{n+1}\in\hat{{\boldsymbol{C}}}({\boldsymbol{X}}_{n+1}))\leq
1-\alpha+\frac{1}{n_{cal}+1}.$

###### Lemma 7

(Van Der Vaart, 1998, Lemma 2.2)

$X_{n}\overset{d}{\to}X$

holds iff, for all bounded, Lipschitz functions $g$,

$\lim_{n\to\infty}\mathbb{E}(g(X_{n}))=\mathbb{E}(g(X))$

###### Lemma 8

For a compact set, $\mathcal{C}$, assume that

$\sup_{x\in\mathcal{C}}|\hat{\mu}_{j}({\boldsymbol{x}})-\mu_{j}({\boldsymbol{x}})|\overset{a.s.}{\to}0\text{
as }n_{tr}\to\infty\ \forall j.$ (6)

Then we have that as $n_{tr}$ and $n_{cal1}$ $\to\infty$, the empirical
probability measure

$\frac{1}{n_{cal1}}\sum_{i=1}^{n_{cal1}}\delta(Y_{i,j}-\hat{\mu}_{j}({\boldsymbol{X}}_{i}))\overset{d}{\to}Y_{j}-\mu_{j}({\boldsymbol{X}}),\quad\forall
j,$

where $\delta(\cdot)$ denotes the Dirac delta probability measure at the
vector enclosed within the parentheses.

Proof  Let $f$ be a bounded Lipschitz function. Then we have for some constant
$K$ that,

$\Big{|}\frac{1}{n_{cal1}}\sum_{i=1}^{n_{cal1}}f(Y_{i,j}-\hat{\mu}_{j}({\boldsymbol{X}}_{i}))-\frac{1}{n_{cal1}}\sum_{i=1}^{n_{cal1}}f(Y_{i,j}-\mu_{j}({\boldsymbol{X}}_{i}))\Big{|}\leq
K\max_{i=1,\ldots,n_{cal1}}|\hat{\mu}_{j}({\boldsymbol{X}}_{i})-\mu({\boldsymbol{X}}_{i})|\overset{a.s.}{\to}0,$

as $n_{tr}\to\infty$, $\forall j,$ by (6). Now, by the strong law of large
numbers,

$\frac{1}{n_{cal1}}\sum_{i=1}^{n_{cal1}}f(Y_{i,j}-\mu_{j}({\boldsymbol{X}}_{i}))\to\mathbb{E}(f(Y_{j}-\mu_{j}({\boldsymbol{X}})))\text{
almost surely, as }n_{cal1}\to\infty,$

hence,

$\frac{1}{n_{cal1}}\sum_{i=1}^{n_{cal1}}f(Y_{i,j}-\hat{\mu}_{j}({\boldsymbol{X}}_{i}))\to\mathbb{E}(f(Y_{j}-\mu_{j}({\boldsymbol{X}})))\text{
as }n_{cal1}\text{ and }n_{tr}\to\infty.$

Then by Lemma 7,
$\frac{1}{n_{cal1}}\sum_{i=1}^{n_{cal1}}\delta(Y_{i,j}-\hat{\mu}_{j}({\boldsymbol{X}}_{i})\overset{d}{\to}Y_{j}-\mu_{j}({\boldsymbol{X}})$.

###### Lemma 9

(Van Der Vaart, 1998, Lemma 21.2) Let $X$ be a random variable with
distribution $F$. Define the quantile function of a cumulative distribution F
as $F^{-1}(p)=\inf\\{x:F(x)\geq p\\}$. Then, for any $0<p<1$,
$F_{n}^{-1}(p)\to F^{-1}(p)$ at all continuity points p of $F^{-1}$ if and
only if $X_{n}\overset{d}{\to}X$.

The proof of Theorem 3

Proof  Denote the oracle conditional quantiles as
$[q_{j,\text{low}}({\boldsymbol{x}}),q_{j,\text{high}}({\boldsymbol{x}})]$ and
the consistent estimators as
$[\hat{q}_{j,\text{low}}({\boldsymbol{x}}),\hat{q}_{j,\text{high}}({\boldsymbol{x}})]$.
Let the first dimension be the reference dimension, that is,
$c_{1}=1,\hat{c}_{1}({\boldsymbol{x}})\equiv 1$.

By consistency, we know that, uniformly in ${\boldsymbol{x}}\in\mathcal{C},\
j=1,\ldots,p$,
$P(|q_{j,\text{low}}({\boldsymbol{x}})-\hat{q}_{j,\text{low}}({\boldsymbol{x}})|>\delta_{\text{low}})\to
0$ for every $\delta_{\text{low}}>0$ as $n_{tr}\to\infty$ and
$P(|q_{j,\text{high}}({\boldsymbol{x}})-\hat{q}_{j,\text{high}}({\boldsymbol{x}})|>\delta_{\text{high}})\to
0$ for every $\delta_{\text{high}}>0$ as $n_{tr}\to\infty$. Hence, for
$j=2,\ldots,p$,

$\hat{c}_{j}({\boldsymbol{x}})=\frac{\hat{c}_{j}({\boldsymbol{x}})}{\hat{c}_{1}({\boldsymbol{x}})}=\frac{\hat{q}_{1,\text{high}}({\boldsymbol{x}})-\hat{q}_{1,\text{low}}({\boldsymbol{x}})}{\hat{q}_{j,\text{high}}({\boldsymbol{x}})-\hat{q}_{j,\text{low}}({\boldsymbol{x}})}$

converges in probability to $c_{j}/c_{1}=c_{j}$ uniformly in
${\boldsymbol{x}}\in\mathcal{C},\ j=2,\ldots,p.$ Hence, for all
$\delta_{j}>0,j=2,\ldots,p$ and uniformly in ${\boldsymbol{x}}$,
$P(|\hat{c}_{j}({\boldsymbol{x}})-c_{j}|>\delta_{j})\to 0$, in probability, as
$n_{tr}\to\infty$.

Define the symmetric set difference, $A\Delta B=(A\cap B^{C})\cup(A^{C}\cap
B)$.

For each $j=1,\ldots,p$, a large enough $n_{tr}$, and a constant $\ell$, by
consistency we have that

$\displaystyle P(Y_{j}$
$\displaystyle\in[q_{j,\text{low}}({\boldsymbol{x}})-\frac{1}{c_{j}}\ell,q_{j,\text{high}}({\boldsymbol{x}})+\frac{1}{c_{j}}\ell]\Delta[\hat{q}_{j,\text{low}}({\boldsymbol{x}})-\frac{1}{\hat{c}_{j}({\boldsymbol{x}})}\ell,\hat{q}_{j,\text{high}}({\boldsymbol{x}})+\frac{1}{\hat{c}_{j}({\boldsymbol{x}})}\ell])$
$\displaystyle\leq P(Y_{j}$
$\displaystyle\in[{q}_{j,\text{low}}({\boldsymbol{x}})-\delta_{\text{low}}-\frac{1}{c_{j}-\delta_{j}}\ell,{q}_{j,\text{low}}({\boldsymbol{x}})+\delta_{\text{low}}-\frac{1}{c_{j}+\delta_{j}}\ell]$
$\displaystyle\cup[{q}_{j,\text{high}}({\boldsymbol{x}})-\delta_{\text{high}}+\frac{1}{c_{j}+\delta_{j}}\ell,{q}_{j,\text{high}}({\boldsymbol{x}})+\delta_{\text{high}}+\frac{1}{c_{j}-\delta_{j}}\ell])$

Because $\delta_{\text{low}}$, $\delta_{\text{high}}$, and $\delta_{j}$ are
arbitrarily small positive numbers, as $n_{tr}\to\infty$

$P(Y_{j}\in[q_{j,\text{low}}({\boldsymbol{x}})-\frac{1}{c_{j}}\ell,q_{j,\text{high}}({\boldsymbol{x}})+\frac{1}{c_{j}}\ell]\Delta[\hat{q}_{j,\text{low}}({\boldsymbol{x}})-\frac{1}{\hat{c}_{j}({\boldsymbol{x}})}\ell,\hat{q}_{j,\text{high}}({\boldsymbol{x}})+\frac{1}{\hat{c}_{j}({\boldsymbol{x}})}\ell])\overset{P}{\to}0.$

By the ICED assumption, the scaled oracle set has the same miscoverage,
$\gamma\in(0,1)$, in each dimension. So using the ICED and the above,

$\eta_{j}=P(Y_{j}\notin[\hat{q}_{j,\text{low}}({\boldsymbol{x}})-\frac{1}{\hat{c}_{j}({\boldsymbol{x}})}\ell,\hat{q}_{j,\text{high}}({\boldsymbol{x}})+\frac{1}{\hat{c}_{j}({\boldsymbol{x}})}\ell])\overset{P}{\to}\gamma,\
j=1,\ldots,p.$

So,

$\max_{j}\\{|\eta_{j}-\frac{1}{p}\sum\limits_{i=1}^{p}\eta_{i}|\\}\overset{P}{\to}0.$

Which is the definition for asymptotic balance. Because $\ell$ is a general
constant, we can let it be the conformal adjustment, so conformal quantile
hyperrectangular regression achieves asymptotic balance under the listed
assumptions.

The proof of Theorem 4

Proof  Note the CQ assumption implies the validity of the ICED assumption. By
the ICED assumption we have that,

$|c_{j}(Y_{j}-\mu_{j}({\boldsymbol{X}})|\overset{\mathcal{D}}{=}G,\
j=1,\ldots,p,$

and by Lemma 8, the empirical probability distribution induced by
$c_{j}\\{Y_{i,j}-\hat{\mu}_{j}({\boldsymbol{X}}_{i})\\},i=1,\ldots,n_{cal1}$,
denoted by $G_{n,j}$, converges in distribution to $G$, as both $n_{tr}$ and
$n_{cal1}\to\infty$, for $j=1,\ldots,p$. Let
$c_{j}\hat{v}_{j}=G_{n,j}^{-1}(\alpha)$ and $c_{j}v_{j}=G^{-1}(\alpha)$, where
$G_{n,j}^{-1}(\alpha)$ and $G^{-1}(\alpha)$ are the empirical and true
$\alpha$ quantiles respectively. Note that $\hat{v}_{j}$ and $v_{j}$ are the
observed initial conformal expansions and the oracle initial conformal
expansions. By assumption CQ and Lemma 9,

$\hat{v}_{j}=\frac{G_{n,j}^{-1}(\alpha)}{c_{j}}\to\frac{G^{-1}(\alpha)}{c_{j}}=v_{j}\text{
as }n_{tr}\mbox{ and }n_{cal1}\to\infty,\ j=1,\ldots,p.$ (7)

Next, let the first dimension be the reference, that is,
$c_{1}=1,\hat{c}_{1}\equiv 1$. Then, our estimates of $\hat{c}_{j}$ are
$\frac{\hat{v}_{1}}{\hat{v}_{j}}$. By (7), for $j=2,\ldots p$,

$\hat{c}_{j}=\frac{\hat{c}_{j}}{\hat{c}_{1}}=\frac{\hat{v}_{1}}{\hat{v}_{j}}\to\frac{v_{1}}{v_{j}}=\frac{c_{j}}{c_{1}}=c_{j}\text{
as }n_{cal1}\to\infty.$

Define $\zeta>0$, $\epsilon>0$, and $\delta>0$. For a large enough $n_{tr}$
and $n_{cal1}$, and some constant, $\ell$, by the ICED assumption and the
consistency of $\hat{\mu}$, we have

$P(Y_{j}\in[\mu_{j}({\boldsymbol{x}})-\frac{1}{c_{j}}(v_{1}+\ell),\mu_{j}({\boldsymbol{x}})+\frac{1}{c_{j}}(v_{1}+\ell)]\Delta[\hat{\mu}_{j}({\boldsymbol{x}})-\frac{1}{\hat{c}_{j}}(\hat{v}_{1}+\ell),\hat{\mu}_{j}({\boldsymbol{x}})+\frac{1}{\hat{c}_{j}}(\hat{v}_{1}+\ell)])$
$\displaystyle\leq P(Y_{j}$
$\displaystyle\in[\mu_{j}({\boldsymbol{x}})-\zeta-\frac{1}{c_{j}-\delta}(v_{1}+\epsilon+\ell),\mu_{j}({\boldsymbol{x}})+\zeta-\frac{1}{c_{j}-\delta}(v_{1}-\epsilon+\ell)]$
$\displaystyle\cup[\mu_{j}({\boldsymbol{x}})-\zeta+\frac{1}{c_{j}-\delta}(v_{1}-\epsilon+\ell),\mu_{j}({\boldsymbol{x}})+\zeta+\frac{1}{c_{j}-\delta}(v_{1}+\epsilon+\ell)])$

Because $\zeta$, $\delta$, and $\epsilon$ are arbitrarily small numbers, as
$n_{tr}\to\infty$ and $n_{cal1}\to\infty$,

$P(Y_{j}\in[\mu_{j}({\boldsymbol{x}})-\frac{1}{c_{j}}(v_{1}+\ell),\mu_{j}({\boldsymbol{x}})+\frac{1}{c_{j}}(v_{1}+\ell)]\Delta[\hat{\mu}_{j}({\boldsymbol{x}})-\frac{1}{\hat{c}_{j}}(\hat{v}_{1}+\ell),\hat{\mu}_{j}({\boldsymbol{x}})+\frac{1}{\hat{c}_{j}}(\hat{v}_{1}+\ell)])\overset{P}{\to}0.$

By the ICED assumption, the scaled oracle set has the same miscoverage,
$\gamma\in(0,1)$, in each dimension. So, using the ICED and the above, as
$n_{tr}\to\infty$ and $n_{cal1}\to\infty$,

$\eta_{j}=P(Y_{j}\notin[\hat{\mu}_{j}({\boldsymbol{x}})-\frac{1}{\hat{c}_{j}}(\hat{v}_{1}+\ell),\hat{\mu}_{j}({\boldsymbol{x}})+\frac{1}{\hat{c}_{j}}(\hat{v}_{1}+\ell)]\overset{P}{\to}\gamma,\
j=1,\ldots,p.$

So,

$\max_{j}\\{|\eta_{j}-\frac{1}{p}\sum\limits_{i=1}^{p}\eta_{i}|\\}\overset{P}{\to}0.$

Which is the definition for asymptotic balance. Because $\ell$ is a general
constant, we can set it to be the final conformal adjustment, so CHR achieves
asymptotic balance under the listed assumptions.