# Spatial modeling of dyadic genetic relatedness data: Identifying factors
associated with M. tuberculosis transmission in Moldova

Joshua L. Warren1, Melanie H. Chitwood2, Benjamin Sobkowiak3,
Valeriu Crudu4, Caroline Colijn3, Ted Cohen2 This work was supported in part
by the National Institute of Allergy and Infectious Diseases (R01 AI147854,
R01 AI137093) and the generous support of the American people through the
United States Agency for International Development (USAID) through the TREAT
TB Cooperative Agreement No. GHN-A-00-08-00004. The contents are the
responsibility of the authors and Subgrantee and do not necessarily reflect
the views of USAID or the United States Government.

(1Department of Biostatistics, Yale University, New Haven, CT 06510, USA
2Department of Epidemiology of Microbial Diseases, Yale University, New Haven,
CT 06510, USA
3Department of Mathematics, Simon Fraser University, Burnaby, BC V5A 1S6, CA
4Microbiology and Morphology Laboratory, Phthisiopneumology Institute,
Chisinau, MLD
)

###### Abstract

Understanding factors that contribute to the increased likelihood of disease
transmission between two individuals is important for infection control.
However, analyzing measures of genetic relatedness is complicated due to
correlation arising from the presence of the same individual across multiple
dyadic outcomes, potential spatial correlation caused by unmeasured
transmission dynamics, and the distinctive distributional characteristics of
some of the outcomes. We develop two novel hierarchical Bayesian spatial
methods for analyzing dyadic genetic relatedness data, in the form of
patristic distances and transmission probabilities, that simultaneously
address each of these complications. Using individual-level spatially
correlated random effect parameters, we account for multiple sources of
correlation between the outcomes as well as other important features of their
distribution. Through simulation, we show the limitations of existing
approaches in terms of estimating key associations of interest, and the
ability of the new methodology to correct for these issues across datasets
with different levels of correlation. All methods are applied to
_Mycobacterium tuberculosis_ data from the Republic of Moldova where we
identify previously unknown factors associated with disease transmission and,
through analysis of the random effect parameters, key individuals and areas
with increased transmission activity. Model comparisons show the importance of
the new methodology in this setting. The methods are implemented in the R
package GenePair.

Keywords: Hierarchical Bayesian methods; Patristic distance; Phylogenetic
tree; Spatial statistics; Transmission probability.

##  1 Introduction

The availability and affordability of whole-genome sequencing technology has
made the use of genomic data increasingly common in epidemiological modeling
studies (Loman and Pallen, 2015; Polonsky et al., 2019). Sequencing data from
bacteriological pathogens provide insight into the spread of infectious
diseases through populations and numerous methods have been developed to infer
individual-to-individual transmission. The simplest approach uses the
difference of single nucleotide polymorphisms (SNPs) between two samples and
assigns any pairs separated by fewer than a threshold number of SNPs to a
transmission cluster (Walker et al., 2014). Other approaches make use of
phylogenetic trees, a reconstruction of the estimated evolutionary history of
a pathogen based on genetic sequencing data (Delsuc et al., 2005). For
example, patristic distances, the summed length of the tree branches
separating two isolates, is a similar measure of genetic distance, but this
approach leverages prior biological assumptions about the substitution rate to
generate a more accurate measure than SNP distance (Poon, 2016). More complex
approaches make use of transmission trees, based on phylogenetic analyses, to
estimate the asymmetric probability of direct transmission between pairs of
cases (Didelot et al., 2017; Klinkenberg et al., 2017; Campbell et al., 2019).

Understanding how potentially modifiable factors contribute to the increased
likelihood of disease transmission between two individuals is important for
infection control. However, while a robust set of methods exists to estimate
novel measures of genetic relatedness between case pairs, advanced statistical
methods to analyze these dyadic data and their associations with other factors
may require modifications/extensions due to the unique features of the
outcomes (e.g., spatial correlation, distributional characteristics).

Fixed effect regression-based approaches that do not account for correlation
between the dyadic outcomes may have serious inferential limitations. For
example, we may expect correlation between dyadic outcomes $Y_{ij}$ and
$Y_{kj}$, describing the genetic relatedness between individuals $i,j$ and
$k,j$ respectively, because individual $j$ is present in both pairs;
particularly if individual $j$ is a major driver of transmission in the
population. Ignoring this potentially positive correlation during analysis can
result in overly optimistic measures of uncertainty for regression parameter
estimates, leading to incorrect conclusions regarding statistical significance
(Hoff, 2003).

Several different analysis frameworks have been developed for modeling network
dependence in dyadic data (Kenny et al., 2020). Random effect regression
models have been widely used in this setting due to their ability to flexibly
characterize the correlation, ability to handle different outcome types, and
relative ease in making statistical inference (Warner et al., 1979; Wong,
1982; Gill and Swartz, 2001; Hoff et al., 2002; Hoff, 2003, 2005, 2021).
However, current models often ignore other important sources of correlation
that may be unique to the infectious disease setting (e.g., spatial
correlation due to unmeasured transmission dynamics) and/or are limited in
their ability to describe non-standard distributions that may be seen when
analyzing novel genetic relatedness outcomes (e.g., zero-inflation), with a
few exceptions. Beck et al. (2006) and Neumayer and Plümper (2010) discuss the
use of spatial lag regression models for analyzing spatially-referenced dyadic
data while Austin et al. (2013) and Ciminelli et al. (2019) integrate
spatially correlated random effects within a network analysis by extending the
latent position distance model of Hoff et al. (2002). With respect to non-
standard distributions, Simpson and Laurienti (2015) introduce a non-spatial
(in terms of the random effect parameters) two-part mixed model for analyzing
the probability of connection in whole-brain network dyadic data.

In this work, we develop hierarchical Bayesian spatial methods for analyzing
dyadic genetic relatedness data in the form of patristic distances and
transmission probabilities by extending the existing random effect frameworks
to better reflect features of the genetic relatedness outcomes. Specifically,
instead of inducing spatial correlation in the data through the latent
position distance components of the model as in Austin et al. (2013) and
Ciminelli et al. (2019), we incorporate spatial structure directly into the
individual-level random effect parameters. We fully investigate the
implications of this choice on the induced correlation structure of the data.
Additionally, we accommodate the distinctive distributional features of the
considered genetic relatedness outcomes during modeling. Finally, we develop
efficient Markov chain Monte Carlo (MCMC) sampling algorithms for several
genetic relatedness outcomes, along with an R package, to facilitate posterior
sampling in future work.

Using a simulation study, we show the importance of these methods for
conducting accurate statistical inference for key regression associations and
the limitations of existing approaches. We also apply our methods to a unique
dataset of _Mycobacterium tuberculosis_ isolates and associated patient data
from the Republic of Moldova and show the benefits of the new methodology with
respect to model fit and uncovering new insights into potentially important
transmission factors. Analyzing the posterior distributions of individual-
specific random effect parameters is shown to be important for understanding
how individuals in the population personally contribute to the transmission
dynamics and the role of spatial versus individual variability in this
process. Given the increasing availability and interest in these types of
data, and the limitations of existing approaches, we anticipate that these
methods will represent important tools for researchers looking to correctly
identify factors associated with genetic relatedness between individuals in
future studies.

In Section 2 we describe the data from the Republic of Moldova while the new
statistical methods are presented in Section 3. Sections 4 and 5 represent the
simulation study and real data application, respectively. We close in Section
6 with conclusions and further discussion.

##  2 Motivating data and application

In this study, we analyze data previously described by Yang et al. (2022). In
that work, the authors performed whole genome sequencing on _Mycobacterium
tuberculosis_ isolates from 2,236 of the 2,770 non-incarcerated adults
diagnosed with culture-positive tuberculosis (TB) in the Republic of Moldova
between January 1, 2018 and December 31, 2019. They constructed a maximum
likelihood phylogenetic tree with RAxML (Stamatakis, 2014) and identified
broad, putative transmission clusters using TreeCluster (Balaban et al., 2019)
with a threshold of 0.001 substitutions per site.

Table 1: Summary of the dyadic genetic relatedness data in the Republic of Moldova study population by transmission probability (TP) status. Means, with interquartile ranges given in parentheses, are shown for continuous variables and percentages are shown for categorical variables. Effect | TP $=0$ ($n=6,958)$ | TP $>0$ ($n=2,744$)
---|---|---
Distance Between Villages (km) | 87.49 (67.92) | 94.76 (78.13)
Same Village (% Yes) | 0.56 | 1.28
Date of Diagnosis Difference (Day) | 225.09 (233.75) | 202.41 (211.00)
Age Difference (Year) | 13.57 (14.00) | 13.33 (15.00)
Sex (%): |  |
Both Male | 56.54 | 58.89
Both Female | 6.11 | 4.63
Mixed Pair | 37.35 | 36.48
Residence Location (%): |  |
Both Urban | 14.62 | 14.18
Both Rural | 37.55 | 38.16
Mixed Pair | 47.83 | 47.67
Working Status (%): |  |
Both Employed | 0.99 | 0.77
Both Unemployed | 80.81 | 80.50
Mixed Pair | 18.19 | 18.73
Education (%): |  |
Both $<$ Secondary | 10.30 | 10.02
Both $\geq$ Secondary | 45.33 | 46.21
Mixed Pair | 44.37 | 43.77

From these analyses, Yang et al. (2022) produced estimates of genetic
relatedness among sequences using two metrics. First, they computed patristic
distance between any pair of isolates within a cluster. Patristic distance is
a measure of genetic relatedness between two sequences in a phylogenetic tree
expressed in substitutions per site. Second, the authors estimated the
probability that one individual infected another individual. These
probabilities are obtained using TransPhylo (Didelot et al., 2017), a Bayesian
approach that augments timed phylogenetic trees with who-infected-whom
information while accounting for the possibility of unsampled individuals.
TransPhylo uses an MCMC framework to obtain a posterior collection of
transmission trees, accounting for the time from infection to infecting others
(generation time) and the time from infection to sampling. For most pairs
$(i,j)$, there is no posterior transmission tree that includes a transmission
event from $i$ to $j$ or from $j$ to $i$. Where there are samples in the
posterior in which $i$ or $j$ infected the other, the posterior probability is
not necessarily symmetric. This can occur, for example, because $i$ was
sampled many months prior to $j$, making it more likely that $i$ infected $j$
than vice versa. For all possible pairs within a transmission cluster,
symmetric estimates of patristic distance (i.e., $Y_{ij}=Y_{ji}$) and
asymmetric estimates of transmission probability (i.e., $Y_{ij}\neq Y_{ji}$
necessarily) are available. To our knowledge, this is the first study to
analyze transmission probabilities in the infectious disease setting.

In addition to estimates of genetic relatedness, demographic data are
available for each TB case. These data include individual characteristics such
as age (in years), sex, education status (less than secondary, secondary or
higher), working status (employed, unemployed), and residence type (urban, not
urban). With these data, we calculate characteristics of the pair of
individuals, such as an indicator for whether the individuals in the pair
reside in the same village, the Euclidean distance between their villages of
residence (in kilometers), the difference between their dates of diagnosis (in
days), and the absolute difference between their ages (in years).

