# MDS-ViTNet: Improving saliency prediction for Eye-Tracking with Vision
Transformer

Polezhaev Ignat12 Goncharenko Igor12 Iurina Natalya1
1 _Sberbank, Moscow, Russia_
2 _MIPT, Dolgoprudny, Russia_
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we present a novel methodology we call MDS-ViTNet (Multi
Decoder Saliency by Vision Transformer Network) for enhancing visual saliency
prediction or eye-tracking. This approach holds significant potential for
diverse fields, including marketing, medicine, robotics, and retail. We
propose a network architecture that leverages the Vision Transformer, moving
beyond the conventional ImageNet backbone. The framework adopts an encoder-
decoder structure, with the encoder utilizing a Swin transformer to
efficiently embed most important features. This process involves a Transfer
Learning method, wherein layers from the Vision Transformer are converted by
the Encoder Transformer and seamlessly integrated into a CNN Decoder. This
methodology ensures minimal information loss from the original input image.
The decoder employs a multi-decoding technique, utilizing dual decoders to
generate two distinct attention maps. These maps are subsequently combined
into a singular output via an additional CNN model. Our trained model MDS-
ViTNet achieves state-of-the-art results across several benchmarks. Committed
to fostering further collaboration, we intend to make our code, models, and
datasets accessible to the public.

##  1 Introduction

The human visual system is remarkably complex. This paper concentrates on the
mechanisms underlying visual attention, predicting how individuals direct
their attention toward specific objects or regions within an image. There are
a number of studies in this field [1, 9-11, 22]. When a person looks at an
image, the foveal vision provides the highest resolution visual information.
Simultaneously, the peripheral vision, although less detailed, continues to
provide long-range visual insights. The eye movements recorded in the
experiments serve as the ground truth for estimating image saliency. This
paper aims to introduce a novel computer vision model that assigns a saliency
value to each pixel of an image, aligning more closely with the ground truth.

The task of predicting human saliency or visual attention maps has widespread
applications in numerous fields, such as marketing (for estimating the amount
of information a person assimilates while watching an advertisement),
robotics, and healthcare.

##  2 Related work

Currently, a variety of methods exist for predicting visual attention maps,
encompassing both neural network-based models and traditional approaches.
Traditional methods [7, 8, 12, 24] rely on low-level features like luminance,
color, texture, and contrast. In contrast, solutions involving neural networks
leverage more complex, higher-level patterns.

Many computer vision models, including eye-tracking, employ convolutional
neural networks (CNNs) [13, 14, 17, 21]. The success of visual saliency
prediction models can be attributed in part to the use of CNNs and to the
emerged of large datasets advancing the research in this field. However, CNN-
based models have their limitations. Each convolutional kernel in a CNN
receives information only from a local subset of image pixels, which restricts
these models’ ability to capture contextual information over a long distances.
Models based on Long Short-Term Memory (LSTM)[5, 18] networks are also
utilized. This approach has proven to be effective in processing local and
long-range visual information, thereby enhancing the accuracy of saliency
prediction.

To tackle the issues mentioned earlier, the TranSalNet model [20], which
integrates Transformers into a CNN-based architecture, was introduced in 2022.
Transformer encoders learn spatial dependencies across long distances via a
multi-head self-attention mechanism. This model has proven to be a state-of-
the-art solution on the largest open-source dataset, SALICON. TranSalNet
incorporates a backbone formed by convolutional ResNet/DenseNet networks. The
CNN generates three feature maps of different dimensions, which are inserted
into the transformer encoders followed by a single decoder, producing visual
salience maps.

After TranSalNet model, the value of transformers for the task of visual
saliency prediction became apparent. The initial discussions about visual
transformers emerged in 2020 [26], challenging the convolutional network
paradigm with a new architecture that surpassed its convolutional counterparts
in accuracy while using fewer parameters. Transformers made a significant
impact in the field of computer vision in 2021, with the introduction of the
ViT (Vision Transformer) [6].

Additionally, 2021 saw a shift in public focus towards the Swin Transformer
models. The Swin (Shifted Windows) Transformer [19], a versatile transformer
tailored for computer vision tasks (meaning it’s suitable for various
applications like semantic segmentation and image classification),
outperformed ViT-based and CNN-based architectures. The model authors made new
architecture more computationally optimal by utilizing localize attention
mechanism. The Shifted Window Multi-Head Attention maintained the network’s
representative capacity, enabling it to compete with contemporary models. This
approach facilitated the development of an architecture that extracts features
from images at different spatial scales, allowing Swin Transformer to excel as
a backbone in segmentation and detection tasks, where Transformers had
previously been less effective.

Figure 1: Example of model work on an image from SALICON

Contributions of this paper are as follows: (i) substitution of the CNN
backbone with visual Swin transformer in the saliency prediction model, (ii)
incorporation of multi-decoder architecture with additional CNN-based decoder
merging layer, thereby improving the accuracy of the model.

As a result of this work, multi-decoder saliency with visual transformer
network (MDS-ViTNet) is proposed and trained on existing datasets. This model
demonstrated excellent performance, establishing a new state-of-the-art (SoTA)
benchmark.

The code of the project can be found on the GitHub:
https://github.com/IgnatPolezhaev/MDS-ViTNet.

##  3 Model Architecture

###  3.1 The ViT encoder

Extensive research over the years has demonstrated the efficacy of utilizing
CNN networks for feature extraction in predicting attention maps. For
instance, TranSalNet [20] has achieved SoTA performance in this domain. In
this study, two models were developed, employing ResNet-50 and DenseNet-161 as
their backbones. Recognizing that TranSalNet employs CNNs as a backbone for
image element recognition and Transformer encoders for identifying elements
typically detected by human peripheral vision, we were motivated to
investigate the potential of Swin Transformers as an alternative backbone to
DenseNet and ResNet. This exploration yielded positive results.

In this paper, we introduce MDS-ViTNet (Multi Decoder Saliency by Vision
Transformer Network), a novel model that combines transformer components with
Swin transformers to capture long-range contextual visual information, a task
conventionally handled by CNNs. We utilized a Swin-T transformer from the
PyTorch library, with 28M parameters, as the model’s backbone. Our
experimental findings indicate that employing Swin transformers as the
backbone is key to capturing the maximum number of essential qualitative
features in images. When these are combined with other types of transformers
as encoders and convolutional decoders, they significantly contribute to the
model’s saliency prediction capabilities, thereby enhancing its perceptual
relevance and performance. The unique approach of our proposed model, which
leverages a combination of transformers for attention prediction, successfully
surpasses the SoTA benchmark in competitions like Salicon.

###  3.2 Transformers in the encoder architecture

Transformers were initially introduced in the context of natural language
processing (NLP) [23]. Since then, this architecture has found widespread
application in various domains, including time series analysis [25], computer
vision [6], and numerous other fields. The popularity of transformers can be
attributed to their ability to handle large contexts. Their capacity to retain
substantial amounts of information in memory is enabled by the self-attention
mechanism. Research has demonstrated their potent capability for capturing
long-range information, a feature that can significantly enhance gaze
prediction tasks.

