# Interpretable pap smear cell representation for cervical cancer screening

Yu Ando  Department of Biomedical Science, Kyungpook National University  Nora
Jee-Young Park  Department of Pathology, School of Medicine, Kyungpook
National University  Department of Pathology, Kyungpook National University
Chilgok Hospital  Gun Oh Chong  Department of Obstetrics and Gynecology,
School of Medicine, Kyungpook National University  Department of Obstetrics
and Gynecology, Kyungpook National University Chilgok Hospital  Seokhwan Ko
Department of Biomedical Science, Kyungpook National University  Donghyeon Lee
Department of Biomedical Science, Kyungpook National University  Junghwan Cho
Corresponding author<EMAIL_ADDRESS>Clinical Omics Institute, Kyungpook
National University  Hyungsoo Han  Department of Biomedical Science, Kyungpook
National University  Clinical Omics Institute, Kyungpook National University

###### Abstract

Screening is critical for prevention and early detection of cervical cancer
but it is time-consuming and laborious. Supervised deep convolutional neural
networks have been developed to automate pap smear screening and the results
are promising. However, the interest in using only normal samples to train
deep neural networks has increased owing to class imbalance problems and high-
labeling costs that are both prevalent in healthcare. In this study, we
introduce a method to learn explainable deep cervical cell representations for
pap smear cytology images based on one class classification using variational
autoencoders.

Findings demonstrate that a score can be calculated for cell abnormality
without training models with abnormal samples and localize abnormality to
interpret our results with a novel metric based on absolute difference in
cross entropy in agglomerative clustering. The best model that discriminates
squamous cell carcinoma (SCC) from normals gives $0.908\pm 0.003$ area under
operating characteristic curve (AUC) and one that discriminates high-grade
epithelial lesion (HSIL) $0.920\pm 0.002$ AUC. Compared to other clustering
methods, our method enhances the V-measure and yields higher homogeneity
scores, which more effectively isolate different abnormality regions, aiding
in the interpretation of our results. Evaluation using in-house and additional
open dataset show that our model can discriminate abnormality without the need
of additional training of deep models.

##  1 Introduction

Cervical cancer screening is an invaluable test for early detection and
prevention of cervical cancer [1]. The Pap smear is a screening method that
examines cervical cells to discern abnormalities [2]. While pap smear
screening has been successful at reducing rates of cervical cancer mortality,
incidence and mortality increased in developing countries owing to the lack of
resources [3]. Developing quantitative and autonomous methods using computer
technology has the potential to improve healthcare worldwide by making the
smear screening process cost-efficient.

Recently, large models trained using supervised deep learning were proposed
and the obtained results have been promising [4, 5, 6]. Despite the positive
outlook, models were limited in providing explanations for their predictions.
Additionally, training these models requires many labeled datasets. In the
medical domain, it is expensive to acquire vast amounts of labeled data due to
the expertise required to interpret and generate data.

One-class classification (OCC) [7] is another machine-learning problem where
the objective is to train a model that can distinguish objects belonging to a
specific class, called the normal class, from all other objects or instances
that are considered anomalous. OCC has gained significant attention in recent
years owing to its potential applications in various domains, including
cybersecurity [8], fault detection [9], and medical tasks [10]. Owing to the
regular screening of cervical cancer [11], we can expect more samples of
normal-cell samples compared with the abnormal ones. This suggests that OCC
may be an invaluable method for screening. The goal is to capture the
intrinsic characteristics of the normal class and identify abnormalities
(deviations from normality).

Notwithstanding promising methods in deep learning, the inherent black-box
nature of deep models prevents the user from knowing how the model made a
prediction or inference. Models that can explain their predictions may gain
the trust of service providers and users and may support those involved with
informed decisions in fields that require informed consent or in cases
entailing high-stake decisions. Deep-learning methods that produce reasons for
their predictions will be more valuable in these domains.

Disentanglement is the property of representation where a unit of the
representation, an entry or entries of a feature vector, correspond to a
factor of variation of data [12]. For example, if a model is trained on images
of human faces, a disentangled representation may give entries that correspond
to age, hair length, facial color, etc. A model that has such property can be
used to interpret results by finding which feature is most important.
Variational autoencoder (VAE) based models [13, 14, 15], were found to have
this property. This is induced by regularizing independence between latent
entries or factors. We focus on these models to find morphological and color
factors in cervical cell image datasets to interpret abnormality.

In this work, we explore an unsupervised method for learning features of
normal cells and using negative log-likelihood to quantify abnormality. To
interpret our latent space, we examine the factors of variation learned by our
model and introduce a metric based on cross-entropy to cluster similar cells.
We test our method against baselines of other unsupervised clustering
algorithms. Our contributions are as follows:

  * •

We show that VAE models are capable of learning interpretable features from
Pap smear dataset

  * •

the trained models can discern abnormality by estimating a Gaussian on latent
feature space using only normal samples

  * •

additional image augmentation during training of generative models can improve
separation between normal and abnormal samples

  * •

a formulation of statistical distance based on cross-entropy with which
agglomerative clustering outperforms conventional clustering methods

  * •

Finally, our model is able to generalize to other datasets containing image of
cervical cells; it is capable of separating normal and abnormal images by
using the pre-trained encoder

##  2 Related Studies

This section reviews related studies concerning automated pap smear screening,
OCC using variational autoencoders, and clustering using statistical
distances.

###  2.1 Automated Pap Smear Screening

