5.1 Introduction

In the previous chapters we have primarily considered the geometric invariance of objects under translation, rotation and scale. In some applications the joint study of size and shape is appropriate, and we use the term size-and-shape in this case. Another common term used for size-and-shape is form, and this is particularly common in biological applications.

In molecule comparisons (see Sections 1.4.4 and 1.4.9) the relative scales of the molecules are known, and so we do not need scale invariance. Hence, molecules are usually compared geometrically in size-and-shape space, rather than shape space.

Another example where size-and-shape is appropriate is the analysis of the microfossil data of Section 1.4.16. It is of interest to examine whether the size is related to shape. The form of the microfossils is represented in the size-and-shape space, which is a product of the positive real line (for size) and shape space.

The definition of size-and-shape of an object was given in Definition 3.9. When carrying out size-and-shape analysis all objects must be recorded on the same scale, which we denote as being commensurate in scale. In many applications we do have scale information so a choice needs to be made as to whether to work in size-and-shape space, or to work in shape space and consider the size variable separately.

5.2 Root mean square deviation measures

One way of measuring size-and-shape differences between sets of points is to consider the root mean square deviation (RMSD) between configurations after Procrustes matching. Let X^o₁ and X^o₂ be two configurations and write X₁ = CX^o₁, X₂ = CX^o₂ for the centred configurations. If we match X₁ to X₂ by Procrustes analysis using rotation and translation only, then denote the fitted configuration as X^P₁. In particular, where and X^T₂X₁ = VDU^T with U, V ∈ SO(m) and D diagonal with positive elements except possibly the smallest singular value (see Lemma 4.2).

The RMSD is:

(5.1)

The RMSD has been used in structural bioinformatics to measure size-and-shape difference (Levitt 1976, p. 77), and it is related to the Riemannian size-and-shape distance d_S that we define below in Equation (5.5), where

Also,

(5.2)

for the two configurations (m-vectors) (x_i) and (y_i)  i = 1, …, k in m dimensions, where is the minimizing rotation and is the minimizing translation from Procrustes analysis.

There are in fact two RMSD measures used in the early bioinformatics literature for comparing size-and-shape (Mardia 2010). As well as the coordinate-based RMSD another measure, denoted by RMSD_D, is obtained using the distances between pairs of atoms δ⁽¹⁾_ij = ||x_i − x_j||,  δ⁽²⁾_ij = ||y_i − y_j||, and the measure is given by:

It has been found from empirical studies by Cohen and Sternberg (1980) that

(5.3)

where c ranges from (0, 0.2). A similar relation can be obtained under an independent isotropic Gaussian model for the points (Mardia 2010) where

(5.4)

The result was suggested in Levitt (1976, p. 77) and it is similar to Equation (5.3). Note that RMSD is now used exclusively in bioinformatics rather than RMSD_D, partly for its computational efficiency.

5.3 Geometry

Consider two k-point configurations in m dimensions, (k × m) matrices.

Result 5.1 The Riemannian distance in size-and-shape space is given by:

(5.5)

where S₁, S₂ are the centroid sizes of X^o₁, X₂^o, and ρ is the Riemannian shape distance of Equation (4.12).

Proof: We initially remove location by pre-multiplying by the Helmert submatrix to give the Helmertized coordinates, X₁ = HX^o₁,  X₂ = HX^o₂ as in Equation (3.7). We have sizes S₁ = ||X₁|| = ||CX^o₁|| and S₂ = ||X₂|| = ||CX^o₂||. The Riemannian distance between the size-and-shape of the configurations is found by minimizing the Euclidean distance over rotations (Le 1988, 1995; Ziezold 1994), that is

(5.6)

We have already seen in Equation (4.2) and Equation (4.12) that

where ρ(X₁, X₂) is Riemannian distance in shape space, hence the result follows.

The Riemannian distance in size-and-shape space SΣ^k_m was derived by Le (1988, 1995), and it is also known as the Procrustes distance or intrinsic size-and-shape distance. The Riemannian size-and-shape distance can be computed in R using the function ssriemdist. For example, from

ssriemdist(macf.dat[,,1],macf.dat[,,2])
[1] 13.4061
ssriemdist(macf.dat[,,1],macm.dat[,,1])
[1] 18.29003

we see that the first and second Macaque female skulls are closer in size-and-shape distance than the first female and first male skulls.

We can also consider classical MDS of the schizophrenia data of Section 1.4.5 using size-and-shape distance:

data(schizophrenia)
distmat <- matrix( 0, 28, 28)
for (i in 1:28){
for (j in 1:28){
rho <- ssriemdist( schizophrenia$x[,,i],schizophrenia$x[,,j] )
distmat[i,j] <- rho
}
}
par(mfrow=c(1,2))
plot(cmdscale(distmat),pch=as.character(schizophrenia$group),
                                xlab="MDS1",ylab="MDS2")
