The Integrated Nested Laplace Approximation for fitting Dirichlet regression models

TL;DR

This paper presents a new Laplace approximation method for Bayesian inference in Dirichlet regression models, enabling analysis of compositional data without data transformation, and introduces an R package extension for practical implementation.

Contribution

The paper develops a novel Laplace approximation approach for Dirichlet regression and provides an R package extension, enhancing Bayesian analysis of compositional data.

Findings

The method performs well in simulation studies.

It effectively models skewed and heteroscedastic compositional data.

The approach is demonstrated on real-world data.

Abstract

This paper introduces a Laplace approximation to Bayesian inference in Dirichlet regression models, which can be used to analyze a set of variables on a simplex exhibiting skewness and heteroscedasticity, without having to transform the data. These data, which mainly consist of proportions or percentages of disjoint categories, are widely known as compositional data and are common in areas such as ecology, geology, and psychology. We provide both the theoretical foundations and a description of how Laplace approximation can be implemented in the case of Dirichlet regression. The paper also introduces the package dirinla in the R-language that extends the R-INLA package, which can not deal directly with Dirichlet likelihoods. Simulation studies are presented to validate the good behaviour of the proposed method, while a real data case-study is used to show how this approach can be applied.

Tables

Table 1. Table 1: Output obtained when summary command is employed in a dirinlaregmodel object.

		mean	sd	0.025quant	0.5quant	0.975quant	mode
y1	intercept	-1.5542	0.2014	-1.9497	-1.5542	-1.1591	-1.5542
	v1	1.0181	0.3690	0.2936	1.0181	1.7420	1.0181
y2	intercept	-1.7601	0.2329	-2.2174	-1.7601	-1.3032	-1.7601
	v2	1.9847	0.3971	1.2051	1.9847	2.7636	1.9847
y1	intercept	-0.0593	0.2500	-0.5501	-0.0593	0.4312	-0.0593
	v3	-1.9470	0.4158	-2.7633	-1.9470	-1.1314	-1.9470
DIC = 1821.0247, WAIC = 1827.6268, LCPO = 913.8319
Number of observations: 100
Number of Categories: 3

Table 2. Table 2: Comparing accuracy between dirinla and R-JAGS computing r ​ a ​ t ​ i ​ o 1 𝑟 𝑎 𝑡 𝑖 subscript 𝑜 1 ratio_{1} and r ​ a ​ t ​ i ​ o 2 𝑟 𝑎 𝑡 𝑖 subscript 𝑜 2 ratio_{2} for β 0 ​ c subscript 𝛽 0 𝑐 \beta_{0c} , β 1 ​ c subscript 𝛽 1 𝑐 \beta_{1c} and β 2 ​ c subscript 𝛽 2 𝑐 \beta_{2c} , c = 1 , … , 4 𝑐 1 … 4 c=1,\ldots,4 .

	$rati{o}_1$		$rati{o}_2$
	dirinla	R-JAGS	dirinla	R-JAGS
${\beta }_{01}$	0.2255	-0.0020	1.0594	1.0026
${\beta }_{02}$	0.2279	0.0237	1.0344	1.0013
${\beta }_{03}$	0.2274	-0.0071	0.9957	1.0042
${\beta }_{04}$	0.2197	0.0075	1.0575	1.0003
${\beta }_{11}$	-0.0056	-0.0219	1.0886	0.9967
${\beta }_{12}$	-0.0116	-0.0180	1.0322	0.9977
${\beta }_{13}$	-0.0108	-0.0099	0.9549	1.0023
${\beta }_{14}$	0.0011	-0.0097	1.0586	1.0011
${\beta }_{21}$	0.0509	-0.0115	0.9318	0.9945
${\beta }_{22}$	0.0560	-0.0006	0.9499	1.0097
${\beta }_{23}$	0.0836	-0.0137	0.9996	0.9890
${\beta }_{24}$	0.0612	-0.0042	1.1022	1.0110

Equations

$$p(\mathbf{y}\mid \mathbf{\alpha })=\frac{1}{\text{B}(\mathbf{\alpha })}\prod _{c=1}^C{y}_c^{{\alpha }_c-1},$$

$${\mathbf{y}}^{\ast }=\frac{\mathbf{y}(N-1)+1/C}{N},$$

$$g({\alpha }_{cn})={\eta }_{cn}^{\ast }={\mathbf{V}}_{n\bullet }^{(c)}{\mathbf{\beta }}_{\bullet c}+{\omega }_{cn},$$

$${\mu }_{cn}=\frac{\exp ({\eta }_{cn}^{\ast \mathbf{\mu }})}{{\sum }_{c=1}^C\exp ({\eta }_{cn}^{\ast \mathbf{\mu }})}=\frac{\exp ({\mathbf{V}}_{n\bullet }^{(c)}{\mathbf{\gamma }}_{\bullet c}+{\omega }_{cn}^{\mathbf{\mu }})}{{\sum }_{c=1}^C\exp ({\mathbf{V}}_{n\bullet }^{(c)}{\mathbf{\gamma }}_{\bullet c}+{\omega }_{cn}^{\mathbf{\mu }})}.$$

$${\tilde{\mathbf{\gamma }}}_{\bullet c}={\tilde{\mathbf{\beta }}}_{\bullet c}-{\tilde{\mathbf{\beta }}}_{\bullet r},$$

$${\omega }_{cn}^{\mathbf{\mu }}={\omega }_{cn}-{\omega }_{rn},$$

$${\bm{\tilde{\eta}}}=\underbrace{\begin{bmatrix}{\bm{\eta}}^{*}_{\bullet 1}\\ \vdots\\ {\bm{\eta}}^{*}_{\bullet N}\end{bmatrix}}_{CN\times 1}\,\$$

$$\tilde{\mathbf{\eta }}=\mathbf{A}\mathbf{x},$$

$$\mathbf{y}\mid \mathbf{x},{\mathbf{\theta }}_1\sim \prod _{n=1}^{N}p(y_n\mid x_n,{\mathbf{\theta }}_1).$$

$$\mathbf{x}\mid {\mathbf{\theta }}_2\sim \mathcal{N}(\mathbf{\mu }({\mathbf{\theta }}_2),{\mathbf{Q}}^{-1}({\mathbf{\theta }}_2))$$

$${\eta }_n={\beta }_0+\sum _j{v}_{nj}{\beta }_j+\sum _k{\omega }_{kn},$$

$$I_n=\int \exp (nf(x))\mathrm{d}x.$$

$$I_n$$

$$=\exp (nf(x_0))\sqrt{\frac{2\pi }{-n{f}^{\prime \prime }(x_0)}}={\tilde{I}}_n.$$

$$p({\gamma }_1)$$

$$\approx {\frac{p(\mathbf{\gamma })}{p_G({\mathbf{\gamma }}_{-1};\mathbf{\mu }({\gamma }_1),\mathbf{Q}({\gamma }_1))}|}_{{\mathbf{\gamma }}_{-1}={\mathbf{\gamma }}_{-1}^{\ast }=\mathbf{\mu }({\gamma }_1)},$$

$$p(x_m\mid \mathbf{y})$$

$$p({\theta }_k\mid \mathbf{y})$$

$${\mathbf{z}}_n^0$$

$$l({\mathbf{y}}_n\mid {\mathbf{\eta }}_n)$$

$${\tilde{\mathbf{\eta }}}^0=\underset{CN\times 1}{\underbrace{\left[\begin{array}{c}{\mathbf{\eta }}_{\bullet 1}^{\ast 0}\\ ⋮\\ {\mathbf{\eta }}_{\bullet N}^{\ast 0}\end{array}\right]}},{\mathbf{g}}_{\tilde{\mathbf{\eta }}}^0=\underset{CN\times 1}{\underbrace{\left[\begin{array}{c}{\mathbf{g}}_1^0\\ ⋮\\ {\mathbf{g}}_N^0\end{array}\right]}},{\mathbf{L}}^0=\underset{CN\times CN}{\underbrace{\left[\begin{array}{ccc}{\mathbf{L}}_1^0& & 0\\ & \ddots & \\ 0& & {\mathbf{L}}_N^0\end{array}\right]}},{\mathbf{H}}_{\tilde{\mathbf{\eta }}}^0=\underset{CN\times CN}{\underbrace{\left[\begin{array}{ccc}{\mathbf{H}}_{{\mathbf{\eta }}_1}^0& & 0\\ & \ddots & \\ 0& & {\mathbf{H}}_{{\mathbf{\eta }}_N}^0\end{array}\right]}},$$