The development of a robust and accurate automated pap smear screening method
was motivated by reducing human errors and examination time [16]. This is
achieved by quantitative assessments of either hand-crafted features [17] or
the use of deep learning [4, 5, 6]. Methods using hand-crafted features are
interpretable in the sense that we can calculate the importance or weights of
features, unlike deep learning wherein the trained model is usually a black
box. A study that is most closely related to ours’ is that of Özbay and Özbay
[18] who proposed an interpretable deep learning framework for cervical cancer
detection using the ’hash’ layer. Unlike their method, by virtue of
variational autoencoder formulation, we are able to determine the qualitative
meaning of each feature dimension by interpreting the images generated by the
decoder along a fixed direction in feature space. This allows comparison
between images by calculating the change in latent space which quantitatively
show how the images are different.

###  2.2 OCC on VAEs

OCC using VAEs has the benefit of learning features for the given class in an
unsupervised manner. This method has been applied to detect deepfakes based on
the reconstruction loss and by comparing the latent features of the input and
the reconstructed image [19]. Another study applied the same method to
intrusion detection in network flow but the anomaly score set to the
reconstruction probability [20]. The advantage of our method is that
additional passes to the decoder are not required as we analyze only the
distribution in the latent space by using the negative likelihood in feature
space as our score. To the best of our knowledge, we are the first to train
our encoder using unsupervised learning and the first to apply OCC using VAEs
in Pap smear anomaly detection.

###  2.3 Clustering using Statistical Distances

Clustering is an unsupervised method used to determine groups of related
samples where interpoint distances are small compared with the points outside
a cluster [21]. In the context of Pap smear classification, several clustering
methods, K-means, and its fuzzy variant [22], have shown relatively good
performance [23, 24]. Unlike previous studies, ours is unique because we do
not use fixed points to cluster our samples. Instead, we utilized non-fixed
points, mapping the data into a multi-dimensional latent space as a
probability distribution. Analyzing data using non-fixed points [25] and
calculating statistical distances to compare class distribution have been
studied [26]. However, to our knowledge, we are the first to apply clustering
on non-fixed points by means of agglomerative clustering on nonfixed points
[27] using information-theory-based pseudometric. We also show that with some
additional assumptions, our pseudometric becomes equal to the Euclidean
distance which provides evidence for comparable but slightly better results
than those obtained using K-Means with the Euclidean distance.

##  3 Methods

In this section, we define our study’s problem about obtaining a score
function on the latent space of variational autoencoders (VAEs) by first
reviewing VAEs in the context of disentanglement. We then define our score
function for OCC and follow up using other methods for interpreting latent
space by qualitatively assessing reconstructions and using agglomerative
clustering method to obtain similar samples from input.

###  3.1 VAEs and Disentanglement

We denote $\theta_{d}$ as the parameters for the decoder, $\theta_{e}$ as the
parameters for the encoder, $q_{\theta_{d}}$ as the approximate posterior
parameterized as a diagonal multivariate normal distribution and
$p_{\theta_{g}}$ as the true data-generating distribution with marginal latent
distribution as a standard normal distribution. Suppose we have samples
$S=\\{\mathbf{x_{i}}\\}_{i=1}^{N}$ generated by hidden factors $\mathbf{z}$.
The evidence lower bound (ELBO) of VAE [13] is

$\mathcal{L}(\theta_{d},\theta_{e})=\mathcal{L}_{recon}(\theta_{d},\theta_{e})-\frac{1}{|S|}\sum_{\mathbf{x}\in
S}\left[D_{KL}(q_{\theta_{d}}(\mathbf{z}|\mathbf{x})\lVert
p_{\theta_{g}}(\mathbf{z}))\right],$ (1)

where
$\mathcal{L}_{recon}(\theta_{d},\theta_{e})=\frac{1}{|S|}\sum_{\mathbf{x}\in
S}E_{q_{\theta_{e}}(\mathbf{z}|\mathbf{x})}(p_{\theta_{d}}(\mathbf{x}|\mathbf{z}))$
and $D_{KL}$ is the Kullback-Leibler divergence [28].

Burgess et al [14] showed that, for curated datasets, adding an extra term
$\beta$ to the ELBO helps disentangle factors of variation[12] which may be
help interpret features. The ELBO for $\beta$-VAE is as follows

$\mathcal{L}_{\beta}(\theta_{d},\theta_{e};\beta)=\mathcal{L}_{recon}(\theta_{d},\theta_{e})-\frac{1}{|S|}\sum_{\mathbf{x}\in
S}\left[\beta D_{KL}(q_{\theta_{d}}(\mathbf{z}|\mathbf{x})\lVert
p_{\theta_{g}}(\mathbf{z}))\right]$ (2)

Chen et al [15] showed that total correlation [29], a term found by
decomposing ELBO[30], is related to the disentangling property found in
$\beta$-VAE. The ELBO for their model, $\beta$-TCVAE is as follows

