To test the variation in landmark position among different specimens and different anatomical regions, we can apply the landmark partition test described in details below. This vignette contains the theoretical description of the test followed by a step by step tutorial of the implementation of the test in landvR.

Landmark partition test

Requirements and data

The packages

This test will be using the landvR package from github and will require devtools, geomorph and dispRity from the CRAN. All these packages can be loaded easily as follows

## Loading the libraries (and installing if necessary)
if(!require(devtools)) install.packages("devtools")
if(!require(geomorph)) install.packages("geomorph")
if(!require(dispRity)) install.packages("dispRity")
if(!require(landvR)) install_github("TGuillerme/landvR")

## Setting a random seed for the repeatability of this tutorial
set.seed(42)

The data

This whole method is based on Procrustes superimpositions that can be easily done through the geomorph package (see the geomorph manual and especially the function gpagen for more details). For this example, we will be using the plethodon dataset from the geomorph R package.

## Loading the dataset
data(plethodon)

## Generating a Procrustes superimposition
procrustes <- gpagen(plethodon$land, print.progress = FALSE)

This results in a dataset of 40 Procrustes superimposition of 12 landmarks each.

Measuring the range of variation for each landmark

Landmark variation can be measured in different ways described in more details in the ?variation.range function manual. For this example, we will only consider the radius of the polar coordinates landmark difference between specimens as the range of variation for each landmark. In other words when one landmark has a big range of variation of one specimen it means that this specimen’s landmark is further away from the mean specimen equivalent landmark.

In more technical details, the range of variation is calculate as the difference in length between a pair of homologous landmarks. This is calculated using a spherical coordinate system expressed by three parameters: (1) the radius ($\rho$, the euclidean distance between the two landmarks); (2) the azimuth angle ($\phi$, the angle on the equatorial plan); and (3) the polar angle ($\theta$, the angle measured on the polar plan). Since we are only interested in the magnitude of difference ($\rho$) regardless of the direction ($\phi$ and $\theta$) we will use only the radius below.

This can be calculated in landvR easily by using the coordinates.difference function. For example, the following will measure the polar coordinates differences for all the specimens from the consensus Procrustes shape:

## Measuring the spherical coordinates differences from the mean
differences_from_mean <- coordinates.difference(coordinates = procrustes$coords,
                                               reference = procrustes$consensus, type = "spherical")

This results in a list of coordinates for each specimen. For example, to look at the differences between the first specimen and the consensus in terms of spherical coordinates, one can call:

## The differences between the first specimen and the mean in terms of spherical coordinates
differences_from_mean[[1]]

##            radius   azimuth
##  [1,] 0.034421355  1.297680
##  [2,] 0.017539211  8.094379
##  [3,] 0.012244155 83.462940
##  [4,] 0.005645058 15.700107
##  [5,] 0.004350167 12.567055
##  [6,] 0.019478593 64.112830
##  [7,] 0.024723159  3.525265
##  [8,] 0.014953078 38.475144
##  [9,] 0.006452463 44.483161
## [10,] 0.004439587 52.440277
## [11,] 0.046264826 25.611556
## [12,] 0.028748953 88.276001

Note that the angles are expressed in degrees here (column $azimuth) but they can be measured in radian using the angle option in coordinates.difference.

The radius of spherical coordinates differences for all the specimens to the mean specimen is equivalent to listing all the vectors between each mean landmark position and each specimens’ one. In small datasets like here, this can be easily visualised using the following:

## Graphical parameters
par(mar = c(0, 0, 0, 0))
## The variation range between the consensus and the first specimen
procrustes.var.plot(procrustes$consensus, procrustes$coords[, , 1])

Here the grey dots are the landmarks’ mean positions and the black arrows are the distances from these positions to the first specimen landmarks. The length of these black arrows is equal to the values in the column $radius in the table above.

PCA-based hypothesis construction

Using this system of coordinates, we can use a principal component analysis (PCA) to select the major axis of variation in terms of landmark coordinates displacement (along the first principal coordinate - PC1) to propose so biological hypothesis.

First we need to ordinate the data

