# Generation of anonymous chest radiographs using latent diffusion models for
training thoracic abnormality classification systems

###### Abstract

The availability of large-scale chest X-ray datasets is a requirement for
developing well-performing deep learning-based algorithms in thoracic
abnormality detection and classification. However, biometric identifiers in
chest radiographs hinder the public sharing of such data for research purposes
due to the risk of patient re-identification. To counteract this issue,
synthetic data generation offers a solution for anonymizing medical images.
This work employs a latent diffusion model to synthesize an anonymous chest
X-ray dataset of high-quality class-conditional images. We propose a privacy-
enhancing sampling strategy to ensure the non-transference of biometric
information during the image generation process. The quality of the generated
images and the feasibility of serving as exclusive training data are evaluated
on a thoracic abnormality classification task. Compared to a real classifier,
we achieve competitive results with a performance gap of only 3.5 % in the
area under the receiver operating characteristic curve.

Index Terms—  Chest Radiography, Synthetic Image Generation, Abnormality
Classification, Patient Privacy

##  1 Introduction

Over the past years, the availability of novel deep learning (DL) [1]
techniques has led to significant advances in the development of diagnostic
algorithms for the detection and classification of thoracic abnormalities [2].
As the training of DL models is data-driven, it typically requires enormous
amounts of training samples to obtain well-performing networks. However,
medical images, e. g., chest radiographs, include personal information and are
subject to certain privacy protection regulations, such as the EU General Data
Protection Regulation (GDPR), which hinders the widespread sharing and, thus,
utilization of acquired data. Therefore, to make large-scale medical image
collections publicly available for research purposes, it is necessary to
strictly comply with applicable objectives and requirements for the robust and
reliable anonymization of sensitive patient data.

Conventional techniques for anonymizing chest radiographs rely on removing
personally identifiable metainformation, e. g., patient names, or replacing
them with pseudonyms. Furthermore, in the image domain itself, critical areas
are typically obscured with black boxes. In one of our recent studies [3],
however, we empirically demonstrate that such approaches do not offer
sufficient means of protecting patient privacy. This is primarily because
chest radiographs inherently contain biometric information (similar to a
fingerprint) that can be exploited to successfully re-identify specific
patients by DL-based linkage attacks. Hence, there is an urgent need for more
sophisticated privacy-enhancing techniques (PETs) in medical imaging.

In this context, the generation of synthetic medical images with privacy
guarantees has attracted attention as a promising solution for anonymization
and overcoming data-sharing limitations. Generative models learn the
probability distribution of a given real dataset and can be used to synthesize
realistic representations by sampling from the learned data distribution.
Provided that synthetic medical images are of sufficient quality, such data
could be used to train DL models at scale. Han et al. [4] utilized
progressively growing generative adversarial networks (PGGANs) [5] to
synthesize chest radiographs and demonstrated the feasibility of GAN-based
data augmentations. Moreover, they confirm the utility of exclusively
synthetic training data in a clinical setting. Recently, diffusion
probabilistic models [6] have become state-of-the-art for tasks such as high-
resolution image synthesis. Chambon et al. [7] leverage a latent diffusion
model (LDM) [8] for the generation of class-conditional chest radiographs.
Nevertheless, the authors do not investigate whether the generated images can
convey relevant class information during the training process of an
abnormality classifier. Furthermore, their dataset is limited to only 2,000
samples.

In this work, we aim to investigate the feasibility of using exclusively
synthetic training datasets for learning and recognizing thoracic
abnormalities in chest radiographs. Therefore, we utilize an LDM [8] to
generate high-quality class-conditional images from a given real data
distribution. To synthesize a fully anonymous dataset from the trained
generative model without transferring biometric patient identifiers, we
propose a rigorous privacy-enhancing sampling strategy that excludes synthetic
images, provided a specific patient identity has been learned from the
original training data. Finally, to quantify the quality of the generated
anonymous dataset and to evaluate its applicability for training a thoracic
abnormality identification system, we conduct multiple classification
experiments and compare the performances of classifiers trained on real versus
synthetic data, respectively. Throughout our work, we leverage a PGGAN [5] as
a synthetic image generation baseline.

##  2 Methods

###  2.1 Dataset

We use the large-scale ChestX-ray14 dataset [9] consisting of 112,120 chest
radiographs from 30,805 individual patients. Accompanying metadata provides
information about the corresponding 14 abnormality labels, including
Atelectasis, Cardiomegaly, Consolidation, Edema, Effusion, Emphysema,
Fibrosis, Hernia, Infiltration, Mass, Nodule, Pleural Thickening, Pneumonia,
and Pneumothorax. Healthy subjects are labeled with an additional class,
indicating that none of the above-mentioned abnormalities are present. To
simplify the class-conditional image generation task in this work, we exclude
all images with multiple abnormality labels. Furthermore, we limit the maximum
number of follow-up scans per patient to 5 to prevent the generative models
from learning patient-specific patterns from over-represented subjects.
Lastly, as the anatomy of young patients may not yet be fully developed, we
only consider images from patients older than 21 years. After this data
reduction procedure, the remaining 56,352 images are split into a training,
validation, and test set by a ratio of 70:10:20. To this end, a patient-wise
splitting strategy is applied to avoid potential patient overlap between the
subsets while maintaining the overall class distribution.

###  2.2 Latent diffusion model

For the class-conditional image generation task, we apply the LDM proposed in
[8] leveraging pre-trained autoencoders. In this approach, input images are
first embedded into a latent space of size 64$\times$64$\times$3 using the
encoder of a vector quantized-variational autoencoder (VQ-VAE) with 32, 64,
and 128 channels in each stage. A diffusion model operates within that lower
space dominated by lower frequencies and performs 1000 denoising steps with a
U-Net (32, 128, and 256 channels). In the final step, the decoder of the
autoencoder increases the spatial resolution to 256$\times$256 pixels while
introducing higher frequencies. The class-conditional information is
incorporated using a trainable lookup table. This is realized by combining the
class embeddings with the diffusion process using cross-attention in the
bottleneck of the U-Net. We closely followed the implementation by Rombach et
al. [8].

Fig. 1: Illustration of the anonymous image generation process using the
proposed privacy-enhancing sampling strategy. We apply generative modeling (1)
techniques to synthesize chest radiographs. We then use patient retrieval (2)
and verification (3) networks to ensure the non-transference of biometric
identifiers. The resulting anonymous dataset maintains the size and class
distribution of the original training set and is used for training an
abnormality classifier (4).

###  2.3 Progressively growing GAN