L_β-TC(θ_d,θ_e;α, β, γ)=L_recon(θ_d, θ_e) - αI_q(z;x)
-βD_KL(q(z)∏_j q(z_j) - γ∑_j D_KL(q(z_j)p(z_j))

where $I$ is the mutual information function. $\alpha=1=\gamma$ was made
constant for all experiments in this study.

In this study, we explore the possibility of interpreting latent factors for
cells in Pap smear images using the aforementioned variations of VAEs. By
training using only normal samples, we hypothesized that color and
morphological factors may be found in the latent space so that we can
interpret results when inputs shows signs of abnormality.

###  3.2 Score Function

The ELBO suggests that our approximate posterior,
$q_{\hat{\theta}_{e}}(\mathbf{z}|\mathbf{x})$, should be a factorized standard
normal distribution after training. However, our distribution of normal
samples may differ from training samples owing to augmentations applied during
training. Instead, we also estimate our distribution of normal samples by a
multivariate normal distribution. We define our scoring function $s$ by the
negative log-likelihood of a multivariate normal distribution

$s(\mathbf{z})=-\log\mathcal{N}(\mathbf{z};\hat{\mu},\hat{V}),$ (3)

where $\hat{\mu}=\frac{1}{N}\sum_{\mathbf{x}}\mu(\mathbf{x};\hat{\phi}_{\mu})$
and
$\hat{V}=\frac{1}{N-1}\sum_{\mathbf{x}}(V(\mathbf{x};\hat{\phi}_{V})+\mu(\mathbf{x};\hat{\phi}_{\mu})\mu(\mathbf{x};\hat{\phi}_{\mu})^{T})-\hat{\mu}\hat{\mu}^{T}$.
$\hat{\phi}_{\mu}$ and $\hat{\phi}_{V}$ are parameters in $\hat{\theta}_{e}$.

###  3.3 Interpreting Latent Space

The formulation of VAEs suggests independent image-generating factors in
latent space. By traversing along the latent entry, we can infer qualitatively
the meaning of each factor in latent space. As we modeled the image
distribution as an augmented distribution, we traversed along each entry
centered at the mean of the nonaugmented data distribution for 10 steps
starting from a value of -4 x standard deviations (with respect to the mean)
to 4 x standard deviations (with respect to the mean).

Another way we can interpret latent space is to identify samples near the
encoded input. The straightforward approach is to cluster samples in latent
space and identify other samples in the cluster associated with the input.
However, a characteristic of variational autoencoders may allows us to use
statistical distances instead of sampling approaches to this problem. The
reparameterization trick[13] for a Gaussian posterior, a differentiable
transformation $g_{\theta_{e}}$, yields

$\mathbf{z}=g_{\theta_{e}}(\mathbf{x},\epsilon)=\mu(\mathbf{x};\theta_{e})+\sigma(\mathbf{x};\theta_{e})\epsilon,$
(4)

where $\epsilon\sim\mathcal{N}(\mathbf{0},\mathbf{1})$. This reveals that the
latent feature vector $\mathbf{z}$ is modeled by $\mu$ with an error term with
variance $\sigma^{2}$. Therefore, instead of using contemporary clustering
methods using the mean vector or a sample from this distribution (adding the
error term), we can use statistical distances. Images in the same cluster
should contain similar image characteristics owing to the proximity of their
location in latent space. By incorporating statistical distances, we
hypothesized that our method clusters similar images more effectively compared
with contemporary methods, such as K-means with Euclidean distance, by
utilizing the error of representation and factoring relative information of
nearby samples.

###  3.4 Cross-entropy-based Referenced Statistical Distance

We introduce a statistical distance that is related to well-known divergences
such as Jenson-Shannon divergence and Jefferys divergence, and an equality to
Mahalanobis distance as a special case which is shown in Appendix A. We call
this Cross-Entropy based Referenced Statistical Distance (CRSD).

Given two probability distribution, $p$ and $q$, the CRSD with respect to a
probability distribution $r$ is defined as

$d_{CRSD}(p,q;r)=|CE(r,q)-CE(r,p)|$ (5)

where $CE(p,q)$ is the cross-entropy of $q$ relative to $p$. Henceforth, for
brevity purposes we will eliminate the acronym. Based on the definition, it
can be inffered that the function is symmetric, $d(p,q;r)=d(q,p;r)$, non-
negative, and zero if $p=q$. However, $d(p,q;r)=0$ does not imply $p=q$ and
this can be easily demonstrated by Gaussian distributions with centered $r$;
$p$ and $q$ have means equal to $\mu$ and $-\mu$. The triangular inequality is
also satisfied and proved using simple algebra.

We also define for multiple reference probability densities: if $R$ is a set
of probability densities, the CRSD with respect to $R$ is defined by

$d(p,q;R)=\sum\limits_{r\in R}|CE(r,q)-CE(r,p)|$ (6)

In this study, we substitute Euclidean distance for CRSD in the agglomerative
clustering algorithm [27] when using non-fixed points. We restrict $R$ to be
finite and in our experiments to define the latent distributions of a given
data related to training, validation, or test fold.

###  3.5 Dataset

An open conventional Pap smear dataset provided by the Center for Recognition
and Inspection of Cells (CRIC) platform[31] was chosen for our experiment.
This dataset contains 400 Pap smear images and 11,534 classified cells which
contains more classified cells than other open datasets [32, 33]. The dataset
contains six classes according to the Bethesda system [34]: negative for
intraepithelial lesion or malignancy (NILM), low-grade epithelial lesion
(LSIL), high-grade epithelial lesion (HSIL), atypical squamous cell of
undetermined significance (ASC-US), atypical squamous cell cannot rule our
HSIL (ASC-H), and squamous cell carcinoma (SCC). As the classified cells were
annotated rather than cropped to separate images, different views of cells
(achieved by rotation and displacement) could be obtained during training
without introducing background elements generated by augmentation on cropped
images. In addition, we have removed some annotations that were close to the
boundary of the image (within 128 pixels in either x or y axis). The number of
cropped images used in our experiment is shown in Table 1.

Dataset | Class | Counts
---|---|---
CRIC | NILM | 5422
HSIL | 1,609
LSIL | 1,287
ASC-H | 894
ASC-US | 563
SCC | 156
SIPAKMED | Superficial-Intermediate | 831
Parabasal | 782
Koliocytotic | 814
Dyskeratotic | 794
Metaplastic | 785
In-house | Normal | 3,336,005
HSIL+ASC-H | 17,970
Table 1: Number of samples per class for each dataset used in this study.

[width=]figures/CellEncodingFigures_3.png

Figure 1: Example images of datasets used in our work. Tile images were
cropped according to location of cells (TOP): A) CRIC Dataset, B) SIPAKMED, C)
In-house normal, D) In-house abnormal. For in-house dataset, we only crop
cells within the bounding box (shown in red). The cropped image were augmented
further to capture the cell morphology which includes the nucleus and its
immediate surroundings. E) Examples of images used for training. We applied
three different augmentation techniques: None, random horizontal-flip and
affine (FA), random horizontal-flip, affine, and bightness adjustment (FAB).
F) Examples of in-house (top) and SIPAKMED (bottom) images used for external
evaluation.

