Gaussian asymptotic limits for the $\alpha$-transformation in the analysis of compositional data

TL;DR

This paper investigates the asymptotic behavior of the $oldsymbol{ extalpha}$-transformation in compositional data analysis, focusing on the limiting case where the transformation parameter approaches zero, with theoretical and numerical insights.

Contribution

It provides a theoretical and numerical analysis of the $oldsymbol{ extalpha}$-transformation's limits in compositional data, especially when the parameter approaches zero, clarifying issues in maximum likelihood estimation.

Findings

As $oldsymbol{ extalpha}$ approaches zero, the distribution's behavior changes significantly.

Maximum likelihood estimates of $oldsymbol{ extalpha}$ near zero lead to large parameter estimates.

Theoretical results explain the observed numerical phenomena in the limiting case.

Abstract

Compositional data consists of vectors of proportions whose components sum to 1. Such vectors lie in the standard simplex, which is a manifold with boundary. One issue that has been rather controversial within the field of compositional data analysis is the choice of metric on the simplex. One popular possibility has been to use the metric implied by logtransforming the data, as proposed by Aitchison [1, 2]; and another popular approach has been to use the standard Euclidean metric inherited from the ambient space. Tsagris et al. [21] proposed a one-parameter family of power transformations, the $\alpha$-transformations, which include both the metric implied by Aitchison's transformation and the Euclidean metric as particular cases. Our underlying philosophy is that, with many datasets, it may make sense to use the data to help us determine a suitable metric. A related possibility is to apply the $\alpha$-transformations to a parametric family of distributions, and then estimate a along with the other parameters. However, as we shall see, when one follows this last approach with the Dirichlet family, some care is needed in a certain limiting case which arises $(\alpha \neq 0)$, as we found out when fitting this model to real and simulated data. Specifically, when the maximum likelihood estimator of a is close to 0, the other parameters tend to be large. The main purpose of the paper is to study this limiting case both theoretically and numerically and to provide insight into these numerical findings.

Tables

Table 1. Table 1: Estimates of the parameters for the Mammals dataset. The columns of the table give estimates of parameters b j subscript 𝑏 𝑗 b_{j} defined in ( 8 ), but labelled by the component name rather than j 𝑗 j . The optimal value for α 𝛼 \alpha is 0.06 0.06 0.06 .

Method	water	protein	fat	lactose	ash
Direct MLE	793.908	671.054	680.629	663.055	604.305
Asymptotic 1	782.297	668.362	677.964	660.235	597.352
Asymptotic 2	793.673	670.856	680.428	662.859	604.127

Table 2. Table 2: Estimates for the East Bay Clams dataset. The columns of the table give estimates of parameters b j subscript 𝑏 𝑗 b_{j} defined in ( 8 ), but labelled by component name rather than j 𝑗 j . The optimal value for α 𝛼 \alpha is 0.28 0.28 0.28 .

Method	dl	dm	ds
Direct MLE	232.218	224.377	222.197
Asymptotic 1	231.568	223.800	221.592
Asymptotic 2	231.990	224.156	221.979

Table 3. Table 3: Estimates for the OECD dataset. The columns of the table give estimates of the parameters b j subscript 𝑏 𝑗 b_{j} defined in ( 8 ), but labelled by component name rather than j 𝑗 j . The optimal value for α 𝛼 \alpha is 0.14 0.14 0.14 .

Method	PCINC	AGR	IND	SER
Direct MLE	212.965	129.782	143.360	136.907
Asymptotic 1	199.034	124.925	139.824	132.928
Asymptotic 2	212.669	129.601	143.161	136.716

Table 4. Table 4: Estimates for the GRTA dataset. The columns of the table give estimates of the parameters b j subscript 𝑏 𝑗 b_{j} defined in ( 8 ), but labelled by component name rather than j 𝑗 j . The optimal value for α 𝛼 \alpha is − 0.04 0.04 -0.04 .

Method	Killed	Seriously injured	Slightly injured
Direct MLE	28366.65	28109.32	25420.99
Asymptotic 1	28305.94	28057.80	25320.62
Asymptotic 2	28366.36	28109.03	25420.73

Equations

$$\mathcal{S}\text{im}(d)=\{\mathbf{x}=(x_1,\dots ,x_D{)}^{\top }\in {\mathbb{R}}^D:x_j\ge 0,j=1,\dots ,D;\sum _{j=1}^D{x}_j=1\},$$

$$w_j=\log (x_j)-D^{-1}\sum _{k=1}^D\log (x_k),j=1,\dots ,D.$$

$${\mathbf{u}}_{\alpha }(\mathbf{x})={\left(\frac{x_1^{\alpha }}{{\sum }_{k=1}^D{x}_k^{\alpha }},\dots ,\frac{x_D^{\alpha }}{{\sum }_{k=1}^D{x}_k^{\alpha }}\right)}^T.$$

