summaries: Authors demonstrated that one can encode the data likelihood function of an HMM using a specialized RNN architecture. Unlike previous work where neurons from different layers were multiplied together, the new encoding strictly followed the classical architecture restrictions of a neural network , i.e. each layer was a weighted sum of the previous layer. Empirically, author showed that the parameter learned by applying gradient descent on the likelihood function is similar to the one that is obtained using the EM algorithm. In addition, authors demonstrated that such formulation enables an application in studying Alzheimer's disease Symptom Progression.

summaries: This paper introduces a novel architecture of recurrent neural network that mimics the working of a standard HMM. In particular, the proposed HMRNN learns model-parameters, which are statistically similar solutions to those of a standard HMM, within a general neural network learning framework. While it is shown the proposed network similarly to the HMM, there are many issues that should be considered critically.

summaries: The authors describe an RNN architecture, the HMRNN, which models the log-likelihood in an HMM.  The authors provide theoretical results showing that the HMRNN objective function does, indeed, correspond to the log-likelihood of an HMM.  Synthetic results are presented which compare the learned parameters between the HMRNN and an HMM trained using the Baum-Welch algorithm.  Finally, HMRNN training is augmented and compared to an HMM trained using the Baum-Welch algorithm for the task of Alzheimer's progression prediction.

abstractText: Hidden Markov models (HMMs) are commonly used for disease progression modeling when the true patient health state is not fully known. Since HMMs may have multiple local optima, performance can be improved by incorporating additional patient covariates to inform estimation. To allow for this, we formulate a special case of recurrent neural networks (RNNs), which we name hidden Markov recurrent neural networks (HMRNNs), and prove that each HMRNN has the same likelihood function as a corresponding discrete-observation HMM. The HMRNN can be combined with any other predictive neural networks that take patient covariate information as input. We show that HMRNN parameter estimates are numerically close to those obtained from via the Baum-Welch algorithm, validating their theoretical equivalence. We then demonstrate how the HMRNN can be combined with other neural networks to improve parameter estimation, using an Alzheimer’s disease dataset. The HMRNN’s solution improves disease forecasting performance and offers a novel clinical interpretation compared with a standard HMM.

authors:
T. Ayer
O. Alagoz
N.K. Stout

title: A POMDP approach to personalize mammography screening

venue: decisions. Operations Research,

authors:
M. Baucum
A. Khojandi
R. Vasudevan

title: Improving deep reinforcement learning with transitional variational autoencoders: A healthcare application

venue: Journal of Biomedical & Health Informatics (forthcoming),

authors:
L.E. Baum
T. Petrie

title: Statistical inference for probabilistic functions of finite state Markov chains

venue: The Annals of Mathematical Statistics,

title: Cerebrospinal fluid protein biomarkers for Alzheimer’s disease

venue: NeuroRx,

title: Alpha-nets: A recurrent ‘neural’ network architecture with a hidden Markov model interpretation

venue: Speech Communication,

authors:
A. Bussone
S. Stumpf
D. O’Sullivan

title: The role of exlpanations on trust and reliance in clinical decision support

authors:
L. Canuet
S. Pusil
M.E. Lopez
R. Bajo
J.A. Pineda-Pardo
P. Cuesta
G. Galvez
J.M. Gatzelu
D. Lourido
G. Garcia-Ribas
F. Maestu

title: Network disruption and cerebrospinal fluid amyloidbeta and phospho-tau levels in mild cognitive impairment

venue: Journal of Neuroscience,

authors:
A. Chopra
T.A. Cavalieri
D.J. Libon

title: Dementia screening tools for the primary care physician

venue: Clinical Geriatrics,

authors:
B. Estebanez
P. del Saz-Orozco
I. Rivas
E. Bauzano
V.F. Muñoz
I. Garcia-Morales

title: Maneuvers recognition in laparoscopic surgery: Artificial Neural Network and hidden Markov model approaches

authors:
M.F. Folstein
S.E. Folstein
P.R. McHugh

title: Mini-mental state: A practical method for grading the cognitive state of patients for the clinician

venue: Journal of Psychiatric Research,

authors:
F.A. Guerrero-Pena
P.D.M. Fernandez
T.I. Ren
M. Yui
E. Rothenberg
A. Cunha

title: Multiclass weighted loss for instance segmentation of cluttered cells

authors:
C. Guihenneuc-Jouyaux
S. Richardson
I.M. Longini Jr.

title: Modeling markers of disease progression by a hidden Markov process: Application to characterizing CD4 cell decline

title: Developing clinical decision support systems for sepsis prediction using temporal and non-temporal machine learning method

venue: Unpublished doctoral dissertation,

authors:
Oskar Hansson
Henrik Zetterberg
Eugeen Vanmechelen
Hugo Vanderstichele
Ulf Andreasson
Elisabet Londos
Anders Wallin
Lennart Minthon
Kaj Blennow

title: Evaluation of plasma aβ40 and aβ42 as predictors of conversion to alzheimer’s disease in patients with mild cognitive impairment

venue: Neurobiology of aging,

authors:
K. Hirao
T. Ohnishi
Y. Hirata
F. Yamashita
T. Mori
Y. Moriguchi
T. Iwamoto

title: The prediction of rapid conversion to Alzheimer’s disease in mild cognitive impairment using regional cerebral blood flow spect

authors:
M. Igl
L. Zintgraf
T.A. Le
F. Wood
S. Whiteson

title: Deep variational reinforcement learning for POMDPs

venue: arXiv preprint,

authors:
C.H. Jackson
L.D. Sharples
S.G. Thompson
S.W. Duffy
E. Couto

title: Multistate markov models for disease progression with classification error

venue: Journal of the Royal Statistical Society: Series D,

authors:
D. Jurafsky
J. Martin

title: Speech and language processing

title: A POMDP approach to age-dependent primary screening policies for cervical cancer in colombia

venue: Unpublished Master’s thesis,

authors:
B.Y. Lim
A.K. Dey

title: Toolkit to support intelligibility in context-aware applications, pp

authors:
Y.Y. Liu
S. Li
F. Li
L. Song
J.M. Rehg

title: Efficient learning of continuous-time hidden Markov models for disease progression

venue: Advances in Neural Information Processing Systems,

authors:
A. Martino
G. Guatteri
A.M. Paganoni

title: Multivariate hidden markov models for disease progression

venue: Statistical Analysis and Data Mining: The ASA Data Science Journal,

authors:
T. Monroe
M. Carter

title: Using the Folstein mini mental state exam (MMSE) to explore methodological issues in cognitive aging research

venue: European Journal of Ageing,

title: Clinical dementia rating: Current version and scoring

venue: rules. Neurology,

authors:
S. Nemati
M.M. Ghassemi
G.D. Clifford

title: Optimal medication dosing from suboptimal clinical examples: A deep reinforcement learning approach

venue: 38th Annual International Conference of IEEE Engineering in Medicine and Biology Society,

authors:
L.A. Rabin
N. Pare
A.J. Saykin
M.J. Brown
H.A. Wishart
L.A. Flashman
R.B. Santulli

title: Differential memory test sensitivity for diagnosing amnestic mild cognitive impairment and predicting conversion to alzheimer’s disease

venue: Aging, Neuropsychology, and Cognition,

authors:
I. Stanculescu
C.K. Williams
Y. Freer

title: Autoregressive hidden Markov models for the early detection of neonatal sepsis

venue: IEEE Journal of Biomedical and Health Informatics,

authors:
M. Stelmach
M. Chlebus

title: Novel multilayer stacking framework with weighted ensemble approach for multiclass credit scoring problem application

authors:
R. Sukkar
E. Katz
Y. Zhang
D. Raunig
B.T. Wyman

title: Disease progression modeling using hidden Markov models

authors:
T. Wessels
C.W. Omlin

title: Refining hidden Markov models with recurrent neural networks

venue: In Proceedings of the IEEE-INNS-ENNS International Joint Conference on Neural Networks,

authors:
S. Yildirim
S.S. Singh
T. Dean
A. Jasra

title: Parameter estimation in hidden markov models with intractable likelihoods using sequential monte carlo

venue: Journal of Computational and Graphical Statistics,

authors:
Z. Zhou
Y. Wang
H. Mamani
D.G. Coffey

title: How do tumor cytogenetics inform cancer treatments? dynamic risk stratification and precision medicine using multi-armed bandits

heading: 1 INTRODUCTION

text:
Hidden Markov models (HMMs; Baum & Petrie, 1966) are commonly used for modeling disease progression, because they allow researchers to conceptualize complex (and noisy) clinical measurements as originating from a smaller set of latent health states. Each latent health state is characterized by an emission distribution that specifies the probabilities of each measurement/observation given that state. This allows HMMs to explicitly account for uncertainty or measurement error, since the system’s true state is not fully observable. Because of their intuitive parameter interpretations and flexibility, HMMs have been used to model biomarker changes in HIV patients (Guihenneuc-Jouyaux et al., 2000), Alzheimer’s disease progression (Liu et al., 2015), breast cancer screening decisions (Ayer et al., 2012), and patient response to blood anticoagulants (Nemati et al., 2016).
Researchers may wish to integrate HMMs with other disease progression models and/or data sources. For instance, researchers in Igl et al. (2018) jointly trained parameters for an HMM and a reinforcement learning policy to maximize patient returns. Other researchers have attempted to learn or initialize HMM parameters based on additional sources of patient data (Gupta, 2019; Zhou et al., 2019). Such modifications typically require multiple estimation steps (e.g., Zhou et al., 2019) or changes to parameter interpretation (e.g., Igl et al., 2018). This is because the standard algorithm for fitting HMMs, the Baum-Welch algorithm (Baum & Petrie, 1966), maximizes the likelihood of a data sequence without consideration of additional covariates.
We introduce Hidden Markov Recurrent Neural Networks (HMRNNs) - neural networks that mimic the computation of hidden Markov models while allowing for substantial modularity with other predictive networks. Unlike past work combining neural networks and HMMs (e.g., Bridle, 1990), HMRNNs are designed to maximize the most commonly-used HMM fit criterion - the likelihood of the data given the parameters. In doing so, our primary contributions are as follows: (1) We prove how recurrent neural networks (RNNs) can be formulated to optimize the same likelihood function as HMMs, with parameters that can be interpreted as HMM parameters (section 3); (2) We empirically demonstrate that our model yields statistically similar parameter solutions compared with the Baum-
Welch algorithm (section 4.1); (3) We demonstrate our model’s utility in a disease progression application, in which combining it with other predictive neural networks improves predictive accuracy and offers unique parameter interpretations not afforded by simple HMMs (section 4.2).