$${\tilde{\mathbf{z}}}^0$$

$${\tilde{\mathbf{z}}}^0\mid \tilde{\mathbf{\eta }}$$

$$p(\mathbf{Y}\mid \mathbf{x},\mathbf{\theta })=p(\mathbf{Y}\mid \tilde{\mathbf{\eta }})\approx p({\tilde{\mathbf{z}}}^0\mid \tilde{\mathbf{\eta }})=p({\tilde{\mathbf{z}}}^0\mid \mathbf{x},\mathbf{\theta }).$$

$$p(\mathbf{\theta },\mathbf{x}\mid \mathbf{y})=p(\mathbf{\theta }\mid \mathbf{y})\cdot p(\mathbf{x}\mid \mathbf{y},\mathbf{\theta }),$$

$$\overline{p}(\mathbf{\theta },\mathbf{x}\mid {\tilde{\mathbf{z}}}^0)=\overline{p}(\mathbf{\theta }\mid {\tilde{\mathbf{z}}}^0)\cdot \overline{p}(\mathbf{x}\mid {\tilde{\mathbf{z}}}^0,\mathbf{\theta }).$$

$${\mathbf{Y}}_{\bullet n}$$

$$\log ({\alpha }_{1n})$$

$$\log ({\alpha }_{2n})$$

$$\log ({\alpha }_{3n})$$

$$rati{o}_1-\text{method}$$

$$rati{o}_2-\text{method}$$

$${\mathbf{Y}}_{\bullet n}$$

$$\log ({\alpha }_{1n})$$

Full text

Keywords: Dirichlet regression, Hierarchical Bayesian models, INLA, multivariate likelihood, random effects

1 Introduction

The use of regression models with multivariate or correlated responses has enormously increased in the last few years. Different R-packages have been developed to deal with them, MCMCglmm (Hadfield, 2010), SabreR (Nowosad and Stepinski, 2018) or mcglm (Bonat, 2018), also packages which uses copulas (Masarotto and Varin, 2017). The use of Multinomial likelihood in the case of multivariate discrete response is one of the most popular (Monyai et al., 2016; Odeyemi et al., 2019; Piccini et al., 2019).

Responses can also be continuous, such as Gaussian (Anderson, 1958) or compositional data (Aitchison and Egozcue, 2005; Hijazi and Jernigan, 2009). Of particular interest are compositional data which consist of proportions or percentages of disjoint categories summing up to one, and play an important role in many fields such as ecology (Kobal et al., 2017; Douma and Weedon, 2019), geology (Buccianti and Grunsky, 2014; Engle and Rowan, 2014), genomics (Tsilimigras and Fodor, 2016; Shi et al., 2016; Washburne et al., 2017) or medicine (Dumuid et al., 2018; Fairclough et al., 2018).

One of the biggest problems one has to face when dealing with models with multivariate or correlated responses is that of performing inference. Their own complexity makes statistical analysis complicated. In the case of compositional data, there are different approaches to deal with these additional complications. One method, due to Aitchison (1986), is based in the idea that “information in compositional vectors is concerned with relative, not absolute magnitudes”, and uses log-ratio analysis to deal with the unit-sum constraint (Aitchison, 1981, 1982, 1983, 1984). Dirichlet regression models (Hijazi and Jernigan, 2009) are another good way of analyzing compositional data. By using appropriate link functions, Dirichlet regression provides a GLM-like framework that relates compositional data with other relevant variables of interest. Beta regression can be considered a special, and effectively univariate, case of the former with only two categories.

Different packages have been implemented in R (R Core Team, 2021) that analyze compositional data using Beta regression and Dirichlet regression, both under the frequentist (Cribari-Neto and Zeileis, 2010; Maier, 2014) and the Bayesian paradigm. In the case of the latter, the largest challenge is implementing the posterior approximation. In particular, it has been implemented in BayesX (Klein et al., 2015), Stan (Sennhenn-Reulen, 2018), BUGS (Van der Merwe, 2019) and R-JAGS (Plummer, 2016). These packages are mainly based on Markov chain Monte Carlo (MCMC) methods, which construct a Markov chain whose stationary distribution converges to the posterior distribution. However, the computational cost of MCMC can be high. On the other hand, the integrated nested Laplace approximation (INLA) methodology (Rue et al., 2009), whose main idea is to approximate the posterior distribution using the Laplace integration method, has become an alternative to MCMC, guaranteeing a higher computational speed for Latent Gaussian models (LGMs).

The INLA methodology is now a well-established tool for Bayesian inference in several research fields, including ecology, epidemiology, econometrics and environmental science, and is implemented in the R-INLA package (Rue et al., 2017). Nevertheless, and spite of its availability for a large number of models, R-INLA does not allow to deal with compositional data when the number of categories is bigger than 2, the reason being that it is constructed for models with univariate responses.

Our objective in this work is twofold. We present an expansion of the INLA method for the particular case of the Dirichlet regression, providing both its theoretical foundations and a description of how it can be implemented for its application, and we introduce the package dirinla in the R-language that allows its practical use. To do so, the remainder of the paper is structured as follows. Section 2 provides the basics of the Dirichlet regression, while Section 3 gives the necessary hints about LGMs and the INLA approach to follow the remainder of the paper. Section 4 depicts the new approach, and Section 5 introduces the dirinla package and how to use it. Simulation studies about the performance of the method introduced is presented in Section 6, followed by an illustration of its use on real data in Section 7. Finally, Section 8 concludes.

2 Dirichlet likelihood

In what follows we present both the Dirichlet distribution and the Dirichlet regression, while introducing some assumptions and the notation that will be used in the rest of the paper.

2.1 Dirichlet distribution

The Dirichlet distribution is the generalization of the widely known Beta distribution, and it is defined by the following probability density,

[TABLE]

where ${\bm{\alpha}}=(\alpha_{1},\ldots,\alpha_{C})$ is known as the vector of shape parameters for each category, $\alpha_{c}>0$ $\forall c$, $y_{c}\in(0,1)$, $\sum_{c=1}^{C}y_{c}=1$, and $\text{B}({\bm{\alpha}})$ is the Multinomial Beta function, which serves as the normalizing constant. The Multinomial Beta function is defined as $\prod_{c=1}^{C}\Gamma(\alpha_{c})/\Gamma(\sum_{c=1}^{C}\alpha_{c})$. The sum of all $\alpha$’s, $\alpha_{0}=\sum_{c=1}^{C}\alpha_{c}$, is usually interpreted as a precision parameter ($\tau$). The Beta distribution is the particular case when $C=2$. In addition, each variable is marginally Beta distributed with $\alpha=\alpha_{c}$ and $\beta=\alpha_{0}-\alpha_{c}$.

Let ${\bm{y}}\sim\mathcal{D}({\bm{\alpha}})$ denote a variable that is Dirichlet distributed. The expected values are $\text{E}(y_{c})=\alpha_{c}/\alpha_{0}$, the variances are $\text{Var}(y_{c})=[\alpha_{c}(\alpha_{0}-\alpha_{c})]/[\alpha_{0}^{2}(\alpha_{0}+1)]$ and the covariances are $\text{Cov}(y_{c},y_{c^{\prime}})=(-\alpha_{c}\alpha_{c^{\prime}})/[\alpha_{0}^{2}(\alpha_{0}+1)]$.

2.1.1 Dealing with zeros and ones

Dirichlet distributions, as it happens with Beta distributions, are defined in the open interval $(0,1)$. Nevertheless, in practice, data may come from the closed interval $[0,1]$. Although converting compositional data into Multinomial data could be an option, there are many situations in real life (see references in the Introduction) where we have to use the continuous nature of the compositional data. In such cases, one option to deal with zeros and ones in Dirichlet distributions is to transform them (Maier, 2014), in a similar way as it happens in Beta distributions (Smithson and Verkuilen, 2006). In particular:

[TABLE]

with $N$ being the number of observations of the $C$ dimensional Dirichlet response. This transformation compresses the data symmetrically around $0.5$, so extreme values are affected more than values lying close to $1/2$. Additionally, as it is pointed out in Maier (2014), if $N\rightarrow\infty$ the compression vanishes, that is, larger data sets are less affected by this transformation. However, for small $N$, the transformation can introduce a noticeable bias. In the INLA software, a new method has been implemented for the Beta distribution that instead treats extreme observations as censored variables, but it is not clear how to extend this to the Dirichlet case.

From now on, we assume that the variable takes values in the open interval $(0,1)$, so that the transformation (2) is not required.

2.2 Dirichlet regression models

Let ${\bm{Y}}$ be a matrix with $C$ rows and $N$ columns denoting $N$ observations for the different categories $C$ of the $C$ dimensional response variable ${\bm{Y}}_{\bullet n}\sim\mathcal{D}({\bm{\alpha}}_{n})$. Let $\eta^{*}_{cn}$ be the linear predictor for the $n$th observation in the $c$th category, so ${\bm{\eta}}^{*}$ is a matrix with $C$ rows and $N$ columns. Let ${\bm{V}}^{(c)}$, $c=1,\ldots,C$, represents a matrix with dimension $N\times J_{c}$ that contains the covariate values for each individual and each category, so ${\bm{V}}^{(c)}_{n\bullet}$ shows the covariate values for the $n$th observation and the $c$th category. Let ${\bm{\beta}}$ be a matrix with $J_{c}$ rows and $C$ columns representing the regression coefficients in each dimension. Finally, let also $\omega_{cn}$ represents a realization of a random effect for the the $c$th category and the $n$th observation. Then, the relationship between the parameters of the Dirichlet distribution and the elements of the linear predictor (including random effects) is set up as:

[TABLE]

where $g(\cdot)$ is the link-function. As $\alpha_{c}>0$ for $c=1,\ldots,C$, log-link $g(\cdot)=\log(\cdot)$ is used. The regression coefficients ${\bm{\beta}}_{\bullet c}$ are a column vector with $J_{c}$ elements.

Dirichlet distributions can also be parametrized in terms of the mean $\mu_{cn}=\frac{\alpha_{cn}}{\sum_{n=1}^{N}\alpha_{cn}}$ and the precision $\tau$. In this case, the relationship between the mean and the linear predictor is stated as:

[TABLE]

The regression coefficients ${\bm{\gamma}}_{\bullet c}$ are now a column vector with $J_{c}$ elements, while $\omega_{cn}^{{\bm{\mu}}}$ represents now a realization of a random effect in this parametrization. Moreover, as one of the categories can be obtained as one less the sum of the rest, we can employ the same Multinomial logit strategy used in Multinomial regression: linear predictor of one category (base category) is set to zero, whereby this category is virtually omitted and becomes the reference.

The equivalence between the two parametrizations comes from the equalities

[TABLE]

and

[TABLE]

where $\tilde{{\bm{\gamma}}}_{\bullet c}$ is a column vector with $J\geq J_{c}$ elements that contains $J_{c}$ non-zero elements and $J-J_{c}$ zero-elements. $J$ dimension comes from a general matrix $\bm{V}$ with dimension $N\times J$ which contains all the covariates values ($J$) for all the observations used in all the categories. $\tilde{{\bm{\beta}}}_{\bullet c}$ is also a column vector with $J\geq J_{c}$ elements, and $\tilde{{\bm{\beta}}}_{\bullet r}$ and $\omega_{rn}$ are the regression coefficients and the realization of a random effect for the base category respectively.

Note that we can include in the formula as many random effects as we consider and from a different nature: temporal, spatial, etc. But here, without loss of generality, we have added just one to show the equivalence between both parametrizations (proof in Appendix A). As the proposal presented in the following Sections relies on the fact that there are no other parameters to estimate in the observation likelihood, from now on we focus on the first parametrization presented.

Equation (3) can be rewritten in a vectorized form. In particular, if

[TABLE]

denotes a restructured linear predictor, being ${\bm{\eta}}^{*}_{\bullet n}$ a column vector representing the linear predictor for the $n$th observation and all the categories, the model in matrix notation is

[TABLE]

where ${\bm{A}}$ is the matrix properly constructed with the covariates values and vectors of 1s for the random effects, and ${\bm{x}}$ a vector formed by the regression coefficients plus the realizations of the random effects. Posteriorly, we will refer to this vector as the Latent Gaussian random field.

3 INLA for Latent Gaussian Models (LGMs)

In this section, we start with a brief explanation about LGMs, the framework where we are going to fit the Bayesian Dirichlet regression (subsection 3.1), followed by the main idea of the Laplace approximation (subsection 3.2) and finishing with the INLA methodology (subsection 3.3).

3.1 LGMs

The popularity of INLA stems from the fact that it allows for fast approximate inference for LGMs, which are a large class of models that include a lot of classically important models (Rue and Held, 2005). LGMs can be written as a three-stage hierarchical model in which observations ${\bm{y}}$ can be assumed to be conditionally independent given a latent Gaussian random field ${\bm{x}}$ and hyperparameters ${\bm{\theta}}_{1}$,

[TABLE]

The versatility of the model class relates to the specification of the latent Gaussian field

[TABLE]

which includes all the latent (non-observable) components of interest, such as fixed effects and random terms, describing the underlying process of the data. The hyperparameters ${\bm{\theta}}=({\bm{\theta}}_{1},{\bm{\theta}}_{2})$ control the latent Gaussian field and/or the likelihood for the data.

The LGMs are a class generalising the large number of related variants of additive and generalized models. If the likelihood $p(y_{n}\mid x_{n},{\bm{\theta}})$ such that “$y_{n}$ only depends on its linear predictor $\eta_{n}$” yields the generalized linear model setup, the set $\{x_{n},n=1,\ldots,N\}$ can be interpreted as $\eta_{n}$, being $\eta_{n}$ the linear predictor which is additive with respect to other effects,

[TABLE]

where $\beta_{0}$ is the intercept, ${\bm{v}}$ represents the fixed covariates with linear effects $\{\beta_{j}\}$, and the terms $\{{\bm{\omega}}_{k}\}$ represent specific Gaussian processes. Each $\omega_{kn}$ is the contribution of the model components ${\bm{\omega}}_{k}$ to the $n$th linear predictor (Rue et al., 2017). If a Gaussian prior is assumed for the intercept and the parameters of the fixed effects, the joint distribution of ${\bm{x}}=\{{\bm{\eta}},\beta_{0},{\bm{\beta}},{\bm{\omega}}_{1},{\bm{\omega}}_{2},\ldots\}$ is a priori Gaussian. This yields the latent field ${\bm{x}}$ in the hierarchical LGM formulation. The hyperparameters ${\bm{\theta}}$ contain the non-Gaussian parameters of the likelihood and the model components. These parameters commonly include variance, scale, or correlation parameters.

In many important cases, the latent field is not only Gaussian, but also sparse Gaussian Markov random field (Rue and Held, 2005, GMRF). A GMRF is a multivariate Gaussian random variable with additional conditional independence properties: $x_{j}$ and $x_{j}^{\prime}$ are conditionally independent given the remaining elements if and only if the $(i,j)$ entry of the precision matrix is [math]. Implementations of the INLA method frequently use this property to speed up computation.

3.2 Laplace Approximation

Laplace approximation (Barndorff-Nielsen and Cox, 1989) is a technique used to approximate integrals of the form

[TABLE]

The main idea is to approximate the target with a scaled Gaussian density that matches the value and the curvature of the target distribution at the mode and evaluate the integral using this Gaussian instead. If $x_{0}$ is the point where $f(x)$ has its maximum, then

[TABLE]