$${\mathbf{u}}_{\alpha }^{-1}(\mathbf{x})=\left(\frac{x_1^{1/\alpha }}{{\sum }_{k=1}^D{x}_k^{1/\alpha }},\dots ,\frac{x_D^{1/\alpha }}{{\sum }_{k=1}^D{x}_k^{1/\alpha }}\right).$$

$$\underset{\alpha \to 0}{\lim }{\alpha }^{-1}\left(D{\mathbf{u}}_{\alpha }(\mathbf{x})-1_D\right)=\mathbf{w}(\mathbf{x}),$$

$${\Delta }_{\alpha }\left(\mathbf{x},\mathbf{y}\right)=\frac{D}{|\alpha |}\Vert {\mathbf{u}}_{\alpha }(\mathbf{x})-{\mathbf{u}}_{\alpha }(\mathbf{y})\Vert =\frac{D}{\left|\alpha \right|}{\left[\sum _{j=1}^D{\left(\frac{x_j^{\alpha }}{{\sum }_{k=1}^D{x}_k^{\alpha }}-\frac{y_j^{\alpha }}{{\sum }_{k=1}^D{y}_k^{\alpha }}\right)}^2\right]}^{1/2}.$$

$${\Delta }_0\left(\mathbf{x},\mathbf{y}\right)=\underset{\alpha \to 0}{\lim }{\Delta }_{\alpha }\left(\mathbf{x},\mathbf{y}\right)={\left[\sum _{j=1}^D{\left(\log \frac{x_j}{{\prod }_{k=1}^D{x}_k^{1/D}}-\log \frac{y_j}{{\prod }_{k=1}^D{y}_k^{1/D}}\right)}^2\right]}^{1/2},$$

$${\Delta }_1\left(\mathbf{x},\mathbf{y}\right)=D{\left[\sum _{j=1}^D{\left(x_j-y_j\right)}^2\right]}^{1/2}$$

$$\left|{\mathbf{J}}_{\alpha }\right|=|\alpha {|}^d\prod _{j=1}^D\frac{x_j^{\alpha -1}}{{\sum }_{k=1}^D{x}_k^{\alpha }}.$$

$$f_{\mathbf{\gamma }}(\mathbf{v})=\frac{\Gamma ({\gamma }_{+})}{{\prod }_{j=1}^D\Gamma ({\gamma }_j)}\prod _{j=1}^D{v}_j^{{\gamma }_j-1},$$

$${\mathbf{u}}_i=(u_{i1},\dots ,u_{iD}{)}^{\top }\sim \text{Dirichlet}(b_1,\dots ,b_D),i=1,\dots ,n,$$

$$u_{ij}=\frac{x_{ij}^{\alpha }}{{\sum }_{k=1}^D{x}_{ik}^{\alpha }},i=1,\dots ,n,j=1,\dots ,D.$$

$$\bar{l}(\alpha,b_{1},\ldots,b_{D})=\sum_{i=1}^{n}\big{[}\log f_{\bf b}({\bf u}_{i})+\log\left|{\bf J}_{\alpha,i}\right|\big{]}$$

$$\bar{l}(\alpha ,b_1,\dots ,b_D)$$

$$+\sum _{i=1}^n\sum _{j=1}^D(\alpha b_j-1)\log (x_{ij})-b_{+}\sum _{i=1}^n\log \left(\sum _{j=1}^D{x}_{ij}^{\alpha }\right),$$

$$y_{ij}=\log (x_{ij}),i=1,\dots ,n;j=1,\dots ,D.$$

$$b_j=\frac{b}{{\alpha }^2}(1+\alpha c_j),j=1,\dots ,D,$$

$$\sum _{j=1}^D{c}_j=0,$$

$$\ell (\alpha ,b,c_1,\dots ,c_D)=\bar{\ell }\left(\alpha ,\frac{b}{{\alpha }^2}(1+\alpha c_1),\dots ,\frac{b}{{\alpha }^2}(1+\alpha c_D)\right).$$

$$\ell (\alpha ,b,c_1,\dots ,c_D)$$

$$+n{\mathcal{C}}_0+\alpha \sum _{i=1}^n{\mathcal{C}}_{1,i}+O({\alpha }^2),$$

$${\mathcal{C}}_0=-\frac{1}{2}\log (D)-D{\bar{y}}_{++},{\mathcal{C}}_{1,i}=-\frac{b}{6}\left(D{\hat{\kappa }}_{3,i}+\sum _{j=1}^D{c}_j^3\right),$$

$${\hat{\kappa }}_{1,i}\equiv {\bar{y}}_{1,i}=D^{-1}\sum _{j=1}^D{y}_{ij},{\hat{\kappa }}_{2,i}=\left(D^{-1}\sum _{j=1}^D{y}_{ij}^2\right)-{\hat{\kappa }}_{1,i}^2$$

