Burnaby's orthogonal subspace computation

Performs Burnaby's multivariate approach for standardizing morphometric data relative to nuisance factors by computing and projecting data onto an orthogonal subspace.

burnaby(x, vars = NULL, tree = NULL, axmat = NULL)

Arguments

x: A matrix with one or more variables as columns and observations as rows, representing the variables the user desires to be rendered independent (orthogonal) to vars.
vars: A matrix with one or more variables as columns and observations as rows, representing the source of variation the user wants to remove from x.
tree: A "phylo" object containing a phylogenetic tree. Tip labels should match the row names from x and vars.
axmat: An optional matrix of axes coefficients to render x orthogonal to. If provided, is used instead of vars to orthogonalize x.

Value

A "burnaby" or "phy_burnaby" object, containing:

$x: a matrix with the scores of observations in the new ordination axes.
$sdev: the standard deviations of the principal components (i.e., the square roots of the eigenvalues of the covariance matrix).
$rotation: a matrix of eigenvector coefficients.
$center: the mean values of the original variables (or their phylogenetic mean, if tree is provided).
$totvar: the sum of the variances of the original data.

Details

Originally devised by Burnaby (1966) to eliminate allometric variation from morphometric data sets by computing a subspace that is orthogonal (i.e., mathematically independent) to an allometric vector, and then projecting the data onto these new set of axes (However, this approach can be used to standardize morphometric data for any type of arbitrary vector(s) -- commonly one or more nuisance factors or covariates introducing variation the user wishes to exclude from the analysis).

In particular, this function is intended to provide a way for computation of ordination (sub)spaces orthogonal to other specific sources of variation, as well as visualization of shape variation associated to their ordination axes, most commonly as part of the mspace %>% proj_* pipeline. If the user wishes to use the Burnaby approach for standardization of morphometric (shape) data, this is better achieved through detrend_shapes by setting method = "orthogonal".

The subspace produced in this way has one less dimension for each axis or variable used for standardization, relative to the original variables. However, the user should be aware that, when analyzing shape variation, a rather large portion of variation in x that is not orthogonal to vars will end up 'hidden' in the lowest ordination axes (i.e., those effectively lost during standardization procedures that remove non-shape variation, and therefore showing null eigenvalues). Therefore, projection of the original shape variables into the orthogonal subspace should be attempted after removing (at least) these null axes.

References

Burnaby, T. P. (1966) Growth-invariant discriminant functions and generalized distances. Biometrics, 22, 96–110.

Claude, J. (2008). Morphometrics with R. Springer Science & Business Media, 316.

Revell, L. J. (2009). Size-correction and principal components for interspecific comparative studies. Evolution, 63, 3258-3268.

Klingenberg, C. P. (2016). Size, shape, and form: concepts of allometry in geometric morphometrics. Development Genes and Evolution, 226(3), 113-137.

Examples

#load packages
library(geomorph)
library(Morpho)

#extract input data: shapes and sizes from S. vacaensis.
data("shells3D")
index <- shells3D$data$species == "vacaensis"
allshapes <- shells3D$shapes
shapes <- allshapes[,,index]
sizes <- log(shells3D$sizes[index])

#construct template for visualization
template <- tps3d(x = shells3D$mesh_meanspec,
                  refmat = allshapes[,,findMeanSpec(allshapes)],
                  tarmat = expected_shapes(shapes))


#find shape subspace orthogonal to size using the Burnaby approach for
#S. vacaensis
burn <- burnaby(x = geomorph::two.d.array(shapes), vars = sizes)

#inspect results
names(burn) #the contents of the resulting object
#> [1] "x"        "rotation" "center"   "sdev"     "totvar"  
exp_var(burn) #interpretation of eigenvalues is complex! not yet implemented
#>         variance cummulative
#> Axis 1   7.18323     7.18323
#> Axis 2   5.23785    12.42108
#> Axis 3   3.45779    15.87887
#> Axis 4   2.09618    17.97505
#> Axis 5   1.66424    19.63929
#> Axis 6   1.17258    20.81187
#> Axis 7   0.84681    21.65868
#> Axis 8   0.56362    22.22230
#> Axis 9   0.49068    22.71299
#> Axis 10  0.33959    23.05257
plot(burn$x) #ordination



#compare shape variation as summarized by different methods

if (FALSE) {

#build morphospace for S. vacaensis using Burnaby's approach (orthogonal
#to size variation)
msp1 <- mspace(shapes, template = template, cex.ldm = 0, bg.models = "gray",
               adj_frame = c(0.9, 0.9), asp.models = 1.1, alpha.models = 0.7,
               FUN = burnaby, vars = sizes) %>%
  proj_shapes(shapes, pch = 16)

#compare against a morphospace made using regular PCA
mspace(shapes, template = template, cex.ldm = 0, bg.models = "gray",
       adj_frame = c(0.9, 0.9), asp.models = 1.1, alpha.models = 0.7) %>%
  proj_shapes(shapes, pch = 16)


#extract sample and background shape models for subspace for S. vacaensis
grid_shapes <- extract_shapes(msp1)$shapes
na_shapes <- extract_shapes(msp1, scores = msp1$projected$scores[,1:2])$shapes

#perform regression of size on shape for S. vacaensis
reg <- lm(two.d.array(shapes) ~ sizes)

#build morphospace for the entire dataset using regular PCA, project both the
#empirical shapes, orthogonal subspace range, and size-standardized shapes
mspace(allshapes, template = template, cex.ldm = 0,
       bg.models = "gray", adj_frame = c(0.9, 0.9), asp.models = 1.1,
       alpha.models = 0.7) %>%
  proj_shapes(shapes, pch = 1) %>%
  proj_shapes(na_shapes, pch = 21, bg = "red") %>%
  proj_axis(reg, lwd = 2) %>%
  proj_groups(grid_shapes, alpha = 0.2, col = "red")
}