Transformers have also gained considerable traction in the field of computer
vision. They exhibit competitive performance alongside convolutional networks,
which have long been dominant in this area. Currently, there is a diverse
array of architectures based on visual Transformers, including DeepViT [28],
Swin Transformer [19], and Swin Transformer V2. These models have attained
significant success in various visual tasks, predominantly due to the
advantages offered by Transformer technology.

In our MDS-ViTNet model, we have implemented similar approaches as in the
field of computer vision with Transformers. Our model has a backbone that
produces six sets of feature maps of different scales. These feature maps are
inserted into the encoder to amplify range and contextual information. The
architecture is illustrated in Figure 2. Within this structure, the feature
sets $x_{2}$, $x_{4}$, $x_{6}$ with spatial dimensions of (w/8, h/8), (w/16,
h/16), (w/32, h/32) are channeled to a first decoder. Meanwhile, the feature
sets $x_{1}$, $x_{3}$, $x_{5}$ with spatial dimensions of (w/4, h/4), (w/8,
h/8), (w/16, h/16) are channeled to a second decoder. More precisely, the
feature sets $x_{1}$ and $x_{2}$ undergo a reduction in their third dimension
to 512, and the sets $x_{3}$ to $x_{6}$ are similarly reduced to 768.

There is no relative or absolute position information in the feature maps, it
is necessary to use position embedding (POS) to provide position understanding
before feeding input data to the transform encoders. Therefore, as in the work
of the TranSalNet colleagues [20], an absolute POS [6] is applied before
feeding the input data into the transformer encoders to form a piecewise
complement of the input data and a trained matrix of the same shape as the
input data. Each transformer encoder contains two identical layers of standard
multi-headed self-attention (MSA) and multi-layer perceptron (MLP) blocks
[link-link (google 16x16)]. The model employs 8-headed attention in transform
coders 1, 2 and 12-headed attention in coders 3 - 6. Layer normalization (LN)
and residual coupling are applied before and after each block, respectively.
The processing in each transform coder can be represented as follows:

In the absence of relative or absolute positional information within the
feature maps, the incorporation of position embedding (POS) is crucial to
incorporate positional understanding prior to the input data being fed into
the transformer encoders. Hence, akin to the approach used by authors of [20],
we apply an absolute POS [6] before inputting the data into the transformer
encoders. This step involves forming a composite of the input data with a
trained matrix, which matches the shape of the input data.

Each transformer encoder within our model comprises two identical layers,
consisting of standard multi-headed self-attention (MSA) and multi-layer
perceptron (MLP) blocks [refer to link-link (google 16x16)]. The model is
designed to use 8-headed attention in transformer encoders 1 and 2, and
12-headed attention in encoders 3 to 6. Layer normalization (LN) and residual
connections are implemented both before and after each block, respectively.

The processing sequence in each transformer encoder can be described as
follows:

$\large z_{0}=Conv_{1\times 1}(x_{i})\oplus POS_{i},i=1,2,3$ (1)

$\large z_{l}\textquotesingle=MSA(LN(z_{l-1}))\oplus z_{l-1},l=1,2$ (2)

$\large z_{l}=MSA(LN(z_{l}\textquotesingle))\oplus
z_{l}\textquotesingle,l=1,2$ (3)

where $z_{l}$ are the output feature maps of the l-th layer in the transformer
encoder, and $x_{i}$ are input feature maps from the ViT encoder.

###  3.3 Multi-decoder architecture

Multidecoder architecture (MDA) model in deep learning refers to a mechanism
that allows the use of multiple decoders to solve a problem. The advantage of
such an architecture is to increase performance and provide faster data
processing speed since independent decoders can work in parallel. Also MDA
contributes to an error reduction, i.e. if one decoder cannot process the
input data correctly, another decoder is able to correct this error. Thus, it
is possible to increase the number of parameters by adding new decoders in
parallel, which will not increase the inference time as much as adding
consecutive layers. MDA is used in many deep learning tasks. For example, in
speech technology [29], due to its multiple decoder and counting head
architecture, the model requires one forward pass at test time on a single
network. MDA has also been mentioned in work on diffusion models [27]. This
work breaks the time interval into multiple stages, where specialized multi-
decoder U-network architectures that combine time-dependent models with a
universal common encoder are used. This approach allows efficient allocation
of computational resources and elimination of inter-stage interference, which
significantly improves the learning efficiency.

In our work, we propose two independent CNN decoders with one common backbone
encoder. The set of features from the transformer-coders (see Fig. 2) fall
into the decoders, which predict based on the attention map. Decoder 1 accepts
maps $x_{2}$, $x_{4}$, and $x_{6}$, and decoder 2 $x_{1}$, $x_{3}$, and
$x_{5}$. The decoders combine these feature maps and restore to the original
resolution of the image. Each decoder uses 7 convolution layers. The $x_{6}$
and $x_{5}$ maps are inputted into decoders 1 and 2, respectively. After the
first and second CNN layers, the significance maps $x_{4}$ ($x_{3}$ for
decoder 2) and $x_{2}$ ($x_{1}$ for decoder 2) are element-wise multiplied
with the preceding feature maps. Batch normalization (BN) and activation
function (ReLU for all convolution layers except the last one; Sigmoid for the
last convolution layer) are applied after each $3\times 3$ convolution
operation ($Conv_{3\times 3}$), where the former is used to promote
convergence and the latter to increase the non-linearity coefficient of the
model.

Figure 2: MDS-ViTNet architecture

###  3.4 An additional model for combining saliency maps

The multi-decoder architecture in our solution predicts two saliency maps. It
was proposed to write an additional CNN-Merge model on convolutional layers
that combines these maps and predicts one. CNN-Merge takes as input $x_{1}$
and $x_{2}$, where $x_{1}$ and $x_{2}$ are two independent importance maps
obtained from the decoder. These maps are concatenated by channel
dimensionality and processed by convolutional layers. The architecture
consists of 7 convolutional layers $Conv_{3\times 3}$. Batch normalization
(BN) and activation function (ReLU) are used after every convolutional layer
except the last one. The Sigmoid activation function comes after the 7th
layer.

###  3.5 Training losses

The question of how to evaluate a model’s ability to predict visual attention
map of an image remains an open question. There are many metrics that compare
the two saliency maps. The paper [3] gives examples of the latter, and has
done a lot of work investigating them. Research in the area of predicting
visual saliency maps [4, 15] has shown that using saliency score metrics to
define a loss function can significantly improve the model performance.

In our work, following research by TranSalNet colleagues [20], we used a
linear combination of different loss functions: Similarity (SIM), Pearson’s
Correlation Coefficient (CC), Kullback-Leibler divergence (KL). These metrics
can be categorized into those that compare the similarity between two
distributions and those based on the pixel-by-pixel difference. KL is the
metric that baed on thethe pixel-by-pixel difference and the others (SIM and
CC) the similarity of the distribution, so in the aggregate loss function KL
is taken with a plus sign and the others with a minus sign. In the TranSalNet
work, a study was conducted on what coefficients to take in a linear
combination of metrics. The minimum loss values on the validation set were
fixed, and the coefficients on the loss functions are assigned so that their
contribution to the total loss is relatively equal. As a result, weights of 1,
2, and 10 were assigned to SIM, CC, and KL, respectively.

