Conformalized Ordinal Classification with Marginal and Conditional Coverage
Subhrasish Chakraborty Chhavi Tyagi Haiyan Qiao Wenge Guo



 















 
    
       
         
      





 


 


 



theoremTheorem
lemmaLemma
corollaryCorollary
propositionProposition
definition
remarkRemark
definitionDefinition


*assumption*

assumption[2]
  
  
  @edefcurrentlabel#1#2




kern[1]#1@kerna











ℝ
ℂ
ℤ
ℕ
𝔼
ℙ
𝖳
Cov
Var
⊥⊥
þ^th
tr
df
dim
col
row
null
rank
nullity
span
sign
supp
diag
aff
conv
dom
ŷ
f̂
μ̂
α̂
β̂
θ̂
𝒜
ℬ
𝒟
ℰ
ℱ
𝒢
𝒦
ℋ
ℐ
ℒ
ℳ
𝒩
𝒫
𝒮
𝒯
𝒲
𝒳
𝒴
𝒵






1]Subhrasish Chakraborty
1]Chhavi Tyagi
2]Haiyan Qiao
1]Wenge GuoAuthor e-mail addresses: sc2325@njit.edu, ct364@njit.edu, hqiao@csusb.edu, wenge.guo@njit.edu
[1]Department of Mathematical Sciences

       New Jersey Institute of Technology
[2]School of Computer Science & Engineering

       California State University San Bernardino