heading: 2 RELATED WORK

text:
A small number of studies have attempted to formally model HMMs in a neural network context. Wessels & Omlin (2000) proposes using neural networks to approximate Gaussian emission distributions in HMMs; however, their method requires pre-training of the HMM. Similar to our work, Bridle (1990) demonstrates how HMMs can be reduced to recurrent neural networks for speech recognition, though it requires that neurons be computed via products (rather than sums), which are not commonly used in modern neural networks. This model also maximizes the mutual information between observations and hidden states; this is a commonly used criterion in speech recognition, but less common than likelihood maximization in other domains (e.g., disease progression modeling). Lastly, Bridle (1990) and Wessels & Omlin (2000) present only theoretical justification, with no empirical comparisons with the Baum-Welch algorithm.
A limited number of studies have also explored connections between neural networks and Markov models in the healthcare domain. For instance, Nemati et al. (2016) employs a discriminative hidden Markov model to estimate ‘hidden states’ underlying patients’ ICU measurements, though these hidden states are not mathematically equivalent to HMM latent states. Estebanez et al. (2012) compares HMM and neural network effectiveness in training a robotic surgery assistant, while Baucum et al. (2020) proposes a generative neural network for modeling ICU patient health based on the mathematical intuition of the HMM. Although these studies showcase the of value of pairing neural networks and Markov models in the healthcare domain, they differ from our approach of directly formulating HMMs as neural networks, which maintains the interpretability of HMMs while allowing for joint estimation of the HMM with other predictive models.
In summary, studies have shown the promise of incorporating elements of HMMs into deep learning tasks, there are no existing methods for optimizing HMM log-likelihood in a neural network context. While past works have also used gradient descent to learn HMM parameters (e.g. Yildirim et al., 2015), we demonstrate how specifically implementing HMMs as neural networks allows additional data sources (e.g., patient covariates) to steer model estimation to better-fitting solutions. We thus develop the first neural network formulation of an HMM that maximizes the observed data likelihood, employs widely-used neural network operations, and compares favorably to the Baum-Welch algorithm when tested on real-world datasets.

heading: 3 METHODS

text: In this section, we briefly review HMM preliminaries, formally define the HMRNN, and prove that it optimizes the same likelihood function as a corresponding HMM.

heading: 3.1 HMM PRELIMINARIES

text:
Formally, an HMM models a system over a given time horizon T , where the system occupies a hidden state xt ∈ S = {1, 2, . . . , k} at any given time point t ∈ {0, 1, . . . , T}; that is, xt = i indicates that the system is in the ith state at time t. For any state xt ∈ S and any time point t ∈ {0, 1, . . . , T}, the system emits an observation according to an emission distribution that is uniquely defined for each state. We consider the case of categorical emission distributions, which are commonly used in healthcare (e.g., Liu et al., 2015; Leon, 2015; Ayer et al., 2012; Stanculescu et al., 2013). These systems emit one of c distinct observations at each time point; that is, for any time t, we observe yt ∈ O, where |O| = c and O = {1, . . . , c}. Thus, an HMM is uniquely defined by a k-length initial probability vector π, k × k transition matrix P , and k × c emission matrix Ψ. Entry i in the vector π is the probability of starting in state i, row i in the matrix P is the state transition probability distribution from state i, and row i of the matrix Ψ is the emission distribution from state i.
HMMs are fit via the Baum-Welch algorithm, which identifies the parameters that (locally) maximize the likelihood of the observed data (Baum & Petrie, 1966). The likelihood of an observation sequence y is a function of an HMM’s initial state distribution (π), transition probability matrix (P ), and emission matrix (Ψ) (Jurafsky & Martin, 2009). Let diag(Ψi) be a k × k diagonal matrix with the ith column of Ψ as its entries - i.e., the probabilities of observation i from each state. We then have
Pr(y|π,P ,Ψ) = π> · diag(Ψy0) · ( T∏
i=1
P · diag(Ψyi)) · 1k×1. (1)
The likelihood function can also be expressed in terms of αt(i), the probability of being in state i at time t and having observed {y0, ..., yt}. We denote αt as the (row) vector of all αt(i) for i ∈ S, with
αt = π > · diag(Ψy0) · ( t∏ i=1 P · diag(Ψyi)) (2)
for t ∈ {1, ..., T}, with α0 = π> · diag(Ψy0). Note that equation 3.1 also implies that αt = αt−1 · P · diag(Ψyt) for t ∈ {1, ..., T}, and that
Pr(y|π,P ,Ψ) = αT · 1k×1. (3)

heading: 3.2 DEFINITION OF HIDDEN MARKOV RECURRENT NEURAL NETWORKS (HMRNNS)

