Robust Morphometric Analysis based on Landmarks. Applications

TL;DR

This paper enhances Procrustes Analysis for morphometric landmark data by introducing robustness to outliers and providing statistical tests for differences, with practical GIS applications.

Contribution

It introduces a robust version of Procrustes Analysis and derives a new tail probability approximation for the Procrustes Statistic under near-normal conditions.

Findings

Robust classification of individuals using the modified Procrustes method.

New statistical test for significant differences between individuals.

Applications demonstrated with GIS tools.

Abstract

Procrustes Analysis is a Morphometric method based on Configurations of Landmarks that estimates the superimposition parameters by least-squares; for this reason, the procedure is very sensitive to outliers. In the first part of the paper we robustify this technique to classify individuals from a descriptive point of view. In the literature there are also classical results, based on the normality of the observations, to test whether there are significant differences between individuals. In the second part of the paper we determine a Von Mises plus Saddlepoint approximation for the tail probability of the Procrustes Statistic when the observations come from a model close to the normal. We conclude the paper with some applications using the Geographical Information System QGIS.

Tables

Table 1. Table 1: Exact and approximate p-values with g = 3 𝑔 3 g=3

$t$	“exact”	approximate
6	$0^{\prime }149$	$0^{\prime }148$
8	$0^{\prime }077$	$0^{\prime }076$
10	$0^{\prime }042$	$0^{\prime }042$
12	$0^{\prime }024$	$0^{\prime }025$
14	$0^{\prime }016$	$0^{\prime }016$
16	$0^{\prime }011$	$0^{\prime }011$
18	$0^{\prime }007$	$0^{\prime }008$

Table 2. Table 2: Exact and approximate p-values with g = 5 𝑔 5 g=5

$t$	“exact”	VOM appr.	VOM+SAD appr.
9	$0^{\prime }1125$	$0^{\prime }1129$	$0^{\prime }1136$
11	$0^{\prime }0538$	$0^{\prime }0539$	$0^{\prime }0545$
13	$0^{\prime }0251$	$0^{\prime }0249$	$0^{\prime }0253$
15	$0^{\prime }0114$	$0^{\prime }0112$	$0^{\prime }0115$
17	$0^{\prime }0050$	$0^{\prime }0049$	$0^{\prime }0051$
19	$0^{\prime }0022$	$0^{\prime }0022$	$0^{\prime }0023$

Equations

$$M=\left(\begin{array}[]{ccc}c_{11}&\cdots&c_{1k}\\ \cdots&\cdots&\cdots\\ c_{p1}&\cdots&c_{pk}\\ \end{array}\right)$$

$$\min d_F(M_1,M_2)=\min ||M_2-\alpha M_1\text{boldmath}\Gamma -{\text{boldmath}1}_p{\beta }^{\prime }||=$$

$$=\sqrt{trace[(M_2-\alpha M_1\text{boldmath}\Gamma -{\text{boldmath}1}_p{\beta }^{\prime }{)}^{\prime }(M_2-\alpha M_1\text{boldmath}\Gamma -{\text{boldmath}1}_p{\beta }^{\prime })]}$$

$$M_c=\left(M_{c1},...,M_{ck}\right)=\left(\frac{1}{p}\sum _{i=1}^p{c}_{i1},...,\frac{1}{p}\sum _{i=1}^p{c}_{ik}\right)$$

$$CS=\sqrt{\sum _{i=1}^p{d}_E^2(l_i,M_c)}=\sqrt{\sum _{i=1}^p\sum _{j=1}^k(c_{ij}-M_{cj}{)}^2}=\sqrt{\sum _{j=1}^{k}p\cdot Var(c_{.j})}$$

$$M_{cs}=\frac{M}{CS}$$

$$\underset{\text{boldmath}\Gamma }{\min }||M_2-M_1\text{boldmath}\Gamma ||=\underset{\text{boldmath}\Gamma }{\min }\sqrt{trace((M_2-M_1\text{boldmath}\Gamma {)}^{\prime }(M_2-M_1\text{boldmath}\Gamma ))}$$

$$d_F(M_1,M_2)=\sqrt{trace{\left(X_2-\beta X_1\text{boldmath}\Gamma \right)}^{\prime }\left(X_2-\beta X_1\text{boldmath}\Gamma \right)}=$$

$$=\sqrt{1-{\left(\sum _{j=1}^k{\lambda }_j\right)}^2}$$

$$X_2^{\prime }{X}_1=UD{V}^{\prime }$$

$$\text{boldmath}\Gamma =V{U}^{\prime }$$

