Multiple scaled symmetric distributions in allometric studies

1 Introduction

Allometry can be roughly devised as a tool to study the joint relation between parts (morphometric variables) in various organisms [1], [2], with the measurements being often log-transformed (for the practical and theoretical reasons why it is often useful to transform data to logarithms, see [3], [4], [5], [6]). So, the statistical models that are commonly considered are inherently multivariate (almost always bivariate), and mainly aim to:

1. reproduce the joint distribution of the log-transformed morphometric variables;

2. summarize the morphometric variation.

Concerning the first aim, the multivariate normal (MN) distribution with mean μ and covariance matrix Σ played, and still plays, a special rule. However, many empirical studies show that the MN is often not appropriate, and models have been introduced that relax the assumption of multivariate normality to an assumption that the resulting distribution is elliptically symmetric (or elliptically contoured or, simply, elliptical), still maintaining the MN model as a special case [7], [8]. This latter assumption is more biologically realistic, given that the distribution of allometric data is often heavy-tailed [9]. A well-known and simple model satisfying these requirements is the MN scale mixture (also known as MN variance mixture [10]); it is a finite/continuous mixture of MN distributions, with the same mean μ , where Σ is scaled by a convenient positive mixing random variable W whose probability mass/density function depends on a parameter vector θ governing the tails behavior (in terms of leptokurtosis) of the unconditional mixture.

As for the second aim, which is more debated in allometry [11], most multivariate studies have used principal component analysis (PCA); in this direction, in the important work by [12] it is suggested to summarize the morphometric variation by the first principal component (PC) of the log-transformed measurements. In the special bivariate case, the problem simplifies in fitting a line to the scatter plot of the log-transformed measurements and making inference about the fitted line, typically testing if the slope equals a given value and constructing confidence intervals for slope and elevation (intercept). Commonly used line-fitting methods in allometry are linear regression, major axis (which is equivalent to the first PC), and standardised major axis. These are least squares methods differing in the direction in which errors from the line are measured [6], [11]. Advantageously, however, major axis methods have a natural extension to more than two dimensions ([6], page 278).

The present paper has four merits. Firstly, concerning the first aim above, by using allometric data we show the existence of situations where, although the joint distribution of the log-transformed morphometric variables can be considered as symmetric and with heavy tails, the assumption of elliptical contours may be rather restrictive, with heavy-tailed distributions accounting for alternative symmetric shapes providing a better fit. To describe these shapes, we consider the family of multiple scaled symmetric (MSS) distributions proposed by Forbes and Wraith [13]. MSS distributions can be roughly considered as an extension of the MN scale mixture based on two key elements: the decomposition of the scale matrix Σ by eigenvalues and eigenvectors matrices Λ and Γ, and the use of the mixing random variable W separately for each dimension (PC) of the space spanned by the columns of Γ. The result is a multivariate heavy-tailed distribution in which the parameter vector θ is allowed to vary across PCs and where the contours are not constrained to be elliptical, although the symmetry is preserved. In allometric terms, these distributions have the advantage to be directly defined in the space of interest, namely the PC space, with each PC being distributed according to a (univariate) normal scale mixture. Moreover, regardless of the mixing distribution considered, and in analogy with the MN scale mixtures, MSS distributions enable (among other) natural/efficient computation of the maximum likelihood (ML) estimates via the expectation-maximization (EM) algorithm and robust ML estimation of μ and Λ based on down-weighting of mild outliers in the PC space. As a second merit of the paper, we propose a new member of the family of MSS distributions obtained by choosing a convenient shifted exponential as mixing distribution, so extending the recently proposed multivariate shifted exponential normal distribution [14] to the MSS framework. As a third merit, we introduce parsimony in the MSS distributions by constraining θ to be tied across PCs, resulting in distributions for the PCs having the same degree of leptokurtosis. In connection with the second aim discussed above, we know that line-fitting methods are quite sensitive to the presence of outliers, in the sense that the fitted line is directed towards them. As a limit case this implies that in certain situations, based on the position of the outlier(s), the resulting fitted line may be flipped into a completely different direction. As a final merit of the paper, via a simulation study and a sensitivity analysis based on allometric bivariate data, we show that the fitted line from MSS distributions, which is a major axis approach, is: 1) more robust to outliers when compared to the classical line-fitting methods discussed above, and 2) comparable with robust versions of major axis and standardised major axis approaches based on Huber’s M-estimator in place of least squares [15].