$${\hat{\kappa }}_{3,i}=\left(D^{-1}\sum _{j=1}^D{y}_{ij}^3\right)-3{\hat{\kappa }}_{1,i}\left(D^{-1}\sum _{j=1}^D{y}_{ij}^2\right)+2{\hat{\kappa }}_{1,i}^3$$

$${\hat{c}}_j=n^{-1}\sum _{i=1}^n(y_{ij}-{\bar{y}}_{i+})+O(\alpha )={\bar{y}}_{+j}-{\bar{y}}_{++}+O(\alpha ),1\le j\le D,$$

$$\hat{b}={\left\{\frac{1}{nd}\sum _{i=1}^n\sum _{j=1}^D{\left(y_{ij}-{\bar{y}}_{+j}-{\bar{y}}_{i+}+{\bar{y}}_{++}\right)}^2\right\}}^{-1}+O(\alpha )$$

$${\mathbf{\Sigma }}^{-}=\frac{1}{{\sigma }^2}\left({\mathbf{I}}_D-D^{-1}{1}_D{1}_D^{\top }\right)$$

$$(Z_1-\bar{Z}-{\mu }_1,\dots ,Z_D-\bar{Z}-{\mu }_D){\mathbf{\Sigma }}^{-}(Z_1-\bar{Z}-{\mu }_1,\dots ,Z_D-\bar{Z}-{\mu }_D{)}^{\top }$$

$$=\frac{1}{{\sigma }^2}\sum _{j=1}^D(Z_j-\bar{Z}-{\mu }_j{)}^2-\frac{1}{D{\sigma }^2}{\left(\sum _{j=1}^D(Z_j-\bar{Z}-{\mu }_j)\right)}^2$$

$$=\frac{1}{{\sigma }^2}\sum _{j=1}^D(Z_j-\bar{Z}-{\mu }_j{)}^2,$$

$$b\sum _{j=1}^D{\left(y_{ij}-{\bar{y}}_{i+}-c_j\right)}^2,$$

$$b^{\ast }=b(1+\alpha {\bar{c}}_{+})and{c}_j^{\ast }=\frac{c_j-{\bar{c}}_{+}}{1+\alpha {\bar{c}}_{+}},j=1,\dots ,D.$$

Full text

Key Words: Dirichlet distribution; log-ratio transformation; manifold; metric; power transformation.

1 Introduction

A compositional data vector has non-negative components which sum to a constant, usually taken to be $1$, as we assume here. Such vectors lie in the standard simplex, $\mathcal{S}\textrm{im}(d)$, which has the form

[TABLE]

where $d=D-1$. Compositional datasets appear naturally in many fields including medicine, geology, archeology, biology and economics; see e.g. [2].

The following two approaches have been widely used in compositional data analysis. One of these approaches, due to Aitchison [1, 2], is to apply the following transformation, known as the centred log-ratio transformation:

[TABLE]

Also, define ${\bf w}=(w_{1},\ldots,w_{D})^{\top}$. Note that, by definition, the components of ${\bf w}$ sum to [math]. Given a sample of compositional vectors ${\bf x}_{1},\ldots,{\bf x}_{n}$, we transform to ${\bf w}_{1},\ldots,{\bf w}_{n}$ using (2). Aitchison then proposed applying standard multivariate methods, such as principal components analysis, to log-transformed compositional vectors, which is tantamount to assuming the standard Euclidean metric is appropriate for the ${\bf w}_{i}$.

Another popular approach is to “do nothing”; i.e. assume that the Euclidean metric inherited from the ambient Euclidean space $\mathbb{R}^{D}$ is appropriate for the ${\bf x}_{i}$, and apply standard multivariate methods directly to the ${\bf x}_{i}$. For advocation of this approach see Baxter [3, 4], Baxter et al. [5] and Baxter and Freestone [6].

For some datasets, the log-ratio approach seems preferable, whereas for other datasets, the “do nothing” approach seems superior. Nevertheless, the choice over which method to use, if either, has occasionally become heated and acrimonious. See the paper by Scealy and Welsh [19] for a careful and detailed discussion of the issues and history of this debate. Our view concerning the choice of metric is that it is not possible to come up with a compelling choice for either method based on purely a priori or theoretical grounds, and that a more pragmatic approach is to make a data-dependent choice of metric.