$$X^{\ast }=X(I_{kp}-x_m^{\prime }{x}_m)$$

$$M_c=\left(M_{c1},...,M_{ck}\right)=\left(\frac{1}{p}\sum _{i=1}^p{c}_{i1},...,\frac{1}{p}\sum _{i=1}^p{c}_{ik}\right)$$

$$CS=\sqrt{\sum _{i=1}^p{d}_E^2(l_i,M_c)}=\sqrt{\sum _{i=1}^p\sum _{j=1}^k(c_{ij}-M_{cj}{)}^2}=\sqrt{\sum _{j=1}^{k}p\cdot Var(c_{.j})}$$

$$MAD=1^{\prime }4826M_e\left\{\left|X_i-M_e(X_i)\right|\right\}$$

$$MS=\sum _{j=1}^{k}MAD(c_{.j})$$

$$M_{rs}=\frac{M}{MS}$$

$$||X_2-X_1|{|}^2=trace\left[(X_2-X_1{)}^{\prime }(X_2-X_1)\right].$$

$$\underset{\text{boldmath}\Gamma }{\min }||X_2-X_1\text{boldmath}\Gamma |{|}^2=\underset{\text{boldmath}\Gamma }{\min }trace\left[(X_2-X_1\text{boldmath}\Gamma {)}^{\prime }(X_2-X_1\text{boldmath}\Gamma )\right].$$

$$G(X_1,X_2)=\underset{\text{boldmath}\Gamma }{\min }||X_2-X_1\text{boldmath}\Gamma |{|}^2.$$

$$X_2=X_1+\eta \text{boldmath}e$$

$$G(X_1,X_2)\approx {\eta }^2{\chi }_g^2$$

$$\frac{X_2-X_1}{\eta }⇝(1-ϵ)N(0,1)+ϵN(0,\nu )$$

$$P_{F}\{G_{s}(X_{1},X_{2})>t\}\simeq g\int_{-\infty}^{\infty}P\{\chi_{g-1}^{2}>t-(\mbox{$\frac{x-\mu}{\sigma}$})^{2}\}\,dF(x)-(g-1)P\{\chi_{g}^{2}>t\}.$$

$$p_g^F=P_F\{G_s(X_1,X_2)>t\}\simeq p_g^{\Phi }+\int \text{mbox}TAIF(x;t;{\chi }_g^2,{\Phi }_{\mu ,\sigma })dF(x)$$

$$\text{mbox}TAIF(x;t;{\chi }_g^2,{\Phi }_{\mu ,\sigma })$$

$$(1-ϵ)N(0,1)+ϵN(0,\nu )$$

$$p_g^F\simeq (1-gϵ)P\{{\chi }_g^2>t\}+gϵ{\int }_{-\infty }^{\infty }P\{{\chi }_{g-1}^2>t-x^2\}d{\Phi }_{0,\nu }(x).$$

$$P_F\left\{G_s(X_1,X_2)>t\right\}\simeq P\left\{{\chi }_g^2>t\right\}-B+B{\int }_{-\infty }^{\infty }\frac{\sqrt{g}}{\sqrt{t}}{e}^{\frac{(t-g)(x-\mu {)}^2}{2t{\sigma }^2}}dF(x)$$

$$P_{\Phi }\left\{\frac{G_s(X_1,X_2)}{g}>t\right\}$$

Full text

Keywords: Robustness; Morphometrics; Von Mises expansion; Saddlepoint approximations; Geographical Information System QGIS

1 Introduction

This paper is about a robust classification problem of $n$ individuals based on their shapes, i.e., using their geometric information. The usual (classical or robust) methods based on a Multivariate Analysis can not extract all the geometric information from the individuals. For this reason, in recent years morphometrics methods based on Configurations of landmarks have been developed. A landmark is a peculiar point whose position is common in all the individuals to classify. For instance, when we classify skulls, the landmarks could be the center of the supraorbital arch, the chin, etc.; or, if we classify projectile points found in an archaeological site, the landmarks could be the ends of the points.

In all the cases, the mathematical (geometric) information that we obtain from individuals is the $k$ coordinates of their $p$ landmarks, $l_{i}=(c_{i1},...,c_{ik})$, $\;i=1,...,p$.

The matrix of landmarks coordinates is called a Configuration. For each individual with $p$ landmarks of dimension $k$ (2 or 3) we shall have a collection of landmark coordinates expressed in $p\times k$ matrix as

[TABLE]

2 Classical Morphometric Analysis from a Descriptive Point of View