The paper is organized as follows. In Section 2 we first recall the MN scale mixtures (Section 2.1) and MSS distributions (Section 2.2), and then we introduce the multiple scaled shifted exponential normal distribution and the parsimonious variants of the MSS distributions (Section 2.3). In Section 3 we focus on the bivariate case by introducing what can be defined as major axis under MSS distributions. To simplify the reading of the paper, we postpone in Appendix A inferential aspects such as ML parameter estimation via the use of the EM algorithm, computation of standard errors for the estimates, formulation of confidence intervals and statistical tests. In Section 4 we describe the results of a simulation study aiming to compare, with and without outliers, the performance of the major axes from our MSS models with the lines fitted by well-established (robust) methods in allometry. In Section 5 we analyze two allometric data sets to demonstrate the better performance of the MSS models over heavy-tailed elliptical distributions which are well-known in the literature and to investigate the robustness of our MSS-based line-fitting methods. We conclude the paper with a brief discussion, as well as with avenues for further research, in Section 6.

2 Methodology

2.3 Special cases of MSS distributions and parsimony

In the MSS family of distributions, Forbes and Wraith [13] introduce the multiple scaled t (MSt) distribution. It is obtained by setting W h ∼ G ( θ h / 2 , θ h / 2 ) . For the MSt, the variance and kurtosis factors for each PC are respectively given by

(8) v ( θ h ) = θ h θ h − 2   a n d   k ( θ h ) = θ h − 2 θ h − 4 ,

being in general Var (Y _h ) = λ _h  v (θ _h ) and Kurt (Y _h ) = 3k (θ _h ), h = 1, …, d; compare with Eqs. (2) and (3).

In this paper we introduce a new member of the MSS family, the multiple scaled SEN (MSSEN) distribution, obtained by assuming W h ∼ S E ( θ h ) . For the MSSEN, the variance and kurtosis factors for each PC are respectively given by

(9) v ( θ h ) = θ h   exp ( θ h ) φ − 1 ( θ h )   a n d   k ( θ h ) = θ h v ( θ h ) 2 [ 1 − v ( θ h ) ] ,

where

φ m ( z ) = E − m ( z )       = ∫ 1 ∞ t m   exp ( − z t ) d t        = z − ( m + 1 ) Γ ( m + 1 , z )

is the Misra function [26], further generalization of the generalized exponential integral function E n ( z ) = ∫ 1 ∞ t − n exp ( − z t ) d t [27] with Γ ( n , z ) = ∫ z ∞ t n − 1 exp ( − t ) d t denoting the upper incomplete gamma function; for further details, see [14]. The MSt and MSSEN are the two special cases of MSS distributions we will consider in this paper.

To evaluate the data configurations favoring the MSt over the MSSEN, and vice versa, it is enough to consider univariate t and SEN distributions. Their main difference relies on the joint behavior of the variance and kurtosis factors v(θ _h ) and k(θ _h ), h = 1, …, d; refer to Figure 1. In particular, λ _h kept fixed, for the t distribution larger values of kurtosis are related to larger values of variance (cf. Figure 1(a)) while, for the SEN distribution, larger values of kurtosis are related to smaller values of variance (cf. Figure 1(b)); see [28].

Figure 1:

Multiplicative factors v (θ _h ) (solid line) and k (θ _h ) (dashed line).