Bearing this in mind, one possibility is to consider a finite-dimensional family of transformations, such as the one-dimensional $\alpha$-transformation family considered in this paper, where each transformation implies a choice of metric. Then use the data to choose the optimal transformation, which in turn determines an optimal metric. It is not so clear how to do this non-parametrically, but if a parametric family of models, such as the Dirichlet in the present context, is assigned to the transformed compositional data, then we can estimate both the distributional parameters and the transformation parameter simultaneously by maximum likelihood. Using this approach we estimate the parameters of the $\alpha$-transformed compositional data model. In effect, the optimal transformation parameter compels the data to resemble the chosen parametric family as closely as possible. However, after carrying out this approach with a number of real and simulated datasets we found that in cases where the maximum likelihood estimator of $\alpha$ is close to [math], there is a tendency for the other (Dirichlet) parameters to be rather large. See Tables 1-4 in Section 5 for examples of this phenomenon. This has motivated us to take a close look at the asymptotics as (i) $\alpha\rightarrow 0$ and (ii) the Dirichlet becomes highly concentrated. This theoretical and numerical investigation is the principal focus of our paper.

The phenomenon mentioned in the previous paragraph potentially arises with analogous families of metrics on other manifolds should we wish to select an optimal metric in a data-dependent way. Consequently, before focusing on compositional data in the remainder of the paper, we briefly discuss the question of choice of metric in the setting of general manifolds in Section 2.

The structure of this paper is as follows. In the next section we discuss the choice of metric in the statistical analysis of manifold-valued data from a broader perspective. Then, in Section 3, we briefly review the $\alpha$-transformation for compositional data. In Section 4 we present some asymptotics for small $\alpha$ when the $\alpha$-transformation is applied to the Dirichlet family, and in Section 5 we present numerical results, both from simulations and the analysis of real datasets. These numerical results reflect the theoretical results in the previous section. Proofs are given in Section 6.

2 Statistics on manifolds and choice of metric

Since the middle of the last century, statistical analysis of data which naturally lies in a manifold has grown enormously in interest and relevance. An important early example was R.A. Fisher’s paper [12] on directional statistics with application to paleomagnetic directional data. Subsequently, the field of directional statistics, in which the sample space consists of the surface of a unit sphere in Euclidean space (or, equivalently, the set of directions or unit vectors) has grown into a well-developed and mature field. See, for example, the books by Mardia [17], Mardia and Jupp [18] and Fisher et al. [13].

In a related but different direction, there has been great interest in statistical shape analysis from the 1980s onwards. A popular and mathematically deep approach to statistical shape analysis was initiated by D.G. Kendall; see the monograph Kendall et al. [16]. Other books on this field include Small [20] and Dryden and Mardia [10, 11]. In the Kendall approach, objects in Euclidean space, usually in two or three dimensions, are represented by the coordinates of a finite set of labelled landmarks. For example, a human face could be represented approximately by landmark coordinates at obvious places including the eyes, ears, nose and mouth. Then the shape is defined as what remains when translation, scale and orientation have been “quotiented out”. Further details may be found in the literature mentioned above.

Whenever analysing manifold data, one question that always needs to be addressed is what metric to use. For example, a choice of metric is required if we wish to define a Frechét mean on the manifold, which is a generalisation of the least-squares mean in Euclidean space. One important dichotomy is between intrinsic and extrinsic metrics. In the former, the metric, often a Riemannian metric, is defined “intrinsically”, without any reference to an ambient space. In contrast, with extrinsic metrics, the manifold of interest is considered to be embedded in an ambient space, typically a Euclidean space endowed with the standard Euclidean metric, and the manifold inherits the metric from the ambient space. From a conceptual point of view, intrinsic metrics are often to be preferred. However, extrinsic metrics are often easier to work with, both from a practical and theoretical point of view; from a practical point of view because intrinsic distance may be more difficult to calculate and from a theoretical point of view because e.g. asymptotic theory tends to be more difficult when intrinsic metrics are used. As an example of the difficulties which can arise, see Hotz and Huckemann [15], who study the non-standard behaviour which arises when the intrinsic, or arc-length, metric is used for circular data.

Even when we have made the choice to use an extrinsic metric, the choice of metric may still not be clear. Two examples (which supplement the case of compositional data considered here) in which a one-parameter family of metrics may be worth considering will now be mentioned.

In [8] a family of power metrics on the space of covariances matrices of a given dimension is considered, with a particular focus on diffusion tensor imaging. In this case, the relevant manifold is the space of positive-definite matrices of given dimension. There is a close analogy with the $\alpha$-transformation considered here; in their setting the powered quantities are the eigenvalues of the covariance matrices, whereas here the powered quantities are components of compositional vectors.

The second example relates to statistical shape analysis. In [9] a family of extrinsic metrics is considered for measuring the distance between two shapes described by labelled landmarks. Here, the family of metrics consists of powers of eigenvalues of certain symmetric non-negative definite matrices, but in this case the relevant matrices do not have full rank.

In both of these examples, the question of how to choose the metric arises, as it does in this paper, where the focus is on compositional data. It is clear that the choice one makes is going to have an effect on the statistical analysis, in some cases a major effect.