As a synthetic image generation baseline, we employ the PGGAN architecture
proposed in [5]. PGGANs are based on the idea of starting the training process
with low-resolution images before gradually increasing the image resolution by
adding layers to the network architecture. This allows the generator to focus
on large-scale image structures at the beginning of the training procedure and
to learn finer high-resolution details in later stages. We use a generator and
discriminator that have the same general structure and grow equally over time.
Both networks incorporate replicated 3-layer blocks consisting of an
upsampling or downsampling layer, a 3$\times$3 convolutional layer, and a
leaky rectified linear unit (LReLU) with a slope of 0.2. We adapt the
originally proposed PGGAN architecture to produce chest radiographs with a
final spatial resolution of 256$\times$256 pixels. The input to the generator
is represented by a 256-dimensional latent noise vector. Its last 15 entries
consist of the conditional vector, which controls the image class to be
generated. For the discriminator, we use a projection-based approach [10] to
incorporate conditional information into the network.

###  2.4 Privacy-enhancing sampling strategy

While synthetic data is often referred to as anonymous per se, generative
models may memorize training examples and reproduce specific patient
identities during inference [11]. To ensure the non-transference of biometric
patterns when generating synthetic chest radiographs, we propose a rigorous
privacy-enhancing sampling strategy (see Fig. 1) that leverages state-of-the-
art patient matching approaches. The strategy is based on the idea of
excluding synthetic images in case patient-specific biometrics have been
reproduced in a generated scan. To realize such an image selection process, we
first apply a patient retrieval network for each generated chest radiograph to
find its top-1 image in terms of patient similarity from the original training
set. Then, we perform a 1-to-1 matching step for the recommended image pairs
(real vs. synthetic) by applying a patient verification network that yields
the probability of whether or not the synthetic image and the recommended real
image belong to the same patient. Based on this probability value, we exclude
all synthetic images that exceed a threshold of $t=0.5$. We employ pre-trained
patient retrieval and verification models proposed in previous work [3]. The
used networks achieve a top-1 precision of 99.6 % and an area under the
receiver operating characteristic curve (AUC) of 99.4 %, respectively.

###  2.5 Thoracic abnormality classification

For the downstream abnormality classification task, we employ CheXNet [12], a
densely connected convolutional network (DenseNet) [13] consisting of 121
layers. Its input layer receives chest radiographs with a resolution of
224$\times$224 pixels. The final classification layer yields a 14-dimensional
output vector indicating the probability of the presence or absence of each
abnormality class appearing in the ChestX-ray14 dataset.

##  3 Experiments

All experiments were conducted using PyTorch and Python. The experimental
setup included the following three steps.

###  3.1 Training of the generative models

The first step of our experimental pipeline consisted of training the
generative models. For the LDM, we first trained the autoencoder with a
combination of the perceptual loss [14] and a patch-based adversarial
objective [15] for 100 epochs using a learning rate of $4.5\cdot 10^{-6}$.
Then, the autoencoder was applied to embed the input images into latent
representations, which were used as inputs for the underlying diffusion model.
The diffusion model was trained with a learning rate of $10^{-6}$ using the
$L_{1}$ loss until no further improvement in the validation loss (250 epochs).

The PGGAN architecture was trained using the Wasserstein loss with a gradient
penalty coefficient of 10. For optimization, we used the Adam [16] optimizer
with a learning rate of 0.001. The training was initiated at an image
resolution of 4$\times$4\. We gradually added layers to both the generator and
the discriminator to increase the resolution by factors of 2 until an image
size of 256$\times$256 was reached. Each stage of the PGGAN was trained until
the generator and discriminator loss values stabilized (100 epochs).

###  3.2 Privacy-enhancing image generation

After training the generative models, we synthesized two anonymous chest X-ray
datasets by leveraging the trained LDM and PGGAN, respectively. For this
purpose, we applied the proposed privacy-enhancing sampling strategy to ensure
the non-transference of biometric identifiers during the image generation
process. Note that the size and class distribution of the real training and
validation set was maintained when creating the synthetic datasets. No
oversampling of under-represented classes was performed. Both the LDM-based
and the GAN-based datasets were used for further experiments as described in
the following step.

Table 1: Comparison of the abnormality identification performance of CheXNet
when using either real or exclusively synthetic training sets. For all
evaluations, the same real testing data is employed. We report the means and
standard deviations of the resulting AUC values after 10 independent training
and testing runs. Best synthetic results are shown in bold.

Training set |  Real |  $\text{Syn}_{\text{PGGAN}}$ |  $\text{Syn}_{\text{LDM}}$
---|---|---|---
Atelectasis |  $81.3\pm 0.8$ |  $70.1\pm 1.2$ |  $\bm{76.2\pm 0.4}$
Cardiomegaly |  $92.9\pm 0.6$ |  $86.4\pm 1.4$ |  $\bm{88.6\pm 0.8}$
Consolidation |  $74.8\pm 1.0$ |  $68.0\pm 2.8$ |  $\bm{75.5\pm 1.1}$
Edema |  $92.8\pm 0.8$ |  $84.4\pm 2.8$ |  $\bm{87.5\pm 1.7}$
Effusion |  $90.7\pm 0.4$ |  $83.2\pm 0.9$ |  $\bm{85.9\pm 0.9}$
Emphysema |  $88.1\pm 0.8$ |  $76.5\pm 1.6$ |  $\bm{83.9\pm 1.0}$
Fibrosis |  $80.8\pm 1.0$ |  $69.4\pm 2.9$ |  $\bm{77.3\pm 0.6}$
Hernia |  $93.5\pm 1.5$ |  $80.6\pm 3.5$ |  $\bm{93.7\pm 1.4}$
Infiltration |  $68.7\pm 0.2$ |  $59.1\pm 0.6$ |  $\bm{63.4\pm 0.5}$
Mass |  $81.0\pm 1.3$ |  $67.7\pm 0.9$ |  $\bm{76.9\pm 0.8}$
Nodule |  $71.4\pm 1.0$ |  $60.8\pm 1.2$ |  $\bm{67.8\pm 0.7}$
Pleural Thickening |  $75.5\pm 1.2$ |  $68.4\pm 1.3$ |  $\bm{73.4\pm 0.8}$
Pneumonia |  $70.8\pm 3.3$ |  $61.5\pm 5.1$ |  $\bm{65.8\pm 2.0}$
Pneumothorax |  $79.9\pm 0.7$ |  $70.1\pm 1.5$ |  $\bm{77.0\pm 0.5}$
Mean |  $81.6\pm 0.4$ |  $71.9\pm 0.8$ |  $\bm{78.1\pm 0.3}$

(a) Infiltration

(b) Nodule

(c) Mass

