Processing math: 49%
+ - 0:00:00
Notes for current slide
Notes for next slide

09: Phylogenetic Comparative Methods and Shape Data

1 / 58

Lecture Outline

  • On Correlated Observations and Errors
  • Multivariate Phylogenetic Comparative Methods
    1. Phylogenetic Least Squares (ANOVA/Regression/Covariation)
    2. Phylogenetic Signal
    3. Phylogenetic Ordination
    4. Comparing Evolutionary Models
2 / 58

Motivation: Correlated Observations

  • Imagine the following scenario:
    • We assemble a large database of data for many mammal species to ask:
  • Is home range size predicted by body size?

  • We might consider using a linear model to investigate: Z=Xβ+ϵ

  • But what potential problem do we face?

3 / 58

Motivation: Correlated Observations

  • Imagine the following scenario:
    • We assemble a large database of data for many mammal species to ask:
  • Is home range size predicted by body size?

  • We might consider using a linear model to investigate: Z=Xβ+ϵ

  • But what potential problem do we face?

  • The observations are not independent (species are related by a phylogeny)!
3 / 58

Comparative Biology Tradition

  • To study adaptation, biologists look comparatively across taxa, and examine trait correlations

  • Species values commonly utilized

4 / 58

Comparative Biology Tradition

  • To study adaptation, biologists look comparatively across taxa, and examine trait correlations

  • Species values commonly utilized

  • But this has the same problem: the observations are not independent!

  • To resolve this, Phylogenetic Comparative Methods are required

4 / 58

Phylogenetic Comparative Methods (PCMs)

  • PCMs condition the data on the phylogeny under an evolutionary model (i.e., account for the phylogeny during the analysis)

  • Empirical Goal: Evaluate evolutionary hypotheses while accounting for (phylogenetic) non-independence

  • This requires a phylogeny, and an evolutionary model of how trait variation is expected to accumulate
5 / 58

Phylogenetic Comparative Methods (PCMs)

  • PCMs condition the data on the phylogeny under an evolutionary model (i.e., account for the phylogeny during the analysis)

  • Empirical Goal: Evaluate evolutionary hypotheses while accounting for (phylogenetic) non-independence

  • This requires a phylogeny, and an evolutionary model of how trait variation is expected to accumulate

  • But to understand PCMs more deeply, let's first go back to something familiar: linear models

5 / 58

Linear Models Revisited

  • Statistical linear models describe patterns in data with both a deterministic component (the relationships between variables) AND a process for how we expect stochastic variation (random error) to be generated 1
  • The Linear Model: Z=Xβ+ϵ

1: For a refreshing perspective on correlated errors in biology see Ives (2022). Methods Ecol. Evol.

6 / 58

Linear Models Revisited

  • Statistical linear models describe patterns in data with both a deterministic component (the relationships between variables) AND a process for how we expect stochastic variation (random error) to be generated 1
  • The Linear Model: Z=Xβ+ϵ

1: For a refreshing perspective on correlated errors in biology see Ives (2022). Methods Ecol. Evol.

  • Xβ describes the predicted values (ˆZ) given the independent variables (X)

  • ϵ captures the stochastic error in the processes generating Z

6 / 58

Linear Models Revisited

  • Statistical linear models describe patterns in data with both a deterministic component (the relationships between variables) AND a process for how we expect stochastic variation (random error) to be generated 1
  • The Linear Model: Z=Xβ+ϵ

1: For a refreshing perspective on correlated errors in biology see Ives (2022). Methods Ecol. Evol.

  • Xβ describes the predicted values (ˆZ) given the independent variables (X)

  • ϵ captures the stochastic error in the processes generating Z

  • To this point (and in most statistics courses) we tend to focus on the deterministic component (Xβ) and whether our model is an adequate description of variation in Z (i.e., by minimizing ϵ)
6 / 58

Linear Models Revisited

  • Statistical linear models describe patterns in data with both a deterministic component (the relationships between variables) AND a process for how we expect stochastic variation (random error) to be generated 1
  • The Linear Model: Z=Xβ+ϵ

1: For a refreshing perspective on correlated errors in biology see Ives (2022). Methods Ecol. Evol.

  • Xβ describes the predicted values (ˆZ) given the independent variables (X)

  • ϵ captures the stochastic error in the processes generating Z

  • To this point (and in most statistics courses) we tend to focus on the deterministic component (Xβ) and whether our model is an adequate description of variation in Z (i.e., by minimizing ϵ)
  • But what about ϵ? What does it really represent and contain?
6 / 58

Thinking about the Error: ϵ

  • When fitting linear models, we describe ϵ as the 'error'

    • It is that portion of Z not explained by the model: ˆZ=Xβ (+ϵ)
  • Statistically, ϵ follows a distribution which describes how we expect the error to be distributed

  • For the linear models used thus far, ϵ is modeled as (assumed to be) iid 1. These are ordinary least squares (OLS) models

1: iid: independent, identically distributed error. This is where the common 'assumptions' of the OLS linear model are embodied.

7 / 58

Thinking about the Error: ϵ

  • When fitting linear models, we describe ϵ as the 'error'

    • It is that portion of Z not explained by the model: ˆZ=Xβ (+ϵ)
  • Statistically, ϵ follows a distribution which describes how we expect the error to be distributed

  • For the linear models used thus far, ϵ is modeled as (assumed to be) iid 1. These are ordinary least squares (OLS) models

1: iid: independent, identically distributed error. This is where the common 'assumptions' of the OLS linear model are embodied.

  • In essence, the random errors are modeled as if drawn independently from a normal distribution: ϵN(0,σ2)

7 / 58