Conformalized Ordinal Classification with Marginal and Conditional Coverage
    [
    April 29, 2024
===========================================================================




Conformal prediction is a general distribution-free approach for constructing prediction sets combined with any machine learning algorithm that achieve valid marginal or conditional coverage in finite samples. Ordinal classification is common in real applications where the target variable has natural ordering among the class labels. In this paper, we discuss constructing distribution-free prediction sets for such ordinal classification problems by leveraging the ideas of conformal prediction and multiple testing with FWER control.  
Newer conformal prediction methods are developed for constructing contiguous and non-contiguous prediction sets based on marginal and conditional (class-specific) conformal p-values, respectively. Theoretically, we prove that the proposed methods respectively achieve satisfactory levels of marginal and class-specific conditional coverages. Through simulation study and real data analysis, these proposed methods show promising performance compared to the existing conformal method. 


KEY WORDS: Conformal prediction, ordinal classification, multiple testing, FWER control, marginal coverage, class-specific conditional coverage



§ INTRODUCTION



Ordinal classification, also known as ordinal regression or ordinal prediction, is a machine learning task that involves predicting a target variable with ordered categories <cit.>. In ordinal classification, the target variable has a natural ordering or hierarchy among its categories, but the intervals between the categories may not be evenly spaced or defined. Unlike regular classification, where the classes are nominal and unordered, ordinal classification takes into account the ordering relationship between the classes. This makes it suitable for situations where the outcome variable has multiple levels of severity, satisfaction ratings, or rankings. Here are a few examples of ordinal classification: customer satisfaction levels, movie ratings, disease severity levels, and education levels <cit.>. 

In ordinal classification, the goal is to learn a model that can accurately predict the ordinal variable's value given a set of input features. The model needs to understand the ordering of the classes and make predictions that respect this order. In the literature, 
some conventional classification algorithms have been adapted or modified to address ordinal classification, for example, ordinal logistic regression, SVM, decision trees, random forest, and neural networks <cit.>. Some alternative methods are also specifically developed for ordinal classification problems by fully exploiting the ordinal structure of the response variables
<cit.>. However, these existing methods can only provide point prediction, which is not adequate in some high
stakes areas such as medical diagnosis and automatic driving. Uncertainty quantification (UQ) techniques aim to go beyond point predictions and provide additional information about the reliability of these predictions. There are various techniques for UQ in machine learning, including Bayesian methods, calibration, and conformal prediction <cit.>. 

Conformal prediction is a unique distribution-free UQ technique that provides a prediction set rather than a point prediction for the true response with guaranteed coverage <cit.>. It can be used as a wrapper with any black-box algorithm. In this paper, we use the conformal prediction technique to construct prediction sets for 
ordinal classification problems. By combining the ideas of conformal prediction and multiple testing, two new conformal prediction methods are introduced for constructing contiguous and non-contiguous prediction sets.
Firstly, the problem of ordinal classification is reformulated as a problem of multiple testing; Secondly, for each constructed hypothesis, the marginal and conditional conformal p-values are respectively calculated; Thirdly, based on these marginal (conditional) conformal p-values, three multiple testing procedures are developed for controlling marginal (conditional) familywise error rate (FWER); Finally, based on the testing outcomes of these procedures, the prediction sets are constructed and proved having guaranteed marginal (conditional) coverage. 

There are almost no works of applying conformal prediction to address ordinal classification in the literature. To our knowledge, <cit.> is the only existing work, in which, a new (split) conformal prediction method is developed for constructing adaptive contiguous prediction region. This method is proved to have guaranteed marginal coverage, however, it cannot guarantee to have more desired conditional coverage. Moreover, it does not work well for high dimensional data. 
Compared to the method introduced in <cit.>, our proposed methods generally show via theoretical and numerical studies performing better in the settings of higher dimensions and in terms of class-specific conditional coverage; especially for the conditional conformal p-values based methods, they are proved to have guaranteed conditional coverage. 

The rest of this paper is structured as follows. In Section 2, we briefly introduce split-conformal prediction and review related works, followed by Section 3 which presents the development
of our proposed conformal methods using the idea of multiple testing. Section 4 provides numerical studies to evaluate the performance of the proposed methods compared to the existing method. Some discussions are presented in Section 5 and all proofs



§ PRELIMINARIES


In this section, we briefly describe the conformal prediction framework and review the related literature.


 §.§ Conformal Prediction


Conformal prediction is a general approach to construct prediction sets combined with a pre-trained classifier. The main advantage of this approach is that it is distribution-free and can work with any black-box algorithm. Conformal prediction is broadly of two types – full conformal prediction and split-conformal prediction. The full conformal prediction uses all the observations to train the black-box algorithms <cit.>. In contrast, split-conformal prediction <cit.> involves splitting the training data into proper training data to train the black-box algorithm and calibration data to calculate the threshold for forming prediction sets. Our proposed methods are based on the split-conformal method.
Consider a multi-class classification problem with feature space 𝒳 and labels 𝒴 = { 1,2,⋯,K }. Given the training observations  (X_i,Y_i)_i=1^2n and a test input X_2n+1, the goal is to find a prediction set C(X_2n+1) : 𝒳↦  2^𝒴 that contains the unknown response Y_2n+1 with enough statistical coverage.

The split-conformal procedure suggests to split 2n observations to n training observations, i.e., (X_i,Y_i)_i=1^n, which are used to train f̂, a black-box classifier such that f̂: 𝒳↦𝒴 and the remaining n observations (X_i,Y_i)_i=n+1^2n for calibration. The central part of this technique involves calculating the conformity scores for each observation, which measures how much the test observation conforms with the calibration observations. There can be several choices of conformity scores for multi-class classification problem, including posterior class probability, cumulative probability, and regularized cumulative probability <cit.>. Given the score function s : 𝒳×𝒴↦ℝ, the conformity score for the i^th calibration observation is defined as s_i = s(X_i,Y_i),  i=n+1,⋯,2n.

For a test input X_2n+1, we compute the conformity score for each class label. Therefore, for a class label y ∈𝒴 the conformity score corresponding to (X_2n+1,y) is s_2n+1 = s(X_2n+1,y). By using the conformity scores obtained for the calibration observations and the test input coupled with a given label y,  we can calculate the conformal p-value to test whether the unknown true label Y_2n+1 corresponding to the test input X_2n+1 is y or not. The (marginal) conformal p-value is defined as,

    p(X_2n+1, y) = ∑_i=n+1^2n𝕀(s_2n+1≤ s_i) + 1/n+1.

The final step involves constructing the prediction set C(X_2n+1) = {y:p(X_2n+1, y) ≥α}, which satisfies

    ℙ(Y_2n+1∈ C(X_2n+1)) ≥ 1 - α,

when the calibration and test observations (X_i,Y_i)_i=n+1^2n+1 are exchangeably distributed, where α∈ (0,1) is a pre-specified mis-coverage level. Equation (<ref>) is called marginal validity of the prediction set C(X_2n+1). It guarantees that the true label Y_2n+1 is contained in the prediction set with 100(1 - α)% confidence. <cit.> introduced another type of conformal p-value which is called as the conditional conformal p-value. Let 𝒟_cal = {n+1, …, 2n} denote the indices of the calibration observations 
(X_i,Y_i)_i=n+1^2n. For a test input X_2n+1 and any class y=1,…,K, the (class-specific) conditional conformal p-value given Y_2n+1 = y is defined as 

    p(X_2n+1|y) = ∑_i∈ℐ_y𝕀{s_i≤ s_2n+1} + 1/n_y + 1,

where ℐ_y = { i ∈𝒟_cal: Y_i = y }, n_y = |ℐ_y| is the size of ℐ_y, and
s_i = s(X_i,y) for i=n+1, …, 2n+1.

In general, the concept of conditional coverage such as object conditional validity and class-specific conditional validity are more relevant to practical applications <cit.>. If it is satisfied, the more desired results are often guaranteed.
Specifically, in classification problems, class-specific conditional validity provides conditional coverage for each given class, which is defined as

    ℙ(Y_2n+1∈ C(X_2n+1)|Y_2n+1 = y) ≥ 1 - α

for any y=1,…,K.
Proposition 1 and 2 below ensure that the marginal and conditional conformal p-values are valid, which result in desired marginal and (class-specific) conditional coverage.

To simplify the notation, we let Z_i = (X_i,Y_i) for i=1, …, 2n+1 and denote the conformity scores of the calibration data, {Z_i}_i=n+1^2n, as s_i's, and the conformity score of the test data, Z_2n+1=(X_2n+1,Y_2n+1), as s_2n+1, where Y_2n+1 is unknown. These notations are used in all propositions and theorems presented in this paper.



Suppose that {Z_i}_i=1^2n+1 where Z_i = (X_i,Y_i) are exchangeable random variables, then the marginal conformal p-values defined below as,


    p(Z_2n+1) = ∑_i=n+1^2n𝕀(s_i ≤ s_2n+1) + 1/n+1
  
is valid in the sense that for any t ∈ [0,1], we have  

    ℙ(p(Z_2n+1) ≤ t) ≤ t.

Moreover, if the conformity scores {s_i}_i=n+1^2n+1 are distinct surely, we have p(Z_2n+1) ∼ U{1/n+1,⋯,1}.




Suppose that {Z_i}_i=1^2n+1 where Z_i = (X_i,Y_i) are exchangeable random variables, then for any y ∈𝒴, given ℐ_y ⊆𝒟_cal and Y_2n+1 = y, the corresponding conditional conformal p-value as defined in equation (<ref>), is conditionally valid in the sense that for any t ∈ [0,1], 

    ℙ(p(X_2n+1|y) ≤ t  | ℐ_y,Y_2n+1 = y) ≤ t.

Moreover, if {s_i}_i∈ℐ_y ∪{2n+1} are distinct surely, we have that conditional on ℐ_y and Y_2n+1 = y,

    p(X_2n+1|y) ∼ U{1/n_y+1,⋯,1}.







 §.§ Related work

The framework of Conformal prediction was introduced by Vladimir Vovk and his collaborators <cit.> and has found many applications in classification problems. <cit.> and <cit.> provided a tutorial introduction and  brief literature review on this field. Several conformal methods have been developed to address binary classification <cit.> and multi-class classification problems <cit.>. Coverage guarantees of all these methods are established under the assumption of exchangeability. Very recently, some new conformal prediction methods have been developed in the settings of non-exchangeability  <cit.>. 

Although various conformal prediction methods have been developed for conventional classification problems, however, to our knowledge, <cit.> is the only reference that is specifically devoted to address ordinal classification problems using conformal prediction methods, in which an adaptive conformal method is developed for constructing contiguous prediction sets for ordinal response and is applied to AI disease rating in medical imaging. In addition, <cit.> is the closely related reference in which newer methods are developed for two types of loss functions specially designed for ordinal classification
in the more general framework of conformal risk control. 





§ METHOD


In  this section, we introduce several new conformal prediction methods for ordinal classification problems, in which there is a natural ordering among the classes labels. For simplicity, we assume a descending order of priority from class 1 to K in the response space 𝒴 = { 1,2,⋯,K }.



 §.§ Problem Formulation




We formulate the ordinal classification problem as a multiple testing problem. Specifically, by using the One-vs-All (OVA) strategy <cit.>,  for each class label, we construct a hypothesis to test whether or not a given test input X_2n+1 belongs to the particular class. The construction of the hypothesis is described as follows,

    H_i : Y_2n+1 = i   vs  H'_i : Y_2n+1≠ i,

for i=1,⋯,K. 
It is easy to see that all these hypotheses are random and there is only one true null. To test each individual hypothesis H_i, we use the corresponding marginal conformal p-value p(X_2n+1, i); to test H_1, …, H_K simultaneously, we consider the following three p-value-based testing procedures:

    
  * Procedure 1 : Test H_1, H_2, ⋯, H_K sequentially. The test is performed as follows.
    
        
  * If p(X_2n+1, 1) ≤α, reject H_1, move to test  H_2 else stop testing;
        
  * For i = 2,⋯, K-1, if p(X_2n+1, i) ≤α, reject H_i, move to test H_i+1 else stop testing;
        
  * If p(X_2n+1, K) ≤α, reject H_K else stop testing.
    
    
  * Procedure 2 : Test H_K, H_K-1,⋯, H_1 sequentially. The test is performed as follows.
    
        
  * If p(X_2n+1, K) ≤α, reject H_K, move to test H_K-1 else stop testing;
        
  * For i = K-1,⋯, 2, if p(X_2n+1, i) ≤α, reject H_i, move to test H_i-1 else stop testing;
        
  * If p(X_2n+1, 1) ≤α, reject H_1 else stop testing. 
    
    
  * Procedure 3 : Single-step procedure with common critical value α. This procedure rejects any hypothesis H_i if and only if p(X_2n+1, i) ≤α.
    


Procedure 1 and 2 are two pre-ordered testing procedures for which Procedure 1 follows the same testing order as that of the K classes  whereas 
Procedure 2 uses the reverse order of these classes <cit.>.  Procedure 3 is actually a conventional Bonferroni procedure for a single true null. 
Since there is only one true null among the K tested hypotheses, by Proposition 1, we have that all these three (marginal) conformal p-value based procedures strongly control family-wise error rate (FWER) at a pre-specified level α <cit.>. For each Procedure i, i = 1, 2, 3 defined above, the index set A_i of the accepted hypotheses is described as follows,

    
  * A_1 = {y_min, y_min+1,⋯, K }, where y_min = min{y∈𝒴:p(X_2n+1, y) > α};
    
  * A_2 = { 1, 2, ⋯, y_max}, where y_max = max{y∈𝒴:p(X_2n+1, y) > α};
    
  * A_3 = {y: p(X_2n+1, y) > α}.







 §.§ Ordinal Prediction Interval

Based on the the acceptance sets A_1 and A_2 of Procedure 1 and 2 given as above, we can obtain a new acceptance region A_12 = A_1 ∩ A_2 = {y_min, …, y_max}, which is used to define the prediction region C(X_2n+1) for the unknown response Y_2n+1. Specifically, the prediction region  C(X_2n+1) consists of the class labels for which the corresponding hypotheses are both accepted by Procedure 1 and 2, resulting in a contiguous set of labels {y_min, …, y_max}. This prediction region is referred to as a prediction interval in this context. The procedure for constructing the prediction interval C(X_2n+1) is summarized in Algorithm <ref> and illustrated in Figure <ref>. 

ruled
Comment/*  */












 §.§ Ordinal Prediction Set

Our second method for constructing ordinal prediction regions is based on Procedure 3. In this method, the prediction region C(X_2n+1) is defined simply using the acceptance region A_3 of Procedure 3, that is, C(X_2n+1) = {y ∈𝒴: p(X_2n+1, y) > α}. Specifically, the prediction region consists of any class labels for which the corresponding hypotheses are not rejected by Procedure 3, resulting in a non-contiguous set of labels. 
This prediction region is referred to as a prediction set in this context. The procedure for constructing the prediction set is detailed in Algorithm <ref> below.


ruled
Comment/*  */





In the following, we present two results regarding the FWER control of Procedure 1-3 and (marginal) coverage guarantees of Algorithm 1-2 introduced as above. 


Suppose that (X_i,Y_i)_i=n+1^2n+1 are exchangeable random variables, then Procedure 1-3 based on marginal conformal p-values, all strongly control the FWER at level α, i.e.,  FWER ≤α.
Specifically, if the conformity scores {s_i}_i=n+1^2n+1 are distinct surely, then for Procedure 3, we also have,
FWER ≥α - 1/n+1.



Suppose that (X_i,Y_i)_i=n+1^2n+1 are exchangeable random variables, then the prediction region C(X_2n+1) determined by Algorithm 1 and 2 both satisfy

    ℙ(Y_2n+1∈ C(X_2n+1)) ≥ 1 - α.

Specifically, for C(X_2n+1) determined by Algorithm 2, if the conformity scores {s_i}_i=n+1^2n+1 are distinct surely, we have 

    ℙ(Y_2n+1∈ C(X_2n+1)) ≤ 1 - α + 1/n+1.






 §.§ Class-specific conditional coverage


 To achieve more desired class-specific conditional coverage for our constructed prediction intervals and prediction sets, in Procedure 1-3 we use (class-specific) conditional conformal p-values p(X_2n+1|y) as described in equation (<ref>), instead of the marginal conformal p-values p(X_2n+1, y) for simultaneously testing H_1, …, H_K formulated in equation (<ref>). It is shown in Proposition 4 below that all these three modified procedures strongly control the conditional familywise error rate (FWER) at level α, i.e., FWER_y = ℙ(V > 0|Y_2n+1 = y) ≤α, for any y∈𝒴, where V is the number of type 1 errors. This result in turn leads to that the prediction regions C(X_n+1|y) constructed by Algorithm 1 and 2 based on p(X_2n+1|y) satisfy more desired (class-specific) conditional coverage, as stated in Theorem 2. 


Under the same exchangeability assumption as in Proposition 2, Procedure 1-3 based on conditional conformal p-values p(X_2n+1|y) all strongly control the conditional FWER at level α, i.e., for any y ∈𝒴, 

    FWER_y = ℙ{reject H_y|Y_2n+1 = y}≤α.

Specifically, if the conformity scores {s_i}_i∈ℐ_y ∪{2n+1} are distinct surely, then for Procedure 3 based on p(X_2n+1|y), we have that for any y ∈𝒴 and ℐ_y ⊆𝒟_cal.

    ℙ(reject H_y|Y_2n+1=y,ℐ_y ) ≥α - 1/n_y + 1.
 



Under the same exchangeability assumption as in Theorem 1, the prediction region C(X_2n+1|y) determined by Algorithm 1 or 2 based on conditional conformal p-values p(X_2n+1|y) satisfies

    ℙ(Y_2n+1∈ C(X_2n+1|y)  |Y_2n+1 = y) ≥ 1 - α
 for any y ∈𝒴.
Specifically, for the prediction set C(X_n+1|y) determined by Algorithm 2 based on p(X_2n+1|y), if the conformity scores {s_i}_i∈ℐ_y ∪{2n+1} are distinct surely, we have

    ℙ(Y_2n+1∈ C(X_2n+1|y)  | Y_n+1=y, ℐ_y ) ≤ 1 - α + 1/n_y+1
 
for any y ∈𝒴 and ℐ_y ⊆𝒟_cal.






§ NUMERICAL STUDY


In this section, we evaluate the performance of our four proposed methods, Ordinal Prediction Interval (OPI) in Algorithm 1 based on marginal conformal p-values (marginal OPI), the OPI based on conditional conformal p-values (conditional OPI), Ordinal Prediction Set (OPS) in Algorithm 2 based on marginal conformal p-values (marginal OPS), and the OPS based on conditional conformal p-values (conditional OPS), in comparison with the existing counterpart developed in Lu et al. (2022), Ordinal Adaptive Prediction Set (OAPS), on simulated data and one real dataset. The comparison is based on the marginal coverage, average set size, and class-specific conditional coverage of the prediction regions for a pre-specified level α. The empirical metric we use to measure the class-specific conditional coverage (CCV) of the above methods is defined as

    CCV = max_y ∈{ 1,2,⋯,K}{ (1-α) - P_y,0},

where P_y is the estimate of ℙ(Y ∈ C(X_2n+1|y) |Y = y)  and C(X_2n+1|y) is the prediction region obtained from Algorithm 1 or 2 using conditional conformal p-values for any y ∈𝒴. Intuitively, the metric measures the maximum of the deviance of the conditional coverage for each of the classes from the desired level of conditional coverage 1 - α. 

In the whole numerical investigations including simulation studies and real data analysis, we use the logistic regression algorithm as the black-box algorithm for our experiments and compute the conformity scores as estimated posterior probabilities of classes. 



 §.§ Simulations


We present the simulation study to evaluate the performance of our proposed methods along with the existing method. We consider two simulation settings below, a Gaussian mixture model and a sparse model. 


    
    
    
  * Gaussian mixture. (X|Y=k) ∼π_1 N_d(μ_k,Σ) + π_2 N_d(μ_k+1,Σ) for k = 1,2,3 and (X|Y = 4) ∼π_1 N(μ_4;Σ) + π_2 N(μ_1;Σ) with π_1 = 0.2, π_2 = 0.8.
   
    
In the above setting, we set d=2, μ_1 = (-1,0), μ_2 = (-1,-1), μ_3 = (0,-1), μ_4 = (1,-1),  and Σ as the equal correlation matrix with correlation = 0.1.

2. Sparse model. The sparse model is generated with different dimensions of feature vector with d = 5, 10, 20, 50,  or  100. The features are generated with X_i ∼ N(0,1) with Cov(X_i,X_j) = 0.5 for i ≠ j. The class labels are generated using the sigmoid function f(x) and the the following decision rule,
    
    Y_i = k  if  k-1/4≤ f(x_i) < k/4 ,   where  k = 1,2,3,4,

where f(x) = 1/(1 + e^- β^' x) with β = (β_1,⋯,β_d)^' and x = (x_1, …, x_d)^'.  
The value of β is set as (β_1,…,β_5) = (1,1,1,-√(2),1), and β_i = 0 for any 5 < i ≤ d. 

The sample size for these two simulation settings is 2,000, out of which 500 samples have been used to train the classifier, 525 observations for calibration, and 975 for validation. The simulations are repeated 500 times, and the results are averaged to obtain the final performance metrics.


Figure <ref> displays the performance of our proposed methods along with the existing method under simulation setting 1. It can be seen from the left panel of Figure <ref> that all these five methods empirically achieve the desired level of marginal coverage. The middle panel of Figure <ref> compares the set sizes of the prediction regions corresponding to these five methods. It can be seen from the figure that the marginal OPS and the conditional OPS use shorter set sizes to attain the proper marginal coverage than the existing OAPS  whereas the  OAPS has shorter set sizes than the marginal OPI and the conditional OPI.  Finally, while the marginal coverage is guaranteed by all the methods, the right panel of Figure <ref> shows their differences in class-specific conditional coverage; 
the existing OAPS exhibits the largest value of CCV compared to the proposed methods and among the four proposed methods, the conditional OPI and conditional OPS exhibit lower values of CCV than the marginal OPI and marginal OPS. 




Figure <ref> shows the performance of our proposed methods along with the existing method  under simulation setting 2 with different dimensions of inputs. It can be seen from the left panel of this figure that all these methods empirically achieve desired marginal coverage for lower dimensions, however, the existing OAPS massively undercovers for higher dimensions and thus loses the control of mis-coverage rate. From the middle and right panels of Figure <ref>, we can also see that the conditional OPI and conditional OPS achieve lower values of CCV than the existing OAPS, although the OAPS has the lower set sizes than our proposed methods
for various dimensions of input. 



 §.§ Application to real data


We also evaluate the performance of our proposed methods on a real dataset, Traffic accident data, which is publicly available on the website of Leeds City Council. The real data consists of 1,908 traffic accidents that occurred in the year 2019. The objective is to predict the severity of the casualties, which are classified into three categories – mild, serious, and fatal based on the features available. In the numerical experiment, 500 observations are used to train the logistic regression model, 35% of the remaining observations are used for calibration, and 65% for validation. Figure <ref> shows that all these methods empirically achieve desired marginal coverage for different levels of mis-coverage, however, the proposed marginal OPI and marginal OPS, and existing OAPS have lower set sizes than 
the conditional OPI and conditional OPS. It is also evident from Figure <ref> that the proposed conditional OPS and conditional OPI both attain desired class-specific conditional coverage unlike the existing method, OAPS, which seems to largely deviate from the desired level of conditional coverage. 







§ CONCLUDING REMARKS


In this paper, we discussed the ordinal classification problem in the framework of conformal prediction and introduced two types of conformal methods, OPI and OPS, for constructing distribution-free contiguous prediction regions and non-contiguous prediction sets, respectively. These methods are developed by leveraging the idea of 
multiple testing with the FWER control and are specifically designed based on marginal conformal p-values and (class-specific) conditional conformal p-values, respectively. Theoretically, it was proved that the proposed methods based on marginal and conditional p-values respectively achieve satisfactory levels of marginal and class-specific conditional coverages. Through some numerical investigations including simulations and real data analysis, our proposed methods show promising results for the settings of higher dimensions and for class-specific conditional coverage. 

This paper discussed constructing valid prediction set for single test input. It would be interesting to discuss how to construct (simultaneous) prediction sets for multiple test inputs with some overall error control such as false discovery rate (FDR) control. 
Another interesting extension might be to relax the conventional distributional assumption we used for classification problems, for which the training data and the test data follow from the same distribution. It will be interesting to see whether the proposed methods can be extended to the settings of distribution shift where the training and test data sets have different distributions. 

plainnat
   






§ PROOFS




 §.§ Proof of Proposition 1



Suppose that {Z_i}_i=1^n are exchangeable based on which the conformal p-values is defined as

    p(Z_n+1) = 1 + ∑_i=1^n𝕀(s_i ≤ s_n+1)/n+1

Suppose, for any given values of conformity scores, v_1, ⋯, v_n+1,  they can be rearranged as 
ṽ_1 < ⋯ < ṽ_k 
with repetitions n_i of ṽ_i such that ∑_i=1^k n_i = n+1. Let E_v denote the event of {s_n+1,⋯,s_2n+1} = {v_1,⋯,v_n+1}. Then, under E_v, for i = 1,⋯, k, we have

    ℙ(s_2n+1=ṽ_i|E_v) = n_i/n+1,

due to the exchangeability of s_i's, 

We also note that under E_v and s_2n+1 = ṽ_i we have from equation (<ref>),

    p(Z_2n+1) = ∑_l=1^i n_l/n+1.

Then, for any t ∈ [0,1] and i = 1, ⋯, k, we have

    ℙ(p(Z_2n+1) ≤ t|E_v,s_2n+1 = ṽ_i) =
    
    0    if  t < ∑_l=1^i n_l/n+1, 
    
    1    o.w.

Thus, for any i = 1, ⋯, k and ∑_l=1^i-1 n_l/n+1≤ t < ∑_l=1^i n_l/n+1, we have 

    ℙ(p(Z_2n+1) ≤ t|E_v) 
        = ∑_l=1^kℙ(p(Z_2n+1) ≤ t|E_v,s_2n+1 = ṽ_l) ·ℙ(s_2n+1 = ṽ_l|E_v) 
        = ∑_l=1^i-1 n_l/n+1≤ t.

By taking the expectation on the above inequality, it follows that the conformal p-value p(Z_2n+1) is marginally valid. 

Specifically, if conformity scores {s_i}_i=n+1^2n+1 are distinct surely,  then k = n+1 and n_i = 1 for i=1, …, n+1. Thus, 


    ℙ(p(Z_2n+1) ≤ t|E_v) = i-1/n+1,  if i-1/n+1≤ t < i/n+1,


that is, 

    p(Z_2n+1)|E_v ∼ U{1/n+1,⋯,1}.

This completes the proof.




 §.§ Proof of Proposition 2



For any given y ∈𝒴, the corresponding (class-specific) conditional conformal p-value is given by

    p(X_2n+1|y) = 1/n_y + 1[ ∑_i∈ℐ_y𝕀{s_i≤ s_2n+1} + 1],

where ℐ_y = {i:(X_i,Y_i) ∈𝒟_cal, Y_i = y}, n_y = |ℐ_y|,
s_i = s(X_i,y) for i ∈ℐ_y, and s_2n+1 = s(X_2n+1,y).
Given I_y and Y_2n+1=y, (X_i)_i ∈ℐ_y∪{X_2n+1} are exchangeably distributed, which is due to the assumption that  (X_i,Y_i)_i=1^2n+1 are exchangeably distributed. Using the similar arguments as in the proof of Proposition 1, for any given values of v_1, ⋯, v_n_y + 1, suppose that they can be arranged as  ṽ_1 < ⋯ < ṽ_l
with repetitions m_i of ṽ_i such that ∑_i=1^l m_i = n_y + 1. 

Let E_v denote the event {s_i}_i ∈ℐ_y∪{s_2n+1} = { v_1,⋯,v_n_y + 1}. Then, given E_v, I_y, and Y_2n+1 = y,  we have 

    ℙ(s_2n+1 = ṽ_i|E_v, I_y,Y_2n+1 = y) = m_i/n_y + 1

for  i = 1,⋯, l and y = 1,⋯, K, due to exchangeability of s_i, i ∈ I_y ∪{2n+1} given I_y,  which in turn is due to exchangeability of (X_i, Y_i)_i=1^2n+1. Note that given E_v,  I_y,  Y_2n+1 = y, and s_2n+1 = ṽ_i, we have from equation (<ref>),


    p(X_2n+1|y) = ∑_j=1^i m_j/n_y+1.

Thus, for any t ∈ [0,1] and i = 1,⋯,l,

    ℙ(p(X_2n+1|y) ≤ t | E_v, I_y, Y_2n+1 = y, s_2n+1 = ṽ_i ) =
    
    0    if  t < ∑_j=1^i m_j/n_y+1, 
    
    1    o.w.

Then, for any given i = 1,⋯,l and ∑_j=1^i m_j/n_y+1≤ t < ∑_j=1^i+1 m_j/n_y+1, we have 

    ℙ (p(X_2n+1|y) ≤ t  | E_v,I_y,Y_2n+1 =y ) 
        = ∑_j=1^lℙ(p(X_2n+1|y) ≤ t  |E_v,I_y,Y_2n+1=y,s_2n+1 = ṽ_i ) ·ℙ(s_2n+1 = ṽ_i |E_v,I_y,Y_2n+1=y ) 
        = ∑_j=1^i m_j/n_y+1≤ t.

By taking expectation, it follows that p(X_2n+1|y) is conditionally valid given Y_2n+1  = y.




 §.§ Proof of Proposition 3



Consider Procedure 1-3 based on marginal conformal p-values. Note that among the tested hypotheses H_1, …, H_K, there is exactly one hypothesis H_Y_2n+1 to be true. Thus, the FWER of Procedure 1-3 are all equal to

    ℙ(reject  H_Y_2n+1) ≤ℙ(p(X_2n+1,Y_2n+1) ≤α) ≤α,

where the last inequality follows by Proposition 1. 

Specifically,  for Procedure 3, if the conformity scores {s_i}_i=n+1^2n+1 are distinct surely, by Proposition 1, we have 

    FWER  = ℙ(reject  H_Y_2n+1) 
             = ℙ{p(X_2n+1,Y_2n+1) ≤α} ≥α - 1/n+1,

the desired result.





 §.§ Proof of Theorem 1



Note that the prediction set derived from Algorithm 1 is given by
C(X_2n+1) = A_1 ∩ A_2.
Thus, by Proposition 1,  
    ℙ(Y_2n+1∈ C(X_2n+1)) ≥ℙ(p(X_2n+1),Y_2n+1) > α) ≥ 1 - α.


Similarly, for Algorithm 2, its prediction set is given by C(X_2n+1) = { y ∈𝒴: p(X_2n+1, y) > α}.
By Proposition 1, it is easy to check that 

    ℙ(Y_2n+1∈ C(X_2n+1))  = ℙ(p(X_2n+1,Y_2n+1) > α) ≥ 1 - α.

Specifically, if the conformity scores {s_i}_i=1^2n+1 are distinct surely, for Algorithm 2, we have

    ℙ(Y_2n+1∈ C(X_2n+1)) = 1 - ℙ(p(X_2n+1,Y_2n+1) ≤α) ≤ 1-α + 1/n+1.

This completes the proof.





 §.§ Proof of Proposition 4


Consider Procedure 1-3 based on conditional conformal p-values. For any y = 1,⋯, K, given Y_2n+1 = y,  the conditional FWER of Procedure 1-3 are all equal to

    ℙ(reject  H_y|Y_2n+1 = y) ≤ℙ(p(X_2n+1|y) ≤α|Y_2n+1 = y) ≤α,

where the inequalities follow the definitions of Procedure 1-3 and Proposition 2. 

Specifically, for Procedure 3,  if the conformity scores {s_i}_i∈ℐ_y ∪{2n+1} are distinct surely, then by Proposition 2, 
the FWER conditional on I_y and Y_2n+1 = y is equal to

    ℙ{reject  H_y   | I_y, Y_2n+1 = y }  = ℙ(p(X_2n+1|y) ≤α|I_y, Y_2n+1 = y) ≥α - 1/n_y + 1.

This completes the proof.






 §.§ Proof of Theorem 2



By using Proposition 4 and the similar arguments as in the proof of Theorem 1, the prediction sets C(X_2n+1|y) derived from Algorithm 1 and 2 based on the conditional conformal p-values p(X_2n+1|y) all satisfy,

    ℙ{ Y_2n+1∈ C(X_2n+1 |y)  | I_y, Y_2n+1 = y } ≥ 1 - α

for any y = 1, …, K.
Specifically, if the conformity scores {s_i}_i∈ℐ_y ∪{2n+1} are distinct surely, for Algorithm 2, we have

    ℙ{ Y_2n+1∈ C(X_2n+1 |y)  | I_y, Y_2n+1 = y }    = ℙ{ p(X_2n+1|y) > α  | I_y, Y_2n+1 = y}
       ≤ 1 - α + 1/n_y + 1,

the desired result.