## Transforming the data into a 2 dimensional array
procrustes_2d_array <- geomorph::two.d.array(procrustes$coords)

## Ordinating the data
ordination <- stats::prcomp(procrustes_2d_array)

We can then measure the difference in terms of spherical coordinates between the two extreme hypothetical specimens along the PC1. This can be done easily using the variation.range function (more on that later):

## Calculating the coordinate differences between the two extreme hypothetical specimens on PC1
PC1_variation <- variation.range(procrustes, axis = 1, ordination = ordination, type = "spherical")

## The range of hypothetical variation
PC1_variation

##            radius   azimuth
##  [1,] 0.028292228 16.270393
##  [2,] 0.041730031 25.174721
##  [3,] 0.038191480 19.448824
##  [4,] 0.003633941 39.462415
##  [5,] 0.009845285 63.332780
##  [6,] 0.012163183 61.621010
##  [7,] 0.023786856 78.288791
##  [8,] 0.011400984 16.949942
##  [9,] 0.028571133 48.377964
## [10,] 0.022508686 17.776472
## [11,] 0.114862307  3.860291
## [12,] 0.045143902 19.522473

As we will see later, this function can be used to select the actual specimens closest to this range. Nonetheless, we can use this information to plot the main variation (on PC1) on real specimens to propose some biological hypothesis based on the landmark variation. To generate the coordinates of the hypothetical specimens, we can use the plotTangentSpace function from geomorph:

## Wrapping the specimens on the tangent space
wrap_PCA <- plotTangentSpace(procrustes$coords, verbose = FALSE)

The two most extreme specimens on the first axis of the PCA are the ones plotted with on the deformed meshes.

We can then plot the coordinates of the hypothetical specimens along the first PC axis and the range of landmark variation using the procrustes.var.plot function in landvR:

## Selecting the specimens
hypothetical_1 <- wrap_PCA$pc.shapes$PC1min
hypothetical_2 <- wrap_PCA$pc.shapes$PC1max

## Plotting the range of variation using a heat colour code and the PC1 variation range
procrustes.var.plot(hypothetical_1, hypothetical_2, col = heat.colors, 
                    col.val = PC1_variation[, "radius"], labels = TRUE)

From this plot, we can hypothesise that the variation in landmark change seems to be concentrated in the the landmarks on the right hand side of the plot (landmarks 2, 11 and 12). This could be tested using the PCA data but would certainly lead to some circularity (i.e. we cannot use the PCA to determine that landmarks 2, 11 and 12 vary the most and then use the same PCA data to test it). We can thus use the observed data from the Procrustes superimposition to test it.

Note that the coordinates differences calculated from either plotTangentSpace (hypothetical_1 and hypothetical_2) and from variation.range are identical:

## Measuring the difference directly from the hypothetical shapes from plotTangentSpace
var_from_geomorph <- coordinates.difference(hypothetical_1, hypothetical_2, type = "spherical")[[1]]

## The results are equal to the ones from the variation.range function (up to the 13th digit!)
all(round(var_from_geomorph, 13) == round(PC1_variation, 13))

## [1] TRUE

Testing the hypothesis

First let’s divide our landmark into two distinct datasets.

## The "front of the skull" partition
land_front <- c(2, 11, 12)

The idea is to see whether the landmark variation is greater in the landmarks in the front of the skull and relatively smaller in the back of the skull. We can measure that using two different statistics:

The displacement difference that is the proxy of the cumulative landmark variation
The Bhattacharyya Coefficient which is the probability of overlap between two distributions

Of course, any other type of statistics can be used in the test that follows. These two however are a pretty intuitive proxy for testing our hypothesis above.

The Displacement difference

This statistic can be measured as follows:

$\Delta_{displacement} = \int_{0}^{n-1} (f_{x} - f_{y})d(x,y)$

Where n is minimum number of comparable landmarks and $f_{x}$ and $f_{y}$ are ranked functions (i.e. $f_{0} \geq f_{1} \geq f_{2} ...$) for the landmarks in the partition and all the landmarks respectively. If one of the functions $f_{x}$ or $f_{y}$ have m elements (with $m > n$) $f^{*}_{z}$, a rarefied estimated of the function with m elements is used instead.