Visualizing ϵ as a Matrix

  • When random error is iid, this means that elements of ϵ are drawn from the same normal distribution: ϵN(0,σ2)
  • Viewed as a matrix, the Gaussian error is: σ2Ω, where σ2 is the error variance and Ω describes the relationships among observations
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
8 / 58

Visualizing ϵ as a Matrix

  • When random error is iid, this means that elements of ϵ are drawn from the same normal distribution: ϵN(0,σ2)
  • Viewed as a matrix, the Gaussian error is: σ2Ω, where σ2 is the error variance and Ω describes the relationships among observations
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
  • The diagonal elements of Ω describe differences in the expected variance for each observation. Here they are all '1', which represents identically distributed
  • The off-diagonal elements of Ω describe expected covariation (correlation) between pairs of observations. Here they are '0', signifying independent observations

  • However, we can make the situation (and thus Ω) more complex

NOTE: In the PCM literature, Ω is often termed C or Σ. In this workshop we use Ω so as not to confuse this matrix with a trait covariance matrix Σ.
8 / 58

Modeling ϵ with Heterogeneous Variance

  • Imagine taking repeated caliper measurements of body size on 4 specimens

    • Error is expected to be proportionately greater when measuring smaller individuals (it is harder to measure smaller specimens as accurately)
9 / 58

Modeling ϵ with Heterogeneous Variance

  • Imagine taking repeated caliper measurements of body size on 4 specimens

    • Error is expected to be proportionately greater when measuring smaller individuals (it is harder to measure smaller specimens as accurately)
  • This means the expected variance differs among observations

smallest next next largest
smallest 2.7 0.0 0.0 0.0
next 0.0 2.1 0.0 0.0
next 0.0 0.0 1.8 0.0
largest 0.0 0.0 0.0 1.2
  • Now Ω models independent observations with heterogeneous random variance (i.e., σ2 varies across observations)
9 / 58

Modeling ϵ with Correlated Observations

  • Now imagine we measure 2 species from each of two genera

    • We expect that values are more similar between congeneric species (due to phylogenetic relatedness)
  • This means the expected covariance is not always zero

10 / 58

Modeling ϵ with Correlated Observations

  • Now imagine we measure 2 species from each of two genera

    • We expect that values are more similar between congeneric species (due to phylogenetic relatedness)
  • This means the expected covariance is not always zero

t1 t2 t3 t4
t1 1.0 0.7 0.0 0.0
t2 0.7 1.0 0.0 0.0
t3 0.0 0.0 1.0 0.7
t4 0.0 0.0 0.7 1.0
  • Now Ω models random variation in non-independent observations
10 / 58

Modeling ϵ with Correlated Observations

  • Now imagine we measure 2 species from each of two genera

    • We expect that values are more similar between congeneric species (due to phylogenetic relatedness)
  • This means the expected covariance is not always zero

t1 t2 t3 t4
t1 1.0 0.7 0.0 0.0
t2 0.7 1.0 0.0 0.0
t3 0.0 0.0 1.0 0.7
t4 0.0 0.0 0.7 1.0
  • Now Ω models random variation in non-independent observations

  • For actual data, two questions remain:

    1. How do we generate Ω?
    2. What do we do with it?
10 / 58

Generating Object Covariance Matrices (Ω)

  • Ω is an N×N matrix describing the expected covariation of the random errors among observations

    • If ϵ is iid, then Ω is an identity matrix
    • If ϵ is NOT iid, then Ω displays some other pattern, which was caused by 'something'
  • Biology informs us of what that something is (phylogeny, space, time)

11 / 58

Generating Object Covariance Matrices (Ω)

  • Ω is an N×N matrix describing the expected covariation of the random errors among observations

    • If ϵ is iid, then Ω is an identity matrix
    • If ϵ is NOT iid, then Ω displays some other pattern, which was caused by 'something'
  • Biology informs us of what that something is (phylogeny, space, time)

  • In this case, generating Ω requires knowledge of the relationships among observations, and a model of how trait variation is expected to accumulate

  • Let's look at this for the case of a phylogeny, and what are termed phylogenetic comparative methods (PCMs)

    • But first, we will describe the model of evolutionary change (Brownian motion)
11 / 58

Brownian Motion and PCMs

  • Brownian motion is a useful NULL model of trait change
    • Trait changes are independent from time step to time step
    • Results in: Δμ=0 and σ2taxa↑∝time

  • We can use this to intuit the expected covariation among species on a phylogeny
12 / 58

From Brownian Motion to Ω

  • Imagine a trait evolving through time
    • From the root of the phylogeny, its value changes following a BM process
    • When speciation occurs, traits in each species evolve via BM

  • The result is that traits tend to be more similar in more closely related species
    • Using BM, we can estimate Ω
13 / 58

Phylogenetic Covariance Matrices (Ω)

  • Under Brownian motion, the expected amount of change in trait values is proportional to time

    • Species evolve 'together' on their parental branch until they speciate
    • Thus, the more time before they speciate, the more similar their traits are expected to be (less time evolving separately means less opportunity to diverge)
  • This gives us the expected object covariance (Ω) under Brownian motion 1

1: Other evolutionary models, such as OU, could also be used.

14 / 58

Phylogenetic Covariance Matrices (Ω): Example

  • What is the expected covariation for species related by a phylogeny under BM?

t4 t5 t1 t2 t3
t4 1.00 0.40 0.21 0.00 0.00
t5 0.40 1.00 0.21 0.00 0.00
t1 0.21 0.21 1.00 0.00 0.00
t2 0.00 0.00 0.00 1.00 0.24
t3 0.00 0.00 0.00 0.24 1.00
  • Now we have Ω for our species

  • The next question is, what do we do with it?

15 / 58