(d) Cardiomegaly

(a) PGGAN

(e) Infiltration

(f) Nodule

(g) Mass

(h) Cardiomegaly

(b) LDM

Fig. 2: Randomly selected images generated by (a) the trained PGGAN and (b)
the LDM, respectively. The positions and areas of induced abnormalities are
highlighted with red arrows or circles.

###  3.3 Thoracic abnormality classification

To assess the extent to which synthetic data can convey relevant class
information during the training process of a thoracic abnormality classifier,
we conducted experiments with both a real training set and with exclusive
synthetic training sets generated by either the trained PGGAN or the LDM
($\text{Syn}_{\text{PGGAN}}$ and $\text{Syn}_{\text{LDM}}$). We trained the
classifiers using the class-wise binary cross-entropy (BCE) loss and
stochastic gradient descent (SGD) with a momentum of 0.9 and weight decay of
0.0001. The initial learning rate was set to 0.01 and was divided by a factor
of 10 if the loss did not improve after an epoch. Model training was conducted
until early stopping with a patience of 3. Each individually trained
classifier was evaluated using the same real testing set. For each scenario,
we performed 10 independent training and testing runs. The network
performances are compared by reporting the means and standard deviations of
the resulting AUC values.

##  4 Results

###  4.1 Quantitative results

The results of the conducted thoracic abnormality classification experiments
are summarized in Table 1. As can be seen, the overall best performance is
achieved when training the classifier with real images from the ChestX-ray14
dataset. In this case, we obtain a mean AUC of 81.6 %. Compared to this
reference, the results significantly decrease for the classifier trained with
the synthetic dataset generated by the PGGAN. The mean AUC of 71.9 % indicates
that the GAN-based images fail to convey relevant class information during the
training of the classification system. In contrast, when using training images
generated by the LDM, we observe competitive classification results with a
mean AUC of 78.1 % on our test set. Thus, the performance is only slightly
below the reference, indicating the ability of the LDM-based chest radiographs
to adequately preserve class-specific patterns. Interestingly, the performance
for pathology classes Consolidation and Hernia is even higher than the
reference.

###  4.2 Qualitative results

Randomly selected chest radiographs generated by the trained PGGAN and the LDM
are shown in Fig. 2. In direct comparison, we perceive better image quality
for the LDM-based chest radiographs. The trained LDM can produce high-
resolution images with realistic patient anatomy. Moreover, as seen in the
figure, each of the presented abnormalities becomes visible in the generated
images. In contrast, the GAN-based chest radiographs show a slightly worse
image quality. This becomes apparent, for instance, by the partially
unrealistic shape of the ribs and lungs, or through isolated artifacts that
appear in the images. Moreover, we observe that the induced class conditions
are a bit less pronounced in the images generated by the PGGAN.

##  5 Discussion and Conclusion

In this paper, we proposed a privacy-enhancing image sampling strategy to
synthesize anonymous chest X-ray datasets. We leveraged an LDM and a PGGAN to
generate high-quality class-conditional images. The obtained quantitative and
qualitative results indicate the competitiveness of diffusion-based chest
radiographs for training a thoracic abnormality classification system. With
our work, we do not want to suggest that synthetic data can replace real data.
In a real-life scenario, this would only be applicable if the quality of
synthetic data is on par with real data. However, if this can be guaranteed,
we hypothesize that our proposed approach allows for significant mitigation of
data-sharing limitations.

In future work, next to further improving the overall image quality, we aim to
exploit the generative models for compensating under-represented classes or
creating infinitely large synthetic datasets. Lastly, we plan to incorporate a
privacy-enhancing mechanism directly into the LDM to ensure the non-
transference of biometric information during the training process.

##  6 Compliance with ethical standards

This research study was conducted retrospectively using human subject data
made available in open access by the National Institutes of Health (NIH)
Clinical Center [9]. Ethical approval was not required as confirmed by the
license attached with the open-access data.

##  7 Acknowledgments

The research leading to these results has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (ERC Grant no. 810316). The authors gratefully acknowledge
the scientific support and HPC resources provided by the Erlangen National
High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-
Universität Erlangen-Nürnberg (FAU). The hardware is funded by the German
Research Foundation (DFG). The authors declare that they have no conflicts of
interest.

## References

  * [1] Andreas Maier, Christopher Syben, Tobias Lasser, and Christian Riess,  “A gentle introduction to deep learning in medical image processing,”  Zeitschrift für Medizinische Physik, vol. 29, no. 2, pp. 86–101, 2019.
  * [2] Sebastian Gündel, Sasa Grbic, Bogdan Georgescu, Siqi Liu, Andreas Maier, and Dorin Comaniciu,  “Learning to Recognize Abnormalities in Chest X-Rays with Location-Aware Dense Networks,”  in Iberoamerican Congress on Pattern Recognition. Springer, 2018, pp. 757–765.
  * [3] Kai Packhäuser, Sebastian Gündel, Nicolas Münster, Christopher Syben, Vincent Christlein, and Andreas Maier,  “Deep learning-based patient re-identification is able to exploit the biometric nature of medical chest X-ray data,”  Scientific Reports, vol. 12, no. 1, pp. 1–13, 2022.
  * [4] Tianyu Han et al.,  “Breaking medical data sharing boundaries by using synthesized radiographs,”  Science Advances, vol. 6, no. 49, pp. eabb7973, 2020.
  * [5] Tero Karras, Timo Aila, Samuli Laine, and Jaakko Lehtinen,  “Progressive Growing of GANs for Improved Quality, Stability, and Variation,”  in International Conference on Learning Representations, 2018.
  * [6] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli,  “Deep Unsupervised Learning using Nonequilibrium Thermodynamics,”  in International Conference on Machine Learning. PMLR, 2015, pp. 2256–2265.
  * [7] Pierre Chambon, Christian Bluethgen, Curtis P Langlotz, and Akshay Chaudhari,  “Adapting Pretrained Vision-Language Foundational Models to Medical Imaging Domains,”  arXiv preprint arXiv:2210.04133, 2022.
  * [8] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer,  “High-Resolution Image Synthesis with Latent Diffusion Models,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2022, pp. 10684–10695.
  * [9] Xiaosong Wang, Yifan Peng, Le Lu, Zhiyong Lu, Mohammadhadi Bagheri, and Ronald M Summers,  “ChestX-ray8: Hospital-scale Chest X-ray Database and Benchmarks on Weakly-Supervised Classification and Localization of Common Thorax Diseases,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 2097–2106.
  * [10] Takeru Miyato and Masanori Koyama,  “cGANs with Projection Discriminator,”  in International Conference on Learning Representations, 2018.
  * [11] Aditya Kuppa, Lamine Aouad, and Nhien-An Le-Khac,  “Towards Improving Privacy of Synthetic DataSets,”  in Annual Privacy Forum. Springer, 2021, pp. 106–119.
  * [12] Pranav Rajpurkar et al.,  “CheXNet: Radiologist-Level Pneumonia Detection on Chest X-Rays with Deep Learning,”  arXiv preprint arXiv:1711.05225, 2017.
  * [13] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger,  “Densely Connected Convolutional Networks,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 4700–4708.
  * [14] Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang,  “The Unreasonable Effectiveness of Deep Features as a Perceptual Metric,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018, pp. 586–595.
  * [15] Patrick Esser, Robin Rombach, and Björn Ommer,  “Taming Transformers for High-Resolution Image Synthesis,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 12873–12883.
  * [16] Diederik P Kingma and Jimmy Ba,  “Adam: A Method for Stochastic Optimization,”  arXiv preprint arXiv:1412.6980, 2014.