text:
An HMRNN is a recurrent neural network whose parameters directly correspond to the initial state, transition, and emission probabilities of an HMM. As such, training an HMRNN optimizes the joint log-likelihood of the N T -length observation sequences given these parameters. Definition 3.1. An HMRNN is a recurrent neural network with trainable parameters π (a k-length stochastic vector), P (a k × k stochastic matrix), and Ψ (a k × c stochastic matrix). It is trained on T + 1 input matrices of size N × c, denoted by Yt for t ∈ {0, 1, . . . , T}, where the n-th row of matrix Yt is a one-hot encoded vector of observation y (n) t for sequence n ∈ {1, . . . , N}. The HMRNN consists of an inner block of hidden layers that is looped T + 1 times (for t ∈ {0, 1, . . . , T}), with each loop containing hidden layers h(t)1 , h (t) 2 , and h (t) 3 , and a c-length input layer h (t) y through which the input matrix Yt enters the model. The HMRNN has a single output unit o(T ) whose value is the joint negative log-likelihood of the N observation sequences under an HMM with parameters π, P , and Ψ; the summed value of o(T ) across all N observation sequences is also the loss which is minimized through any neural network optimizer (e.g., gradient descent).
Layers h(t)1 , h (t) 2 , h (t) 3 , and o (T ) are defined in the following equations. Note that the block matrix in equation (5) is a c× (kc) block matrix of c 11×kvectors , arranged diagonally, while the block matrix in equation (6) is a (kc)× k row-wise concatenation of c k × k identity matrices.
h (t) 1 =
{ π>, t = 0,
h (t−1) 3 P , t > 0.
(4)
h (t) 2 = ReLu ( h (t) 1 [diag(Ψ1) . . . diag(Ψc)] + Yt [ 11×k . . . 01×k . . . . . . . . .
01×k . . . 11×k
] − 1n×(kc) ) (5)
h (t) 3 = h (t) 2 [Ik . . . Ik] > (6)
o(T ) = − log(h(T )3 1k×1). (7)
Fig. 1 outlines the structure of the HMRNN. Intuitively, operations within each recurrent block mimic matrix multiplication by diag(Ψyt) (i.e., h (t) 3 = h (t) 1 diag(Ψyt)), while connections between blocks mimic multiplication by P . In each block, layer h(t)1 contains k units that represent the joint probability of being in each state 1–k and all observations up to time t - note that this is equivalent to the α values in traditional HMM notation. Layer h(t)2 , expands each unit in h (t) 1 into c units via connections with weights Ψi,j for i ∈ {1, . . . , k} and j ∈ {1, . . . , c}, resulting in k · c units; this is equivalent to multiplying h(t)1 (in row-vector form) by a column-wise concatenation of diag(Ψj)
for j ∈ {1, . . . , c}. The resulting units represent the probabilities of all possible state/outcome combinations at time t.
Each column in Yt is connected to all units in h (t) 2 that correspond to that column’s observation, with connection weights set to 1. A bias of −1 and subsequent ReLu activation are then applied to layer h (t) 2 ; this leaves the units that correspond to Yt unchanged, while all other units (i.e., probabilities for non-occurring observations) are made negative by the −1 bias, then forced to zero through the ReLu activation. Thus, layer h(t)2 identifies the joint probability of being in each state and observing Yt. Layer h(t)3 then sums across all c units for each of the k states (all of which are zero except for those corresponding to Yt), yielding k unitsrepresenting the probabilities of being in each state given all previous observations and the observation at time t, i.e., h(t)3 = h (t) 1 diag(Ψyt).
We then apply a fully-connected layer of weights to transform h(3)3 to h (t+1) 1 , equivalent to matrix multiplication by P , i.e., h(t+1)1 = h (t) 3 P , forcing the rows of P to sum to 1. Lastly, activations for h (T ) 3 are summed into unit o
(T ) and subject to a negative logarithmic activation function, yielding the negative log-likelihood of the data given the model’s parameters. Note that the hidden layers may suffer from underflow for long sequences. This can be addressed by normalizing layer h(t)3 to sum to 1 at each time point, then simply subtracting the logarithm of the normalization term (i.e., the log-sum of the activations) from network’s output − log(o(T )). The HMRNN is a special case of an RNN, due to its time-dependent layers and shared weights between time points. Note that its use of recurrent blocks (each of which containing three layers) differs from many RNNs that use only a single layer at each time point; this distinction allows the HMRNN to mimic HMM computations of an HMM at each time point.

heading: 3.3 PROOF OF HMM/HMRNN EQUIVALENCE

text:
We now formally establish that the HMRNN’s output unit, o(T ), is the negative log-likelihood of an observation sequence under an HMM with parameters π, P , and Ψ. We prove this for the case of N = 1 and drop notational dependence on n (i.e., we write y(1)t as yt), though extension to N > 1 is trivial since the log-likelihood of multiple independence sequences is simply the sum of their individual log-likelihoods. We first rely on the following lemma.
Lemma 3.1. If all units in h(t)1 (j) are between 0 and 1 (inclusive), then h (t) 3 = h (t) 1 diag(Ψyt).
Proof. Let h(t)1 (j) and h (t) 3 (j) represent the jth units of layer h (t) 1 and h (t) 3 , respectively. Showing h (t) 3 = h (t) 1 diag(Ψyt) is equivalent to showing that h (t) 3 (j) = Ψj,yth (t) 1 (j) for j ∈ {1, .., k}.
To show this, recall that h(t)2 contains k×c units, which we index with a tuple (l,m) for l ∈ {1, . . . , c} and m ∈ {1, . . . , k}. The connection matrix between h(t)1 and h (t) 2 is [ diag(Ψ1) . . . diag(Ψc)]. Thus, the connection between units h(t)1 (j) and h (t) 2 (l,m) is Ψj,l when j = m, and equals 0 otherwise. Also recall that matrix Yt enters the model through a c-length input layer h (t) y , where the jth unit is 1 when yt = j, and equals 0 otherwise. This layer is connected to h (t) 2 by a c× (kc) diagonal block matrix of c (1× k) row vectors of ones. Thus, the connection between the jth unit of this input layer and unit (l,m) of h(t)2 is 1 when j = l, and equals 0 otherwise. Lastly, a bias of −1 is added to all units in h(t)2 , which is then subject to a ReLu activation, resulting in the following expression for each unit in h(t)2 :
h (t) 2 (l,m) = ReLu(Ψm,l · h (t) 1 (m) + h (t) y (l)− 1). (8)
Because h(t)y (l) is 1 when yt = l, and equals 0 otherwise, then if all units in h (t) 1 are between 0 and 1, we have
h (t) 2 (l,m) =
{ ReLu(Ψm,l · h(t)1 (m)) = Ψm,l · h (t) 1 (m), j = yt,
ReLu(Ψm,l · h(t)1 (m)− 1) = 0, otherwise . (9)
The connection matrix between h(t)2 and h (t) 3 is a (kc)× k row-wise concatenation of k × k identity matrices; thus, the connection between h(t)2 (l,m) and h (t) 3 (j) is 1 if j = m, and 0 otherwise. Hence,
h (t) 3 (j) = c∑ j=0 h (t) 2 (l, j) = Ψj,yt · h (t) 1 (j). (10)
Thus, h(t)3 = h (t) 1 diag(Ψyt).
Theorem 3.1. An HMRNN with parameters π (1× k stochastic vector), P (k× k stochastic matrix), and Ψ (k×c stochastic matrix), and with layers defined as in equations (4-7), produces output neuron o(T ) for sequence n ∈ {1, . . . , N} whose value is the negative log-likelihood of a corresponding HMM as defined in equation (1).
Proof. Note that, based on Lemma 3.1 and equation 4, h(t)3 = h (t−1) 3 · P · diag(Ψyt) for t ∈ {1, ..., T}, assuming that h(t)1 (j) ∈ [0, 1] for j ∈ {1, .., k}. Since αt = αt−1 · P · diag(Ψyt), then if h(t−1)3 = αt−1, then h (t) 1 (j) ∈ [0, 1] for j ∈ {1, .., k} and therefore h (t) 3 = αt. We show the initial condition that h(0)3 = α0, since h (0) 1 = π > implies that h(0)3 = π > · diag(Ψy0) = α0. Therefore, by induction, h(T )3 = αT , and o (T ) = − log(αT · 1k×1), which is the logarithm of the HMM likelihood from equation 3.

heading: 4 EXPERIMENTS AND RESULTS

text: In this section, we compare HMRNNs to Baum-Welch through computational experiments with synthetic data and a case study of Alzheimer’s disease patients. The experiment in section 4.1 demonstrates that Baum-Welch and and the HMRNN yield similar parameter solutions and predictive accuracy. Section 4.2 demonstrates how augmenting the HMRNN with additional neural networks yields HMM parameters with improved predictive ability on real-world data.

heading: 4.1 EMPIRICAL VALIDATION OF HMRNN

text:
We demonstrate that an HMRNN trained via gradient descent yields statistically similar solutions to the Baum-Welch algorithm. We show this with synthetically-generated observations sequences for which the true HMM parameters are known, allowing us to validate Theorem 3.1 empirically.
We simulate systems with state spaces S = 1, 2, ..., k that begin in state 1, using k = 5, 10, or 20 states. These state sizes are consistent with disease progression HMMs, which often involve less than 10 states (Zhou et al., 2019; Jackson et al., 2003; Sukkar et al., 2012; Martino et al., 2020). We assume that each state ‘corresponds’ to one observation, implying the same number of states and observations (c = k). The probability of correctly observing a state (P (yt = xt)) is ψii, which is the diagonal of Ψ and is the same for all states. We simulate systems with ψii = 0.6, 0.75, and 0.9.
We test three variants of the transition probability matrix P . Each is defined by their same-state transition probability pii, which is the same for all states. For all P the probability of transitioning to higher states increases with state membership; this is known as ‘increasing failure rate’ and is a common property for Markov processes. As pii decreases, the rows of P stochastically increase, i.e., lower values of pii imply a greater chance of moving to higher states. We use values of pii = 0.4, 0.6, and 0.8, for 27 total simulations (k = {5, 10, 20} ×Ψii = {0.6, 0.75, 0.9} × pii = {0.4, 0.6, 0.8}). For each of the 27 simulations, we generate 100 trajectories of length T = 60; this time horizon might practically represent one hour of data collected each minute or two months of data collected each day. Initial state probabilities are fixed at 1 for state 1 and 0 otherwise. Transition parameters are initialized based on the observed number of transitions in each dataset, using each observation as a proxy for its corresponding state. Since transition probabilities are initialized assuming no observation error, the emission matrices are correspondingly initialized using ψii = 0.95 (with the remaining 0.05 distributed evenly across all other states). In practice, gradient descent on the HMRNN rarely yields parameter values that are exactly zero. To facilitate comparability between the HMRNN and Baum-Welch, we post-process all HMRNN results with one iteration of the Baum-Welch algorithm, which forces low-probability entries to zero. For Baum-Welch and HMRNN, training ceased when all parameters ceased to change by more than 0.001. For each simulation, we compare Baum Welch’s and the HMRNN’s average Wasserstein distance between the rows of the estimated and ground truth P and Ψ matrices. This serves as a measure of each method’s ability to recover the true data-generating parameters. We also compare the Baum-Welch and HMRNN solutions’ log-likelihoods using a separate hold-out set of 100 trajectories.
Across all simulations, the average Wasserstein distance between the rows of the true and estimated transition matrices was 0.191 for Baum-Welch and 0.178 for HMRNN (paired t-test p-value of 0.483). For the emission matrices, these distances were 0.160 for Baum-Welch and 0.137 for HMRNN (paired t-test p-value of 0.262). This suggests that Baum-Welch and the HMRNN recovered the ground truth parameters with statistically similar degrees of accuracy. This can be seen in Figure 2, which presents the average estimated values of pii and ψii under each model. Both models’ estimated pii values are, on average, within 0.05 of the ground truth values, while they tended to estimate ψii values of around 0.8 regardless of the true ψii. Note that, while Baum-Welch was slightly more accurate at estimating pii and ψii, the overall distance between the ground truth and estimated parameters did not significantly differ between Baum-Welch and the HMRNN.
For each simulation, we also compute the log-likelihood of a held-out set of 100 sequences under the Baum-Welch and HMRNN parameters, as a measure of model fit. The average holdout loglikelihoods under the ground truth, Baum-Welch, and HMRNN parameters are -9250.53, -9296.03, and -9303.27, respectively (paired t-test p-value for Baum-Welch/HMRNN difference of 0.440). Thus, Baum-Welch and HMRNN yielded similar degrees of model fit on held-out sequence data.

heading: 4.2 ALZHEIMER’S DISEASE SYMPTOM PROGRESSION APPLICATION

text:
We also demonstrate how incorporating additional data into an HMRNN can improve parameter fit and offer novel clinical interpretations. We test our HMRNN on clinical data from n = 426 patients with mild cognitive impairment (MCI), collected over the course of three (n = 91), four (n = 106), or five (n = 229) consecutive annual clinical visits (Initiative). Given MCI patients’ heightened risk of Alzheimer’s, modeling their symptom progression is of considerable clinical interest (Hirao et al., 2005; Hansson et al., 2010; Rabin et al., 2009). We analyze patients’ overall cognitive functioning based on the Mini Mental Status Exam (MMSE; Folstein et al., 1975).
MMSE scores range from 0 to 30, with scores of 27−30 indicating no cognitive impairment, 24−26 indicating borderline cognitive impairment, and 17−23 indicating mild cognitive impairment (Chopra et al., 2007; Monroe & Carter, 2012). Scores below 17 were infrequent (1.2%) and were treated as scores of 17 for analysis. We use a 3-state latent space S = {0, 1, 2}, with xt = 0 representing
Estimated vs. Ground Truth Parameters, Baum-Welch vs. HMRNN
‘no cognitive impairment,’ xt = 1 representing ‘borderline cognitive impairment,’ and xt = 2 representing ‘mild cognitive impairment.’ The observation space is O = {0, 1, 2}, using yt = 0 for scores of 27− 30, yt = 1 for scores of 24− 26, and yt = 2 for scores of 17− 23. To showcase the benefits of the HMRNN’s flexibility, we augment the HMRNN through two substantive modifications. First, the initial state probabilities in the augmented HMRNN are predicted from patients’ gender, age, degree of temporal lobe atrophy (Hua et al., 2008), and amyloid-beta 42 levels (Aβ42, a relevant Alzheimer’s biomarker (Canuet et al., 2015; Blennow, 2004)), using a single-layer neural network. Second, at each time point, the probability of being in the most impaired state, h(1)t (2), is used to predict concurrent scores on the Clinical Dementia Rating (CDR, Morris, 1993), a global assessment of dementia severity, allowing another clinical metric to inform estimation. We use a single connection and sigmoid activation to predict patients’ probability of receiving a CDR score above 0.5 (corresponding to ‘mild dementia’).
We compare HMRNN and Baum-Welch parameter solutions’ ability to predict patients’ final MMSE score categories from their initial score categories, using 10-fold cross-validation. We evaluate performance using weighted log-loss L (Guerrero-Pena et al., 2018; Stelmach & Chlebus, 2020), i.e., the log-probability placed on each final MMSE score category averaged across score categories. This metric accounts for class imbalance and rewards models’ confidence in their predictions, an important component of medical decision support (Bussone et al., 2015; Lim & Dey, 2010). We also report p̄, the average probability placed on patients’ final MMSE scores (computed from L). We train all models using a relative log-likelihood tolerance of 0.001% (we do not use parameter convergence since the number of parameters differs between models). Converge runtimes for Baum-Welch and the HMRNN are 2.89 seconds and 15.24 seconds, respectively.
Model results appear in Table 1. Note that the HMRNN’s weighted log-loss L is significantly lower than Baum-Welch’s (paired t-test p-value= 2.396× 10−6), implying greater predictive performance. This is supported by Figure 3, which shows p̄, the average probability placed on patients’ final MMSE scores by score category. Note that error bars represent marginal sampling error and do not represent statistical comparisons between Baum-Welch and HMRNN. Interestingly, the HMRNN yields lower transition probabilities and lower MMSE accuracy (i.e., lower diagonal values of Ψ) than Baum-Welch, suggesting that score changes are more likely attributable to testing error as opposed to true state changes.

heading: 5 DISCUSSION

text:
We outline a flexible approach for HMM estimation using neural networks. The HMRNN produces statistically similar solutions to the Baum-Welch algorithm in its standard form, yet can be combined
with other neural networks to improve predictive accuracy when additional data is available. In our Alzheimer’s case study, augmenting an HMRNN with two predictive networks improves the predictive performance of its parameter solution compared with Baum-Welch. Interestingly, the HMRNN yields a clinically distinct parameter interpretation compared with Baum-Welch, predicting relatively poor diagnostic accuracy for the ‘borderline’ and ‘mild’ cognitive impairment states of the MMSE. This suggests that fewer diagnostic categories might improve the MMSE’s utility, which is supported by existing MMSE research (e.g., Monroe & Carter, 2012), and suggests the HMRNN might also be used to improve the interpretability of HMM parameter solutions.
In addition to demonstrating the HMRNN’s utility in a practical setting, we also make a theoretical contribution by formulating discrete-observation HMMs as a special case of RNNs and by proving coincidence of their likelihood functions. Unlike past approaches, our formulation relies only on matrix multiplication and nonlinear activations, and is designed for generalized use by optimizing for maximum likelihood, which is widely used for HMM training.
Future work may formally assess the time complexity of the HMRNN formulation. The HMRNN converges slower than Baum-Welch for the Alzheimer’s case study. Yet packages for training neural networks are designed to efficiently handle large datasets, and therefore the HMRNN may converge faster than Baum-Welch for problems with larger sample sizes. Yet since sequence lengths in healthcare are often considerably shorter than in other domains that employ HMMs (e.g., speech analysis), runtimes will likely be reasonable for many healthcare datasets. We also limited our case study to disease progression; future work might explore the HMRNN in other healthcare domains.


summaries: This manuscript presents a novel dimensionality reduction method, called Factorized Linear Discriminant Analysis. The method starts from a real problem in neurobiology, and tries to link expression levels of neural genes to phenotypes. In particular, the main goal of the proposed technique is to find linear projections of the genes expression which vary maximally with one phenotypical aspect and minimally with the others. The approach is evaluated using a synthetic example and a real case study involving Drosophila T4/T5 cells.

summaries: This manuscript describes a generalization of ANOVA that is intended to be used in the interpretation of single-cell RNA-seq data. The method requires specification of an orthogonal, discrete categorization of cells, nominally by phenotype.  The method then linearly factorizes the observed gene expression values into features and their interactions, relative to the phenotypic categories.  The factors can then be used to help interpret the categories, especially in conjunction with a regularizer to reduce the number of genes involved in the factors.

summaries: The paper studies a very important problem in gene data analysis. The proposed method is technically sound. The method is intuitive in its idea and easy to implement. The results are interpretable. And according to the experimental evaluations, the proposed method is consistent to existing biological observations and could further identify unknown genetic targets. Therefore, it, potentially, has insightful scientific implications.

abstractText: A central goal in neurobiology is to relate the expression of genes to the structural and functional properties of neuronal types, collectively called their phenotypes. Single-cell RNA sequencing can measure the expression of thousands of genes in thousands of neurons. How to interpret the data in the context of neuronal phenotypes? We propose a supervised learning approach that factorizes the gene expression data into components corresponding to individual phenotypic characteristics and their interactions. This new method, which we call factorized linear discriminant analysis (FLDA), seeks a linear transformation of gene expressions that varies highly with only one phenotypic factor and minimally with the others. We further leverage our approach with a sparsity-based regularization algorithm, which selects a few genes important to a specific phenotypic feature or feature combination. We applied this approach to a single-cell RNA-Seq dataset of Drosophila T4/T5 neurons, focusing on their dendritic and axonal phenotypes. The analysis confirms results obtained by conventional methods but also points to new genes related to the phenotypes and an intriguing hierarchy in the genetic organization of these cells.

authors:
Yoshua Bengio
Aaron Courville
Pascal Vincent

title: Representation Learning: A Review and New Perspectives

authors:
Hung-I. Harry Chen
Yufang Jin
Yufei Huang
Yidong Chen

title: Detection of high variability in gene expression from single-cell RNA-seq profiling

venue: BMC Genomics,

authors:
John P. Cunningham
Zoubin Ghahramani

title: Linear Dimensionality Reduction: Survey, Insights, and Generalizations

venue: Journal of Machine Learning Research,

authors:
Thomas Euler
Silke Haverkamp
Timm Schubert
Tom Baden

title: Retinal bipolar cells: Elementary building blocks of vision

venue: Nature Reviews Neuroscience,

authors:
Fangxiang Feng
Xiaojie Wang
Ruifan Li

title: Cross-modal Retrieval with Correspondence Autoencoder

venue: In Proceedings of the 22nd ACM International Conference on Multimedia, MM

title: The Use of Multiple Measurements in Taxonomic Problems

venue: Annals of Eugenics,

title: The Correlation between Relatives on the Supposition of Mendelian Inheritance

venue: Royal Society of Edinburgh],

authors:
Rohan Gala
Nathan Gouwens
Zizhen Yao
Agata Budzillo
Osnat Penn
Bosiljka Tasic
Gabe Murphy
Hongkui Zeng
Uygar Sümbül

title: A coupled autoencoder approach for multi-modal analysis of cell types

venue: Advances in Neural Information Processing Systems

authors:
Benyamin Ghojogh
Fakhri Karray
Mark Crowley

title: Eigenvalue and Generalized Eigenvalue Problems: Tutorial. arXiv:1903.11240 [cs, stat], March 2019

authors:
Anton Arkhipov
Amy Bernard
John W. Phillips
Hongkui Zeng
Christof Koch

title: Classification of electrophysiological and morphological neuron types in the mouse visual cortex

venue: Nat. Neurosci.,

authors:
Josef Sulc
Susan M. Sunkin
Aaron Szafer
Michael Tieu
Amy Torkelson
Jessica Trinh
Herman Tung
Wayne Wakeman
Katelyn Ward
Grace Williams
Zhi Zhou
Jonathan T. Ting
Anton Arkhipov
Uygar Sümbül
Ed S. Lein
Christof Koch
Zizhen Yao
Bosiljka Tasic
Jim Berg
Gabe J. Murphy
Hongkui Zeng

title: Integrated Morphoelectric and Transcriptomic Classification of Cortical GABAergic Cells

authors:
G.E. Hinton
R.R. Salakhutdinov

title: Reducing the Dimensionality of Data with

venue: Neural Networks. Science,

authors:
Aapo Hyvärinen
Juha Karhunen
Erkki Oja

title: Independent Component Analysis. WileyInterscience, New York, 1st edition

title: Principal Component Analysis

authors:
Theofanis Karaletsos
Serge Belongie
Gunnar Rätsch

title: Bayesian representation learning with oracle constraints. arXiv:1506.05011 [cs, stat], March 2016

venue: Comment: 16 pages,

authors:
Durk P Kingma
Shakir Mohamed
Danilo Jimenez Rezende
Max Welling

title: Semi-supervised Learning with Deep Generative Models

authors:
Tejas D Kulkarni
William F. Whitney
Pushmeet Kohli
Josh Tenenbaum

title: Deep Convolutional Inverse Graphics Network

authors:
Yerbol Z Kurmangaliyev
Juyoun Yoo
Samuel A LoCascio
S Lawrence Zipursky

title: Modular transcriptional programs separately define axon and dendrite connectivity. eLife, 8:e50822, November 2019

title: Decomposition using maximum autocorrelation factors

venue: Journal of Chemometrics,

authors:
Evan Z. Macosko
Anindita Basu
Rahul Satija
James Nemesh
Karthik Shekhar
Melissa Goldman
Itay Tirosh
Allison R. Bialas
Nolan Kamitaki
Emily M. Martersteck
John J. Trombetta
David A. Weitz
Joshua R. Sanes
Alex K. Shalek
Aviv Regev
Steven A. McCarroll

title: Highly parallel genome-wide expression profiling of individual cells using nanoliter

venue: droplets. Cell,

title: Discriminant Analysis and Statistical Pattern Recognition

authors:
B. Moghaddam
Ming-Hsuan Yang

title: Learning gender with support faces

venue: IEEE Transactions on Pattern Analysis and Machine Intelligence,

authors:
Charles R. Gerfen
Christof Koch
Amy Bernard
Chinh Dang
Allan R. Jones
Hongkui Zeng

title: A mesoscale connectome of the mouse

venue: brain. Nature,

authors:
Shristi Pandey
Karthik Shekhar
Aviv Regev
Alexander F. Schier

title: Comprehensive Identification and Spatial Mapping of Habenular Neuronal Types Using Single-Cell RNA-Seq

venue: Curr. Biol.,

authors:
Yi-Rong Peng
Karthik Shekhar
Wenjun Yan
Dustin Herrmann
Anna Sappington
Gregory S. Bryman
Tavé van Zyl
Michael Tri. H. Do
Aviv Regev
Joshua R. Sanes

title: Molecular Classification and Comparative Taxonomics of Foveal and Peripheral Cells in Primate

venue: Retina. Cell,

authors:
Filipe Pinto-Teixeira
Clara Koo
Anthony Michael Rossi
Nathalie Neriec
Claire Bertet
Xin Li
Alberto Del-Valle-Rodriguez
Claude Desplan

title: Development of Concurrent Retinotopic Maps in the Fly Motion Detection Circuit

authors:
Karl Ridgeway
Michael C Mozer

title: Learning Deep Disentangled Embeddings With the F-Statistic Loss

authors:
Joshua R. Sanes
Richard H. Masland

title: The types of retinal ganglion cells: Current status and implications for neuronal classification

venue: Annu. Rev. Neurosci.,

authors:
Karthik Shekhar
Sylvain W. Lapan
Irene E. Whitney
Nicholas M. Tran
Evan Z. Macosko
Monika Kowalczyk
Xian Adiconis
Joshua Z. Levin
James Nemesh
Melissa Goldman
Steven A. McCarroll
Constance L. Cepko
Aviv Regev
Joshua R. Sanes

title: COMPREHENSIVE CLASSIFICATION OF RETINAL BIPOLAR NEURONS

venue: BY SINGLE-CELL TRANSCRIPTOMICS. Cell,

title: General Intelligence,” Objectively Determined and Measured

venue: The American Journal of Psychology,

authors:
Kean Ming Tan
Zhaoran Wang
Han Liu
Tong Zhang

title: Sparse generalized eigenvalue problem: Optimal statistical rates via truncated Rayleigh flow

venue: Journal of the Royal Statistical Society: Series B (Statistical Methodology),

authors:
Bosiljka Tasic
Vilas Menon
Thuc Nghi Nguyen
Tae Kyung Kim
Tim Jarsky
Zizhen Yao
Boaz Levi
Lucas T. Gray
Staci A. Sorensen
Tim Dolbeare
Darren Bertagnolli
Jeff Goldy
Nadiya Shapovalova
Sheana Parry
Changkyu Lee
Kimberly Smith
Amy Bernard
Linda Madisen
Susan M. Sunkin
Michael Hawrylycz
Christof Koch
Hongkui Zeng

title: Adult mouse cortical cell taxonomy revealed by single cell transcriptomics

authors:
Tae Kyung Kim
Aaron Szafer
Nick Dee
Peter Groblewski
Ian Wickersham
Ali Cetin
Julie A. Harris
Boaz P. Levi
Susan M. Sunkin
Linda Madisen
Tanya L. Daigle
Loren Looger
Amy Bernard
John Phillips
Ed Lein
Michael Hawrylycz
Karel Svoboda
Allan R. Jones
Christof Koch
Hongkui Zeng

title: Shared and distinct transcriptomic cell types across neocortical

venue: areas. Nature,

authors:
Nicholas M. Tran
Karthik Shekhar
Irene E. Whitney
Anne Jacobi
Inbal Benhar
Guosong Hong
Wenjun Yan
Xian Adiconis
McKinzie E. Arnold
Jung Min Lee
Joshua Z. Levin
Dingchang Lin
Chen Wang
Charles M. Lieber
Aviv Regev
Zhigang He
Joshua R. Sanes

title: Single-cell profiles of retinal neurons differing in resilience to injury reveal neuroprotective genes

authors:
Weiran Wang
Raman Arora
Karen Livescu
Jeff Bilmes

title: On Deep Multi-View Representation Learning: Objectives and Optimization

authors:
Laurenz Wiskott
Terrence J. Sejnowski

title: Slow Feature Analysis: Unsupervised Learning of Invariances

venue: Neural Computation,

authors:
Hongkui Zeng
Joshua R. Sanes

title: Neuronal cell-type classification: Challenges, opportunities and the path forward

authors:
Zhifei Zhang
Yang Song
Hairong Qi

title: Age Progression/Regression by Conditional Adversarial Autoencoder

heading: 1 INTRODUCTION

text:
The complexity of neural circuits is a result of many different types of neurons that specifically connect to each other. Each neuronal type has its own phenotypic traits, which together determine the role of the neuronal type in a neural circuit. Typical phenotypic descriptions of neurons include features such as dendritic and axonal laminations, electrophysiological properties, and connectivity (Sanes & Masland, 2015; Zeng & Sanes, 2017; Gouwens et al., 2019). However, the genetic programs behind these phenotypic characteristics are still poorly understood.
Recent progress in characterizing neuronal cell types and investigating their gene expression, especially with advances in high-throughput single-cell RNA-Seq (Zeng & Sanes, 2017), provides an opportunity to address this challenge. With massive data generated from single-cell RNA-Seq, we now face a computational problem: how to factorize the high-dimensional data into gene expression modules that are meaningful to neuronal phenotypes? Specifically, given phenotypic descriptions of neuronal types, such as their dendritic stratification and axonal termination, can one project the original data into a low-dimensional space corresponding to these phenotypic features and their interactions, and further extract genes critical to each of these components?
Here we propose a new analysis method named factorized linear discriminant analysis (FLDA). Inspired by multi-way analysis of variance (ANOVA) (Fisher, 1918), this method factorizes data into components corresponding to phenotypic features and their interactions, and seeks a linear transformation that varies highly with one specific factor but not with the others. The linear nature of this approach makes it easy to interpret, as the weight coefficients directly inform the relative importance of each gene to each factor. We further introduce a sparse variant of the method, which constrains the number of genes contributing to each linear projection. We illustrate this approach by applying FLDA to a single-cell transcriptome dataset of T4/T5 neurons in Drosophila (Kurmangaliyev et al., 2019), focusing on two phenotypes: dendritic location and axonal lamination.

heading: 2 FACTORIZED LINEAR DISCRIMINANT ANALYSIS (FLDA)

text:
Suppose that we are given gene expression data of single neurons which are typically very highdimensional. These cells are classified into cell types, as a result of clustering in the highdimensional space and annotations based on prior knowledge or verification outcome (Macosko et al., 2015; Tasic et al., 2016; Shekhar et al., 2016; Tasic et al., 2018; Peng et al., 2019). We know the phenotypic traits of each neuronal type, therefore each type can also be jointed defined by the phenotypic features. We want to find an interpretable low-dimensional embedding in which certain dimensions represent factors of phenotypic features or their interactions. This requires that variation along one of the axes in the embedding space causes the variation of only one factor. In reality, this is hard to satisfy due to noise in the data, and we relax the constraint by letting data projected along one axis vary largely with one factor while minimally with the others. In addition, we ask that cells classified as the same type are still close to each other in the embedding space, while cells of different types are far apart.
As a start, let us consider only two phenotypic features of neurons, dendritic stratification, and axonal termination, both of which can be described with discrete categories, such as different regions or layers in the brain (Oh et al., 2014; Euler et al., 2014; Sanes & Masland, 2015; Kurmangaliyev et al., 2019). Suppose that each cell type can be jointly represented by its dendritic location indexed as i and axonal lamination indexed as j, with the number of cells within each cell type nij . This representation can be described using a contingency table (Figure 1A,B). Note here that we allow the table to be partially filled.
Let xijk(k ∈ 1, 2, ...nij) represent the expression values of g genes in each cell (xijk ∈ Rg)). How to find linear projections yijk = uTxijk and zijk = vTxijk that are aligned with features i and j respectively (Figure 1C)? We first asked whether we could factorize, for example, yijk, with respect to components depending on features i and j. Indeed, motivated by the linear factor models used in multi-way ANOVA and the idea of partitioning variance, we constructed an objective function as the following, and found u∗ that maximizes the objective (see detailed analysis in Appendix A):
u∗ = arg max u∈Rg
uTNAu uTMeu (1)
When we have a complete table, and there are a levels for the feature i and b levels for the feature j, we have
NA = MA − λ1MB − λ2MAB (2)
where MA, MB , and MAB are the covariance matrices explained by the feature i, the feature j, and the interaction of them. λ1 and λ2 are hyper-parameters controlling the relative weights of MB
and MAB with respect to MA. Me is the residual covariance matrix representing noise in gene expressions. Formal definitions of these terms are the following:
MA = 1
a− 1 a∑ i=1 b∑ j=1 (mi. −m..)(mi. −m..)T (3)
MB = 1
b− 1 a∑ i=1 b∑ j=1 (m.j −m..)(m.j −m..)T (4)
MAB = 1
(a− 1)(b− 1) a∑ i=1 b∑ j=1 (mij −mi. −m.j + m..)(mij −mi. −m.j + m..)T (5)
Me = 1
N − ab a∑ i=1 b∑ j=1 [ 1 nij nij∑ k=1 (xijk −mij)(xijk −mij)T ] (6)
where
m.. = 1
ab a∑ i=1 b∑ j=1 mij (7)
mi. = 1
b b∑ j=1 mij (8)
m.j = 1
a a∑ i=1 mij (9)
in which
mij = 1
nij nij∑ k=1 xijk (10)
An analogous expression provides the linear projection v∗ for the feature j, and w∗ for the interaction of both features i and j. Similar arguments can be applied to the scenario of a partial table to find u∗ or v∗ as the linear projection for the feature i or j (see Appendix B for mathematical details).
Note that NA is symmetric and Me is positive definite. Therefore the optimization problem is a generalized eigenvalue problem (Ghojogh et al., 2019). When Me is invertible, u∗ is the eigenvector associated with the largest eigenvalue of M−1e NA. In general, if we want to embed xijk into a ddimensional subspace aligned with the feature i (d < a), we can take the eigenvectors with the d largest eigenvalues of M−1e NA, which we call the top d factorized linear discriminant components (FLDs). Since multi-way ANOVA can handle contingency tables with more than two dimensions, our analysis can be easily generalized to more than two features.

