Ordinal losses for classification of cervical cancer risk

Fixing CE with an ordinal loss term

A possible fix for CE is to add a regularization term that penalizes the deviations from the unimodal setting. Defining 1(x) as the indicator function of x and ReLU(x) = x1(x > 0) = max(0, x), a tentative solution for an order-aware loss could be (1)${\displaystyle \text{CO}\left({\mathbf{y}}_n,{\hat{\mathbf{y}}}_n\right)=\text{CE}\left({\mathbf{y}}_n,{\hat{\mathbf{y}}}_n\right)+\lambda \sum _{k=1}^{K-1}1\left(k\ge k_n^{\star }\right)\text{ReLU}\left({\hat{y}}_{n\left(k+1\right)}-{\hat{y}}_{n\left(k\right)}\right)+\lambda \sum _{k=1}^{K-1}1\left(k\le k_n^{\star }\right)\text{ReLU}\left({\hat{y}}_{n\left(k\right)}-{\hat{y}}_{n\left(k+1\right)}\right),}$ where λ ≥ 0 controls the relative importance of the extra terms favoring unimodal distributions. Predicted probability values are expected to decrease monotonously as we depart left and right from the true class. The added terms penalize any deviation from this expected unimodal distribution, with a penalty proportional to the difference of the consecutive probabilities. The additional terms, although promoting uni-modality, still allow flat distributions. A generalization of the previous idea is to add a margin of δ > 0 to the ReLU, imposing that the difference between consecutive probabilities is at least δ. This leads us to a second CE loss, CO2, suitable for ordinal classes: (2)${\displaystyle \text{CO2}\left({\mathbf{y}}_n,{\hat{\mathbf{y}}}_n\right)=\text{CE}\left({\mathbf{y}}_n,{\hat{\mathbf{y}}}_n\right)+\lambda \sum _{k=1}^{K-1}1\left(k\ge k_n^{\star }\right)\text{ReLU}\left(\delta +{\hat{y}}_{n\left(k+1\right)}-{\hat{y}}_{n\left(k\right)}\right)+\lambda \sum _{k=1}^{K-1}1\left(k\le k_n^{\star }\right)\text{ReLU}\left(\delta +{\hat{y}}_{n\left(k\right)}-{\hat{y}}_{n\left(k+1\right)}\right).}$ A value of δ = 0.05 has been empirically found to provide a sensible margin. This loss is aligned with the proposal present in Belharbi et al. (2019).

Beyond CO2: ordinal entropy loss function

In CO2, the CE term by itself is only trying to maximize the probability estimated in the true output class (while ignoring the remaining probabilities); the ordinal terms are promoting unimodality but not penalizing (almost) flat distributions. This also explains why the ordinal terms by themselves (especially the version without margin) are not enough to promote strong learning: the model could converge to solutions where the predicted probability in the true class is only slightly above the neighbouring probabilities, which will not, most likely, provide a strong generalization for new observations.

However, the extreme nature of CE, ignoring almost everything in the predicted distribution ${\hat{\mathbf{y}}}_n$ is equivalent to assuming that the perfect probability distribution is one on the true class and zero everywhere else. This assumes a strong belief and dependence on the chosen one-hot encoding, which is often a crude approximation to the true probability class distribution. Seldom, for a fixed observation x_n, the class is deterministically known; rather, we expect a class distribution with a few non-zero values. This is particularly true for observations close to the boundaries between classes. A softer assumption is that the distribution should have a low entropy, only.

This leads us to propose the ordinal entropy loss, HO2, for ordinal data as (3)${\displaystyle \text{HO2}\left({\mathbf{y}}_n,{\hat{\mathbf{y}}}_n\right)=\text{H}\left({\hat{\mathbf{y}}}_n\right)+\lambda \sum _{k=1}^{K-1}1\left(k\ge k_n^{\star }\right)\text{ReLU}\left(\delta +{\hat{y}}_{n\left(k+1\right)}-{\hat{y}}_{n\left(k\right)}\right)+\lambda \sum _{k=1}^{K-1}1\left(k\le k_n^{\star }\right)\text{ReLU}\left(\delta +{\hat{y}}_{n\left(k\right)}-{\hat{y}}_{n\left(k+1\right)}\right),}$ where H(p) denotes the entropy of the distribution p.