As we have mentioned before, we shall use the shape of the individuals in their classification. Shape is a property of an object invariant under scaling, rotation and translation; otherwise, for instance, an object and itself with double size could be classified into two different groups.

There are many morphometric methods; see for instance C or DM . In this paper we shall consider Superimposition Methods; namely, Procrustes Analysis, obtaining the Procrustes coordinates with it, adapting the Configurations to a common (local) reference system and matching them at the common center. For these reasons, a Local Coordinate Reference System is needed and a Geographical Information System will be very useful.

A common graphical representation of a Configuration is a scatter plot of its landmarks coordinates. Joining them with segments we obtaining a polygon as, for instance, in Fig. 1, where the landmarks coordinates are the vertices of the polygon.

As we have said above, to classify individuals we have first to remove the effect of Size (scale), Location (translation) and Orientation (rotation) to standardize them and match them in a common center (the centroid of the polygon) in order to make them comparable.

To apply the Procrustes superimposition method we have to estimate by least-squares the superimposition parameters $\alpha$, $\beta$ and ${\boldmath\Gamma}$ (scale, translation and rotation) in order to minimize the full Procrustes distance $d_{F}$ between Configurations $M_{1}$ and $M_{2}$, i.e.,

[TABLE]

[TABLE]

where $\alpha$ is a scalar representing the Size, $\beta$ is a vector of $k$ values corresponding to a Location parameter formed by the centroid coordinates, ${\boldmath 1}_{p}$ is a column vector of dimension $p\times 1$ and ${\boldmath\Gamma}$ a $k\times k$ square rotation matrix.

The idea that we pursue with this transformation is to match both Configurations, i.e., a superimposition of $M_{1}$ onto $M_{2}$.

2.1 Removing the Size Effect

The first step we must take in Procrustes Analysis to standardize Configurations is to remove the Size effect. If, as usual, we consider as center the centroid-mean of dimension $k$ (sample mean by columns) defined by

[TABLE]

and easily computed with R as (RCT )

apply(M,2,mean)

the Centroid-mean Size is defined as

[TABLE]

being $d_{E}^{2}(l_{i},M_{c})$ the square of the Euclidean distance between the $i$th landmark $l_{i}$ and the centroid-mean $M_{c}$. Hence, the Centroid-mean Size depends on the sample variance and so, it will be very sensitive to outliers. This size can be computed as

sqrt(sum(apply(M,2,var)*(p-1)))

The coordinates of a scaled Configuration are now calculated dividing the original coordinates by $CS$

[TABLE]

In Fig. 2 we see the previous Configuration (with red centroid) and the scaled to Centroid-mean Size equal to 1 (the Configuration with green centroid).

2.2 Removing Location by Translation

We remove the Location effect translating the Configuration matrix so that its centroid is the new origin. We do this with the R sentence

scale(M,scale=F)

In Fig. 3 we have the previous Configurations and the centered one (with blue centroid).

2.3 Removing Orientation by Rotation

After the effect of Size and Location have been removed, we estimate (by least-squares) the rotation matrix ${\boldmath\Gamma}$ minimizing the distance between Configurations $M_{1}$ and $M_{2}$, i.e.,

[TABLE]

where ${\boldmath\Gamma}$ is a $k\times k$ square rotation matrix, a matrix that must be determined in order to maximize the correlation between the two sets of landmarks, i.e., to minimize the distance between landmarks. More precisely:

If $M_{1}$ and $M_{2}$ are two Configurations and $X_{1}$ and $X_{2}$ the corresponding centered Configurations scaled to unit Centroid-mean Size, the (full) Procrustes distance is defined as

[TABLE]

[TABLE]

where $\lambda_{j}$ are the diagonal elements of matrix $D$ in the singular-value decomposition of $X_{2}^{\prime}X_{1}$:

[TABLE]

But in fact, the previous problem is a known mathematical issue: If we have the previous singular-value decomposition of $X_{2}^{\prime}X_{1}$, the rotation matrix we are looking for is

[TABLE]

In Fig. 4 we have two Configurations, before and after rotated one of them, according to the previous method.

2.4 More than two Configurations (Generalized Procrustes Analysis)

In the previous sections we have done, in three steps, what is called a classical Partial Procrustes Analysis because we have compared, from a descriptive point of view, two Configurations.

If we have more than two Configurations we have to do what is called a Generalized Procrustes Analysis in three steps:

All the Configurations must be standardized, i.e., their centroids matched at a common origin and scaled to unit size. 2. 2.