heading: 3 SPARSITY-BASED REGULARIZATION OF FLDA

text:
For this domain-specific application in neurobiology, there is particular interest in finding a small group of genes that best determine one of the phenotypic features. This leads to finding axes that have only a few non-zero components. To identify such a sparse solution, we solved the following optimization problem:
u∗ = arg max u∈Rg
uTNAu uTMeu subject to ||u||0 ≤ l (11)
from which the number of non-zero elements of u∗ is less or equal to l.
This is known as a sparse generalized eigenvalue problem, which has three unique challenges, as listed in Tan et al. (2018): first, when the data are very high-dimensional, Me can be singular and non-invertible; second, because of the normalization term uTMeu, many solutions for sparse eigenvalue problems cannot be applied directly; finally, this problem involves maximizing a convex objective over a nonconvex set, which is NP-hard.
To solve it, we used truncated Rayleigh flow (Rifle), a method specifically developed to solve sparse generalized eigenvalue problems. The algorithm of Rifle is composed of two steps (Tan et al., 2018): first, to obtain an initial vector u0 that is close to u∗. We used the solution from the non-sparse FLDA as an initial estimate of u0; second, iteratively, to perform a gradient ascent step on the objective function, and then execute a truncation step that preserves the l entries of u with the largest values and sets the remaining entries to 0. Pseudo-code for this algorithm can be found in Appendix C.
As proved in Tan et al. (2018), if there is a unique sparse leading generalized eigenvector, Rifle will converge linearly to it with the optimal statistical rate of convergence. The computational complexity of the second step is O(lg + g) for each iteration, therefore Rifle scales linearly with g, the dimensionality of the original data. Based on the theoretical proof, to guarantee convergence, the hyperparameter η was selected to be sufficiently small such that ηλmax(Me) < 1, where λmax(Me) is the largest eigenvalue of Me. In our case, the other hyperparameter l, indicating how many genes to be preserved, was empirically selected based on the design of a follow-up experiment. As mentioned later in Results, we chose l to be 20, a reasonable number of candidate genes to be tested in a biological study.