# An efficient counting method for the colored triad census

Jeffrey Lienert  Laura Koehly  Felix Reed-Tsochas  Christopher Steven Marcum
Social Network Methods Section, National Human Genome Research Institute,
National Institutes of Health  CABDyN Complexity Centre, Saïd Business School,
University of Oxford  Corresponding Author

###### Abstract

The triad census is an important approach to understand local structure in
network science, providing comprehensive assessments of the observed
relational configurations between triples of actors in a network. However,
researchers are often interested in combinations of relational and categorical
nodal attributes. In this case, it is desirable to account for the label, or
color, of the nodes in the triad census. In this paper, we describe an
efficient algorithm for constructing the colored triad census, based, in part,
on existing methods for the classic triad census. We evaluate the performance
of the algorithm using empirical and simulated data for both undirected and
directed graphs. The results of the simulation demonstrate that the proposed
algorithm reduces computational time many-fold over the naïve approach. We
also apply the colored triad census to the Zachary karate club network
dataset. We simultaneously show the efficiency of the algorithm, and a way to
conduct a statistical test on the census by forming a null distribution from
$1,000$ realizations of a mixing-matrix conditioned graph and comparing the
observed colored triad counts to the expected. From this, we demonstrate the
method’s utility in our discussion of results about homophily, heterophily,
and bridging, simultaneously gained via the colored triad census. In sum, the
proposed algorithm for the colored triad census brings novel utility to social
network analysis in an efficient package.

###### keywords:

triad census , labeled graphs , simulation

††journal: XXXX

##  1 Introduction

The triad census is an important approach towards understanding local network
structure. Holland1976 first presented the 16 isomorphism classes of
structurally unique triads possible in a directed network without loops. To
conduct a triad census, one simply counts each occurrence of these structures,
without respect to the labeling of the nodes (here we use node label, color,
characteristic, and attribute interchangeably). This is useful insofar as
specific triads, or combinations thereof, may relate to underlying social
processes giving rise to an observed network. For example, bridges (triads
with one null dyad and two non-null dyads) may be important in navigating
social networks [Granovetter1973], and certain triads may be more or less
favorable based on structural balance theory (e.g. the 300 is balanced but the
201 is not, see Figure 1) [Cartwright1956]. Moreover, a variant of the triad
census, motif analysis, investigates the statistics of various triad
configurations (motifs), and has found wide application in biology
[milo2002network].

Also important to network structure are nodal characteristics and how they
relate to tie formation or dissolution. This has been the subject of research
on homophily (individuals having similar attributes with those to whom they
are connected) [McPherson2014]. However, homophily is an observed phenomenon,
not a process. The processes giving rise to homophily are varied, often
confound the relationship between networks and outcomes, and are difficult to
tease apart [Shalizi2012]. Methodological advances, such as stochastic actor-
oriented models can disentangle these effects to some extent
[snijders1996stochastic]. Other analyses have attempted to disentangle the
processes leading to homophily from structural processes, such as triadic
closure [goodreau&2009bfff]. Additionally, the coloring of nodes in a network
has been an important question for many graph theorists and indeed represents
a major topic in this field [Jensen2011].

Although nodal characteristics and the triad census are important, they have
rarely been examined fully in conjunction. Yet, there are a few cases where
specific colored triads have been studied. For example, gould1989structures
study brokerage based on triad structure and group membership simultaneously.
This same approach has been used to study brokerage in dynamic networks
[SPIRO2013130]. As well, a study by Marcum2015 examined specific colored
triads based on generational membership within families; in this work the
authors showed that inter-generational ties were observed in different
quantities than expected based on the underlying null model. None of the past
research evaluated the full census of colored triads, rather, researchers have
focused instead on specific colored triads that were _a priori_ expected to be
relevant to the processes at hand. As a result, these foundational works were
not exhaustive with respect to all alternatives. In other words, previous
research examining a subset of colored triads likely had an amount of false
negatives due to not examining every colored triad; this could be addressed by
censusing the colored triads.

The examination of node characteristics together with local structure is
important as it provides opportunity to simultaneously study the occurrence of
triadic structure, nodal attributes, and the interactions between them. For
instance, certain colored triads may be forbidden, such as three-cycles
between strict heterosexuals in mixed-orientation sexual contact networks
[marcum2016fourcycle]. Impermissible triads would be categorized the same as
those that were not observed due to chance in a triad census, potentially
missing important social processes or constraints at play in this type of
network. Only by incorporating node coloring into the triad census can this
pattern be fully elucidated.

Based on this methodological gap in the literature, we develop a method to
census the colored triads for any one-mode binary network with arbitrary
number of colors. Due to the large numbers of isomorphism classes of size 3 as
the number of colors increases, this method requires computational efficiency
in addition to mathematical accuracy. As well, one is often interested in
forming a null distribution with which to compare observed colored triad
counts. If the null distribution cannot be analytically solved, one would
likely census the colored triads of many simulated networks, further
increasing the need for the algorithm to be computationally efficient.

Current efficient methods for the triad census exploit the sparseness of
networks [Batagelj2001], and scale sub-quadratically (as the number of edges
increases the time to run the algorithm is faster than the number of edges
squared). However, methods that exploit network sparseness by inferring the
number of null triads do not work in the colored case because they do not
explicitly interrogate every triad, and there are variations within the null
triads due to the coloring. Therefore, we extend the methodology of Moody1998,
which is based on matrix algebra and interrogates every triad; his method
scales sub-quadratically with the number of nodes.