Data pre-processing

Given that all images from the Herlev dataset are of different sizes, all images were resized to 224 × 224 pixels; however, before the resize of cytological images, a zero-padding must be done to avoid the loss of essential information regarding cells shape. The last pre-processing step was to apply the same normalization as used by ImageNet (Simonyan & Zisserman, 2014).

Since the Herlev database has a relatively small number of observations (917), the training dataset was augmented by a series of random transformations: 10% of width and height shift, 10% of zoom, image rotation, horizontal and vertical flips, and color saturation. These transformations are illustrated in Fig. 3.

Convolutional neural networks

A convolutional neural network (CNN) is a neural network that successively applies convolutions of filters to the image. These filters are learned and consist of quadrilateral patches that are convolved across the whole input image—unlike previous fully-connected networks, only local inputs are connected at each layer. Typically, each convolution is intertwined with downsampling operations, such as max-pooling, that successively reduce the size of the original image.

The final layers are fully-connected and then the final output is processed by a soft-max for multi-class problems or a logistic function for binary classification. Dropout was used to reduce overfitting by constraining these fully-connected layers (Srivastava et al., 2014).

Network architectures

Two different models were trained and tested in this work for multi-class (4-class and 7-class) classification of Pap smear cells images (Fig. 4). Both models were trained and tested with nine different convolutional network architectures: AlexNet (Krizhevsky, Sutskever & Hinton, 2012), GoogLeNet (Szegedy et al., 2015), MobileNet_V2 (Howard et al., 2017), ResNet18 (He et al., 2016), ResNeXt50_32X4D (Xie et al., 2017), ShuffleNet_V2_X1_0 (Zhang et al., 2017), SqueezeNet1_0 (Iandola et al., 2016), VGG-16 (Simonyan & Zisserman, 2014), and Wide_ResNet50_2 (Zagoruyko & Komodakis, 2016). The goal of testing these different architectures is to evaluate how well the proposed loss behaves in a wide range of architectures. These nine different architectures were chosen as they are often used in the literature and came pre-trained with PyTorch on ImageNet (https://pytorch.org/docs/stable/torchvision/models.html). The last block of each architecture was replaced by the following layers: dropout with p = 20%, 512-unit dense layer with ReLU, dropout with p = 20%, a 256-wide dense layer with ReLU, followed by K output neurons.

A brief introduction of each architecture is now presented. AlexNet, based on LeNet, formalized the Convolutional Neural Network (CNN) as is known today: a series of convolutions intertwined by downsampling blocks. Max-pooling was used for downsampling and ReLU was used as the activation function. It became famous for winning ImageNet, the first CNN to do so (Krizhevsky, Sutskever & Hinton, 2012). The following ImageNet competitions were also won by other CNNs–VGG and GoogLeNet–which were evolutions on top of AlexNet that consist mostly of a much higher number of parameters (Simonyan & Zisserman, 2014; Szegedy et al., 2015). Then, MobileNet (Howard et al., 2017) introduced hyperparameters to help the user choose between latency and accuracy trade-offs. An attempt was then made at curbing the number of parameters with ShuffleNet (Zhang et al., 2018) by approximating convolution operators using fewer parameters.

Finally, an attempt was made at curbing the number of parameters, which had been exploding, while keeping the accuracy of these early CNNs with SqueezeNet (Iandola et al., 2016).

In another line of research, ResNet (He et al., 2016) introduced residual blocks whose goal was to make the optimization process easier for gradient descent. Each residual block learns a = f(x) + x instead of a = f(x). Given that weights are initialized randomly around zero and most activation functions are also centred in zero (an exception would be the logistic activation function), then, in expectation, all neurons output zero before any training. Therefore, when using residual blocks, at time = 0, a = x, i.e., activations produce the identity function. This greatly helps gradient descent focus on finding improvements (residuals) on top of the identity function. While this model allowed for deeper neural networks, each per cent of improved accuracy required nearly doubling the number of layers, which motivated WideResNet (Zagoruyko & Komodakis, 2016) and ResNeXt (Xie et al., 2017) to improve the residual architecture to improve learning time.

Training

The weights of the architectures previously mentioned are already initialized by pre-training on ImageNet. Adam was used as the optimizer and started with a learning rate of 10⁻⁴. The learning rate is reduced by 10% whenever the loss is stagnant for 10 epochs. The training process is completed after 100 epochs.

The dataset was divided into 10 different folds using stratified cross-validation, in order to maintain the class ratios. Therefore, the results are the average and deviation of these 10 folds. In the case of the proposed loss, the hyperparameter λ is tuned by doing nested k-fold cross-validating using the training set (with k = 5) in order to create an unbiased validation set.

Evaluation metrics

The most popular classification metric is accuracy (Acc). For N observations, taking k_i and ${\hat{k}}_i$ to be the label and prediction of the n-th observation, respectively, then $\text{Acc}=\frac{1}{N}{\sum }_{n=1}^{N}1\left({\hat{k}}_n^{\star }=k_n^{\star }\right),$ where 1 is the indicator function.

However, this metric treats all class errors as the same, whether the error is between adjacent classes or between classes in the extreme. If we have K classes represented by a set $\mathcal{C}=\left\{{\mathcal{C}}^{\left(1\right)},{\mathcal{C}}^{\left(2\right)},\dots ,{\mathcal{C}}^{\left(K\right)}\right\}$, then accuracy will treat an error between ${\mathcal{C}}^{\left(1\right)}$ and ${\mathcal{C}}^{\left(2\right)}$ with the same magnitude as an error between ${\mathcal{C}}^{\left(1\right)}$ and ${\mathcal{C}}^{\left(K\right)}$ which is clearly worse. As an illustration, in a medical setting, a misdiagnosis between Stage II and Stage III of a disease, while bad, is not as bad as a misdiagnosis between Healthy and Stage III. For that reason, a popular metric for ordinal classification is the Mean Absolute Error (MAE), $\text{MAE}=\frac{1}{N}{\sum }_i|k_i^{\star }-{\hat{k}}_i^{\star }|.$ This metric is not perfect since it treats an ordinal variable as a cardinal variable. An error between classes ${\mathcal{C}}^{\left(1\right)}$ and ${\mathcal{C}}^{\left(3\right)}$ will be treated as two times worse than an error between classes ${\mathcal{C}}^{\left(1\right)}$ and ${\mathcal{C}}^{\left(2\right)}$. Naturally, the assumption of cardinality is not always warranted.

To evaluate the models’ performance, we also used a metric specific for ordinal classification, Uniform Ordinal Classification Index (UOC) which considers accuracy and ranking in the performance assessment and is also robust against imbalanced classes (Silva, Pinto & Cardoso, 2018). The better the performance, the lower the UOC.

By combining a quality assessment (accuracy) with a quantity assessment (MAE) and also with a specific metric for ordinality (UOC) we hope to provide a balanced view of the performance of the methods.

The two other metrics used are the AUC of ROC or AUROC (Area Under the Receiver Operating Characteristic) and Kendall’s τ rank correlation coefficient. AUROC measures how well-calibrated are the probabilities produced by the model. This first metric is used in the binary classification context (two classes) and is extended for multi-class by comparing each class against the rest (one vs rest strategy) and performing an overall average, known as macro averaging. On the other hand, Kendall’s Tau is a non-parametric evaluation of relationships between columns of ranked data, so it is a measure of ordinal association between data. The τ correlation coefficient returns a value that ranges from −1 to 1, with 0 being no correlation and 1 perfect correlation.