3 Review of the $\alpha$-transformation

For the purpose of this paper it is convenient to define the $\alpha$-transformation of a compositional vector $\mathbf{x}\in\mathcal{S}\textrm{im}(d)$ by

[TABLE]

This transformation, which has a slightly different but mathematically equivalent form to that given in [21] and [22], is also known as the compositional power transformation [2].

Note that the transformation ${\bf u}_{\alpha}$ given by (3) defines a bijection of the interior of $\mathcal{S}\textrm{im}(d)$ with inverse

[TABLE]

If one excludes the boundary of the simplex, which corresponds to compositional vectors that have one or more components equal to zero, then the $\alpha$-transformation (3) and its inverse (4) are well defined for all $\alpha\in\mathbb{R}$. The motivation for transformation (3) is that the case $\alpha=0$ corresponds to the log-ratio approach, whereas $\alpha=1$ corresponds to the "do nothing" approach; see the discussion in the Introduction. We define the case $\alpha=0$ in terms of the limit $\alpha\rightarrow 0$; then, rescaling by $\alpha$ as follows, we obtain

[TABLE]

where ${\bf 1}_{D}$ is the $D$-vector of ones and ${\bf w}(\mathbf{x})$ is the centred log-ratio transformation (2). For the case $\alpha=1$, (3) is just the identity transformation of the simplex.

The $\alpha$-transformation (3) leads naturally to a simplicial distance measure $\Delta_{\alpha}\left({\bf x},{\bf y}\right)$, which we call the $\alpha$-metric, between observations ${\bf x},{\bf y}\in\mathcal{S}\textrm{im}(d)$, defined in terms of the Euclidean distance $\|\cdot\|$ between transformed observations, i.e.,

[TABLE]

The special case

[TABLE]

is the Euclidean distance on the log-ratio transformed data, whereas

[TABLE]

is just Euclidean distance multiplied by $D$.

The choice of the value of $\alpha$ depends upon the context of the analysis. Maximum likelihood estimation, assuming the transformation (3) or (4) has been applied to a parametric family of distributions such as the Dirichlet, requires the Jacobian determinant of (3), which is given by (see [23] for the derivation)

[TABLE]

4 Gaussian asymptotic limits as $\alpha\rightarrow 0$

In this section, we examine the $\alpha$-transformation defined in (3), applied to a sample ${\bf x}_{1},\ldots,{\bf x}_{n}$ of compositional data on $\mathcal{S}\textrm{im}(d)$. The $\alpha$-transformation along with the assumption of a $D$-dimensional Dirichlet distribution for the transformed data result in a parametric family which has $D+1$ parameters in total (i.e. one additional parameter). The particular focus in this Section is to see what happens when we apply maximum likelihood estimation to the transformed $D+1$ parameter family in the limit where the Dirichlet becomes highly concentrated and Gaussian and $\alpha\rightarrow 0$. This question arose when fitting this transformed model to datasets where the MLE of $\alpha$ turns out to be small, while the MLE of the Dirichlet parameters tended to be large. For numerical evidence of this phenomenon, see Tables 1-4 below.

4.1 The Dirichlet case

Recall that the Dirichlet density on the simplex at a point ${\bf v}=(v_{1},\ldots,v_{D})^{\top}\in\mathcal{S}\textrm{im}(d)$ is

[TABLE]

where $\gamma_{+}=\sum_{j=1}^{D}\gamma_{j}$. Suppose that

[TABLE]

where

[TABLE]

Then the log-likelihood which results when the $\alpha$-transformation is applied to ${\bf x}_{1},...,{\bf x}_{n}$ is given by

[TABLE]

which can be rewritten after expansion and re-organization as

[TABLE]

where $b_{+}=\sum_{j=1}^{D}b_{j}$. Also, define

[TABLE]

We now state our first result about the asymptotic behavior of the log-likelihood.

Proposition 1. Define $y_{ij}$ as in (7) and

[TABLE]

where

[TABLE]

and define

[TABLE]

Then as $\alpha\rightarrow 0$, the log-likelihood $\ell(\alpha,b,c_{1},\ldots,c_{D})$ for $b$ and the $c_{j}$ satisfies

[TABLE]

where $\bar{y}_{i+}=D^{-1}\sum_{j=1}^{D}y_{ij}$,   $\bar{y}_{++}=(nD)^{-1}\sum_{i=1}^{n}\sum_{j=1}^{D}y_{ij}$,

[TABLE]

and

[TABLE]

and

[TABLE]

are the first three sample cumulants of the $y_{i1},\ldots,y_{iD}$.

Proposition 1 is proved in Section 6. We present several remarks below concerning Proposition 1.

