Dean Adams, Iowa State University
Taxonomic Diversity
Taxonomic Diversity
Morphological Diversity
Taxonomic Diversity
Morphological Diversity
How do we characterize patterns, and hypothesize processes?
Trait correlations often used to study coevolution and adaptation
Species values commonly utilized
The problem?
Species are not independent of one another
Phylogenetic comparative methods condition the data on the phylogeny during the analysis
Allows one to assess trait covariation while accounting for the non-independence due to shared evolutionary history
The conceptual development of PCMs
Phylogenetic regression/anova/association models
Phylogenetic Generalized Least Squares methods
Phylogenetic PLS
Phylogenetic ordination
Phylomorphospace
Phylogenetic PCA
Phylogenetically Aligned Component Analysis (PACA)
Exploring and modeling evolutionary processes
Phylogenetic signal
Evolutionary rates
Evolutionary models
PCMs condition the data on the phylogeny under an evolutionary model
Data conditioning: Account for phylogenetic non-independence
Evolutionary model: How trait variation is expected to accumulate
Most PCMs use GLS (generalized least squares) as a model:
\[\small\mathbf{Y}=\mathbf{X{\hat{\beta}}+\epsilon}\] Here, \(\small\epsilon\) is not iid but is \(\small\sim\mathcal{N}(0,\textbf{V})\): containing expected covariation between taxa as described by the phylogeny (in matrix \(\small\mathbf{V}\))
\(\small\mathbf{V}\) is the phylogenetic covariance matrix
Describes the amount of evolutionary time species share via common ancestry (and thus how similar they are expected to be)
\[\small\mathbf{Y}=\mathbf{X{\hat{\beta}}+\epsilon}\]
Model design (\(\small\mathbf{X}\)) describes the type of analysis
Parameters of model (and model significance) obtained in various ways
Evaluate \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) in a phylogenetic context
The workhorse of PCMs
ANOVA and regression models that account for phylogeny
Requires a model of evolutionary change: typically Brownian motion (BM)
ANOVA and regression models that account for phylogeny
Requires a model of evolutionary change: typically Brownian motion (BM)
\(\Delta\mu = 0\)
\(\sigma^2\) (variance among taxa) \(\uparrow\) \(\propto\) time
ANOVA and regression models that account for phylogeny
Requires a model of evolutionary change: typically Brownian motion (BM)
\(\Delta\mu = 0\)
\(\sigma^2\) (variance among taxa) \(\uparrow\) \(\propto\) time
Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)
Several implementations possible (yield identical results if implemented properly)
Estimate contrast scores between pairs of taxa (tips or nodes)
Use contrasts for analyses (OLS solution)
Coefficients found as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}_{pic} \mathbf{X}_{pic}\right )^{-1}\left ( \mathbf{X}^{T}_{pic} \mathbf{Y}_{pic}\right )\)
What is the association between Y and X?
What is the association between Y and X?
## Analysis of Variance Table
##
## Response: pic.y
## Df Sum Sq Mean Sq F value Pr(>F)
## pic.x 1 14.3519 14.3519 19.285 0.00461 **
## Residuals 6 4.4651 0.7442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pic.x
## 0.8846971
Use GLS model with non-independent error structure
Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)
Use GLS model with non-independent error structure
Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)
\(\textbf{V}\) is the phylogenetic covariance matrix
Describes expected covariance among taxa due to shared evolutionary history (typically under BM)
Coefficients found as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}^{-1}\mathbf{Y}\right )\)
## Denom. DF: 6
## numDF F-value p-value
## (Intercept) 1 12.87792 0.0115
## X 1 19.28544 0.0046
## (Intercept) X
## -2.3296557 0.8846971
Identical results to PICs!
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]
Weighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} (v_1+v_2) & v_{12} & 0 \\ v_{12} & (v_1+v_2) & 0 \\ 0 & 0 & 1 \end{array} \right) \]
In PGLS, the weights are the phylogenetic distances, which describe the phylogenetic non-independence
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]
Weighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} (v_1+v_2) & v_{12} & 0 \\ v_{12} & (v_1+v_2) & 0 \\ 0 & 0 & 1 \end{array} \right) \]
In PGLS, the weights are the phylogenetic distances, which describe the phylogenetic non-independence
Utilize OLS transformation of GLS models
Phylogenetic transformation matrix: \(\small\mathbf{T}=\left ( \mathbf{U}\mathbf{W}^{-1/2} \mathbf{U}^{T}\right )^{-1}\)
U and W are eigenvectors and eigenvalues of V
Transformed data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\)
Transformed design matrix: \(\small\tilde{\mathbf{X}}=\mathbf{TX}\)
Coefficients solved as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{X}}\right )^{-1}\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{Y}}\right )\)
##
## 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)
## X 1 14.3519 14.3519 0.76271 19.285 2.9475 0.002 **
## Residuals 6 4.4651 0.7442 0.23729
## Total 7 18.8170
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ X, data = df, Cov = df$Cov, print.progress = FALSE)
## (Intercept) X
## -2.3296557 0.8846971
Identical results to PICs & PGLS!
One can also use linear mixed models (PGLMM: phylogenetic generalized linear mixed models) to evaluate trait associations
Here the phylogeny is considered a random effect (sense Lynch 1991)
Results from MCMC runs are analogous (though not always identical) to the linear model solutions
## (Intercept) X
## -2.3296557 0.8846971
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -1.4721864 -5.1694902 1.813651 98 0.38775510
## X 0.8097609 0.4361165 1.224440 98 0.01020408
PIC, PGLS, and Phylogenetic transform yield identical parameter estimates
\[\tiny \hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}^{-1}\mathbf{Y}\right )\]
\[\tiny =\left ( \mathbf{X}^{T}_{pic} \mathbf{X}_{pic}\right )^{-1}\left ( \mathbf{X}^{T}_{pic} \mathbf{Y}_{pic}\right )\]
\[\tiny =\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{X}}\right )^{-1}\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{Y}}\right )\]
Model significance typically accomplished using parametric procedures
F-ratios compared to parametric F-distribution
Optimize \(\log{L}\) for model
\[\tiny \log{L}=\log{ \left(\frac{exp{\left({-\frac{1}{2}{\left({\mathbf{Y}-E(\mathbf{Y})}\right)^{T}}} \mathbf{V}^{-1}\left({\mathbf{Y}-E(\mathbf{Y})}\right)\right)}} {\sqrt{\left({2\pi}\right)^{Np}\times\begin{vmatrix} \mathbf{V} \end{vmatrix}}}\right)}\]
PROBLEM: Parametric significance testing suffers from Rao’s paradox
Power reduces as p-dimensions increase
Another solution required
Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)
Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\)
Permute residuals from reduced model (\(\small\tilde{\mathbf{E}}_{R}\)), add to predicted values (\(\small\hat{\mathbf{Y}}_{R}\)), assess model, repeat
Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)
Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\)
Permute residuals from reduced model (\(\small\tilde{\mathbf{E}}_{R}\)), add to predicted values (\(\small\hat{\mathbf{Y}}_{R}\)), assess model, repeat
Appropriate for ANOVA, regression and more complex PGLS models with high-dimensional data (Adams & Collyer. Evol. 2018)
RULE: Transform first, permute second!
Are body dimensions associated with overall size across Plethodon salamander species?
Are body dimensions associated with overall size across Plethodon salamander species?
PGLS.reg<-procD.pgls(Y ~ size, phy = phy, data = gdf, iter = 999, print.progress = FALSE)
PGLS.reg$aov.table## Df SS MS Rsq F Z Pr(>F)
## size 1 214.21 214.214 0.35852 19.561 3.0528 0.001 **
## Residuals 35 383.29 10.951 0.64148
## Total 36 597.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Do body dimensions differ by salamander ecotype?
PGLS.reg<-procD.pgls(Y ~ gp, phy = phy, data = gdf, iter = 999, print.progress = FALSE)
PGLS.reg$aov.table## Df SS MS Rsq F Z Pr(>F)
## gp 1 23.44 23.435 0.03922 1.4288 0.77209 0.236
## Residuals 35 574.07 16.402 0.96078
## Total 36 597.50
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
Assess correlation between two blocks of variables in a phylogenetic context
PLS of evolutionary covariance (rate) matrix
\[\small\mathbf{R}=\frac{\left(\mathbf{Y}-E(\mathbf{Y})\right)^{T}\mathbf{V}^{-1} \left(\mathbf{Y}-E(\mathbf{Y})\right)} {N-1}\]
SVD of \(\small\mathbf{R}_{12}=\mathbf{U}_{12}\mathbf{D}\mathbf{V}^{T}_{12}\)
Assess correlation between two blocks of variables in a phylogenetic context
PLS of evolutionary covariance (rate) matrix
\[\small\mathbf{R}=\frac{\left(\mathbf{Y}-E(\mathbf{Y})\right)^{T}\mathbf{V}^{-1} \left(\mathbf{Y}-E(\mathbf{Y})\right)} {N-1}\]
SVD of \(\small\mathbf{R}_{12}=\mathbf{U}_{12}\mathbf{D}\mathbf{V}^{T}_{12}\)
Equivalently found from PLS of \(\small\Sigma\) from \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\)
Significance from RRPP of \(\small\tilde{\mathbf{Y}}\)
Do head traits correlate with body traits (same data)?
PLS.Y <- phylo.integration(A = Y[,2:3], A2 = Y[,c(1,4:6)], phy= plethtree, print.progress = FALSE)
PLS.Y$r.pls; PLS.Y$P.value## [1] 0.6133857
## obs
## 0.001
Is there evolutionary integration across species between cranial and mandible morphology?
Yes
Visual inspection of phylomorphospaces useful for hypothesis generation
The degree to which phenotypic similarity associates with phylogenetic relatedness
Kappa statistic: Blomberg et al. Evol. (2003)
\[\small K=\frac{\left(\mathbf{Y}-E(\mathbf{Y})\right)^{T} \left(\mathbf{Y}-E(\mathbf{Y})\right)} {\left(\mathbf{Y}-E(\mathbf{Y})\right)^{T}\mathbf{V}^{-1} \left(\mathbf{Y}-E(\mathbf{Y})\right)} / \frac{tr(\mathbf{V})-N(\mathbf{1}^{T}\mathbf{V1})^{-1}}{N-1}\]
Univariate response data only!
\(\small{K}_{mult}\): A multivariate generalization of Kappa (Adams. Syst. Biol. 2014)
Based on distances between objects in original and phylo-transformed dataspaces
\[\small\ K_{mult}= \frac{\mathbf{D}_{\mathbf{Y},E(\mathbf{Y})}^{T}\mathbf{D}_{\mathbf{Y},E(\mathbf{Y})}} {\mathbf{PD}_{\tilde{\mathbf{Y}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Y}},0}} / \frac{tr(\mathbf{V})-N(\mathbf{1}^{T}\mathbf{V1})^{-1}}{N-1}\]
\(\small{K}_{mult}\): A multivariate generalization of Kappa (Adams. Syst. Biol. 2014)
Based on distances between objects in original and phylo-transformed dataspaces
\[\small\ K_{mult}= \frac{\mathbf{D}_{\mathbf{Y},E(\mathbf{Y})}^{T}\mathbf{D}_{\mathbf{Y},E(\mathbf{Y})}} {\mathbf{PD}_{\tilde{\mathbf{Y}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Y}},0}} / \frac{tr(\mathbf{V})-N(\mathbf{1}^{T}\mathbf{V1})^{-1}}{N-1}\]
Where \(\small\mathbf{D}_{\mathbf{Y},E(\mathbf{Y})}\) is the distance from each object to the phylogenetic mean, and \(\small\mathbf{PD}_{\tilde{\mathbf{Y}},0}\) is the distance between each object in the phylogenetically-transformed space (found as: \(\small\tilde{\mathbf{Y}}=\mathbf{P}(\mathbf{Y-E(\mathbf{Y})})\)) and the origin
Do body dimensions in Plethodon display phylogenetic signal?
##
## Call:
## physignal(A = Y, phy = plethtree, iter = 999, print.progress = FALSE)
##
##
##
## Observed Phylogenetic Signal (K): 0.4342
##
## P-value: 0.002
##
## Effect Size: 2.6861
##
## Based on 1000 random permutations
Yes they do.
What evolutionary model best describes trait variation?
Fit data to phylogeny under differing evolutionary models
What evolutionary model best describes trait variation?
Fit data to phylogeny under differing evolutionary models
Model comparisons of:
1: Evolutionary rate models: BM1, BMM, etc.
2: Evolutionary ‘mode’ models: BM1, OU1, OUM, etc.
“How fast, as a matter of fact, do animals evolve in nature?” (Simpson, 1944)
Is there evidence that one trait evolves faster than another?
Obtain \(\small\mathbf{R}_0\) and \(\small{logL}\):
\[\tiny\mathbf{R}=\frac{(\mathbf{Y}-E(\mathbf{Y}))^{T}\mathbf{C}^{-1}(\mathbf{Y}-E(\mathbf{Y}))}{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)\]
\[\tiny\log{L}=-\frac{np}{2}ln{2\pi}-\frac{n}{2}ln{\begin{vmatrix}{(\mathbf{R\otimes{C}})}\end{vmatrix}}-\frac{1}{2}tr{(\mathbf{Y}-E(\mathbf{Y}))^{T}(\mathbf{R\otimes{C}})^{-1}(\mathbf{Y}-E(\mathbf{Y}))}\]
Is there evidence that one trait evolves faster than another?
Obtain \(\small\mathbf{R}_0\) and \(\small{logL}\):
\[\tiny\mathbf{R}=\frac{(\mathbf{Y}-E(\mathbf{Y}))^{T}\mathbf{C}^{-1}(\mathbf{Y}-E(\mathbf{Y}))}{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)\]
\[\tiny\log{L}=-\frac{np}{2}ln{2\pi}-\frac{n}{2}ln{\begin{vmatrix}{(\mathbf{R\otimes{C}})}\end{vmatrix}}-\frac{1}{2}tr{(\mathbf{Y}-E(\mathbf{Y}))^{T}(\mathbf{R\otimes{C}})^{-1}(\mathbf{Y}-E(\mathbf{Y}))}\]
Estimate \(\small\mathbf{R}_C\) and \(\small{logL}_C\) where diagonal of \(\mathbf{R}\) is contrained to be equal:
\[\tiny\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}\)
Is there evidence for multiple evolutionary rates on the phylogeny?
Procedure: Fit data to phylogeny under single-rate and multi-rate models. Estimate \(\small\mathbf{R}\) and \(\small\log{L}\) & compare \(\small\log{L}\) between models (LRT, AIC, simulation, etc.)
Is there evidence for multiple evolutionary rates on the phylogeny?
Procedure: Fit data to phylogeny under single-rate and multi-rate models. Estimate \(\small\mathbf{R}\) and \(\small\log{L}\) & compare \(\small\log{L}\) between models (LRT, AIC, simulation, etc.)
Works well for univariate \(\mathbf{Y}\), but Type I error \(\small\uparrow\) with p-dimensions
\(\small\log{L_{mult}}\) not a general solution for high-dimensional data
Generalizes \(\small\sigma^{2}\) to estimate the net evolutionary rate in multivariate data
Procedure: Phylogenetic-transform data \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\)
Estimate \(\small\sigma^{2}_{mult}=\frac{\mathbf{PD}_{\tilde{\mathbf{Y}},0}^{T}\mathbf{PD}_{\tilde{\mathbf{Y}},0}}{N}\) for regimes in BMM
Obtain rate ratio: \(\small{RR}=\frac{\sigma^{2}_{mult.max}}{\sigma^{2}_{mult.min}}\); evaluate via simulation or permutation (RRPP).
Displays appropriate type I error and power (breaks Rao’s paradox)
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, etc.)
Unfortunately, ALL current implementations for multivariate OU models display high misspecification rates
More work is needed in this area