heading: 4 RELATED WORK

heading: 4.1 LINEAR DIMENSIONALITY REDUCTION

text:
FLDA is a method for linear dimensionality reduction (Cunningham & Ghahramani, 2015). Formally, linear dimensionality reduction is defined as the following: given n data points each of g dimensions, X = [x1,x2, ...,xn] ∈ Rg×n, and a choice of reduced dimensionality r < g, optimize an objective function f(.) to produce a linear projection P ∈ Rr×g , and Y = PX ∈ Rr×n is the low-dimensional transformed data.
State-of-the-art methods for linear dimensionality reduction include principal component analysis (PCA), factor analysis (FA), linear multidimensional scaling (MDS), linear discriminant analysis (LDA), canonical correlations analysis (CCA), maximum autocorrelation factors (MAF), slow feature analysis (SFA), sufficient dimensionality reduction (SDR), locality preserving projections (LPP), and independent component analysis (ICA) (Cunningham & Ghahramani, 2015). Some of them are obviously unsuitable to the problem we want to solve, for example, MAF and SFA are developed for data with temporal structures (Larsen, 2002; Wiskott & Sejnowski, 2002), and LPP focuses on the local structures of the data instead of their global organization (He & Niyogi). The remaining approaches can be roughly grouped for either unsupervised or supervised linear dimensionality reduction, and we discuss them separately.