We have to define a consensus or average Configuration as a reference, called Mean Shape because, in fact, it is the sample mean of all the Configurations, by vertices, i.e., the mean of the vertices $(i,j)$ (homologous coordinates) for the $n$ Configurations. The Mean Shape Configuration can be computed as

apply(M,c(1,2),mean)

Finally, we perform a (partial) Procrustes superimposition between the Mean Shape and each Configuration.

2.5 Configuration Projection onto the Tangent Space

The shape space defined by the previous Procrustes superimposition method is non-Euclidean and corresponds to a curved surface. This means that the distance between two landmarks is not the length of the segment joining them and hence, we cannot apply traditional statistics to the Procrustes coordinates of the $n$ individuals.

From a convenient point of view, the $n$ individuals are aligned in the $n\times kp$ matrix $X$ and then projected onto the (Kendall) tangent space, where the vectorized Mean Shape $x_{m}$ (i.e., a vector of dimension $kp\times 1$) is the contact point between spaces. The projected (tangent) coordinates are obtained in a matrix $X^{*}$ as

[TABLE]

where $I_{kp}$ is the $\;\;kp\times kp\;\;$ identity matrix.

Then, we can apply the usual statistical techniques to these projected coordinates, for example, classifying the resulting observations with the scores of their Principal Components.

Example 1

In paper HD , 59 gorilla skulls were considered. We know, in the example but not in a real case, that 30 of them are female and 29 male. In their paper, 8 landmarks were considered. If we represent these data in a scatter plot we obtain Fig. 5 where no apparent classification between males and females is observed.

If we make the four previous steps of the Generalized Procrustes analysis and conclude with a Principal Component analysis of the scores, we obtain Fig. 6 where we cannot appreciate the two groups very clearly although the vertical bar at PC1=0 is the usual classification rule taken for this example.

3 Robust Morphometric Analysis from a Descriptive Point of View

Let us consider again the data of Example 1 plus a Configuration, that we call a, for which we replace the coordinates of the third landmark as in the following diagram:

[,1] [,2] [,1] [,2] [1,] 36 187 [1,] 36 187 [2,] 59 -31 [2,] 59 -31 [3,] 0 0 ----> [3,] 30 30 [4,] 0 36 a = [4,] 0 36 [5,] 12 102 [5,] 12 102 [6,] 38 171 [6,] 38 171 [7,] 91 103 [7,] 91 103 [8,] 100 19 [8,] 100 19

If we give the four previous steps to perform a Generalized Procrustes analysis, we obtain the classification given in Fig. 7 where all the individuals are in one group except outlier a.

But in Fig. 8 we see that a is in the bulk of the data and also the mean shape in Fig. 9. Hence, no apparent solution is clear.

3.1 Removing the Size Effect in a Robust Way

We propose, instead of using the centroid-mean

[TABLE]

as before, that essentially is a sample mean computed with

apply(M,2,mean)

to use the median (or the trimmed-mean) by columns with the following two R sentences,

apply(M,2,median)

apply(M,2,mean,trim = .2)

obtaining in this way a more robust centroid. Now, instead of considering the Centroid-mean Size $CS$ that, as we saw before, is essentially the variance

[TABLE]

computed with

sqrt(sum(apply(M,2,var)*(p-1)))

we propose to use the Median Absolute Deviation $MAD$, defined as

[TABLE]

obtaining with it what we call the Centroid-median Size

[TABLE]

computed as

sum(apply(M,2,mad))

and that satisfies the Size invariance property $\,MS(aM)=aMS(M)\,$ for any positive scalar $a$.

In this way we obtain a robust size measure. For instance, considering the following two configurations A and B

A > B [,1] [,2] [,1] [,2] [1,] 2 0 [1,] 20 0 [2,] 1 1 [2,] 1 1 [3,] 0 0 [3,] 0 0 [4,] 0 -1 [4,] 0 -1 [5,] 2 -2 [5,] 2 -2

that differ in just a wrong digit in the first landmark of Configuration B, the classical Centroid-mean Size is very sensitive:

sqrt(sum(apply(A,2,var)*4)) [1] 3.03315 sqrt(sum(apply(B,2,var)*4)) [1] 17.44706

but not the the Centroid-median Size:

sum(apply(A,2,mad)) [1] 2.9652 sum(apply(B,2,mad)) [1] 2.9652

And what is more importante, this new size measure keeps the relative size of the Configurations avoiding a possible masking effect. For instance, in Fig. 10, if we divide both Configurations by the classical Centroid-mean Size we obtain Fig. 11 and they would probably be classified in the same final group. Nevertheless, standardizing them with the new robust Centroid-median Size, we see in Fig. 12 that the differences in size between them, remain.

