How to threshold multiclass probability prediction to get confusion matrix?

Question

Lets say my multinomial logistic regression predict that a chance of a sample belonging to a each class is A=0.6, B=0.3, C=0.1 How do I threshold this values to get just binary prediction of a sample belonging to a class, taking in to an account imbalances of classes. I know what I would do if it's just a binary decision (threshold based on classes prevalence), or if the classes are balanced (classify to a class with highest probability). My end goal is to get 3x3 confusion matrix

Answer 1

According to @cangrejo's answer: https://stats.stackexchange.com/a/310956/194535, suppose the original output probability of your model is the vector $v$, and then you can define the prior distribution:

$\pi=(\frac{1}{\theta_1}, \frac{1}{\theta_2},..., \frac{1}{\theta_N})$, for $\theta_i \in (0,1)$ and $\sum_i\theta_i = 1$, where $N$ is the total number of labeled classes, $i$ is the class index.

Take $v' = v \odot \pi$ as the new output probability of your model, where $\odot$ denotes an element-wise product.

Now, your question can be reformulate to this: Finding the $\pi$ which optimize the metrics you have specified (eg. roc_auc_score) from the new output probability model. Once you find it, the $\theta s (\theta_1, \theta_2, ..., \theta_N)$ is your optimal threshold for each classes.

The Code part:

---

1. Create a proxyModel class which takes your original model object as an argument and return a proxyModel object. When you called predict_proba() through the proxyModel object, it will calculate new probability automatically based on the threshold you specified:

   class proxyModel():
       def __init__(self, origin_model):
           self.origin_model = origin_model
   
       def predict_proba(self, x, threshold_list=None):
           # get origin probability
           ori_proba = self.origin_model.predict_proba(x)
   
           # set default threshold
           if threshold_list is None:
               threshold_list = np.full(ori_proba[0].shape, 1)
   
           # get the output shape of threshold_list
           output_shape = np.array(threshold_list).shape
   
           # element-wise divide by the threshold of each classes
           new_proba = np.divide(ori_proba, threshold_list)
   
           # calculate the norm (sum of new probability of each classes)
           norm = np.linalg.norm(new_proba, ord=1, axis=1)
   
           # reshape the norm
           norm = np.broadcast_to(np.array([norm]).T, (norm.shape[0],output_shape[0]))
   
           # renormalize the new probability
           new_proba = np.divide(new_proba, norm)
   
           return new_proba
   
       def predict(self, x, threshold_list=None):
           return np.argmax(self.predict_proba(x, threshold_list), axis=1)
   

2. Implement a score function:

   def scoreFunc(model, X, y_true, threshold_list):
       y_pred = model.predict(X, threshold_list=threshold_list)
       y_pred_proba = model.predict_proba(X, threshold_list=threshold_list)
   
       ###### metrics ######
       from sklearn.metrics import accuracy_score
       from sklearn.metrics import roc_auc_score
       from sklearn.metrics import average_precision_score
       from sklearn.metrics import f1_score
   
       accuracy = accuracy_score(y_true, y_pred)
       roc_auc = roc_auc_score(y_true, y_pred_proba, average='macro')
       pr_auc = average_precision_score(y_true, y_pred_proba, average='macro')
       f1_value = f1_score(y_true, y_pred, average='macro')
   
       return accuracy, roc_auc, pr_auc, f1_value
   
   

3. Define weighted_score_with_threshold() function, which takes the threshold as input and return weighted score:

   def weighted_score_with_threshold(threshold, model, X_test, Y_test, metrics='accuracy', delta=5e-5):
       # if the sum of thresholds were not between 1+delta and 1-delta, 
       # return infinity (just for reduce the search space of the minimizaiton algorithm, 
       # because the sum of thresholds should be as close to 1 as possible).
       threshold_sum = np.sum(threshold)
   
       if threshold_sum > 1+delta:
           return np.inf
   
       if threshold_sum < 1-delta:
           return np.inf
   
       # to avoid objective function jump into nan solution
       if np.isnan(threshold_sum):
           print("threshold_sum is nan")
           return np.inf
   
       # renormalize: the sum of threshold should be 1
       normalized_threshold = threshold/threshold_sum
   
       # calculate scores based on thresholds
       # suppose it'll return 4 scores in a tuple: (accuracy, roc_auc, pr_auc, f1)
       scores = scoreFunc(model, X_test, Y_test, threshold_list=normalized_threshold)    
   
       scores = np.array(scores)
       weight = np.array([1,1,1,1])
   
       # Give the metric you want to maximize a bigger weight:
       if metrics == 'accuracy':
           weight = np.array([10,1,1,1])
       elif metrics == 'roc_auc':
           weight = np.array([1,10,1,1])
       elif metrics == 'pr_auc':
           weight = np.array([1,1,10,1])
       elif metrics == 'f1':
           weight = np.array([1,1,1,10])
       elif 'all':
           weight = np.array([1,1,1,1])
   
       # return negatitive weighted sum (because you want to maximize the sum, 
       # it's equivalent to minimize the negative sum)
       return -np.dot(weight, scores)
   

4. Use optimize algorithm differential_evolution() (better then fmin) to find the optimal threshold:

   from scipy import optimize
   
   output_class_num = Y_test.shape[1]
   bounds = optimize.Bounds([1e-5]*output_class_num,[1]*output_class_num)
   
   pmodel = proxyModel(model)
   
   result = optimize.differential_evolution(weighted_score_with_threshold, bounds, args=(pmodel, X_test, Y_test, 'accuracy'))
   
   # calculate threshold
   threshold = result.x/np.sum(result.x)
   
   # print the optimized score
   print(scoreFunc(model, X_test, Y_test, threshold_list=threshold))
   
   

Answer 2

This was helpful, thanks! But it is not applicable during model training. When it comes to using this method after training the model (after finding the hyperparameters relevant to the model), it's valid; only that there has to be some way of standardizing this to avoid loss of generality and for it to be applicable on test data.