Remark 1. When $\alpha$ is small, the MLE obtained by maximising (10) is given by

[TABLE]

where $\bar{y}_{+j}=n^{-1}\sum_{i=1}^{n}y_{ij}$, and

[TABLE]

Remark 2. The limit distribution is logistic normal with isotropic covariance structure (see Remark 3 below for clarification). The simple covariance structure in the limit as $\alpha\rightarrow 0$ is a consequence of the restricted nature of the covariance structure of the Dirichlet distribution. Note that the limit distribution has the same number of parameters as the $D$-dimensional Dirchlet distribution, namely $D$.

Remark 3. Although the likelihood in (10) clearly has a Gaussian character, it may not be immediately clear how it arises, given that when we sum $y_{ij}-\bar{y}_{i+}$ over $j$ we get [math]. To clarify, suppose we have independent Gaussian variables $Z_{j}\sim N(\mu_{j},\sigma^{2})$, $j=1,\ldots,D$, where it is assumed that $\sum_{j=1}^{D}\mu_{j}=0$. Define $\bar{Z}=D^{-1}\sum_{j=1}^{D}Z_{j}$ to be the sample mean. The covariance matrix of $(Z_{1}-\bar{Z},\ldots,Z_{D}-\bar{Z})^{\top}$ is ${\boldsymbol{\Sigma}}=\sigma^{2}\left({\bf I}_{D}-D^{-1}{\bf 1}_{D}{\bf 1}_{D}^{\top}\right)$, where ${\bf I}_{D}$ is the $D\times D$ identity matrix and ${\bf 1}_{D}$ is the $D$-vector of ones. The Moore-Penrose inverse of ${\boldsymbol{\Sigma}}$ is given by

[TABLE]

and the relevant quadratic form is given by

[TABLE]

because we have assumed that $\mu_{1}+\cdots+\mu_{D}=0$. Note that if we put $b=1/\sigma^{2}$, $c_{j}=\mu_{j}$ and $Z_{j}=y_{ij}$, $1\leq j\leq D$, then (16) is equal to

[TABLE]

which agrees with the inner sum of the second term in (10).

Remark 4. The assumption (9) only plays a cosmetic role in the proof provided that we are prepared to allow $b$ and the $c_{j}$ to depend weakly on $\alpha$. Specifically, suppose that $\bar{c}_{+}=D^{-1}\sum_{j=1}^{D}c_{j}\neq 0$. Then define

[TABLE]

It is easily checked that

[TABLE]

and also $b^{\ast}$ and the $c_{j}^{\ast}$ only depend weakly on $\alpha$ in the sense that $b^{\ast}\rightarrow b$ and $c_{j}^{\ast}\rightarrow c_{j}-\bar{c}_{+}$, $1\leq j\leq D$, as $\alpha\rightarrow 0$.

Remark 5. Note that the $b_{j}$ defined in (8) have a rather special structure in that they are asymptotically equal as $\alpha\rightarrow 0$. A natural question is: what happens to the asymptotics in Proposition 1 if the $b_{j}$ do not coalesce as $\alpha\rightarrow 0$, specifically if the $b_{j}$ are of the form $b_{j}=\gamma_{j}/\alpha^{2}$. It turns out that, in contrast to Proposition 1, when the $b_{j}$ do not coalesce the log-likelihood does not converge to a finite limit as $\alpha\rightarrow 0$, as we shall see later.

4.2 The role of asymptotic normality and the general case

The asymptotic regime considered in the previous section is quite limited because the $b_{j}$ in (8) satisfy

[TABLE]

However, it seems that we need this strong assumption to hold if we are to have convergence to a finite log-likelihood as indicated in Proposition 1. One important point we have not discussed yet is the fact that the Dirichlet distribution in Proposition 1 is converging to a multivariate normal. Consider the following elementary result.

Proposition 2. Suppose ${\bf u}^{(\alpha)}\sim\textrm{Dirchlet}({\boldsymbol{\gamma}}/\alpha^{2})$, where ${\boldsymbol{\gamma}}=(\gamma_{1},\ldots,\gamma_{D})^{\top}$. Then as $\alpha\rightarrow 0$,

[TABLE]

where ${\boldsymbol{\Sigma}}=\gamma_{+}^{-1}\left\{\textrm{diag}(\bar{\boldsymbol{\gamma}})-\bar{\boldsymbol{\gamma}}\bar{\boldsymbol{\gamma}}^{\top}\right\}$, $\bar{\boldsymbol{\gamma}}=\gamma_{+}^{-1}{\boldsymbol{\gamma}}$ and $\gamma_{+}=\sum_{j=1}^{D}\gamma_{j}$. *Moreover, the convergence of densities also occurs uniformly in growing balls of the form * ${\bf v}\in\mathcal{B}_{\alpha^{-1+\epsilon}}(0)$ as $\alpha\rightarrow 0$ for any fixed $\epsilon\in(0,1)$, in the sense that, on the set $\sum_{j=1}^{D}v_{j}=0$,