$\int_{0}^{n-1}f^*_{z}d(z) = \frac{\sum_1^p\int f^*_{zi}}{s}$

Where s is the number of rarefaction replicates. s is chosen based on the Silverman’s rule of thumb for choosing the bandwidth of a Gaussian kernel density estimator multiplied by 1000 with a result forced to be 100 $\leq p \leq$ 1000 (Silverman 1986). This rule of thumb is implemented in R in the bw.nrd0 (a bandwith selector for Kernel Density estimations)

This statistic is implemented in landvR in the area.diff function that simply intakes two distributions:

## Two random distributions
distribution_1 <- rnorm(10)
distribution_2 <- rnorm(10)

## Their displacement difference
area.diff(distribution_1, distribution_2)

## [1] 6.195767

Probability of overlap between size differences (Bhattacharyya Coefficient)

The Bhattacharyya Coefficient calculates the probability of overlap of two distributions [Bhattacharyya (1943); Guillerme2016146]. When it is equal to zero, the probability of overlap of the distributions is also zero, and when it is equal to one, the two distributions are entirely overlapping. It forms an elegant and easy to compute continuous measurement of the probability of similarity between two distributions. The coefficient is calculated as the sum of the square root of the relative counts shared in n bins among two distributions:

$BC=\sum_{i=1}^{n} \sqrt{{\sum{a_i}}\times{\sum{b_i}}}$

Where ${a_i}$ and ${b_i}$ are the number of counts in bin i for the distributions a and b respectively divided by the total number of counts in the distribution a and b respectively. n was determined using the Silverman’s rule of thumb (see above).

This statistic is implemented in the dispRity package in the bhatt.coeff function that intakes the same arguments as area.diff:

## Their probability of overlap
bhatt.coeff(distribution_1, distribution_2)

## [1] 0.6146264

Random permutation test

Now that we have two statistics (note that one could be sufficient), we can finally test our hypothesis using a random permutation test. This test measures for on statistic between two populations (e.g. the displacement difference between the landmarks variation in the front or the back of the skull) come from the same global statistical population (H0) or different ones (H1).

First we measured the statistic between the landmark partition of interest and all the other landmarks (including the ones from the partition). Second, we generated 1000 statistics by randomly sampling the same number of landmarks as in the partition in the whole distributions and compared them again to the full distribution. This resulted in 1000 null comparisons (i.e. assuming the null hypothesis that the statistic in the partition is the same as the statistic in the total distribution). We then calculated the p value based on:

$p=\frac{\sum_1^B\left(random_{i} >= observed\right)+1}{B + 1}$

Where B is the number of random draws (1000 bootstrap pseudo-replicates), $random_{i}$ is the $i^{th}$ statistic from the comparison of the $i^{th}$ randomly drawn partition and the distribution and $observed$ is the observed statistic between the partition and the distribution.

Selecting the range of variation

First we will need to calculate the range of variation between or specimen of interest. In fact, the range of variation measured above was based on the hypothetical PC1 extremes, we thus need to have some kind of variation based on the actual Procrustes specimens to test our hypothesis. In practice, we need to select two specimens on which to apply the coordinates.difference function detailed above. We could for example decide to consider:

the range between two species using their mean shapes
the range within one species using their most different shapes
the range within one species using the 0.95 confidence interval (CI) between their most different shapes

Note that there are many other ways to select this range of variation but only these three will be detailed here:

Range between two species

For this, we can use a simple function select.procrustes that will calculate the mean shape for different subsets of a dataset. The plethodon dataset contains information on two sets of species stored in plethodon$species. We can use this data to select the Procrustes means for each species.

The select.procrustes function intakes a Procrustes dataset, a list of species and a function to use for summarising each groups (here we want the mean of each group).

## Making the list of species
species_list <- list("Jord" = which(plethodon$species == "Jord"),
                    "Teyah" = which(plethodon$species == "Teyah"))

## Selecting the procrustes means for each species
species_means <- select.procrustes(procrustes, factors = species_list)

## The list of both species means
species_means