This paper (1) presents the colored triad census and its computational
complexity, (2) shows that this approach can be used on large networks (tested
for up to $10,000$ nodes) with up to 10 colors in relatively efficient time,
and (3) uses the method many times to create null distributions of colored
triad censuses to form the basis of conditional uniform graph tests. We
illustrate the benefits of an analysis incorporating the colored triad census
using a well-known dataset, Zachary’s Karate Club [zachary1977information].

##  2 Algorithm

Since the original appearance of the triad census in 1976, a number of papers
have explored how to compute the triad census of a network in an efficient
manner. Although, for sparse networks, sub-quadratic methods (in terms of
number of nodes) exist for calculating the triad census (e.g. Batagelj2001),
we use the quadratic algorithm presented by Moody1998 here. This is because
the more efficient methods avoid interrogating null triads directly by taking
advantage of the sparseness of graphs, the subsequent large number of null
(003) triads, and the known number of total triads. Instead, they interrogate
all triads with at least one edge, and then subtract that count from the total
number of triads in the network to arrive at the number of null triads. This
is insufficient in the colored triad census as there are differently-colored
null triads, and the count of each cannot therefore be algebraically
determined. For example, if there are two colors, four different null colored
triads exist (0-3 nodes color A). The exact breakdown of the null triad into
the four colored triads cannot be determined without interrogating each null
triad, thereby losing the efficiency gained when not considering colors.
Moody’s algorithm does not employ this limiting shortcut, and we therefore use
it as a basis for our colored triad census algorithm. Additionally, because
many networks are sparse, we can leverage computational techniques for
increasing the efficiency of sparse matrix operations [duff2017direct],
further reducing the computational complexity of our method.

Moody1998 showed that the count of each of the 16 triad isomorphism classes
could be derived by using matrix algebra on the adjacency matrix of the graph
and its derivatives. To review, let $\mathbf{A}$ be the adjacency matrix of a
network, and $A_{ij}=1$ when a tie exists from node $i$ to node $j$. Let $E$
be the symmetrized matrix $A$, formed by making any edge in $A$ reciprocal via
$E_{ij}=\max(A_{ij},A_{ji})$. The complement of $E$, $\bar{E}$, is formed by
subtracting the complete network adjacency matrix from $E$, so that $E_{ij}=1$
if and only if there is neither a tie from $i$ to $j$ nor a tie from $j$ to
$i$. Next, we have $M$, the mutual matrix of $A$, and is made by removing any
asymmetric edges from $A$, or $M_{ij}=\min(A_{ij},A_{ji})$. Finally, $C$ is
the matrix of only asymmetric edges, and is calculated by $C=A-M$. Therefore,
$C_{ij}=1\iff A_{ij}=1\;\&\;A_{ji}=0$. Based on these matrices, Moody
demonstrates how to calculate the number of each of the 16 isomorphism classes
for the case of unlabeled graphs (or, equivalently, for a graph consisting of
nodes of the same single color). Generally, this was done by multiplying
(either through dot-product or element-wise multiplication) the three matrices
corresponding to the relevant edges in the triad of interest. There were two
triads ($111U$ and $111D$) that were not directly amenable to this process and
were calculated via addition and subtraction of other triad types,
respectively.

To extend this work to the case of multiple colors, we introduce the out-
coloring and in-coloring matrices, $K^{r}$ and $K^{r^{\prime}}$, respectively,
where $r$ is the focal color of matrix $K$. Here, the in-coloring matrix is
the transpose of the out-coloring matrix. The out-coloring matrix is
calculated by evaluating the color of the nodes row-wise, such that rows
indexing nodes of the focal color are composed in the following way:

$K^{r}_{i\bullet}=\begin{cases}1&\text{if }R(i)=r\\\ 0&\text{if }R(i)\neq
r\end{cases}$ (1)

Where $R(i)$ is a function returning the color of node $i$. As above, the in-
coloring matrix is the transpose of the out-coloring matrix in Eq. 1.

Our algorithm works by using the in- and out-coloring matrices to evaluate and
“switch on” edges that have nodes of the focal colors at the ends (or tails)
of edges in the adjacency matrix $A$ of the network. We adapt the triad census
nomenclature of Holland1976 by appending the colors after the name of the
triad. The colors are ordered from the top node proceeding clockwise in Figure
1. We have arbitrarily adapted the orientation of the triads from the triad
census figure in Holland1976 for computational reasons. The orientation is
important here because triads with the same orientation may no longer be
isomorphic when color is introduced. Figure 1 makes it possible to count
unambiguously and name only unique colored triads. Therefore, $T_{102-123}$ is
the triad consisting of 1 symmetric dyad and 2 null dyads, where the top node
is of color $1$, the bottom-right node is of color $2$ and the bottom-left
node is of color $3$. This is distinct from the $T_{102-312}$ triad because
the coloring of the nodes is not identical from the previous triad.

Figure 1: The 16 isomorphism classes of triads and their orientation used here
with respect to the color numbering. When colors are added to these triads,
they are labeled starting from the top node and proceeding clockwise.

Following this, the general formula for an arbitrary triad “T” with an
arbitrary coloring triplet is:

$T=\operatorname{Tr}\big{(}(K^{1}\times H(T,1,2)\times
K^{2^{\prime}})(K^{2}\times H(T,2,3)\times K^{3^{\prime}})(K^{3}\times
H(T,3,1)\times K^{1^{\prime}})\big{)}$ (2)

In the above, “$\times$” refers to element-wise multiplication, and “Tr” is
the trace function. For an arbitrary triad, “T” has a color triplet
$r_{1},r_{2},r_{3}$. $H(T,i,j)$ is a function returning the matrix specific to
the type of edge between nodes $i$ and $j$ in triad “$T$”. For example, in a
$102$ triad, the first edge from the top node going clockwise is a symmetric
edge from node one to node two (Figure 1). $H(T_{102},1,2)$ in this case would
be the matrix $E$ for the symmetric matrix, and the sandwiching color matrices
would turn the proper edges on and off if nodes one and two were of the
specified colors. If the edge is an asymmetric one, and the direction of the
edge in Figure 1 is counter-clockwise, then $C^{\prime}$ is used instead of
$C$ to force the edge to go in the proper direction.

At this point, there are redundant triads due to certain colored triads being
isomorphic. For instance, the $T_{003-122}$ is isomorphic with $T_{003-221}$
and $T_{003-212}$, and would be triple-counted. These are removed by checking
for isomorphisms based on matrix row and column permutations of the triad. If
two colored matrices are identical after such row and/or column permutations,
then they are isomorphic, and one is removed. We arbitrarily decide to discard
the triad whose coloring triplet name comes second alphanumerically. It should
be noted that removing in this way is computationally expensive, particularly
as the number of colors and nodes grows large. We therefore shorten this
process by performing it once for 1 to 10 colors and storing the unique
isomorphism classes. This leaves only unique isomorphism classes of colored
triads, which can then be accessed in linear time.