$\large L_{CC}(\hat{m},m)=\frac{cov(\hat{m},m)}{\sigma(\hat{m})\sigma(m)}$ (4)

$\large L_{SIM}(\hat{m},m)=\sum_{i}min(\hat{m}_{i},m_{i})$ (5)

$\large
L_{KL}(\hat{m},m)=\sum_{i}m_{i}\log(\varepsilon+\frac{m_{i}}{\varepsilon+\hat{m}_{i}})$
(6)

where $\hat{m}$, $m$ \- predicted and true saliency maps,
$\varepsilon=2.2204\times 10^{-16}$ \- regularization constant.

###  3.6 Evaluation metrics

In evaluating the quality of visual attention map prediction, various metrics
are employed, broadly classified into two categories based on their alignment
with human judgments. The first category, ’perception-based metrics’, includes
measures like CC (Pearson’s Correlation Coefficient) and SIM (Similarity). The
second category, ’non-perception-based metrics’, comprises AUC (Area Under the
Curve) and KLD (Kullback-Leibler Divergence). It’s important to note that non-
perception-based metrics, while not directly measuring perception, are still
valuable as they can detect specific gaze behavior patterns, such as
topography. Metrics like AUC and KL are particularly suited for detection
tasks as they penalize incorrectly identified objects. Others, like SIM, are
adept at estimating the relative importance of different regions within an
image. In the context of eye-tracking, the AUC metric is commonly used to
compare prediction models.

In our study, we developed a custom metric to assess the prediction quality of
visual importance maps using AUC. This metric involves selecting a range of
threshold values from 0 to 1, corresponding to the range of values in the
attention map. For each threshold, True Positives (TP), True Negatives (TN),
False Positives (FP), and False Negatives (FN) are calculated point by point.
Then, the True Positive Rate (TPR) and False Positive Rate (FPR) for each
threshold are computed, leading to the calculation of the area under the
curve. This process enables the metric to evaluate the model’s predictive
ability in terms of accurately identifying true position fixations.

Finally, we compare MDS-ViTNet model with TranSalNet using the metrics SIM,
CC, KL, and AUC, providing an evaluation of its performance in predicting
saliency maps.

##  4 Datasets and data

###  4.1 SALICON and CAT2000 dataset

In this work, two public datasets, SALICON and CAT2000, were used to train the
model. The SALICON [16] is a large dataset on the popular Microsoft Common
Objects in Context (MS COCO) image database. MS COCO is a large-scale image
dataset that emphasizes non-iconic images and objects in context. It provides
a rich set of specific annotations for image recognition, segmentation and
captioning. The SALICON dataset was collected in 2017. It contains 10,000
training images and 5,000 validation images. A set of 5,000 test images
without attention maps is also provided.

The CAT2000 dataset [2] contains 2,000 images from 20 different categories
with eye tracking data from 24 observers. The entire dataset was not used in
the paper, only part of it. This is because some categories are very different
from the main distribution of the SALICON dataset, i.e. there were sketches,
drawings or monochrome images.

###  4.2 Augmentation

Augmentation was used in this paper to extend the dataset. The augmentation
was chosen so that the images do not differ from what a person sees in
reality. That is, the image should not be inverted, monochrome or large
changes in color. Otherwise, with large and irregular changes, a person’s
attention will be distorted, for example, in reality some objects cannot be
upside down.

Thus, in this paper, augmentation in the form of mirrored horizontal display,
small Gaussians blur and changes in contrast, brightness sharpness and
saturation were chosen. The images are further normalized and converted to a
size of (288, 384).

##  5 Experiments

###  5.1 Learning a multi-decoder model

The model was trained for 15 epochs on an A100 80Gb graphics card. The number
of model parameters is 111.5M. 4 images were taken as the batch size. The
8-bit Adam optimizer of bitsandbytes library was used in training. This
optimizer was chosen to reduce the video memory consumption in the training in
order to increase the size of the batches. StepLR scheduler was also used to
reduce the learning rate. There was a stopping criterion in training if 5
epochs in a row loess on validation does not improve the value. The best
quality was achieved at the 9th epoch. The results of the model comparison are
summarized in the table 1.

###  5.2 Learning a CNNMerge model

The CNN-Merge model was also trained for 15 epochs on an A100 80Gb graphics
card. The number of convolution model parameters is about 400K. The size of
the batches is 64. AdamW optimizer with learning rate parameter 1e-5 c
scheduler StepLR with learning rate halving every 10 epochs. The CNN-Merge
training also utilized the stopping criterion. The best quality was achieved
at epoch 30. The results of the model comparison are summarized in the table
1.

Model | AUC $\nearrow$ | KL $\searrow$ | CC $\nearrow$ | SIM $\nearrow$
---|---|---|---|---
MDS-ViTNet | 0.8684 | 0.2171 | 0.8949 | 0.7855
MDS-ViTNet decoder 1 | 0.8650 | 0.2215 | 0.8911 | 0.7817
MDS-ViTNet decoder 2 | 0.8676 | 0.2127 | 0.8980 | 0.7887
TranSalNet DenseNet | 0.8486 | 0.2218 | 0.8875 | 0.7757
Table 1: Table of model comparison on validation sample

##  6 Conclusion

The conducted research successfully addressed the eye-tracking problem. A
thorough analysis of TranSalNet model was undertaken. This led to a
comprehensive training cycle for our model, which included compiling an
extensive dataset and implementing several imrpovements of the architecture. A
key decision was the adoption of the Swin Transformer as an encoder, which has
demonstrated superior performance in computer vision tasks compared to
traditional encoders based on convolutional layers. Additionally, a multi-
decoder architecture was employed, wherein each decoder independently predicts
saliency maps. These maps are further processed by an auxiliary CNN model.
This additional processing step contributed to an improvement in the quality
of the final output. Finally, this approach enabled us to achieve a state-of-
the-art result.

## Acknowledgments

We are grateful to Mark Potanin for his fruitful comments and suggestions. We
also thank our colleagues Pavel Abramov and Andrey Mamontov for helpful
discussions.