heading: 4.2 UNSUPERVISED METHODS FOR LINEAR DIMENSIONALITY REDUCTION

text:
Unsupervised linear dimensionality reduction, including PCA (Jolliffe, 2002), ICA (Hyvärinen et al., 2001), FA (Spearman, 1904) and more, project data into a low-dimensional space without using supervision labels. They are not suitable to solve our problem, due to the specific characteristics of the gene expression data. The dimensionality of gene expression data is usually very high, with tens of thousands of genes, and expressions of a fair number of them can be very noisy. These noisy genes cause large variance among individual cells, but in an unstructured way. Without supervision signals from the phenotypic features, unsupervised methods tend to select these genes
to construct the low-dimensional space, which offers neither the desired alignment nor a good separation of cell type clusters. To illustrate this, we performed PCA on the gene expression data and compared it with FLDA. Briefly, we solved the following objective to find the linear projection:
u∗ = arg max u∈Rg
uTXXTu
uTu (12)
The outcome of this comparison is shown in Results.

heading: 4.3 SUPERVISED METHODS FOR LINEAR DIMENSIONALITY REDUCTION

text: Supervised linear dimensionality reduction, represented by LDA (Fisher, 1936; McLachlan, 2004) and CCA (Hotelling, 1936; Wang et al., 2016), can overcome the above issue. By including supervised signals of phenotypic features, we can devalue genes whose expressions are non-informative about the phenotypes.

heading: 4.3.1 LINEAR DISCRIMINANT ANALYSIS (LDA)

text:
We name our method FLDA because its objective function has a similar format as that of LDA. LDA also models the difference among data organized in pre-determined classes. Formally, LDA solves the following optimization problem:
u∗ = arg max u∈Rg
uTΣbu uTΣeu (13)
where Σb and Σe are estimates of the between-class and within-class covariance matrices respectively.
Different from FLDA, the representation of these classes is not explicitly formulated as a contingency table composed of multiple features. The consequence is that, when applied to the example problem in which neuronal types are organized into a two-dimensional contingency table with phenotypic features i and j, in general, axes from LDA are not aligned with these two phenotypic features.
However, in the example above, we can perform two separate LDAs for the two features. This allows the axes from each LDA to align with its specific feature. We call this approach “2LDAs”. There are two limitations of this approach: first, it discards information about the component depending on the interaction of the two features which cannot be explained by a linear combination of them; second, it explicitly maximizes the segregation of cells with different feature levels which sometimes is not consistent with a good separation of cell type clusters. Detailed comparisons between LDA, “2LDAs” and FLDA can be found in Results.

heading: 4.3.2 CANONICAL CORRELATION ANALYSIS (CCA)

text:
CCA projects two datasets Xa ∈ Rg×n and Xb ∈ Rd×n to Ya ∈ Rr×n and Yb ∈ Rr×n, such that the correlation between Ya and Yb is maximized. Formally, it tries to maximize this objective:
(u∗,v∗) = arg max u∈Rg,v∈Rd
uT (XaX T a ) − 12XaX T b (XbX T b ) − 12v
(uTu)− 1 2 (vTv)− 1 2
(14)
To apply CCA to our problem, we need to set Xa to be the gene expression matrix, and Xb to be the matrix of d phenotypic features (d = 2 for two features as illustrated later). In contrast with FLDA, CCA finds a transformation of gene expressions aligned with a linear combination of phenotypic features, instead of a factorization of gene expressions corresponding to individual phenotypic features. This difference is quantified and shown in Results.

heading: 5 EXPERIMENTS

heading: 5.1 DATASETS

text:
In order to quantitatively compare FLDA with other linear dimensionality reduction methods, such as PCA, CCA, LDA, and the “2LDAs” approach, we created synthetic datasets. Four types of cells, each containing 25 examples, were generated from a Cartesian product of two features i and j, organized in a 2x2 complete contingency table. Expressions of 10 genes were generated for these cells, in which the levels of Genes 1-8 were correlated with either the feature i, the feature j, or the interactions of them, and the levels of the remaining 2 genes were purely driven by noise (Figure 2A). Details of generating the data can be found in Appendix D.
To illustrate FLDA in analyzing single-cell RNA-Seq datasets for real problems of neurobiology, and demonstrate the merit of our approach in selecting a few important genes for each phenotype, we used a dataset of Drosophila T4/T5 neurons (Kurmangaliyev et al., 2019). T4 and T5 neurons are very similar in terms of general morphology and physiological properties, but they differ by the location of their dendrites in the medulla and lobula, two distinct brain regions. T4 and T5 neurons each contain four subtypes, with each pair of the four laminating their axons in a specific layer in the lobula plate (Figure 3A). Therefore, we can use two phenotypic features to describe these neurons: the feature i indicates the dendritic location at the medulla or lobula; the feature j describes the axonal lamination at one of the four layers (a/b/c/d) (Figure 3B). In this experiment, we focused on the dataset containing expression data of 17492 genes from 3833 cells collected at a defined time during brain development.

heading: 5.2 DATA PREPROCESSING

text: The T4/T5 neuron dataset was preprocessed as previously reported (Shekhar et al., 2016; Peng et al., 2019; Tran et al., 2019). Briefly, transcript counts within each column of the count matrix (genes×cells) were normalized to sum to the median number of transcripts per cell, resulting in the normalized counts Transcripts-per-median or TPMgc for Gene g in Cell c. We used the logtransformed expression data Egc = ln (TPMgc + 1) for further analysis. We adopted a common approach in single-cell RNA-Seq studies that is based on fitting a relationship between mean and coefficient of variation (Chen et al., 2016; Pandey et al., 2018) to select highly variable genes, and performed FLDA on the expression data with only these genes. We preprocessed the data with PCA and kept principal components (PCs) explaining ∼99% of the total variance before running FLDA but not the sparse version of the algorithm. In the experiment below, we set the hyper-parameters λs in Equation (2) to 1.

heading: 5.3 METRICS

text:
We included the following metrics to evaluate our method: A signal-to-noise ratio (SNR) measures how well each discriminant axis separates cell types compared with noise estimated from the vari-
ance within cell type clusters. The explained variance (EV) for each discriminant axis measures how much variance of the feature i or j is explained among the total variance explained by that axis. The mutual information (MI) between each discriminant axis and each feature quantifies how ”informative” an axis is to a specific feature. Built on the calculation of MI, we included the modularity score which measures whether each discriminant axis depends on at most one feature (Ridgeway & Mozer, 2018). The implementation details of these metrics can be found in Appendix E.

heading: 6 RESULTS