Hence, instead of dividing the coordinates of the Configuration by the classical Centroid-mean Size $CS$, we propose to divide the configuration $M$ (the coordinates) by the robust Centroid-median Size to distinguish between individuals in a better way, avoiding a possible masking effect,

[TABLE]

3.2 Removing Location in a Robust Way

In the same way as we have removed the location effect in a classical way, translating the Configuration matrix so that its centroid-mean was the new origin, with the sentence

scale(M,scale=F) = scale(M,scale=F,center=apply(M,2,mean))

subtracting the mean of each column to the whole column, in the robust version we subtract the median with the sentence

scale(M,scale=F,center=apply(M,2,median))

being the centroid-median the new origin.

After robustifying with respect scale and location we keep the classical rotation matrix for the robust coordinates. These three steps are in our new R function rpgpa1 that can be obtained on request from the authors. We conclude the process with the same coordinates projection formula than before.

3.3 More than two Configurations

If there are more than two Configurations, the key point in the robustification process is the selection of a robust mean shape, that in the classical Morphometric analysis is the sample mean of the Configuration coordinates. In our robust version we propose to choose as mean shape the median of the Configuration coordinates, obtaining Fig. 13 for the gorillas example (after doing a classical Principal Component analysis of the scores).

Considering the $0^{\prime}2$-trimmed mean as mean shape we obtain Fig. 14. Finally, considering the $0^{\prime}1$-trimmed mean as mean shape we obtain Fig. 15.

These three options are in new R function rpgpa2 that can be composed with rpgpa1.

4 Classical Morphometric Analysis from an Inferential Point of View

Instead of considering a descriptive morphometric analysis it is more interesting to test if there are significant differences between two Configurations. From a classical point of view, we have the following result in LC and S : If $X_{1}$ and $X_{2}$ are two scaled and centered Configurations with $p\times k$ landmarks, the Residual Distance between Configurations $X_{1}$ and $X_{2}$ is defined as

[TABLE]

As saw in the previous sections, the $k\times k$ square rotation matrix ${\boldmath\Gamma}$ is determined such that the Procrustes distance between these two Configurations $X_{1}$ and $X_{2}$ (i.e., between landmarks) is minimal

[TABLE]

This minimum (i.e., after matching, i.e., after translation, rotation and scaling) that we obtain is called the Procrustes statistic:

[TABLE]

Under the null hypothesis $H_{0}$ that there is no systematic difference between Configurations $X_{1}$ and $X_{2}$, i.e., they belong to the same group, or more precisely, that they are of the form

[TABLE]

where all the elements of Configuration $\;{\boldmath e}\;$ are i.i.d. $N(0,1)$, then

[TABLE]

i.e., $G_{s}(X_{1},X_{2})=G(X_{1},X_{2})/\eta^{2}\approx\chi^{2}_{g}\,$, where $\,g=kp-k(k+1)/2-1\,$ obtaining so, a way to compute the tail probabilities (p-values) for testing $H_{0}$. It must be $p>(k+1)/2+1/k$ and obviously an integer.

5 Robust Morphometric Analysis from an Inferential Point of View

The standard normality of the landmarks is a very difficult assumption to assume and check. For this reason we shall use robust methods for testing $H_{0}$ assuming that the $p\times k$ landmarks of $e$ follow, not a standard normal distribution but a contaminated normal model:

[TABLE]

In this section we are going to compute the tail probabilities (p-values), assuming this contaminated model, using a VOM+SAD approximation.

We use this scale contaminated normal mixture model because the Configurations are matched at the common centroid that is the new origin and equal to 0, being the contamination in the scale the source of contamination in the observations.

5.1 Von Mises Approximations for the p-value of the Procrustes Statistic

In order to test the null hypothesis $H_{0}$ that there is no systematic difference between the standardized Configurations $X_{1}$ and $X_{2}$, using the Procrustes statistic $G_{s}(X_{1},X_{2})$ that follow a $\chi^{2}_{g}$ distribution under a normal model, we have the following result.

Proposition 1

Let $G_{s}(X_{1},X_{2})$ be the Procrustes statistic, that follows a $\chi_{g}^{2}$ distribution when the underlying model is a normal distribution, $\Phi_{\mu,\sigma}$. If the previous null hypothesis $H_{0}$ holds, the von Mises (VOM) approximation for the functional tail probability (if $F$ is close to the normal $\Phi_{\mu,\sigma}$) is

