In TGuillerme/dispRity: Measuring Disparity
This vignette provides advanced details on the mini chains MCMCglmm pipeline and provides code and details to reproduce the results presented in the paper associated with this vignette.
Due to the amount of data and computational time required to reproduce the results in the paper (two CPU years and 350GB of data!) we only demonstrate the method on the order Charadriiformes (shorebirds).
Important note on the structure of this vignette
This vignette is intended to provide general details on mcmcmcglmmm that will be written in plain text (like this sentence).
## This vignette also contains reproducible examples written in R snippets
this_is_a_reproducible_R_snippet <- TRUE
And finally, details about specific parameters used in the paper for repeatability will be displayed in notes like this one.
Installation and requirements
## Installing the dispRity package from CRAN
install.packages("dispRity")
## Loading the packages
library(dispRity)
Projection analyses
With the variance-covariance matrices obtained from the mcmcmcglmmm analyses, we can use linear algebra to make a mathematical interpretation of @endler2005's elaboration and innovation patterns (Fig 1. in @endler2005).
In fact, if we can interpret the variance-covariance matrix as estimated evolutionary space in the trait-space that contains all the observed traits of a group and takes evolutionary history into account, we can then use the first eigenvector of that variance-covariance matrix as the major evolutionary axis, i.e. the axis in multidimensional space along which most trait values co-evolve.
In other words, if the variance-covariance matrix is a multidimensional ellipse, its major axis is its major evolutionary axis.
Then, if species evolve parallel to this axis, they are elaborating on trait combinations; conversely, if species evolve orthogonaly to this axis, they are innovating new traits combinations (innovating).
% NC: Might need to re-write the above when GT finishes re-writing these sections of the paper/when we decide on our terminology.
Mathematically we can use linear algebra to project and reject each species in the trait-space along this major axis of evolution Figure 1.
The detailed operations are described in the supplementary materials.
source("test.rotation.R")
Figure 1: {#fig1}
Consider a space with five elements named \"a\" to \"e\" (in grey) and one focal axis in this space defined as the vector from point1: \"d\" to point2: \"e\" (the thin solid grey line linking the grey \"d\" and \"e\"). We can then define all points from \"d\" (point1) as vectors (the dashed grey lines). We can then rotate the whole space so that the vector \"d\" to \"e\" becomes the reference coordinate vector (0,1) (the thick grey line) and rotate and translate all the other vectors accordingly (the solid grey lines linking the coloured letters \"a\" to \"c\"). We can then project each of these vectors on the reference vector and measure their distance on the reference vector (the projection; the coloured solid lines) and their distance from the reference vector (the rejection; the coloured dotted lines). We can then also measure the angles between each vector and their respective projections (the angles between a grey solid line to a coloured letter and the corresponding coloured vector). Note that if this angle is equal to 90 degrees, the projection is equal to 0 (point \"c\") or if that angle is > 90 degrees, the projection is negative (point \"b\").
In practice, we can use these projections and rejections of vectors or points onto a specific other vector as a way to measure the elaboration and innovation of species within a group or of entire groups.
The first instance, measuring the elaboration and innovation of species within a group, is rather straightforward from Figure 2: the idea is to project each point (defined as the dashed vectors in Figure 2) onto a vector of interest (e.g. either the group's own main axis of variance-covariance or the entire phylogeny's one).
We measured innovation and elaboration using two main approaches (see Figure 2):
- By clades (i.e innovation/elaboration between clades), where we calculated the projections of the major axis of each clade onto the global phylogenetic major axis.
- By species (i.e innovation/elaboration within clades), where we calculated the projections of each species onto 1) the global phylogenetic major axis and 2) the phylogenetic major axis of their respective clade.
source("../../R/covar.utilities_fun.R")
## Get projection lines (from a centred matrix)
get.proj.line <- function(point, base) {
point <- t(point)
base <- t(base)
proj_x <- ( (point[1]*base[1] + point[2]*base[2]) / (base[1]*base[1] + base[2]*base[2]) ) * base[1];
proj_y <- ( (point[1]*base[1] + point[2]*base[2]) / (base[1]*base[1] + base[2]*base[2]) ) * base[2];
return(rbind(cbind(proj_x, proj_y), point))
}
set.seed(42)
## Space plot
cor_matrix <- matrix(cbind(1, 0.8, 0.8, 1), nrow = 2)
space <- space.maker(10, 2, rnorm, cor.matrix = cor_matrix)
lim <- c(floor(range(space)[1]), ceiling(range(space)[2]))
par(mfrow = c(2, 2))
## Plot original space
plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "black", main = "Original space")
abline(v = 0, col = "grey", lwd = 0.5)
abline(h = 0, col = "grey", lwd = 0.5)
## Plotting the ellipse
lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "blue", lwd = 1)
## Plotting the major axis
lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "blue", lty = 1)
## Plotting a phylogenetic major axis
lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "orange")
# legend("topright", pch = c(19, NA, NA, NA), lty = c(NA, 1,1,1), lwd = c(NA, 1, 3, 3), legend = c("Group of species", "95% CI ellipse for the group", "Major axis for the group", "Major axis for the phylogeny"), col = c("black", "blue", "blue", "orange"))
legend("bottomright", legend = "A")
## Projections of the groups
## Translate the major axis
plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "grey", main = "Groups' projection\non the phylogeny")
abline(v = 0, col = "grey", lwd = 0.5)
abline(h = 0, col = "grey", lwd = 0.5)
## Plotting the ellipse
lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "grey", lwd = 1)
## Plotting the major axis
lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "grey", lty = 1)
lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "orange")
legend("bottomright", legend = "B")
## Do the translation
major_axis <- get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75)
point_to_translate <- c(-1.2, 0)
# lines(rbind(major_axis[2,], point_to_translate), lty = 3, col = "grey")
arrows(x0 = major_axis[2,1], y0 = major_axis[2,2], x1 = point_to_translate[1], y1 = point_to_translate[2], col = "black", lty = 3, length = 0.1)
move_x <- point_to_translate[1] - major_axis[2, 1]
move_y <- point_to_translate[2] - major_axis[2, 2]
axis_trans <- major_axis
axis_trans[,1] <- axis_trans[,1] + move_x
axis_trans[,2] <- axis_trans[,2] + move_y
lines(axis_trans, lwd = 3, col = "blue", lty = 1)
## Project the translation
matrix <- axis_trans
point1 <- c(-1.2, 0)
point2 <- c(2, 0.75)
base_vector <- rbind(point1, point2)
## Get all the space (with the two last rows being the base vectors)
space_cent <- rbind(matrix, base_vector)
## Centre all the space
space_cent <- space_cent - rep(point1, rep.int(nrow(space_cent), ncol(space_cent)))
## Re-attribute the centred variables
matrix_cent <- space_cent[1:nrow(matrix), ]
base_vector_cent <- space_cent[-c(1:nrow(matrix)), ]
## Projected translation
cent_proj <- get.proj.line(matrix_cent[1,], base_vector_cent[2, ])
cent_proj[, 1] <- (cent_proj[, 1] + point1[1])
lines(cent_proj, col = "darkgreen")
## Extend the phylogeny line
extend_base <- rbind(base_vector[1, ], cent_proj[1, ])
lines(extend_base, lwd = 1, col = "orange", lty = 2)
## Projections per species
## Translate the major axis
plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "grey", main = "Species' projection\non the phylogeny")
abline(v = 0, col = "grey", lwd = 0.5)
abline(h = 0, col = "grey", lwd = 0.5)
## Plotting the ellipse
lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "grey", lwd = 1)
## Plotting the major axis
lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "grey", lty = 1)
lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "orange")
legend("bottomright", legend = "C")
## Do the projection
matrix <- space
point1 <- c(-1.2, 0)
point2 <- c(2, 0.75)
base_vector <- rbind(point1, point2)
## Get all the space (with the two last rows being the base vectors)
space_cent <- rbind(matrix, base_vector)
## Centre all the space
space_cent <- space_cent - rep(point1, rep.int(nrow(space_cent), ncol(space_cent)))
## Re-attribute the centred variables
matrix_cent <- space_cent[1:nrow(matrix), ]
base_vector_cent <- space_cent[-c(1:nrow(matrix)), ]
## Plot the projections
for(i in 1:nrow(matrix)) {
cent_proj <- get.proj.line(matrix_cent[i,], base_vector_cent[2, ])
cent_proj[, 1] <- (cent_proj[, 1] + point1[1])
lines(cent_proj, col = "darkgreen")
}
## Species projection onto their own axis
plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "grey", main = "Species' projection\non their own axis")
abline(v = 0, col = "grey", lwd = 0.5)
abline(h = 0, col = "grey", lwd = 0.5)
lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "grey")
## Plotting the ellipse
lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "blue", lwd = 1)
## Plotting the major axis
lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "blue", lty = 1)
legend("bottomright", legend = "D")
## Do the projection
matrix <- space
base_vector <- get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75)
point1 <- base_vector[1,]
point2 <- base_vector[2,]
## Get all the space (with the two last rows being the base vectors)
space_cent <- rbind(matrix, base_vector)
## Centre all the space
space_cent <- space_cent - rep(point1, rep.int(nrow(space_cent), ncol(space_cent)))
## Re-attribute the centred variables
matrix_cent <- space_cent[1:nrow(matrix), ]
base_vector_cent <- space_cent[-c(1:nrow(matrix)), ]
## Extend the base vector a little bit
extend_base <- base_vector
extend_base[1, ] <- extend_base[1,1] + 0.05
lines(extend_base, lwd = 1, col = "blue", lty = 2)
## Plot the projections
for(i in 1:nrow(matrix_cent)) {
cent_proj <- get.proj.line(matrix_cent[i,], base_vector_cent[2, ])
cent_proj[, 1] <- (cent_proj[, 1] + point1[1])
cent_proj[, 2] <- (cent_proj[, 2] + point1[2])
lines(cent_proj, col = "darkgreen")
}
Figure 2: {#fig2}
Different types of elaboration/innovation analyses for a group of species. A) Given a group of 10 species (in black) with two traits, we can calculate the 95% confidence interval variance-covariance ellipse and its major axis (in blue). We also calculate another vector (in orange), for example, the major axis of another variance-covariance matrix (e.g. from another group or the whole phylogeny). B) We can calculate the projection of the group's major axis on the other major axis by translating it (black dotted arrow) and calculating its elaboration relative to the orange vector (orange dashed line) or its relative innovation (green line). C) We can also calculate the elaboration and innovation scores for each individual species relative to the orange major axis (e.g. the phylogenetic major axis); D) or relative to their own major axis.
Note that the examples depicted in Figure 2 are not scaled for clarity.
In the algorithm implemented in dispRity, the results are always effectively scaled as described in Figure 1 and in the detailed algorithm.
In practice, this makes the vector onto which elements are projected of length one.
Furthermore, because the projections happen based on the first point of the major axis of interest (which can be arbitrary), we advise that you always centre the elaboration results on 0.5 and rescale to 1 (effectively removing 0.5 from the scaled projection value and rescaling it to 1) and to make them absolute.
This makes the elaboration and innovation results more intuitive to interpret: an elaboration score comprised in $[0, 1]$ is distributed onto the major axis, a score of $[1, \infty)$ is beyond the major axis (in any n direction).
If the major axis is, for example, outside the 95% confidence interval from the phylogenetic variance-covariance matrix, this denotes a species that is especially elaborative.
Similarly for innovation, any innovation score comprised in $[0, 1]$ indicates an innovation that is within the major axis' variance-covariance 95% confidence interval; a score of $[1, \infty)$ denotes a particularly innovative species.
Reproducible example
Here we illustrate the elaboration and innovation analysis pipeline by applying it to a data set from @cooney2017 focused on the bird order Charadriiformes (shorebirds).
Data
This order contains 359 species divided into three clades: gulls, sandpipers, and plovers that have respectively 159, 102 and 98 species each. %NC: Check numbers here as you had sandpipers twice so I'm not sure which was the number for them and which was for plovers.
For each of these species, we used the 3D beak dataset from @cooney2017 that is an ordination of the shapes of the beaks of 8748 birds (the resulting shape space used here does not contain any information about the beak size - i.e. centroid size).
We extracted the 359 charadriiform species from this space where the first three dimensions account for more than 99% of the variance in the dataset.
For the rest of the example analysis here, we used these 359 species with three dimensions as a trait-space of (hereafter the shapespace).
## Loading the Charadriiformes data
data(charadriiformes)
## Extracting the tree
tree <- charadriiformes$tree
## Extracting the data column that contains the clade assignments
data <- charadriiformes$data[, "clade"]
## Changing the levels names (the clade names) to colours
levels(data) <- c("orange", "blue", "darkgreen")
data <- as.character(data)
## Matching the data rownames to the tip order in the tree
data <- data[match(ladderize(tree)$tip.label, rownames(charadriiformes$data))]
## Matching the tip colours (labels) to their descending edges in the tree
## (and making the non-match edges grey)
clade_edges <- match.tip.edge(data, tree, replace.na = "grey")
## Plotting the results
plot(ladderize(tree), show.tip.label = FALSE, edge.color = clade_edges, main = "Charadriiformes")
legend("bottomleft", lty = c(1,1,1), col = c("blue", "darkgreen", "orange"), legend = c("plovers", "sandpipers", "gulls"))
axisPhylo()
Figure 3: {#fig3}
The Charadriiformes data used in this pipeline illustration. The axis units are millions of years ago.
Projections analysis
Using the MCMCglmm.subsets we can visualise the variance-covariance matrices calculated previously for each group in different ways (Figure 4).
## Creating a dispRity object containing the covar data
covar_data <- MCMCglmm.subsets(data = charadriiformes$data,
posteriors = charadriiformes$posteriors,
group = MCMCglmm.levels(charadriiformes$posteriors)[1:4],
rename.groups = c("gulls", "plovers", "sandpipers", "phylogeny"))
par(mfrow = c(2,2))
covar.plot(covar_data, points = TRUE, points.cex = 0.5, col = c("orange", "blue", "darkgreen"))
legend("topleft", legend = "A")
covar.plot(covar_data, points = TRUE, points.cex = 0.5, ellipses = mean, major.axes = mean, legend = TRUE, col = c("orange", "blue", "darkgreen", "grey"), apply.to.VCV = TRUE, legend.x = "bottomright", scale = "phylogeny")
legend("topleft", legend = "B")
covar.plot(covar_data, points = FALSE, ellipses = TRUE, n = 100, col = c("orange", "blue", "darkgreen", "grey"))
legend("topleft", legend = "C")
covar.plot(covar_data, points = FALSE, major.axes = TRUE, n = 100, col = c("orange", "blue", "darkgreen", "grey"))
legend("topleft", legend = "D")
Figure 4: {#fig4}
Visualisation of the posteriors of the model: A) the first two dimensions of the input trait-space; B) the same space with the average ellipses and average major axis scaled to the size of the phylogenetic (grey) ellipse; C) the distribution of 100 random ellipses (unscaled); D) the distribution of 100 random major axes (unscaled).
We can then use this "dispRity" object to calculate which of the three Charadriiformes clades is the most elaborative and innovative (Figure 2B) and then the distribution of species absolute (Figure 2C) and relative (Figure 2D) elaboration and innovation scores on the 1000 posterior samples.
In THE_PAPER we calculated the results based on the 4000 posterior samples.
This can be done automatically using the dispRity.covar.projections function:
## Measure the group differences compared to the phylogeny
group_differences <- dispRity.covar.projections(covar_data,
type = "groups",
base = "phylogeny",
output = c("position", "distance"))
## Measure the individual absolute differences
species_differences_abs <- dispRity.covar.projections(covar_data,
type = "elements",
base = "phylogeny",
output = c("position", "distance"))
## Measure the individual relative differences
species_differences_rel <- dispRity.covar.projections(covar_data,
type = "elements",
output = c("position", "distance"))
## Plotting the group_differences
par(mfrow = c(1, 2))
plot(group_differences[[1]], col = c("orange", "blue", "darkgreen", "grey"), ylab = c("elaboration"), xlab = "")
plot(group_differences[[2]], col = c("orange", "blue", "darkgreen", "grey"), ylab = c("innovation"), xlab = "")
## Creating the colouring vector
colour_vector <- c(rep("orange", size.subsets(species_differences_abs[[1]])[1]),
rep("blue", size.subsets(species_differences_abs[[1]])[2]),
rep("darkgreen", size.subsets(species_differences_abs[[1]])[3]))
## Plotting the group correlations
par(mfrow = c(1, 2))
plot(species_differences_abs, specific.args = list(correlation.plot = c("position", "distance")),
col = colour_vector, ylab = "innovation", xlab = "elaboration", pch = 19, main = "Absolute correlation")
plot(species_differences_rel, specific.args = list(correlation.plot = c("position", "distance")),
col = colour_vector, ylab = "innovation", xlab = "elaboration", pch = 19, main = "Relative correlation")
Supplementary details {#projections-supplementary}
Detailed projection and rejection operations
In more detail, we can define the major axis from the variance-covariance matrices and then project and reject each element of interest in the space onto this axis.
The following steps are generalised to $n$ dimensions and detailed below (as well as the algorithm used in dispRity to perform the transformations):
Major axis definition
For $O_{n}$, the unit hypersphere matrix of $n$ dimensions and a radius composed of the two identity matrices $I_{n}$ and $-I_{n}$ so that:
\begin{equation}
O_{n} =
\begin{pmatrix}
1 & 0 & \cdots & 0 \
0 & 1 & \cdots & 0 \
\vdots & \vdots & \ddots & \vdots \
0 & 0 & \cdots & 1 \
-1 & 0 & \cdots & 0 \
0 & -1 & \cdots & 0 \
\vdots & \vdots & \ddots & \vdots \
0 & 0 & \cdots & -1 \
\end{pmatrix}
\end{equation}
In other words, $O_{n}$ is the matrix representing the hypersphere of $n$ dimensions and of radius $1$ that fits in the centre of the trait-space;
And $O'_{n}$ is the scaled matrix hypersphere to the 95\% confidence interval size using the $\chi^2$ distribution:
$$O'{n} = O{n} \sqrt{\chi^2(0.95)}$$
Then, for the variance-covariance matrix $VCV_{n}$ of $n$ dimensions obtained from the posterior distribution of the mcmcmcglmmmm:
\begin{equation}
VCV_{n} =
\begin{pmatrix}
\sigma(a) & \sigma(a,b) & \cdots & \sigma(a,n) \
\sigma(a,b) & \sigma(b) & \cdots & \sigma(b,n) \
\vdots & \vdots & \ddots & \vdots \
\sigma(n,a) & \sigma(n,b) & \cdots & \sigma(n) \
\end{pmatrix}
\end{equation}
and the eigenvectors v and the eigenvalues $\lambda$ satisfying the following eigen decomposition:
$$VCV_{n} v = \lambda v$$
We can get $M_{n}$, the matrix containing all the edge coordinates of the 0.95 CI hypersphere from $VCV_{n}$ using the transposition of the cross product between the eigenvectors v and the product of the scaled 0.95 CI unit sphere $O'_{n}$ and the eigenvalues $\lambda$:
$$M_{n} = [(O'_{n}\sqrt{\lambda}) \times v]^{\text{T}}$$
Finally, we can centre the matrix $M_{n}$ on the estimated solution of each GLMM (model$Sol, in MCMCglmm) corresponding to the estimation of the position of the variance-covariance matrix in the trait-space.
$M_{n}$ then contains all the major axes of the 0.95 hyper-ellipse fitting the variance-covariance matrix.
We can then define the first row of the matrix, $M_{1,n}$, as the major axis, the second row, $M_{2,n}$, as the second major axis (the minor axis in a 2D ellipse), etc.
The detailed procedure was adapted from 李哲源's post on Stack Overflow and implemented in dispRity::axis.covar.
Measuring projection and rejection
Once we have defined a major axis, we can project any elements in the trait-space onto that axis.
Specifically, for any elements in the trait-space, we can define it as the vector $\vec{a}$ with one set of coordinates in $n$ dimensions:
\begin{equation}
\vec{a} =
\begin{bmatrix}
x \
y \
\cdots \
n \
\end{bmatrix}
\end{equation}
And for any major axis that we can define as a vector $\vec{b}$ as a set of pairs of coordinates in $n$ dimensions:
\begin{equation}
\vec{b} =
\begin{bmatrix}
x_{1} & x_{2} \
y_{1} & y_{2} \
\cdots & \cdots \
n_{1} & n_{2} \
\end{bmatrix}
\end{equation}
We can then calculate $\vec{a_{1}}$, the orthogonal projection of $\vec{a}$ onto $\vec{b}$ using:
\begin{equation}
\vec{a_{1}} = \frac{\vec{a} \cdot \vec{b}}{\|\vec{b}\|}
\end{equation}
With $\|\vec{b}\| = \sqrt{\vec{b} \cdot \vec{b}}$ being the norm of $\vec{b}$.
And $\vec{a_{2}}$, the rejection of $\vec{a}$ onto $\vec{b}$:
\begin{equation}
\vec{a_{2}} = \vec{a} - \vec{a_{1}}
\end{equation}
Generalisation of projection onto any vector in a set space
Using this, we can generalise the procedure so as to calculate the projection and rejection for any element within a trait-space $TS_{m,n}$:
\begin{equation}
TS_{m,n} =
\begin{bmatrix}
x_{1} & x_{2} & \cdots & x_{m} \
y_{1} & y_{2} & \cdots & y_{m} \
\vdots & \vdots & \ddots & \vdots \
n_{1} & n_{2} & \cdots & n_{m} \
\end{bmatrix}
\end{equation}
And any major axis defined as a vector $\vec{B}$:
\begin{equation}
B =
\begin{bmatrix}
x_{1} & x_{2}\
y_{1} & y_{2}\
\vdots & \vdots \
n_{1} & n_{2} \
\end{bmatrix}
\end{equation}
By using the linear transformation $f_{\vec{B}}$ of the trait-space $TS$ moving $\vec{B}$ onto $TS$'s first axis unit vector $\vec{\hat{\imath}}$:
$$f_{\vec{B}}(TS) = \left( \frac{TS - [Bx_{1}, By_{1}, \cdots, Bn_{1}]^{\text{T}}}{\|\vec{B}\|} \right) \cdot R_{\vec{B}}$$
With $R_{\vec{B}}$ being the rotation matrix of the vector $\vec{B}$ onto $\vec{\hat{\imath}}$:
\begin{equation}
R_{\vec{B}} = I_{\vec{B}} - \vec{B}\vec{B}^\text{T} - \vec{\hat{\imath}}\vec{\hat{\imath}}^\text{T} + [\vec{B} \vec{\hat{\imath}}]
\begin{bmatrix}
cos(\theta) & -sin(\theta)\
sin(\theta) & cos(\theta)\
\end{bmatrix} [\vec{B} \vec{\hat{\imath}}]^\text{T}
\end{equation}
Where $\theta$ is:
\begin{equation}
\theta = acos \left(\frac{\vec{B} \cdot \vec{\hat{\imath}}}{\|\vec{B}\| \cdot \|\vec{\hat{\imath}}\|} \right)
\end{equation}
Or $\theta = acos (B_x)$ since both $\|\vec{B}\|$ and $\|\vec{\hat{\imath}}\|$ are equal to 1 and $\|\vec{\hat{\imath}}\|$ is the unit vector on the first axis.
Algorithm for calculating the projection/rejection of any element in a defined space
In practice we followed this procedure and applied a modification of this implementation (see @aguilera2004 for the formal generalisation of this algorithm in $n$ dimensions) using the following algorithm implemented in dispRity::projections (@guillermedisprity):
- In the trait-space, define $\vec{B}$ as the base vector (typically $\vec{B}$ is defined as the pair of coordinates from the major axis described above).
- Centre the trait-space on the origin of $\vec{B}$ so that the first set of coordinates of $\vec{B}$ are 0.
- Scale the trait-space to the norm of $\vec{B}$ so that the norm of $\vec{B}$ is now 1.
- Rotate the trait-space using the rotation matrix $R_{\vec{B}}$ to satisfy the linear transformation $\vec{B} \rightarrow \vec{\hat{\imath}}$ (with $\vec{\hat{\imath}}$ being the first unit vector of the trait-space - typically the x axis unit vector).
- Project/reject every element in the trait-space on $\vec{B}$ (that is now $\vec{\hat{\imath}}$). In practice, the first coordinate (x) of each element is now its projection onto $\vec{B}$.
References
TGuillerme/dispRity documentation built on Dec. 21, 2024, 4:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
This vignette provides advanced details on the mini chains MCMCglmm pipeline and provides code and details to reproduce the results presented in the paper associated with this vignette.
Due to the amount of data and computational time required to reproduce the results in the paper (two CPU years and 350GB of data!) we only demonstrate the method on the order Charadriiformes (shorebirds).
Important note on the structure of this vignette
This vignette is intended to provide general details on mcmcmcglmmm that will be written in plain text (like this sentence).
## This vignette also contains reproducible examples written in R snippets this_is_a_reproducible_R_snippet <- TRUE
And finally, details about specific parameters used in the paper for repeatability will be displayed in notes like this one.
Installation and requirements
## Installing the dispRity package from CRAN install.packages("dispRity")
## Loading the packages library(dispRity)
Projection analyses
With the variance-covariance matrices obtained from the mcmcmcglmmm analyses, we can use linear algebra to make a mathematical interpretation of @endler2005's elaboration and innovation patterns (Fig 1. in @endler2005).
In fact, if we can interpret the variance-covariance matrix as estimated evolutionary space in the trait-space that contains all the observed traits of a group and takes evolutionary history into account, we can then use the first eigenvector of that variance-covariance matrix as the major evolutionary axis, i.e. the axis in multidimensional space along which most trait values co-evolve.
In other words, if the variance-covariance matrix is a multidimensional ellipse, its major axis is its major evolutionary axis.
Then, if species evolve parallel to this axis, they are elaborating on trait combinations; conversely, if species evolve orthogonaly to this axis, they are innovating new traits combinations (innovating).
% NC: Might need to re-write the above when GT finishes re-writing these sections of the paper/when we decide on our terminology.
Mathematically we can use linear algebra to project and reject each species in the trait-space along this major axis of evolution Figure 1. The detailed operations are described in the supplementary materials.
source("test.rotation.R")
Figure 1: {#fig1}
Consider a space with five elements named \"a\" to \"e\" (in grey) and one focal axis in this space defined as the vector from point1: \"d\" to point2: \"e\" (the thin solid grey line linking the grey \"d\" and \"e\"). We can then define all points from \"d\" (point1) as vectors (the dashed grey lines). We can then rotate the whole space so that the vector \"d\" to \"e\" becomes the reference coordinate vector (0,1) (the thick grey line) and rotate and translate all the other vectors accordingly (the solid grey lines linking the coloured letters \"a\" to \"c\"). We can then project each of these vectors on the reference vector and measure their distance on the reference vector (the projection; the coloured solid lines) and their distance from the reference vector (the rejection; the coloured dotted lines). We can then also measure the angles between each vector and their respective projections (the angles between a grey solid line to a coloured letter and the corresponding coloured vector). Note that if this angle is equal to 90 degrees, the projection is equal to 0 (point \"c\") or if that angle is > 90 degrees, the projection is negative (point \"b\").
In practice, we can use these projections and rejections of vectors or points onto a specific other vector as a way to measure the elaboration and innovation of species within a group or of entire groups. The first instance, measuring the elaboration and innovation of species within a group, is rather straightforward from Figure 2: the idea is to project each point (defined as the dashed vectors in Figure 2) onto a vector of interest (e.g. either the group's own main axis of variance-covariance or the entire phylogeny's one).
We measured innovation and elaboration using two main approaches (see Figure 2):
- By clades (i.e innovation/elaboration between clades), where we calculated the projections of the major axis of each clade onto the global phylogenetic major axis.
- By species (i.e innovation/elaboration within clades), where we calculated the projections of each species onto 1) the global phylogenetic major axis and 2) the phylogenetic major axis of their respective clade.
source("../../R/covar.utilities_fun.R") ## Get projection lines (from a centred matrix) get.proj.line <- function(point, base) { point <- t(point) base <- t(base) proj_x <- ( (point[1]*base[1] + point[2]*base[2]) / (base[1]*base[1] + base[2]*base[2]) ) * base[1]; proj_y <- ( (point[1]*base[1] + point[2]*base[2]) / (base[1]*base[1] + base[2]*base[2]) ) * base[2]; return(rbind(cbind(proj_x, proj_y), point)) } set.seed(42) ## Space plot cor_matrix <- matrix(cbind(1, 0.8, 0.8, 1), nrow = 2) space <- space.maker(10, 2, rnorm, cor.matrix = cor_matrix) lim <- c(floor(range(space)[1]), ceiling(range(space)[2])) par(mfrow = c(2, 2)) ## Plot original space plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "black", main = "Original space") abline(v = 0, col = "grey", lwd = 0.5) abline(h = 0, col = "grey", lwd = 0.5) ## Plotting the ellipse lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "blue", lwd = 1) ## Plotting the major axis lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "blue", lty = 1) ## Plotting a phylogenetic major axis lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "orange") # legend("topright", pch = c(19, NA, NA, NA), lty = c(NA, 1,1,1), lwd = c(NA, 1, 3, 3), legend = c("Group of species", "95% CI ellipse for the group", "Major axis for the group", "Major axis for the phylogeny"), col = c("black", "blue", "blue", "orange")) legend("bottomright", legend = "A") ## Projections of the groups ## Translate the major axis plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "grey", main = "Groups' projection\non the phylogeny") abline(v = 0, col = "grey", lwd = 0.5) abline(h = 0, col = "grey", lwd = 0.5) ## Plotting the ellipse lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "grey", lwd = 1) ## Plotting the major axis lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "grey", lty = 1) lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "orange") legend("bottomright", legend = "B") ## Do the translation major_axis <- get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75) point_to_translate <- c(-1.2, 0) # lines(rbind(major_axis[2,], point_to_translate), lty = 3, col = "grey") arrows(x0 = major_axis[2,1], y0 = major_axis[2,2], x1 = point_to_translate[1], y1 = point_to_translate[2], col = "black", lty = 3, length = 0.1) move_x <- point_to_translate[1] - major_axis[2, 1] move_y <- point_to_translate[2] - major_axis[2, 2] axis_trans <- major_axis axis_trans[,1] <- axis_trans[,1] + move_x axis_trans[,2] <- axis_trans[,2] + move_y lines(axis_trans, lwd = 3, col = "blue", lty = 1) ## Project the translation matrix <- axis_trans point1 <- c(-1.2, 0) point2 <- c(2, 0.75) base_vector <- rbind(point1, point2) ## Get all the space (with the two last rows being the base vectors) space_cent <- rbind(matrix, base_vector) ## Centre all the space space_cent <- space_cent - rep(point1, rep.int(nrow(space_cent), ncol(space_cent))) ## Re-attribute the centred variables matrix_cent <- space_cent[1:nrow(matrix), ] base_vector_cent <- space_cent[-c(1:nrow(matrix)), ] ## Projected translation cent_proj <- get.proj.line(matrix_cent[1,], base_vector_cent[2, ]) cent_proj[, 1] <- (cent_proj[, 1] + point1[1]) lines(cent_proj, col = "darkgreen") ## Extend the phylogeny line extend_base <- rbind(base_vector[1, ], cent_proj[1, ]) lines(extend_base, lwd = 1, col = "orange", lty = 2) ## Projections per species ## Translate the major axis plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "grey", main = "Species' projection\non the phylogeny") abline(v = 0, col = "grey", lwd = 0.5) abline(h = 0, col = "grey", lwd = 0.5) ## Plotting the ellipse lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "grey", lwd = 1) ## Plotting the major axis lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "grey", lty = 1) lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "orange") legend("bottomright", legend = "C") ## Do the projection matrix <- space point1 <- c(-1.2, 0) point2 <- c(2, 0.75) base_vector <- rbind(point1, point2) ## Get all the space (with the two last rows being the base vectors) space_cent <- rbind(matrix, base_vector) ## Centre all the space space_cent <- space_cent - rep(point1, rep.int(nrow(space_cent), ncol(space_cent))) ## Re-attribute the centred variables matrix_cent <- space_cent[1:nrow(matrix), ] base_vector_cent <- space_cent[-c(1:nrow(matrix)), ] ## Plot the projections for(i in 1:nrow(matrix)) { cent_proj <- get.proj.line(matrix_cent[i,], base_vector_cent[2, ]) cent_proj[, 1] <- (cent_proj[, 1] + point1[1]) lines(cent_proj, col = "darkgreen") } ## Species projection onto their own axis plot(space, pch = 19, xlab = "Trait 1", ylab = "Trait 2", xlim = lim, ylim = lim, col = "grey", main = "Species' projection\non their own axis") abline(v = 0, col = "grey", lwd = 0.5) abline(h = 0, col = "grey", lwd = 0.5) lines(x = c(-1.2, 2), y = c(0, 0.75), lwd = 3, col = "grey") ## Plotting the ellipse lines(ellipse::ellipse(cor_matrix, level = 0.75), col = "blue", lwd = 1) ## Plotting the major axis lines(get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75), lwd = 3, col = "blue", lty = 1) legend("bottomright", legend = "D") ## Do the projection matrix <- space base_vector <- get.one.axis(list(VCV = cor_matrix, Sol = c(0,0)), level = 0.75) point1 <- base_vector[1,] point2 <- base_vector[2,] ## Get all the space (with the two last rows being the base vectors) space_cent <- rbind(matrix, base_vector) ## Centre all the space space_cent <- space_cent - rep(point1, rep.int(nrow(space_cent), ncol(space_cent))) ## Re-attribute the centred variables matrix_cent <- space_cent[1:nrow(matrix), ] base_vector_cent <- space_cent[-c(1:nrow(matrix)), ] ## Extend the base vector a little bit extend_base <- base_vector extend_base[1, ] <- extend_base[1,1] + 0.05 lines(extend_base, lwd = 1, col = "blue", lty = 2) ## Plot the projections for(i in 1:nrow(matrix_cent)) { cent_proj <- get.proj.line(matrix_cent[i,], base_vector_cent[2, ]) cent_proj[, 1] <- (cent_proj[, 1] + point1[1]) cent_proj[, 2] <- (cent_proj[, 2] + point1[2]) lines(cent_proj, col = "darkgreen") }
Figure 2: {#fig2}
Different types of elaboration/innovation analyses for a group of species. A) Given a group of 10 species (in black) with two traits, we can calculate the 95% confidence interval variance-covariance ellipse and its major axis (in blue). We also calculate another vector (in orange), for example, the major axis of another variance-covariance matrix (e.g. from another group or the whole phylogeny). B) We can calculate the projection of the group's major axis on the other major axis by translating it (black dotted arrow) and calculating its elaboration relative to the orange vector (orange dashed line) or its relative innovation (green line). C) We can also calculate the elaboration and innovation scores for each individual species relative to the orange major axis (e.g. the phylogenetic major axis); D) or relative to their own major axis.
Note that the examples depicted in Figure 2 are not scaled for clarity.
In the algorithm implemented in dispRity, the results are always effectively scaled as described in Figure 1 and in the detailed algorithm.
In practice, this makes the vector onto which elements are projected of length one.
Furthermore, because the projections happen based on the first point of the major axis of interest (which can be arbitrary), we advise that you always centre the elaboration results on 0.5 and rescale to 1 (effectively removing 0.5 from the scaled projection value and rescaling it to 1) and to make them absolute.
This makes the elaboration and innovation results more intuitive to interpret: an elaboration score comprised in $[0, 1]$ is distributed onto the major axis, a score of $[1, \infty)$ is beyond the major axis (in any n direction).
If the major axis is, for example, outside the 95% confidence interval from the phylogenetic variance-covariance matrix, this denotes a species that is especially elaborative.
Similarly for innovation, any innovation score comprised in $[0, 1]$ indicates an innovation that is within the major axis' variance-covariance 95% confidence interval; a score of $[1, \infty)$ denotes a particularly innovative species.
Reproducible example
Here we illustrate the elaboration and innovation analysis pipeline by applying it to a data set from @cooney2017 focused on the bird order Charadriiformes (shorebirds).
Data
This order contains 359 species divided into three clades: gulls, sandpipers, and plovers that have respectively 159, 102 and 98 species each. %NC: Check numbers here as you had sandpipers twice so I'm not sure which was the number for them and which was for plovers. For each of these species, we used the 3D beak dataset from @cooney2017 that is an ordination of the shapes of the beaks of 8748 birds (the resulting shape space used here does not contain any information about the beak size - i.e. centroid size). We extracted the 359 charadriiform species from this space where the first three dimensions account for more than 99% of the variance in the dataset. For the rest of the example analysis here, we used these 359 species with three dimensions as a trait-space of (hereafter the shapespace).
## Loading the Charadriiformes data data(charadriiformes)
## Extracting the tree tree <- charadriiformes$tree ## Extracting the data column that contains the clade assignments data <- charadriiformes$data[, "clade"] ## Changing the levels names (the clade names) to colours levels(data) <- c("orange", "blue", "darkgreen") data <- as.character(data) ## Matching the data rownames to the tip order in the tree data <- data[match(ladderize(tree)$tip.label, rownames(charadriiformes$data))] ## Matching the tip colours (labels) to their descending edges in the tree ## (and making the non-match edges grey) clade_edges <- match.tip.edge(data, tree, replace.na = "grey") ## Plotting the results plot(ladderize(tree), show.tip.label = FALSE, edge.color = clade_edges, main = "Charadriiformes") legend("bottomleft", lty = c(1,1,1), col = c("blue", "darkgreen", "orange"), legend = c("plovers", "sandpipers", "gulls")) axisPhylo()
Figure 3: {#fig3}
The Charadriiformes data used in this pipeline illustration. The axis units are millions of years ago.
Projections analysis
Using the MCMCglmm.subsets we can visualise the variance-covariance matrices calculated previously for each group in different ways (Figure 4).
## Creating a dispRity object containing the covar data covar_data <- MCMCglmm.subsets(data = charadriiformes$data, posteriors = charadriiformes$posteriors, group = MCMCglmm.levels(charadriiformes$posteriors)[1:4], rename.groups = c("gulls", "plovers", "sandpipers", "phylogeny"))
par(mfrow = c(2,2)) covar.plot(covar_data, points = TRUE, points.cex = 0.5, col = c("orange", "blue", "darkgreen")) legend("topleft", legend = "A") covar.plot(covar_data, points = TRUE, points.cex = 0.5, ellipses = mean, major.axes = mean, legend = TRUE, col = c("orange", "blue", "darkgreen", "grey"), apply.to.VCV = TRUE, legend.x = "bottomright", scale = "phylogeny") legend("topleft", legend = "B") covar.plot(covar_data, points = FALSE, ellipses = TRUE, n = 100, col = c("orange", "blue", "darkgreen", "grey")) legend("topleft", legend = "C") covar.plot(covar_data, points = FALSE, major.axes = TRUE, n = 100, col = c("orange", "blue", "darkgreen", "grey")) legend("topleft", legend = "D")
Figure 4: {#fig4}
Visualisation of the posteriors of the model: A) the first two dimensions of the input trait-space; B) the same space with the average ellipses and average major axis scaled to the size of the phylogenetic (grey) ellipse; C) the distribution of 100 random ellipses (unscaled); D) the distribution of 100 random major axes (unscaled).
We can then use this "dispRity" object to calculate which of the three Charadriiformes clades is the most elaborative and innovative (Figure 2B) and then the distribution of species absolute (Figure 2C) and relative (Figure 2D) elaboration and innovation scores on the 1000 posterior samples.
In THE_PAPER we calculated the results based on the 4000 posterior samples.
This can be done automatically using the dispRity.covar.projections function:
## Measure the group differences compared to the phylogeny group_differences <- dispRity.covar.projections(covar_data, type = "groups", base = "phylogeny", output = c("position", "distance")) ## Measure the individual absolute differences species_differences_abs <- dispRity.covar.projections(covar_data, type = "elements", base = "phylogeny", output = c("position", "distance")) ## Measure the individual relative differences species_differences_rel <- dispRity.covar.projections(covar_data, type = "elements", output = c("position", "distance"))
## Plotting the group_differences par(mfrow = c(1, 2)) plot(group_differences[[1]], col = c("orange", "blue", "darkgreen", "grey"), ylab = c("elaboration"), xlab = "") plot(group_differences[[2]], col = c("orange", "blue", "darkgreen", "grey"), ylab = c("innovation"), xlab = "")
## Creating the colouring vector colour_vector <- c(rep("orange", size.subsets(species_differences_abs[[1]])[1]), rep("blue", size.subsets(species_differences_abs[[1]])[2]), rep("darkgreen", size.subsets(species_differences_abs[[1]])[3])) ## Plotting the group correlations par(mfrow = c(1, 2)) plot(species_differences_abs, specific.args = list(correlation.plot = c("position", "distance")), col = colour_vector, ylab = "innovation", xlab = "elaboration", pch = 19, main = "Absolute correlation") plot(species_differences_rel, specific.args = list(correlation.plot = c("position", "distance")), col = colour_vector, ylab = "innovation", xlab = "elaboration", pch = 19, main = "Relative correlation")
Supplementary details {#projections-supplementary}
Detailed projection and rejection operations
In more detail, we can define the major axis from the variance-covariance matrices and then project and reject each element of interest in the space onto this axis.
The following steps are generalised to $n$ dimensions and detailed below (as well as the algorithm used in dispRity to perform the transformations):
Major axis definition
For $O_{n}$, the unit hypersphere matrix of $n$ dimensions and a radius composed of the two identity matrices $I_{n}$ and $-I_{n}$ so that:
\begin{equation} O_{n} = \begin{pmatrix} 1 & 0 & \cdots & 0 \ 0 & 1 & \cdots & 0 \ \vdots & \vdots & \ddots & \vdots \ 0 & 0 & \cdots & 1 \ -1 & 0 & \cdots & 0 \ 0 & -1 & \cdots & 0 \ \vdots & \vdots & \ddots & \vdots \ 0 & 0 & \cdots & -1 \ \end{pmatrix} \end{equation}
In other words, $O_{n}$ is the matrix representing the hypersphere of $n$ dimensions and of radius $1$ that fits in the centre of the trait-space;
And $O'_{n}$ is the scaled matrix hypersphere to the 95\% confidence interval size using the $\chi^2$ distribution:
$$O'{n} = O{n} \sqrt{\chi^2(0.95)}$$
Then, for the variance-covariance matrix $VCV_{n}$ of $n$ dimensions obtained from the posterior distribution of the mcmcmcglmmmm:
\begin{equation} VCV_{n} = \begin{pmatrix} \sigma(a) & \sigma(a,b) & \cdots & \sigma(a,n) \ \sigma(a,b) & \sigma(b) & \cdots & \sigma(b,n) \ \vdots & \vdots & \ddots & \vdots \ \sigma(n,a) & \sigma(n,b) & \cdots & \sigma(n) \ \end{pmatrix} \end{equation}
and the eigenvectors v and the eigenvalues $\lambda$ satisfying the following eigen decomposition:
$$VCV_{n} v = \lambda v$$
We can get $M_{n}$, the matrix containing all the edge coordinates of the 0.95 CI hypersphere from $VCV_{n}$ using the transposition of the cross product between the eigenvectors v and the product of the scaled 0.95 CI unit sphere $O'_{n}$ and the eigenvalues $\lambda$:
$$M_{n} = [(O'_{n}\sqrt{\lambda}) \times v]^{\text{T}}$$
Finally, we can centre the matrix $M_{n}$ on the estimated solution of each GLMM (model$Sol, in MCMCglmm) corresponding to the estimation of the position of the variance-covariance matrix in the trait-space.
$M_{n}$ then contains all the major axes of the 0.95 hyper-ellipse fitting the variance-covariance matrix.
We can then define the first row of the matrix, $M_{1,n}$, as the major axis, the second row, $M_{2,n}$, as the second major axis (the minor axis in a 2D ellipse), etc.
The detailed procedure was adapted from 李哲源's post on Stack Overflow and implemented in dispRity::axis.covar.
Measuring projection and rejection
Once we have defined a major axis, we can project any elements in the trait-space onto that axis. Specifically, for any elements in the trait-space, we can define it as the vector $\vec{a}$ with one set of coordinates in $n$ dimensions:
\begin{equation} \vec{a} = \begin{bmatrix} x \ y \ \cdots \ n \ \end{bmatrix} \end{equation}
And for any major axis that we can define as a vector $\vec{b}$ as a set of pairs of coordinates in $n$ dimensions:
\begin{equation} \vec{b} = \begin{bmatrix} x_{1} & x_{2} \ y_{1} & y_{2} \ \cdots & \cdots \ n_{1} & n_{2} \ \end{bmatrix} \end{equation}
We can then calculate $\vec{a_{1}}$, the orthogonal projection of $\vec{a}$ onto $\vec{b}$ using:
\begin{equation} \vec{a_{1}} = \frac{\vec{a} \cdot \vec{b}}{\|\vec{b}\|} \end{equation}
With $\|\vec{b}\| = \sqrt{\vec{b} \cdot \vec{b}}$ being the norm of $\vec{b}$. And $\vec{a_{2}}$, the rejection of $\vec{a}$ onto $\vec{b}$:
\begin{equation} \vec{a_{2}} = \vec{a} - \vec{a_{1}} \end{equation}
Generalisation of projection onto any vector in a set space
Using this, we can generalise the procedure so as to calculate the projection and rejection for any element within a trait-space $TS_{m,n}$:
\begin{equation} TS_{m,n} = \begin{bmatrix} x_{1} & x_{2} & \cdots & x_{m} \ y_{1} & y_{2} & \cdots & y_{m} \ \vdots & \vdots & \ddots & \vdots \ n_{1} & n_{2} & \cdots & n_{m} \ \end{bmatrix} \end{equation}
And any major axis defined as a vector $\vec{B}$:
\begin{equation} B = \begin{bmatrix} x_{1} & x_{2}\ y_{1} & y_{2}\ \vdots & \vdots \ n_{1} & n_{2} \ \end{bmatrix} \end{equation}
By using the linear transformation $f_{\vec{B}}$ of the trait-space $TS$ moving $\vec{B}$ onto $TS$'s first axis unit vector $\vec{\hat{\imath}}$:
$$f_{\vec{B}}(TS) = \left( \frac{TS - [Bx_{1}, By_{1}, \cdots, Bn_{1}]^{\text{T}}}{\|\vec{B}\|} \right) \cdot R_{\vec{B}}$$
With $R_{\vec{B}}$ being the rotation matrix of the vector $\vec{B}$ onto $\vec{\hat{\imath}}$:
\begin{equation} R_{\vec{B}} = I_{\vec{B}} - \vec{B}\vec{B}^\text{T} - \vec{\hat{\imath}}\vec{\hat{\imath}}^\text{T} + [\vec{B} \vec{\hat{\imath}}] \begin{bmatrix} cos(\theta) & -sin(\theta)\ sin(\theta) & cos(\theta)\ \end{bmatrix} [\vec{B} \vec{\hat{\imath}}]^\text{T} \end{equation}
Where $\theta$ is:
\begin{equation} \theta = acos \left(\frac{\vec{B} \cdot \vec{\hat{\imath}}}{\|\vec{B}\| \cdot \|\vec{\hat{\imath}}\|} \right) \end{equation}
Or $\theta = acos (B_x)$ since both $\|\vec{B}\|$ and $\|\vec{\hat{\imath}}\|$ are equal to 1 and $\|\vec{\hat{\imath}}\|$ is the unit vector on the first axis.
Algorithm for calculating the projection/rejection of any element in a defined space
In practice we followed this procedure and applied a modification of this implementation (see @aguilera2004 for the formal generalisation of this algorithm in $n$ dimensions) using the following algorithm implemented in dispRity::projections (@guillermedisprity):
- In the trait-space, define $\vec{B}$ as the base vector (typically $\vec{B}$ is defined as the pair of coordinates from the major axis described above).
- Centre the trait-space on the origin of $\vec{B}$ so that the first set of coordinates of $\vec{B}$ are 0.
- Scale the trait-space to the norm of $\vec{B}$ so that the norm of $\vec{B}$ is now 1.
- Rotate the trait-space using the rotation matrix $R_{\vec{B}}$ to satisfy the linear transformation $\vec{B} \rightarrow \vec{\hat{\imath}}$ (with $\vec{\hat{\imath}}$ being the first unit vector of the trait-space - typically the x axis unit vector).
- Project/reject every element in the trait-space on $\vec{B}$ (that is now $\vec{\hat{\imath}}$). In practice, the first coordinate (x) of each element is now its projection onto $\vec{B}$.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.