From Ordinary Least Squares to Generalized Least Squares

  • Linear models with iid error: ϵN(0,σ2), are called Ordinary Least Squares methods

    • These include all the methods we have discussed previously
    • ANOVA, regression, ANCOVA, factorial ANOVA, multiple regression, pairwise comparisons, trajectory analysis, are part of this paradigm
16 / 58

From Ordinary Least Squares to Generalized Least Squares

  • Linear models with iid error: ϵN(0,σ2), are called Ordinary Least Squares methods

    • These include all the methods we have discussed previously
    • ANOVA, regression, ANCOVA, factorial ANOVA, multiple regression, pairwise comparisons, trajectory analysis, are part of this paradigm
  • With Ω we move from iid to correlated error

  • Linear models with correlated error: ϵN(0,σ2Ω), are called Generalized Least Squares methods

    • With GLS, we can perform versions of ANOVA, regression, ANCOVA, multiple regression, pairwise comparisons, trajectory analysis, etc.
    • These are the SAME statistical models as before, but now account for the correlations among observations as described by Ω
    • Because Ω describes covariation of the errors, we must adjust our linear model algebra to accommodate it
  • Phylogenetic Comparative Methods (PCMs) are GLS approaches

16 / 58

PCMs: General Model

  • Most PCMs use GLS model: Z=Xβ+ϵ

    • Xβ describes the predicted values (ˆZ) given the independent variable (X)
    • ϵN(0,σ2Ω) captures the stochastic error in the processes generating Z given Ω
  • Model design (X) describes the type of analysis 1

  • The Phylogenetic Comparative Method Toolkit is a set of procedures for evaluating different hypotheses

1: reviewed in: Adams and Collyer (2018;2019)

17 / 58

PCMs: An Incomplete Historical Walk

  • The conceptual development of PCMs

18 / 58

Phylogenetic Comparative Method Toolkit

19 / 58

1: Phylogenetic Least Squares: Regression (Concept)

  • Evaluate Z=Xβ+ϵ in a phylogenetic context

  • The workhorse of PCMs

20 / 58

1: Phylogenetic Least Squares: Regression and ANOVA

  • Z=Xβ+ϵ, where ϵN(0,σ2Ω)

  • Parameter estimation: 1 ˆβ=(XTΩ1X)1(XTΩ1Z).

  • Also expressed as: ˆβ=(˜XT˜X)1(˜XT˜Z)

1: Using Dean's favorite equation!

2: Grafen (1989). Phil. Trans. Roy. Soc.

21 / 58

1: Phylogenetic Least Squares: Regression and ANOVA

  • Z=Xβ+ϵ, where ϵN(0,σ2Ω)

  • Parameter estimation: 1 ˆβ=(XTΩ1X)1(XTΩ1Z).

  • Also expressed as: ˆβ=(˜XT˜X)1(˜XT˜Z)

1: Using Dean's favorite equation!

2: Grafen (1989). Phil. Trans. Roy. Soc.

  • Model evaluation steps:
ˆβ ˆZ Residuals Transformed Residuals RSS
ˆβF XFˆβF (ZˆZF) \small(\tilde{\mathbf{Z}} - \tilde{\mathbf{X}}_F\boldsymbol{\hat{\beta}}_F) \tilde{\mathbf{E}}_F^T \tilde{\mathbf{E}}_F
\small\hat{\mathbf{\beta}}_R \small\mathbf{X}_R\hat{\mathbf\beta}_R \small(\mathbf{Z}-\hat{\mathbf{Z}}_R) \small(\tilde{\mathbf{Z}} - \tilde{\mathbf{X}}_R\boldsymbol{\hat{\beta}}_R) \tilde{\mathbf{E}}_R^T \tilde{\mathbf{E}}_R
  • Obtain F and compare to RRPP empirical sampling distribution to assess significance
  • EVERYTHING is the same as before, only the phylogeny is incorporated in the analysis!
21 / 58

1: A Visual of Phylogenetic Least Squares

22 / 58

1: Phylogenetic Least Squares Regression: Example

  • Is there evolutionary allometry of head shape (does shape covary with size across species while accounting for phylogeny)?
##
## No curves detected; all points appear to be fixed landmarks.

23 / 58

1: Phylogenetic Least Squares Regression: Example (Cont.)

  • Yes, there is significant evolutionary allometry across the Plethodon phylogeny
pgls.reg <- procD.pgls(f1 = shape~svl, phy=plethtree, data=gdf, print.progress = FALSE)
summary(pgls.reg)
##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Generalized Least-Squares (via OLS projection)
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## svl 1 0.0006581 0.00065811 0.07046 3.0323 1.9503 0.028 *
## Residuals 40 0.0086814 0.00021704 0.92954
## Total 41 0.0093395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: procD.lm(f1 = shape ~ svl, iter = iter, seed = seed, RRPP = TRUE,
## SS.type = SS.type, effect.type = effect.type, int.first = int.first,
## Cov = Cov, data = data, print.progress = print.progress)
24 / 58

1: Phylogenetic Least Squares Regression: Visualization

  • Now let's visualize these results
plot.out <- plot(pgls.reg, type = "regression",
predictor = svl, reg.type = "RegScore", pch=19, col = "black")
abline(lm(plot.out$RegScore~svl))

preds <- shape.predictor(pgls.reg$GM$pgls.fitted,
x= svl, Intercept = TRUE, predmin = min(svl), predmax = max(svl))

25 / 58

1: Phylogenetic Least Squares ANOVA: Example

  • Does head shape in Plethodon differ between high and low elevation species?

26 / 58

1: Phylogenetic Least Squares ANOVA: Example (Cont.)

  • No, head shape does not differ across elevational groups when phylogeny is considered