[TABLE]

Proof

The von Mises (VOM) approximation for the functional tail probability is (if $F$ is close to the normal $\Phi_{\mu,\sigma}$)

[TABLE]

where TAIF is the Tail Area Influence Function defined in FR .

Replacing the normal model by the contaminated normal model $\Phi^{\epsilon}=(1-\epsilon)\,\Phi_{\mu,\sigma}+\epsilon\,\delta_{x}\;$ and computing the derivative at $\epsilon=0$ we obtain that

[TABLE]

integrating now, we obtain the result. ∎

Considering a scale contaminated normal (SCN) model

[TABLE]

the VOM approximation is

[TABLE]

In Table 1 appear the Exact values (obtained through a simulation of 100.000 samples) and the VOM approximations when $\epsilon=0^{\prime}05$, $\nu=2$ and $g=3$.

To obtain the previous numerical results we had to deal with numerical integration. Sometimes, we would like to have analytic expressions of $p_{g}^{F}$ to value the effect of contamination $\epsilon$, etc. For this reason and for controlling the relative error of the approximation, in the next section we shall compute the Saddlepoint approximation for the p-value of the Procrustes Statistic.

5.2 Saddlepoint Approximations for the p-value of the Procrustes Statistic

Using Lugannani and Rice formula, LR , for the sample mean of $g$ independent square normal variables, we obtain the VOM+SAD approximation given in the next result.

Proposition 2

Let $G_{s}(X_{1},X_{2})$ be the Procrustes statistic, that follows a $\chi_{g}^{2}$ distribution when the underlying model is a normal distribution, $\Phi_{\mu,\sigma}$. If the null hypothesis $H_{0}$ holds, the saddlepoint approximation of the von Mises expansion, VOM+SAD approximation, for the functional tail probability (if $F$ is close to the normal $\Phi_{\mu,\sigma}$) is

[TABLE]

where $B=\frac{g\,\sqrt{g}}{\sqrt{\pi}\,(t-g)}\,e^{-(t-g-g\cdot\log(t/g))/2}$.

Proof

If $G_{s}(X_{1},X_{2})$ follows a $\chi_{g}^{2}$ distribution, and $Y_{1},...,Y_{g}$ are $g$ independent gamma distributions $\gamma(1/2,1/2)$ with moment generating function $M$ and cumulant generating function $K=\log M$, it is, following LR , D or J ,

[TABLE]

where $\Phi_{s}$ and $\phi_{s}$ are the cumulative distribution and density functions of the standard normal distribution.

If $K$ is the cumulant generating function, that is the functional of $\Phi_{\mu,\sigma}$,

[TABLE]

and $z_{0}$ is the (functional) saddlepoint, i.e., it is the solution of the equation $K^{\prime}(z_{0})=t$, the functionals that appear in (3) are

[TABLE]

[TABLE]

As we saw before, the VOM approximation for the tail probability depends on the TAIF. To obtain the TAIF of $\;G_{s}(X_{1},X_{2})/g\;$ at $\;\Phi_{\mu,\sigma}\;$ we have to replace the model $\Phi_{\mu,\sigma}$ by the contaminated model $\Phi^{\epsilon}=(1-\epsilon)\Phi_{\mu,\sigma}+\epsilon\,\delta_{x}$ in all the functionals in the right side of (3) that depend on $\Phi_{\mu,\sigma}$ , and then to obtain the derivative at $\epsilon=0$; this process is represented with a dot over the functional. Since $\;\phi^{\prime}_{s}(w)=-\phi_{s}(w)\,w\;\;$ and $\;\phi_{s}(w)\leq 1\;$, we obtain that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

because the functionals $w_{1},\stackrel{{\scriptstyle\bullet}}{{w}}_{1},r_{1}$ and $\stackrel{{\scriptstyle\bullet}}{{r}}_{1}$ do not depend on $g$. Since

[TABLE]

it will be

[TABLE]

Hence, we have to compute the influence functions $\stackrel{{\scriptstyle\bullet}}{{K}}(z_{0})$ and $\stackrel{{\scriptstyle\bullet}}{{z}}_{0}$. To do this, because

[TABLE]

from the saddlepoint equation, $\;K^{\prime}(z_{0})=t\;$, we obtain

[TABLE]

Replacing again the model by the contaminated model $\Phi^{\epsilon}=(1-\epsilon)\,\Phi_{\mu,\sigma}+\epsilon\,\delta_{x}\;$ before obtaining the derivative at $\epsilon=0$, and making the change of variable $(u-\mu)/\sigma=y$, we obtain