In models (6) and (7) there is a θ _h -parameter for each PC. This means that, as the dimension d grows, the total number of free parameters to be estimated can quickly become large leading to overfitting. To introduce parsimony, we allow the d tailedness parameters to be tied across PCs, i.e. we assume θ ₁ = ⋯ =  θ _d . This results in normal scale mixtures, for the PCs, having the same degree of leptokurtosis. When this constraint is applied to MSt and MSSEN distributions, we will simply write C-MSt and C-MSSEN, respectively.

3 Line-fitting from bivariate MSS distributions

As emphasized in Section 1, allometry mainly focuses on the bivariate case (d = 2), where X = ( X 1 , X 2 ) ⊤ . In this case, we would like to make inference about the elevation (intercept) α ∈ ℝ and the slope β ∈ ℝ of the line describing the relationship between X₁ and X₂. Operationally, we are often interested in testing if β equals a particular value, and estimating confidence intervals for α and β. Therefore, it is important: 1) to evaluate how the MSS distributions specialize in the bivariate case, and 2) parametrize their pdf in terms of α and β.

First of all, it is important to underline that the line determined under MSS distributions belongs to the class of the major axis (or first PC axis) methods; as such, the line we are interested in is the axis of Y₁.

There are different ways to parametrize the MSS distributions, in the bivariate case, in terms of α and β. We decided to express Γ and μ₂ (mean of X₂) in terms of α and β in the following way

where μ₁ is the mean of X₁. According to this choice, the pdf of the MSt distribution becomes

while the pdf of the MSSEN distribution becomes

where Ψ = (α, β, μ₁, λ₁, λ₂, θ₁, θ₂)′.

Figure 2 shows, via isodensities, some possible shapes of the bivariate C-MSt and C-MSSEN distributions by fixing α, β, and μ₁ (it is equivalent to fix the mean vector μ and orientation matrix Γ), and by varying λ₁, λ₂, and the common θ so to produce PCs with given variances and kurtoses; refer to Section 2.3. In detail, in all the plots, μ₁ = 0, α = 0, and β = 1; note that in this case μ₂ = 0. Variances and kurtoses of the two PCs, along with the values of λ₁, λ ₂ , and θ, are given in Table 1. Figure 2 clearly shows how the shape of the MSS distributions is not constrained to be elliptical, although the symmetry is preserved. In particular, the choices made for θ (refer to Table 1) produce “smoothed” rhomboidal (Figure 2(a) and (b)) and starred (Figure 2(c) and (d)) contours.

Figure 2:

Examples of contour plots of bivariate MSS distributions with common θ where μ₁ = 0, α = 0, and β = 1. The values of λ _h , h = 1, 2, and θ are given in Table 1.

4 Simulation study: comparison of line-fitting methods

In this section, we investigate the behaviour of the fitted lines from our four bivariate MSS models – MSSEN, MSt, C-MSSEN, and C-MSt – through a large-scale simulation study performed using R [29]. We further provide a comparison with the following classical line-fitting methods: ordinary least-squares (OLS) regression, OLS bisector, major axis (MA; also known as orthogonal regression), standardized major axis (SMA; also known as reduced major-axis), robust MA (R-MA), and robust SMA (R-SMA), the last two approaches being based on Huber’s M-estimator; see [6], [11], [15] for details. Emphasis is given to the evaluation of the impact of outliers on the fitted slope of the competing methods.

To generate the data we consider a bivariate normal N 2 ( μ , Σ ) . We fix μ = 0 and Σ = ΓΛΓ ^⊤ , where

with λ∈(0, 1) governing the axis-ratio of the elliptical scatter: as λ approaches to 1, the ellipse becomes more and more similar to a circle while, on the other side, as λ approaches to 0, the ellipse degenerates on its major axis. With the notation of Section 3, this would be equivalent to say that we generate data from a bivariate normal distribution whose elevation, slope, and mean of X₁ are α = 0, β = 1, and μ₁ = 0, respectively. If we want outliers to be included in the sample, we generate them from a bivariate normal N 2 ( μ out , Σ out ) , where μ _out = q ⋅ (1, − 1)^⊤ and Σ_out = 0.01 ⋅ Σ, with q being a quantity governing the distance of the outliers from the bulk of the data (represented by μ  = 0).