pgls.aov <- procD.pgls(f1 = shape~elev, phy=plethtree, data=gdf, print.progress = FALSE)
summary(pgls.aov)
##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Generalized Least-Squares (via OLS projection)
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## elev 1 0.0004220 0.00042201 0.04519 1.893 1.3479 0.1
## Residuals 40 0.0089175 0.00022294 0.95481
## Total 41 0.0093395
##
## Call: procD.lm(f1 = shape ~ elev, iter = iter, seed = seed, RRPP = TRUE,
## SS.type = SS.type, effect.type = effect.type, int.first = int.first,
## Cov = Cov, data = data, print.progress = print.progress)
27 / 58

1: Phylogenetic Least Squares ANOVA: Visualization

  • Now let's visualize these results
plot.res <- gm.prcomp(shape,phy=plethtree)
plot(plot.res,phylo = FALSE, pch=21, bg=gdf$elev, cex=2)
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))

preds <- shape.predictor(arrayspecs(pgls.aov$pgls.fitted, 11,2), x= pgls.aov$X[,-1],
Intercept = TRUE, Low = Low, High = High)

28 / 58

1: Phylogenetic ANOVA: Group Aggregation

  • How groups are distributed on the phylogeny can make a difference: Adams and Collyer. Evol. (2018)

  • Too few group 'state' changes across phylogeny results in lower power
29 / 58

1: Phylogenetic Least Squares: Implementations

  • Three implementations are found in the literature

  • All yield equivalent results if performed properly
30 / 58

1B: Phylogenetic Partial Least Squares ^1

  • A related approach assesses the covariation between blocks of variables while accounting for phylogenetic non-independence
  • PLS of evolutionary covariance (i.e., rate) matrix
    1. Calculate: \mathbf{R}=\frac{(\mathbf{Z}-E(\mathbf{Z}))^{T}\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))}{N-1}
    2. Decompose: \small\mathbf{R}_{12}=\mathbf{U}_{12}\mathbf{D}\mathbf{V}^{T}_{12}
    3. Obtain r_{PLS} from correlation of singular values, \mathbf{P}_1 and \mathbf{P}_2
    4. Significance evaluated using RRPP

1: Adams and Felice. PLoS One (2014); Adams and Collyer. Syst. Biol. (2018)

2: Use of PICs yields same r_{PLS} (Klingenberg and Marugan-Loban. Syst. Biol. 2013); but significance tests based on permuting PICs is incorrect (see Adams and Collyer 2015. Evolution).

31 / 58

1B: Phylogenetic Partial Least Squares ^1

  • A related approach assesses the covariation between blocks of variables while accounting for phylogenetic non-independence
  • PLS of evolutionary covariance (i.e., rate) matrix
    1. Calculate: \mathbf{R}=\frac{(\mathbf{Z}-E(\mathbf{Z}))^{T}\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))}{N-1}
    2. Decompose: \small\mathbf{R}_{12}=\mathbf{U}_{12}\mathbf{D}\mathbf{V}^{T}_{12}
    3. Obtain r_{PLS} from correlation of singular values, \mathbf{P}_1 and \mathbf{P}_2
    4. Significance evaluated using RRPP

1: Adams and Felice. PLoS One (2014); Adams and Collyer. Syst. Biol. (2018)

2: Use of PICs yields same r_{PLS} (Klingenberg and Marugan-Loban. Syst. Biol. 2013); but significance tests based on permuting PICs is incorrect (see Adams and Collyer 2015. Evolution).

  • This may be obtained directly from the phylo-transformed data:
    1. Decompose: \tilde{\mathbf{Z}}_1^T\tilde{\mathbf{Z}}_2=\mathbf{U}_{12}\mathbf{D}\mathbf{V}^{T}_{12}
    2. Obtain r_{PLS} from correlation of singular values, \mathbf{P}_1 and \mathbf{P}_2 from phylogenetically-transformed data
    3. Significance evaluated using RRPP
31 / 58

1B: Phylogenetic Partial Least Squares: Example

  • Is there phenotypic integration between the cranium and mandible in Plethodon salamanders?
land.gps<-c("A","A","A","A","A","B","B","B","B","B","B")
PLS.Y <- phylo.integration(A = gdf$shape, partition.gp = land.gps, phy= plethtree, print.progress = FALSE)
summary(PLS.Y)
##
## Call:
## phylo.integration(A = gdf$shape, phy = plethtree, partition.gp = land.gps,
## print.progress = FALSE)
##
##
##
## r-PLS: 0.843
##
## Effect Size (Z): 4.7366
##
## P-value: 0.001
##
## Based on 1000 random permutations
32 / 58

1B: Phylogenetic Partial Least Squares: Visualization

  • Now let's visualize these results
plot(PLS.Y)

33 / 58

2: Phylogenetic Signal

  • The degree to which phenotypic similarity associates with phylogenetic relatedness

34 / 58

2: Phylogenetic Signal Summary Measures

  • Kappa statistic: Blomberg et al. Evol. (2003)

\small{K=\frac{(\mathbf{Z}-E(\mathbf{Z}))^T(\mathbf{Z}-E(\mathbf{Z}))}{(\mathbf{Z}-E(\mathbf{Z}))^T\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))}/\frac{tr(\mathbf\Omega)-N(\mathbf{1}^T\mathbf{\Omega1})^{-1}}{N-1}}

35 / 58

2: Phylogenetic Signal Summary Measures

  • Kappa statistic: Blomberg et al. Evol. (2003)

\small{K=\frac{(\mathbf{Z}-E(\mathbf{Z}))^T(\mathbf{Z}-E(\mathbf{Z}))}{(\mathbf{Z}-E(\mathbf{Z}))^T\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))}/\frac{tr(\mathbf\Omega)-N(\mathbf{1}^T\mathbf{\Omega1})^{-1}}{N-1}}

  • \lambda statistic: Pagel Nature (1999)