## References

  * (1) Ali Borji  “Saliency Prediction in the Deep Learning Era: Successes and Limitations”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 43, 2019, pp. 679–700
  * (2) Ali Borji and Laurent Itti  “CAT2000: A Large Scale Fixation Dataset for Boosting Saliency Research”  In _ArXiv_ abs/1505.03581, 2015
  * (3) Zoya Bylinskii et al.  “What Do Different Evaluation Metrics Tell Us About Saliency Models?”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 41, 2016, pp. 740–757
  * (4) Zhaohui Che et al.  “How is Gaze Influenced by Image Transformations? Dataset and Model”  In _IEEE Transactions on Image Processing_ 29, 2019, pp. 2287–2300
  * (5) Marcella Cornia, Lorenzo Baraldi, Giuseppe Serra and Rita Cucchiara  “Predicting Human Eye Fixations via an LSTM-Based Saliency Attentive Model”  In _IEEE Transactions on Image Processing_ 27, 2016, pp. 5142–5154
  * (6) Alexey Dosovitskiy et al.  “An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale”  In _ArXiv_ abs/2010.11929, 2020
  * (7) Erkut Erdem and Aykut Erdem  “Visual saliency estimation by nonlinearly integrating features using region covariances.”  In _Journal of vision_ 13 4, 2013, pp. 11
  * (8) Shu Fang et al.  “Learning Discriminative Subspaces on Random Contrasts for Image Saliency Analysis”  In _IEEE Transactions on Neural Networks and Learning Systems_ 28, 2017, pp. 1095–1108
  * (9) Ke Gu et al.  “Automatic Contrast Enhancement Technology With Saliency Preservation”  In _IEEE Transactions on Circuits and Systems for Video Technology_ 25, 2015, pp. 1480–1494
  * (10) Jungong Han, Eric J. Pauwels and Paul M. Zeeuw  “Fast saliency-aware multi-modality image fusion”  In _Neurocomputing_ 111, 2013, pp. 70–80
  * (11) Yiyuan Han, Bing Han and Xinbo Gao  “Human scanpath estimation based on semantic segmentation guided by common eye fixation behaviors”  In _Neurocomputing_ 453, 2020, pp. 705–717
  * (12) Jonathan Harel, Christof Koch and Pietro Perona  “Graph-Based Visual Saliency”  In _Neural Information Processing Systems_ , 2006
  * (13) Kaiming He, X. Zhang, Shaoqing Ren and Jian Sun  “Deep Residual Learning for Image Recognition”  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015, pp. 770–778
  * (14) Gao Huang, Zhuang Liu and Kilian Q. Weinberger  “Densely Connected Convolutional Networks”  In _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 2261–2269
  * (15) Sen Jia  “EML-NET: An Expandable Multi-Layer NETwork for Saliency Prediction”  In _ArXiv_ abs/1805.01047, 2018
  * (16) Ming Jiang, Shengsheng Huang, Juanyong Duan and Qi Zhao  “SALICON: Saliency in Context”  In _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015, pp. 1072–1080
  * (17) Alex Krizhevsky, Ilya Sutskever and Geoffrey E. Hinton  “ImageNet classification with deep convolutional neural networks”  In _Communications of the ACM_ 60, 2012, pp. 84–90
  * (18) Nian Liu and Junwei Han  “A Deep Spatial Contextual Long-Term Recurrent Convolutional Network for Saliency Detection”  In _IEEE Transactions on Image Processing_ 27, 2016, pp. 3264–3274
  * (19) Ze Liu et al.  “Swin Transformer: Hierarchical Vision Transformer using Shifted Windows”  In _2021 IEEE/CVF International Conference on Computer Vision (ICCV)_ , 2021, pp. 9992–10002
  * (20) Jianxun Lou et al.  “TranSalNet: Towards perceptually relevant visual saliency prediction”  In _Neurocomputing_ 494, 2021, pp. 455–467
  * (21) Karen Simonyan and Andrew Zisserman  “Very Deep Convolutional Networks for Large-Scale Image Recognition”  In _CoRR_ abs/1409.1556, 2014
  * (22) Kaikai Song, Ting Yao, Qiang Ling and Tao Mei  “Boosting image sentiment analysis with visual attention”  In _Neurocomputing_ 312, 2018, pp. 218–228
  * (23) Ashish Vaswani et al.  “Attention is All you Need”  In _Neural Information Processing Systems_ , 2017
  * (24) Dirk B. Walther and Christof Koch  “Modeling attention to salient proto-objects”  In _Neural networks : the official journal of the International Neural Network Society_ 19 9, 2006, pp. 1395–407
  * (25) Qingsong Wen et al.  “Transformers in Time Series: A Survey”  In _International Joint Conference on Artificial Intelligence_ , 2022
  * (26) Bichen Wu et al.  “Visual Transformers: Token-based Image Representation and Processing for Computer Vision”  In _ArXiv_ abs/2006.03677, 2020
  * (27) Huijie Zhang et al.  “Improving Efficiency of Diffusion Models via Multi-Stage Framework and Tailored Multi-Decoder Architectures”  In _ArXiv_ abs/2312.09181, 2023
  * (28) Daquan Zhou et al.  “DeepViT: Towards Deeper Vision Transformer”  In _ArXiv_ abs/2103.11886, 2021
  * (29) Junzhe Zhu, Raymond A. Yeh and Mark A. Hasegawa-Johnson  “Multi-Decoder DPRNN: High Accuracy Source Counting and Separation”  In _ArXiv_ abs/2011.12022, 2020


# An exact test for significance of clusters in binary data

James Mathews, Cameron Crowe, Rami Vanguri *,
Margaret Callahan*, Travis Hollmann*, Saad Nadeem* Memorial Sloan-Kettering
Cancer Center

###### Abstract

Unsupervised clustering of feature matrix data is an indispensible technique
for exploratory data analysis and quality control of experimental data.
However, clusters are difficult to assess for statistical significance in an
objective way. We prove a formula for the distribution of the size of the set
of samples, out of a population of fixed size, which display a given
signature, conditional on the marginals (frequencies) of each individual
feature comprising the signature. The resulting “exact test for coincidence”
is widely applicable to objective assessment of clusters in any binary data.
We also present a software package implementing the test, a suite of
computational verifications of the main theorems, and a supplemental tool for
cluster discovery using Formal Concept Analysis.

###### Contents

  1. 1 Introduction
  2. 2 Theory
    1. 2.1 Setup
    2. 2.2 Binary matrix configurations
    3. 2.3 Incidence statistic, its PMF, and CDF
    4. 2.4 CDF generating function and incomplete beta function
  3. 3 Software implementation
    1. 3.1 Python package
    2. 3.2 Command-line tool
    3. 3.3 Web application
    4. 3.4 Testing
  4. 4 Related work
  5. A Formal Concept Analysis bicluster identification
  6. B Figures

##  1 Introduction

A typical visualization of a binary data matrix is a hierarchically-clustered
heatmap, with dendrograms in which the higher-level clusters are recursively
comprised of smaller clusters, the hierarchy being computed with an
agglomeration strategy involving a distance function defined pairwise between
samples (or features) to be clustered. In favorable cases a cluster may appear
at some level of the hierachy which is especially characteristic of an
important underlying state or measure, i.e. an outcome. For example,
likelihood of favorable response to some medical treatment.

But it is often difficult to decide whether a cluster found this way, or any
other way, could just as easily have occurred by random chance. This is
obviously a primary concern in the unsupervised context, where outcomes which
might guide cluster assessment are not present. It is also a concern in the
supervised context, due to the possibility of overfitting or multiple-
hypothesis false discovery.

