Phylogenetic Factor Analysis

TL;DR

This paper introduces phylogenetic factor analysis (PFA), a Bayesian method that simplifies high-dimensional trait evolution modeling by identifying independent factors, improving computational efficiency and model fit in evolutionary studies.

Contribution

The paper presents PFA, a novel Bayesian approach that infers evolutionary factors, handles mixed data types, and incorporates phylogenetic uncertainty with efficient Gibbs sampling.

Findings

PFA outperforms multivariate diffusion in computational speed.

PFA provides better model fit in evolutionary case studies.

PFA effectively integrates missing data and discrete traits.

Abstract

Phylogenetic comparative methods explore the relationships between quantitative traits adjusting for shared evolutionary history. This adjustment often occurs through a Brownian diffusion process along the branches of the phylogeny that generates model residuals or the traits themselves. For high-dimensional traits, inferring all pair-wise correlations within the multivariate diffusion is limiting. To circumvent this problem, we propose phylogenetic factor analysis (PFA) that assumes a small unknown number of independent evolutionary factors arise along the phylogeny and these factors generate clusters of dependent traits. Set in a Bayesian framework, PFA provides measures of uncertainty on the factor number and groupings, combines both continuous and discrete traits, integrates over missing measurements and incorporates phylogenetic uncertainty with the help of molecular sequences. We develop Gibbs samplers based on dynamic programming to estimate the PFA posterior distribution, over three-fold faster than for multivariate diffusion and a further order-of-magnitude more efficiently in the presence of latent traits. We further propose a novel marginal likelihood estimator for previously impractical models with discrete data and find that PFA also provides a better fit than multivariate diffusion in evolutionary questions in columbine flower development, placental reproduction transitions and triggerfish fin morphometry.

Figures

Tables

Table 1. Table 1 : Log marginal likelihood estimates for the number K 𝐾 K of independent factors driving evolution under a phylogenetic factor analysis (PFA) and a latent multivariate Brownian diffusion (LMBD) model in Aquilegia , and Poeciliidae and multivariate Brownian diffusion (MBD) in Balistidae . The K = 2 𝐾 2 K=2 model for Aquilegia , the K = 3 𝐾 3 K=3 and K = 4 𝐾 4 K=4 model for Poeciliidae and the K = 5 𝐾 5 K=5 model for Balistidae achieve the highest marginal likelihoods.

		Log marginal
	Model	likelihood
Aquilegia	$K=1$	$-385.4$
	$K=2$	$-366.9$
	$K=3$	$-374.3$
	LMBD	$-391.1$
	$K=2$	$-536.0$
Poeciliidae	$K=3$	$-500.7$
	$K=4$	$-501.0$
	$K=5$	$-505.9$
	LMBD	$-592.3$
Balistidae	$K=4$	$-15622.0$
	$K=5$	$-15603.5$
	$K=6$	$-15610.4$
	MBD	$-15673.2$

Table 2. Table 2 : Precision 𝚲 𝚲 \boldsymbol{\Lambda} posterior mean and 95% Bayesian credible interval estimates under the latent factor model for the traits in Aquilegia , in Poeciliidae and in Balistidae . The PFA model explains all of the continuous traits in these models better than a N ​ ( 0 , 1 ) 𝑁 0 1 N(0,1) distribution on the standardized traits.

		Posterior	95% Bayesian
	Trait	mean	credible interval
	Orientation	$2.1$	$[1.0,3.3]$
	Spur length	$4.4$	$[2.0,7.1]$
	Blade length	$3.0$	$[1.4,4.8]$
Aquilegia	Sepal length	$2.6$	$[1.3,4.1]$
	Spur chroma	$4.2$	$[1.8,6.9]$
	Spur hue	$6.2$	$[2.6,10.5]$
	Spur brightness	$2.7$	$[1.2,4.3]$
	Blade chroma	$2.3$	$[1.1,3.7]$
	Blade hue	$2.1$	$[1.0,3.2]$
	Blade brightness	$3.3$	$[1.4,.6]$
	Matrotrophy index	$14.3$	$[5.6,23.2]$
Poeciliidae ($K=3$)	Gonopodium length	$9.3$	$[4.3,16.1]$
	Male body length	$3.5$	$[2.4,4.6]$
	Male body weight	$2.8$	$[1.9,3.7]$
	Female body length	$10.5$	$[5.7,15.5]$
	Female body weight	$15.1$	$[8.0,24.3]$
	Matrotrophy index	$13.8$	$[5.5,22.7]$
Poeciliidae ($K=4$)	Gonopodium length	$9.1$	$[4.4,15.5]$
	Male body length	$3.5$	$[2.3,4.8]$
	Male body weight	$2.8$	$[1.9,3.8]$
	Female body length	$10.5$	$[5.8,15.5]$
	Female body weight	$14.7$	$[8.2,22.5]$

Equations