The number of unique isomophism classes for a given number of colors can be
shown for each of the 16 ismorphism classes in the triad census. The 16
classes separate into four types of colored triads, depending on how many
structurally-distinct positions there are in the triad (e.g. the two ends of
the edge in a 102 triad are not structurally-distinct from one another, but
are distinct from the node with no edges). The calculation for the number of
each isomorphism class for arbitrary number of colors ($k$) is shown in Table
1. Each combinatoric term in each row (together with their respective leading
permutation coefficients) counts the number of colored triads when there are
three, two, or one unique color(s), respectively. For example, in a network
with three colors,the ‘300’ and ‘003’ classes have only one accessible
permutation when there are three colors present in the triad (i.e.
$\binom{3}{3}$), six ways when there are two colors (i.e. $2\,\binom{3}{2}$),
and one way when there is one color in the triad (i.e. $\binom{3}{1}$).

Isomorphism classes | Number of colored triads
---|---
$300$ and $003$ |  ${{k}\choose{3}}+2{{k}\choose{2}}+{{k}\choose{1}}$
$030C$ | $2{{k}\choose{3}}+2{{k}\choose{2}}+{{k}\choose{1}}$
$102$ $021D$ $021U$ $201$ $120D$ and $120U$ | $3{{k}\choose{3}}+4{{k}\choose{2}}+{{k}\choose{1}}$
$012$ $021C$ $111D$ $111U$ $030T$ $120C$ and $210$ | $6{{k}\choose{3}}+6{{k}\choose{2}}+{{k}\choose{1}}$
Table 1: Expression for the number of isomorphism classes within a triad
class. $k$ is the number of colors

If these numbers are summed over the 16 isomorphism classes, the total number
of colored isomorphism classes of triads for $k$ colors is returned.
Similarly, the same can be done for undirected triads, solely summing over the
$4$ triads observed in the undirected case. Table 2 reports the total number
of colored triads for undirected and directed networks over a range of $k$.
Clearly, the number of isomorphism classes grows quite quickly as $k$
increases.

Number of colors |  Number of directed colored triads |  Number of undirected colored triads
---|---|---
1 | 16 | 4
2 | 104 | 20
3 | 328 | 56
4 | 752 | 120
5 | 1440 | 220
6 | 2456 | 364
7 | 3864 | 560
8 | 5728 | 816
9 | 8112 | 1140
10 | 11080 | 1540
Table 2: The number of colored triad isomorphism classes for directed and
undirected networks for $k$ ranging from 1 to 10.

The algorithm implemented as an R package is publicly available and is linked
to this paper via github: https://github.com/jlienert/ColoredTriadCensus.

##  3 Algorithmic Performance

If a naïve implementation of matrix multiplication is used, this algorithm
runs with computational complexity $O(N^{3}*3^{k})$. It scales with the number
of nodes cubed ($N^{3}$) because of the matrix multiplication involved in the
algorithm. However, many software packages use algorithms that reduce the
complexity of matrix multiplication to $O(N^{2.38})$ [davie_stothers_2013].
Furthermore, by taking advantage of methods for matrix multiplication using
sparse matrices (as appropriate due to the sparse nature of most social
networks), this complexity is reduced to something closer to $O(N^{2})$
[yuster2005fast]. The exact benefit gained by using sparse matrix
multiplication varies based on how sparse the matrix is. This ranges from the
nearly-optimal $O(N^{2})$ when very few edges exist, to worse than the
optimized algorithm when many edges exist. The scaling with $3^{k}$ comes from
the number of distinct colored triads the algorithm needs to evaluate, and the
number of isomorphism classes scales in such a manner.

To test the efficiency of the algorithm, we apply it to networks ranging in
size from $n=10$ to $n=10,000$ with the number of colors ranging from $k=3$ to
$k=10$, all holding the average degree constant at 6 by creating Erdös-Rényi
graphs with an edge probability of $\frac{6}{N-1}$. This reflects the average
number of ties participants enumerate in social networks surveys
[marsden2003interviewer]. The runtime of the algorithm with these parameters
can be seen in Figure 2. In general, increasing $K$ results in constant
increases in $\log(runtime)$, which is what we expect based on the theoretical
computational complexity. As expected, we also observe a super-quadratic
increase in $\log(runtime)$ as $N$ increases. Although it is super-linear, it
is still below the curve that would exist if we used matrix multiplication not
optimized for sparse matrices (dotted line in Figure 2). This difference shows
the expected time saved by using sparse matrix methods. Finally, we observe
changes in the rank-order and decreases in runtime going from $10$ to $100$
nodes. This is also due to the computational time involved in initializing the
sparse matrices and storing and operating on sparse matrices, and as such is
not unexpected. Additionally, because the average degree was held constant,
the smaller networks are much more dense, and therefore are actually less
efficient than if they used standard matrix multiplication methods. To be
perfectly optimized, therefore, the algorithm would use standard matrix
multiplication for small networks, and switch to sparse methods for larger
networks. However, the gains would be minimal, generally under 10 seconds, and
would require additional logical steps to check for network size, further
minimizing the gain. We therefore use sparse matrix methods for all network
sizes.

Figure 2: Runtime of the algorithm on networks ranging from size $10$ to
$10,000$ nodes in orders of magnitude, and from one to ten colors.
Additionally, the dashed line represents the computational time that would be
expected using standard matrix multiplication methods for $k=10$. These
runtimes were generated using virtual PCs including 1 dual core CPU and 10GBs
of RAM.

##  4 Empirical Use and Example

To show the empirical value of this algorithm, we use the Zachary karate club
social network [zachary1977information]. This is a well-known historical
network that describes the social relationships between 34 members of a
university karate club. Ties exist between members if they overlapped in at
least one of eight contexts representing undirected relations. These relations
varied in terms of likely strength of the association. Likely at the weak end
of the spectrum is being enrolled in the same class at the university, while
likely at the strong end is being a student-teacher at the studio.
Additionally, three ties are specific to activities with a part-time
instructor. Member “factions” were identified as a node attribute, taking one
of five mutually exclusive values: strongly associated with the president,
weakly associated with the president, neutral, weakly associated with the
part-time instructor, or strongly associated with the part-time instructor.
These are labeled ”Zs”, ”Zw”, ”N”, ”Hw”, and ”Hs”, respectively. These labels
can be placed on an ordinal scale from -2 (Zs) to 2 (Hs) to quantify members’
direction and strength of alignment. This undirected network with five colors
represents a case that is rich in the number of colored triads (220) for
detailed conclusions to be drawn using the proposed algorithm (which is
general to both undirected and directed networks).