Figure 1: Patristic distances (substitutions per site; panel 1), transmission
probabilities that are $>0$ (panel 2), and village locations (panel 3) from
the largest putative cluster in the Republic of Moldova data analysis. In
panel 3, gray, blue, orange, and red points represent villages with one, two,
three, and five cases, respectively.

In this work, we analyze data from individuals in the largest putative
transmission cluster. After removing six individuals with missing covariates,
this cluster includes 99 individuals resulting in 4,851 symmetric patristic
distance pairs and 9,702 asymmetric transmission probability pairs.
Characteristics of the dyadic data are shown in Table 1 while in Figure 1 we
display distributions of the different genetic relatedness outcomes along with
the village locations for the included individuals.

##  3 Methods

We develop models to analyze spatially-referenced dyadic genetic relatedness
outcomes (i.e., patristic distances, transmission probabilities) while
accounting for multiple sources of correlation between responses and important
features of the outcomes. The methods are designed to yield accurate
statistical inference for the primary regression associations of interest
while also identifying individuals/locations that play a more important role
in the transmission dynamics within the population.

###  3.1 Patristic distances

We model the patristic distance (log scale) between individuals $i$ and $j$ as
a function of pair- and individual-level covariates as well as spatially-
referenced, individual-specific random effect parameters such that

$\ln\left(P_{ij}\right)=\textbf{x}_{ij}^{\text{T}}\boldsymbol{\beta}+\left(\textbf{d}_{i}+\textbf{d}_{j}\right)^{\text{T}}\boldsymbol{\gamma}+\theta_{i}+\theta_{j}+\epsilon_{ij},$
(1)

$i=1,\ldots,n-1,\ j=i+1,\ldots,n,$ where
$\epsilon_{ij}|\sigma^{2}_{\epsilon}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,\sigma^{2}_{\epsilon}\right)$
is the error term; $n$ is the total number of individuals in the study;
$P_{ij}>0$ is the symmetric patristic distance between a unique pair of
individuals $i$ and $j$ (i.e., $P_{ij}=P_{ji}$ for all $i\neq j$);
$\textbf{x}_{ij}$ is a vector of covariates describing differences between
individuals $i$ and $j$ (e.g., spatial distance), with $\boldsymbol{\beta}$
the corresponding vector of regression parameters; $\textbf{d}_{i}$ is a
vector of covariates specific to individual $i$ where the impact of the
covariates on the response, described by the $\boldsymbol{\gamma}$ vector, is
assumed to be the same across all individuals; and $\theta_{i}$ is a
spatially-referenced, individual-specific random effect parameter which
describes individual $i$’s role in transmission events within the population
(both as infector and infectee). Small $\theta_{i}$ values indicate that
across all outcome pairs involving individual $i$, the patristic distance is
smaller on average, suggestive of an increased likelihood of being in
transmission pairs.

We allow for the possibility of spatial correlation between responses by
modeling the random effect parameters using a spatially-referenced Gaussian
process such that

$\displaystyle\begin{split}&\theta_{i}=\eta\left\\{h\left(\textbf{s}_{i}\right)\right\\}+\zeta_{i},\
i=1,\ldots,n,\\\
&\boldsymbol{\eta}^{\text{T}}=\left\\{\eta\left(\textbf{s}_{1}^{*}\right),\ldots,\eta\left(\textbf{s}_{m}^{*}\right)\right\\}|\phi,\tau^{2}\sim\text{MVN}\left\\{\boldsymbol{0}_{m},\tau^{2}\Sigma\left(\phi\right)\right\\},\text{
and}\\\
&\Sigma\left(\phi\right)_{ij}=\text{Corr}\left\\{\eta\left(\textbf{s}_{i}^{*}\right),\eta\left(\textbf{s}_{j}^{*}\right)\right\\}=\exp\left\\{-\phi\left\|\textbf{s}_{i}^{*}-\textbf{s}_{j}^{*}\right\|\right\\}\end{split}$
(2)

where $\theta_{i}$ is decomposed into two pieces; one that is purely spatial,
$\eta\left(\textbf{s}_{i}^{*}\right)$, and one that is non-spatial,
$\zeta_{i}$. This allows for each parameter to be individual-specific and not
driven solely by spatial location. In the case where individuals are co-
located, the function $h\left(.\right)$ maps the spatial location of an
individual to an entry within a smaller set of $m<n$ unique locations such
that
$h\left(\textbf{s}_{i}\right)\in\left\\{\textbf{s}_{1}^{*},\ldots,\textbf{s}_{m}^{*}\right\\}$,
where it is possible that
$h\left(\textbf{s}_{i}\right)=h\left(\textbf{s}_{j}\right)=\textbf{s}^{*}_{k}$
for some $i\neq j$. When all individuals have a unique location,
$h\left(\textbf{s}_{i}\right)=\textbf{s}^{*}_{i}$ for all $i$ and therefore,
$m=n$.

The vector of purely spatial random effect parameters, $\boldsymbol{\eta}$, is
modeled using a Gaussian process centered at zero (i.e., $\textbf{0}_{m}$ is
an $m$ length vector of zeros) with correlation structure defined by the
Euclidean distances between spatial locations (i.e.,
$\left\|\textbf{s}^{*}_{i}-\textbf{s}^{*}_{j}\right\|$), where $\phi>0$
controls the level of spatial correlation between the parameters and
$\tau^{2}$ the total variability of the spatial process. Small values of
$\phi$ indicate strong spatial correlation even at larger distances. Finally,
$\zeta_{i}|\sigma^{2}_{\zeta}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,\sigma^{2}_{\zeta}\right)$
represent the individual-specific parameters that account for the possibility
that two people could be in very similar spatial locations but have vastly
different patterns of behavior that impact their likelihood of being
transmitted to and/or transmitting to others.

####  3.1.1 Prior distributions

We assign weakly informative prior distributions to the remaining model
parameters when appropriate. The regression parameters are specified as
$\beta_{j},\gamma_{k}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,100^{2}\right)$
for $j=1,\ldots,p_{x}$ and $k=1,\ldots,p_{d}$, where $p_{x}$ and $p_{d}$ are
the lengths of the $\textbf{x}_{ij}$ and $\textbf{d}_{i}$ vectors,
respectively; the variance parameters as
$\sigma^{2}_{\epsilon},\tau^{2},\sigma^{2}_{\zeta}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{Inverse
Gamma}\left(0.01,0.01\right)$; and the spatial correlation parameter as
$\phi\sim\text{Gamma}\left(1.00,1.00\right)$. We scale the spatial distances
used in the analysis to allow the prior distribution for $\phi$ to be
minimally informative (i.e., ranging from relatively weak to strong spatial
correlation at both short and long spatial distances).

###  3.2 Transmission probabilities

Next, we introduce a method for analyzing transmission probabilities which,
unlike patristic distances, contain potentially rich information regarding the
direction of transmission. As a result, the outcomes are not necessarily
symmetric as they were in Section 3.1 (i.e., $T_{ij}\neq T_{ji}$ necessarily).
Here, we introduce a framework for analyzing transmission probabilities that
includes similar inferential goals as the model in (1, 2) while also
accounting for the asymmetry of the outcome as well as other important data
features (e.g., zero-inflation).

Specifically, we model the probability that individual $j$ infected individual
$i$ (i.e., $T_{ij}$) as a function of individual- and pair-specific covariates
and individual/location-specific random effect parameters using a mixed-type
distribution. This specification includes a binary component to account for
the large proportion of exact zero transmission probabilities (see Figure 1),
and a continuous piece to model the non-zero probabilities. We define the
probability density function (pdf) for a transmission probability,
$f_{t_{ij}}\left(t\right)$, as

$T_{ij}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}f_{t_{ij}}\left(t\right)=\left(1-\pi_{ij}\right)^{1\left(t=0\right)}\left[\frac{\pi_{ij}}{t\left(1-t\right)}f_{w_{ij}}\left\\{\ln\left(\frac{t}{1-t}\right)\right\\}\right]^{1\left(t>0\right)},\
t\in\left[0,1\right)$ (3)

where all pairs are now included in the analysis (i.e., $i=1,\ldots,n$;
$j=1,\ldots,n$; $i\neq j$); $1\left(.\right)$ is an indicator function taking
the value of one when the input statement is true and the value of zero
otherwise; the transmission probabilities can be exactly equal to zero; and
$\left(1-\pi_{ij}\right)$ describes the probability that this event occurs.

We use a logistic regression framework to connect these underlying
probabilities with covariates and random effect parameters such that

$\ln\left(\frac{\pi_{ij}}{1-\pi_{ij}}\right)=\textbf{x}_{ij}^{\text{T}}\boldsymbol{\beta}_{z}+\textbf{d}_{j}^{\text{T}}\boldsymbol{\gamma}_{z}^{\left(g\right)}+\textbf{d}_{i}^{\text{T}}\boldsymbol{\gamma}_{z}^{\left(r\right)}+\theta_{zj}^{\left(g\right)}+\theta_{zi}^{\left(r\right)}+\nu_{zi}\nu_{zj}$
(4)

where $\textbf{x}_{ij}$ and $\textbf{d}_{i}$ were previously described in
Section 3.1. Because we now observe two responses for each pair (i.e.,
$T_{ij}$ and $T_{ji}$), we are able to separate the included parameters into
the different roles of the individuals within a pair; specifically, we can
estimate the infector or “giver” ($g$) and the infectee or “receiver” ($r$)
terms. For example, $T_{ij}$ describes the probability that individual $j$
transmits to $i$, so $\textbf{d}_{j}^{\text{T}}\boldsymbol{\gamma}_{z}^{(g)}$
from (4) represents the impact of the giver’s (individual $j$’s) covariates on
the probability of transmission while
$\textbf{d}_{i}^{\text{T}}\boldsymbol{\gamma}_{z}^{(r)}$ describes the impact
of the receiver’s (individual $i$’s) characteristics. Similarly, now each
individual has two different additive random effect parameters,
$\boldsymbol{\theta}_{zi}^{\text{T}}=\left(\theta_{zi}^{(g)},\theta_{zi}^{(r)}\right)$;
one for the giver and receiver roles, respectively. The $\nu_{zi}\nu_{zj}$
interaction term accounts for correlation between observations from the same
pair of individuals in different roles (i.e.,
$\nu_{zi}\nu_{zj}=\nu_{zj}\nu_{zi}$), where
$\nu_{zi}|\sigma^{2}_{\nu_{z}}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,\sigma^{2}_{\nu_{z}}\right)$,
$i=1,\ldots,n$ (Hoff, 2005).

In order to understand if the magnitude of a non-zero transmission probability
is also impacted by covariates and random effect parameters, the pdf in (3)
includes a separate regression framework for the non-zero probabilities.
Specifically, $f_{w_{ij}}\left(w\right)$ is a pdf introduced to model the non-
zero probabilities on the logit scale such that