\small{logL=log(\frac{(\mathbf{Z}-E(\mathbf{Z}))^T\mathbf\Omega^\lambda{^{-1}}(\mathbf{Z}-E(\mathbf{Z})))}{\sqrt{(2\pi)^{Np}\times|\mathbf\Omega^\lambda|}}}

  • Branch-length scaling during ML fit of data to phylogeny
35 / 58

2: Phylogenetic Signal Summary Measures

  • Kappa statistic: Blomberg et al. Evol. (2003)

\small{K=\frac{(\mathbf{Z}-E(\mathbf{Z}))^T(\mathbf{Z}-E(\mathbf{Z}))}{(\mathbf{Z}-E(\mathbf{Z}))^T\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))}/\frac{tr(\mathbf\Omega)-N(\mathbf{1}^T\mathbf{\Omega1})^{-1}}{N-1}}

  • \lambda statistic: Pagel Nature (1999)

\small{logL=log(\frac{(\mathbf{Z}-E(\mathbf{Z}))^T\mathbf\Omega^\lambda{^{-1}}(\mathbf{Z}-E(\mathbf{Z})))}{\sqrt{(2\pi)^{Np}\times|\mathbf\Omega^\lambda|}}}

  • Branch-length scaling during ML fit of data to phylogeny
  • Both are for univariate response data only!
35 / 58

2: Multivariate Phylogenetic Signal ^1

  • A multivariate generalization of Kappa has been proposed: \small{K}_{mult}
    • Approach based on the distance-covariance equivalency
    • Uses distances between species in the original and phylo-transformed dataspaces

K_{mult}=\frac{\mathbf{D}_{\mathbf{Z},E(\mathbf{Z})}^{T}\mathbf{D}_{\mathbf{Z},E(\mathbf{Z})}}{\mathbf{PD}_{\tilde{\mathbf{Z}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Z}},0}}/\frac{tr(\mathbf\Omega)-N(\mathbf{1}^{T}\mathbf{\Omega1})^{-1}}{N-1}

  • \small\mathbf{D}_{\mathbf{Z},E(\mathbf{Z})} is the distance from each object to the phylogenetic mean
  • \small\mathbf{PD}_{\tilde{\mathbf{Z}},0} is the distance between each object in the phylogenetically-transformed space (found as: \small\tilde{\mathbf{Z}}=\mathbf{P}(\mathbf{Z-E(\mathbf{Z})})) and the origin
36 / 58

2: Multivariate Phylogenetic Signal ^1

  • A multivariate generalization of Kappa has been proposed: \small{K}_{mult}
    • Approach based on the distance-covariance equivalency
    • Uses distances between species in the original and phylo-transformed dataspaces

K_{mult}=\frac{\mathbf{D}_{\mathbf{Z},E(\mathbf{Z})}^{T}\mathbf{D}_{\mathbf{Z},E(\mathbf{Z})}}{\mathbf{PD}_{\tilde{\mathbf{Z}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Z}},0}}/\frac{tr(\mathbf\Omega)-N(\mathbf{1}^{T}\mathbf{\Omega1})^{-1}}{N-1}

  • \small\mathbf{D}_{\mathbf{Z},E(\mathbf{Z})} is the distance from each object to the phylogenetic mean
  • \small\mathbf{PD}_{\tilde{\mathbf{Z}},0} is the distance between each object in the phylogenetically-transformed space (found as: \small\tilde{\mathbf{Z}}=\mathbf{P}(\mathbf{Z-E(\mathbf{Z})})) and the origin

This is the same as:

K_{mult}=\frac{tr\left((\mathbf{Z}-E(\mathbf{Z}))^T(\mathbf{Z}-E(\mathbf{Z})\right)}{tr\left((\mathbf{Z}-E(\mathbf{Z}))^T\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))\right)}/\frac{tr(\mathbf\Omega)-N(\mathbf{1}^T\mathbf{\Omega1})^{-1}}{N-1}

1: Adams (2014). Syst. Biol.

36 / 58

2: Multivariate Phylogenetic Signal: Example

Does head shape in in Plethodon display phylogenetic signal?

PS.shape <- physignal(A=shape,phy=plethtree,iter=999, print.progress = FALSE)
summary(PS.shape)
##
## Call:
## physignal(A = shape, phy = plethtree, iter = 999, print.progress = FALSE)
##
##
##
## Observed Phylogenetic Signal (K): 0.392
##
## P-value: 0.001
##
## Based on 1000 random permutations
##
## Use physignal.z to estimate effect size.
plot(PS.shape)

37 / 58

2: Phylogenetic Signal: An Effect Size ^1

  • Optimizing \Omega by \lambda in the equation of Kappa provides a direct link between the two statistics ^2
  • This formulation lends itself nicely to evaluation via RRPP

    • RRPP many times, obtaining logL

    • Calculate Effect size as:

    \small{Z=\frac{\theta_{obs}-\mu_\theta}{\sigma_\theta}} where \theta_{obs}=f(logL_\Omega)

1: See Collyer et al. (2022) Methods Ecol. Evol.

2: Blomberg et al. (2003) Evol.

38 / 58

2: Phylogenetic Signal: An Effect Size ^1

  • Optimizing \Omega by \lambda in the equation of Kappa provides a direct link between the two statistics ^2
  • This formulation lends itself nicely to evaluation via RRPP

    • RRPP many times, obtaining logL

    • Calculate Effect size as:

    \small{Z=\frac{\theta_{obs}-\mu_\theta}{\sigma_\theta}} where \theta_{obs}=f(logL_\Omega)

1: See Collyer et al. (2022) Methods Ecol. Evol.

2: Blomberg et al. (2003) Evol.

  • One can now compare effect sizes to evaluate the strength of phylogenetic signal across datasets:

\small{Z_{12}=\frac{(\theta_{obs_1}-\mu_{\theta_1})-(\theta_{obs_2}-\mu_{\theta_2})}{\sqrt{\sigma^2_{\theta_1}+\sigma^2_{\theta_2}}}}

38 / 58

2: Comparing Phylogenetic Signal: Example

Greater phylogenetic signal in SA:V ratios!

From Collyer et al. (2022). Methods Ecol. Evol.
39 / 58

3: Phylogenetic Ordination

  • Ordinations that incorporate phylogenetic relatedness

  • Intent is to provide visual insights into macroevolutionary patterns

  • Three approaches:

40 / 58

3: Phylogenetic Ordination

  • Ordinations that incorporate phylogenetic relatedness

  • Intent is to provide visual insights into macroevolutionary patterns

  • Three approaches:

  1. Phylomorphospace: project phylogeny into PCA space (using ancestral states)

  2. Phylogenetic PCA (pPCA): account for phylogeny in PCA computations

    • Rotates data 'away' from phylogenetic signal (i.e., pPC_1 is uncorrelated with phylogenetic signal)
  3. Phylogenetically-Aligned Components Analysis (PACA)

    • Rotates data 'towards' phylogenetic signal (i.e., PACA_1 maximizes phylogenetic signal)
40 / 58

3: Phylogenetic Ordination ^1

  • Recall the general ordination equation (Shape Statistics I lecture):

\mathbf{I}^T_{n \times n} \mathbf{Z}=\mathbf{UDV}^T

where: \mathbf{A} = \mathbf{I}_{n \times n}. This is PCA, and finds the principal directions of variation in the data independent of any other association.

1: See Collyer and Adams (2021). Methods Ecol. Evol.

41 / 58

3: Phylogenetic Ordination ^1

  • Recall the general ordination equation (Shape Statistics I lecture):

\mathbf{I}^T_{n \times n} \mathbf{Z}=\mathbf{UDV}^T

where: \mathbf{A} = \mathbf{I}_{n \times n}. This is PCA, and finds the principal directions of variation in the data independent of any other association.

1: See Collyer and Adams (2021). Methods Ecol. Evol.

But what if we do:

\mathbf{T^T_\Omega}\mathbf{Z} =\mathbf{UDV}^T

where: \mathbf{A}=\mathbf{T_\Omega} and \mathbf\Omega^{-1}=\mathbf{TT^T}. Data is aligned to directions independent of phylogenetic signal (the least phylogenetic signal in the first component).

41 / 58

3: Phylogenetic Ordination ^1

  • Recall the general ordination equation (Shape Statistics I lecture):

\mathbf{I}^T_{n \times n} \mathbf{Z}=\mathbf{UDV}^T

where: \mathbf{A} = \mathbf{I}_{n \times n}. This is PCA, and finds the principal directions of variation in the data independent of any other association.

1: See Collyer and Adams (2021). Methods Ecol. Evol.

But what if we do:

\mathbf{T^T_\Omega}\mathbf{Z} =\mathbf{UDV}^T

where: \mathbf{A}=\mathbf{T_\Omega} and \mathbf\Omega^{-1}=\mathbf{TT^T}. Data is aligned to directions independent of phylogenetic signal (the least phylogenetic signal in the first component).

Or if instead we perform:

\mathbf{\Omega}^T\mathbf{Z} = \mathbf{UDV}^T

where \mathbf{A}=\mathbf{\Omega}. Data is aligned to directions that maximize phylogenetic signal (most phylogenetic signal in the first few components).

41 / 58

3: Phylogenetic Ordination ^1

  • Recall the general ordination equation (Shape Statistics I lecture):

\mathbf{I}^T_{n \times n} \mathbf{Z}=\mathbf{UDV}^T

where: \mathbf{A} = \mathbf{I}_{n \times n}. This is PCA, and finds the principal directions of variation in the data independent of any other association.

1: See Collyer and Adams (2021). Methods Ecol. Evol.

But what if we do:

\mathbf{T^T_\Omega}\mathbf{Z} =\mathbf{UDV}^T

where: \mathbf{A}=\mathbf{T_\Omega} and \mathbf\Omega^{-1}=\mathbf{TT^T}. Data is aligned to directions independent of phylogenetic signal (the least phylogenetic signal in the first component).

Or if instead we perform:

\mathbf{\Omega}^T\mathbf{Z} = \mathbf{UDV}^T

where \mathbf{A}=\mathbf{\Omega}. Data is aligned to directions that maximize phylogenetic signal (most phylogenetic signal in the first few components).

Mathematically the difference between these phylogenetic ordination procedures is that they align the data relative to different matrices in \mathbf{A}.

41 / 58

3: Phylogenetic Ordination: Phylomorphospace (PCA)

  • Phylomorphospace is a PCA with the phylogeny superimposed (via projection of ancestral states)

\mathbf{I}^T_{n \times n} \mathbf{Z}=\mathbf{UDV}^T

plot.pca <- gm.prcomp(shape,phy=plethtree)
plot(plot.pca,phylo = TRUE, pch=21, bg=gdf$elev, cex=2, phylo.par = list(tip.labels = FALSE, node.labels = FALSE) )
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))

42 / 58

3: Phylogenetic Ordination: Phylogenetic PCA (pPCA)

  • Aligns data to directions independent of phylogenetic signal

\mathbf{T^T_\Omega}\mathbf{Z} =\mathbf{UDV}^T

plot.ppca <- gm.prcomp(shape,phy=plethtree, GLS = TRUE, transform = FALSE)
plot(plot.ppca,phylo = TRUE, pch=21, bg=gdf$elev, cex=2, phylo.par = list(tip.labels = FALSE, node.labels = FALSE) )
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))