We initially ran the colored triad census on the social network using the
faction as the nodal attribute. This provided the basis for our empirical
observed colored triad census. To determine whether these triads were observed
more or less often than expected by chance, we constructed a null model. As
the choice of null model can have important ramifications for the null
distribution of triads, we chose a model where edge formation is a function of
the probability of ties between nodes of specific attributes
[faust2010puzzle]. The null model is a mixing-matrix conditioned uniform
random graph distribution based on probabilities of edges between nodes of
particular color combinations [newman2003mixing]. This matrix comprises
empirical probabilities of ties between groups, with the diagonal representing
within-group tie probabilities. Observations of significantly over- or under-
represented colored triads are the result of network effects beyond homophily
and heterophily. Networks are then generated from this matrix via a Bernoulli
random graph process á la ER1959. This null model therefore conditions on
graph size, the distribution of node factions, and the probability of ties
within and between factions. By generating networks from the null model, we
can observe whether colored triad counts deviate from that expected based on
the marginal distribution of faction mixing. Because we condition on the above
parameters, if we observe statistical deviations in our colored triad census,
it indicates that the structure of the network is dependent on parameters
other than those on which we conditioned.

Moreover, for any triad, the expected number and variance can be calculated
assuming each tie follows a Binomial distribution (which is a reasonable
assumption for most binary social network data). The observed number can then
be compared to these numerical results and a p-value extracted from an exact
Binomial test. This equates to the following probability, expectation, and
variance for an example colored triad:

$\displaystyle\begin{split}P(T)=&P(A_{ij}=1|R(i)=r_{1},R(j)=r_{2})\times
P(A_{ij}=1|R(i)=r_{2},R(j)=r_{3})\\\
\times&P(A_{ij}=1|R(i)=r_{3},R(j)=r_{1})\end{split}$ (3)
$E(T)=P(T)\times\prod_{r=1}^{L(T)}{{\sum K^{r}_{\bullet 1}}\choose{S(T,r)}}$
(4) $V(T)=E(T)\times(1-P(T))$ (5)

The probability of $T$, $P(T)$ in Equation 3 is based on the mixing-matrix of
the three colors ($r$) involved in the triad $T$. As is standard for the
mixing-matrix approach, this continues to assume that all edges in the graph
are independent. For the expected value of a specific triad, we multiply the
probability of a single one of those triads by the total number of colored
triplets that exist in the graph. In Equation 4, the expectation of the triad,
$L(T)$ returns the number of unique colors in $T$ and $\sum K^{r}_{\bullet 1}$
is the number of nodes of color $r$ in the graph. Also, we take the nodes one,
two, or three at a time depending on how many times that color repeats in $T$,
represented by $S(T,r)$. This expectation therefore follows a binomial
distribution, and it’s variance follows accordingly in Equation 5.

However, to show that this method also works for null distributions that are
not analytically solvable, we construct a null distribution based on simulated
draws from the null model. As the number of trials increases, the simulated
null distribution of the colored triad census should asymptotically approach
the analytical solution shown above. For each of $1,000$ trials, we draw
random networks from the null distribution, and run the triad census on all
these networks. Comparing our observed count to the null distribution then
allows us to get an approximate p-value for a conditional uniform graph test,
and test the over- or under-representation of each colored triad. We now turn
to these results.

###  4.1 Results

Figure 3 is a heatmap of the approximate p-values associated with each
binomial exact test against the null for each triad, clustered by the triad
and the colored triplet as returned by the proposed algorithm. We use a
clustering algorithm to group color triplets with similar profiles across the
types of triads. This assists with identifying trends across different colored
triads, leading to conclusions that would likely be missed if all the colored
triads were individually examined. We find particular importance in three
branch cutpoints in the clustering algorithm on the color triplets. The first
branch in the clustering algorithm (A in Figure 3) separates four color
triplets, comprising 16 colored triads, with a pattern of over-observed 003
and 102 triads, and under-observed 201 and 300 triads. These results show that
these color triplets are those that are less clustered than expected by
chance. The color triplets all contain nodes of two factions with the first
two nodes being $Hs$, that is, those strongly aligned with the part-time
instructor. This indicates that those who are so aligned are likely to form
ties to one another, but not to members of other factions. The only exception
in this group is that two $Hs$ nodes are more likely to form a tie from one of
the $Hs$ members to a $Hw$ member, but even in this case the complete triad
(003) is still observed less than expected by chance. This particular result
is, perhaps, unsurprising, since $Hs$ and $Hw$ members are close in alignment,
more so than with those aligning with the president. Therefore, given the
tendency towards homophily they are likely to overlap, though less strongly
than members of the same faction; hence, the under-observed $T_{003-HsHsHw}$.

Figure 3: Heatmap of colored triads and their corresponding p-value of how
often they were observed in the empirical networks relative to the null
distribution. The columns separate triads based on the MAN configuration, and
the rows separate triads based on the triplet of colors. Standard clustering
algorithms were used to create the dendrograms. White space indicates
redundant isomorphism classes. Gray boxes are either those with 0 triads
observed in the network or in any of the networks of the null distribution,
and therefore have an undefined pseudo p-value, or those with a pseudo p-value
of $0.5$. The three labels correspond to three breakpoints in the clustering
that separate meaningful groups. (A) is a group of four color triplets
exhibiting homophily between $Hs$ nodes. (B) is a group of 21 colored triplets
exhibiting low clustering between heterogeneous nodes. (C) is a group of 6
colored triplets that show potential significant amounts of bridging.

The second branching point in the clustering (B in Figure 3) separates the
group of color triplets that are over-observed for the 003 triad, under-
observed for the 102 and 201 triads, and observed about as much as expected
for the 300 triads. All the triplets in question have nodes of different
factions in the first and second position. Because the edge in the 102 triad
is between the first and second node in the triplet (Figure 1), this means
that these are all triplets where the first edge is less likely than expected
by chance, and the lack of formation of the first edge subsequently hampers
the formation of the edge between the second and third nodes in the triplet
(201 triad). The first two nodes of these triplets are often (e.g., 16 out of
21) two factions at least a distance of two away (e.g. $N$ and $Hs$),
indicating members of a faction are not likely to overlap with members who are
too disparate from their faction. Put another way, this pattern of triads
shows a lack of faction heterophily.