## $Jord
##              [,1]         [,2]
##  [1,]  0.15972600 -0.018694066
##  [2,]  0.19292481 -0.088683758
##  [3,] -0.03102401 -0.002242806
##  [4,] -0.28009869 -0.094375946
##  [5,] -0.31098111 -0.062441507
##  [6,] -0.32577543 -0.037323461
##  [7,] -0.32151306  0.037915739
##  [8,] -0.18816791  0.101728921
##  [9,]  0.01989675  0.104643792
## [10,]  0.18915661  0.080485056
## [11,]  0.35711201  0.064196747
## [12,]  0.53874402 -0.085208711
## 
## $Teyah
##              [,1]        [,2]
##  [1,]  0.14493869 -0.03177919
##  [2,]  0.19365475 -0.10139889
##  [3,] -0.03637704 -0.01161697
##  [4,] -0.28354985 -0.08436582
##  [5,] -0.31047223 -0.05322464
##  [6,] -0.32622497 -0.02684087
##  [7,] -0.31363236  0.04219763
##  [8,] -0.18832063  0.09896532
##  [9,]  0.02328874  0.09306291
## [10,]  0.18977919  0.06939520
## [11,]  0.34715000  0.05883370
## [12,]  0.55976570 -0.05322839

We can then use the coordinates.difference function as detailed above to get the range of difference between both species means:

## Species means differences
species_means_difference <- coordinates.difference(species_means[[1]], species_means[[2]],
                                                   type = "spherical")

Range within one species using the most extreme shapes

When looking at variation within one single species (or a pool of two species as in the example here), we can calculate the range of landmark variation (here the two individuals with the most different radii) using the following algorithm:

Calculate the radii for all n landmarks between the mean Procrustes superimposition and each specimen’s superimposition.
Rank each set of n radii and measure the area under the resulting curve.
Select the specimen with the highest n radii area. This is now the “X” Procrustes (i.e. the specimen with the most different shape compared to the mean).
Calculate the radii for all n landmarks between the “X” Procrustes and the remaining individual ones.
Repeat step 2 and 3. The resulting specimen with the highest n radii area is now the “Y” Procrustes (i.e. the specimen with the most different shape compared to the “X” Procrustes).

These two “X” and “Y” Procrustes superimpositions are not the biggest/smallest, widest/narowest, etc. features per se but rather the two extremes of the overall distribution of landmark variation ($\rho$). Also, note that we use the radii here as an optimisation for differentiating our specimens but any other measurement can be used.

This algorithm is implemented in the variation.range function that we used a bit above (but without the ordination arguments). This function simply intakes the set of landmarks, the type of difference to measure ("spherical") and what to use for the optimisation ("radius"):

## Measuring the range of variation between the two most extreme specimens
all_spec_difference <- variation.range(procrustes, type = "spherical", what = "radius",
                                       return.ID = TRUE)

Range within one species using the 0.95 CI

One of the problems using the algorithm above can occur when the distribution contains outliers that could artificially increase this range of variation. One easy way to correct for that is to consider the same algorithm described above but rather than selecting the specimen “X” or “Y” the furthest away from respectively the mean or “X”, we can use the one the furthest away excluding the last 5% of variation. This is done easily through the CI argument in the variation.range function:

## Measuring the range of variation between the 0.95 CI extreme specimens
ci95_spec_difference <- variation.range(procrustes, type = "spherical", what = "radius", CI = 0.95,
                                       return.ID = TRUE)

Note that the resulting specimens should be different in both cases:

## The two most extreme specimens
all_spec_difference$min.max

## [1] 34 14

## The 0.95 CI extreme specimens
ci95_spec_difference$min.max

## [1]  7 27

Running the permutation test

We can use the rand.test function from landvR that’s based on a modification from the as.randtest function from the ade4 package (Thioulouse et al. 1997)). This function allows to pass the statistic as a function to the test argument as follows. We can use the argument test.parameter to calculate the p value as outlined above (note that by default, these tests are two sided).

## The random permutation test for the displacement difference between species
sp_mean_test_area <- rand.test(species_means_difference[[1]][, "radius"], subset = land_front,
                               test = area.diff, test.parameter = TRUE)