To check for the generalizability of our model, we used SIPAKMED [33] and an
in-house dataset to compare scores between abnormal and normal cells. The
SIPAKMED dataset contained information on cell morphological information which
we used to determine the cell coordinates. We have retrospectively collected
100 normal and 50 abnormal patients including HSIL and ASC-H cells under the
approval Daegu Joint Institutional Review Board (DGIRB-2021-11-002). Our data
consisted of whole-slide images of liquid-based pap smears without precise
bounding boxes for individual cells. Whole-slide images containing abnormal
regions were annotated. To obtain bounding boxes for individual cells, we used
the Segment Anything Model (SAM) to detect cell locations [35]. From the 100
normal whole slide images and 26 abnormal whole slide images, we applied SAM
to 1024 by 1024 tiled images to obtain our cell images. To obtain only
abnormal cells from abnormal whole-slide images, we detected cells only within
manually annotated regions containing abnormal cells. The dimension of all
cell images were is 256 by 256; these images were center cropped to 64x64
before they were input into the model. This limits the field of view of our
model to focus on the nucleus of a single cell and its immediate surroundings.
Figure 1 shows the cropped images which were used as input.

###  3.6 Experiments

Our study consists largely of four parts: training of VAEs and the estimation
of normal sample distribution in latent space, abnormality scoring,
interpretation method for latent representations, and clustering. An overview
is shown in Figure 2.

We trained VAE, $\beta$-VAE, and $\beta$-TCVAE using normal samples (NILM) in
the CRIC dataset to learn disentangled features and to calculate a simple
scoring function to discern abnormal samples. Five-fold cross-validation was
performed on a two class classification task (NILM vs. rest) to identify the
best hyperparameters. The hyperparameters were varied using $\beta=1,4,16$,
and latent dimensions 8, 32, and 128.

[width=]figures/CellEncodingFigures_5_1.png

Figure 2: Overview of Interpretable Cervical Cancer Screening using VAE. (a)
Variants of VAE are trained on the CRIC dataset. After training, we estimate
parameters of multivariate guassian distribution. (b) To screen for
abnormality, we use negative log-likelihood in latent space using our
estimated distribution obtained in (a). (c) Generative latent variables are
characterized by their generating factor assesed qualitatively after training
the model. Images were generated by traversing in orthogonal directions with
equal steps. (d) Agglomerative clustering is performed with a statistical
pseudo-distance to cluster samples containing both normal and abnormal
samples. Similar samples, with respect to latent space, will be close to each
other.

Additionally, image augmentation techniques were varied to determine if
additional image preprocessing had a positive effect on model performance. In
generative models, severe augmentations tend to distort the learned data
distribution; thus we selected plausible augmentations, such as rotation,
translation, and brightness adjustment, that mimic possible variations in the
usual data acquisition process using an optical microscope. Three methods were
tested as hyperparameters in the cross-validation process: none, random flip,
rotation and translation (FA), and random flip, rotation, translation and
brightness (FAB). Example images are shown in Figure 1. For each abnormal
class, scores were then compared against NILM to assess the per-class
performance of our model. Finally, clustering algorithms were compared using
our best model. Clustering algorithms included agglomerative clustering [27],
K-means [36], spectral clustering[37], and DBSCAN [38]. We chose homogeneity,
completeness and $V$-measure[39] as a metric for comparison. We conducted
tests using 100 clusters and other parameters were set to default using sci-
kit learn implementation [40].

We experimented on additional external datasets to see if our model could
score abnormality with comparable or better performance. The model was fixed,
only to encode images to latent space. The distribution of normal images was
estimated using the same procedure as with the CRIC dataset. The datasets used
in the futher analysis were the SIPAKMED dataset and in house dataset obtained
from Kyungpook National University Chilgok Hospital. Unlike the SIPAKMED
dataset, the in-house dataset consisted of whole slide images which required
additional preprocessing in order to obtain cropped images of individual
cells. Bounding boxes containing HSIL were also given. We used Segment
Anything model (SAM)[35] to generate masks of cells. The center of the contour
was used as a reference point to crop a 256x256 image. For abnormal cells, we
applied a mask generated by bounding boxes from which candidate abnormal cells
were obtained using SAM. In total, we found 3,336,005 normal cell crops and
17,970 abnormal cell crops. During our evaluation process, we used repeated
random sampling of normal cells and abnormal cells. 6,000 normal cells were
sampled of which 3,000 was used for estimating the Gaussian and the
abnormality score was calculated on the other 3,000 normal cells and 3,000
abnormal cells. For SIPAKMED, we estimated Gaussian on different folds of
normal data (5-fold split on training datasets). Abnormal data were split into
5-folds each evaluated against normal data to calculate the receiver operating
characteristic (ROC) curve.