[TABLE]

uniformly for ${\bf v}\in\mathcal{B}_{\alpha^{-1+\epsilon}}(0)$, where in the line above, $f$ is the Dirichlet density.

The proof of Proposition 2 is similar to that of Proposition 1 but is somewhat simpler and so is omitted.

We now consider the asymptotic behaviour of $\bf y$ as $\alpha\rightarrow 0$ in the context of a Dirichlet distribution with parameters which do not coalesce as in Proposition 1.

Corollary 1. Under the assumptions of Proposition 2,

[TABLE]

where ${\bf y}=(y_{1},\ldots,y_{D})^{\top}$, $y_{j}=\log(x_{j})$ and $\bar{y}=D^{-1}\sum_{j=1}^{D}y_{j}$, $u_{j}=x_{j}^{\alpha}/\sum_{k=1}^{D}x_{k}^{\alpha}$ and ${\bf u}=(u_{1},\ldots,u_{D})^{\top}\sim\text{Dirichlet}(\gamma_{1}/\alpha^{2},\ldots,\gamma_{D}/\alpha^{2})$.

Proof of Corollary 1. If $u_{j}=x_{j}^{\alpha}/\sum_{k=1}^{D}x_{k}^{\alpha}$, then elementary calculations show that as $\alpha\rightarrow 0$,

[TABLE]

where $y_{j}=\log(x_{j})$ and $\bar{y}=D^{-1}\sum_{j=1}^{D}y_{j}$. Consequently, from (18),

[TABLE]

and so Corollary 1 follows because, from Proposition 2, ${\bf u}\approx\mathcal{N}_{D}(\bar{\boldsymbol{\gamma}},\alpha^{2}{\boldsymbol{\Sigma}})$. $\square$

So for ${\bf u}$ to attain a general limiting mean vector $\bar{\boldsymbol{\gamma}}$ in the limit as $\alpha\rightarrow 0$, it is necessary for the expectation vector of ${\bf y}-\bar{y}{\bf 1}_{D}$ to go to infinity as $\alpha\rightarrow 0$.

Finally, we go to a more general case in which we do not assume that ${\bf u}\equiv{\bf u}^{(\alpha)}$ is Dirichlet but just assume that as $\alpha\rightarrow 0$,

[TABLE]

Since the ${\bf u}$ are defined on the simplex $\mathcal{S}\textrm{im}(d)$ the vector ${\boldsymbol{\eta}}=(\eta_{1},\ldots,\eta_{D})^{\top}$ has non-negative components which satisfy $\sum_{j=1}^{D}\eta_{j}=1$ and the covariance matrix ${\boldsymbol{\Omega}}$ has an eigenvector ${\bf 1}_{D}$ with corresponding eigenvalue [math]. In this case, we obtain the same limiting result as in (17), but with the more general covariance matrix ${\boldsymbol{\Omega}}$ replacing ${\boldsymbol{\Sigma}}$, and with the same requirement that the mean vector of ${\bf y}-\bar{y}{\bf 1}_{D}$ goes to infinity if we wish for ${\bf u}$ to obtain a general limiting mean vector.

5 Numerical results

We first show the asymptotic results of the previous section on synthetic data and then apply the $\alpha$-transformation on real datasets.

5.1 Numerical Simulations

As discussed earlier, a finite log-likelihood is obtained under the assumption in (8) for $\bf{b}$. In contrast, when the Dirichlet’s parameters do not converge to the same value as $\alpha\to 0$ the log-likelihood is expected to be infinite. Figure 1 presents this behaviour as it is encoded by the asymptotic limit for the mean value ${\bf y}-\bar{y}$. Blue lines correspond to the former case with $b=1$ and ${\bf c}=[0.1,0.3,-0.4]^{T}$. As expected from Proposition 1, the mean value of ${\bf y}-\bar{y}$ converges to a finite value which is however different at each coordinate due to the different values of ${\bf c}$’s elements. Red lines correspond to the case with ${\bf b}=[1.1,1.3,0.6]^{T}$ where the mean value diverges as implied by Corollary 1. Evidently, the rate of divergence is inverse proportional to $\alpha$.

Remark 6. The estimation of ${\bf y}$ requires the exponentiation of ${\bf u}$ with $1/\alpha$ which leads to numerical instability in the general case (i.e., the red lines) when $\alpha$ is below $0.001$. In order to ensure numerical stability, we calculate ${\bf y}$ with the following equivalent formula

[TABLE]

where $j^{*}=\operatorname*{\arg\!\max}_{j}b_{j}$.

5.2 Real data examples

Three real datasets are used for illustration purposes.