As an example, in recent work of the authors[9], a subtype of several types of
cancers (including lung and uterus cancers) was identified which exhibited a
molecular signature defined by about 10 genes, the PSGs. Network analysis
methods implicated the gene subset, but initially confidence concerning its
actual significance was low. Pearson correlation analysis was inconclusive due
to the presence of outliers. The rarity of the subtype displaying the full
signature added to this uncertainty. Ultimately Kaplan-Meier analysis did show
that the PSG+ phenotype confers a poor prognosis, confirming the biological
significance of this subtype, but we still lacked an objective basis for any
claim of statistical significance of the signature/subtype itself. The exact
test we introduce in this article turns out to provide such a basis, as
described in Figure 1.

##  2 Theory

###  2.1 Setup

Let $\overline{M}$ be a binary matrix of shape $(N,K)$. We call the $K$
columns _features_ and the $N$ rows _samples_. Given a $k$-element subset $F$
of the feature set, let $M$ denote the restriction of $\overline{M}$ to the
columns $F$, and let $v:=(v_{1},\ldots,v_{k})$ denote the corresponding column
sums. Let $S$ denote the set of samples which have all of the features $F$.
That is,

$\displaystyle S=\\{\thinspace s\quad|\quad M(s,f)=1\quad\forall f\in F\\}$

A pair $(F,S)$ obtained as above may be called a _maximal bicluster_ or a
_formal concept_. We shall use the term _signature_ , emphasizing the feature
set $F$, and call $S$ the set of samples _displaying signature_ $F$. In
appendix A we explain how to identify, in practice, many examples of $(F,S)$
for which $S$ is non-empty and relatively large. Of course, any other
signature discovery method may be used instead.

We propose to assess the significance of a given signature finding in terms of
the size $|S|$, under the intuition that simultaneous display of multiple
features by a large set of samples indicates a non-trivial relation between
the features. We call this size the _incidence_ or _intersection_ statistic,
and denote it $I$.

###  2.2 Binary matrix configurations

We are concerned with binary matrices $M$. If $M$ has $k$ columns, it will be
convenient to do some calculations in the ring of formal power series
$T:=\mathbb{Z}[[t_{1},\ldots,t_{k}]]$. This is because of the correspondence
between:

  1. 1.

Multiplicity-free monomials in $T$, i.e. elements of the form
$t_{J}:=\prod\limits_{j\in J}t_{j}$ for some $J\subset\\{1,\ldots,k\\}$

  2. 2.

Subsets $J\subset\\{1,\ldots,k\\}\thinspace$ (i.e. $J\in\mathscr{P}_{k}$)

  3. 3.

Possible rows $r=(r^{1},\ldots,r^{k})$ of $M$

The correspondence is

$\displaystyle t_{J}\quad\longleftrightarrow\quad
J\quad\longleftrightarrow\quad r=(r^{1},\ldots r^{k}),\thinspace
r^{j}=\begin{cases}1\mbox{ if }j\in J\\\ 0\mbox{ if }j\notin J\end{cases}$

Denote by $\mathscr{F}(k)$ the set defined by any of these 3 equivalent
descriptions. (Here $\mathscr{F}$ stands for ”features”.)

Symmetrically, if $M$ has $n$ rows, we consider the ring
$W:=\mathbb{Z}[[s_{1},\ldots,s_{n}]]$, and the 3 sets in correspondence:

  1. 1.

Multiplicity-free monomials in $W$, i.e. elements of the form
$s_{U}:=\prod\limits_{u\in U}s_{u}$ for some $U\subset\\{1,\ldots,n\\}$

  2. 2.

Subsets $U\subset\\{1,\ldots,n\\}\thinspace$ (i.e. $U\in\mathscr{P}_{n}$)

  3. 3.

Possible columns $c=(c^{1},\ldots,c^{n})$ of $M$

Denote this set by $\mathscr{S}(n)$. (Here $\mathscr{S}$ stands for
”samples”).

In these terms, the set of all $M$ is naturally identified with
$(\mathscr{F}(k))^{n}$ and with $(\mathscr{S}(n))^{k}$ by regarding $M$ as an
$n$-tuple of rows or, respectively, as a $k$-tuple of columns.

We will also call the matrices $M$ _configurations_ , writing

$(\mathscr{F}(k))^{n}\cong(\mathscr{S}(n))^{k}=:\mathscr{C}$
$(J_{1},\ldots,J_{n})\quad\longleftrightarrow\quad(U_{1},\ldots,U_{k})\quad\longleftrightarrow
M$

In counting configurations satisfying certain conditions, we will appeal to
the notation introduced above for corresponding elements in lieu of explicit
notation for the bijection functions.

###  2.3 Incidence statistic, its PMF, and CDF

Define integers $a(n,v)$, for integers $n\geq 0$ and $v=(v_{1},\ldots,v_{k})$
with $v_{j}\geq 0\thinspace\forall j$, by the generating function:

$\displaystyle(f(t)-t_{1}\cdots t_{k})^{n}$
$\displaystyle=:\sum_{v}a(n,v)t_{1}^{v_{1}}\cdots t_{k}^{v_{k}}$
$\displaystyle f(t):$ $\displaystyle=(1+t_{1})\cdots(1+t_{k})$

The following counting theorem is the underlying fact needed to prove a
formula for the probability mass function (PMF) of the incidence statistic.

######  Theorem 1.

  1. 1.

$a(n,v)$ is the number of configurations in which the mutual intersection of
the $U_{j}$ is empty, that is $\cap_{j=1}^{j=k}U_{j}=\emptyset$, and such that
$|U_{j}|=v_{j}$ for each $j$.

  2. 2.