If $nf(x)$ is interpreted as the sum of log-likelihoods and $x$ as the unknown parameter, the Gaussian approximation will be very accurate as $n\rightarrow\infty$ under appropriate regularity conditions.

If we are interested in computing a marginal distribution $p(\gamma_{1})$ from a joint distribution $p({\bm{\gamma}})$, the Laplace approximation of the integral $\int p({\bm{\gamma}})\,\mathrm{d}{\bm{\gamma}}_{-1}$ can be expressed as follows:

[TABLE]

where the first equality holds for any valid ${\bm{\gamma}}^{*}_{-1}$, and the mean ${\bm{\mu}}(\gamma_{1})$ and precision ${\bm{Q}}(\gamma_{1})$ are the parameters of the multivariate Gaussian density derived from the derivatives of $\log p({\bm{\gamma}})$ with respect to ${\bm{\gamma}}_{-1}$, for fixed $\gamma_{1}$. If the posterior is close to a Gaussian density, the results will be more accurate than if the posterior is very non-Gaussian. In this context, unimodality is necessary since the integrand is being approximated with a Gaussian at the mode ${\bm{\gamma}}^{*}_{-1}={\bm{\mu}}(\gamma_{1})$.

3.3 INLA

The main idea of INLA approach is to approximate the posteriors of interest: the marginal posteriors for the latent field, $p(x_{m}\mid{\bm{y}})$, and the marginal posteriors for the hyperparameters, $p(\theta_{k}\mid{\bm{y}})$. These posteriors can be written as

[TABLE]

The nested formulation is used to compute $p(x_{m}\mid{\bm{y}})$ by approximating $p(x_{m}\mid{\bm{\theta}},{\bm{y}})$ and $p({\bm{\theta}}\mid{\bm{y}})$, and then using numerical integration to integrate out ${\bm{\theta}}$. Similarly, $p(\theta_{k}\mid{\bm{y}})$ can be computed by approximating $p({\bm{\theta}}\mid{\bm{y}})$ and integrating out ${\bm{\theta}}_{-k}$.

The marginal posterior distributions in (12) and (13) are computed using the Laplace approximation presented in subsection 3.2. In Rue et al. (2009) it is shown that the nested approach yields a very accurate approximation if applied to LGMs.

All this methodology can be used through R with the R-INLA package. For more details about R-INLA we refer the reader to Blangiardo and Cameletti (2015); Zuur et al. (2017); Wang et al. (2018); Krainski et al. (2018); Moraga (2019); Gómez-Rubio (2020), where practical examples and code guidelines are provided.

However, and despite the advantages of R-INLA implementation, there are some limitations when we deal with multivariate responses. The Multinomial case is solved, as the Multinomial likelihood can be approximated using the Poisson trick. It consists of transforming the Multinomial likelihood into a Poisson likelihood with additional parameters (Baker, 1994). However, there is no method available when we deal with compositional data. In what follows, we propose an expansion of the INLA method for a Dirichlet response variable.

4 Inference in Dirichlet likelihoods

The INLA methodology is a tool that allows us to deal with a wide range of LGMs. However, when a multivariate response is required and several linear predictors are needed to explain it, the implemented R-INLA methodology has some limitations. In the particular case of Dirichlet likelihoods, the main idea to incorporate them in the R-INLA is first to approximate the effect of the log likelihood on the posterior using the Laplace approach and then convert the multivariate response data into observations that R-INLA can deal with. The remainder of the Section presents both the theoretical fundamentals to approximate the effect of the log-likelihood function $\log p({\bm{Y}}\mid{\bm{x}},{\bm{\theta}})$ in the posterior using the Laplace approximation that provides the conditioned independent Gaussian pseudo-observations, and then an algorithmic representation of the method.

4.1 Fundamentals of the approximation

Let ${\bm{\eta}}_{n}:={\bm{\eta}}^{*}_{\bullet n}$ denote the linear predictor corresponding to the $n$th observation ${\bm{y}}_{n}:={\bm{Y}}_{\bullet n}$. If $l({\bm{y}}\mid{\bm{x}})$ represents $-\log p({\bm{y}}\mid{\bm{x}})$ for any ${\bm{y}}$ and ${\bm{x}}$, then $l({\bm{y}}_{n}\mid{\bm{\eta}}_{n})=-\log p({\bm{y}}_{n}\mid{\bm{\eta}}_{n})$ is the log-likelihood function expressed for the $n$th observation, being ${\bm{y}}_{n}$ and ${\bm{\eta}}_{n}$ vectors with $C$ components. Moreover, if ${\bm{\eta}}^{0}_{n}$ is a vector with dimension $C$, we express the gradient of $l({\bm{y}}_{n}\mid{\bm{\eta}}_{n})$ in ${\bm{\eta}}^{0}_{n}$ as ${{\bm{g}}^{0}_{\bm{\eta}}}_{n}=\nabla_{{\bm{\eta}}_{n}}(l)({\bm{\eta}}^{0}_{n},{\bm{y}}_{n})$, and the Hessian of $l({\bm{y}}_{n}\mid{\bm{\eta}}_{n})$ in ${\bm{\eta}}^{0}_{n}$ as ${\bm{H}}^{0}_{{\bm{\eta}}_{n}}$. Depending on which is more computationally convenient, ${\bm{H}}^{0}_{{\bm{\eta}}_{n}}$ can be either the true Hessian ($\nabla^{2}_{{\bm{\eta}}_{n}}(l)({\bm{\eta}}^{0}_{n},{\bm{y}}_{n})$) or the expected Hessian ($\text{E}_{{\bm{y}}_{n}\mid{{\bm{\eta}}_{n}}}(\nabla^{2}_{{\bm{\eta}}_{n}}(l)({\bm{\eta}}^{0}_{n},{\bm{y}}_{n}))$) in ${\bm{\eta}}^{0}_{n}$. Let ${\bm{L}}^{0}_{n}$ be the result of applying the Cholesky factorization to ${\bm{H}}^{0}_{{\bm{\eta}}_{n}}$, ${\bm{H}}^{0}_{{\bm{\eta}}_{n}}={\bm{L}}^{0}_{n}({\bm{L}}^{0}_{n})^{T}$.

Theorem 4.1**.**

If the Laplace approximation method is applied to $l({\bm{y}}_{n}\mid{\bm{\eta}}_{n})$ with respect to ${\bm{\eta}}^{0}_{n}$, then the vector

[TABLE]

is conditionally independent Gaussian distributed

[TABLE]

i.e., ${\bm{z}}^{0}_{n}\mid{\bm{\eta}}_{n}\sim\mathcal{N}(({\bm{L}}^{0}_{n})^{T}{\bm{\eta}}_{n},{\bm{I}}_{d})$ and $z^{0}_{cn}\mid{\bm{\eta}}_{n}\sim\mathcal{N}([({\bm{L}}^{0}_{n})^{T}{\bm{\eta}}_{n}]_{c},1)$, and the constant value of the expression is $l({\bm{y}}_{n}\mid{\bm{\eta}}^{0}_{n})-\frac{1}{2}({\bm{g}}^{0}_{{\bm{\eta}}_{n}})^{T}({\bm{H}}^{0}_{{\bm{\eta}}_{n}})^{-1}{\bm{g}}^{0}_{{\bm{\eta}}_{n}}$.

Proof.

For proof of the theorem see Appendix B. ∎

Theorem 4.1 allows us to convert the observation vector ${\bm{y}}_{n}$ into Gaussian conditionally independent pseudo-observations ${\bm{z}}^{0}_{n}$. More importantly, this theorem can be expanded to multiple observations. In particular, if we denote

[TABLE]

then the following proposition stands.

Proposition 4.2**.**

The matrix

[TABLE]

is conditionally independent Gaussian distributed by columns,

[TABLE]

Proof.

For proof of this proposition see Appendix B. ∎

This approximation has been constructed for a generic ${\bm{\tilde{\eta}}}^{0}$, but, as we are interested in building a Gaussian approximation of the effect of the likelihood on the posterior distribution, ${\bm{\tilde{\eta}}}^{0}$ has been chosen as the posterior mode of $l({\bm{\tilde{\eta}}}\mid{\bm{Y}})$. Then:

[TABLE]

The model posterior is factorized as:

[TABLE]

and the approximation is factorized as:

[TABLE]

Note that the approximation is constructed for a generic ${\bm{\tilde{\eta}}}$, in other words, this approximation is conditionally dependent on the linear predictor. The linear predictor can be formed by fixed effects or random effects. Here, in the implementation of the algorithm, we focus on the case where fixed effects and iid random effects are added to the linear predictor.

4.2 The algorithm

In what follows, we depict the different steps to compute the marginal posterior distributions of the latent field, $p({\bm{x}}_{m}\mid{\bm{Y}})$, and the marginal posterior distribution of the hyperparameters $p({\bm{\theta}}_{k}\mid{\bm{Y}})$. To obtain them, it is necessary to numerically find the mode of the posterior distribution of the linear predictor ${\bm{\tilde{\eta}}}^{0}$. This can be done by means of an iterative method in a similar way to inlabru (Bachl et al., 2019). We define a functional $f(\overline{p}_{{\bm{u}}})$ of the posterior distribution at ${\bm{u}}$ which generates a latent field configuration. Our aim is to find an invariant point of the functional, so that ${\bm{x}}^{0}=f(\overline{p}_{{\bm{x}}^{0}})$. The choice for $f(\cdot)$ is the joint conditional mode $f(\overline{p}_{{\bm{x}}})=$ arg max${}_{{\bm{x}}}\overline{p}_{{\bm{x}}}({\bm{x}}\mid{\bm{\tilde{z}}}^{0},{\bm{\theta}}^{0})$, being ${\bm{\theta}}^{0}=$ arg max${}_{{\bm{\theta}}}\overline{p}_{{\bm{\theta}}}({\bm{\theta}}\mid{\bm{\tilde{z}}}^{0})$.

The final algorithm is:

Let ${\bm{x}}^{1}$ and ${\bm{\theta}}^{1}$ be initial candidate points for the latent variables and the hyperparameters. 2. 2.