We consider four experimental factors: the sample size (n ∈ {100, 200}), the proportion of simulated outliers (p ∈ {0, 0.05, 0.10}), λ ∈ {0.2, 0.3}, and q ∈ {8, 10}. In Figure 3 we show an example of generated dataset, with outliers, in the case n = 200, p = 0.05, λ = 0.3, and q = 8. The contamination scheme is not chosen randomly, but it is similar to a classical contamination scheme used to evaluate robustness for MA methods (see, e.g., [30], [31], [32], [33]). We generate 100 data sets for each combination of (n, λ) when p = 0, and for each combination of (n, p, λ, q) when p > 0; the difference is due to the fact that, when p = 0, there is no reason to consider the experimental factor q. This yields a total of 20 scenarios (2² = 4 in the case p = 0, and 2⁴ = 16 otherwise) and 2000 generated data sets (100 × 2² = 400 in the case p = 0, and 100 × 2⁴ = 1600 for p > 0). For each generated data set, we fit the 10 competing models/lines discussed above. The MA, SMA, R-MA and R-SMA methods are estimated using the smatr package for R [34].

Figure 3:

Example of simulated data, with outliers, in the case n = 200, p = 0.05, λ = 0.3, and q = 8.

For comparison’s sake, for each scenario and each line-fitting method, we display the box-plots of the estimated slopes over the 100 replications. While in the 4 scenarios without outliers we expect all the methods working well and comparably, in the 16 scenarios with outliers we hope the estimated slope is not flipped into the opposite direction.

Figure 4 reports the obtained results for the scenarios without outliers. A horizontal dashed gray line is placed at β = 1 to have, as a benchmark, the true but unknown value of β in the data generating process (DGP). Note that, we use a larger than expected range for the y-axis simply to have the same range for all the box-plots of this simulation study (compare with Figures 4 and 5); this should make easier the comparison of results referred to different scenarios. Apart from the OLS regression, the competing methods perform comparably, with a variability of the estimates that decreases either as n grows (from 100 to 200) or as λ decreases (from 0.3 to 0.2). The different behavior of the OLS estimates is due to the fact that the DGP violates the zero conditional mean assumption of the OLS regression (see, e.g., [35], Chapter 2).

Figure 4:

Box-plots, over 100 replications, of the estimated slopes in the scenarios without outliers (p = 0). Each plot refers to a different pair (n, λ).

Figure 5:

Box-plots, over 100 replications, of the estimated slopes in the scenarios with q = 8. Each plot refers to a different triplet (n, λ, p).

In Figures 5 and 6 we report the obtained results for the scenarios with q = 8 and q = 10, respectively. Regardless of the considered scenario, the estimated slopes from our MSS distributions, as those from R-MA and R-SMA, are robust to the presence of outliers, with a behavior very similar to the one analyzed in the scenarios without contamination (compare with Figure 4); similarity is also in terms of variability of the estimates which is decreasing when either n increases or λ decreases. The only difference between the estimated slopes from our MSS distributions, and those from R-MA and R-SMA, is a slight larger variability for the estimates from the latter methods, especially when λ grows. Interestingly enough, the proportion of outliers p does not seem to have an impact on the estimates behavior. Instead, the bisector, MA, and SMA methods of line-fitting are strongly affected by outliers, with a median estimated β (over the 100 replications) being always close to −1, the opposite of the true β. The OLS regression provides box-plots that are mainly located on values of β between 0 and −1. Although the OLS-estimates are also in this case very far from the true value β = 1, they are slightly better than those from bisector, MA, and SMA methods.

Figure 6:

Box-plots, over 100 replications, of the estimated slopes in the scenarios with q = 10. Each plot refers to a different triplet (n, λ, p).

5 Real data analyses

In this section we illustrate two applications of our models on real allometric data. The first one (Section 5.1) considers bivariate data; one of the aims is to evaluate how a single outlier affects the behavior of the competing line-fitting methods already considered in the simulation study of Section 4. The second analysis (Section 5.2) aims to fill one of the gaps in Section 4, showing an application of our models on data with more than two measurements.