## The random permutation test for the probability of overlap between species
sp_mean_test_bc <- rand.test(species_means_difference[[1]][, "radius"], subset = land_front,
                            test = bhatt.coeff, test.parameter = TRUE)

## The random permutation test for the displacement difference within both species
all_sp_test_area <- rand.test(all_spec_difference$range[, "radius"], subset = land_front,
                              test = area.diff, test.parameter = TRUE)
## The random permutation test for the probability of overlap within both species
all_sp_test_bc <- rand.test(all_spec_difference$range[, "radius"], subset = land_front,
                            test = bhatt.coeff, test.parameter = TRUE)

## The random permutation test for the displacement difference within both species
ci95_sp_test_area <- rand.test(ci95_spec_difference$range[, "radius"], subset = land_front,
                               test = area.diff, test.parameter = TRUE)
## The random permutation test for the probability of overlap within both species
ci95_sp_test_bc <- rand.test(ci95_spec_difference$range[, "radius"], subset = land_front,
                             test = bhatt.coeff, test.parameter = TRUE)

Each test can be summarised as a S3 randtest object by simply calling them as follows:

sp_mean_test_area

## Monte-Carlo test
## Call: rand.test(distribution = species_means_difference[[1]][, "radius"], 
##     subset = land_front, test = area.diff, test.parameter = TRUE)
## 
## Observation: 0.01230886 
## 
## Based on 100 replicates
## Simulated p-value: 0.1089109 
## Alternative hypothesis: two-sided 
## 
## Normal residuals      Random mean  Random variance 
##     1.484526e+00     9.682568e-04     5.835753e-05

The Observation value is the observed statistic in the subset. For example, this can be the displacement difference or the probability of overlap in landmark displacement between the two tested specimens in the tested subset (here the front of the skull). This value is then compared to 100 replicated random values that are equivalent of measuring this statistic over any randomly picked landmark displacements. The p value will indicate the type I error when rejecting the null hypothesis (that is: the landmarks in the subset are not different than the other).

To facilitate interpreting these results, it is possible to plot them:

## Graphical parameters
par(mfrow = c(3,2))

## Plotting one result
plot(sp_mean_test_area, main = "Species mean displacement difference")
## Adding the p value
legend("topright", legend = paste("p =", round(sp_mean_test_area$pvalue, 3)), bty = "n")

## Plotting the other results
plot(sp_mean_test_bc, main = "Species mean Bhattacharyya Coefficient")
legend("topright", legend = paste("p =", round(sp_mean_test_bc$pvalue, 3)), bty = "n")

plot(all_sp_test_area, main = "Overall range displacement difference")
legend("topright", legend = paste("p =", round(all_sp_test_area$pvalue, 3)), bty = "n")
plot(all_sp_test_bc, main = "Overall range Bhattacharyya Coefficient")
legend("topright", legend = paste("p =", round(all_sp_test_bc$pvalue, 3)), bty = "n")

plot(ci95_sp_test_area, main = "Overall range displacement difference")
legend("topright", legend = paste("p =", round(ci95_sp_test_area$pvalue, 3)), bty = "n")
plot(ci95_sp_test_bc, main = "Overall range Bhattacharyya Coefficient")
legend("topright", legend = paste("p =", round(ci95_sp_test_bc$pvalue, 3)), bty = "n")

References

Bhattacharyya, A. 1943. “On a Measure of Divergence Between Two Statistical Populations Defined by Their Probability Distributions.” Bulletin of the Calcutta Mathematical Society 35: 99–109.

Silverman, B.W. 1986. Density Estimation for Statistics and Data Analysis. Chapman & Hall/Crc Monographs on Statistics & Applied Probability. Taylor & Francis. http://books.google.ie/books?id=e-xsrjsL7WkC.

Thioulouse, Jean, Daniel Chessel, Sylvain Dole, Jean-Michel Olivier, and others. 1997. “ADE-4: A Multivariate Analysis and Graphical Display Software.” Statistics and Computing 7 (1). Springer: 75–83.

Testing the variation in landmark position

Thomas Guillerme

2019-03-12

What does it do?