Locate a good candidate in the latent variables for the computation of the new pseudo-observations. Compute the mode (${\bm{x}}^{2}$) in ${\bm{x}}$ of $p({\bm{x}}\mid{\bm{Y}},{\bm{\theta}}^{1})$ by means of a quasi-Newton method (Dennis and Moré, 1977) with line search strategy and Armijo conditions (Nocedal and Wright, 2006) in $-\log(p({\bm{x}}\mid{\bm{Y}},{\bm{\theta}}^{1})\propto l({\bm{Y}}\mid{\bm{x}},{\bm{\theta}}^{1})+l({\bm{x}}\mid{\bm{\theta}}^{1})$, being $p({\bm{x}}\mid{\bm{\theta}}^{1})$ multivariate Gaussian as we are in the LGMs context. ${\bm{\tilde{\eta}}}^{1}$ can be easily calculated from the expression ${\bm{\tilde{\eta}}}={\bm{A}}{\bm{x}}$. 3. 3.

Calculate the conditionally independent Gaussian pseudobservations ${\bm{\tilde{z}}}^{1}$. In this case, at the modal configuration established ${\bm{x}}^{1}$, the Hessian matrix ${\bm{H}}^{1}_{{\bm{\tilde{\eta}}}}$ is computed. If the submatrix corresponding to the $n$th individual ${\bm{H}}^{1}_{\tilde{{\bm{\eta}}}_{n}}$ is not positive definite, the expected Hessian is used instead to guarantee a positive definite ${\bm{H}}^{1}_{{\bm{\tilde{\eta}}}}$. Following the approximation previously presented, the Cholesky factorization is computed in ${\bm{H}}^{1}_{{\bm{\tilde{\eta}}}}={\bm{L}}^{1}({\bm{L}}^{1})^{T}$. The gradient (${\bm{g}}^{1}_{{\bm{\tilde{\eta}}}}$) is also calculated in ${\bm{\tilde{\eta}}}^{1}$. According to the equation (16), the scale and rotation of the original observations are done to get the pseudo-observations ${\bm{\tilde{z}}}^{1}$. 4. 4.

Call R-INLA. As pseudo-observations ${\bm{\tilde{z}}}^{1}$ are conditionally independent Gaussian observations, we are able to call R-INLA. If in Step 2 algorithm has converged for a given tolerance, the posterior distributions obtained here are the one that we were looking for: $\overline{p}({\bm{\theta}}\mid{\bm{\tilde{z}}}^{1})$ and $\overline{p}({\bm{x}}\mid{\bm{\tilde{z}}}^{1})$. Else, repeat from step 1 defining ${\bm{x}}^{3}=f(\overline{p}_{{\bm{x}}^{2}})$ as the new starting point for the latent variables, and ${\bm{\theta}}^{2}=$ arg max${}_{{\bm{\theta}}}\overline{p}_{{\bm{\theta}}}({\bm{\theta}}\mid{\bm{\tilde{z}}}^{1})$ the new starting point for the hyperparameters.

After depicting the complete method, we focus on an implementation in R (R Core Team, 2021) of this approximation.

5 The R-package dirinla

In what follows we present dirinla, an R (R Core Team, 2021) package developed to fit Dirichlet regression models. This package can be installed and upgraded via the repository https://github.com/inlabru-org/dirinla. To show how it works, we present an example of a Dirichlet regression without random effect, although, they can be included in the model, as done in Section 6. Then, this Section is divided in three parts: the first one presents the necessary commands to perform a simulation from a Dirichlet regression model; the second one is devoted to show how to fit those models; and the last one depicts how to predict using the package (An extended version of this example is available in the vignette of the R-package dirinla). In particular, we firstly illustrate how to simulate 100 data points from a Dirichlet regression model with three different categories and one different covariate per category:

[TABLE]

being the parameters that compose the latent field $\beta_{01}=-1.5$, $\beta_{02}=-2$, $\beta_{03}=0$ (the intercepts), and $\beta_{11}=1$, $\beta_{12}=2.3$, $\beta_{13}=-1.9$ (the slopes). Note that covariates are different for each category. This could be particularized for a situation where all of them are the same.

For simplicity, covariates are simulated from a Uniform distribution on (0,1). To posteriorly fit the model, and following the structure of LGMs, Gaussian prior distributions are assigned with precision $10^{-4}$ to all the elements of the Gaussian field.

5.1 Data simulation

This subsection presents an example of how simulation can be conducted using the functions of dirinla.

First, we simulate the covariates from a Uniform(0,1):

R> N <- 100R> V <- as.data.frame(matrix(runif((3) * N, 0, 1), ncol = 3))R> names(V) <- paste0(’v’, 1:3)

We then define the formula that we want to fit to keep the values of the different categories in a list. This object will be used to construct the ${\bm{A}}$ matrix. We use the function formula_list() from the package dirinla.

R> formula <- y ~ 1 + v1 | 1 + v2 | 1 + v3R> names_cat <- formula_list(formula)

The values for the parameters composing the latent field are assigned to conduct the simulation. As we have previously depicted, $\beta_{01}=-1.5$, $\beta_{02}=-2$, $\beta_{03}=0$ are the intercepts, and $\beta_{11}=1$, $\beta_{12}=2.3$, $\beta_{13}=-1.9$ are the slopes:

R> x <- c(-1.5, 1, -2, 2.3, 0, -1.9)

We call the function data_stack_dirich() of the package dirinla to construct the ${\bm{A}}$ matrix presented in previous sections. This function uses the inla.stack() structure of the package R-INLA. As a consequence, the returning object is an inla.stack object. Observe that the arguments are the response variable y (in this case it has not been generated yet), the names of the categories covariates, a matrix with the values of the covariates data, the number of categories d and the number of observations N. The sparse matrix ${\bm{A}}$ is then computed.

R> mus <- exp(x) / sum(exp(x))R> C <- length(names_cat)R> A_construct <-+    data_stack_dirich(y = as.vector(rep(NA, N * C)),+                      covariates = names_cat,+                      data       = V,+                      d          = C,+                      n          = N) The next step is to construct the linear predictor as ${\bm{\tilde{\eta}}}={\bm{A}}{\bm{x}}$ using the parameters fixed in the latent field. Using the exponential transformation it is easy to get the parameters ${\bm{\alpha}}$ of the Dirichlet distribution:

R> eta <- A_construct %*% xR> alpha <- exp(eta)

The last stage is to generate the response variable using the function rdirichlet() from DirichletReg (Maier, 2014). The output is a matrix with the response variable summing their rows up to one.

R> y <- rdirichlet(N, alpha)

5.2 Fitting the model

We now present the functions of dirinla needed to fit Dirichlet regression models, the main one being dirinlareg. This function carries out all the steps presented in Section 4, and its use (for the simulated data from previous subsection) is as follows:

R> model.inla <- dirinlareg(+    formula  = y ~ 1 + v1 | 1 + v2 | 1 + v3 ,+    y        = y,+    data.cov = V,+    prec     = 0.0001,+    verbose  = FALSE)

where we have to specify the model formula, the response variable ${\bm{Y}}$ in a matrix format, the data.frame with the covariates data.cov, and the precision of the Gaussian prior (prec) for the latent field ${\bm{x}}$. If we want to follow the process step by step, we can add the instruction verbose = TRUE.

Once the model is fitted, we can summarize the posterior distribution of the fixed effects by means of the function summary applied to the object generated. This object belongs to the dirinlaregmodel class. Three model selection criteria are also displayed: Deviance Information Criterion (Spiegelhalter et al., 2002, DIC), Watanabe-Akaike information criteria (Gelman et al., 2014, WAIC), and the mean of the logarithm of the conditional predictive ordinate (Gneiting and Raftery, 2007, LCPO). Lastly, the number of observations and the number of categories are also depicted (Table 1).

Using the implemented plot method, we obtain the marginal posterior distributions of the latent field for the different categories (Figure 1). These marginals and a summary of each one are stored in model.inla$marginals_fixed and model.inla$summary_fixed (see example.R from the supplementary code to see the details of the code). The plot method also depicts the posterior predictive distribution in the simplex (Figure 2).

Finally, the posterior distribution for the scale parameters of the Dirichlet $\alpha$ can also be computed. In particular, model.inla$marginals_fixed and model.inla$summary_fixed provide the marginals and a summary for each category. Mean parameters can also be obtained via model.inla$marginals_means or model.inla$summary_means, and similarly, model.inla$marginals_precision or model.inla$summary_precision provides the precision parameters.

5.3 Prediction

In most cases, practitioners want to be able to predict the composition of a new observation. The package also provides a function predict to compute posterior predictive distributions for new individuals. To show how it works, we now present how to predict for a value of v1 = 0.2, v2 = 0.5, and v3 = -0.1:

R> model.prediction <-+  predict(model.inla,+          data.pred.cov = data.frame(v1 = 0.2 ,+                                     v2 = 0.5,+                                     v3 = -0.1))

The resulting object also belongs to the dirinlaregmodel class. In a similar way as above, the elements summary_predictive_alphas and marginals_predictive_alphas describe the posterior predictive distribution for the scale parameters of the Dirichlet ${\bm{\alpha}}$, obtained for the new values of the covariates. In addition, means ($\mu$) and precisions ($\tau$) are available via summary_predictive_means, marginals_predictive_means, summary_predictive_precision and marginals_predictive_precision

6 Simulation studies

This section provides a comparison of the performance of the INLA approach for Dirichlet regression models using the dirinla package with the widely used method for Bayesian inference using MCMC algorithms, R-JAGS (Plummer, 2016). The comparison was performed in four different simulated scenarios, and consists of comparisons of computational times and method accuracy.

In all the cases to compare dirinla with R-JAGS, we employed three different methods to make inference: a standard application of the R-JAGS package with a number of iterations enough to achieve a given effective sample size; the INLA methodology through the dirinla package; and a “long” application of the R-JAGS package (for simplicity called long R-JAGS from now on), in this case with a large amount of iterations to get really good representation of the posterior distributions. All computations were performed on a computer with a processor Intel(R) Core(TM) i5-10500 CPU @ 3.10GHz and 32 Gb RAM memory.

To compare the different ways of inference, we implemented four different strategies. The first one was to plot the approximate marginal posterior distributions of each parameter jointly with the real value. The second one consisted on comparing computational times between the three different methods. The remaining two strategies involved the computation of two different ratios:

[TABLE]

where method refers to the method we want to compare with long R-JAGS. It can be either dirinla or R-JAGS. $\phi$ is the parameter of interest. $E(\cdot)$ represents the computed mean of the marginal posterior distribution, and $SD(\cdot)$ is the computed standard deviation of the posterior distribution. In what follows, we make a brief description of the four different scenarios with the aim of extending the explanation at a later stage.

Simulation 1: the first scenario comprised the simulation of data coming from a Dirichlet regression with just intercepts in the model, with different data sizes. 2. 2.

Simulation 2: the second scenario is similar to the previous one, but adding some complexity to the model. In particular, by adding a covariate per category, and checking with different data sizes. 3. 3.

Simulation 3: in the third scenario, we show how the method behaves when random effects are added to the model. The simulations and model include random effects levels. 4. 4.

Simulation 4: the objective of this last scenario is to check how our INLA approach for Dirichlet regression behaves when the number of categories increases. We simulated data from a Dirichlet regression with no slopes but increasing the number of categories.

6.1 Simulation 1

Our first setting is based on a Dirichlet regression with four categories and one parameter per category, the intercept, that is:

[TABLE]

Five different datasets of sizes $N=50,100$, $500,1000,10000$ with this structure were simulated letting $\beta_{0c},c=1,\ldots,4$ to be $-2.4$, $1.2$, $-3.1$ and $1.3$, respectively. We used vague prior distributions for the latent field ($\bm{x}=$ {$\beta_{0c},c=1,\ldots,4$}). In particular, $p(x_{m})\sim$ $\mathcal{N}(0,\tau=0.0001)$, $m=1,\ldots,4$. As the response values are not close to 0 and 1, no transformation was needed.

As above mentioned, for each simulated dataset, we employed three different methods to make inference: a standard application of the R-JAGS package with $2000$ iterations, a burnin of $200$, a thin number of $5$ and $3$ chains; the INLA methodology through the dirinla package; and long R-JAGS, using $1000000$ iterations with a burnin of $100000$, a thin number of $5$ and $3$ chains. To compare methodologies, the strategies presented at the begining of the section were used.

As seen in Figure 3 (the rest of the plots are shown in the Simulation 1 of the supplementary material in https://jmartinez-minaya.github.io/supplementary/CODA/tests.html), the posterior densities have similar shape, with the new method tending to agree more with the long R-JAGS result than the more variable short run R-JAGS does, illustrating how our method reduces estimator variability, at the potential cost of a generally small bias.

With respect the computational effort needed to get those results, Figure 4 displays that, dirinla methodology has a faster computational speed for given data sizes, for the given R-JAGS chain lengths.

In Figure 5, we see that $ratio_{1}$-dirinla is so closed to 0 as happen with $ratio_{1}$-R-JAGS, and $ratio_{2}$-dirinla is always close to 1, similar to $ratio_{2}$-R-JAGS, meaning that using dirinla we obtain as good approximations as with R-JAGS reducing considerably the computational cost.

6.2 Simulation 2

The second setting is based on a Dirichlet regression with a different covariate per category:

[TABLE]

Again, we simulated five different datasets of sizes $N=50,100$, $500,1000,10000$. We set values for $\beta_{0c}$ and $\beta_{1c}$, $c=1,\ldots,4$ to $-1.5,1,-3,1.5,2,-3,-1,5$ respectively, and we simulated covariates from a Uniform distribution with mean in the interval $(0,1)$. We assigned vague prior distributions for the latent field ($\bm{x}=${$\beta_{0c},\beta_{1c},c=1,\ldots,4$}). Particularly, $p(x_{m})\sim$ $\mathcal{N}(0,\tau=0.0001)$, $m=1,\ldots,8$. As the data generated did not present zeros and ones, we did not use any transformation.

As in the previous simulation, we employed the same three different inference methods with the same configurations and the four strategies used in the previous simulation. In Figure 6 (see Simulation 2 in https://jmartinez-minaya.github.io/supplementary/CODA/tests.html for further plots), the posterior distributions are similar in both $\beta_{0c}$ and $\beta_{1c}$, $c=1,\ldots,4$.

Figure 7 depicts that dirinla has a faster computational speed for given data sizes and for given R-JAGS chain lengths. In Figure 8, we show that $ratio_{1}$-dirinla and $ratio_{1}$-R-JAGS are so close to 0 in all cases, and $ratio_{2}$-dirinla and $ratio_{2}$-R-JAGS are always close to 1. Again, we conclude that using dirinla we get similar approximations to R-JAGS, while substantially reducing the computational cost.

6.3 Simulation 3

The third setting is based on a Dirichlet regression with a different covariate per category without intercept and adding two shared independent random effects.

[TABLE]

We simulated four different datasets of sizes $N=50,100$, $500,1000$. We set values for $\beta_{1c}$ for $c=1,\ldots,4$ to $-1.5,2,1,-3$ respectively, and we simulated covariates from a Uniform distribution on the interval $(-1,1)$. Random effects ${\bm{\omega}}_{1}$ and ${\bm{\omega}}_{2}$ were simulated from Gaussian distributions with mean 0 and standard deviations $\sigma_{1}=1/2$ (precision $\tau_{1}=4$) and $\sigma_{2}=1/3$ (precision $\tau_{2}=9$) varying the levels of the factor ($I$), in particular, they were set to $I=2,5,10,25$. The $i_{n}$ sub-index assigns each individual $n$ to a level of the factor.

As we are in the context of Bayesian LGMs, we established Gaussian prior distributions for the latent field, in this case, formed by the parameters corresponding to the fixed effects and the random effects. In particular, we assigned Gaussian prior distributions with mean 0 and precision $0.0001$ to {$\beta_{1c},c=1,\ldots,4$}, and Gaussian priors distribution with mean [math] and precisions $\tau_{1}$ and $\tau_{2}$ for the two shared random effects ${\bm{\omega}}_{1}$ and ${\bm{\omega}}_{2}$. Two types of priors for the $\tau_{1}$ and $\tau_{2}$ parameters were employed. Half-Gaussian priors with location [math] and precision parameter $1$ were used for R-JAGS, long R-JAGS and dirinla. For dirinla, and additional model with a PC-prior(1, 0.01) (Simpson et al., 2017) was also used. The generated data did not contain zeros and ones, so we did not use any transformation. As these models have two hyperparameters, we increased the number of iterations of the R-JAGS method to $20000$ and burnin to $2000$ to achieve a given effective sample size of the MCMC method.

In Figures 9 and 10 (see Simulation 3 in the supplementary material, https://jmartinez-minaya.github.io/supplementary/CODA/tests.html for the complete simulation summaries), we display the posterior distribution of the parameters and the hyperparameters with 2 levels in the factors obtained with the different methods and the different priors. In the case of the precision parameters, we also depict the priors used for them: PC-prior (pc) and Half Gaussian (hn). We observe that the shapes of the posterior densities are similar in all the cases.

Figure 11 shows that dirinla provides higher computational efficiency. And, in view of Figure 12, $ratio_{1}$-dirinla and $ratio_{1}$-R-JAGS are both closed to 0. Similar is the behaviour of $ratio_{2}$, in both cases, it is closed to 1, proving that when we include random effects in the model, we can also get similar approximations by decreasing the computational cost.

If we look at the rest of the simulations conducted (See Simulation 3 in the supplementary material, https://jmartinez-minaya.github.io/supplementary/CODA/tests.html), we observe that large $N$ and small $J$ shows the influence of the different priors. And large $J$ display the difficulties of dirinla in approaching the posterior of the hyperparameters.

6.4 Simulation 4

This last simulation is based on a Dirichlet regression with just one parameter per category and without covariates as in Simulation 1. The idea is that for a fixed data size, we increase the number of categories and then simulate from all those models and compare inferential methods. In particular, we used $N=100$ and $C=5,10,15,20$ and $30$. We also employed R-JAGS, dirinla with the same configuration as in simulations 1 and 2. However, in the case of long R-JAGS, iterations were reduced to 100000 and burnin to 10000.

After comparing computational times needed (Figure 13), we also plot (see Figure 14) the measures $ratio_{1}$ and $ratio_{2}$ computed for the Dirichlet means $\mu_{c}$, $c=1,\ldots,C$ to see how accurate is our method in comparison with R-JAGS. Once more, dirinla shows a higher computational speed than R-JAGS and a similar accuracy.

7 Real example: Arabidopsis thaliana

After validating the use of the package and the approximation in simulated examples, this Section shows an application of the INLA approach for Dirichlet regression models in a real setting. In particular, we worked with a collection of 301 accessions of the annual plant Arabidopsis thaliana in the Iberian Peninsula using the four genetic clusters (gc) infered in Martínez-Minaya et al. (2019): gc1, gc2, gc3 and gc4; categories whose values are summing up to one. Here, we depict how these four gcs can be modeled using Dirichlet regression in terms of two covariates: the annual mean temperature (BIO1) and the annual precipitation (BIO12) (Martínez-Minaya et al., 2019). The aim is explain the distribution of each gc using those two climatic covariates. The complete dataset was downloaded from the repository https://zenodo.org/record/2552025#.YtbLgHZByUl.

The Dirichlet regression model that relates the proportion of each gc, i.e., the multivariate response ${\bm{Y}}_{\bullet n}$, to the covariates is

[TABLE]

Remark that due to the different scale of the covariates, we have scaled both to have mean 0 and standard deviation 1. We assigned vague prior distributions for the latent field, in particular $p(x_{m})\sim$ $\mathcal{N}(0,\tau=0.01)$, $m=1,\ldots,12$. As the data did not include zeros and ones, we did not do any transformation. In a similar approach as in the previous Section, the number of iterations used in R-JAGS were $20000$ with a burnin of $2000$, a thinning of $5$ and $3$ chains, while in the case of the long R-JAGS we used $1000000$ of iterations with a burnin of $100000$, a thinning of $5$ and $3$ chains.

Figure 15 displays the marginal posterior distributions for $\beta_{0c}$, $\beta_{1c}$ and $\beta_{2c}$, $c=1,\ldots,4$. In most cases distributions obtained with R-JAGS perfectly match with those obtained using dirinla. R-JAGS took $405.152$ seconds, dirinla $8.517$ seconds and long R-JAGS $9886.677$.

As in the previous examples, we also computed the measures $ratio_{1}$ and $ratio_{2}$ for dirinla and R-JAGS. Results are depicted in Table 2 and, we can observe that approaches done for dirinla seems to have a good behaviour in comparison with R-JAGS not only in accuracy but also in computational cost.

8 Concluding remarks and future work

In this paper the INLA methodology is extended to fit models with Dirichlet response. In particular, we present both the calculations and a package to make inference and prediction for Dirichlet regression. The main idea underneath the proposed method is to approximate the multivariate likelihoods with univariate ones that can be fitted by R-INLA, in particular, Gaussian likelihoods. This idea is similar to the one proposed for modeling Multinomial likelihood in R-INLA, where using the Poisson trick (Baker, 1994) to reparametrize the model we just need to fit independent Poisson observations. Simpson et al. (2016) use a similar strategy, constructing a Poisson approximation to the true log-Gaussian Cox process likelihood and enabling to make inference on a regular lattice over the observation window by counting the number of points in each cell. This technique has been already implemented in the R package inlabru (Bachl and Lindgren, 2018; Bachl et al., 2019).

It is widely known that there exist a close connection between the Bayesian integrated likelihood and an adjusted profile likelihood (Lee et al., 2018). In particular, the integrated likelihood is approximately the adjusted profile likelihood in the case of orthogonal parameters (Cox and Reid, 1987), we just need to use the Laplace integral approximation to proof it. For the case of Dirichlet response that we deal with in this paper, we propose a quadratic approximation of the likelihood, allowing us to measure the effect of the likelihood in the posterior. However, for the case of the adjusted profile likelihood, we have not found a similar approximation for the Dirichlet regression. Maybe, the method we have developed can help to do approximations using the adjusted profile likelihood of Dirichlet response.

Finally, all examples presented here are focused on models that include fixed and i.i.d. random effects. As we are converting the multivariate observations to conditionally independent Gaussian observations that only depend on the linear predictor, and not directly on the individual latent variables, it should be expected that more complex random effects can be incorporated in the model structures. In particular, the R-INLA method itself uses the same hierarchical model building technique to separate general complex structured random effect components (including spatial, temporal, and spatio-temporal models) from the observation models, that are only linked to the combined linear predictor values. For such more general models, the computational scaling cost benefits in comparison with R-JAGS should be even more pronounced, due to the sparse matrix algebra taken advantage of by R-INLA.

SUPPLEMENTARY MATERIAL

R-package dirinla:

R-package dirinla containing code to perform the Dirichlet regression using R-INLA described in the article. The package also contains the necessary code to reproduce the real and simulation examples presented in the article. It can be installed from https://github.com/inlabru-org/dirinla.

Simulations:

the document test.html (also in https://jmartinez-minaya.github.io/supplementary/CODA/tests.html) contains a complete scenario of simulations by checking the performance of the dirinla R-package.

Appendix

Appendix A Equivalence between parametrizations

In this Appendix, we present the proof for the equivalence between the two parametrizations on Dirichlet regression models presented in the equations (3) and (4).

As mentioned in the text, we just need to construct a matrix ${\bm{V}}$ with dimension $N\times J$ which contains all the observed values for all the observations of all the covariates using in the different categories. Note that if the used covariates are the same in all the categories ${\bm{V}}$ is equal to ${\bm{V}}^{(c)}$. We rewrite equation (3) as:

[TABLE]

being ${\bm{\tilde{\beta}}}_{\bullet c}$ a $J$ column vector having $J_{c}$ non zero values, and being zero the rest of the elements; and equation (4)

[TABLE]

being ${\bm{\tilde{\gamma}}}^{c}$ a $J$ column vector having $J_{c}$ non zero values and being zero the rest of the elements. $\omega_{cn}^{{\bm{\mu}}(c)}$ represents a realization of a random effect for this parametrization.

In the alternative parametrization (28), as we have $C$ categories for which the means must always sum up to $1$, we employ a Multinomial logit strategy as in Multinomial regression, where the linear predictor of one category is set to zero, whereby it is virtually omitted and becomes the reference. Let $\eta_{1n}^{*{\bm{\mu}}}$ be the reference, then:

[TABLE]

Rewriting the equation (29) reveals that

[TABLE]

so that

[TABLE]

and

[TABLE]

as we wanted to proof.

Appendix B Likelihood approximation effect

This Appendix presents both the proof for the Theorem 4.1 and the Proposition 4.2, which is the expansion of the theorem to multiple observations. Let ${\bm{\eta}}_{n}:={\bm{\eta}}_{\bullet n}$ denote the linear predictor corresponding to the $n$th observation ${\bm{y}}_{n}:={\bm{Y}}_{\bullet n}$; let also define $l({\bm{y}}\mid{\bm{x}})=-\log p({\bm{y}}\mid{\bm{x}})$ for any ${\bm{y}}$ and ${\bm{x}}$, and so, let $l({\bm{y}}_{n}\mid{\bm{\eta}}_{n})=-\log p({\bm{y}}_{n}\mid{\bm{\eta}}_{n})$ denote the log-likelihood function expressed for the $n$th observation, being ${\bm{y}}_{n}\in\mathbb{R}^{C}$ and ${\bm{\eta}}_{n}\in\mathbb{R}^{C}$. Using the Taylor series expansion in vector ${\bm{\eta}}^{0}_{n}$, we obtain the approximation:

[TABLE]

where ${\bm{g}}^{0}_{{\bm{\eta}}_{n}}=\nabla_{{\bm{\eta}}_{n}}(l)({\bm{\eta}}^{0}_{n},{\bm{y}}_{n})$ and ${\bm{H}}^{0}_{{\bm{\eta}}_{n}}$ is either the true Hessian ($\nabla^{2}_{{\bm{\eta}}_{n}}(l)({\bm{\eta}}^{0}_{n},{\bm{y}}_{n})$) or the expected Hessian ($\text{E}_{{\bm{y}}_{n}\mid{\bm{\eta_{n}}}}(\nabla^{2}_{{\bm{\eta}}_{n}}(l)({\bm{\eta}}^{0}_{n},{\bm{y}}_{n}))$). $C_{1}$ is a constant whose value is $l({\bm{y}}_{n}\mid{\bm{\eta}}^{0}_{n})-\frac{1}{2}({\bm{g}}^{0}_{{\bm{\eta}}_{n}})^{T}({\bm{H}}^{0}_{{\bm{\eta}}_{n}})^{-1}{\bm{g}}^{0}_{{\bm{\eta}}_{n}}$. Considering now the Cholesky factorization of ${\bm{H}}^{0}_{{\bm{\eta}}_{n}}$, ${\bm{H}}^{0}_{{\bm{\eta}}_{n}}={\bm{L}}^{0}_{n}({\bm{L}}^{0}_{n})^{T}$, expression (B) can be rewritten as follows:

[TABLE]

Defining

[TABLE]

a conditionally Gaussian approximation is constructed.

[TABLE]

Thus, ${\bm{z}}^{0}_{n}\mid{\bm{\eta}}_{n}\sim\mathcal{N}(({\bm{L}}^{0}_{n})^{T}{\bm{\eta}}_{n},{\bm{I}}_{d})$, i.e., $z^{0}_{nc}\mid{\bm{\eta}}_{n}\sim\mathcal{N}([({\bm{L}}^{0}_{n})^{T}{\bm{\eta}}_{n}]_{c},1)$. The observation vector ${\bm{y}}_{n}$ has been converted into conditionally independent Gaussian pseudo-observations ${\bm{z}}^{0}_{n}$. This approximation can be expanded to the $N$ observations.

To present a proof for Proposition (4.2), we rewrite $l({\bm{Y}}\mid{\bm{\eta}})$ for all the observations $N$ as in equation (34).

[TABLE]

Using the notation

[TABLE]

[TABLE]

equation (B) can be rewritten as follows:

[TABLE]

Defining

[TABLE]

we obtain $p({\bm{\tilde{z}}}^{0}\mid{\bm{\tilde{\eta}}})$,

[TABLE]

and the observation matrix ${\bm{Y}}$ has been turned into Gaussian conditionally independent pseudo-observations ${\bm{\tilde{z}}}^{0}$, a likelihood which R-INLA can deal with.

Appendix C Calculus for the Dirichlet likelihood

In this Appendix, we present all the calculus required for the particular case of the Dirichlet likelihood. For the sake of simplicity, we present those required just for one observation. We start with presenting the likelihood in terms of the linear predictor, we continue with the gradient, followed by the Hessian and finishing with the expected Hessian.

C.1 Likelihood

The density function corresponding to Dirichlet distribution has been depicted in expression (1). Using the Dirichlet regression displayed in Equation (3), we know that $\alpha_{c}=\exp(\eta_{c}),\ c=1,\ldots,C$. Then the density function or the likelihood for just one observation can be expressed as:

[TABLE]

Taking logarithms and using the definition of the B function, the next expression is obtained:

[TABLE]

C.2 Gradient

Here, the gradient of the log likelihood is calculated.

[TABLE]

where $c=1\ldots,C$ and $\phi$ is the digamma function.

C.3 Hessian

The second derivatives are calculated for the log likelihood. Let $c$ and $d$ two naturals such as $1\leq c\leq C$ and $1\leq d\leq C$, then

[TABLE]

where $\phi$ is the digamma function and $\phi^{1}$ is the trigamma function.

C.4 Expected Hessian

The expected second derivatives are calculated for the log likelihood. Thus

[TABLE]

where $\phi$ is the digamma function and $\phi^{1}$ is the trigamma function.