$f_{w_{ij}}\left(w\right)\equiv\text{N}\left(\textbf{x}_{ij}^{\text{T}}\boldsymbol{\beta}_{w}+\textbf{d}_{j}^{\text{T}}\boldsymbol{\gamma}_{w}^{\left(g\right)}+\textbf{d}_{i}^{\text{T}}\boldsymbol{\gamma}_{w}^{\left(r\right)}+\theta_{wj}^{\left(g\right)}+\theta_{wi}^{\left(r\right)}+\nu_{wi}\nu_{wj},\
\sigma^{2}_{\epsilon}\right)$ (5)

where $\sigma^{2}_{\epsilon}$ represents the variance of the distribution and
the remaining terms in (5) have been previously described. The “$w$”
subscripts in (5) serve to differentiate these parameters from those used by
the binary model in (4) (i.e., “$z$” subscripts).

In total, each individual in the analysis has four additive random effect
parameters,
$\boldsymbol{\theta}_{i}^{\text{T}}=\left(\boldsymbol{\theta}_{zi}^{\text{T}},\boldsymbol{\theta}_{wi}^{\text{T}}\right)$,
representing the residual (i.e., after adjustment for known risk factors and
repeated pair correlation) likelihood of transmitting ($g$) or being
transmitted to ($r$) in the binary ($z$) and positive ($w$) transmission
probability regressions. Large values of the $z$ subscript random effect
parameters suggest an increased chance of a non-zero transmission probability,
while large values of the $w$ subscript parameters indicate an increasingly
positive transmission probability.

As in (2), the introduced random effect parameters serve to adjust for pair-
and proximity-based correlation in the outcomes. Unlike in Simpson and
Laurienti (2015), we anticipate that the collection of parameters
corresponding to the same individual may themselves be correlated and
introduce a multivariate model as a result. The model for one set of the
parameters is similar to (2) and is given as (for $i=1,\ldots,n)$

$\displaystyle\theta_{zi}^{(g)}=\eta_{z}^{(g)}\left\\{h\left(\textbf{s}_{i}\right)\right\\}+\zeta_{zi}^{(g)}$

where
$\zeta_{zi}^{(g)}|\sigma^{2}_{\zeta^{(g)}_{z}}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,\sigma^{2}_{\zeta^{(g)}_{z}}\right)$
once again account for individual variability, and the remaining parameters,
$\theta^{(r)}_{zi}$, $\theta^{(g)}_{wi}$, $\theta^{(r)}_{wi}$, are defined
similarly. To account for cross-covariance and spatial correlation among the
parameters, we specify a multivariate Gaussian process for the mean parameters
such that

$\displaystyle\begin{split}&\boldsymbol{\eta}|\Omega,\phi\sim\text{MVN}\left\\{\boldsymbol{0}_{4m},\Sigma\left(\phi\right)\otimes\Omega\right\\}\text{
where}\\\
&\boldsymbol{\eta}^{\text{T}}=\left\\{\boldsymbol{\eta}\left(\textbf{s}^{*}_{1}\right)^{\text{T}},\ldots,\boldsymbol{\eta}\left(\textbf{s}^{*}_{m}\right)^{\text{T}}\right\\}\text{
and}\\\
&\boldsymbol{\eta}\left(\textbf{s}^{*}_{i}\right)^{\text{T}}=\left\\{\eta_{z}^{(g)}\left(\textbf{s}^{*}_{i}\right),\eta_{z}^{(r)}\left(\textbf{s}^{*}_{i}\right),\eta_{w}^{(g)}\left(\textbf{s}^{*}_{i}\right),\eta_{w}^{(r)}\left(\textbf{s}^{*}_{i}\right)\right\\}.\end{split}$
(6)

The complete collection of mean parameters across all $m$ unique spatial
locations is denoted by $\boldsymbol{\eta}$ (similar to (2));
$\boldsymbol{\eta}\left(\textbf{s}^{*}_{i}\right)$ represents the collection
of mean parameters across each component of the model and different roles,
specific to unique location $\textbf{s}_{i}^{*}$; $\Sigma\left(\phi\right)$
was previously described in (2); $\otimes$ is the Kronecker product; and
$\Omega$ represents a four-by-four unstructured covariance matrix describing
the cross-covariance among the set of four random effect parameters specific
to a unique spatial location.

####  3.2.1 Prior distributions

We complete the model specification by assigning prior distributions to the
introduced model parameters. As in Section 3.1.1, the regression parameters
are specified as $\beta_{zj}$, $\beta_{wj}$, $\gamma_{zk}^{(g)}$,
$\gamma_{zk}^{(r)}$, $\gamma_{wk}^{(g)}$, $\gamma_{wk}^{(r)}$
$\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{N}\left(0,100^{2}\right)$
for $j=1,\ldots,p_{x}$ and $k=1,\ldots,p_{d}$; the variance parameters as
$\sigma^{2}_{\epsilon}$, $\sigma^{2}_{\zeta_{z}^{(g)}}$,
$\sigma^{2}_{\zeta_{z}^{(r)}}$, $\sigma^{2}_{\zeta_{w}^{(g)}}$,
$\sigma^{2}_{\zeta_{w}^{(r)}}$, $\sigma^{2}_{\nu_{z}}$, $\sigma^{2}_{\nu_{w}}$
$\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{Inverse
Gamma}\left(0.01,0.01\right)$; the spatial correlation parameter as
$\phi\sim\text{Gamma}\left(1.00,1.00\right)$; and the cross-covariance matrix
as $\Omega^{-1}\sim\text{Wishart}\left(5,I_{4}\right)$, a weakly informative
choice resulting in uniform cross-correlations a priori (Gelman et al., 2013).

###  3.3 Induced correlation structure

The inclusion of individual-specific and spatially correlated random effect
parameters in the dyadic outcome models detailed in Sections 3.1 and 3.2
results in a positive correlation between the responses whose magnitude varies
depending on (i) if the two pairs share a common individual and (ii) the
geographic distance between the people in the pairs. To better understand the
induced correlations, we consider two different cases specifically for
patristic distances; when there is, and is not, a shared individual between
the pairs. Derivations for the transmission probability model are similar but
more complicated due to the mixed-type distribution and two regression
frameworks used. Therefore, in Figure S1 of the Supplementary Material, we
present simulation-based correlation estimates for this model over a range of
spatial variance/correlation settings; see Section S1 for full details.

First, we derive the correlation between two dyadic patristic distances that
consist of entirely different individuals, $\ln\left(P_{ij}\right)$ and
$\ln\left(P_{kl}\right)$ where $i\neq j\neq k\neq l$. The variance of one of
the responses (assuming every individual in the study has a unique spatial
location) is given as

$\displaystyle\begin{split}&\text{Var}\left\\{\ln\left(P_{ij}\right)\right\\}=\text{Var}\left(\theta_{i}\right)+\text{Var}\left(\theta_{j}\right)+\text{Var}\left(\epsilon_{ij}\right)+2\text{Cov}\left(\theta_{i},\theta_{j}\right)\\\
&=\text{Var}\left\\{\eta\left(\textbf{s}^{*}_{i}\right)+\zeta_{i}\right\\}+\text{Var}\left\\{\eta\left(\textbf{s}^{*}_{j}\right)+\zeta_{j}\right\\}+\text{Var}\left(\epsilon_{ij}\right)+2\text{Cov}\left\\{\eta\left(\textbf{s}^{*}_{i}\right)+\zeta_{i},\eta\left(\textbf{s}^{*}_{j}\right)+\zeta_{j}\right\\}\\\
&=2\tau^{2}\left(1+\exp\left\\{-\phi\left\|\textbf{s}^{*}_{i}-\textbf{s}^{*}_{j}\right\|\right\\}\right)+2\sigma^{2}_{\zeta}+\sigma^{2}_{\epsilon},\end{split}$

due to the spatial dependence between the $\theta_{i}$ parameters. The
covariance between these responses is given as

$\displaystyle\begin{split}&\text{Cov}\left\\{\ln\left(P_{ij}\right),\ln\left(P_{kl}\right)\right\\}=\text{E}\left\\{\left(\theta_{i}+\theta_{j}\right)\left(\theta_{k}+\theta_{l}\right)\right\\}\\\
&=\text{E}\left[\left\\{\eta\left(\textbf{s}^{*}_{i}\right)+\zeta_{i}+\eta\left(\textbf{s}^{*}_{j}\right)+\zeta_{j}\right\\}\left\\{\eta\left(\textbf{s}^{*}_{k}\right)+\zeta_{k}+\eta\left(\textbf{s}^{*}_{l}\right)+\zeta_{l}\right\\}\right]\\\
&=\tau^{2}\sum_{p_{1}\in\left\\{i,j\right\\}}\sum_{p_{2}\in\left\\{k,l\right\\}}\exp\left\\{-\phi\left\|\textbf{s}^{*}_{p_{1}}-\textbf{s}^{*}_{p_{2}}\right\|\right\\},\end{split}$

and is simply a function of spatial distances between the individuals in the
pairs. This suggests that if spatial correlation is negligible then the
covariance/correlation between outcome pairs without a shared individual is
effectively equal to zero since
$\exp\left\\{-\phi\left\|\textbf{s}^{*}_{i}-\textbf{s}^{*}_{j}\right\|\right\\}\approx
0$ for all $i,j$. In the case of strong spatial correlation (i.e,
$\exp\left\\{-\phi\left\|\textbf{s}^{*}_{i}-\textbf{s}^{*}_{j}\right\|\right\\}\approx
1$), the variance and covariance become
$4\tau^{2}+2\sigma^{2}_{\zeta}+\sigma^{2}_{\epsilon}$ and $4\tau^{2}$,
respectively. This yields a correlation between the dyadic outcomes of $0$ and
$\frac{4\tau^{2}}{4\tau^{2}+2\sigma^{2}_{\zeta}+\sigma^{2}_{\epsilon}}$ for
the weak and strong spatial correlation settings, respectively.

Following similar derivations, the covariance between observations from two
pairs that include one of the same individuals, say $\ln\left(P_{ij}\right)$
and $\ln\left(P_{ik}\right)$, is given as

$\tau^{2}\left(1+\exp\left\\{-\phi\left\|\textbf{s}^{*}_{i}-\textbf{s}^{*}_{j}\right\|\right\\}+\exp\left\\{-\phi\left\|\textbf{s}^{*}_{i}-\textbf{s}^{*}_{k}\right\|\right\\}+\exp\left\\{-\phi\left\|\textbf{s}^{*}_{j}-\textbf{s}^{*}_{k}\right\|\right\\}\right)+\sigma^{2}_{\zeta}.$