$a(n,v)=\sum\limits_{m=0}^{m=n}(-1)^{n+m}\binom{n}{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$

###### Proof.

(1) By expansion, $f(t)$ consists of the sum of all the monomials in
$\mathscr{F}(k)$. So $f(t)-t_{1}\cdots t_{k}$ is the sum of all the monomials
except $t_{1}\cdots t_{k}$. Before collecting terms with the same monomial
part, the terms of $(f(t)-t_{1}\cdots t_{k})^{n}$ are labelled by ordered
$n$-tuples of elements of $\mathscr{F}(k)\backslash\\{t_{1}\cdots t_{k}\\}$.
That is, by certain elements of $\mathscr{C}$. Thus the notation we have
introduced for elements of $\mathscr{C}$ may be brought to bear. In
particular, the monomial part of a given term is

$t_{1}^{|U_{1}|}\cdots t_{k}^{|U_{k}|}$

It follows that the coefficient of $t_{1}^{v_{1}}\cdots t_{k}^{v_{k}}$ is the
number of configurations, in which no $J_{i}$ equals the whole set
$\\{1,\ldots,k\\}$ (due to the missing element $t_{1}\ldots t_{k}$), such that
$|U_{j}|=v_{j}$ for all $j$. The condition that no $J_{i}$ be equal to the
whole set is equivalent to the mutual intersection of $U_{j}$ being empty.

(2) We apply the binomial theorem $1+k$ times:

$\displaystyle(f(t)-t_{1}\cdots t_{k})^{n}$
$\displaystyle=\sum\limits_{m=0}^{m=n}(-1)^{n-m}\binom{n}{m}(f(t))^{m}(t_{1}^{n-m}\cdots
t_{k}^{n-m})$
$\displaystyle=\sum\limits_{m=0}^{m=n}(-1)^{n+m}\binom{n}{m}(1+t_{1})^{m}\cdots(1+t_{k})^{m}(t_{1}^{n-m}\cdots
t_{k}^{n-m})$
$\displaystyle=\sum\limits_{m=0}^{m=n}(-1)^{n+m}\binom{n}{m}\left(\sum\limits_{u=0}^{u=m}\binom{m}{u}t_{1}^{u}\right)\cdots\left(\sum\limits_{u=0}^{u=m}\binom{m}{u}t_{k}^{u}\right)(t_{1}^{n-m}\cdots
t_{k}^{n-m})$
$\displaystyle=\sum\limits_{m=0}^{m=n}(-1)^{n+m}\binom{n}{m}\left(\sum\limits_{v}\prod\limits_{j=1}^{j=k}\binom{m}{v_{j}}t_{1}^{v_{1}}\cdots
t_{k}^{v_{k}}\right)(t_{1}^{n-m}\cdots t_{k}^{n-m})$
$\displaystyle=\sum\limits_{v}\sum\limits_{m=0}^{m=n}(-1)^{n+m}\binom{n}{m}\prod\limits_{j=1}^{j=k}\binom{m}{v_{j}}t_{1}^{n-m+v_{1}}\cdots
t_{k}^{n-m+v_{k}}$
$\displaystyle=\sum\limits_{v}\sum\limits_{m=0}^{m=n}(-1)^{n+m}\binom{n}{m}\prod\limits_{j=1}^{j=k}\binom{m}{v_{j}-(n-m)}t_{1}^{v_{1}}\cdots
t_{k}^{v_{k}}$
$\displaystyle=\sum\limits_{v}\sum\limits_{m=0}^{m=n}(-1)^{n+m}\binom{n}{m}\prod\limits_{j=1}^{j=k}\binom{m}{n-v_{j}}t_{1}^{v_{1}}\cdots
t_{k}^{v_{k}}$

∎

The proof above is clarified somewhat by the observation that $f(t)$ can be
expressed as a specialization of the power series in $W\otimes_{\mathbb{Z}}T$,

$g(s,t):=\prod\limits_{u,j}(1+s_{u}t_{j})\quad,$

namely $f(t)=g(\mathbf{1},t)$ where $\mathbf{1}=(1,\ldots,1)$.

######  Theorem 2.

Fix integers $i\geq 0$, $v=(v_{1},\ldots,v_{k})$, $v_{j}\geq 0$, and $n>0$.
Consider the $n\times k$ configurations $M$ in which:

  1. 1.

$|U_{j}|=v_{j}$ for each $j$.

  2. 2.

The cardinality of the intersection of the $U_{j}$ is exactly $i$, that is
$|\cap_{j=1}^{j=k}U_{j}|=i$.

The number of such configurations is given by the formula:

$\binom{n}{i}\sum\limits_{m=0}^{m=n-i}(-1)^{n-i+m}\binom{n-i}{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$

###### Proof.

The indicated set of configurations is partitioned equally into $\binom{n}{i}$
sets, according to which $i$-element sample subset is the mutual intersection,
denoted $X$. By construction the reduced configuration matrix, not involving
the elements of $X$, must consist of $k$ features with sample sets of sizes
$(v_{1}-i,\ldots,v_{k}-i)$ and with no intersection. Thus the size of each
part of the partition is $a(n-i,(v_{1}-i,\ldots,v_{k}-i))$. The number of
configurations is therefore

$\binom{n}{i}a(n-i,(v_{1}-i,\ldots,v_{k}-i))$

The result follows from the formula for $a$ given in Theorem 1.2. ∎

The null assumption we make for our test is the one that is made implicitly in
a standard permutation test, namely the uniform distribution on the subset of
$\mathscr{C}$ defined by $|U_{j}|=v_{j}$, given $v=(v_{1},\ldots,v_{k})$. Note
that this entails that we do _not_ assume $M$ is comprised of $n$ independent
and identically distributed (iid) samples. Also, despite the fact that $M$
appears to be $n$ samples from a set of binary discrete variables, it is
definitely not $n$ samples of Bernoulli variables; for example, the variance
of the number of positives is 0 for each feature, rather than $np(1-p)$ for
some positivity rate $p$.

Under this assumption the incidence statistic $I$ is an integer-valued random
variable. The following corollary provides a formula for its PMF.

######  Corollary 3.

Consider $n$ samples observed with $k$ binary features of respective
frequencies $v_{1},\ldots v_{k}$. The probability of observing exactly $i$
samples positive for all $k$ features is:

$p(I=i)=\frac{\binom{n}{i}\sum\limits_{m=0}^{m=n-i}(-1)^{n-i+m}\binom{n-i}{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}}{\prod\limits_{j=1}^{j=k}\binom{n}{v_{j}}}$

By summing over several values of $i$ in Corollary 3, one can compute a value
of the cumulative distribution function (CDF) of $I$. This is (one minus) the
$p$-value for the proposed “exact test for coincidence”.

The next theorem provides an alternative, more closed-form calculation of the
CDF, with significantly decreased computational complexity compared with
direct summation of PMF values, namely $O(n)$ rather than $O(n^{2})$.

The proof of this theorem depends on two basic lemmas.

######  Lemma 4.

$\binom{a}{b}\binom{b}{c}=\binom{a-c}{a-b}\binom{a}{c}$

###### Proof.

$\displaystyle\frac{a!}{(a-b)!b!}\cdot\frac{b!}{(b-c)!c!}=\frac{1}{(a-b)!(b-c)!}\cdot\frac{a!}{c!}=\frac{(a-c)!}{(a-b)!(b-c)!}\cdot\frac{a!}{(a-c)!c!}$

∎

######  Lemma 5.

$\sum\limits_{h=0}^{h=l}(-1)^{h}\binom{g}{h}=(-1)^{l}\binom{g-1}{l}$

###### Proof.

By induction. Base case $g=1$:

$\displaystyle(-1)^{0}\binom{1}{0}$ $\displaystyle=1=(-1)^{0}\binom{0}{0}$
$\displaystyle\binom{1}{0}-\binom{1}{1}$
$\displaystyle=0=(-1)^{1}\binom{0}{1}$

Now assume the formula holds (for all $l$) for a fixed $g\geq 0$.

$\displaystyle\sum\limits_{h=0}^{h=l}(-1)^{h}\binom{g+1}{h}$
$\displaystyle=\sum\limits_{h=0}^{h=l}(-1)^{h}\left(\binom{g}{h}+\binom{g}{h-1}\right)$
$\displaystyle=\sum\limits_{h=0}^{h=l}(-1)^{h}\binom{g}{h}+\sum\limits_{h=0}^{h=l}(-1)^{h}\binom{g}{h-1}$
$\displaystyle=\sum\limits_{h=0}^{h=l}(-1)^{h}\binom{g}{h}+\sum\limits_{h=1}^{h=l}(-1)^{h}\binom{g}{h-1}$
$\displaystyle=\sum\limits_{h=0}^{h=l}(-1)^{h}\binom{g}{h}-\sum\limits_{h=0}^{h=l-1}(-1)^{h}\binom{g}{h}$
$\displaystyle=(-1)^{l}\binom{g}{l}$

∎

######  Theorem 6.

$\displaystyle\sum\limits_{u=i}^{u=n}p(I=u)=\frac{N}{D}$

where

$\displaystyle
N:=\sum\limits_{m=\text{\emph{max}}\\{n-v_{j}\\}}^{m=n-i}(-1)^{m}\binom{n}{m}\left((-1)^{\text{\emph{max}}\\{n-v_{j}\\}}\binom{n-m-1}{n-\text{\emph{max}}\\{n-v_{j}\\}}+(-1)^{n-i}\binom{n-m-1}{i-1}\right)\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$
$\displaystyle D:=\prod\limits_{j=1}^{j=k}\binom{n}{v_{j}}$

###### Proof.

First note that $p(I=u)=0$ if $u>\text{min}\\{v_{j}\\}$, so the sum stops at
$u=\text{min}\\{v_{j}\\}$. We apply the formula for $p(I=u)$:

$\displaystyle\sum\limits_{u=i}^{u=\text{min}\\{v_{j}\\}}\binom{n}{u}a(n-u,(v_{1}-u,\ldots,v_{k}-u))=\sum\limits_{u=i}^{u=\text{min}\\{v_{j}\\}}\binom{n}{u}\sum\limits_{m=0}^{m=n-u}(-1)^{n-u+m}\binom{n-u}{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$
$\displaystyle=(-1)^{n}\sum\limits_{m=0}^{m=\infty}(-1)^{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}\sum\limits_{u=i}^{u=\text{min}\\{v_{j}\\}}(-1)^{u}\binom{n}{n-u}\binom{n-u}{m}$
$\displaystyle=(-1)^{n}\sum\limits_{m=0}^{m=\infty}(-1)^{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}\sum\limits_{u=i}^{u=\text{min}\\{v_{j}\\}}(-1)^{u}\binom{n-m}{u}\binom{n}{m}$
$\displaystyle=(-1)^{n}\sum\limits_{m=0}^{m=\infty}(-1)^{m}\binom{n}{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}\sum\limits_{u=i}^{u=\text{min}\\{v_{j}\\}}(-1)^{u}\binom{n-m}{u}$
$\displaystyle=(-1)^{n}\sum\limits_{m=0}^{m=\infty}(-1)^{m}\binom{n}{m}\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}\left((-1)^{\text{min}\\{v_{j}\\}}\binom{n-m-1}{\text{min}\\{v_{j}\\}}-(-1)^{i-1}\binom{n-m-1}{i-1}\right)$
$\displaystyle=(-1)^{n}\sum\limits_{m=0}^{m=\infty}(-1)^{m}\binom{n}{m}\left((-1)^{\text{min}\\{v_{j}\\}}\binom{n-m-1}{\text{min}\\{v_{j}\\}}+(-1)^{i}\binom{n-m-1}{i-1}\right)\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$
$\displaystyle=\sum\limits_{m=0}^{m=\infty}(-1)^{m}\binom{n}{m}\left((-1)^{n-\text{min}\\{v_{j}\\}}\binom{n-m-1}{\text{min}\\{v_{j}\\}}+(-1)^{n-i}\binom{n-m-1}{i-1}\right)\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$
$\displaystyle=\sum\limits_{m=0}^{m=\infty}(-1)^{m}\binom{n}{m}\left((-1)^{\text{max}\\{n-v_{j}\\}}\binom{n-m-1}{n-\text{max}\\{n-v_{j}\\}}+(-1)^{n-i}\binom{n-m-1}{i-1}\right)\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$
$\displaystyle=\sum\limits_{m=\text{max}\\{n-v_{j}\\}}^{m=n-i}(-1)^{m}\binom{n}{m}\left((-1)^{\text{max}\\{n-v_{j}\\}}\binom{n-m-1}{n-\text{max}\\{n-v_{j}\\}}+(-1)^{n-i}\binom{n-m-1}{i-1}\right)\prod\limits_{j=1}^{j=m}\binom{m}{n-v_{j}}$

∎

Plots of the PMF/CDFs for some values of the parameters are shown in Figure 2.
The behavior of the test in an example case is illustrated in Figure 3.

###  2.4 CDF generating function and incomplete beta function

The generating function for the values of CDF$(i)$, that is with $i$ and $n$
fixed and the $v=(v_{1},\ldots,v_{k})$ variable, is nearly expressible as the
regularized incomplete beta function $I_{x}(a,b)$ with certain arguments,
establishing a strong analogy to the binomial distribution. The number of
configurations with up to $i$ incidence statistic is given by the generating
function:

$\displaystyle\sum\limits_{v}$
$\displaystyle\sum\limits_{u=0}^{u=i}\binom{n}{u}a(n-u,(v_{1}-u,\ldots,v_{k}-u))t^{v}=\sum\limits_{u=0}^{u=i}\binom{n}{u}(f(t)-t_{1}\cdots
t_{k})^{n-u}(t_{1}\cdots t_{k})^{u}$
$\displaystyle=f(t)^{n}\sum\limits_{u=0}^{u=i}\binom{n}{u}\left(1-\frac{t_{1}\cdots
t_{k}}{f(t)}\right)^{n-u}\left(\frac{t_{1}\cdots t_{k}}{f(t)}\right)^{u}$
$\displaystyle=f(t)^{n}\thinspace I_{1-\frac{t_{1}\cdots
t_{k}}{f(t)}}(n-i,i+1)$

The last equation above is a ”formal” application of the expression for the
CDF of a binomial distribution with $n$ trials, that is,

$\displaystyle\sum\limits_{u=0}^{u=i}\binom{n}{u}p^{u}(1-p)^{n-u}=I_{1-p}(n-i,i+1)$

except that instead of the usual real parameter $p\in[0,1]$ of such a
distribution, $p$ must be permitted to be equal to the power series
$\frac{t_{1}\cdots t_{k}}{f(t)}$ which tabulates information across all of the
different values of the parameters $v=(v_{1},\ldots,v_{k})$.

The total number of configurations is given by the generating function
$(f(t))^{n}$, so the generating function for CDF$(i)$ is the ratio:

$\displaystyle f(t)^{n}I_{1-\frac{t_{1}\cdots
t_{k}}{f(t)}}(n-i,i+1)\thinspace//\thinspace f(t)^{n}$

Here the double division symbol $//$ means the coefficient-wise ratio of the
multi-dimensional series represented by the respective generating functions.
Thus, despite the analogy with the binomial distribution, the generating
function for CDF$(i)$ is not literally equal to $I_{1-\frac{t_{1}\cdots
t_{k}}{f(t)}}(n-i,i+1)$.

##  3 Software implementation

###  3.1 Python package

A Python package coincidencetest is released on PyPI. It contains a self-
contained module, with no dependencies beyond the standard library, that
calculates the $p$-value for the test.

###  3.2 Command-line tool

A command-line tool is distributed with coincidencetest that bundles together
a basic, lightweight signature discovery algorithm as well as test evaluation
on an input binary matrix file. This may be run in a non-interactive context
on a remote server or as part of a pipeline.

###  3.3 Web application

A simple GUI performs signature discovery and evaluation in real-time after
user upload of a binary matrix file. A screenshot is shown in Figure 4.

###  3.4 Testing

The Python package contains a test suite which verifies the $p$-value formulas
(i.e. the PMF and CDF) against brute-force enumerations for several small
values of the parameters, furnishing rigorous computational evidence for the
main theorems in addition to the proofs.

##  4 Related work

The test turns out to specialize to the Fisher exact test[2] in the case of 2
features, $|F|=2$. The incidence statistic and the frequencies of each feature
provide the same information as a $2\times 2$ integer contingency table, and
the formula for the probability value agrees with ours in this case.

The Fisher exact test has been generalized to larger, $r\times c$ contingency
tables[10]. Whether such tables are regarded as pertaining to 2 categorical
variables with $r$ and $c$ categories respectively, or as pertaining to pairs
of binary variables, one from a list of $r$ variables and one from a list of
$c$ variables, contingency table methods are second-order in that they only
involve interactions between pairs of variables. Much work on exact inference
generally has focused on contingency tables, with multi-dimensional
generalizations appearing in the literature up to order 3 (e.g. $I\times
J\times K$ tables[1]).

By contrast our test is inherently higher-order, depending, albeit in a simple
way, on the mutual interaction of all $k$ variables. As for other higher-order
methods, an investigation of the joint distribution of Bernoulli variables
under certain constraints has been published[8], and this may yield a test
with comparable domain of applicability as our test. However, as indicated in
section 2.3, the Bernoulli context involves a different null assumption.

In Good[5] a very similar generating function to our $g(s,t)$ is identified as
a tabulation of the number of _contingency tables_ (not binary matrix
configurations) with fixed column and row sums. The function is
$g(-s,t)^{-1}=g(s,-t)^{-1}$ (c.f. page 1166 item 5.6 and page 1182,
“$f(\mathbf{z})$”). This connection may help to explain the appearance of the
beta function in the generating function for the CDF of the incidence
statistic.

##  Appendix A Formal Concept Analysis bicluster identification

Formal Concept Analysis (FCA)[3] studies a binary data matrix, called a
_formal context_ , in terms of a lattice of certain patterns found in the
matrix. The patterns are known as _(formal) concepts_. Such a concept consists
of a bicluster $(F,S)$, defined as a set of features $F$ and a set of samples
$S$ for which the submatrix along $(F,S)$ consists of all 1s, which is maximal
in two senses: (1) $S$ cannot be enlarged without reducing $F$, and (2) $F$
cannot be enlarged without reducing $S$.

The containment relations of the sets $F$ (respectively $S$) confer a partial
ordering or lattice structure on the set of all concepts, which turns out to
be complete. The maximality condition amounts to a closure condition on the
sets $F$ (respectively $S$), and the whole apparatus can be formulated as a
Galois correspondence between two closure systems on the full feature set and
full sample set.

A straightforward recursive algorithm can be used to enumerate _all_ concepts
in a given context[4]. This algorithm applies to any finite closure system,
and it works by computing the closure of the union of any pair of previously-
found closed sets.

In practice, however, data sets of intermediate size or larger furnish too
many concepts for a complete enumeration to provide a useful direction of
attention to important subsamples or signatures. The present work is partly
motivated by this problem, as it can be used to filter signatures by
significance.

##  Appendix B Figures

Figure 1: (Above) The dichotomized expression of several PSG genes on 510 lung
tumor samples from the TCGA-LUAD project. (Below) The same expression matrix,
with the 19 samples that are positive for all features grouped together on the
left. The number of positives for each feature are respectively 101, 105, 106,
73, 69, 104. The exact test for coincidence yields $p=5.1\cdot 10^{-56}$,
suggesting that the PSG+ phenotype is highly statistically significant. The
loci of the PSG genes are very near to each other, so this is not too
surprising; it is likely that gene amplification events near this locus were
the cause of the observation. Figure 2: (Blue) The probability mass functions
for the incidence statistic at several values of the set sizes $v$ and the
ambient set size $n$. (Red) The cumulative distribution functions. Figure 3:
Illustration of an application of the exact test for coincidence. Figure 4: A
screenshot of the in-browser GUI. Figure 5: A portion of the lattice of
feature subsets in cell data extracted from 29-channel multiplexed mass
cytometry of breast tumor Tissue Micro Arrays (TMA)[6]. The original data are
available on Zenodo[7]. For signature discovery, a random subsample was taken
from the cell table. Green indicates lower $p$-value, and pink indicates
higher $p$-value. The node size and vertical placement are proportional to the
frequency of the sample set displaying the given signature. Only signatures
with frequency between 60 and 400 (out of 2280) are shown.

## References

  * [1] Alan Agresti  “A survey of exact inference for contingency tables”  In _Statistical science_ 7.1 Institute of Mathematical Statistics, 1992, pp. 131–153
  * [2] R.. Fisher  “On the Interpretation of chi-2 from Contingency Tables, and the Calculation of P”  In _Journal of the Royal Statistical Society_ 85.1 [Wiley, Royal Statistical Society], 1922, pp. 87–94  URL: http://www.jstor.org/stable/2340521
  * [3] Bernhard Ganter, Gerd Stumme and Rudolf Wille  “Formal concept analysis: foundations and applications”  springer, 2005
  * [4] Bernhard Ganter, Sergei Obiedkov, Sebastian Rudolph and Gerd Stumme  “Conceptual exploration”  Springer, 2016
  * [5] Irving J Good  “On the application of symmetric Dirichlet distributions and their mixtures to contingency tables”  In _The Annals of Statistics_ 4.6 Institute of Mathematical Statistics, 1976, pp. 1159–1189
  * [6] Hartland W Jackson et al.  “The single-cell pathology landscape of breast cancer”  In _Nature_ 578.7796 Nature Publishing Group, 2020, pp. 615–620
  * [7] Hartland W. Jackson et al.  “The Single-Cell Pathology Landscape of Breast Cancer”  Zenodo, 2019  DOI: 10.5281/zenodo.3518284
  * [8] Nikolai Kolev, Ekaterina T Kolkovska and José Alfredo López-Mimbela  “Joint probability generating function for a vector of arbitrary indicator variables”  In _Journal of computational and applied mathematics_ 186.1 Elsevier, 2006, pp. 89–98
  * [9] James C Mathews et al.  “Functional network analysis reveals an immune tolerance mechanism in cancer”  In _Proceedings of the National Academy of Sciences_ 117.28 National Acad Sciences, 2020, pp. 16339–16345
  * [10] Daniel Zelterman, Ivan Siu-Fung Chan and Paul W Mielke Jr  “Exact tests of significance in higher dimensional tables”  In _The American Statistician_ 49.4 Taylor & Francis, 1995, pp. 357–361