plot(schizophrenia$group,centroid.size(schizophrenia$x))

From Figure 5.1 we see that there is a clear difference between the two groups, with the control group being larger in centroid size. Note the difference in the multi-dimensional scaling (MDS) plot compared with Figure 4.3 where shape distances were used (with scale invariance).

Instead of the Helmertized coordinates we could work with the centred configurations X₁ = CX^o₁ = H^THX^o₁, X₂ = CX^o₂ = H^THX^o₂. Note the similarity between Equation (5.5) and the law of cosines. In the m = 2 case complex notation can be used, and if Z₁ and Z₂ are the pre-shapes, then we see that ρ can be regarded as the angle between the complex vectors S₁Z₁ and S₂Z₂. Equation (5.6) is the sum of squares from ordinary partial Procrustes matching, as seen later in Section 7.2.3 and Equation (7.8).

The differential geometric properties of the space are as in the following (Le 1988). The size-and-shape space is a cone with a warped-product metric. The vertex of the cone corresponds to complete coincidence of the points and the unit section (corresponding to unit centroid size) is the shape space Σ^k_m. The singularities in the space are the vertex of the cone and all the points of the rays which meet the unit section at a singularity of Σ^k_m. As in the pure shape case the singularity set has Lebesgue measure zero if dealing with continuous distributions, and so can easily be dealt with in practical applications.

On a historical note Herbert Ziezold presented work in 1974 at the European Meeting of Statisticians which also contained metrics for comparing the size-and-shape of point configurations based on Fréchet (1948) means, and the proceedings were published in 1977. The size-and-shape distance of Equation (5.5) was studied by Ziezold (1977), who considered the m-dimensional case and also studied the 2D case in detail.

5.4 Tangent coordinates for size-and-shape space

Let X and p be Helmertized or centred versions of the configurations X^o, p^o. The Procrustes tangent coordinates of X in the tangent space to the size-and-shape space at pole p are given by:

(5.7)

where is the Procrustes rotation of X onto p of Equation (4.6). Note that ||V|| = d_S(X^o, p^o), and these tangent coordinates are also the inverse exponential map coordinates and also the Procrustes residual coordinates. So, for size-and-shape space this is a very natural choice of tangent coordinates.

Note that V satisfies the following conditions:

1. V is centred V^T1 = 0.
2. p^TV is symmetric: p^TV = V^Tp.
3. The eigenvalues of the symmetric matrix p^T((1 − ||V||²)^1/2p + V) are optimally signed and non-degenerate, see Equation (4.4) and Equation (4.7).

The number of linear constraints on V is m + m(m − 1)/2. The Procrustes size-and-shape tangent coordinates are available in the procGPA routine in the R shapes library, by using the option scale=FALSE. For example,

V <- procGPA( digit3.dat , scale=FALSE )$tan

gives the Procrustes size-and-shape tangent coordinates for the digit 3 data of Section 1.4.2 with the pole as the sample Procrustes size-and-shape mean (see Section 6.7).

5.5 Geodesics

Minimal geodesics in size-and-shape space have a simple form. In particular, consider the minimal geodesic which passes through the size-and-shapes of Helmertized/centred configurations X and p corresponding to original configurations X^o and p^o. This geodesic passes through p in the direction of W = V/||V|| and is given by:

where are size-and-shape tangent coordinates of X with pole p given by Equation (5.7). The length of the minimal geodesic between X and p is ||V|| = d_S(X^o, p^o). Geodesics can be useful, although in practice more complicated models are often required to described changes on manifolds.

5.6 Size-and-shape coordinates

In order to represent size-and-shape we can use any choice of shape coordinates together with any size variable, for example centroid size or its logarithm. A common approach in morphometrics is to use shape tangent coordinates jointly with log centroid size (Mitteroecker et al., 2004).

5.6.1 Bookstein-type coordinates for size-and-shape analysis

Following Mardia (2009) we indicate how to construct Bookstein-type coordinates for size-and-shape (form) analysis which have been used in registering backbones in structural bioinformatics (e.g. see Killian et al. 2007). Let X(k × 3) be the configuration matrix with rows x_i. First note that there are six unknowns in the rotation matrix A (3 Eulerian angles) and a translation vector, but x₁ and x₂ are not enough to determine the form coordinates, since the distance between x₁ and x₂ is fixed, so we have to use x₃ for the one remaining constraint. We first use the point x₁ as the origin so that

Then use y₂ (i.e. x₂) to fix the co-latitude and longitude of the points, that is let (θ₂, ϕ₂) be the polar coordinates of y₂ (z-axis is the north pole), then

where

Now rotate the new coordinates around the coordinates u₃ (or x₃), that is if

then the Bookstein type coordinates for the form in 3D (angle ϕ₃ with new x-axis) are:

where

In bioinformatics, these are termed bond--angle--torsion (BAT) coordinates (e.g. see Killian et al. 2007) and can be given a tree structure.

For 2D, the first two new coordinates are (0, 0) and (d, 0) using x₁ and x₂. That is

where now S(ϕ) is a 2 × 2 rotation matrix

such that

5.6.2 Goodall–Mardia QR size-and-shape coordinates

Let us consider now general m ≥ 2 dimensions for k point configurations where k > m. Given a k × m matrix of landmark coordinates the location information can be removed by pre-multiplying by a suitable matrix, say X_H = HX, where H is the Helmert submatrix, defined in Equation (2.10). A QR decomposition of X_H is given by:

where the (k − 1) × m lower triangular matrix T contains size-and-shape coordinates. Note that T is invariant under the original location and rotation of the configuration. The matrix T has zero entries above the leading diagonal and therefore a total of (k − 1)m − m(m − 1)/2 size-and-shape coordinates.

To obtain the shape coordinates we divide by the centroid size to give W = T/||T|| and so there are q = (k − 1)m − m(m − 1)/2 − 1 shape coordinates. Note that ||T|| = ||X_H|| is the centroid size. This coordinate system was developed by Goodall and Mardia (1992, 1993).

Slight variants of this coordinate system include using a different translation matrix from H [such as B from Equation (3.14)] or dividing by a different scale such as the (1, 1)th element of T (a baseline size). In particular, if k = 3 and m = 2 and B is used to remove location, then the QR decomposition leads to the size-and-shape variables given by:

where A is given by Equation (2.8). Removing size by dividing by t₁₁ we have:

where u^B₃ = t₂₁/t₁₁ and v^B₃ = t₂₂/t₁₁ are Bookstein coordinates of Equation (2.5). Hence, for m = 2 (and m = 3) the QR decomposition is closely related to Bookstein coordinates. As is the case for Bookstein coordinates, one should be careful in interpreting correlations between QR shape variables.

For the general number of m ≥ 2 dimensions the Goodall–Mardia polar shape coordinates are obtained from the QR decomposition as in the following. If X_H = HX = TΓ is the QR decomposition and

then write:

(5.8)

and let ϕ be the generalized spherical polar coordinates of v. The size-and-shape of X is therefore represented by {||T||, ρ, ϕ} and the shape of X is represented by {ρ, ϕ}.

Goodall--Mardia polar coordinates were used by Ball et al. (2008) to construct Ornstein--Uhlenbeck processes for modelling cell shape motion, where a drift was specified to a specific reference shape.

Note that this spherical representation is different from the coordinates for Kendall’s shape sphere for k = 3, from Equation (2.14). Also, ρ is the Riemannian distance from the shape of X to the matrix μ which contains zeros apart from the (1, 1)th entry, which is 1.

5.7 Allometry

Definition 5.1 Allometry involves the study of the relationships between shape and size, and in particular the manner in which shape depends on size.

Traditional methods in allometry involve the fitting of linear or non-linear regression equations between size and/or shape measures, see, for example, Sprent (1972). The notion of allometry was introduced by Huxley (1924, 1932).

We can also investigate allometry using our geometrical framework, for example using linear regression of shape coordinates on size. A shape coordinate is chosen as a response variable and a choice of size variable is used as the explanatory variable, and we could use standard normal-based inference. Simple exploratory plots of the data are useful and the variables might need to be transformed to model a linear relationship.

Example 5.1 Consider the microfossil data of Section 1.4.16. We have actually added 0.5 to the variable U^B to tie in with Bookstein’s original analysis of the data with the baseline sent to (0, 0) and (1, 0), rather than and (Bookstein 1986). There are n = 21 triangles in the dataset. Scatter plots are show in Figure 5.2 and Figure 5.3. We regress and V = V^B on √S and then on log S. There is a slight suggestion of non-linearity in the residual plot from V on √S and perhaps slightly less of a parabolic structure from V regressed on log S, see Figure 5.4. We fit the regression lines

and the fitted values (with standard errors) are , , and . Hence there is strong evidence for allometry between V and log S, as a standard test of H₀: β₂ = 0 would be rejected, with a t statistic of t = 6.6. Examination of a plot of the residuals (Figure 5.4) versus fitted values does not show much structure, although there is possibly a suggestion that one of the observations is an outlier (the observation with size 65 and U = 0.35, V = 0.66).

From Mitteroecker et al. (2004) allometry can also be investigated by appending the logarithm of centroid size to the shape tangent coordinates, and principal components analysis then provides an estimate of the allometric equation from the coefficients of the first principal component. This approach is commonly used in geometric morphometrics. Another example investigating allometry is given in Section 9.5, and further discussion of allometry is given by Klingenberg (1996); Bookstein (2013a); and Cardini et al. (2015) among many others.