##  4 Results

###  4.1 Hyperparameter Tuning

In supervised learning, additional augmentations tend to enhance the
performances of deep models but this may not directly translate to generative
models because augmented samples distort the target data distribution to be
learned. Therefore, our experiments were conducted on relatively plausible
augmentations: rotation of cells, systematic error in centering a cell in an
image, and brightness configuration which depends on the light source of the
microscope. Figure 3 shows that for any hyperparameter combination, using the
full augmentation (FAB) yields the best area under the receiver operating
characteristic curve (AUROC) on average. The best model was found to be
$\beta$-TCVAE with $\beta=4$ and latent dimension of 8. Further analysis was
performed using this model.

[width=0.9]figures/hyperparameter_3.png

Figure 3: Cross-fold AUROC) ($\pm$ standard deviation.) for each
hyperparameter (latent dimension, $\beta$), data augmentation (None, FA or
FAB), and model configuration. Applying extra augmentation (in particular the
combination of horizontal, rotation, and brightness augmentations), AUROC
increased in all hyperparameteric combination cases. The best model was TCVAE
with $\beta$=4 and a latent dimension of eight.

###  4.2 Latent Feature Space

Morphological and color-related features were learnt by the model using
unsupervised training. Figure 4 shows features learnt by one of the models.
The reconstructed images were obtained by calculating the mean of latent
features of normal samples and traversing from -4 $\times$ standard deviations
to 4 $\times$ standard deviations from the mean in 10 steps. The top seven
rows correspond to morphology of the cell image concerning the locations,
shapes and sizes of the nucleus. The bottom two rows correspond to color. D-31
is a brightness adjustment and D-26 is a hue adjustment from gold-brownish to
blue.

[width=0.6]figures/Traversal.png

Figure 4: Unsupervised Morphological and Color Features. Features learned by a
model, $\beta$-VAE with $\beta=4$ and latent dimension of 32 trained with FAB
augmentation, are shown. Color (D-26), brightness (D-31), and several
morphological features were discovered by the model.

Pairwise KDE plots are shown for selected two latent features in Figure 5.
Each abnormal class contains different image characteristics which is evident
by differences in densities between classes. Notice that our estimated normal
sample distribution is not a standard Gaussian. KDE plots for the other
dimensions are in Supplementary Figure A1.

[width=0.8]figures/cric_distribution_figure.png

Figure 5: Pairwise KDE for latent dimensions 2 and 3 with generated images in
this plane. KDE plots (left) show that different classes have different
distributions. The classes are NILM (magenta), LSIL+ASC-US (green), HSIL+ASC-H
(blue), and SCC (red). With reference to the decoded images, we can compare
distributions of different classes. The images on the right are decoded images
in the latent plane (latent features 2 and 3). The rows correspond to latent
dimension 3 with increasing value from left to right and the columns
correspond to latent dimension 2 increasing in value from bottom to top. KDE
plots for the other dimensions are in Supplementary Figure A1.

###  4.3 Class-wise Screening Detection

We also calculated the performances of our scoring function by comparing each
abnormal class against normal samples. Table 2 lists accuracy, AUC, F1 scores,
sensitivity and specificity outcomes for each abnormal class. The model
performed well at discerning SCC and HSIL from normals followed by ASC-H,
LSIL, and ASC-US. The AUC score corresponds to severity of the lesion; AUC
scores are ordered from high grades HSIL, SCC, and ASC-H, to low grades LSIL
and ASC-US.

0.8!  class Accuracy AUC F1 score Sensitivity Specificity SCC 0.806±0.029
0.908±0.003 0.883±0.020 0.800±0.035 0.870±0.036 HSIL 0.849±0.002 0.920±0.002
0.850±0.006 0.829±0.027 0.871±0.026 ASC-H 0.742±0.005 0.808±0.006 0.787±0.007
0.731±0.018 0.761±0.020 LSIL 0.647±0.006 0.689±0.005 0.684±0.018 0.672±0.044
0.613±0.043 ASC-US 0.449±0.030 0.549±0.013 0.492±0.055 0.358±0.058 0.724±0.055

Table 2: Classwise Metrics for Best Model. We calculated metrics for each
abnormal class compared to normal dataset. The model was good at discerning
SCC and HSIL from normals followed by ASC-H, LSIL, and ASC-US.

###  4.4 Clustering Performance

Cross-validation (5-fold) on the abnormal dataset was performed to compare
different unsupervised algorithms using best model. Table 3 shows that
agglomerative methods yield higher $V$ measures than other methods with the
agglomerative clustering with statistical metric (Agg-SM) as the best
performing algorithm. Homogeniety for aggomerative clustering with Euclidean
distance (Agg-EM) and KMeans were similar but KMeans had lower $V$ measure.
Spectral clustering and DBSCAN had very low $V$ measures compared to other
algorithms but DBSCAN had highest completeness score suggesting that the
algorithm favors smaller number of clusters.