Under negligible spatial correlation, this expression approaches
$\tau^{2}+\sigma^{2}_{\zeta}$ due to the fact that the same individual is
represented in both pairs; while for strong spatial dependency, it approaches
$4\tau^{2}+\sigma^{2}_{\zeta}$. This yields correlations of
$\frac{\tau^{2}+\sigma^{2}_{\zeta}}{2\tau^{2}+2\sigma^{2}_{\zeta}+\sigma^{2}_{\epsilon}}$
and
$\frac{4\tau^{2}+\sigma^{2}_{\zeta}}{4\tau^{2}+2\sigma^{2}_{\zeta}+\sigma^{2}_{\epsilon}}$
for the weak and strong spatial correlation settings, respectively. Both of
these quantities are larger than the comparable versions for pairs that do not
share an individual.

These findings, along with the results shown in Figure S1 of the Supplementary
Material, suggest that there is increased correlation between outcomes when at
least one individual is shared between the pairs. In the transmission
probability model, there is also extremely high correlation between outcomes
when the pairs include both individuals but in different roles, with spatial
correlation/variability having almost no impact. However, even when pairs have
no shared individuals, strong correlation can still exist depending on the
strength/magnitude of spatial correlation/variability in the data. This may be
an important feature to consider in the infectious disease setting where
unmeasured transmission dynamics could results in residual spatial correlation
between individuals.

##  4 Simulation study

We design a simulation study to investigate the implications of ignoring
multiple forms of correlation between dyadic genetic relatedness outcomes when
making inference on regression parameters of interest. Additionally, we aim to
determine if a common Bayesian model comparison tool can be used to identify
datasets that contain non-negligible levels of correlation and also
differentiate the sources of the correlation. In order to avoid repetition of
text, the process described throughout Section 4 is specifically for the
patristic distances model. However, we carried out the same steps for
transmission probabilities using the corresponding equations from Section 3.2.
These full details are provided in Section S2 of the Supplementary Material.

###  4.1 Data generation

We simulate data from the model in (1) under three different scenarios. In
Setting 1, we specify that there is no unmeasured correlation by setting
$\theta_{i}=0$ for all $i$. In Setting 2, we simulate data from (1) and (2)
with $\eta\left(\textbf{s}_{i}^{*}\right)=0$ for all $i$, resulting in only
non-spatial variability in the random effect parameters. In Setting 3, we also
simulate data from (1) and (2) with
$\phi=-\ln\left(0.05\right)/\max\left\\{\left\|\textbf{s}_{i}^{*}-\textbf{s}_{j}^{*}\right\|;i<j\right\\}$
and $\zeta_{i}=0$ for all $i$, resulting in spatially correlated random effect
parameters with correlation that decreases to $0.05$ at the maximum distance
observed in the data (i.e., the effective range).

When simulating data from these models, we use results from our data
application in the Republic of Moldova (Section 5) to ensure that we are
working with realistic outcomes. Specifically, we use the same sample size
($n=99$), same covariates ($\textbf{x}_{ij}$, $\textbf{d}_{i}$), and choose
the true parameter values needed to simulate from (1) and (2) based on
posterior estimates obtained from the data application. In Table S1 of the
Supplementary Material, the specific values used in each simulation study are
given.

When simulating data from Settings 2 and 3, we define the total variance of
the random effect process as $\tau^{2}+\sigma^{2}_{\zeta}$ and use the values
in Table S1 of the Supplementary Material to estimate this quantity. In
Setting 2, all of this variability is assumed to be non-spatial while in
Setting 3 it is attributed entirely to spatial correlation. We generate a new
vector of $\boldsymbol{\theta}$ parameters for each dataset from the setting-
specific model. We also create a unique set of spatial locations, with no co-
located individuals (i.e., $m=n$), for each dataset. In total, we simulate 100
datasets from each setting.

###  4.2 Competing models

We apply three competing models to every dataset. The first model (i.e.,
Fixed) represents a simplified fixed effects regression form of (1) where
$\theta_{i}=0$ for all $i$ (i.e., no random effect parameters). As Fixed
matches the data generated from Setting 1, we expect it to perform well
overall in that setting. However, in Settings 2 and 3, it will likely struggle
in estimating the associations of interest as it ignores all correlation. The
second model (i.e., Non-spatial) also represents a variant of (1) and (2)
where $\eta\left(\textbf{s}^{*}_{i}\right)=0$ for all $i$ (matching data
generation Setting 2). Therefore, Non-spatial accounts for correlation due to
the nature of dyadic responses, but ignores the potential for spatial
correlation in the data. It is currently unknown how inference for the
regression parameters is impacted when dyadic correlation is accounted for but
spatial correlation is ignored. The final competing model (i.e., Spatial) is
our newly developed model in Section 3.1 which can address both sources of
correlation simultaneously.

We monitor several pieces of information collected from the analyzed datasets
to compare the methods. First, we calculate the mean absolute error (MAE) for
every regression parameter in the model using the posterior mean as the point
estimate. Next, we calculate 95% quantile-based equal-tailed credible
intervals (CIs) for each regression parameter and monitor how often this
interval includes the true value (ideally around 95% of the time) and its
length. Finally, we formally compare the models using Watanabe Akaike
information criterion (WAIC), a metric that balances model fit and complexity
where smaller values suggest that a model is preferred (Watanabe, 2010).

###  4.3 Results

We apply each method to each dataset and collect 20,000 posterior samples
after removing the first 5,000 iterations prior to convergence of the model.
Additionally, we thin the remaining samples by a factor of two to reduce
posterior autocorrelation, resulting in 10,000 samples with which to make
posterior inference. The priors from Sections 3.1.1 and 3.2.1 were used other
than when we encountered convergence issues when applying Spatial in the
transmission probabilities framework to a few of the datasets generated from
Setting 3. In that case, $\sigma^{2}_{\nu_{z}}\sim\text{Inverse
Gamma}\left(100.00,100.00\right)$ or $\text{Inverse
Gamma}\left(1000.00,1000.00\right)$ was used to stabilize estimation of the
$\nu_{zi}$ parameters.

The full results are shown in Table 2 for all models where we report the
average MAE, average empirical coverage (EC), and average CI length across all
regression parameters and simulated datasets, as well as the average
difference in WAIC values, with respect to Spatial, across all simulated
datasets.

Table 2: Simulation study results. $\Delta$ WAIC = WAIC $-$ WAIC${}_{\text{new}}$. Averages across the 100 simulated datasets are reported with standard errors given in parentheses. Bold entries indicate the “best” value within a setting and outcome type across models. MAE results for the patristic distances model are multiplied by 100 for presentation purposes. |  | Patristic |  | Trans. Probs.
---|---|---|---|---
Metric | Setting | Standard | Non-spatial | Spatial |  | Standard | Non-spatial | Spatial
MAE | 1 |  1.02 (0.03) | 1.03 (0.03) | 1.04 (0.03) |  | 1.62 (0.20) | 1.53 (0.19) |  1.40 (0.16)
| 2 | 6.77 (0.23) | 6.21 (0.23) |  6.19 (0.23) |  | 0.40 (0.00) | 0.35 (0.01) |  0.34 (0.01)
| 3 | 5.56 (0.20) | 4.84 (0.19) |  2.31 (0.08) |  | 0.43 (0.03) | 0.61 (0.24) |  0.35 (0.08)
EC | 1 |  0.94 (0.01) | 0.97 (0.00) | 0.97 (0.00) |  | 0.92 (0.01) | 0.92 (0.01) |  0.93 (0.01)
| 2 | 0.41 (0.02) | 0.93 (0.01) |  0.94 (0.01) |  | 0.44 (0.01) | 0.93 (0.00) |  0.94 (0.00)
| 3 | 0.39 (0.01) | 0.92 (0.01) |  0.96 (0.01) |  | 0.47 (0.01) | 0.90 (0.01) |  0.94 (0.01)
CI Length | 1 | 0.05 (0.00) | 0.06 (0.00) | 0.06 (0.00) |  | 4.33 (0.41) | 4.26 (0.40) | 4.28 (0.39)
| 2 | 0.10 (0.00) | 0.28 (0.00) | 0.30 (0.00) |  | 0.38 (0.00) | 1.61 (0.01) | 1.62 (0.01)
| 3 | 0.08 (0.00) | 0.21 (0.00) | 0.13 (0.00) |  | 0.52 (0.09) | 1.90 (0.48) | 1.31 (0.18)
$\Delta$ WAIC | 1 | -31.74 (0.97) | -2.95 (0.07) | – |  | 11.28 (1.95) | -0.49 (1.74) | –
| 2 | 6138.68 (56.90) | 0.57 (0.02) | – |  | 9013.90 (87.72) | 26.04 (22.47) | –
| 3 | 3983.08 (98.48) | 4.89 (0.10) | – |  | 8349.95 (298.65) | 111.87 (10.67) | –

As expected, all competing models perform similarly in terms of MAE, EC, and
CI length in Setting 1, with Spatial having improved MAE results for the
transmission probabilities framework. For the patristic distances model, WAIC
also correctly favors Fixed on average as it most closely matches the way the
data were generated. However, with respect to WAIC, Fixed is slightly
outperformed by Spatial in this setting for transmission probabilities.

In Settings 2 and 3 where correlation is present, Spatial consistently
outperforms Fixed and Non-spatial across all metrics. Fixed displays troubling
behavior in these settings, particularly with respect to EC. The 95% CIs are
only capturing the true parameter values around 39-47% of the times and are
much shorter on average than those from the competing methods. This suggests
that failing to account for correlation may result in CIs that are too narrow,
possibly leading to an inflated type I error rate.

The MAE results suggest that regression parameter inference may also suffer
when correlation is ignored and/or mischaracterized. For example, in Setting 3
where the data exhibit spatial correlation, Spatial produces substantially
smaller MAEs than Non-spatial (and Fixed), showing the importance of
accounting for spatial correlation. However, in Setting 2 where there is no
spatial correlation, Non-spatial and Spatial perform almost identically with
respect to MAE.

The WAIC results in Settings 2 and 3 are decisively in favor of Spatial over
Fixed, with differences between Spatial and Non-spatial increasing under
spatially correlated data. When combined with results from Setting 1, this
provides evidence to suggest that WAIC may be a useful tool for
differentiating datasets in terms of the presence and composition of
correlation.

##  5 _Mycobacterium tuberculosis_ in the Republic of Moldova

We apply the methods developed in Sections 3.1 and 3.2 to better understand
factors related to _Mycobacterium tuberculosis_ transmission dynamics in the
Republic of Moldova. Specifically, $\textbf{x}_{ij}$ from (1, 4, 5) includes
the pair-specific covariates described in Section 2 and Table 1 (i.e.,
intercept, same village indicator, distance between villages, difference in
diagnosis dates, difference in ages) while $\textbf{d}_{i}$ includes the
individual-specific covariates (i.e., age, sex, education, working status,
residence type). In addition to the new methods, we also present results from
the competing approaches detailed in Section 4.2 (i.e., Fixed and Non-
spatial). We use WAIC to compare the different models, identify the source of
correlation, and ultimately determine the need for the new methodology.