text:
To quantitatively compare the difference between FLDA and other alternative models including PCA, CCA, LDA, and “2LDAs”, we measured the proposed metrics from analyses of the synthesized datasets (Figure 2A). Given that the synthesized data were organized in a 2x2 contingency table, each LDA of the “2LDAs” approach could find only one dimension for the specific feature i or j. Therefore, as a fair comparison, we only included the corresponding dimensions in FLDA (FLDi and FLDj) and the top two components of PCA, CCA, and LDA. The overall SNR values normalized by that of LDA and the overall modularity scores were plotted for data generated with different noise levels (Figure 2B,C). The performance of PCA is the worst among all these models because the unsupervised approach cannot prevent the noise from contaminating the signal. The supervised approaches in general have good SNRs, but LDA and CCA suffer from low modularity scores. This is expected because LDA maximizes the separation of cell type clusters but overlooks the alignment of the axes to the feature i or j, and CCA maximizes the correlation to a linear combination of phenotypic features instead of individual ones. By contrast, “2LDAs” achieves the highest modularity scores but has the worst SNR among the supervised approaches, because it tries to maximize the separation of cells with different feature levels, which is not necessarily consistent with maximizing the segregation of cell types. Both the SNR value and the modularity score of FLDA
are close to the optimal, as it not only considers the alignment of axes to different features but also constrains the variance within cell types. A representative plot of the EV and MI metrics of these models is shown in Figure 5, reporting good alignment of axes to either the feature i or j in FLDA and ‘2LDAs”, but not in the others.
A question of significance in neurobiology is whether the diverse phenotypes of neuronal cell types are generated by combinations of modular transcriptional programs, and if so, what is the gene signature for each of the programs. To illustrate the ability of our approach in addressing this problem, we applied FLDA to the dataset of Drosophila T4/T5 neurons. The T4/T5 neurons could be organized in a 2x4 contingency table, therefore, FLDA was able to project the expression data into a subspace of seven dimensions, with one FLD aligned with dendritic location i (FLDi), three FLDs aligned with axonal termination j (FLDj1−3), and the remaining three representing the interaction of both phenotypes (FLDij1−3). We ranked these axes based on their SNR metrics and found that FLDj1, FLDi, and FLDj2 have much higher SNRs than the rest (Figure 3C). Indeed, data representations in the subspace consisting of these three dimensions show a clear separation of the eight neuronal types (Figure 3D). As expected, FLDi teases apart T4 from T5 neurons, whose dendrites are located at different brain regions (Figure 3E). Interestingly, FLDj1 separates T4/T5 neurons into two groups, a/b vs c/d, corresponding to the upper or lower lobula place, and FLDj2 divides them into another two, a/d vs b/c, indicating whether their axons laminate at the middle or lateral part of the lobula plate (Figure 3E,F). Unexpectedly, among these three dimensions, FLDj1 has a much higher SNR than FLDi and FLDj2, whose SNR values are similar. This suggests a hierarchical structure in the genetic organization of T4/T5 neurons: they are first separated into either a/b or c/d types, and subsequently divided into each of the eight subtypes. In fact, this exactly matches the sequence of their cell fate determination during development, as revealed in a previous genetic study (Pinto-Teixeira et al., 2018). Finally, the last discriminant axis of the axonal feature FLDj3 separates the group a/c from b/d, suggesting its role in fine-tuning the axonal depth within the upper or lower lobula plate (Figure 3G).
To seek gene signatures for the discriminant components in FLDA, we applied the sparsity-based regularization to constrain the number of genes with non-zero weight coefficients. Here we set the number to 20, a reasonable number of candidate genes that might be tested in a follow-up biological study. We extracted a list of 20 genes each for the axis of FLDi or FLDj1. The relative importance of these genes to each axis is directly informed by their weight values (Figure 4A,C). Side-by-side, we plotted expression profiles of these genes in the eight neuronal types (Figure 4B,D). For both axes, the genes critical in separating cells with different feature levels are differentially expressed in corresponding cell types. We compared our gene lists with those obtained using conventional methods which were reported in Kurmangaliyev et al. (2019). Consistent with the report, we found indicator genes for dendritic location such as TfAP -2 , dpr2, CG34155, and CG12065, and those for axonal lamination such as klg, bi, pros. In addition, we found genes that were not reported in this previous study. For example, our results suggest that the genes Thor and pHCl -1 are important to the dendritic phenotype, and Lac and Mip are critical to the axonal phenotype. These are promising genetic targets to be tested in biological experiments. Lastly, FLDA allowed us to examine the component that depends on the interaction of both features and identify its gene signature, which provides clues to transcriptional regulation of gene expressions in the T4/T5 neuronal types (Figures 6 and 7).
As a supervised approach, FLDA depends on correct phenotype labels to extract meaningful information. But if the phenotypes are annotated incorrectly, can we use FLDA to raise a flag? We propose a perturbation analysis of FLDA to address this question built on the assumption that among possible phenotype annotations, the projection of gene expression data based on correct labels leads to better metric measurements than incorrect ones. As detailed in Appendix F, we generated three kinds of incorrect labels for the dataset of T4/T5 neurons, corresponding to three common scenarios of mislabeling: the phenotypes of a cell type were mislabeled with those of another type; a singular phenotypic level was incorrectly split into two; two phenotypic levels are incorrectly merged into one. FLDA was applied to gene expressions of T4/T5 neurons but with these perturbed annotations. Proposed metrics such as the SNR value and modularity score were plotted in Figure 8. Indeed, the projection of gene expressions with correct annotation leads to the best SNR value and modularity score compared with incorrect annotations. This implies that this type of perturbation analysis is a useful practice in general: it raises the confidence that the original annotation is correct if FLDA on the perturbed annotations produces lower metric scores.

heading: 7 DISCUSSION

text:
We developed FLDA, a novel supervised linear dimensionality reduction method for understanding the relationship between high-dimensional gene expression patterns and cellular phenotypes. We illustrate the power of FLDA by analyzing a gene expression dataset from Drosophila T4/T5 neurons that are labeled by two phenotypic features, each with multiple levels. The approach allowed us to identify new genes for each of the phenotypic features that were not apparent under conventional methods. Furthermore, we found a hierarchical relationship in the genetic organization of these cells. These findings point the way for new biological experiments.
The approach is motivated by multi-way ANOVA, and thus it generalizes easily to more than two features. Future applications in neurobiology include the analysis of phenotypic characteristics such as electrophysiology and connectivity (Zeng & Sanes, 2017; Gouwens et al., 2019; 2020). More generally FLDA can be applied to any labeled data set for which the labels form a Cartesian product of multiple features. For example, this would include face images that can be jointly labeled by the age, gender, and other features of a person (Moghaddam & Ming-Hsuan Yang, 2002; Zhang et al., 2017).
FLDA factorizes gene expression data into features and their interactions, and finds a linear projection of the data that varies with only one factor but not the others. This provides a modular representation aligned with the factors (Bengio et al., 2012). Ridgeway & Mozer (2018) argued that modularity together with explicitness could define disentangled representations. Our approach is linear, which presents an explicit mapping between gene expressions and phenotypic features, therefore our approach can potentially serve as a supervised approach to disentanglement (Kingma et al., 2014; Kulkarni et al., 2015; Karaletsos et al., 2016).
Compared with other non-linear embedding methods for cell types (Hinton & Salakhutdinov, 2006; Feng et al., 2014; Gala et al., 2019), the linear nature of FLDA makes it extremely easy to interpret the low-dimensional representations, as the weight vector directly informs the relative importance of each gene. To allow the selection of a small set of critical genes, we leveraged our approach with sparse regularization. This makes FLDA especially useful to experimentalists who can take the list of genes and test them in a subsequent round of experiments.

heading: 8 APPENDIX

heading: 8.1 A. OBJECTIVE FUNCTIONS

text: Here we derive the objective functions used in our analysis. Again if xijk(k ∈ 1, 2, ...nij) represents the expression values of g genes in each cell (xijk ∈ Rg)), we seek to find a linear projection yijk = u Txijk that is aligned with the feature i.

heading: 8.1.1 INSPIRATION FROM ANOVA

text:
We asked what is the best way to factorize yijk. Inspired by multi-way ANOVA (Fisher, 1918), we identified three components: one depending on the feature i, another depending on the feature j, and the last one depending on the interaction of both features. We therefore followed the procedures of ANOVA to partition sums of squares and factorize yijk into these three components.
Let us first assume that all cell types defined by i and j contain the same number of cells. With cell types represented by a complete contingency table (Figure 1A), yijk can be linearly factorized using the model of two crossed factors. Formally, the linear factorization is the following:
yijk = µ+ αi + βj + (αβ)ij + ijk (15)
where yijk represents the coordinate of the kth cell in the category defined by i and j; µ is the average level of y; αi is the component that depends on the feature i, and βj is the component that depends on the feature j; (αβ)ij describes the component that depends on the interaction of both features i and j; ijk ∼ N (0, σ2) is the residual of this factorization. Let us say that the features i and j fall into a and b discrete categories respectively. Then without loss of generality, we can require:
a∑ i=1 αi = 0 (16)
b∑ j=1 βj = 0 (17)
a∑ i=1 (αβ)ij = b∑ j=1 (αβ)ij = 0 (18)
Corresponding to these, there are three null hypotheses:
H01 : αi = 0 (19)
H02 : βj = 0 (20)
H03 : (αβ)ij = 0 (21)
Here we want to reject H01 while accepting H02 and H03 in order that yijk is aligned with the feature i.
Next, we partition the total sum of squares. If the number of cells within each cell type category is n, and the total number of cells is N , then we have
a∑ i=1 b∑ j=1 n∑ k=1 (yijk − ȳ...)2 = bn a∑ i=1 (ȳi.. − ȳ...)2 + an b∑ j=1 (ȳ.j. − ȳ...)2
+ n a∑ i=1 b∑ j=1 (ȳij. − ȳi.. − ȳ.j. + ȳ...)2 + a∑ i=1 b∑ j=1 n∑ k=1 (yijk − ȳij.)2 (22)
where ȳ is the average of yijk over the indices indicated by the dots. Equation (22) can be written as
SST = SSA + SSB + SSAB + SSe (23)
with each term having degrees of freedomN−1, a−1, b−1, (a−1)(b−1), andN−ab respectively. Here SSA, SSB , SSAB , and SSe are partitioned sum of squares for the factors αi, βj , (αβ)ij , and the residual.
ANOVA rejects or accepts a null hypothesis by comparing its mean square (the partitioned sum of squares normalized by the degree of freedom) to that of the residual. This is done by constructing F-statistics for each factor as shown below:
FA = MSA MSe = SSA a−1 SSe N−ab
(24)
FB = MSB MSe = SSB b−1 SSe N−ab
(25)
FAB = MSAB MSe =
SSAB (a−1)(b−1)
SSe N−ab
(26)
Under the null hypotheses, the F-statistics follow the F-distribution. Therefore, a null hypothesis is rejected when we observe the value of a F-statistic above a certain threshold calculated from the F-distribution. Here we want FA to be large enough so that we can reject H01, but FB and FAB to be small enough for us to accept H02 and H03. In other words, we want to maximize FA while minimizing FB and FAB . Therefore, we propose maximizing an objective L:
L = FA − λ1FB − λ2FAB (27)
where λ1 and λ2 are hyper-parameters determining the relative weights of FB and FAB compared with FA.

heading: 8.1.2 OBJECTIVE FUNCTIONS UNDER A COMPLETE CONTINGENCY TABLE