[TABLE]

i.e.,

[TABLE]

In a similar way, we obtain that

[TABLE]

Also it is

[TABLE]

Therefore, from (4), it will be

[TABLE]

where

[TABLE]

From (1), we obtain now the VOM+SAD approximation for the p-value of the test statistic $G_{s}(X_{1},X_{2})/g$,

[TABLE]

and from this, we obtain the approximation (2) for the test statistic $G_{s}(X_{1},X_{2})$. ∎

If $F$ is the location contaminated normal mixture (LCN),

[TABLE]

the VOM+SAD approximation is

[TABLE]

In Table 2 appear the Exact values (obtained through simulation of 100.000 samples), the VOM and the VOM+SAD approximations when $\epsilon=0^{\prime}01$, $\theta=1$ and $g=5$.

Corollary 1

To test the null hypothesis $H_{0}$ that there is no systematic difference between the standardized Configurations $X_{1}$ and $X_{2}$ with $p\times k$ landmarks (i.e., they belong to the same group) using the Procrustes statistic $G(X_{1},X_{2})$ and assuming that the error difference between Configurations

[TABLE]

follow a scale contamination normal model $\,(1-\epsilon)N(0,1)+\epsilon N(0,\nu)\,$, the VOM+SAD approximation for the tail probability (p-value) is

[TABLE]

[TABLE]

where $g=kp-k(k+1)/2-1$. It must be $p>(k+1)/2+1/k$ and obviously an integer.

Then, if $k=2$, it is $g=2p-4$ and $p>2$. And if $k=3$, it is $g=3p-7$ and $p\geq 3$.

6 Applications

In this section we are going to consider the following Example in which we make two comparisons using the previous theory.

Example 2

We are going to consider two test to check if there are or not significance differences between two arrows of Notch tips and bay leaves, of Solutrense period, arrows that were found in caves of Asturias (Spain). We shall make this analysis using a photo of the “Museo Arqueológico de Asturias” (Oviedo), included in QGIS as a raster layer, Fig. 16.

In this figure we see large differences among the arrows except in two pairs of arrows: Arrows 1 and 3 and arrows 5 and 6, Fig. 17. Hence, we are going to test the null hypothesis of no significance differences between arrows 1 and 3, and then, with another test, we shall check the null hypothesis of no significance differences between arrows 5 and 6.

To do this, we first create the polygons in QGIS marking the landmarks with the mouse. We consider $p=7$ landmarks. Also, with QGIS we export the coordinates of the landmarks, that are:

Arrow 3

punta3<-matrix(c( 151.77884,-794.21946, 384.34151,-714.48369, 533.84608,-706.17788, 543.81305,-756.01273, 587.00326,-794.21946, 583.68094,-842.39315, 384.34151,-849.03780), ncol=2,byrow = T)

Arrow 1

punta1<-matrix(c( 157.59291,-1934.60710, 444.97392,-1841.58204, 640.99102,-1848.22668, 650.95799,-1891.41689, 734.01609,-1917.99548, 735.67725,-1966.16918, 428.36230,-1977.79731), ncol=2,byrow = T)

Arrow 6

punta6<-matrix(c( 1170.17428,-1072.54821, 1423.18465,-971.34406, 1550.56229,-974.83386, 1557.54188,-1039.39513, 1606.39906,-1074.29311, 1608.14396,-1156.30337, 1410.97036,-1170.26256), ncol=2,byrow = T)

Arrow 5

punta5<-matrix(c( 1119.57220,-1510.51789, 1327.21520,-1383.14025, 1533.11329,-1386.63005, 1529.62350,-1465.15051, 1639.55214,-1564.60976, 1637.80724,-1618.70164, 1349.89889,-1625.68123), ncol=2,byrow = T)

Comparison between Arrows 1 and 3

After removing the effect of Size (scale), Location (translation) and Orientation (rotation) to standardize the individuals, we match them at the common centroid obtaining the polygons of Fig. 18.

The minimum Residual Distance between configurations (arrows), i.e., the value of the Procrustes statistic for testing the null that “No significance differences exist between arrows 1 and 3” is $0.01567681$:

tamapunta1<-sqrt(sum(apply(punta1,2,var)(7-1))) spunta1<-scale(punta1/tamapunta1,scale=F) tamapunta3<-sqrt(sum(apply(punta3,2,var)(7-1))) spunta3<-scale(punta3/tamapunta3,scale=F) library(shapes) (procdist(spunta1,spunta3,type="partial"))^2 [1] 0.01567681