$$Z_{ij}=\left\{\begin{array}[]{l l}(Y_{ij}-\hat{Y}_{j})/\hat{\sigma}_{j}&\text{ if trait }j\text{ is continuous }\\ Z_{ij}&\text{ if trait }j\text{ is binary or ordinal, }\end{array}\right.$$

$${\gamma }_{j(c-1)}<Z_{ij}\le {\gamma }_{jc}\text{ given }{Y}_{ij}=c$$

$$\mathbf{Z}\sim \text{MN}\left({\mathbf{\mu }}_0,{\mathbf{\Psi }}_{\mathcal{F}}+{\kappa }_0^{-1}\mathbf{J},\mathbf{\Sigma }\right),$$

$$p(\mathbf{Z}|\mathcal{F},\mathbf{\Sigma })=\frac{\text{mbox}exp\left\{-\frac{1}{2}\text{mbox}tr\left[{\mathbf{\Sigma }}^{-1}{\left(\mathbf{Z}-1{\mathbf{\mu }}_0^t\right)}^t{\left({\mathbf{\Psi }}_{\mathcal{F}}+{\kappa }_0^{-1}\mathbf{J}\right)}^{-1}\left(\mathbf{Z}-1{\mathbf{\mu }}_0^t\right)\right]\right\}}{{\left(2\pi \right)}^{NP/2}{\left|\mathbf{\Sigma }\right|}^{N/2}{\left|{\mathbf{\Psi }}_{\mathcal{F}}+{\kappa }_0^{-1}\mathbf{J}\right|}^{P/2}},$$

$$\mathbf{Z}$$

$$\mathbf{F}$$

$$\mathbf{ϵ}$$

$$L_{kj}$$

$${\Lambda }_j$$

$$p(K,\mathbf{F},\mathbf{L},\mathbf{\Lambda },\mathbf{\gamma },\mathcal{F}|\mathbf{Y},\mathbf{S})$$

$$\times p(\mathbf{L}|K)\times p(\mathbf{\Lambda })\times p(\mathbf{\gamma })\times p(K)$$

$$=\left(\int p(\mathbf{Y}|\mathbf{Z},\mathbf{\gamma })p(\mathbf{Z}|K,\mathbf{F},\mathbf{L},\mathbf{\Lambda })\text{mbox}d\mathbf{Z}\right)p(\mathbf{F}|K,\mathcal{F})\times p(\mathcal{F}|\mathbf{S})$$

$$\times p(\mathbf{L}|K)\times p(\mathbf{\Lambda })\times p(\mathbf{\gamma })\times p(K),$$

$${\left(L_{1j},\dots ,L_{k^{\prime }j}\right)}^t|\mathbf{Z},\mathbf{F},\mathbf{\Lambda }\sim \text{MVN}\left({\mathbf{M}}_j^{\left(\mathbf{L}\right)},{\mathbf{V}}_j^{\left(\mathbf{L}\right)}\right)\text{ for }j=1,\dots ,P,$$

$${\mathbf{M}}_j^{\left(\mathbf{L}\right)}$$

$${\mathbf{V}}_j^{\left(\mathbf{L}\right)}$$

$${\Lambda }_j|\mathbf{Z},\mathbf{F},\mathbf{L}\sim \Gamma \left({\alpha }_{\mathbf{\Lambda }}+\frac{N}{2},{\beta }_{\mathbf{\Lambda }}+\frac{1}{2}{\mathbf{e}}_j^t{\left(\mathbf{Z}-\mathbf{FL}\right)}^t\left(\mathbf{Z}-\mathbf{FL}\right){\mathbf{e}}_j\right),$$

$${\mathbf{F}}_{i\cdot }^t|\mathbf{Z},{\mathbf{F}}_{\text{-}i\cdot },\mathbf{L},\mathbf{\Lambda },\mathcal{F}\sim \text{MVN}\left({\mathbf{M}}_i^{\left(\mathbf{F}\right)},{\mathbf{V}}_i^{\left(\mathbf{F}\right)}\right)\text{ for }i=1,\dots ,N,$$

$${\mathbf{M}}_i^{\left(\mathbf{F}\right)}$$

$${\mathbf{V}}_i^{\left(\mathbf{F}\right)}$$

$$\frac{p(K=k|\mathbf{Y},\mathbf{S})}{p(K=k^{\prime }|\mathbf{Y},\mathbf{S})}=\frac{p(\mathbf{Y},\mathbf{S}|K=k)}{p(\mathbf{Y},\mathbf{S}|K=k^{\prime })}\frac{p(K=k)}{p(K=k^{\prime })}.$$

$$q(\beta ,\mathbf{Y},\mathbf{S},\mathbf{\theta })=1(\mathbf{Y}|\mathbf{Z},\mathbf{\gamma })\times h(\mathbf{Y},\mathbf{S},\mathbf{\theta }{)}^{\beta }\times \hat{p}(\mathbf{\theta }{)}^{1-\beta }.$$

$$\int \left\{\underset{\beta \to 0}{\lim }q(\beta ,\mathbf{Y},\mathbf{S},\mathbf{\theta })\right\}\text{mbox}d\mathbf{\theta }=1.$$

$$\int \left\{\underset{\beta \to 0}{\lim }q(\beta ,\mathbf{Y},\mathbf{S},\mathbf{\theta })\right\}\text{mbox}d\mathbf{\theta }=\int 1(\mathbf{Y}|\mathbf{Z},\mathbf{\gamma })\times \hat{p}(\mathbf{\theta })\text{mbox}d\mathbf{\theta }.$$

$${\int }_{\mathbf{\Omega }}\hat{p}(\mathbf{\theta })\text{mbox}d\mathbf{\theta }<1,$$

$$q^{\ast }\left(\beta ,\mathbf{Y},\mathbf{S},\mathbf{\theta }\right)=\left\{1-\left[1-1(\mathbf{Y}|\mathbf{Z},\mathbf{\gamma })\right]\beta \right\}\times h(\mathbf{Y},\mathbf{S},\mathbf{\theta }{)}^{\beta }\times \hat{p}(\mathbf{\theta }{)}^{1-\beta }.$$

$$\int \left\{\underset{\beta \to 0}{\lim }{q}^{\ast }\left(\beta ,\mathbf{Y},\mathbf{S},\mathbf{\theta }\right)\right\}\text{mbox}d\mathbf{\theta }=\int \hat{p}(\mathbf{\theta })\text{mbox}d\mathbf{\theta }=1,$$

$$\frac{\partial }{\partial \beta }\log q^{\ast }\left(\beta ,\mathbf{Y},\mathbf{S},\mathbf{\theta }\right)=-\frac{1-1(\mathbf{Y}|\mathbf{Z},\mathbf{\gamma })}{1-\left[1-1(\mathbf{Y}|\mathbf{Z},\mathbf{\gamma })\right]\beta }+\log h(\mathbf{Y},\mathbf{S},\mathbf{\theta })-\log \hat{p}(\mathbf{\theta }),$$

$$p(\mathbf{Z}|\mathbf{F},\mathbf{L},\mathbf{\Lambda })p(\mathbf{F})$$

$$\times \exp \left\{-\frac{1}{2}\text{mbox}tr\left[\mathbf{F}{\mathbf{F}}^t{\left({\mathbf{\Psi }}_{\mathcal{F}}+{\kappa }_0^{-1}\mathbf{J}\right)}^{-1}\right]\right\}.$$

$$q^{\ast }\left({\mathbf{F}}_{i.}|\beta ,{\mathbf{e}}_i\mathbf{Z},{\mathbf{F}}_{-i.},\mathbf{L},\mathbf{\Lambda }\right)\propto l({\mathbf{e}}_i\mathbf{Z}|{\mathbf{F}}_{i.},\mathbf{L},\mathbf{\Lambda }{)}^{\beta }\hat{p}({\mathbf{F}}_{i.}|{\mathbf{F}}_{-i.}).$$

$$q^{\ast }({\mathbf{F}}_{i.}|\beta ,$$

$$\propto \exp \left\{-\frac{1}{2}\beta \left({\mathbf{e}}_i^t\mathbf{Z}-{\mathbf{F}}_{i\cdot }\mathbf{L}\right)\mathbf{\Lambda }{\left({\mathbf{e}}_i^t\mathbf{Z}-{\mathbf{F}}_{i\cdot }\mathbf{L}\right)}^t\right\}$$

$$\times \exp \left\{-\frac{1}{2}{\left({\mathbf{F}}_{i\cdot }^t-{\mathbf{\mu }}_{{\mathbf{F}}_{\text{-}\mathbf{i}\cdot }}\right)}^t{\mathbf{\Lambda }}_{{\mathbf{F}}_{\text{-}\mathbf{i}\cdot }}\left({\mathbf{F}}_{i\cdot }^t-{\mathbf{\mu }}_{{\mathbf{F}}_{\text{-}\mathbf{i}\cdot }}\right)\right\}$$

$$\propto \exp \{-\frac{1}{2}{\mathbf{F}}_{i.}\left(\beta \mathbf{L\Lambda }{\mathbf{L}}^t+{\mathbf{\Lambda }}_{{\mathbf{F}}_{\text{-}\mathbf{i}\cdot }}\right){\mathbf{F}}_{i\cdot }^t-2{\mathbf{F}}_{i\cdot }\left(\beta \mathbf{L\Lambda }{\mathbf{Z}}^t{\mathbf{e}}_i+{\mathbf{\Lambda }}_{{\mathbf{F}}_{\text{-}\mathbf{i}\cdot }}{\mathbf{\mu }}_{{\mathbf{F}}_{\text{-}\mathbf{i}\cdot }}\right)\}$$

$$\propto \exp \left\{-\frac{1}{2}{\left({\mathbf{F}}_{i.}^t-\mathbf{M}{\left(\beta \right)}_i^{\left(\mathbf{F}\right)}\right)}^t{\left(\mathbf{V}{\left(\beta \right)}_i^{\left(\mathbf{F}\right)}\right)}^{-1}\left({\mathbf{F}}_{i.}^t-\mathbf{M}{\left(\beta \right)}_i^{\left(\mathbf{F}\right)}\right)\right\},$$

Full text

Version dated:

Phylogenetic Factor Analysis

Max R. Tolkoff1, Michael E. Alfaro2, Guy Baele3, Philippe Lemey3,

and Marc A. Suchard1,4,5

1Department of Biostatistics, Jonathan and Karin Fielding School of Public Health, University of California, Los Angeles, United States

2*Department of Ecology and Evolutionary Biology, University of California, Los Angeles, United States

3Department of Microbiology and Immunology, Rega Institute, KU Leuven, Leuven, Belgium

4Department of Biomathematics, David Geffen School of Medicine at UCLA, University of California, Los Angeles, United States

5Department of Human Genetics, David Geffen School of Medicine at UCLA, Universtiy of California, Los Angeles, United States*

Corresponding author: Marc A. Suchard, Departments of Biostatistics, Biomathematics, and Human Genetics, University of California, Los Angeles, 695 Charles E. Young Dr., South, Los Angeles, CA 90095-7088, USA; E-mail: msuchard@ucla.edu

Abstract

Phylogenetic comparative methods explore the relationships between quantitative traits adjusting for shared evolutionary history. This adjustment often occurs through a Brownian diffusion process along the branches of the phylogeny that generates model residuals or the traits themselves. For high-dimensional traits, inferring all pair-wise correlations within the multivariate diffusion is limiting. To circumvent this problem, we propose phylogenetic factor analysis (PFA) that assumes a small unknown number of independent evolutionary factors arise along the phylogeny and these factors generate clusters of dependent traits. Set in a Bayesian framework, PFA provides measures of uncertainty on the factor number and groupings, combines both continuous and discrete traits, integrates over missing measurements and incorporates phylogenetic uncertainty with the help of molecular sequences. We develop Gibbs samplers based on dynamic programming to estimate the PFA posterior distribution, over three-fold faster than for multivariate diffusion and a further order-of-magnitude more efficiently in the presence of latent traits. We further propose a novel marginal likelihood estimator for previously impractical models with discrete data and find that PFA also provides a better fit than multivariate diffusion in evolutionary questions in columbine flower development, placental reproduction transitions and triggerfish fin morphometry.

(Keywords: Bayesian inference; comparative methods; morphometrics; phylogenetics)

1 Introduction

Phylogenetic comparative methods revolve around uncovering relationships between different characteristics or traits of a set of organisms over the course of their evolution. One way to gain insight into these interactions is to analyze unadjusted correlations between traits across taxa. However, as insightfully noted by Felsenstein (1985), unadjusted analyses introduce the inherent challenge that any association uncovered may reflect the shared evolutionary history of the organisms being studied, and hence their similar traits values, rather than processes driving traits to co-vary over time. Thus, studies to identify co-varying evolutionary trait processes must simultaneously adjust for shared evolutionary history.

There have been many attempts to accomplish this goal. Felsenstein (1985) and Ives and Garland Jr. (2010) are two such important examples, but they rely on a known evolutionary history described by a fixed phylogenetic tree and consider univariate evolutionary processes giving rise to only single traits. Felsenstein (1985) treats continuous traits as undergoing conditionally independent, Brownian diffusion down the branches of the phylogenetic tree and Ives and Garland Jr. (2010) posit a regression model where the tree determines the error structure in the univariate outcome model. Huelsenbeck and Rannala (2003) adapt the Brownian diffusion description in a Bayesian framework with the goal of drawing simultaneous inference on both the tree from molecular sequence data as well as the correlations of interest related to a small number of traits through a multivariate Brownian diffusion process. Lemey et al. (2010) extend the multivariate process by relaxing the strict Brownian assumption along distinct branches in the tree using a scale mixture of normals representation. Cybis et al. (2015) jointly model molecular sequence data and multiple traits using a multivariate latent liability formulation to combine both continuous and discrete observations and determine their correlation structure while adjusting for shared ancestry. This method is effective, but inference remains computationally expensive and estimates of the high-dimensional correlation matrix between latent traits is often difficult to interpret when addressing scientifically relevant questions. Additional frequentist methods include, Revell (2009) who use a phylogenetically adjusted principal components analysis, Adams (2014) who use a phylogenetic least squares analysis, and Clavel et al. (2015) who also use a multivariate diffusion method. All of these methods, however require large matrix inversions which make them ill suited to adaptations to full Bayesian inference, or bootstrapping to provide measures of uncertainty.

One way to alleviate these problems lies with dimension reduction through exploratory factor analysis (Aguilar and West 2000). Factor analysis is the inferred decomposition of observed data into two matrices, a factor matrix representing a set of underlying unobserved characteristics of the subject which give rise to the observed characteristics and a loadings matrix which explains the relationship between the unobserved and observed characteristics. Another form of dimension reduction through matrix decomposition is an eigen decomposition known as a principal components analysis (PCA). Santos (2009) provides a method for constructing PCA adjusted for evolutionary history. This method, however, has the same problems typically associated with PCA, namely that it is not invariant to the scaling of the data and the elimination of the smaller components necessitates some information loss. In a frequentist setting, the author also provides no approach for simultaneous inference on the phylogenetic tree that is rarely known without error (Huelsenbeck and Rannala 2003). In addition, there lacks a reasonable prescription for measuring uncertainty about which traits contribute to which principle components. Rai and Daume (2008) design a factor analysis method which uses a Kingman coalescent to construct a dendrogram across a factor analysis for genetic data. While this is similar to the idea we will employ, this specific method uses a dendrogram between, rather than within, factors and is thus ill suited to handle the important problem we tackle in this paper. Namely, researchers often seek to identify a small number of relatively independent evolutionary processes, each represented by a factor changing over the tree, that ultimately give rise to a large number of observed, dependent traits.

To formulate such a phylogenetic factor analysis (PFA) model, we begin with usual Bayesian factor analysis, as posited by Lopes and West (2004) and Quinn (2004), which represents underlying latent characteristics of a group of organisms through a factor matrix and maps those latent characteristics to observed characteristics via a loadings matrix. In a standard factor analysis, the underlying factors for each species would be assumed to be independent of each other, however this does nothing to adjust for evolutionary history. Vrancken et al. (2015) describe how a high-dimensional Brownian diffusion can be used to describe the relationship between all of these observed traits, however the signal strength of the results of analyzing this model can be quite poor. By using independent Brownian diffusion priors on our factors, our PFA model groups traits into a parsimonious number of factors while successfully adjusting for phylogeny. Scientifically, these diffusions represent independent evolutionary processes. We use Markov chain Monte Carlo (MCMC) integration in order to draw inference on our model through a Metropolis-within-Gibbs approach. This facilitates both a latent data representation (Cybis et al. 2015) for integrating discrete and continuous traits and a natural method to handle missing data relevant to our problems. We further rely on path sampling methods (Gelman and Meng 1998) to determine the appropriate number of factors (Ghosh and Dunson 2009). Since the latent, probit model necessitates the use of hard thresholds, we now have introduced an inherent difficulty in path sampling. In order to get around this difficulty, we employ a novel method which relies on softening the threshold necessitated by the probit model slowly over the course of the path. We additionally develop a novel method by which to handle identifiability issues inherent to factor analysis by taking advantage of the fact that correlated elements in the loadings matrix tend to be correlated across the MCMC chain.

We show that our PFA method performs superiorly to a high-dimensional Brownian diffusion in both signal strength and, when we are inferring large numbers of latent traits, speed using the examples of the evolution of the flower genus Aquilegia, as well as the reproduction of the fish family Poeciliidae that involves trait measurements missing at random. Lastly, we explore the dorsal, anal and pectoral fin shapes of the fish family Balistidae in order to explore this method’s ability to handle situations where the number of traits are large compared to the number of species and to explore the simultaneous inference on our method along with the evolutionary history of these organisms with the aid of sequence data. The PFA model and its inference tools will be released in the popular phylogenetic inference package BEAST (Drummond et al. 2012).

2 Methods

2.1 Phenotypic Trait Evolution

Consider a collection of $N$ biological entities (taxa). From each taxon $i=1,\ldots,N$, we observe a $P$-dimensional measurement $\mathbf{Y}_{i}=(Y_{i1},\ldots,Y_{iP})$ of traits and, if available, a molecular sequence $\mathbf{S}_{i}$. We organize these phenotypic traits into an $N\times P$ matrix $\mathbf{Y}=(\mathbf{Y}_{1},\ldots,\mathbf{Y}_{N})^{\prime}$ and an aligned sequence matrix $\mathbf{S}$. These taxa are related to each other through an evolutionary history ${\cal{F}}$, informed through $\mathbf{S}$, and we are interested in learning about the evolutionary processes along this history that give rise to observed traits $\mathbf{Y}$.

The history ${\cal{F}}$ consists of a tree topology $\tau$ and a series of branch lengths $\mathbf{B}$. The tree topology is a bifurcating directed acyclic graph with a single generating point called the root, representing the most recent common ancestor of the given taxa, and with end points, each of which corresponds to a different taxon. The branch lengths correspond to edge weights of the graph, reflecting the evolutionary time before bifurcations. The history ${\cal{F}}$ may be known and fixed, or unknown and jointly inferred using $\mathbf{Y}$ and $\mathbf{S}$. For further details on constructing the sequence-informed prior distribution $p({\cal{F}}\,|\,\mathbf{S})$ and integrating over ${\cal{F}}$ when unknown, see, e.g., Suchard et al. (2001) or Drummond et al. (2012).

In order to simultaneously model continuous, binary and ordinal traits, we adapt a latent data representation through the partially observed, standardized matrix $\mathbf{Z}$ with entries

[TABLE]

where $\hat{Y}_{j}$ is the mean of trait $j$ across taxa, $\hat{\sigma}_{j}$ is its standard deviation for $j=1,\ldots,P$ and, more importantly, $Z_{ij}\in\mathbb{R}$ is an unknown random variable that satisfies the restrictions

[TABLE]

and $c\in\{1,\ldots,m_{j}\}$ for $m_{j}$-valued binary/ordinal data for trait $j$.

For identifiability, latent trait cut-points $\boldsymbol{\gamma}_{j}=(\gamma_{j0},\dots,\gamma_{jm_{j}})$ take on the restrictions $\gamma_{j0}=-\infty$, $\gamma_{j1}=0$ and $\gamma_{jm_{j}}=\infty$ or are otherwise random and jointly inferred. Grouping cut-points for all binary or ordinal traits into $\boldsymbol{\gamma}$, Cybis et al. (2015) suggest assuming that differences between the small number of successive, random cut-points are a priori exponentially distributed with mean $\frac{1}{2}$ to define their density $p(\boldsymbol{\gamma})$. Cybis et al. (2015) also discuss in detail how to treat categorical data in this sort of analysis. Since we do not use examples which contain non-ordered categorical data we elect not to describe those methods in these sections, but we will mention that they are implemented in BEAST and are easily adapted to fit the methods described in this paper.

In order to uncover the biological relationships amongst traits in $\mathbf{Z}$ while controlling for evolutionary history, previous work relies on a Gaussian process generative model induced through considering conditionally independent Brownian diffusion along each branch in ${\cal{F}}$ (Felsenstein 1985). In a multivariate setting, a $P\times P$ variance matrix $\boldsymbol{\Sigma}$ and unobserved, $P$-dimensional root trait value $\boldsymbol{\mu}_{R}$ characterize the process. Pybus et al. (2012) identify that analytic integration of $\boldsymbol{\mu}_{R}$ is possible by assuming that $\boldsymbol{\mu}_{R}$ is a priori multivariate normally distributed with a fixed hyperprior mean $\boldsymbol{\mu}_{0}$ and variance equal to $\kappa_{0}^{-1}\boldsymbol{\Sigma}$, where $\kappa_{0}$ is a fixed hyperprior sample-size. Consequentially, given ${\cal{F}}$ and $\boldsymbol{\Sigma}$, the latent traits $\mathbf{Z}$ are distributed according to a matrix-normal (MN)

[TABLE]

where $\boldsymbol{\Psi}_{{\cal{F}}}+\kappa_{0}^{-1}\mathbf{J}$ is the across-taxa (row) variance and a deterministic function of phylogeny ${\cal{F}}$, $\boldsymbol{\Sigma}$ is the across-trait (column) variance, and $\mathbf{J}$ is a $N\times N$ matrix of ones (Vrancken et al. 2015). Traits $\mathbf{Z}$ have density function

[TABLE]

where $\mbox{tr}\left[\cdot\right]$ is the trace operator and $\mathbf{1}$ is a $N$-dimensional column vector of ones. Tree variance matrix $\boldsymbol{\Psi}_{{\cal{F}}}$ contains diagonal elements that are equal to the sum of the adjusted branch lengths in ${\cal{F}}$ between the root node and taxon $i$, and off-diagonal elements $(i,i^{\prime})$ that are equal to the sum of the adjusted branch lengths between the root node and the most recent common ancestor of taxa $i$ and $i^{\prime}$, where the adjusted branch lengths represent a function of wall time and a branch rate accounting for variation in evolution rate over the course of the tree. For our diffusion model, we scale our tree such that from the root to the most recent tip we say that the process has undergone one diffusion unit.

Placing a conjugate prior distribution on $\boldsymbol{\Sigma}$, such as $\boldsymbol{\Sigma}^{-1}\sim\text{Wishart}_{\nu}(\boldsymbol{\Lambda}_{R_{0}})$ where $\nu$ is the hyperprior degrees of freedom and $\boldsymbol{\Lambda}_{R_{0}}$ is the hyperprior belief on the structure of the inverse of the variance matrix $\boldsymbol{\Sigma}$, enables inference about its posterior distribution, shedding light on how the evolution of these traits relate to each other. Such inference often requires repeated evaluation of density (4), especially when the phylogeny ${\cal{F}}$ or variance $\boldsymbol{\Sigma}$ is random. This evaluation suggests a computational order ${\cal O}\left(N^{3}+P^{3}\right)$, arising from the inversion of the $N\times N$ variance matrix $\boldsymbol{\Psi}_{{\cal{F}}}+\kappa_{0}^{-1}\mathbf{J}$ and $P\times P$ variance matrix $\boldsymbol{\Sigma}$. One easily avoids the latter by parameterizing the model in terms of $\boldsymbol{\Sigma}^{-1}$ (Lemey et al. 2010). To address the former, Pybus et al. (2012) provide an ${\cal O}\left(NP^{2}\right)$ dynamic programming algorithm to evaluate (4) without inversion of the across-taxa variance matrix, similar to Freckleton (2012). This advance certainly makes for more tractable inference under these diffusion models as $N$ grows large, but the quadratic dependence on $P$ still hampers their use for high-dimensional traits. Inference can often be slow, taking as long as a day for problems with a dozen traits and about 30 taxa to mix properly (Cybis et al. 2015). Finally, direct inference on $\boldsymbol{\Sigma}$ can often fail to produce a coherent and interpretable conclusion about the number of independent evolutionary processes generating the traits if the matrix cannot be reordered to form approximately separated blocks especially if the signal is too weak to produce many statistically significant cells.

2.2 Factor Analysis

To infer potentially low dimensional evolutionary structure among traits, we rely on dimension reduction via a phylogenetic factor analysis (PFA). This model builds on the premise that a small, but unknown number $K\ll\text{min}\left(N,P\right)$ of a priori independent univariate Brownian diffusion processes along ${\cal{F}}$ provides a more parsimonious description of the covariation in $\mathbf{Z}$ than a $P$-dimensional multivariate diffusion. We parameterize the PFA in terms of an $N\times K$ factor matrix $\mathbf{F}=(\mathbf{F}_{1},\ldots,\mathbf{F}_{K})$ whose $K$ columns $\mathbf{F}_{k}=\left(F_{1k},\ldots,F_{Nk}\right)^{t}$ for $k=1,\ldots,K$ represent the unobserved independent realizations of univariate diffusion at each of the $N$ tips in ${\cal{F}}$, a $K\times P$ loadings matrix $\mathbf{L}=\{L_{kj}\}$ that relates the independent factor columns to $\mathbf{Z}$, and an $N\times P$ model error matrix $\boldsymbol{\epsilon}$, such that

[TABLE]

To inject information about and control for shared evolutionary history ${\cal{F}}$, we specify that

[TABLE]

where $\mathbf{I}_{\left(\cdot\right)}$ is the identity matrix of appropriate dimension and the residual column precision $\mathbf{\Lambda}$ is a diagonal matrix with entries $(\Lambda_{1},\ldots,\Lambda_{P})$. Lastly, since $K$ is unknown, we place a reasonably conservative zero-truncated-Poisson prior on it, such that $p(K=1)=1/2$.

To better appreciate the details of the PFA model, we briefly compare it to a typical Bayesian factor analysis. Typical factor analyses assume that all entries of $\mathbf{F}$ are independent and identically distributed (iid) as $\text{N}\left(0,1\right)$, normal random variables with mean [math] and variance $1$. In PFA, the shared evolutionary history ${\cal{F}}$ specifies the correlation structure within the $N$ entries of column $\mathbf{F}_{k}$. Often, one refers to a given column as a “factor.” Across factors, the column variance remains $\mathbf{I}_{K}$ to reflect our assertion that the underlying evolutionary processes generating $\mathbf{F}_{k}$ are independent of each other. Note that in this model the number of parameters undergoing Brownian Diffusion is assumed to be of dimension $K$ as opposed to of dimension $P$ in the previous model.

To complete model specification of the loadings $\mathbf{L}$ and residual error $\boldsymbol{\epsilon}$, we assume

[TABLE]

otherwise $\Lambda_{j}=1$ to preserve identifiability under the scale-free latent model for discrete traits. Here, $\Gamma\left(\alpha_{\mathbf{\Lambda}},\beta_{\mathbf{\Lambda}}\right)$ signifies a gamma distributed random variable with hyperparameter scale $\alpha_{\mathbf{\Lambda}}$ and rate $\beta_{\mathbf{\Lambda}}$.

Without further restrictions on $\mathbf{L}$, any factor analysis remains over-specified.

For example, given an orthogonal $K\times K$ matrix $\mathbf{T}$, one may rotate $\mathbf{F}$ in one direction and $\mathbf{L}$ in the other and arrive at the same data likelihood, since $\mathbf{F}\mathbf{L}=\mathbf{F}\mathbf{T}\mathbf{T}^{t}\mathbf{L}$. To address this identifiability issue, we fix lower triangular entries $L_{kj}=0$ for $k>j$ (Geweke and Zhou 1996; Aguilar and West 2000). It is also standard practice to apply the restriction $L_{kk}>0$, since otherwise $\mathbf{F}\mathbf{L}=(\text{-}\mathbf{F})(\text{-}\mathbf{L})$. While the constraint yields an identifiable posterior distribution with respect to $\mathbf{F}$ and $\mathbf{L}$, we do not pursue it here because it introduces bias into our scientific inference on $\mathbf{L}$ and, instead, search for an alternative.

The diagonal and upper triangular entries $L_{kj}$ for $k\leq j$ of the loadings $\mathbf{L}$ inform the magnitude and effect-direction that the evolutionary process captured in factor $\mathbf{F}_{k}$ contributes to trait $j$. It is possible, and we would argue likely, that $\mathbf{F}_{k}$ has little or no influence on the trait arbitrarily labeled $k$, such that most of the posterior mass of $L_{kk}$ lies around and close to [math]. Artificially restricting $L_{kk}>0$ forces all of this mass above [math], signifying a positive association with prior, and hence posterior, probability $1$.

To combat this bias, we recouch these identifiability conditions as a label switching problem in a mixture model and propose a post hoc relabeling algorithm (Stephens 2000). We require $K$ sign constraints, one for each column-row outer-product in forming $\mathbf{F}\mathbf{L}$, for posterior identification. In our prior, we modify Equation (7) to further assign one non-zero entry $L_{kj}>0$ per row, but do not specify which one; this assignment mirrors the mixture model labeling. Hence, we allow the data, not an arbitrary decision, to determine which entry per row reflects a positive association with probability $1$, decreasing potential bias.

Recalling that continuous traits are standardized in $\mathbf{Z}$ to have mean [math] and variance $1$ affords several benefits. First, we can posit a $\mathbf{0}$-matrix mean for $\mathbf{F}$ in Equation (6) without loss of information. But, more importantly, when we draw inference on $\mathbf{\Lambda}$, we can interpret traits which have precision elements that demonstrate considerable posterior mass at or below 1 to be described insufficiently by the model, since the factors provide no insight beyond a random normal model. A third advantage is that standardization helps us select reasonable scales for the non-zero entries in $\mathbf{L}$, namely that these have variance $1$, and hyperparameters for $\mathbf{\Lambda}$, specifically that $\frac{\alpha_{\mathbf{\Lambda}}}{\beta_{\mathbf{\Lambda}}}=1$. In practice, $\alpha_{\mathbf{\Lambda}}=\frac{1}{3}$ and $\beta_{\mathbf{\Lambda}}=\frac{1}{3}$ for analyses in this paper. While these hyperparameter choices are by no means perfect we feel that, under the paradigm of data scaling, they are reasonable and generalizable across a variety of problems.

This model is a simplified form of the item factor analysis models that are described by Quinn (2004) in the political science literature and Beguin and Glas (2001) in the psychology literature with a tree as a prior on the factors instead of an independent normal distribution. In fact, the methods for treating binary and ordinal data described in Quinn (2004) are the same as those described in Cybis et al. (2015), making for a convenient adaptation of this factor analysis model to phylogenetics using existing software in BEAST.

2.3 Inference

Given the trait measurements $\mathbf{Y}$ and aligned sequences $\mathbf{S}$, we strive to learn about the joint posterior distribution of the number of evolutionary processes $K$, factors $\mathbf{F}$, loadings $\mathbf{L}$, column precisions $\mathbf{\Lambda}$, latent trait cut-points $\boldsymbol{\gamma}$ and evolutionary history ${\cal{F}}$

[TABLE]

where $p(\mathbf{Y}\,|\,\mathbf{Z},\boldsymbol{\gamma})\propto\mathbf{1}(\mathbf{Y}\,|\,\mathbf{Z},\boldsymbol{\gamma})$ is the indicator function that the restrictions in Equation (2) hold. We accomplish this inference through MCMC, using a random-scan Metropolis-within-Gibbs scheme (Liu et al. 1995) for fixed $K$ and a modification of path sampling to then estimate the marginal posterior $p(K\,|\,\mathbf{Y},\mathbf{S})$. For fixed $K$, our Metropolis-within-Gibbs scheme employs transition kernels described in Cybis et al. (2015) and references therein to integrate over the evolutionary history ${\cal{F}}$ and unobserved, latent traits $Z_{ij}$ and cut-points $\boldsymbol{\gamma}_{j}$ where trait $j$ is discrete.

Here, we focus on transition kernels within the scheme to integrate over the factors $\mathbf{F}$, loadings $\mathbf{L}$ and residual column precision $\mathbf{\Lambda}$. Lopes and West (2004) derive full conditional distributions for the columns of $\mathbf{L}$ and diagonals of $\mathbf{\Lambda}$ under a traditional factor analysis. These full conditional distributions do not change under a PFA and we use them for Gibbs sampling. Specifically, for column $j$ of $\mathbf{L}$, the first $k^{\prime}=\text{min}\left(j,K\right)$ entries are non-zero and, given all other random variables, distributed according to a multivariate normal (MVN)

[TABLE]

parameterized in terms of its mean

[TABLE]

where $\mathbf{F}_{1:k^{\prime}}=\left(\mathbf{F}_{1},\ldots,\mathbf{F}_{k^{\prime}}\right)$ is the first $k^{\prime}$ columns of $\mathbf{F}$ and $\mathbf{e}_{j}$ is the unit-vector in the direction of trait $j$. Further,

[TABLE]

if trait $j$ is continuous. Appendix A provides derivations of these full conditional distributions. Gibbs sampling all columns of $\mathbf{L}$ carries a computation order ${\cal O}\left(NK^{2}P\right)$, arising from the matrix multiplication of $\mathbf{F}_{1:k^{\prime}}^{t}\mathbf{F}_{1:k^{\prime}}$ for each trait. The matrix inversion is not rate-limiting here since $N\gg K$. Likewise, Gibbs sampling $\mathbf{\Lambda}$ remains very light-weight at ${\cal O}\left(NKP\right)$, stemming from the sparse multiplication of $\mathbf{F}\mathbf{L}\mathbf{e}_{j}$ for each trait. While we write that the order of both Gibbs samplers depend on $P$ to be clear that we must iterate over all traits, the astute reader has already recognized the conditional independence of updates between traits, such that we may execute updates for each trait in parallel.

The traditional Gibbs sampler for $\mathbf{F}$ fails in the phylogenetic setting for more than a handful of taxa, since determining the full conditional distribution of $\mathbf{F}$ requires inverting the matrix $\left(\boldsymbol{\Psi}_{{\cal{F}}}+\kappa_{0}^{-1}\mathbf{J}\right)$. As mentioned previously, but worth repeating, this task stands as prohibitive with a computational order ${\cal O}\left(N^{3}\right)$ and presents a major challenge for PFA.

We circumvent this difficulty by exploiting the structure of the phylogenetic tree ${\cal{F}}$. Probability models on directed, acyclic graphs lend themselves well to dynamic programming for determining marginalized data likelihoods, such as Felsenstein’s pruning algorithm for sequence data (Felsenstein 1973) and related work for Brownian diffusion (Pybus et al. 2012), and conditional predictive distributions, like those obtained for (ancestral) sequence reconstruction.

In extending these conditional distributions to Brownian diffusion, first let $\mathbf{F}_{i\cdot}=(F_{i1},\ldots,F_{iK})$ identify row $i$ of $\mathbf{F}$, more specifically all latent factor values attributed to taxon $i$, and let $\mathbf{F}_{\text{-}i\cdot}$ concatenate the remaining rows. Given that $\mathbf{F}$ is matrix-normally distributed with an across-taxa (row) variance that depends on the phylogeny ${\cal{F}}$, Cybis et al. (2015) provide a tree-traversal-based algorithm to determine $p(\mathbf{F}_{i\cdot}\,|\,\mathbf{F}_{\text{-}i\cdot},{\cal{F}})$ that remains a multivariate normal distribution. The algorithm requires first a post-order tree-traversal to determine the joint distribution of all tip-values descendent to each internal node and then a pre-order tree-traversal back to taxon $i$ to compute its prior conditional mean $\boldsymbol{\mu_{\mathbf{F}_{\text{-}i\cdot}}}$ and precision $\boldsymbol{\Lambda_{\mathbf{F}_{\text{-}i\cdot}}}$. Since the across-factor (column) variance on $\mathbf{F}$ is diagonal, the dynamic programming algorithm runs quickly in ${\cal O}\left(NK\right)$. Using this result, we determine the full conditional distribution

[TABLE]

with mean

[TABLE]

where $\mathbf{e}_{i}$ is the unit-vector in the direction of taxon $i$. Appendix A delivers a derivation of this full conditional distribution. The evaluation of this full conditional distribution runs in ${\cal O}\left(K^{2}P\right),$ where the term $\mathbf{L}\mathbf{\Lambda}\mathbf{L}^{t}$ is rate limiting.

Employing Equations (13) - (15), we can cycle over $i$ to fabricate a tractable Gibbs sampler for $\mathbf{F}$ with total computational order ${\cal O}\left(N^{2}K+NK^{2}P\right)$. It is fruitful to compare this work with the rate-limiting step for inference under the non-sparse model. Here, sampling the precision matrix $\boldsymbol{\Sigma}^{-1}$ carries a computational cost of ${\cal O}\left(NP^{2}\right)$. From these bounds, it is clear that increasing numbers of taxa $N$ should limit PFA, while increasing numbers of traits $P$ should limit the non-sparse model from a computational work per MCMC iteration perspective. However, per-iterative arguments ignore the posterior correlation between model parameters and its influence on MCMC mixing times.

Finally, to maintain identifiability with respect to $\mathbf{F}$ and $\mathbf{L}$ in the posterior, we propose a simple post hoc relabeling algorithm (Stephens 2000). We sample $(\mathbf{F}^{(m)},\mathbf{L}^{(m)})$ from $p(K,\mathbf{F},\mathbf{L},\mathbf{\Lambda},\boldsymbol{\gamma},{\cal{F}}\,|\,\mathbf{Y},\mathbf{S})$ for MCMC iteration $m=1,\ldots,M$ assuming a sign-unconstrained prior. From this unconstrained sample, we select for each row $k$ in $\mathbf{L}$ the column element with the fewest number of sign changes between iterations. Assume for row $k$, this is column $j_{k}$. We then constrain our sample by multiplying $\mathbf{F}_{k}^{(m)}$ and row $k$ of $\mathbf{L}^{(m)}$ by the sign of $L_{kj_{k}}^{(m)}$. No further sample reweighing is necessary because $p(\mathbf{F}\,|\,K,{\cal{F}})=p(\text{-}\mathbf{F}\,|\,K,{\cal{F}})$ is also invariant to reflection.

Model Selection

To estimate the marginal posterior density $p(K\,|\,\mathbf{Y},\mathbf{S})$, we rely on a variant of path sampling that we equip to successfully integrate latent variable $\mathbf{Z}$ when traits are discrete. We employ our variant to approximate each marginal likelihood $p(\mathbf{Y},\mathbf{S}\,|\,K=k)$ for $k=1,\ldots,S$, where $S$ is a relatively small number such as $\min\{P,10\}$, after which we approximate $p(\mathbf{Y},\mathbf{S}\,|\,K>S)=0$. Then, invoking Bayes theorem, $p(K=k\,|\,\mathbf{Y},\mathbf{S})\propto p(\mathbf{Y},\mathbf{S}\,|\,K=k)p(K=k)$. Moreover, through this approach, we can address the model selection problem of how many independent factors do the data support through Bayes factors (Jeffreys 1935):

[TABLE]

Lopes and West (2004) and Ghosh and Dunson (2009) have been strong proponents of Bayes factors to determine the optimal number of factors in a traditional factor analysis, where Lopes and West (2004) employ a simple harmonic mean estimator (Newton and Raftery 1994) to estimate their marginal likelihoods. This estimator performs poorly in highly structured phylogenetic models and path sampling has largely supplanted it (Baele et al. 2012).

Path sampling is an MCMC-based integration technique to estimate marginal likelihoods, such as $p(\mathbf{Y},\mathbf{S}\,|\,K)$. The technique constructs a series of power posteriors (Friel and Pettitt 2008) at various temperatures $\beta\in[0,1]$, where $\beta=1$ corresponds to a joint density $l(\mathbf{Y},\mathbf{S},\mathbf{Z},\mathbf{F},\mathbf{L},\mathbf{\Lambda},\boldsymbol{\gamma}\,|\,K)$ proportional, but with an unknown constant, to $p(\mathbf{Y},\mathbf{S}\,|\,K)$ and $\beta=0$ yields a normalized density $\hat{p}(\mathbf{Z},\mathbf{F},\mathbf{L},\mathbf{\Lambda},{\cal{F}},\boldsymbol{\gamma}\,|\,K)$ that does not depend on the data, often a combination of the prior and other working distributions, see e.g. (Baele et al. 2016). The usual power posterior path is $q(\beta,\mathbf{Y},\mathbf{S},\boldsymbol{\theta})=l(\mathbf{Y},\mathbf{S},\boldsymbol{\theta})^{\beta}\times\hat{p}(\boldsymbol{\theta})^{1-\beta}$, where $\boldsymbol{\theta}$ is the set of all parameters in the model we are considering. For example, in PFA, $\boldsymbol{\theta}=\left\{\mathbf{Z},\mathbf{F},\mathbf{L},\mathbf{\Lambda},{\cal{F}},\boldsymbol{\gamma}\right\}.$

In latent models with discrete traits, however, the support of the latent variable $\mathbf{Z}$ changes when the data are observed (Heaps et al. 2014). In particular, our unnormalized joint density $l(\mathbf{Y},\mathbf{S},\boldsymbol{\theta})$ is zero for values of $\mathbf{Z}$ that are incompatible with $\mathbf{Y}$ because $p(\mathbf{Y}\,|\,\mathbf{Z},\boldsymbol{\gamma})=0$, therefore a trait $Z_{ij}$ only has support over $(\gamma_{i(c-1)},\gamma_{ic}]$ if $Y_{ij}=c$, while $\hat{p}(\cdot)$ places non-zero density over all possible values $Z_{ij}\in(-\infty,\infty)$. Our working distribution, for example, assumes $Z_{ij}\sim N(0,1)$ when $Z_{ij}$ is random. If we factor $l(\mathbf{Y},\mathbf{S},\boldsymbol{\theta})$ into a support condition $\mathbf{1}(\mathbf{Y}\,|\,\mathbf{Z},\boldsymbol{\gamma})$ and the remaining likelihood $h(\mathbf{Y},\mathbf{S},\boldsymbol{\theta})$, then the standard path used in this scenario (Heaps et al. 2014) is

[TABLE]

For the power posterior method to yield the marginal likelihood $p(\mathbf{Y}\,|\,K),$ it is necessary (Friel and Pettitt 2008) that

[TABLE]

Plugging (17) into (18), we find

[TABLE]

If we define $\boldsymbol{\Omega}$ as the region where $\mathbf{1}(\mathbf{Y}\,|\,\mathbf{Z},\boldsymbol{\gamma})=1,$ then we see that

[TABLE]

since $\boldsymbol{\Omega}\subsetneq$ the support of $\boldsymbol{\theta}.$ While it is theoretically possible to construct $\hat{p}(\boldsymbol{\theta})$ such that it is normalized to 1 over $\boldsymbol{\Omega}$, previous attempts to do so have failed. Alternatively, Heaps et al. (2014) attempt to approximate such a distribution by fixing $\boldsymbol{\gamma}$ and ignoring the corresponding integral.

We posit an exact solution by proposing a new path that relies on a softening threshold. Consider the modified path

[TABLE]

Following from (18), we find that

[TABLE]

by construction.

Lastly, in order to adapt the power posterior method, at each step in the series we need to compute the derivative of $\log q^{*}\left(\beta,\mathbf{Y},\mathbf{S},\boldsymbol{\theta}\right)$ with respect to $\beta$. From Equation (21), we see that

[TABLE]

and observe that there is no singularity at $\beta=1$ since, at that point in the path, latent variable $\mathbf{Z}$ only assumes values in $\boldsymbol{\Omega}$, such that $\mathbf{1}(\mathbf{Y}\,|\,\mathbf{Z},\boldsymbol{\gamma})=1$.

3 Empirical Examples

3.1 Columbine Flower Development

Columbine genus Aquilegia flowers have attracted at least three different pollinators across their evolutionary history: bumblebees (Bb), hawkmoths (Hm) and hummingbirds (Hb). Whittall and Hodges (2007) question the role that these pollinators play in the tempo of columbine flower evolution, tracked through the color, length and orientation of different anatomical floral features, and are particularly interested in how transitions between pollinators relate to spur length. Cybis et al. (2015) take up this question by examining $P=12$ different traits for $N=30$ monophyletic populations from the genus Aquilegia that include $10$ continuously valued traits, a binary trait that indicates presence or absence of anthocyanin pigment and a final ordinal trait indicating the primary pollinator for that population. Whittall and Hodges (2007) propose a Bb-Hm-Hb ordering and we use the fixed phylogenetic tree the authors employ in their analysis. Through fitting a latent multivariate Brownian diffusion (LMBD) model parameterized in terms of a $12\times 12$ variance matrix $\boldsymbol{\Sigma}$, Cybis et al. (2015) find the data strongly support the proposed ordering over alternative orderings. We return to the relationship between pollinator and the other traits and test whether a PFA returns a better understanding of the evolutionary factors driving their interrelated change compared to an LMBD model.

Under our PFA, the most probable number of independent evolutionary processes is $K=2$, with a log Bayes factor $>7$ over the neighboring $K=1$ or $K=3$ factor parameterizations (Table 1). Further, the PFA with $K=2$ is favored over the LMBD model with a log Bayes factor $>24$ when assuming equal prior probabilities over these two models.

The PFA has high explanatory power for all continuous traits (Table 2) and Figure 1 presents our inference on the relationships between traits under the PFA with $K=2$ and compares these findings to inference under the LMBD model. The first evolutionary process $\mathbf{F}_{1}$ approximately partitions the traits into two groups. One group includes: orientation, blade brightness, spur brightness, sepal length, blade length, pollinator type, spur hue, spur length, blade hue, and expected trait values increase (displayed loadings entries $L_{kj}$ in purple) as the factor grows over the phylogeny. The other group includes: blade chroma, anthocyanins pigment presence and, with less posterior probability, spur chroma, and expected trait values decrease (green) as the factor grows. A possible exception to the $\mathbf{F}_{1}$ partitioning is the pollinator trait, where we estimate only a 92% posterior probability that this cell has the same sign as its posterior mean.

Ignoring the uncertainty in pollinator trait inclusion for the moment, this partitioning recapitulates the block structure that Cybis et al. (2015) report using an LMBD model and an arbitrary thresholding on the posterior mean estimates of the individual pairwise correlation entries in $\boldsymbol{\Sigma}$. However, in Figure 1 we quantify the LMBD uncertainty by shading our inference using the same probability measure as we do for our PFA model. Taking correlation uncertainty into consideration we see that, for example the LMBD model would assert that there is no correlation between blade chroma and spur hue. The PFA model by contrast offers the more nuanced assessment that these traits are related through two independent underlying processes, one process of which has a positive association between these traits, the other of which has a negative association.

In addition to improved uncertainty quantification in the block structure of traits, our PFA returns a second independent evolutionary process $\mathbf{F}_{2}$ that relates pollinator with spur length and, in addition, spur and blade chroma and hue, with posterior probability approaching $1$. The existence of two distinct processes, one of which directly connects pollinator and spur length, sheds additional insight into the original hypothesis that Whittall and Hodges (2007) pose. The LMBD model fails to pick up on this, in addition to returning a worse fit to the data.

3.2 Transitions to Placental Reproduction

The freshwater fish Poeciliidae represent a family of model organisms in which one can study the transition from non-placental to placental reproduction and the evolutionary pressures associated with placental introduction. Pollux et al. (2014) define a matrotrophy index to be the log-ratio of the dry weight of newborn fish to the dry weight of eggs at fertilization as a proxy measure of how reliant a fish species is on its placenta for reproduction. Using phylogenetic generalized least squares (PGLS) (Ives and Garland Jr. 2010), Pollux et al. (2014) find that Poeciliidae dichromatism, courtship behavior, superfetation, and a sexual selection index are all correlated over evolutionary history with the matrotrophy index. Unlike PFA, PGLS as used by Pollux et al. (2014) does not adjust for potential evolutionary relationships between the traits. Failure to do so can lead to false positive measures of association between individual traits and the matrotrophy index.

Pollux et al. (2014) collect from the literature or measure 14 life-history traits and compile from GenBank or sequence 28 different genes across Poeciliidae species. In our analysis, we only use $P=11$ traits since three of the original traits are functions of the included ones. Of these traits, five are discrete-valued: dimorphic coloration (dichromatism), courtship behavior, superfetation, the presence or absence of ornamental display traits and a count composite of the presence or absence of three other male behaviors (sexual selection index). Six are continuous-valued: log weight and log length for males and females, gonopodium length and matrotrophy index. Considering species with at least one trait measurement, there are $N=98$ taxa, for which we assume the same fixed phylogenetic tree that Pollux et al. (2014) estimate and similarly condition on in their PGLS analysis. Importantly, 182 trait measurements remain missing. We treat these measurements as missing-at-random in our PFA and do not need to further prune the tree or impute values that may further introduce bias.

Pollux et al. (2014) find that dichromatism, courtship behavior, superfetation, and sexual selection index are all correlated with the matrotrophy index. Figure 2 shows that this concurs with the results of a $K=2$ factor PFA. This small model fit also highlights a weakness of traditional factor analysis assumptions that fix the diagonal elements of the loadings matrix to be positive. In particular, dichromatism is unrelated to the other traits in the second factor, while the positivity constraint would have forced its inclusion. However, the most probable number of independent evolutionary processes is $K=3$ or $K=4$, with a log Bayes factor in favor $K=3$ over $K=2$ of $35.3$ and a log Bayes factor in favor of $K=4$ over $K=5$ of $4.9$ (Table 1). Since a log Bayes factor of only $0.3$ separates the $K=3$ and $K=4$ models, we include both models in our results, and the data strongly support these PFA models over the LMBD model (log Bayes factor $\approx$ 92).

Loadings for the independent evolutionary process factors $\mathbf{F}^{\left(3\right)}_{k}$ and $\mathbf{F}^{\left(4\right)}_{k}$ under the $K=3$ and $K=4$ PFA models, respectively, recapitulate a negative assocation between the matrotrophy index and dichromatism, courtship behavior, and sexual selection index, and a positive association with superfetation (Figure 2, first loading). However, unlike in Pollux et al. (2014), the PFA does not recover with high posterior probability a relationship between matrotrophy index and gonopodium length nor with body weights and lengths, suggesting that these were false positive findings. For both PFA models, second independent processes $\mathbf{F}^{\left(3\right)}_{2}$ and $\mathbf{F}^{\left(4\right)}_{2}$ drive dichromatism, courtship behavior, ornamental display traits and sexual selection index positively and superfetation and gonopodium length negatively. Both models also identify similar third independent processes $\mathbf{F}^{\left(3\right)}_{3}$ and $\mathbf{F}^{\left(4\right)}_{3}$ relating body lengths and weights. It is perhaps surprising that these size measurements are unrelated to any of the other reproductive characteristics. The only marked difference between the $K=3$ and $K=4$ factor models exists in the presence of a fourth evolutionary process $\mathbf{F}^{\left(4\right)}_{4}$ in the $K=4$ factor model that controls the presence or absence of superfetation independently of all other traits.

The precision elements $\boldsymbol{\Lambda}$ for both the $K=3$ and $K=4$ factor models are all significantly greater than 1 and therefore indicate that, for both models, our PFA provides good insight into the relationship of the continuous traits (Table 2). Further, the precision elements are in broad agreement between the $K=3$ and $K=4$ factor models, as we expect due to the negligible difference in marginal likelihoods.

Frequentist-based factor analysis is only identifiable if the number of parameters inferred for a variance/covariance matrix is greater than the number of parameters that need to be inferred for the factor analysis. Interestingly, our PFA model produces interpretable results in spite of the fact that the correlation model has 66 free parameters as opposed to 333 free parameters for the $K=3$ factor model, and 436 free parameters for the $K=4$ factor model.

3.3 Triggerfish Fin Shape

The fish family Ballistidae, commonly know as triggerfish, live mostly in reefs; however, the particular part of the reef in which they live can vary. This variability affects not only their diet, but also their mobility needs that fin shapes well reflect (Dornburg et al. 2011). To model shape changes through evolution, phylogenetic morphometrics often relies heavily on principle components analysis (PCA) (Revell 2009; Polly et al. 2013). However, deterministic data reduction via PCA can introduce bias (Uyeda et al. 2015) and, more importantly, inference of principal components while simultaneously adjusting for an uncertain evolutionary history remains a continuing challenge. PFA offers an alternative approach.

For $N=24$ triggerfish species, Dornburg et al. (2011) sequence and align 12S (833 nucleotides, nt) and 16S (563 nt) mitochondrial genes and RAG1 (1471 nt), rhodopsin (564 nt) and Tmo4C4 (575 nt) nuclear genes, and Dornburg et al. (2008) digitally photograph and mark 13 semi-landmark Cartesian coordinates for pectoral, dorsal and anal fins, generating $P=78$ measurements per species. Among these morphometric measurements, the species Balistapus undulatus is missing dorsal and anal fins landmarks, and the species Rhinecanthus assasi lacks pectoral fin landmarks. For these, we assume the missing data are missing at random.

To accommodate phylogenetic uncertainty within $p({\cal{F}}\,|\,\mathbf{S})$, we concatenate gene alignments into $\mathbf{S}$ and model nucleotide sequence substitution along the unknown evolutionary history ${\cal{F}}$ through the Hasegawa et al. (1985) continuous-time Markov chain with unknown transition:transversion rate ratio $\kappa$ and stationary distribution $\boldsymbol{\pi}$. We incorporate across-site rate variation using a discretized, one-parameter Gamma distribution (Yang 1994) with unknown shape $\alpha$ and proportion $p_{\text{\tiny inv}}$ of invariant sites. To specify prior $p({\cal{F}},\kappa,\boldsymbol{\pi},\alpha,p_{\text{\tiny inv}})$, we make relatively uninformative choices, documented in the BEAST extensible markup language (XML) file in the Supplementary Material.

These triggerfish sequences and traits favor the $K=5$ factor model with a log Bayes factor of $18.5$ over the $K=4$ factor model and $6.9$ over the $K=6$ factor model (Table 1). Further, these data favor the $K=5$ factor model over the multivariate Brownian diffusion (MBD) model with a log Bayes factor of $69.7$. Even if this support were equivocal, we caution against using a MBD to model these traits. The unknown variance matrix $\Sigma$ carries $P(P+1)/2=3081$ degrees-of-freedom that dwarfs the $N\times P=1872$ possible measurements.

For 2 of the 5 factors in the $K=5$ model, Figure 3 demonstrates how fin shape changes as a function of latent factor values. We vary $\mathbf{F}_{1}$ and $\mathbf{F}_{3}$ between $-1$ and $1$ that approximates their highest posterior density range over their reconstructed evolutionary history. For $\mathbf{F}_{1}$, increasing values lead to dorsal and anal fins that become less pointed and more rounded. For $\mathbf{F}_{3}$, increasing values lead to a counterclockwise rotation of the dorsal fin. Our credible band decreases in size as the factor value gets closer to 0 since the standard deviation of the posterior inference on our loadings is multiplied by these factor values as well.

We also include the corresponding maximum clade credibility (MCC) tree, colored by factor value, with purple representing positive values and green representing negative values for the first factor $\mathbf{F}_{1},$ and the blue representing positive factor values and orange representing negative factor values for $\mathbf{F}_{3}$ in figure 4. This tree shows us that the species Balistes polylepsis and Balistes vetula, have negative factor values for $\mathbf{F}_{1}$, but those species as well as the rest of the clade with the genus Balistes and species Pseudobalistes fuscus have positive factor values for $\mathbf{F}_{3},$ whereas the clade containing the genus Rhinecanthus has negative factor values for $\mathbf{F}_{1},$ but a close to 0 factor value for $\mathbf{F}_{3}$. Conversely, the genus Xanthichthys has a negative factor value for $\mathbf{F}_{3},$ and a closer to 0 factor value for $\mathbf{F}_{1}.$ We also display posterior clade probabilities for those clades with probability $<99\%.$

For brevity, we have only considered two factors in this section. We selected $\mathbf{F}_{1}$ and $\mathbf{F}_{3}$ since these factors relate distinctive information, however we include the results for the remaining factors in the supplementary information. We additionally include our inference on the precision elements as well as our results on the inference on the other aspects of our tree model in the supplementary material.

Lastly, PFA facilitates ancestral shape reconstruction. Figure 5 depicts inferred pectoral, dorsal and anal fin shapes for ancestors of Xanthichthys mento and Balistes capriscus at arbitrary points into their evolutionary past. We choose reconstructions at the most recent common ancestor (MRCA) of all 24 species in our study and $1/4$, $1/2$ and $3/4$ of the expected sequence substitution distance between the MRCA and both contemporaneous species. Typically, high aspect ratio fins, or long fins with a small area, are associated with swimming quickly over large distances. The diet of Xanthichthys mento consists mostly of plankton and swims above reefs and has a high aspect ratio, perhaps reflecting a need to hunt down more evasive prey. We see that these low aspect ratio dorsal and anal fins arose from a moderate MRCA which flatten as the species evolved. The pectoral fin rotated clockwise as this species evolved. By contrast, Balistes capriscus has low aspect ratio dorsal and anal fins, reflecting the fact that it swims more towards the reef floors which may be more useful in navigating the complex habitat. This species evolved from a species with a moderate aspect ratio in its dorsal and anal fins which became broader and more pointed as it evolved. However, the aspect ratio increases again about $3/4$ of the way through its evolution. The pectoral fin rotated counterclockwise as it evolved.

This ancestral reconstruction can provide new insights into the trajectories of shape change that could be further investigated with biomechanical and fluid dynamic models.

4 Computational Aspects

To draw posterior inference, we simulate MCMC chains of between 200M and 1B steps, subsampling every 10K steps to eliminate unnecessary overhead and ensure the rate-limiting computation remains the PFA and L/MDB transition kernels. For path sampling, we employ 100 path points based on the quantiles of a beta $\beta\left(0.3,1\right)$ random variable (Xie et al. 2010), with warm-started chains of 10M steps at each point. In our examples, the PFA chains generate draws three- to five-fold faster than the L/MBD chains. Further, with the relatively large ratio of latent to non-latent traits in the Aquilegia example, we find an approximately 27-fold larger median effective sample size (ESS) across $\mathbf{L}$, $\mathbf{F}$ and $\gamma$ than in the latent components of $\boldsymbol{\Sigma}$, demonstrating both faster and more efficient sampling.

5 Discussion

This paper merges traditional factor analysis with phylogenetics to provide a new inference tool for comparative studies. The key connection rests on modeling each factor independently as a Brownian diffusion along a phylogeny. The tool we provide not only serves as a dimension reduction technique in the face of high-dimensional traits, but directly addresses the principal scientific questions that many comparative studies raise – specifically, how many independent evolutionary processes are driving these traits? Set in a Bayesian framework, we succeed in inferring these processes for combinations of discrete and continuous traits through model selection, while simultaneously accounting for missing measurements and possible phylogenetic uncertainty.

To make inference under PFA practical, we develop two new MCMC integration techniques. While we rely on previously proposed Gibbs samplers for integrating the loading matrix $\mathbf{L}$ and residual trait precisions $\mathbf{\Lambda}$, we require an original algorithm based on dynamic programming to integrate the factors $\mathbf{F}$ along the phylogeny efficiently. Second, we extend path sampling through a softening threshold to handle discrete traits, in which their latent support depends on the path location $\beta$. Such changing support previously has limited marginal likelihood estimation across many Bayesian models with latent random variables to combine discrete and continuous observations.

In examples involving columbine flower and fish families Poeciliidae and Balistidae evolution, inference under the PFA is notably quicker under the presence of latent traits, more interpretable and consistently favored via model selection over competing LMBD / MBD models. Interestingly, this success even holds in the Poeciliidae example, where one might expect an LMBD model to outperform. Here, the number of parameters inferred in the variance matrix is small relative to the number of parameters that form a PFA. The Poeciliidae and Balistidae examples also demonstrate our Bayesian approach’s ability to integrate missing data if we make a simple missing-at-random assumption.

Unlike many univariate comparative methods, the PFA simultaneously adjusts for correlation between all traits. This advantage reveals that some previously identified trait relationships in Poeciliidae evolution may be spurious. Further, as demonstrated in the columbine flower example, the inferred factors and their associated loadings probabilistically cluster traits into independent processes that provide additional scientific insight, often hard to discern from the correlation matrix that a LMBD model provides.

An important computational limitation of PFA arises when the number of taxa $N$ is much greater than the number of traits $P$. For the PFA, computational cost of our current MCMC integration scales as ${\cal O}\left(N^{2}K+NK^{2}P\right)$, while the cost is ${\cal O}\left(NP^{2}\right)$ for the LMBD / MBD models. Nonetheless, the Poeciliidae example carries $N/P\approx 9$ and, still, the PFA model integrates about $3\times$ more efficiently due to the example’s large ratio of latent traits. For larger $N/P$ ratios, we are currently devising algorithms that remain linear in $N$ as future work.

Arguably, PFA reaches its greatest potential when the number of traits stands large relative to the number of taxa – the reputed “large $P$, small $N$” setting. This setting arises commonly in the field of geometric morphometrics where very long series of Cartesian, (semi-) land-mark coordinate measurements define the shape of the organism. In our Balistidae example, the PFA identifies a number of independent evolutionary processes driving pectoral, dorsal and anal fin shapes. With the help of sequence data, the PFA also simultaneously infers the phylogeny and reconstructs ancestral shapes. We believe that morphometrics stands poised as a prime beneficiary of PFA.

One potential extension of this method comes from Lemey et al. (2010), where they place different diffusion rates on different branches. Additionally we can adapt the methods in Gill et al. (2016), which allows us to incorporate inference on drift in our factors whose direction changes at different points in the evolutionary process. Both of these methods are implemented in BEAST and are therefore easily adapted.

6 Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013) under Grant Agreement no. 278433-PREDEMICS and ERC Grant agreement no. 260864 and the National Institutes of Health (R01 AI107034, R01 AI117011 and R01 HG006139) and the National Science Foundation (DMS 1264153).

A Phylogenetic factor analysis Gibbs sampling

While the Gibbs samplers for a standard factor analysis are known and well documented (Lopes and West 2004), there are two aspects of our phylogenetic model that differ sufficiently to require a fresh look at how to draw posterior inference. First, our prior on $\mathbf{F}$ is based on a phylogenetic tree and therefore requires particular consideration in order to produce an efficient Gibbs sampler. Second, our inference on $K$ uses a path sampling approach where we need to infer $\mathbf{L},$ $\mathbf{F},$ and $\mathbf{\Lambda}$ at each point along the path $q^{*}\left(\beta,\mathbf{Y},\mathbf{S},\boldsymbol{\theta}\right),$ and deriving a Gibbs sampler that works for any point in the path $\beta$ will aid this process.

Sampling factors

In a standard Bayesian factor analysis, the prior on each element $F_{ij}$ is $N(0,1),$ and so the entire matrix $\mathbf{F}$ can be Gibbs sampled efficiently in a single step (Lopes and West 2004). For the phylogenetic factor analysis model, the prior on the factors is defined by Brownian motion on a phylogenetic tree as defined in (6). Thus the conditional density of $\mathbf{F}|\mathbf{Z},\mathbf{L},\mathbf{\Lambda}$ in our model is proportional to

[TABLE]

This expression does not appear to represent a distribution from which we can easily sample, principally stemming from the fact that $\mathbf{\Lambda}$ is a between-column precision and $\boldsymbol{\Psi}_{{\cal{F}}}+\kappa_{0}^{-1}\mathbf{J}$ is a between-row precision.

Fortunately, Cybis et al. (2015) devise a pre-order tree-traversal algorithm to determine the conditional distribution $\mathbf{F}_{i.}^{t}|\mathbf{F}_{-i.}$ of the factors at a single tip given all other tip values. This distribution is multivariate normal $\text{MVN}(\boldsymbol{\mu_{\mathbf{F}_{\text{-}i\cdot}}},\boldsymbol{\Lambda_{\mathbf{F}_{\text{-}i\cdot}}})$ with conditional mean $\boldsymbol{\mu_{\mathbf{F}_{\text{-}i\cdot}}}$ and conditional precision $\boldsymbol{\Lambda_{\mathbf{F}_{\text{-}i\cdot}}}$. Further, in order to numerically estimate $\mathbf{F}$ at any point along the path $q^{*}\left(\beta,\mathbf{Y},\mathbf{S},\boldsymbol{\theta}\right),$ we define

[TABLE]

Substituting in the appropriate densities and completing the square, we find that this path is proportional to

[TABLE]

where

[TABLE]

Equation (LABEL:eq:MVNkernel) is proportional to the density of a $\text{MVN}\left(\mathbf{M}\left(\beta\right)^{\left(\mathbf{F}\right)}_{i},\mathbf{V}\left(\beta\right)^{\left(\mathbf{F}\right)}_{i}\right)$; therefore, in order to sample $\mathbf{F}$ at a particular point in the path $\beta,$ we can draw a row $\mathbf{F}_{i.}$ from the distribution $\text{MVN}\left(\mathbf{M}\left(\beta\right)^{\left(\mathbf{F}\right)}_{i},\mathbf{V}\left(\beta\right)^{\left(\mathbf{F}\right)}_{i}\right)$.

Sampling loadings

The loadings matrix can be Gibbs sampled using the same method described by Lopes and West (2004) with an additional adaptation for use in path sampling. For the examples provided in this paper, we place a $N(0,1)$ prior on each cell in the loadings matrix; however, in this section we prove the Gibbs Sampler for a generic $N(\mu,\lambda)$ prior. To begin, we again define for a point on the path $\beta,$

[TABLE]

Plugging in the proper values for the sampling density and priors, rearranging and completing the square, we find that

[TABLE]

where $\mathbf{L}_{.j}=\left(L_{1j},\ldots,L_{k^{\prime}j}\right)$, $\mathbf{1}$ is a matrix of 1’s with the same dimensions as $\mathbf{L}$,

[TABLE]

Hence we find the expression in (LABEL:eq:loaddensity) is proportional to a product of independent $\text{MVN}\left(\mathbf{M}\left(\beta\right)^{\left(\mathbf{L}\right)}_{j},\mathbf{V}\left(\beta\right)^{\left(\mathbf{L}\right)}_{j}\right)$ densities. Therefore, if we wish to numerically sample a loadings column $\mathbf{L}_{.j}$ at a point on the path $\beta$ then we can sample from the distribution $\text{MVN}\left(\mathbf{M}\left(\beta\right)^{\left(\mathbf{L}\right)}_{j},\mathbf{V}\left(\beta\right)^{\left(\mathbf{L}\right)}_{j}\right).$ Since the densities across columns are independent, we may sample from them in parallel.

Sampling residual precision

We wish to sample $\mathbf{\Lambda}$ at any point in our path $q^{*}\left(\beta,\mathbf{Y},\mathbf{S},\boldsymbol{\theta}\right).$ Let $\mathbf{\Lambda}_{c}$ be a matrix equivalent to $\mathbf{\Lambda}$ with rows and columns corresponding to discrete traits removed. We then say that $\mathbf{\Lambda}_{c}=\left(\Lambda_{\left(1\right)},\ldots,\mathbf{\Lambda}_{\left(P^{\prime}\right)}\right)^{t}$ where $\Lambda_{\left(j\right)}$ models continuous trait $j$ and $P^{\prime}$ is the number of continuous traits in our model. If we define $\mathbf{L}_{c}$ and $\mathbf{Z}_{c}$ as the matrices $\mathbf{L}$ and $\mathbf{Z}$ with the columns corresponding to discrete traits removed, then we can say $\mathbf{Z}_{c}\sim\text{MVN}(\mathbf{F}\mathbf{L}_{c},\mathbf{\Lambda}_{c}).$ Our prior on $\Lambda_{\left(j\right)}$ is i.i.d. for different values of $j$ and has distribution $\Gamma(\alpha_{\mathbf{\Lambda}},\beta_{\mathbf{\Lambda}}).$ For an arbitrary point $\beta$ in our path $q^{*}\left(\beta,\mathbf{Y},\mathbf{S},\boldsymbol{\theta}\right),$ we then define

[TABLE]

with density

[TABLE]

The expression in (LABEL:eq:gamRes) is proportional to the density of a gamma $\Gamma\left(\alpha_{\mathbf{\Lambda}}+\frac{\beta N}{2},\beta_{\mathbf{\Lambda}}+\frac{\beta}{2}\,\mathbf{e}_{j}^{t}\left(\mathbf{Z}-\mathbf{F}\mathbf{L}\right)^{t}\left(\mathbf{Z}-\mathbf{F}\mathbf{L}\right)\mathbf{e}_{j}\right)$ random variable, and therefore we can sample from this gamma distribution in order to sample $\Lambda_{\left(j\right)}$ at a given point in the path.