• •

Mammals. This dataset is taken from Hartigan [14] and describe the percentage composition of 24 mammals’ milk on the basis of 5 different constituents (water, protein, fat, lactose and ash). A 4-group solution has been considered as the optimal grouping of such data.

• •

East Bay Clams. From the many colonies of clams in East Bay, 20 were selected at random and from each a sample of clams was taken [2]. Each sample was sieved into three size ranges, large, medium, small; then each size range was sorted by shell colour, dark, light. For each colony the proportions of clams in each colour-size combination was estimated and the corresponding compositions are recorded.

• •

OECD. The dataset contains percentages of the labour force and the per capita income of 20 European OECD countries in 1960. The data can be downloaded from the DASL library.

• •

GRTA. This dataset contains information on the number of Greek road traffic accidents and persons injured from January 2010 up to August 2017 on a monthly basis. Three quantities of interest are measured, Killed, seriously injured and slightly injured.

For each of these three datasets, we apply the $\alpha$-transformation (3) and maximise the Dirichlet log-likelihood. This is the same as maximizing the profile log-likelihood of $\alpha$ presented in (6). For a range of values of $\alpha$ we apply the $\alpha$-transformation (3) to the compositional data and maximize the sum of the Dirichlet log-likelihood in (6) and the logarithm of the Jacobian in (5) of the $\alpha$-transformation given in (3). Figure 2 presents these profile log-likelihoods. Tables 1–4 show the resulting estimates of the Dirichlet parameters obtained from the MLE applied to the $\alpha$-transformed data and based on the two asymptotic forms. Asymptotic 1 corresponds to $\textrm{Dirichlet}\left(b({\bf 1}_{D}+\alpha{\bf c})/\alpha^{2}\right)$ model with $b$ and ${\bf c}$ being estimated from (15) and (14), respectively while Asymptotic 2 corresponds to the general case $\textrm{Dirichlet}\left({\bf b}/\alpha^{2}\right)$. The estimation of the parameter vector ${\bf b}$ is based again on MLE. Therefore the results for Asymptotic 2 are expected to be similar with the results when direct MLE is applied.

6 Proofs

Proof of Proposition 1. Let us focus on the log-likelihood for a single observation $i$. At the end of the proof we will sum over $i$ to obtain the log-likelihood (10).

After re-organization, the log-likelihood for observation $i$ may be written

[TABLE]

We now make use of Lemma 1 and Lemma 2 which are stated and proved below. Noting that

[TABLE]

and substituting (22) and (25) into (20), we obtain, as $\alpha\rightarrow 0$,

[TABLE]

where $\mathcal{C}_{0}$ and $\mathcal{C}_{1,i}$ are defined in (11). Finally, sum (21) over $i$ to obtain (10). $\square$

Lemma 1 and Lemma 2 are now stated and proved.

Lemma 1. With the $b_{j}$ defined as in (8) with the $c_{j}$ subject to (9), the following result holds: as $\alpha\rightarrow 0$,

[TABLE]

Proof of Lemma 1. This involves no more than repeated application of Stirling’s formula

[TABLE]

Making use of (9) and applying Stirling’s formula gives

[TABLE]

Similarly,

[TABLE]

Summing (24) over $j=1,\ldots,D$ and subtracting from (23) yields (22) after some further calculations. $\square$

Lemma 2. As $\alpha\rightarrow 0$,

[TABLE]

where $\hat{\kappa}_{1,i}$, $\hat{\kappa}_{2,i}$ and $\hat{\kappa}_{3,i}$ are defined in (12) and (13).

Proof of Lemma 2. As $\alpha\rightarrow 0$,

[TABLE]

where $\hat{\kappa}_{i,1}$, $\hat{\kappa}_{i,2}$ and $\hat{\kappa}_{i,3}$ are defined in (12) and (13). Since, by (9),

[TABLE]

(25) follows easily when we multiply (26) by $(bD)/\alpha^{2}$ and collect terms. $\square$

7 Discussion and Conclusions

Our numerical and theoretical results show that some care is needed when applying a finite-dimensional family of transformations to a parametric model, and then fitting the resulting larger model by maximum likelihood. The particular focus of study here is on using the $\alpha$-transformed Dirichlet distribution as the parametric model. The numerical results in Section 5 and the theoretical results Proposition 1 and Corollary 1 in Section 4 suggest that when the maximum likelihood estimator of $\alpha$ is close to [math], there is a tendency for the other parameter estimates to be rather large. In some of these cases with $\hat{\alpha}$ small one may be better off using the log-ratio transformation, provided there are no sample data components which are zero or very close to zero.

Acknowledgements

This work was partially supported by EPSRC grant EP/K022547/1, for which we are grateful. Partial results from this research were obtained when the second author was a PhD student at the University of Nottingham.