According to the considerations of Section 1, both the analyses fill another gap with the simulation study of Section 4, i.e. the evaluation of the ability of our models to reproduce the joint distribution of the log-transformed morphometric variables. We also compare our models with the following multivariate elliptical (heavy-tailed) distributions: Mt, MTIN, MSEN, MLN, MCN, MSVG, MSH, SNIG, and MSGH given in Section 2.1. With MLN we denote the multivariate leptokurtic normal distribution introduced by [36].

In analogy with Section 4, the analysis is entirely conducted in R. Parameters are estimated based on the ML approach. We adopt the EM algorithm, illustrated in Appendix A.1, for MSS distributions. For them, we start the algorithm from the M-step, providing initial values for w ˙ i h and Γ ˙ (refer to Appendix A.1.2). In detail, we fix Γ ˙ as the eigenvector matrix of the sample covariance matrix, so that Γ ˙ does not depend on the considered distribution, while we differentiate w ˙ i h according to the considered model: we fix w ˙ i h = 1 for MSt and C-MSt, and w ˙ i h = 2 for MSSEN and C-MSSEN, i = 1, …, n and h = 1, …, d. The code for density evaluation, random number generation, and fitting for the considered MSS distributions is available at http://docenti.unict.it/punzo/Rcode.htm . As concerns the competing models, we use the WML.MLN() function, of the code available at http://docenti.unict.it/punzo/Rcode.htm , to fit the MLN, the mtin.ML.ECME() function of the mtin package^[1] to fit the MTIN, the msen.ML.EM() function of the msen package¹ to fit the MSEN, the CNmixt() function of the ContaminatedMixt package [37] to fit the MCN, and the ghyp package [38] to fit the remaining models.

Having the competing models a differing number of parameters, their goodness-of-fit comparison is made, as usual, via the Akaike information criterion (AIC) [39] and the Bayesian information criterion (BIC) [40] that, in our formulation, need to be maximized. Moreover, we use likelihood-ratio (LR) tests not only to make inference on the parameters (refer to Appendix A.3) but also to: 1) compare the MN distribution (null model) with each of the competing distributions (alternative models), and 2) for testing MSS-parsimony (refer to the last part of Section 2.3), that is to compare the constrained (null) and unconstrained (alternative) MSS models based on the same mixing random variable. In the sequel, the obtained p-value will be always compared to the classical 5% significance level.

5.1 Anthropometric measurements of male twins

The first analysis considers the m.twins data set accompanying the Flury package [41] for R. The data set contains six anthropometric measurements for n = 89 pairs of monozygotic/dizygotic male twins measured in the 1950s at the University of Hamburg, Germany. For details on these data, see [42]. We focus on the logarithm of d = 2 of the available measurements: the chest circumference of first and second twin. The variables are expressed in centimeters. The scatter plot of the log-transformed data is displayed in Figure 7.

Figure 7:

Scatter plot of the of the log-transformed m.twins data.

The first aim of this analysis is to evaluate the best model for the bivariate distribution of the data. The scatter in Figure 7 seems to be symmetric enough and this is corroborated by the Mardia test of symmetry (p-value = 0.444), as implemented by the mvn() function of the MVN package [43]. Table 2 presents the model comparison. AIC and BIC provide the same ranking for the best four models. They are all MSS distributions, with the MSSEN being the overall best, followed by the other unconstrained model (MSt). The top-four models have the same ranking induced by the p-values of the LR test of bivariate normality. In particular, the null hypothesis of bivariate normality is rejected for the unconstrained MSS models only (0.003 for the MSSEN and 0.004 for the MSt). The better performance of the unconstrained MSS models is also corroborated by the p-values from the LR test of MSS-parsimony given in the last column of Table 2: the MSSEN should be preferred to the C-MSSEN (p-value = 0.003), and the MSt should be preferred to the C-MSt (p-value = 0.007). Figure 8 shows the scatter plot of the log-transformed data where isodensities of the selected MSSEN model are superimposed. Summarizing, the results in Table 2, as well as the plot in Figure 8, underline the need for types of symmetry other than elliptical.