All models are fit in the Bayesian setting using MCMC sampling techniques,
with the full conditional distributions detailed in Section S3 of the
Supplementary Material. For each method, we collect 10,000 samples from the
joint posterior distribution after removing the first 50,000 iterations prior
to convergence and thinning the remaining 200,000 posterior samples by a
factor of 20 to reduce posterior autocorrelation. For each parameter,
convergence was assessed using Geweke’s diagnostic (Geweke, 1991) while
effective sample size was calculated to ensure we collected sufficient post-
convergence samples to make accurate statistical inference. Neither tool
suggested any issues of concern. We present posterior means and 95% quantile-
based equal-tailed CIs when discussing posterior inference.

###  5.1 Patristic distances

Results from the patristic distance analyses are shown in Figure 2 and Table
3. From Figure 2 it is clear that WAIC favors Spatial over the other
approaches. In the simulation study, WAIC was shown to consistently identify
the correct data generating setting, suggesting that there is likely non-
negligible and spatially structured correlation in these data. Consequently,
Fixed should not be used as it is likely to underestimate regression parameter
uncertainty, resulting in CIs that are often too narrow, while Non-spatial may
result in regression parameter estimates with inflated MAE. This can partly be
seen in Figure 2 where the Fixed CIs are much shorter on average. The point
estimates seen in Figure 2 are generally consistent across all methods in this
case.

Figure 2: Results from the patristic distance analyses in the Republic of
Moldova. Posterior means and 95% quantile-based equal-tailed credible
intervals are shown for the raw regression parameters from each of the
competing methods. Red lines indicate that the interval excludes zero. WAIC
values are provided along with the expected number of parameters/model
complexity term in parentheses (i.e., p${}_{\text{WAIC}}$), with smaller
values of WAIC preferred.

Posterior inference from Spatial for the exponentiated regression parameters
in Table 3 suggest that the indicator of whether the pair of individuals was
in the same village is the only association whose CI excludes one. Patristic
distances from pairs of individuals from the same village are around 45% (40%,
49%) smaller on average than those from pairs of individuals in different
villages.

Table 3: Results from the Spatial patristic distance analyses in the Republic of Moldova. Posterior means and 95% quantile-based equal-tailed credible intervals are shown for the exponentiated regression parameters (i.e., ratio of expected patristic distances per specified change in covariate value). The displayed credible intervals for the bolded entries exclude one. Effect | Estimate (95% CI)
---|---
Distance Between Villages (50 km) | 1.00 (0.99, 1.01)
Same Village (Yes vs. No) | 0.55 (0.51, 0.60)
Date of Diagnosis Difference (1/2 Year) | 1.01 (1.00, 1.01)
Age Difference (10 Years) | 1.00 (0.99, 1.01)
Age (10 Years) | 1.00 (0.95, 1.04)
Sex: |
Mixed Pair vs. Both Female | 0.99 (0.88, 1.11)
Both Male vs. Both Female | 0.99 (0.78, 1.24)
Residence Location: |
Mixed Pair vs. Both Rural | 0.98 (0.88, 1.08)
Both Urban vs. Both Rural | 0.97 (0.78, 1.19)
Working Status: |
Mixed Pair vs. Both Unemployed | 0.98 (0.84, 1.14)
Both Employed vs. Both Unemployed | 1.00 (0.72, 1.35)
Education: |
Mixed Pair vs. Both $\geq$ Secondary | 1.05 (0.95, 1.17)
Both $<$ Secondary vs. Both $\geq$ Secondary | 1.12 (0.90, 1.37)

###  5.2 Transmission probabilities

Results from the transmission probability analyses are shown in Figure 3 and
Table 4. Similar to the patristic distance analyses, WAIC also favors Spatial
over the competing approaches in this setting. Once again, Fixed may be
producing CIs that are too narrow while the Non-spatial point estimates may
have inflated MAEs. Graphical results from Spatial show several more
significant associations than observed with the competing approaches. Recall
that in the transmission probabilities framework, Non-spatial specifies that
$\boldsymbol{\eta}\left(\textbf{s}^{*}_{i}\right)=\textbf{0}_{4}$ for all $i$,
meaning that it accounts for neither spatial correlation nor cross-covariance
between the parameters. This could explain the relatively large differences
observed between results from the two methods.

Figure 3: Results from the transmission probability analyses in the Republic
of Moldova. Posterior means and 95% equal-tailed quantile-based credible
intervals are shown for the binary component (top row) and continuous
component (bottom row) raw regression parameters from each of the competing
methods. Red lines indicate that the interval excludes zero. WAIC values are
provided along with the expected number of parameters/model complexity term in
parentheses (i.e., p${}_{\text{WAIC}}$), with smaller values of WAIC
preferred.

For the pair-specific covariates in the binary regression model, odds ratio
results in Table 4 from Spatial suggest that pairs of individuals in the same
village and with similar dates of diagnosis are more likely to have non-zero
transmission probabilities on average. Results from the continuous regression
model are similar, with the addition of a significant finding for the distance
between villages. Pairs of individuals whose villages are closer together
geographically are more likely to have larger non-zero transmission
probabilities.

Table 4: Results from the Spatial transmission probability analyses in the Republic of Moldova. Posterior means and 95% quantile-based equal-tailed credible intervals are shown for the exponentiated regression parameters (i.e., odds ratios). The displayed credible intervals for the bolded entries exclude one. Effect | Binary |  | Continuous
---|---|---|---
Distance Between Villages (50 km) | 0.99 (0.89, 1.10) |  | 0.93 (0.87, 1.00)
Same Village (Yes vs. No) | 5.41 (2.22, 11.11) |  | 10.84 (5.54, 19.12)
Date of Diagnosis Difference (1/2 Year) | 0.77 (0.68, 0.86) |  | 0.75 (0.69, 0.82)
Age Difference (10 Years) | 1.03 (0.93, 1.13) |  | 0.93 (0.86, 1.01)
Giver Role: |  |  |
Age (10 Years) | 1.28 (1.04, 1.57) |  | 0.94 (0.86, 1.03)
Sex: |  |  |
Male vs. Female | 0.82 (0.46, 1.36) |  | 1.15 (0.87, 1.51)
Residence Location: |  |  |
Urban vs. Rural | 0.35 (0.11, 0.82) |  | 1.23 (0.90, 1.65)
Working Status: |  |  |
Employed vs. Unemployed | 4.69 (1.23, 12.26) |  | 1.05 (0.71, 1.49)
Education: |  |  |
$<$ Secondary vs. $\geq$ Secondary | 3.01 (1.49, 5.61) |  | 1.08 (0.82, 1.39)
Receiver Role: |  |  |
Age (10 Years) | 1.28 (1.01, 1.62) |  | 1.00 (0.90, 1.10)
Sex: |  |  |
Male vs. Female | 0.65 (0.33, 1.19) |  | 1.20 (0.90, 1.59)
Residence Location: |  |  |
Urban vs. Rural | 0.25 (0.06, 0.66) |  | 1.27 (0.90, 1.72)
Working Status: |  |  |
Employed vs. Unemployed | 6.18 (1.32, 18.66) |  | 1.20 (0.80, 1.73)
Education: |  |  |
$<$ Secondary vs. $\geq$ Secondary | 3.70 (1.60, 7.66) |  | 1.15 (0.86, 1.48)

For the individual-specific covariates, overall the giver and receiver results
are very similar within both the binary and continuous regression models,
suggesting that the impact of the included factors do not differ based on what
role the individual is serving in the pair. For the binary component of the
model we see that older, employed, and less formally educated individuals
living in rural areas are more likely to be in pairs that have non-zero
transmission probabilities. There were no statistically significant
individual-level associations identified for the continuous component of the
model.

###  5.3 Random effect parameter analyses

In Section S4 of the Supplementary Material, we present two additional
analyses for the estimated random effect parameters, $\theta_{i}$ and
$\left(\theta_{zi}^{(g)},\theta_{zi}^{(r)},\theta_{wi}^{(g)},\theta_{wi}^{(r)}\right)^{\text{T}}$,
from the fitted patristic distances and transmission probabilities models,
respectively. First, we present several summaries of the estimated random
effect parameters/hyperparameters with respect to the observed dyadic genetic
distance outcomes, finding an intuitive connection between the random effect
parameters and outcomes as well as the clear presence of spatially structured
variability in the data. Results suggest that the frameworks have the ability
to identify key individuals associated with increased transmission activity.
Next, we use Bayesian kriging techniques (Banerjee et al., 2014) to predict
the random effect parameters at new/unobserved spatial locations for both
outcomes across the Republic of Moldova. Posterior predicted mean and standard
deviation maps are presented which allow us to determine areas of increased
residual transmission risk for both outcomes.

##  6 Discussion

In this work, we presented innovative hierarchical Bayesian spatial methods
for modeling two types of dyadic genetic relatedness data, patristic distances
and transmission probabilities. The models account for multiple sources of
correlation (i.e., dyadic and spatial) and important features of each outcome
(e.g., zero-inflation). Through simulation, we showed that these approaches
perform as well as Fixed (i.e., regression while ignoring correlation) and
Non-spatial (i.e., no spatial or cross-covariance) when applied to
uncorrelated data, and greatly outperform them under correlation in terms of
estimating and quantifying uncertainty in the regression parameters. Under any
type of correlation, Fixed produces CIs that are too narrow and as a result
should not be used for estimating associations between genetic relatedness
measures and other factors.

When applying the models to _Mycobacterium tuberculosis_ data from the
Republic of Moldova, we found significant associations between genetic
relatedness and spatial proximity, dates of diagnosis, and multiple
individual-level factors; many of which represent new insights in this
setting. Analysis of the random effect parameters identified individuals and
geographic areas that are generally associated with higher levels of
transmission, a potentially useful aspect of the model with respect to
infection control as heterogeneity in infectiousness among TB patients has
been investigated in previous studies (Ypma et al., 2013; Melsew et al.,
2019). Additionally, WAIC showed that the correlation between outcomes was at
least partially spatially structured, further motivating the development of
the new methodology.

The unique data in our study represent a major strength of the work, however,
a number of factors inherent to _Mycobacterium tuberculosis_ make the study of
transmission dynamics difficult. First, TB epidemics are relatively slow
moving, often occurring over years (Pai et al., 2016), and it is difficult to
observe these processes in a two-year study. Second, the time between
infection and disease onset varies greatly, from weeks to years (Borgdorff et
al., 2011). Third, not all TB cases are ‘culture positive’, meaning it is more
difficult to culture, and, by extension, sequence isolates for every known TB
case (Cruciani et al., 2004; Nguyen et al., 2019). Taken together, these
factors make it difficult to infer direct transmission events. Even within a
putative transmission cluster, the probabilities of direct transmission among
case pairs may be low. For this reason, it can be difficult to fit a model to
these types of data. We address this challenge by including a binary component
in our model specification (i.e., an indicator of a non-zero transmission
probability between case pairs). Finally, as is standard in most
epidemiological analyses, our data do not include information on individuals
who were exposed but not infected. Therefore, the results should be
interpreted conditional on both individuals being infected.