###  4.5 External Validation

!  Algorithm h c V Agg-SM 0.393$\pm$0.0117 0.151$\pm$0.00259 0.218$\pm$0.005
Agg-EM 0.282$\pm$0.00591 0.0940$\pm$0.00174 0.141 $\pm$ 0.00265 KMeans
0.275$\pm$0.006 0.0725$\pm$0.0019 0.114$\pm$0.003 Spectral Clustering
0.0526$\pm$0.0031 0.0950$\pm$0.0044 0.0676$\pm$0.0025 DBSCAN
0.000234$\pm$0.000413 0.746$\pm$0.414 0.000463$\pm$0.000818

Table 3: Unsupervised Clustering Algorithm Comparisons. Homogeneity (h),
completeness (c) and V-measure (V) are calculated for each algorithm with
number of clusters used. Our method, Agg-SM, gives best V-measure followed by
K-Means. Homogeneity was similar between agglomerative clustering using
Euclidean distance (Agg-EM) and K-Means.

[width=]figures/auroc_figures.png

Figure 6: ROC curves using SIPAKMED and in-house datasets

Finally, we evaluated the AUROC of discerning abnormality from normals by
using the score function for additional datasets: SIPAKMED and the in-house
dataset. Figure 6 shows receiver operating characteristic curves of each
abnormal class for each dataset. On average, the AUC values for dyskeratotic,
metaplastic, and koliocytotic cells were 0.97, 0.89, and 0.9 respectively. The
AUC value for HSIL+ASC-H on in-house data was 0.83. Our method shows good
generalizing performance outcomes when used on external datasets.

##  5 Discussion

VAEs enable the disentanglement of latent features that allow modeling of the
latent space of cervical cells in Pap smear images. Experiments showed that we
were successful in obtaining features corresponding to cytoplasmic color,
nuclear shape, and image brightness in an unsupervised manner while obtaining
respectable detection for SCC and HSIL. More importantly, a simple Gaussian
estimation on the latent feature space of normal examples produced these
results without the requirement of additional highly parameterized model
training.

Adding additional augmentations during the training of a generative model
improved the discerning capability. However, as the generative model learns
the data distribution, applying image augmentation can distort the original
data distribution [41]. This is also shown in Supplementary Figure A1, where
the normal class (NILM) is not represented by a centerted gaussian
distribution for each latent entry. We taken account of this change in latent
distribution by fitting a Gaussian distribution on the normal samples rather
than using the assumed prior of the data distribution.

A comparison of our results with those of other studies showed that there is a
trade-off using limited number of labeled samples for deep feature learning.
Table 4 shows the performance outcomes of other models with supervised
training. Despite training with only normal examples, our model was able to
discern abnormality with 0.7 accuracy in binary classification task. This may
be due to low AUC scores for low grade lesions (LSIL and ASC-US) as seen from
Table 2, which contributes to overall lower performance.

!  Sensitivity Specificity Accuracy F1 Training Method Diniz et al. (Random
Forest)[42] 0.958 0.958 0.958 0.958 Supervised Diniz et al. (Ensemble)[5] 0.96
0.96 0.96 0.96 Supervised CVM-Cervix [43] 0.916 0.9172 0.852 Supervised Ours
($\beta-tcvae$-b4d8) 0.759 0.683 0.70 0.58 OCC

Table 4: Binary classification results. OCC with VAEs perform worse than
supervised training.

This study was associated with some limitations. While our model can detect
abnormalities for most classes, it does not discriminate between abnormal
classes which supervised models do effectively. To overcome this, we devised a
querying method using agglomerative clustering. Images contained in some
single cluster will have similar features. Although we cannot label a
particular cluster as some particular abnormality, related images can be
queried and morphological and image related features can be assessed
quantitatively to assist pathologists examining cells in Pap smears. Despite
the fact that we have shown our model was effective at abnormality detection
on different datasets to the training dataset, external data will contain
different latent distributions depending on acquisition method as seen by
comparing latent spaces for CRIC, in-house and SIPAKMED ( Supplementary Figure
A1, A2, and A3 respectively).

##  6 Conclusion

VAEs provide means to interpret latent features by examining the reconstructed
output and give a simple formulation for scoring abnormality. OCC of normal
cells in Pap smear images can predict cell abnormality without using abnormal
labels. Performance was shown to increase when applying additional image
augmentation during training. Along with interpreting latent features
directly, a cross-entropy based metric was introduced for agglomerative
clustering. The method consistently outperformed other common unsupervised
clustering algorithms. We have shown that, although binary classification
performance is lower than supervised methods, using OCC can discern HSIL and
SCC from normals with high accuracy. We also demonstated a method for
explainability that may contribute to the development of explainable
artificial inteligence (AI) and AI trust in healthcare.

##  7 Acknowledgements

This work was supported by the Brain Pool program funded by the Ministry of
Science and ICT through the National Research Foundation of Korea
[2020H1D3A2A02102040, 2022H1D3A2A01096490]; Ministry of Education
[2021R1I1A3056903]

##  Appendix A Relation of CSRD to other divergences

The distance is equal to the Kullback-Leibler, Jefferys, and Jenson-Shannon
divergences depending on which reference probability distribution or a set of
probabilities is chosen. We denote the divergences as $D_{KL}$, $D_{J}$, and
$D_{JSD}$, respectively. The following are some relations of our formulation
to well-known divergences [44] d(p,q;p) &= D_KL(pq)
d(p,q;{p,q}) = 2 D_J(p,q)
d(p,q; {p,q}) ≥4 D_JSD(p,q) If $p$, $q$, are Gaussian distributions with equal
covariances, then we express the squared Mahalonobis distance as