43 / 58

3: Phylogenetic Ordination: Phylogenetically Aligned Components (PACA)

  • Aligns data to directions that maximize phylogenetic signal

\mathbf{\Omega^T}\mathbf{Z} =\mathbf{UDV}^T

plot.paca <- gm.prcomp(shape,phy=plethtree, align.to.phy = TRUE)
plot(plot.paca,phylo = TRUE, pch=21, bg=gdf$elev, cex=2, phylo.par = list(tip.labels = FALSE, node.labels = FALSE) )
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))

44 / 58

3: Phylogenetic Ordination: Side By Side

45 / 58

3: Phylogenetic Ordination: Comments

  • Methods provide contrasting views of the data
    • Phylomorphospace is a 'pure' view of the space
    • pPCA and PACA rotate relative to the phylogeny (albeit differently)

  • Utilizing several may be helpful
46 / 58

4: Comparing Evolutionary Models

  • What evolutionary model best describes trait variation?
  • Models fit as: \mathbf{Z=1\beta+\epsilon}, where \epsilon\small\sim\mathcal{N}(0,\sigma^2\mathbf\Omega)

    • Differing evolutionary models are embodied in different \mathbf\Omega

47 / 58

4: Comparing Evolutionary Models

  • What evolutionary model best describes trait variation?
  • Models fit as: \mathbf{Z=1\beta+\epsilon}, where \epsilon\small\sim\mathcal{N}(0,\sigma^2\mathbf\Omega)

    • Differing evolutionary models are embodied in different \mathbf\Omega

  • Model comparisons of:
  1. Evolutionary rate models: BM1, BMM, etc.
  2. Evolutionary 'mode' models: BM1, OU1, OUM, etc.
  • Major issue: multivariate generalizations for some evolutionary model comparisons remain underdeveloped (see Adams and Collyer (2018;2019))
47 / 58

4: Comparing Evolutionary Rates Among Clades

Is there evidence for multiple evolutionary rates on the phylogeny?

  • Fit data to phylogeny under single-rate and multi-rate models, where \small\log{L}=-\frac{np}{2}ln{2\pi}-\frac{n}{2}ln{\begin{vmatrix}{(\mathbf{\sigma^2\Omega})}\end{vmatrix}}-\frac{1}{2}tr{(\mathbf{Z}-E(\mathbf{Z}))^{T}(\mathbf{\sigma^2\Omega})^{-1}(\mathbf{Z}-E(\mathbf{Z}))} Perform LRT, AIC, simulation, or permutation-based (RRPP) evaluation
48 / 58

4: Comparing Multivariate Evolutionary Rates Among Clades

  • Approach has been extended to multivariate traits (using distance-based procedure of Adams 2014)
  • \sigma^2_{mult}=\mathbf{PD}_{\tilde{\mathbf{Z}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Z}},0}/N-1
    • RRPP used for model evaluation
49 / 58

4: Comparing Multivariate Evolutionary Rates Among Clades

  • Approach has been extended to multivariate traits (using distance-based procedure of Adams 2014)
  • \sigma^2_{mult}=\mathbf{PD}_{\tilde{\mathbf{Z}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Z}},0}/N-1
    • RRPP used for model evaluation

Does head shape in high-elevation Plethodon species evolve more rapidly than in low-elevation species?

ER<-compare.evol.rates(A=gdf$shape, phy=plethtree,gp=gdf$elev,iter=999, method = 'permutation',print.progress = FALSE)
summary(ER)
##
## Call:
##
##
## Observed Rate Ratio: 1.0133
##
## P-value: 0.9645
##
## Effect Size: -1.5892
##
## Based on 1000 random permutations
##
## The rate for group High is 1.00093880316058e-05
##
## The rate for group Low is 1.01425751777744e-05

49 / 58

4: Comparing Evolutionary Rates Among Traits

  • Is there evidence that one trait evolves faster than another?
  • Calculate evolutionary rate matrix:

\small{\mathbf{R}=\frac{(\mathbf{Z}-E(\mathbf{Z}))^{T}\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))}{N-1}=\left( \begin{array}{ccc} \sigma^2_{1} & \sigma^2_{1,2} & \sigma^2_{1,3} \\ \sigma^2_{2,1} & \sigma^2_{2} & \sigma^2_{2,3} \\ \sigma^2_{3,1} & \sigma^2_{3,2} & \sigma^2_{3} \end{array} \right)}

  • Obtain \small{logL}: \small\log{L}=-\frac{np}{2}ln{2\pi}-\frac{n}{2}ln{\begin{vmatrix}{(\mathbf{R\otimes{\Omega}})}\end{vmatrix}}-\frac{1}{2}tr{(\mathbf{Z}-E(\mathbf{Z}))^{T}(\mathbf{R\otimes{\Omega}})^{-1}(\mathbf{Z}-E(\mathbf{Z}))}

1: Adams (2013). Syst. Biol.; Denton and Adams (2015). Evol.

50 / 58

4: Comparing Evolutionary Rates Among Traits

  • Is there evidence that one trait evolves faster than another?
  • Calculate evolutionary rate matrix:

\small{\mathbf{R}=\frac{(\mathbf{Z}-E(\mathbf{Z}))^{T}\mathbf\Omega^{-1}(\mathbf{Z}-E(\mathbf{Z}))}{N-1}=\left( \begin{array}{ccc} \sigma^2_{1} & \sigma^2_{1,2} & \sigma^2_{1,3} \\ \sigma^2_{2,1} & \sigma^2_{2} & \sigma^2_{2,3} \\ \sigma^2_{3,1} & \sigma^2_{3,2} & \sigma^2_{3} \end{array} \right)}

  • Obtain \small{logL}: \small\log{L}=-\frac{np}{2}ln{2\pi}-\frac{n}{2}ln{\begin{vmatrix}{(\mathbf{R\otimes{\Omega}})}\end{vmatrix}}-\frac{1}{2}tr{(\mathbf{Z}-E(\mathbf{Z}))^{T}(\mathbf{R\otimes{\Omega}})^{-1}(\mathbf{Z}-E(\mathbf{Z}))}