While we use patristic distance as an outcome representative of symmetric
dyadic genetic relatedness data in this work, the framework we have created
can be readily adapted to SNP distances by modifying the likelihood in (1) to
accommodate discrete count data. In our R package GenePair (available at:
https://github.com/warrenjl/GenePair), we have developed software for fitting
the negative binomial and binary (i.e., genetically clustered or not) versions
of the patristic distances model, along with the two versions described in
Sections 3.1 and 3.2.

Future work in this area could focus on alternative techniques for accounting
for the correlation caused by analyzing dyadic outcomes while also
incorporating spatial correlation. Our current methods are based on the idea
of shared spatially correlated random effect parameters between pairs of data
including the same individual, which was shown to induce an intuitive
correlation structure between observations in Section 3.3. Additionally, these
parameters were shown to be useful in identifying key individuals/areas that
drive transmission in the study. However, there are undoubtedly other
approaches for achieving these same goals.

A meta-regression approach could also be used in future work for genetic
relatedness outcomes that are estimated in preceding analyses and include
measures of uncertainty (something not applicable in our data). The introduced
random effects frameworks could still be used for defining the true but
unobserved outcomes in this model, with a first stage that treats the observed
outcome as an estimate of the true value. A joint framework that combines the
model which estimates the genetic relatedness outcomes with our models for
characterizing variability and exploring associations in the outcomes could
also represent an important extension. However, the computational burden of
this sort of approach would likely be great and could possibly require the
addition of more efficient posterior sampling techniques.

Overall, we recommend the use of the newly developed methods when the goal of
a study is to estimate associations between spatially-referenced genetic
relatedness and other variables of interest. Even under the most simplistic,
and likely unrealistic assumption of no correlation between dyadic responses,
the new methodology performs well. When correlation is present, competing
approaches can yield potentially overly optimistic insights about the
significance of the associations and/or poorly estimate the associations of
interest and should not be used for decision making.

## References

  * Austin et al. (2013) Austin, A., C. Linkletter, and Z. Wu (2013).  Covariate-defined latent space random effects model.  Social Networks 35(3), 338–346.
  * Balaban et al. (2019) Balaban, M., N. Moshiri, U. Mai, X. Jia, and S. Mirarab (2019, 08).  Treecluster: Clustering biological sequences using phylogenetic trees.  PLoS ONE 14(8), 1–20.
  * Banerjee et al. (2014) Banerjee, S., B. P. Carlin, and A. E. Gelfand (2014).  Hierarchical Modeling and Analysis for Spatial Data.  CRC press.
  * Beck et al. (2006) Beck, N., K. S. Gleditsch, and K. Beardsley (2006).  Space is more than geography: Using spatial econometrics in the study of political economy.  International Studies Quarterly 50(1), 27–44.
  * Borgdorff et al. (2011) Borgdorff, M. W., M. Sebek, R. B. Geskus, K. Kremer, N. Kalisvaart, and D. van Soolingen (2011).  The incubation period distribution of tuberculosis estimated with a molecular epidemiological approach.  International Journal of Epidemiology 40(4), 964–970.
  * Campbell et al. (2019) Campbell, F., A. Cori, N. Ferguson, and T. Jombart (2019, 03).  Bayesian inference of transmission chains using timing of symptoms, pathogen genomes and contact data.  PLoS Computational Biology 15, 1–20.
  * Ciminelli et al. (2019) Ciminelli, J. T., T. Love, and T. T. Wu (2019).  Social network spatial model.  Spatial Statistics 29, 129–144.
  * Cruciani et al. (2004) Cruciani, M., C. Scarparo, M. Malena, O. Bosco, G. Serpelloni, and C. Mengoli (2004).  Meta-analysis of BACTEC MGIT 960 and BACTEC 460 TB, with or without solid media, for detection of mycobacteria.  Journal of Clinical Microbiology 42(5), 2321–2325.
  * Delsuc et al. (2005) Delsuc, F., H. Brinkmann, and H. Philippe (2005).  Phylogenomics and the reconstruction of the tree of life.  Nature Review Genetics 6, 361–375.
  * Didelot et al. (2017) Didelot, X., C. Fraser, J. Gardy, and C. Colijn (2017).  Genomic infectious disease epidemiology in partially sampled and ongoing outbreaks.  Molecular Biology and Evolution 34(4), 997–1007.
  * Gelman et al. (2013) Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin (2013).  Bayesian Data Analysis.  CRC press.
  * Geweke (1991) Geweke, J. (1991).  Evaluating the Accuracy of Sampling-based Approaches to the Calculation of Posterior Moments, Volume 196.  Federal Reserve Bank of Minneapolis, Research Department Minneapolis, MN, USA.
  * Gill and Swartz (2001) Gill, P. S. and T. B. Swartz (2001).  Statistical analyses for round robin interaction data.  Canadian Journal of Statistics 29(2), 321–331.
  * Hoff (2003) Hoff, P. D. (2003).  Dynamic Social Network Modeling and Analysis: Workshop Summary and Papers, Chapter Random effect models for network data, pp. 303–312.  Washington DC: National Academies Press.
  * Hoff (2005) Hoff, P. D. (2005).  Bilinear mixed-effects models for dyadic data.  Journal of the american Statistical association 100(469), 286–295.
  * Hoff (2021) Hoff, P. D. (2021).  Additive and multiplicative effects network models.  Statistical Science 36(1), 34–50.
  * Hoff et al. (2002) Hoff, P. D., A. E. Raftery, and M. S. Handcock (2002).  Latent space approaches to social network analysis.  Journal of the american Statistical Association 97(460), 1090–1098.
  * Kenny et al. (2020) Kenny, D. A., D. A. Kashy, and W. L. Cook (2020).  Dyadic Data Analysis.  Guilford Publications.
  * Klinkenberg et al. (2017) Klinkenberg, D., J. A. Backer, X. Didelot, C. Colijn, and J. Wallinga (2017, 05).  Simultaneous inference of phylogenetic and transmission trees in infectious disease outbreaks.  PLoS Computational Biology 13, 1–32.
  * Loman and Pallen (2015) Loman, N. and M. Pallen (2015).  Twenty years of bacterial genome sequencing.  Nature Review Microbioly 13, 787–794.
  * Melsew et al. (2019) Melsew, Y. A., M. Gambhir, A. C. Cheng, E. S. McBryde, J. T. Denholm, E. L. Tay, and J. M. Trauer (2019).  The role of super-spreading events in Mycobacterium tuberculosis transmission: evidence from contact tracing.  BMC Infectious Diseases 19(1), 1–9.
  * Neumayer and Plümper (2010) Neumayer, E. and T. Plümper (2010).  Spatial effects in dyadic data.  International Organization 64(1), 145–166.
  * Nguyen et al. (2019) Nguyen, M.-V. H., N. S. Levy, S. D. Ahuja, L. Trieu, D. C. Proops, and J. M. Achkar (2019, 02).  Factors associated with sputum culture-negative vs culture-positive diagnosis of pulmonary tuberculosis.  JAMA Network Open 2(2), e187617–e187617.
  * Pai et al. (2016) Pai, M., M. A. Behr, D. Dowdy, K. Dheda, M. Divangahi, C. C. Boehme, A. Ginsberg, S. Swaminathan, M. Spigelman, H. Getahun, D. Menzies, and M. Raviglione (2016).  Tuberculosis.  Nature Reviews Disease Primers 2(1), 16076.
  * Polonsky et al. (2019) Polonsky, J. A., A. Baidjoe, Z. N. Kamvar, A. Cori, K. Durski, W. J. Edmunds, R. M. Eggo, S. Funk, L. Kaiser, P. Keating, et al. (2019).  Outbreak analytics: A developing data science for informing the response to emerging pathogens.  Philosophical Transactions of the Royal Society B 374(1776), 20180276.
  * Poon (2016) Poon, A. F. Y. (2016).  Impacts and shortcomings of genetic clustering methods for infectious disease outbreaks.  Virus Evolution 2(2).  vew031.
  * Simpson and Laurienti (2015) Simpson, S. L. and P. J. Laurienti (2015).  A two-part mixed-effects modeling framework for analyzing whole-brain network data.  NeuroImage 113, 310–319.
  * Stamatakis (2014) Stamatakis, A. (2014, 01).  RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies.  Bioinformatics 30(9), 1312–1313.
  * Walker et al. (2014) Walker, T. M., M. K. Lalor, A. Broda, L. S. Ortega, M. Morgan, L. Parker, S. Churchill, K. Bennett, T. Golubchik, A. P. Giess, C. Del Ojo Elias, K. J. Jeffery, I. C. J. W. Bowler, I. F. Laurenson, A. Barrett, F. Drobniewski, N. D. McCarthy, L. F. Anderson, I. Abubakar, H. L. Thomas, P. Monk, E. G. Smith, A. S. Walker, D. W. Crook, T. E. A. Peto, and C. P. Conlon (2014).  Assessment of Mycobacterium tuberculosis transmission in Oxfordshire, UK, 2007-12, with whole pathogen genome sequences: An observational study.  The Lancet Respiratory Medicine 2(4), 285–292.
  * Warner et al. (1979) Warner, R. M., D. A. Kenny, and M. Stoto (1979).  A new round robin analysis of variance for social interaction data.  Journal of Personality and Social Psychology 37(10), 1742\.
  * Watanabe (2010) Watanabe, S. (2010).  Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  Journal of Machine Learning Research 11, 3571–3594.
  * Wong (1982) Wong, G. Y. (1982).  Round robin analysis of variance via maximum likelihood.  Journal of the American Statistical Association 77(380), 714–724.
  * Yang et al. (2022) Yang, C., B. Sobkowiak, V. Naidu, A. Codreanu, N. Ciobanu, K. S. Gunasekera, M. H. Chitwood, S. Alexandru, S. Bivol, M. Russi, J. Havumaki, P. Cudahy, H. Fosburgh, C. J. Allender, H. Centner, D. M. Engelthaler, N. A. Menzies, J. L. Warren, V. Crudu, C. Colijn, and T. Cohen (2022).  Phylogeography and transmission of M. tuberculosis in Moldova: A prospective genomic analysis.  PLoS Medicine 19(2), e1003933.
  * Ypma et al. (2013) Ypma, R. J., H. K. Altes, D. Van Soolingen, J. Wallinga, and W. M. Van Ballegooijen (2013).  A sign of superspreading in tuberculosis: Highly skewed distribution of genotypic cluster sizes.  Epidemiology 24(3), 395–400.


# Cascaded Self-Representation for Outlier Detection

###### Abstract

Many methods based on sparse and low-rank representation been developed along
with guarantees of correct outlier detection. Self-representation states that
a point in a subspace can always be expressed as a linear combination of other
points in the subspace. A suitable Markov Chain can be defined on the self-
representation and it allows us to recognize the difference between inliers
and outliers.However, the reconstruction error of self-representation that is
still informative to detect outlier detection, is neglected. Inspired by the
gradient boosting, in this paper, we propose a new outlier detection framework
that combines a series of weak ”outlier detectors” into a single strong one in
an iterative fashion by constructing multi-pass self-representation. At each
stage, we construct a self-representation based on elastic-net and define a
suitable Markov Chain on it to detect outliers. The residual of the self-
representation is used for the next stage to learn the next weaker outlier
detector. Such a stage will repeat many times. And the final decision of
outliers is generated by the previous all results. Experimental results on
image and speaker datasets demonstrate its superiority with respect to state-
of-the-art sparse and low-rank outlier detection methods.

Index Terms—  Outlier Detection, Self-Representation, Sparse Coding

##  1 Introduction

Outlier detection, also called anomaly detection, identifies unusual data
patterns that are different from the majority of data. These unexpected
patterns are always called anomalies or outliers, the ability to detect
anomalies has significant and critical help, and anomalies often provides
useful information in various application domains, such as intrusion
detection, fraud detection, fault detection, suspicious transaction detection
and abnormal moving activity detection, etc. Anomaly detection can be roughly
categorized into three ways: statistical, algebraic and self-representation
based.

RANdom SAmple Consensus (RANSAC) [1], as a statistical method for outlier
detection, employs a sampling strategy to find a subspace that fits as many
data as possible iteratively. In this processing, data points are removed from
the dataset to find a new subspace until a given threshold on the percents of
inliers is reached. Although theoretically, it is able to find correct
subspace in noisy data even in the presence of outliers, the computational
cost is still a challenge. Its variants like [2] and [3] still face the
similar problem.

Algebraic methods have been to robustly learn the subspaces by penalizing the
sum of an unsquared distance of points to the closest subspace [4, 5] compared
with Principal Component Analysis (PCA) minimizing the sum of square error.
Thus they are robust to outliers because of reducing the contributions from
large residuals arising from outliers. However, the optimization problem is
usually nonconvex and a good initialization is extremely important for finding
the optimal solution. Recently the PCA problem with robustness have been
solved to corrupted entries [6], which has led to many recent methods for PCA
with robustness to outliers. A prominent advantage of convex optimization
techniques is that they are guaranteed to correctly identify outliers under
certain conditions. Nonetheless, these methods typically model a unique inlier
subspace, e.g., by a low-rank constraint in Outlier Pursuit [7], and therefore
cannot deal with multiple inlier subspaces since the union of multiple
subspaces could be high-dimensional.

An alternative to algebraic approaches are methods that exploit the self-
expressiveness of subspaces. In most of these approaches, the goal becomes
finding a clustering that renders a suitably defined affinity matrix,
typically obtained by solving a sparsification [8] or rank minimization
problem [9]. These methods work well in the presence of noise and offer
theoretical recovery guarantees. On the other hand, handling outliers requires
augmenting the objective function with additional regularization terms, whose
parameters must be usually hand-tuned to obtain good performance and typically
recovery guarantees are lost. Based on self-representation, [10] design a
random walk process and identify outliers as those whose probabilities tend to
zero, which is parameter-free. However, the residual of self-representation
which is neglected is still informative and contains clues to make outlier
detection better.

In this paper, we introduce a novel cascade architecture that aims to extend
the random walk based outlier detection onto hierarchical self-representations
in a principled way. To achieve the information from the reconstruction error
in a linear combination, we present a multi-stages cascade framework for self-
representation based on sparse coding with Elastic Net [11]. The residuals at
a layer is computed by the difference between the original input and the
aggregated reconstructions of the previous layers. And then the residuals at a
layer can be represented with self-representation just like the original
inputs did. Multi self-representation can be used to detect outliers with a
random work based method, the final decision of outliers are fused by a linear
combination of different results from different layers.

##  2 Related work

Motivated by the observation that outliers do not have sparse representations,
[12] declare a point as an outlier if its $\ell 1$ norm of representation is
above a threshold. However, this $\ell 1$-thresholding strategy is not robust
to outliers that is not good enough to discriminate outliers since their
representation vectors may have many small small $\ell 1$-norms. Low-Rank
Representation (LRR) [9] employs sparsity assumption (i.e. $\ell_{2,1}$-norm)
on the reconstruction error of low rank self-representation to learn a robust
model, and the column-sparse matrix rather than the low-rank matrix can be
used to indicate the outliers.

In [13], a graph-based outlier detection framework is proposed. The main idea
is to represent the underlying dataset as a weighted undirected graph and then
use the random work on the graph to find outliers. The random walk model is
designed to find nodes that are most “central” to the graph and thus the
outliers are far from the ”central” can be efficiently detected. [10] extends
graph construction from kernel based similarity to sparse representation based
self-representation that inliers express themselves with only other inliers
when they lie in a union of low dimensional subspaces. To avoid the non-
convergence because of directed graph, researcher choose to calculate the
average of the first $T$ t-step probability distributions.

Fig. 1: The Framework of Cascaded Self-Representation for Outlier Detection

##  3 Method

In this section, we present cascaded self-representation based outlier
detection method shown in Fig. 1. We first describe the data self-
representation and a random walk algorithm on the graph based on the
representation to identify the sets of inliers and outliers. Then, we
introduce a cascaded architecture that calculates self-representation from the
resulting residuals of the previous layers. Finally we ensemble the all the
results into one set of scores to identify the inliers and outliers.

###  3.1 Elastic Net based Self-Representation

Given unlabeled data points $\\{X\\}_{i=1,...,N}$ drawn from multiple linear
subspaces $\\{S\\}_{i=1,...,K}$, one can express a point in a subspace as a
linear combination of other points in the same subspace, and the data contains
inliers and outliers. We can get the self-representation coefficient matrix by
the optimization problem:

$\quad\min\limits_{C}\frac{\gamma}{2}\|X-XC\|^{2}_{2}+\lambda\|C\|_{1}+\frac{1-\lambda}{2}\|C\|_{2}^{2}\quad
s.t.\;C_{ii}=0\\\ $ (1)

where $C$ is the self-representation coefficient matrix, $\gamma>0$ and
$\lambda\in[0,1)$ , the optional diagonal constraint on $C$ prevents trivial
solutions for sparsity-inducing norms, such as the $\ell 1$ norm. Compared
with the $\ell 1$ norm constraint, A mixture of $\ell 1$ and $\ell 2$ norms is
able to balance the subspace preserving and connectedness properties. However,
the residual term $\|X-XC\|_{2}^{2}$, also called reconstruction error, is
ignored by the later random walk based outlier detection.

###  3.2 Graph-based Outlier Detection

With a Self-Representation $C$, we can get the directed weighted graph $G$:
The vertices of $G$ correspond to the data points $X$, and the edges are given
by the (weighted) adjacency matrix $A=|C|^{T}$. In this graph $G$, inliers
only have connections with the inliers while the outliers have edges with both
inliers and outliers, which we can use random walk to detect the outliers. The
transition probability from $x_{i}$ to $x_{j}$ at the next time is given by
$p_{ij}=a_{ij}/d_{i}$ with $d_{i}=\sum_{j}a_{ij}$, where $p_{ij}\in P$ and
$a_{ij}\in A$ . By this definition, if the starting point of a random walk is
an inlier then it will never escape the set of inliers as there is no edge
going from any inlier to any outlier. In contrast, a random walk starting from
an outlier will likely end up in an inlier state since once it enters any
inlier it will never return to an outlier state. Thus, by using different data
points to initialize random walks, outliers can be identified by observing the
final probability distribution of the state of the random walks. A uniform
distribution $\pi^{(0)}$ = [$\frac{1}{N}$,…,$\frac{1}{N}$] is defined as the
initial probability distribution, we define the formula:

$\hat{\pi}^{(T)}=\frac{1}{T}\sum_{t=1}^{T}\pi^{(0)}P^{t}$ (2)

where the average of the t-step can ensure the convergence . It is expected
that eventually the inliers have high probabilities states and outliers have
low probabilities, so we can use the threshold value $\epsilon$ to identify
the ouliers if $\hat{\pi}^{(T)}\leq\epsilon$.

###  3.3 Outlier Detection by Cascaded Self-Representation

As mentioned above, the method of self-representation based outlier detection
often formulates the problem at hand using a linear model with regularization
by $\ell 1$ and $\ell 2$-norms, but the reconstruction error is ignored that
hinders exploiting self-representation in their full potential. In comparison,
our approach exploits a recursive way where we encode self-representation with
the resulting residuals of the stage in previous stages. In a single stage, we
represent the residual as a linear combination of other points, where we keep
the same as single stage. Let $\hat{X}^{n}$ denote the estimated $n$-th stage
and $X$ denote unlabeled data points, then the overall process can be
described as:

$X=\hat{X}^{0}+\hat{X}^{1}+...\hat{X}^{n}\quad
s.t.\;\hat{X}^{i}=(X-\sum_{j=0}^{i-1}\hat{X}^{j})C^{i},\;C^{i}_{jj}=0$ (3)

where $\hat{X}^{0}$ is a matrix with zero entries which is used to easily
formulate our motivation with the Equ. 3.

The flow diagram of our cascade framework is showed as Fig.1 and the algorithm
is Alg. 1. Given the data $X$, we first get the self-representation $C^{1}$
and the transition probability, then run Equ. 2 on it to get the outlier score
$S^{1}$, for the $\hat{X}^{1}$, it can be solve by $\hat{X}^{1}=XC^{1}$. And
then we can use the $X-\hat{X}^{1}$ do again the self-representation and
random-walk to get the outlier score $S^{2}$, we know that the outlier score
$S^{1}$ get from the random-walk’s initial probability distribution is
$\pi^{(0)}$, however, the outlier score $S^{2}$ solved by the same process’s
initial probability distribution is $S^{1}$. The rest can be done in the same
manner to get the $\hat{C}^{2}$, $\hat{C}^{3}$,…,$\hat{C}^{n}$, and the
residuals $X-\sum_{j=0}^{i-1}\hat{X}^{j}$, and the $S^{3}$, $S^{4}$,…,$S^{n}$.
After all is done, do the score fusion of the $S^{1}$,…,$S^{n}$ to get the
final Outlier Score $S$, and identify the Outliers if $S\leq\epsilon$.

0: Unlabeled data points $X$, the threshold value $\epsilon$

0: Outlier score $S$ of each points and Outliers by $\epsilon$

1: Initial probability distribution $\pi^{(0)}=[\frac{1}{N},...,\frac{1}{N}]$

2: for $i=1$ to $n$ do

3: if $i=1$ then

4: $C^{1}$ $\leftarrow$ solved by Equ. 1

5: Outlier Score $S^{1}$ $\leftarrow$ solved by Equ. 2 with $\pi^{(0)}$

6: $\hat{X}^{1}$ $\leftarrow$ solved by Equ. 3

7: end if

8: Update probability distribution = Outlier Score $S^{(i-1)}$

9: $C^{i}$ $\leftarrow$ solved by Equ. 1 with residuals
$X-\sum_{j=0}^{i-1}\hat{X}^{j}$

10: Outlier Score $S^{i}$ $\leftarrow$ solved by Equ. 2 with $S^{(i-1)}$

11: $\hat{X}^{i}$ $\leftarrow$ solved by Equ. 3

12: end for

13: Score fusion result $S$ $\leftarrow$ from $S^{1},S^{2},..,S^{n}$

14: Identify the Outliers if $S\leq\epsilon$

15: return Outlier score $S$ and Outliers

Algorithm 1 Outlier Detection by Cascaded Self-Representation

##  4 Experiments

In this section, we use three different image datasets and one speaker dataset
to evaluate the proposed method for outlier detection. In particular, we
analyzed the performance results in details and compared with state-of-the-art
techniques.

###  4.1 Experimental setup

We construct our method in matlab, and evaluated it on four publicly available
databases, i.e., the Extended Yale B [14] , the Caltech-256 [15], the Coil-100
object image dataset [16] and the TIMIT Small dataset [17]. For our method we
set the number of iterations T to be 1000, and the number of stages is 3 since
larger stages do not contribute better performance significantly.

We compare with 3 other representative methods that are designed for detecting
outliers in one or multiple subspaces: LRR [9] , $\ell$-thresholding [12] and
R-graph [10]. All other methods are implemented according to the description
in their respective papers. We call our method Outlier Detection by Cascaded
Self-Representation (ODCSR) in this section.

For each outlier detection method, we use two metric to evaluate its
performance. One is the area under the curve (AUC) as a metric of performance
in terms of the ROC, the other is the F1-score, which is the harmonic mean of
precision and recall. Notice that a numerical value for each data point that
indicates its “outlierness” and a threshold value for determining inliers and
outliers are required.

###  4.2 Extended Yale B Dataset

The Extended Yale B dataset is a popular benchmark for subspace clustering. It
contains frontal face images of 38 individuals each under 64 different
illumination conditions. Following the experimental setup of [8], we down-
sampled the original face images from $192\times 168$ to $42\times 42$ pixels,
which makes it computationally feasible for the baselines.

We designed a series of experiments to test the ability of the methods to deal
with multiple inlier groups. For example, we randomly choose 1 or 3 individual
from the 38 subjects and use all 64 images of the subjects as the inliers.
Then the images from the remaining 37 or 35 subjects as outliers with at most
one image from each subject. The overall data set has 35% or 15% outliers. The
results of this experiment are reported in Table 1. The last column have the
results of our methods, note that in each row the best result is typeset in
bold. We can see that our cascade subspace clustering method is the best
method.

Table 1: Results on the Extended Yale B dataset | LRR |  $\ell 1$-thresholding | R-graph | ODCSR(ours)
---|---|---|---|---
Inliers:from one subject, Outliers:35% from other subjects
AUC | 0.857 | 0.844 | 0.986 | 0.990
F1 | 0.797 | 0.763 | 0.951 | 0.956
Inliers:from three subjects, Outliers:15% from other subjects
AUC | 0.807 | 0.848 | 0.985 | 0.986
F1 | 0.509 | 0.545 | 0.878 | 0.886

###  4.3 Caltech-256 Dataset

For Caltech-256, which contains 256 object categories with a total of 30,607
images, and each category has at least 80 images. In our experiments, images
in Caltech-256 are represented by a 4,096-dimensional feature vector extracted
from the last fully connected layer of the 16-layer VGG network [18]. Like the
above, images in $n\in\\{1,3,5\\}$ randomly chosen categories and the first
150 images in it are used as inliers, and we randomly pick a certain number of
images from the “clutter” category as outliers such that there are 50%
outliers in each experiment.

The results of this experiment are reported in Table 2. The last column have
the results of our methods, note that in each row the best result is typeset
in bold, and we are the best method.

Table 2: Results on the Caltech-256 dataset | LRR |  $\ell 1$-thresholding | R-graph | ODCSR(ours)
---|---|---|---|---
Inliers:from one category, Outliers:50% from 257-clutter
AUC | 0.907 | 0.772 | 0.948 | 0.983
F1 | 0.893 | 0.772 | 0.914 | 0.946
Inliers:from three categories, Outliers:50% from 257-clutter
AUC | 0.479 | 0.810 | 0.929 | 0.984
F1 | 0.671 | 0.782 | 0.880 | 0.947
Inliers:from five categories, Outliers:50% from 257-clutter
AUC | 0.337 | 0.774 | 0.913 | 0.984
F1 | 0.667 | 0.762 | 0.858 | 0.952

###  4.4 Coil-100 Dataset

The Coil-100 dataset contains 7,200 images of 100 different objects. Each
object has 72 images taken at pose intervals of 5 degrees, with the images
being of size $32\times 32$. For Coil-100, we randomly pick $n\in\\{1,4,7\\}$
categories as inliers and pick at most one image from each of the remaining
categories as outliers. The results of this experiment are reported in Table
3. And we have the best performance.

Table 3: Results on the Coil-100 dataset | LRR |  $\ell 1$-thresholding | R-graph | ODCSR(ours)
---|---|---|---|---
Inliers:from one subject, Outliers:50% from other subjects
AUC | 0.847 | 0.991 | 0.997 | 0.999
F1 | 0.872 | 0.978 | 0.900 | 0.995
Inliers:from four subjects, Outliers:25% from other subjects
AUC | 0.687 | 0.992 | 0.996 | 0.998
F1 | 0.541 | 0.941 | 0.970 | 0.981
Inliers:from seven subjects, Outliers:15% from other subjects
AUC | 0.628 | 0.991 | 0.996 | 0.997
F1 | 0.366 | 0.897 | 0.955 | 0.963

###  4.5 TIMIT Small Dataset

The TIMIT dataset is composed of 6,300 phrases (10 phrases per speaker),
spoken by 438 males (70%) and 192 females (30%). In our experiment, we used
the same 40 speakers dataset as reported by these earlier attempts (here
called TIMIT Small). we randomly pick $n\in\\{3,5,7\\}$ speakers as inliers
and pick at most one phrase from each of the remaining speakers as outliers.
The results of this experiment are reported in Table 4. And we can see that
our method is still the best.

Table 4: Results on the TIMIT Small dataset | LRR |  $\ell 1$-thresholding | R-graph | ODCSR(ours)
---|---|---|---|---
Inliers:from three speakers, Outliers:15% from other speakers
AUC | 0.867 | 0.873 | 0.959 | 0.972
F1 | 0.698 | 0.660 | 0.837 | 0.867
Inliers:from five speakers, Outliers:15% from other speakers
AUC | 0.847 | 0.839 | 0.944 | 0.951
F1 | 0.630 | 0.598 | 0.774 | 0.814
Inliers:from seven speakers, Outliers:15% from other speakers
AUC | 0.813 | 0.826 | 0.936 | 0.944
F1 | 0.536 | 0.557 | 0.749 | 0.759

##  5 Conclusion

In this paper, we have proposed a general framework for outlier detection in a
manner of cascade subspace clustering. Our architecture consists of multi-
stages, each one is random walk outlier detection based on self-representation
with the residual of last layers. Eventually, each result from different paths
is fused to make the final decision of outliers. Our experiments have
demonstrated that our cascade method provides significant improvement over the
state-of-the-art outlier detection solutions in terms of AUC and F1 on several
datasets.

## References

  * [1] Martin A Fischler and Robert C Bolles,  “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,”  Communications of the ACM, vol. 24, no. 6, pp. 381–395, 1981.
  * [2] Phil Torr and Andrew Zisserman,  “Robust computation and parametrization of multiple view relations,”  in Computer Vision, 1998. Sixth International Conference on. IEEE, 1998, pp. 727–732.
  * [3] Philip H. S. Torr,  “Bayesian model estimation and selection for epipolar geometry and generic manifold fitting,”  International Journal of Computer Vision, vol. 50, no. 1, pp. 35–61, 2002.
  * [4] Gilad Lerman, Teng Zhang, et al.,  “Robust recovery of multiple subspaces by geometric lp minimization,”  The Annals of Statistics, vol. 39, no. 5, pp. 2686–2715, 2011.
  * [5] Teng Zhang, Arthur Szlam, Yi Wang, and Gilad Lerman,  “Hybrid linear modeling via local best-fit flats,”  International journal of computer vision, vol. 100, no. 3, pp. 217–240, 2012.
  * [6] John Wright, Arvind Ganesh, Shankar Rao, Yigang Peng, and Yi Ma,  “Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization,”  in Advances in neural information processing systems, 2009, pp. 2080–2088.
  * [7] Huan Xu, Constantine Caramanis, and Sujay Sanghavi,  “Robust pca via outlier pursuit,”  in Advances in Neural Information Processing Systems, 2010, pp. 2496–2504.
  * [8] Ehsan Elhamifar and René Vidal,  “Sparse subspace clustering,”  in Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on. IEEE, 2009, pp. 2790–2797.
  * [9] Guangcan Liu, Zhouchen Lin, and Yong Yu,  “Robust subspace segmentation by low-rank representation,”  in Proceedings of the 27th international conference on machine learning (ICML-10), 2010, pp. 663–670.
  * [10] Chong You, Daniel P Robinson, and René Vidal,  “Provable selfrepresentation based outlier detection in a union of subspaces,”  in Proceedings of the 2017 IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 1–10.
  * [11] Hui Zou and Trevor Hastie,  “Regularization and variable selection via the elastic net,”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 67, no. 2, pp. 301–320, 2005.
  * [12] Mahdi Soltanolkotabi, Emmanuel J Candes, et al.,  “A geometric analysis of subspace clustering with outliers,”  The Annals of Statistics, vol. 40, no. 4, pp. 2195–2238, 2012.
  * [13] HDK Moonesinghe and Pang-Ning Tan,  “Outrank: a graph-based outlier detection framework using random walk,”  International Journal on Artificial Intelligence Tools, vol. 17, no. 01, pp. 19–36, 2008.
  * [14] Athinodoros S. Georghiades, Peter N. Belhumeur, and David J. Kriegman,  “From few to many: Illumination cone models for face recognition under variable lighting and pose,”  IEEE Trans. Pattern Anal. Mach. Intell., vol. 23, no. 6, pp. 643–660, June 2001.
  * [15] Gregory Griffin, Alex Holub, and Pietro Perona,  “Caltech-256 object category dataset,”  2007\.
  * [16] Sameer A Nene, Shree K Nayar, Hiroshi Murase, et al.,  “Columbia object image library (coil-20),”  1996\.
  * [17] Feliks Hibraj, Sebastiano Vascon, Thilo Stadelmann, and Marcello Pelillo,  “Speaker clustering using dominant sets,”  arXiv preprint arXiv:1805.08641, 2018.
  * [18] Karen Simonyan and Andrew Zisserman,  “Very deep convolutional networks for large-scale image recognition,”  arXiv preprint arXiv:1409.1556, 2014.