$d(p,q;\\{p,q\\})=(\mu_{p}-\mu_{q})^{T}\Sigma^{-1}(\mu_{p}-\mu_{q}).$ (7)

Additionally, if $\Sigma=I$, we are then essentially calculating distances
between latent distributions by the Euclidean norm between the means, i.e.
$d(p,q;\\{p,q\\})=\lVert\mu_{p}-\mu_{q}\rVert_{2}^{2}$. Therefore, a
connection is made to our method using standard clustering algorithms using
the mean vectors of the latent representations to our method of using
statistical distances; the former is a special case of the latter with
assumptions on the covariances and the choice of the reference distributions.

## References

  * [1] E. E. Jansen, N. Zielonke, A. Gini, A. Anttila, N. Segnan, Z. Vokó, U. Ivanuš, M. McKee, H. J. de Koning, I. M. de Kok, et al., Effect of organised cervical cancer screening on cervical cancer mortality in europe: a systematic review, European Journal of Cancer 127 (2020) 207–223.
  * [2] S. L. Bedell, L. S. Goldstein, A. R. Goldstein, A. T. Goldstein, Cervical cancer screening: past, present, and future, Sexual medicine reviews 8 (1) (2020) 28–37.
  * [3] H. Sung, J. Ferlay, R. L. Siegel, M. Laversanne, I. Soerjomataram, A. Jemal, F. Bray, Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries, CA: a cancer journal for clinicians 71 (3) (2021) 209–249.
  * [4] H. Basak, R. Kundu, S. Chakraborty, N. Das, Cervical cytology classification using pca and gwo enhanced deep features selection, SN Computer Science 2 (5) (2021) 369.
  * [5] D. N. Diniz, M. T. Rezende, A. GC Bianchi, C. M. Carneiro, E. JS Luz, G. JP Moreira, D. M. Ushizima, F. NS de Medeiros, M. JF Souza, A deep learning ensemble method to assist cytopathologists in pap test image classification, Journal of Imaging 7 (7) (2021) 111.
  * [6] A. Manna, R. Kundu, D. Kaplun, A. Sinitca, R. Sarkar, A fuzzy rank-based ensemble of cnn models for classification of cervical cytology, Scientific Reports 11 (1) (2021) 14538.
  * [7] M. M. Moya, M. W. Koch, L. D. Hostetler, One-class classifier networks for target recognition applications, NASA STI/Recon Technical Report N 93 (1993) 24043.
  * [8] Y. Kim, H. K. Kim, Anomaly detection using clustered deep one-class classification, in: 2020 15th Asia Joint Conference on Information Security (AsiaJCIS), IEEE, 2020, pp. 151–157.
  * [9] J. Lee, Y. C. Lee, J. T. Kim, Fault detection based on one-class deep learning for manufacturing applications limited to an imbalanced database, Journal of Manufacturing Systems 57 (2020) 357–366.
  * [10] L. Gao, L. Zhang, C. Liu, S. Wu, Handling imbalanced medical image data: A deep-learning-based one-class classification approach, Artificial intelligence in medicine 108 (2020) 101935.
  * [11] S.-C. Lim, C. W. Yoo, Current status of and perspectives on cervical cancer screening in korea, Journal of pathology and translational medicine 53 (4) (2019) 210–216.
  * [12] Y. Bengio, A. Courville, P. Vincent, Representation learning: A review and new perspectives, IEEE transactions on pattern analysis and machine intelligence 35 (8) (2013) 1798–1828.
  * [13] D. P. Kingma, M. Welling, Auto-encoding variational bayes, arXiv preprint arXiv:1312.6114 (2013).
  * [14] C. P. Burgess, I. Higgins, A. Pal, L. Matthey, N. Watters, G. Desjardins, A. Lerchner, Understanding disentangling in beta-vae, arXiv preprint arXiv:1804.03599 (2018).
  * [15] R. T. Chen, X. Li, R. B. Grosse, D. K. Duvenaud, Isolating sources of disentanglement in variational autoencoders, Advances in neural information processing systems 31 (2018).
  * [16] B. Chitra, S. Kumar, Recent advancement in cervical cancer diagnosis for automated screening: a detailed review, Journal of Ambient Intelligence and Humanized Computing (2022) 1–19.
  * [17] P. Wang, L. Wang, Y. Li, Q. Song, S. Lv, X. Hu, Automatic cell nuclei segmentation and classification of cervical pap smear images, Biomedical Signal Processing and Control 48 (2019) 93–103.
  * [18] E. Özbay, F. A. Özbay, Interpretable pap-smear image retrieval for cervical cancer detection with rotation invariance mask generation deep hashing, Computers in Biology and Medicine 154 (2023) 106574.
  * [19] H. Khalid, S. S. Woo, Oc-fakedect: Classifying deepfakes using one-class variational autoencoder, in: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops, 2020, pp. 656–657.
  * [20] S. Zavrak, M. İskefiyeli, Anomaly-based intrusion detection from network flow features using variational autoencoder, IEEE Access 8 (2020) 108346–108358.  doi:10.1109/ACCESS.2020.3001350.
  * [21] C. M. Bishop, N. M. Nasrabadi, Pattern recognition and machine learning, Vol. 4, Springer, 2006.
  * [22] J. C. Bezdek, FUZZY-MATHEMATICS IN PATTERN CLASSIFICATION., Cornell University, 1973.
  * [23] M. Sharma, S. K. Singh, P. Agrawal, V. Madaan, Classification of clinical dataset of cervical cancer using knn, Indian Journal of Science and Technology 9 (28) (2016) 1–5.
  * [24] S. N. Sulaiman, N. A. Mat-Isa, N. H. Othman, F. Ahmad, Improvement of features extraction process and classification of cervical cancer for the neuralpap system, Procedia Computer Science 60 (2015) 750–759.
  * [25] M. Löffler, J. M. Phillips, Shape fitting on point sets with probability distributions, in: Algorithms-ESA 2009: 17th Annual European Symposium, Copenhagen, Denmark, September 7-9, 2009. Proceedings 17, Springer, 2009, pp. 313–324.
  * [26] X. Lu, Information mandala: Statistical distance matrix with clustering, IEEE Access 9 (2021) 56563–56577.
  * [27] D. Müllner, Modern hierarchical, agglomerative clustering algorithms, arXiv preprint arXiv:1109.2378 (2011).
  * [28] T. M. Cover, Elements of information theory, John Wiley & Sons, 1999.
  * [29] S. Watanabe, Information theoretical analysis of multivariate correlation, IBM Journal of research and development 4 (1) (1960) 66–82.
  * [30] M. D. Hoffman, M. J. Johnson, Elbo surgery: yet another way to carve up the variational evidence lower bound, in: Workshop in Advances in Approximate Bayesian Inference, NIPS, Vol. 1, 2016.
  * [31] M. T. Rezende, R. Silva, F. d. O. Bernardo, A. H. Tobias, P. H. Oliveira, T. M. Machado, C. S. Costa, F. N. Medeiros, D. M. Ushizima, C. M. Carneiro, et al., Cric searchable image database as a public platform for conventional pap smear cytology data, Scientific data 8 (1) (2021) 151.
  * [32] E. Martin, J. Jantzen, Pap-smear classification, Technical University of Denmark (2003) 2227–1899.
  * [33] M. E. Plissiti, P. Dimitrakopoulos, G. Sfikas, C. Nikou, O. Krikoni, A. Charchanti, Sipakmed: A new dataset for feature and image based classification of normal and pathological cervical cells in pap smear images, in: 2018 25th IEEE International Conference on Image Processing (ICIP), IEEE, 2018, pp. 3144–3148.
  * [34] D. Solomon, D. Davey, R. Kurman, A. Moriarty, D. O’Connor, M. Prey, S. Raab, M. Sherman, D. Wilbur, T. Wright Jr, et al., The 2001 bethesda system: terminology for reporting results of cervical cytology, Jama 287 (16) (2002) 2114–2119.
  * [35] A. Kirillov, E. Mintun, N. Ravi, H. Mao, C. Rolland, L. Gustafson, T. Xiao, S. Whitehead, A. C. Berg, W.-Y. Lo, et al., Segment anything, arXiv preprint arXiv:2304.02643 (2023).
  * [36] M. Ahmed, R. Seraj, S. M. S. Islam, The k-means algorithm: A comprehensive survey and performance evaluation, Electronics 9 (8) (2020) 1295.
  * [37] A. Ng, M. Jordan, Y. Weiss, On spectral clustering: Analysis and an algorithm, Advances in neural information processing systems 14 (2001).
  * [38] M. Ester, A density-based algorithm for discovering clusters in sarge spatial databases with noise, in: Pro. of 2nd Int. Conf. Knowledge Discovery and Data Mining, 1996, 1996, pp. 291–316.
  * [39] A. Rosenberg, J. Hirschberg, V-measure: A conditional entropy-based external cluster evaluation measure, in: Proceedings of the 2007 joint conference on empirical methods in natural language processing and computational natural language learning (EMNLP-CoNLL), 2007, pp. 410–420.
  * [40] L. Buitinck, G. Louppe, M. Blondel, F. Pedregosa, A. Mueller, O. Grisel, V. Niculae, P. Prettenhofer, A. Gramfort, J. Grobler, R. Layton, J. VanderPlas, A. Joly, B. Holt, G. Varoquaux, API design for machine learning software: experiences from the scikit-learn project, in: ECML PKDD Workshop: Languages for Data Mining and Machine Learning, 2013, pp. 108–122.
  * [41] H. Jun, R. Child, M. Chen, J. Schulman, A. Ramesh, A. Radford, I. Sutskever, Distribution augmentation for generative modeling, in: International Conference on Machine Learning, PMLR, 2020, pp. 5006–5019.
  * [42] D. N. Diniz, M. T. Rezende, A. G. Bianchi, C. M. Carneiro, D. M. Ushizima, F. N. de Medeiros, M. J. Souza, A hierarchical feature-based methodology to perform cervical cancer classification, Applied Sciences 11 (9) (2021) 4091.
  * [43] W. Liu, C. Li, N. Xu, T. Jiang, M. M. Rahaman, H. Sun, X. Wu, W. Hu, H. Chen, C. Sun, et al., Cvm-cervix: A hybrid cervical pap-smear image classification framework using cnn, visual transformer and multilayer perceptron, Pattern Recognition 130 (2022) 108829.
  * [44] J. Lin, Divergence measures based on the shannon entropy, IEEE Transactions on Information Theory 37 (1) (1991) 145–151.  doi:10.1109/18.61115.