sd(spunta1-spunta3) [1] 0.03502222 (1)

Because of (1), choosing $\eta=0.03502222$, we shall obtain a standard normal distribution for $\,(X_{2}-X_{1})/\eta\,$ and hence

[TABLE]

being $g=2p-4=10$. Then, the p-value of this classical test will be

[TABLE]

[TABLE]

accepting the null hypothesis of no significance differences between Arrows 1 and 3.

Nevertheless, using the Mahalanobis distance we can conclude that the errors do not follow a multivariate normal distribution,

dipuntas2<-mahalanobis(spunta1-spunta3,colMeans(spunta1-spunta3),

• var(spunta1-spunta3))

ks.test(dipuntas2,"pchisq",7)

   One-sample Kolmogorov-Smirnov test

data: dipuntas2 D = 0.76612, p-value = 8.265e-05 alternative hypothesis: two-sided

Hence, to assume a common $\eta$ for all the $c_{ij}$ such that

[TABLE]

is unrealistic. It is better to consider a model

[TABLE]

and to use the VOM+SAD approximation (5), programmed into the R function   apro3$(g,\nu,\epsilon,t)\;$ to compute the p-value. Obtaining from the data $\eta=0.04020902$ and $\nu=0.032261$, we have

[TABLE]

[TABLE]

because

apro3(10,0.032261,0.1,22.95901) [1] 0.006318776

Then, because

[TABLE]

we reject, in a more robust way, the null hypothesis of no significance differences between Arrows 1 and 3.

Comparison between Arrows 5 and 6

After removing the effect of Size (scale), Location (translation) and Orientation (rotation) to standardize the individuals we match them at the common centroid having the polygons of Fig. 19.

The minimum Residual Distance between configurations (arrows), i.e., the value of the Procrustes statistic for testing the null that “No significance differences exist between arrows 5 and 6” is $0.03711933$,

tamapunta5<-sqrt(sum(apply(punta5,2,var)(7-1))) spunta5<-scale(punta5/tamapunta5,scale=F) tamapunta6<-sqrt(sum(apply(punta6,2,var)(7-1))) spunta6<-scale(punta6/tamapunta6,scale=F) library(shapes) (procdist(spunta5,spunta6,type="partial"))^2 [1] 0.03711933

sd(spunta5-spunta6) (2) [1] 0.05343598

Because of (2), if $\eta=0.05343598$ we shall obtain a standard normal distribution for $\,(X_{2}-X_{1})/\eta\,$. Then,

[TABLE]

and hence, the p-value of this classical test will be

[TABLE]

[TABLE]

accepting the null hypothesis of no significance differences between Arrows 5 and 6.

Nevertheless, using the Mahalanobis distance we can conclude that the errors do not follow a multivariate normal distribution,

dipuntas<-mahalanobis(spunta6-spunta5,colMeans(spunta6-spunta5),

• var(spunta6-spunta5))

ks.test(dipuntas,"pchisq",7)

    One-sample Kolmogorov-Smirnov test

data: dipuntas D = 0.76677, p-value = 8.093e-05 alternative hypothesis: two-sided

Then, to assume a common $\eta$ for all the $c_{ij}$ such that

[TABLE]

is unrealistic. It is better to consider a model

[TABLE]

and to use the VOM+SAD approximation (5), programmed as the R function   apro3$(g,\nu,\epsilon,t)\;$ to compute p-values. From the data we obtain $\eta=0.06834322$ and $\nu=0.0389347$ and hence,

[TABLE]

[TABLE]

because

apro3(10,0.0389347,0.1,7.947111) [1] 0.5405565

Then, because

[TABLE]

we finally accept, in a more robust way, the null hypothesis of no significance differences between Arrows 5 and 6.

7 Conclusions

Classical Morphometric Analysis based on Landmarks is reviewed from a descriptive and inferential point of view. Because both are based on sample means and least-squares they are not robust.

We first robustify the descriptive measures proposing robust ones. Then we consider a Contaminated Normal Model distribution instead of a classical Normal one to make robust inferences. Namely, for this mixture model we obtain an von Mises approximation for the p-value of a test for the null hypothesis of no significance differences between two individuals based on their shapes.

We also obtain a very accurate saddlepoint approximation of this von Mises approximation. We conclude the paper with some applications using QGIS as Geographical Information System.

Acknowledgements.

This work is partially supported by Grant HAR2015-68876-P from Ministerio de Economía y Competitividad (Spain).