text:
When the numbers of cells within categories defined by i and j (nij) are not all the same, the total sum of squares cannot be partitioned as in Equation (22). However, if we only care about distinctions between cell types instead of individual cells, we can use the mean value of each cell type cluster (ȳij.) to estimate the overall average value (ỹ...), and the average value of each category i (ỹi..) or j (ỹ.j.). Therefore, Equation (22) can be modified as the following:
a∑ i=1 b∑ j=1 [ 1 nij nij∑ k=1 (yijk − ỹ...)2] = b a∑ i=1 (ỹi.. − ỹ...)2 + a b∑ j=1 (ỹ.j. − ỹ...)2
+ a∑ i=1 b∑ j=1 (ȳij. − ỹi.. − ỹ.j. + ỹ...)2 + a∑ i=1 b∑ j=1 [ 1 nij nij∑ k=1 (yijk − ȳij.)2]
(28)
where
ȳij. = ∑nij k=1 yijk nij
(29)
ỹi.. =
∑b j=1 ȳij.
b (30)
ỹ.j. = ∑a i=1 ȳij. a
(31)
ỹ... =
∑a i=1 ∑b j=1 ȳij.
ab (32)
If we describe Equation (28) as:
S̃ST = S̃SA + S̃SB + S̃SAB + S̃Se (33)
then following the same arguments, we want to maximize an objective function in the following format:
L =
S̃SA a−1 − λ1 S̃SB b−1 − λ2 ˜SSAB (a−1)(b−1)
S̃Se N−ab
(34)

heading: 8.1.3 OBJECTIVE FUNCTIONS UNDER A PARTIAL CONTINGENCY TABLE

text:
When we have a representation of a partial table, we can no longer separate out the component that depends on the interaction of both features. Therefore, we use another model, a linear model of two nested factors, to factorize yijk, which has the following format:
yijk = µ+ αi + βj(i) + ijk (35)
Note that we now have βj(i) instead of βj + (αβ)ij . In this model, we identify a primary factor, for instance, the feature denoted by iwhich falls into a categories, and the other (indexed by j) becomes a secondary factor, the number of whose levels bi depends on the level of the primary factor. We merge the component depending on the interaction of both features into that of the secondary factor as βj(i).
Similarly, we have
a∑ i=1 bi∑ j=1 [ 1 nij nij∑ k=1 (yijk − ỹ...)2] = a∑ i=1 [ bi∑ j=1 (ỹi.. − ỹ...)2]
+ a∑ i=1 bi∑ j=1 (ȳij. − ỹi..)2 + a∑ i=1 bi∑ j=1 [ 1 nij nij∑ k=1 (yijk − ȳij.)2]
(36)
which can be written as
S̃ST = S̃SA + S̃SB + S̃Se (37)
with degrees of freedom N − 1, a− 1, M − a, and N −M for each of the terms, where M is:
M = a∑ i=1 bi (38)
Therefore, we want to maximize the following objective:
L = S̃SA a−1 − λ S̃SB M−a
S̃Se N−M
(39)

heading: 8.2 B. FLDA WITH A PARTIAL CONTINGENCY TABLE

text:
Here we provide the mathematical details of FLDA under the representation of a partial table. When we have a partial table, if the feature i is the primary feature with a levels, and the feature j is the secondary feature with bi levels, then NA in Equation (1) is defined as follows:
NA = MA − λMB|A (40)
where
MA = 1
a− 1 a∑ i=1 bi∑ j=1 (mi. −m..)(mi. −m..)T (41)
MB|A = 1
M − a a∑ i=1 bi∑ j=1 (mij −mi.)(mij −mi.)T (42)
and M is defined as in Equation (38). Correspondingly, Me in Equation (1) is defined as:
Me = 1
N −M a∑ i=1 b∑ j=1 [ 1 nij nij∑ k=1 (xijk −mij)(xijk −mij)T ] (43)
and
m.. = 1
M a∑ i=1 bi∑ j=1 mij (44)
mi. = 1
bi bi∑ j=1 mij (45)
The remaining mathematical arguments are the same as those for the complete table. In this scenario, because we don’t observe all possible combinations of features i and j, we cannot find the linear projection for the interaction of both features.

heading: 8.3 C. PSEUDO-CODE FOR RIFLE

text:
We show pseudo-code for the Rifle algorithm as follows:
procedure RIFLE(NA,Me,u0, l, η) . η is the step size t = 1 . t indicates the iteration number while not converge do . Converge when ut ' ut−1
ρt−1 ← uTt−1NAut−1
uTt−1Meut−1
C ← I + ( ηρt−1 )(NA − ρt−1Me) ut ← Cut−1||Cut−1||2 Truncate ut by keeping the top l entries of u with the largest values and setting the rest
entries to 0 ut ← ut||ut||2 t← t+ 1
end while return ut
end procedure

heading: 8.4 D. IMPLEMENTATION DETAILS OF DATA SYNTHESIS

text:
To quantitatively compare FLDA with alternative approaches, we synthesized data of four cell types, each of which contained 25 cells. The four cell types were generated from a Cartesian product of two features i and j, where i ∈ {0, 1} and j ∈ {0, 1}. Expressions of 10 genes were generated for each cell. The expression value of the hth gene in the kth cell of the cell type ij, ghijk was defined as the following:
g1ijk = i+ ijk (46)
g2ijk = j + ijk (47)
g3ijk = i ∧ j + ijk (48)
g4ijk = i ∨ j + ijk (49)
g5ijk = 2i+ ijk (50)
g6ijk = 2j + ijk (51)
g7ijk = 2i ∧ j + ijk (52)
g8ijk = 2i ∨ j + ijk (53)
g9ijk = ijk (54)
g10ijk = 2 + ijk (55)
where
i ∧ j = {
1, if i = 1, j = 1 0, otherwise
(56)
and
i ∨ j = {
0, if i = 0, j = 0 1, otherwise
(57)
were interactions of the two features. Here ijk was driven by Gaussian noise, namely,
ijk ∼ N (0, σ2) (58)
We synthesized datasets of 5 different σ values (σ ∈ {0.2, 0.4, 0.6, 0.8, 1.0}). This was repeated 10 times and metrics for each σ value were calculated as the average across the 10 repeats.

heading: 8.5 E. IMPLEMENTATION DETAILS OF THE METRICS USED IN THE STUDY

text: We measured the following metrics in our experiments:

heading: 8.5.1 SIGNAL-TO-NOISE RATIO (SNR)

text:
Because we care about the separation of cell types, we define the SNR metric as the ratio of the variance between cell types over the variance of the noise, which is estimated from within-cluster variance. For the entire embedding space, given q cell types, if the coordinate of each cell is indicated by c, then we define the overall SNR metric as the following:
SNRoverall = tr(Σqp=1np(c̄p. − c̄..)(c̄p. − c̄..)T ))
tr(Σqp=1Σ np k=1(cpk − c̄p.)(cpk − c̄p.)T )
(59)
where c̄p. is the center of each cell type cluster, and c̄.. is the center of all data points.
Let c denote the embedded coordinate along a specific dimension. The SNR metric for that axis is therefore:
SNR = Σqp=1np(c̄p. − c̄..)2
Σqp=1Σ np k=1(cpk − c̄p.)2
(60)

heading: 8.5.2 EXPLAINED VARIANCE (EV)

text:
We want to know whether the variation of a specific dimension is strongly explained by that of a specific feature. Therefore, we measure, for each axis, how much of the total explained variance is explained by the variance of the feature i or j. Formally, given the embedded coordinate yijk, we calculate the EV as the following:
EVi =
∑a i=1 ∑b j=1 nij(ȳi.. − ȳ...)2∑a
i=1 ∑b j=1 ∑nij k=1(yijk − ȳ...)2
(61)
EVj =
∑a i=1 ∑b j=1 nij(ȳ.j. − ȳ...)2∑a
i=1 ∑b j=1 ∑nij k=1(yijk − ȳ...)2
(62)
where ȳ is the average of yijk over the indices indicated by the dots.

heading: 8.5.3 MUTUAL INFORMATION (MI)

text:
The MI between a discriminant axis u and a feature quantifies how much information of the feature is obtained by observing data projected along that axis. It is calculated as the MI between data representations along the axis y = uTX and feature labels of the data f , where X is the original gene expression matrix:
I(y,f) = H(y) +H(f)−H(y,f) = − ∑ y∈Y p(y) log2 p(y)− ∑ f∈F p(f) log2 p(f)− ∑ y∈Y ∑ f∈F p(y, f) log2 p(y, f) (63)
Here H indicates entropy. To calculate H(y) and H(y,f), we discretize y into 10 bins.

heading: 8.5.4 MODULARITY

text:
Ridgeway & Mozer (2018) argued that in a modular representation, each axis should depend on at most a single feature. Following the arguments in their paper, the modularity score is computed as follows: we first calculate the MI between each feature and each axis (mif denotes the MI between one axis i and one feature f ). If an axis is perfectly modular, it will have high mutual information for only one feature and zeros for the others, we therefore compute a template tif as the following:
tif = { θi, if f = arg maxgmig 0, otherwise
(64)
where θi = maxgmig . We then calculate the deviation from the template as:
δi =
∑ f (mif − tif )2
θ2i (N − 1) (65)
where N is the number of features. The modularity score for the axis i is 1− δi. The mean of 1− δi over i is defined as the overall modularity score.

heading: 8.6 F. IMPLEMENTATION DETAILS OF ANNOTATION PERTURBATION

text:
To evaluate the effect of mislabeling phenotypic levels, we made use of the dataset of T4/T5 neurons, and generated three kinds of perturbation to the original labels:
First, we switched the phenotype labels of T4a neurons with one of the seven other types (T4b, T4c, T4d, T5a, T5b, T5c, T5d). In this scenario, phenotype labels of two cell types were incorrect, but the number of cell type clusters was the same. We had two levels of the dendritic phenotypes (T4/T5), and four levels of the axonal phenotypes (a/b/c/d). Therefore we kept one dimension for the dendritic feature, and three dimensions for the axonal feature.
Second, we merged the axonal phenotypic level a with another level (b/c/d), as an incorrect new level (a+b/a+c/a+d). In this scenario, we had three axonal phenotypes, therefore we kept two dimensions for the axonal feature.
Third, we randomly split each of the four axonal lamination labels (a/b/c/d) into two levels. For instance, among neurons with the original axonal level a, some of them were labeled with a level a1, and the others were labeled with a level a2. In this scenario, we had eight axonal phenotypes (a1/a2/b1/b2/c1/c2/d1/d2), and we kept seven dimensions for the axonal feature.
We performed FLDA on the dataset of T4/T5 neurons but with these perturbed annotations. Metrics from each of the perturbed annotations were measured and compared with those from the original annotation.

heading: 9 ADDITIONAL INFORMATION

title: FACTORIZED LINEAR DISCRIMINANT ANALYSIS FOR PHENOTYPE-GUIDED REPRESENTATION LEARNING OF NEURONAL GENE EXPRESSION DATA