The third branch point (unlabeled) is primarily singling out the group of
color triplets that were not observed in the network, and we cannot draw
conclusions about their prevalence. The fourth branch point (C in Figure 3),
however, distinguishes a group of five triplets that are under-observed for
the 102 triad and over-observed for the 201 triad. This means that the edge
between the first two nodes is less likely than expected by chance, but once
that edge does occur, the second edge occurs more often than expected by
chance. All these triplets begin with a $Zs$ member, and the 201 triad in this
case is effectively a bridging tie between it and another. Interestingly, the
bridging node is anything other than an $Hs$ (whom are primarily consigned to
this role in branch $A$, as discussed above). The third node was another $Zs$
member in four of five triplets. This indicates that $Zs$ members of the
karate club did not often overlap members of other factions, but when they
did, provided it was not with an $Hs$, that second person also often
overlapped with another $Zs$.

Although the above examples show homophily and bridging, analyzing the full
colored triad census allows us to draw further conclusions by looking at other
colored triads. In particular, the homophily has mostly been a story of the
$Hs$ nodes, and the bridging primarily about the $Zs$ nodes. The 300 triad of
both of these factions, when comprising three nodes of the same faction, are
observed more often than expected by chance in both cases, which has different
implications on the previously-noted results. For the $Hs$ nodes, homophily is
strengthened, as not only do $Hs$ nodes not often overlap with members of
other factions, they also very strongly overlap with one another. This may
partially be an artifact of the types of overlap, as stated before, three of
the overlap activities involve direct participation in the part-time
instructor’s studio, but there are no corresponding groups for the president.
This means that those who are $Hs$ or $Hw$ may have more opportunity to
overlap with one another due solely to the structure of the data. On the other
hand, the triplet of all $Zs$ members also has an over-observed 300 triad.
Although there are other triads that seem to indicate bridging between $Zs$
members (C in Figure 3), given that $Zs$ members are also densely connected to
one-another, the practical effect of these potential bridging ties is reduced.
Observing this joint effect of homophily and bridging ties was possible only
through the complete colored triad census. Neither a standard triad census nor
a brokerage analysis would have revealed the intricacies of these results.

In sum, it is clear from these results that the colored triad census allows
one to examine multiple trends simultaneously that are often done in isolated
analyses, including homophily, heterophily, and brokerage. Importantly, it
also allows for generalizations based on the clustering of various triads or
color triplets, as well as specific results based on individual triads. In
this manner, the colored triad census can yield results on multiple structural
levels simultaneously, all while examining local structure, nodal attributes,
and their interaction—that is, net of all alternatives involving mixtures of
node coloring and triadic configurations.

##  5 Limitations

There are some limitations to this method. First, it is only computationally
efficient relative to existing methods (including brute force counting).
Networks of $10,000$ nodes or more will take over a day to run using the
proposed algorithm for the colored triad census. However, this is an easily
parallelizable process (by partitioning the separate algebraic steps, for
example), and so the real time necessary to run the analysis can be greatly
reduced by taking advantage of this feature. The time needed for the
parallelized colored triad census is approximately inversely proportional to
the number of computational cores used in the calculation (plus some
overhead). Second, the interpretation and visualization of these results is
complicated, particularly as the number of colors increases. Examining all of
the triads simultaneously reduces the likelihood of missing interesting
results because a specific colored triad was excluded. However, the sheer
number of colored triads means that making complete sense of results can be
difficult. Even if the results are carefully examined for all colored triads,
it is conceivable that one might miss an important result out of the $11,080$
colored triads in a directed, 10-color network, no matter how meticulous the
examiner’s eye. However, use of standard clustering algorithms and heatmaps
(as in Figure 3) may help to ease interpretation of the results at both a
coarse- (general groups of triads) or fine-grained (individual colored triads)
perspective. That said, we recognize that the interpretation is not
straightforward and that this is a first effort at understanding these
results, but we believe that having an algorithm to efficiently calculate the
colored triad census will spur additional work towards interpreting and using
the results. As a result better approaches therein will emerge with time and
use.

##  6 Conclusions

In this paper, we have extended the matrix algebra methods of Moody1998 to
calculate the colored triad census for any network, directed or undirected,
with an arbitrary number of colors in a relatively computationally efficient
manner. We have shown a number of mathematical results regarding the colored
triad census, including a generalized equation for an arbitrary colored triad,
the number of isomorphism classes for arbitrary numbers of colors, and the
expectation and variances for colored triads. We analyzed an empirical social
network using our algorithm, and calculated approximate p-values for each
colored triad, based on an analytic exact binomial test for less complex null
distributions, or approximately through simulation for more complex null
distributions. We have also shown the type of conclusions that can be drawn
from these results, observing results that would not be feasible with many
other currently available methods.

One additional benefit of this method is that it can be directly used as a
counting tool for sufficient statistics in network inference models, such as
exponential random graphs (ERGM). The colored triad census essentially allows
one to simultaneously evaluate the effect of local structure and node
attribute on network structure in an ERGM, building off previous work where
researchers explicated the ERGMs capacity for including the triad census
[Yaveroglu2015]. We believe that the colored triad census is a useful
technique with an efficient implementation that can be widely-applicable in
social networks research, showing the continued importance of the triad census
even in this era of stochastic models for complex networks.

##  7 Acknowledgements

The authors would like to thank the two anonymous reviewers and the editors
for their contributions to this manuscript. This work utilized the
computational resources of the NIH HPC Biowulf cluster (http://hpc.nih.gov).
This research was funded via National Human Genome Research Institute,
National Institutes of Health (Grant number ZIA HG200335) and the Oxford
Martin School, University of Oxford (Grant number LC1213-006).

##  8 References

##  Appendix A Variable and Functional Definitions

Variable or function notation |  Description of variable or function
---|---
$A$ | Adjacency matrix
$E$ | Symmetrized adjacency matrix
$\bar{E}$ | Complement of symmetrized adjacency matrix
$M$ | Adjacency matrix including only mutual ties
$C$ | Adjacency matrix including only asymmetric ties
$K^{r}$ | Coloring matrix for color $r$
$R(i)$ | Function returning the color of node $i$
$H(T,i,j)$ | Function returning the matrix of the edge in triad T
| between nodes i and j
$T$ | An arbitrary colored triad, with a MAN configuration,
| and colored triplet $r_{1},r_{2},r_{3}$
$L(T)$ | A function returning the number of unique colors for
| a given colored triad
$S(T,r)$ | Function returning the number of times color $r$
| appears in colored triad $T$
$P(T)$ | The probability of observing triad $T$
$E(T)$ | The expectation of triad $T$ under a binomial model
$V(T)$ | The variance of triad $T$ under a binomial model
Table 3: List of variables, constants, and functions defined in this
manuscript.