1: Adams (2013). Syst. Biol.; Denton and Adams (2015). Evol.

  • Estimate \small\mathbf{R}_C and \small{logL}_C where diagonal of \mathbf{R} is constrained to be equal:

\small\mathbf{R}_C=\left( \begin{array}{ccc} \sigma^2_{p} & & \\ \sigma^2_{2,1} & \sigma^2_{p} & \\ \sigma^2_{3,1} & \sigma^2_{3,2} & \sigma^2_{p} \end{array} \right)

  • Compare \small{logL}_C against \small{logL}
50 / 58

4: Comparing Evolutionary Rates Among Traits: Example

51 / 58

4: Comparing Evolutionary Rates Among Multivariate Traits

  • Approach has been extended to multivariate traits (using distance-based procedure of Adams 2014)
  • \sigma^2_{mult}=\mathbf{PD}_{\tilde{\mathbf{Z}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Z}},0}/N-1
    • Rate-ratio used as test measure: R=\frac{max(\sigma^2_{mult})}{min(\sigma^2_{mult})}
    • Permutation tests used for statistical evaluation
EMR <- compare.multi.evol.rates(A=gdf$shape, phy=plethtree, gp=c(rep(1,5),rep(2,6)), print.progress = FALSE)
summary(EMR)
##
## Call:
##
##
## Observed Rate Ratio: 1.0826
##
## P-value: 0.975
##
## Effect Size: -1.943
##
## Based on 1000 random permutations
##
## The rate for group 1 is 9.67170500545254e-06
##
## The rate for group 2 is 1.04710160170648e-05

52 / 58

4: Comparing Multivariate Evolutionary Rates: Other Approaches ^1

  • Likelihood-based model comparisons comparing rates for multivariate data have also been proposed

  • logL (various authors): Methods work so long as N>>p.

    • When this is not the case, |\mathbf{R\otimes\Omega}| is singular, and the logL cannot be completed.

1: see Adams and Collyer (2018;2019)

53 / 58

4: Comparing Multivariate Evolutionary Rates: Other Approaches ^1

  • Likelihood-based model comparisons comparing rates for multivariate data have also been proposed

  • logL (various authors): Methods work so long as N>>p.

    • When this is not the case, |\mathbf{R\otimes\Omega}| is singular, and the logL cannot be completed.

1: see Adams and Collyer (2018;2019)

  • PCL (Pairwise Composite Likelihood: Goolsby 2016)

    • Method obtains set of logL for pairs of traits and sums them

    • Method has high type I error rate which differs by the rotation of the data!

53 / 58

4: Comparing Evolutionary Models: Evolutionary 'Modes'

  • Compare models that describe the manner in which trait variation accumulates over evolutionary time

  • BM1, BMM, OU (Ornstein-Uhlenbeck), Early Burst (EB), ACDC, etc.

  • Fit data to phylogeny under differing evolutionary models and compare fit (LRT, AIC, simulation, RRPP, etc.)

54 / 58

4: Comparing Evolutionary 'Modes': Univariate Example

  • How did Anolis body size groups evolve?
    • 5 models: BM, \small{OU}_1, \small{OU}_3 \small{OU}_4 (3 group+anc), \small{OU}_{LP} (3 gp + history of colonization)

  • \small{OU}_{LP} (3 gp + col. hist.) best explains body size evolution
55 / 58

4: Comparing Evolutionary 'Modes': Multivariate Data

  • Several methods proposed for comparing multivariate evolutionary models: logL_{mult}, PCL, \sum{logL_{indiv}}

  • Unfortunately, nearly all current implementations for multivariate OU models display high misspecification rates

  • More analytical development needed here
56 / 58

4: Comparing Evolutionary 'Modes': Recent Developments

  • New developments in logL_{mult}
  • Bartoszek et al. (2024) updated algorithms in mvSlouch ^1

  • Misspecification levels now appropriate, and method rotation invariant

  • Provides a path forward for multivariate OU models (though still requires n>p and no redundant dimensions)

1: took our suggestions to heart!

57 / 58

Phylogenetic Comparative Methods: Conclusions

  • Phylogenetic comparative methods provide useful tools to account for phylogenetic non-independence among observations

  • Hypotheses one can evaluate with shape data include:

    • phylogenetic ANOVA

    • Phylogenetic regression

    • Phylogenetic covariation (PLS)

    • Phylogenetic ordination (PCA, pPCA, PACA)

    • Phylogenetic signal (development regarding \lambda_{mult} ongoing)

    • Comparing evolutionary rates

    • Comparing evolutionary modes (some advance: requires more analytical development)

58 / 58

Lecture Outline

  • On Correlated Observations and Errors
  • Multivariate Phylogenetic Comparative Methods
    1. Phylogenetic Least Squares (ANOVA/Regression/Covariation)
    2. Phylogenetic Signal
    3. Phylogenetic Ordination
    4. Comparing Evolutionary Models
2 / 58
Paused

Help

Keyboard shortcuts

, , Pg Up, k Go to previous slide
, , Pg Dn, Space, j Go to next slide
Home Go to first slide
End Go to last slide
Number + Return Go to specific slide
b / m / f Toggle blackout / mirrored / fullscreen mode
c Clone slideshow
p Toggle presenter mode
t Restart the presentation timer
?, h Toggle this help
Esc Back to slideshow