Table 2:

Model comparison, for the m.twins data, in terms of number of parameters (#par), log-likelihood, AIC, BIC, and p-values from the LR tests of multivariate normality and MSS-parsimony. Rankings induced by the considered criteria are also displayed.

Model	# par.	Log-Lik.	AIC	Ranking	BIC	Ranking	Normality LR p-value	Ranking	MSS-parsimony LR p-value
MN	5	213.822	417.644	12	405.201	5	–	–	–
Mt	6	215.653	419.306	6	404.374	6	0.176	5	–
MTIN	6	215.594	419.188	7	404.256	7	0.183	6	–
MSEN	6	214.711	417.421	13	402.489	12	0.346	12	–
MLN	6	214.978	417.957	11	403.025	11	0.282	11	–
MCN	7	216.943	419.885	5	402.465	13	0.210	7	–
MSVG	6	215.093	418.186	10	403.254	10	0.260	10	–
MSH	6	215.264	418.529	9	403.597	9	0.230	9	–
SNIG	6	215.360	418.721	8	403.789	8	0.215	8	–
MSGH	7	215.653	417.306	14	399.886	14	0.400	13	–
MSSEN	7	225.347	436.693	1	419.273	1	0.003	1	0.003
MSt	7	224.689	435.377	2	417.957	2	0.004	2	0.007
C-MSSEN	6	216.552	421.103	4	406.172	4	0.099	4	–
C-MSt	6	217.492	422.983	3	408.051	3	0.055	3	–

Figure 8:

Scatter plot of the log-transformed m.twins data with superimposed isodensities from the ML-fitted MSSEN model.

The second aim of the analysis is to compare the behavior of the line-fitting methods considered in Section 4. Table 3 reports the estimated elevation α ˆ and slope β ˆ , and approximate 95% confidence intervals (CIs). For the bisector method and for all our models, these intervals are based on the asymptotic normality of the estimators and are computed as α ˆ ∓ 1.96 ⋅ se ˆ ( α ˆ ) and β ˆ ∓ 1.96 ⋅ se ˆ ( β ˆ ) ; see Appendix A.3 for details. Standard errors are computed as described in [11] for the bisector method, and using the score based-approach illustrated in Appendix A.2 for our models. For the remaining methods, approximate 95% CIs are those provided by the smatr package (see [15] for details). Moreover, because of interest in allometry, we also evaluate whether the chest circumferences of the first and second twin are directly proportional. If this happens, then the log-transformed variables will be linearly related with a slope of 1. Hence, we simply need to test the null hypothesis H₀ : β = 1. We do that by checking if the value of the slope under the null is included in the approximate 95% CI. For our models only, we also show the p-values from the Wald and LR tests outlined in Appendix A.3.

Table 3:

Estimates and 95% confidence intervals (CIs) for the parameters α and β of the competing lines fitted to the m.twins data in Figure 7. The last two columns report the p-values of the Wald and LR tests for H₀:β = 1.

Methods		95% CI			95% CI		H0:β = 1 (p-value)
	α ˆ	2.5%	97.5%	β ˆ	2.5%	97.5%	Wald	LR
OLS	0.2583	−0.0494	0.5659	0.9398	0.8683	1.0113
MA	0.0089	−0.3177	0.3356	0.9978	0.9245	1.0768
SMA	0.0084	−0.2993	0.3161	0.9979	0.9289	1.0720
R-MA	0.0464	−0.2069	0.2996	0.9893	0.9321	1.0500
R-SMA	0.0449	−0.2003	0.2901	0.9896	0.9343	1.0483
Bisector	0.0084	−0.2468	0.2635	0.9979	0.9385	1.0573
MSSEN	0.0820	−0.1519	0.3159	0.9810	0.9263	1.0356	0.4942	0.2974
MSt	0.0733	−0.2043	0.3509	0.9831	0.9184	1.0477	0.6074	1.0000
C-MSSEN	0.0601	−0.2471	0.3673	0.9861	0.9147	1.0575	0.7019	0.6295
C-MSt	0.0706	−0.2275	0.3687	0.9837	0.9144	1.0529	0.6440	0.5766

As we can note by the obtained results, OLS regression works differently from the other competing methods, providing a larger α ˆ and, in line with the simulation results of Section 4, a smaller β ˆ . Moreover, 95% CIs for elevation and slope do not contain the values α = 0 and β = 1; these values would represent a direct proportionality between the chest circumferences of the first and second twin (β = 1) with unitary proportionality constant (α = 0). For the remaining methods we have the following considerations. They work comparably providing an estimated slope close to 1 and an estimated elevation close to 0. The 95% CIs for elevation and slope contain the reference values α = 0 and β = 1, respectively. As regards our models, the p-values from the Wald and LR tests are always greater than 0.05.

The final aim of the analysis is to evaluate how a single outlier affects the behavior of the competing line-fitting methods. We create a “perturbed” version of the m.twins data by adding a single outlier at ( x ‾ 1 − 1.6 , x ‾ 2 + 1.6 ) ⊤ , being x ‾ 1 and x ‾ 2 the sample means of the log-transformed chest circumferences of the first and second twin, respectively. The scatter plot of the perturbed log-transformed data is displayed in Figure 9.

Figure 9:

Scatter plot of the of the log-transformed m.twins data with the inclusion of an outlier displayed as •.

Table 4 reports the summary results about the estimated elevation and slope for each competing method, with corresponding lines displayed in Figure 10 for OLS, bisector, MA, and SMA, and in Figure 11 for the remaining approaches.

Table 4:

Estimates and 95% confidence intervals (CIs) for the parameters α and β of the competing lines fitted to the m.twins data with a single outlier (Figure 9). The last two columns report the p-values of the Wald and LR tests for H₀:β = 1.

Methods		95% CI			95% CI		H0:β = 1 (p-value)
	α ˆ	2.5%	97.5%	β ˆ	2.5%	97.5%	Wald	LR
OLS	5.6369	4.7729	6.5009	−0.3075	−0.5089	−0.1061
MA	8.5935	5.7935	11.3935	−0.9976	−2.1832	−0.4551
SMA	8.6006	7.7361	9.4651	−0.9993	−1.2208	−0.8179
R-MA	0.0462	−0.2177	0.3102	0.9895	0.9299	1.0528
R-SMA	0.0447	−0.2102	0.2997	0.9899	0.9324	1.0509
Bisector	8.6020	8.5509	8.6532	−0.9996	−1.0164	−0.9827
MSSEN	0.0824	−0.1499	0.3146	0.9809	0.9266	1.0351	0.4895	0.4962
MSt	0.0707	−0.2248	0.3662	0.9837	0.9148	1.0525	0.6418	1.0000
C-MSSEN	0.0826	−0.1427	0.3079	0.9808	0.9283	1.0334	0.4744	0.5505
C-MSt	0.0711	−0.2157	0.3579	0.9836	0.9168	1.0503	0.6293	1.0000

Figure 10:

Scatter plot of the of the log-transformed perturbed m.twins data and fitted lines from non-robust competing methods.

Figure 11:

Scatter plot of the of the log-transformed perturbed m.twins data and fitted lines from robust competing methods.

From these plots we realize that, while the lines from R-MA, R-SMA, and our MSS models are not affected at all by the outlier, the lines from the other methods are severely dragged towards that point, and this is especially true for the bisector, MA and SMA methods where the outlier lies on the fitted line. We can obtain similar insights by comparing Table 3 with Table 4.

5.2 Anthropometric measurements of female twins

The second analysis considers the f.twins data set always included in the Flury package. The data set contains d = 6 anthropometric measurements for n = 79 pairs of monozygotic/dizygotic female twins measured in the 1950s at the University of Hamburg, Germany. The measurements, expressed in centimeters, for each pair of twins are: stature (STA1 and STA2), hip width (HIP1 and HIP2), and chest circumference (CHE1 and CHE2) for each of the two sisters. For details on these data, see [42].

The matrix of pairwise scatter plots, for the log-transformed data, is displayed in Figure 12. These data can be considered as symmetric according to the Mardia test of symmetry (p-value = 0.205).

Figure 12:

Matrix of scatter plots of the log-transformed f.twins data.

Table 5 presents the model comparison. The best four models for the AIC are the MSS ones, with the MSSEN and C-MSSEN occupying the first and second position, respectively. As concerns the BIC, the best two models are the constrained MSS distributions (the C-MSSEN is the overall best), followed by the MSSEN. This podium is confirmed by LR test of multivariate normality, with the models on the podium being the only ones leading to the rejection of the null hypothesis. The best performance of the constrained MSS models, suggested by the BIC, is corroborated by the p-values from the LR test of MSS-parsimony: the C-MSSEN should be preferred to the MSSEN (p-value = 0.206), and the C-MSt should be preferred to the MSt (p-value = 0.640).

Table 5:

Model comparison, for the f.twins data, in terms of number of parameters (#par), log-likelihood, AIC, BIC, and p-values from the LR tests of multivariate normality and MSS-parsimony. Rankings induced by the considered criteria are also displayed.

Model	# par.	Log-Lik.	AIC	Ranking	BIC	Ranking	Normality LR p-value	Ranking	MSS-parsimony LR p-value
MN	27	763.103	1472.207	7	1408.232	4	–	–	–
Mt	28	763.675	1471.350	12	1405.006	11	0.450	11	–
MTIN	28	763.690	1471.380	11	1405.035	10	0.444	10	–
MSEN	28	764.802	1473.605	5	1407.260	5	0.192	5	–
MLN	28	763.845	1471.690	10	1405.345	9	0.389	9	–
MCN	29	764.286	1470.573	13	1401.859	13	0.554	12	–
MSVG	28	764.271	1472.543	6	1406.198	6	0.280	6	–
MSH	28	764.059	1472.119	8	1405.774	7	0.328	7	–
SNIG	28	763.900	1471.800	9	1405.455	8	0.372	8	–
MSGH	29	764.271	1470.543	14	1401.829	14	0.558	13	–
MSSEN	33	778.249	1490.499	1	1412.307	3	0.019	3	0.206
MSt	33	773.235	1480.470	4	1402.278	12	0.119	4	0.640
C-MSSEN	28	771.059	1486.118	2	1419.773	1	0.005	1	–
C-MSt	28	769.845	1483.691	3	1417.346	2	0.009	2	–

6 Conclusions and discussion

In allometric studies, the joint distribution of log-transformed morphometric variables is often assumed to be elliptical. Moreover, with reference to the bivariate case, we know the importance of using a good line-fitting method, preferably robust in the presence of outliers. Concerning the first point, by considering real allometric data, and by using for the first time multiple scaled symmetric (MSS) distributions in this context, we showed the existence of cases where the symmetry involved by the available data is different from the elliptical one. As further merits of the paper, we introduced a new member of this class of models, the multiple scaled shifted exponential normal distribution, and proposed a simple way to add parsimony for MSS distributions. Concerning the second issue, by using artificial and real data we showed that the fitted line from MSS distributions is: 1) more robust to outliers when compared to classical line-fitting methods considered in allometry, and 2) comparable with robust versions of major axis and standardized major axis approaches based on Huber’s M-estimator.

As an open point for further research, it could be interesting to consider MSS distributions in common principal component analysis (CPCA), generalization of standard PCA to several groups under the rigid mathematical assumption of equality of all principal components across groups. From an allometric point of view, this will allow us to extend our inferential results to the case where the user is interested in comparing several slopes or elevations, and testing for shift along the axis amongst several groups (refer to [6]).