Regression models for compositional data: General log-contrast formulations, proximal optimization, and microbiome data applications

TL;DR

This paper introduces a versatile convex optimization framework for linear log-contrast regression models tailored for compositional data, with applications in microbiome research, offering exact solutions and convergence guarantees.

Contribution

It presents a general convex optimization model for log-contrast regression that encompasses previous methods, along with a proximal algorithm with proven convergence.

Findings

The proposed method performs well on soil microbiome data.

The approach effectively models gut microbiome compositions.

The algorithm guarantees exact solutions with rigorous convergence.

Abstract

Compositional data sets are ubiquitous in science, including geology, ecology, and microbiology. In microbiome research, compositional data primarily arise from high-throughput sequence-based profiling experiments. These data comprise microbial compositions in their natural habitat and are often paired with covariate measurements that characterize physicochemical habitat properties or the physiology of the host. Inferring parsimonious statistical associations between microbial compositions and habitat- or host-specific covariate data is an important step in exploratory data analysis. A standard statistical model linking compositional covariates to continuous outcomes is the linear log-contrast model. This model describes the response as a linear combination of log-ratios of the original compositions and has been extended to the high-dimensional setting via regularization. In this contribution, we propose a general convex optimization model for linear log-contrast regression which includes many previous proposals as special cases. We introduce a proximal algorithm that solves the resulting constrained optimization problem exactly with rigorous convergence guarantees. We illustrate the versatility of our approach by investigating the performance of several model instances on soil and gut microbiome data analysis tasks.

Equations

$$y=X\overline{b}+\overline{o}+Se,\text{with}{C}^{\top }\overline{b}=0,$$

$${\eta }_i={\beta }_1{\xi }_{i,1}-{\beta }_1{\xi }_{i,2}.$$

$${\eta }_i={\beta }_1\log {\zeta }_{i,1}-{\beta }_1\log {\zeta }_{i,2}={\beta }_1\log \frac{{\zeta }_{i,1}}{{\zeta }_{i,2}},$$

$$\underset{\begin{array}{c}\begin{array}{c}b\in {\mathbb{R}}^p\\ {\sum }_{k=1}^p{\beta }_k=0\end{array}\end{array}}{\mathrm{minimize}}\frac{1}{2n}\Vert Xb-y{\Vert }_2^2+\lambda \Vert b{\Vert }_1,$$

$$\underset{\begin{array}{c}\begin{array}{c}b\in {\mathbb{R}}^p\\ C^{\top }b=0\end{array}\end{array}}{\mathrm{minimize}}\frac{1}{2n}\Vert Xb-y{\Vert }_2^2+\lambda \Vert b{\Vert }_1,$$

$$\underset{\begin{array}{c}\begin{array}{c}b\in {\mathbb{R}}^p\\ {\sum }_{k=1}^p{\beta }_k=0\end{array}\end{array}}{\mathrm{minimize}}\frac{1}{2n}\Vert Xb-y{\Vert }_2^2+\lambda \Vert J_{\alpha }b{\Vert }_1$$

$$\underset{\begin{array}{c}b\in {\mathbb{R}}^p,o\in {\mathbb{R}}^n\end{array}}{\mathrm{minimize}}\frac{1}{2n}\Vert Xb-y-o{\Vert }_2^2+{\lambda }_1\Vert b{\Vert }_1+{\lambda }_2\Vert o{\Vert }_1,\text{where}{C}^{\top }b=0,$$

$$h_{\rho }:\mathbb{R}\to \mathbb{R}:u\mapsto \{\begin{array}{ll}\rho |u|-\frac{{\rho }^2}{2},& \text{if}|u|>\rho ;\\ \frac{|u{|}^2}{2},& \text{if}|u|⩽\rho ,\end{array}$$

$$\underset{\begin{array}{c}\begin{array}{c}b\in {\mathbb{R}}^p\\ C^{\top }b=0\end{array}\end{array}}{\mathrm{minimize}}\frac{1}{2n}\sum _{i=1}^n{h}_{\rho }(x_{i}b-{\eta }_i)+{\lambda }_1\Vert b{\Vert }_1.$$

$$\begin{array}[]{rl}\widetilde{\varphi}_{i}\colon\mathbb{R}\times\mathbb{R}^{n_{i}}&\to\left]-\infty,+\infty\right]\\[5.69054pt] (\sigma_{i},u_{i}){}{}&\mapsto\begin{cases}\sigma_{i}\varphi_{i}(u_{i}/\sigma_{i}),&\text{if}\;\>\sigma_{i}>0;\\ \displaystyle{\sup_{u\in\text{\rm dom}\,\varphi_{i}}}\big{(}\varphi_{i}(u+u_{i})-\varphi_{i}(u)\big{)},&\text{if}\;\>\sigma_{i}=0;\\ {+\infty},&\text{if}\;\;\sigma_{i}<0\end{cases}\end{array}$$

$$X=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_N\end{array}\right]\text{and}y=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_N\end{array}\right].$$

$$E=\big{\{}{b\in\mathbb{R}^{p}}~{}\big{|}~{}{C^{\top}b=0}\big{\}}$$

$$\underset{\begin{subarray}{c}{s\in D,\>b\in E}\end{subarray}}{\text{minimize}}\;\;\displaystyle\sum_{i=1}^{N}\widetilde{\varphi}_{i}\big{(}\sigma_{i},X_{i}b-y_{i}\big{)}+\displaystyle\sum_{i=1}^{M}\psi_{i}\big{(}L_{i}b\big{)}.$$

$$\begin{array}[]{rl}\widetilde{\varphi}\colon\mathbb{R}\times\mathbb{R}^{P}&\to\left]-\infty,+\infty\right]\\[5.69054pt] (\sigma,u){}{}&\mapsto\begin{cases}\dfrac{\sigma}{2}+\dfrac{\|u\|_{2}^{2}}{\sigma},&\text{if}\;\>\sigma>0;\\ 0,&\text{if}\;\>\sigma=0\;\>\text{and}\;\>u=0;\\ {+\infty},&\text{otherwise.}\end{cases}\end{array}$$

$$\gamma t^3+2(2\sigma +3\gamma )t-8\Vert u{\Vert }_2=0,$$

$$\text{\rm prox}_{\gamma\widetilde{\varphi}}(\sigma,u)=\begin{cases}\bigg{(}\sigma+\dfrac{\gamma}{2}\bigg{(}\dfrac{t^{2}}{2}-1\bigg{)},u-\gamma p\bigg{)},&\text{if}\;\;4\gamma\sigma+\|u\|_{2}^{2}>2\gamma^{2};\\ (0,0),&\text{if}\;\;4\gamma\sigma+\|u\|_{2}^{2}\leqslant 2\gamma^{2}.\end{cases}$$

$$\begin{array}[]{rl}\widetilde{\varphi}\colon\mathbb{R}\times\mathbb{R}&\to\left]-\infty,+\infty\right]\\[5.69054pt] (\sigma,u)&\mapsto\begin{cases}\dfrac{(1-\rho^{2})\sigma}{2}+\rho|u|,&\text{if}\;\>|u|>\sigma\rho\;\text{and}\;\sigma>0;\\[5.69054pt] \dfrac{\sigma}{2}+\dfrac{|u|^{2}}{2\sigma},&\text{if}\;\>|u|\leqslant\sigma\rho\;\text{and}\;\sigma>0;\\ \rho|u|,&\text{if}\;\>\sigma=0;\\ {+\infty},&\text{if}\;\;\sigma<0.\end{cases}\end{array}$$

$$\text{\rm prox}_{\gamma\widetilde{\varphi}}(\sigma,u)=\Bigg{(}0,\bigg{(}1-\dfrac{\gamma\rho}{|u|}\bigg{)}u\Bigg{)}.$$

$$\text{\rm prox}_{\gamma\widetilde{\varphi}}(\sigma,u)=\Bigg{(}\sigma+\dfrac{\gamma}{2}\big{(}\rho^{2}-1\big{)},\bigg{(}1-\dfrac{\gamma\rho}{|u|}\bigg{)}u\Bigg{)}.$$

$$\gamma t^3+(2\sigma +\gamma )t-2|u|=0.$$

$$\text{\rm prox}_{\gamma\widetilde{\varphi}}(\sigma,u)=\begin{cases}\big{(}\sigma+\gamma(t^{2}-1)/2,u-\gamma t\text{\rm sign}(u)\big{)},&\text{if}\;\;2\gamma\sigma+|u|^{2}>\gamma^{2};\\ \big{(}0,0\big{)},&\text{if}\;\;2\gamma\sigma+|u|^{2}\leqslant\gamma^{2}.\end{cases}$$

$$A=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_N\\ L_1\\ ⋮\\ L_M\end{array}\right],Q=A^{\top }(\mathrm{Id}+{AA}^{\top }{)}^{-1},\text{and}W=\mathrm{Id}-C(C^{\top }C{)}^{-1}{C}^{\top },$$

$$\begin{array}[]{rl}\mathsf{g}\colon&\mathbb{R}^{N}\times\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{N}}\times\mathbb{R}^{m_{1}}\times\cdots\times\mathbb{R}^{m_{M}}\to\left]-\infty,+\infty\right]\\ &\big{(}s,u_{1},\ldots,u_{N},v_{1},\ldots,v_{M}\big{)}\mapsto{\displaystyle\sum_{i=1}^{N}\widetilde{\varphi}_{i}(\sigma_{i},u_{i}-y_{i})+\displaystyle\sum_{i=1}^{M}\psi_{i}(v_{i})},\end{array}$$

$$\{\begin{array}{l}{s}_k=({\sigma }_{1,k},\dots ,{\sigma }_{N,k})\in {\mathbb{R}}^N\\ h_{s,k}=({\eta }_{1,k},\dots ,{\eta }_{N,k})\in {\mathbb{R}}^N\\ h_{b,k}=(h_{1,k},\dots ,h_{N,k},h_{N+1,k},\dots ,h_{N+M,k})\\ \text{hfill}\in {\mathbb{R}}^{n_1}\times \cdots \times {\mathbb{R}}^{n_N}\times {\mathbb{R}}^{m_1}\times \cdots \times {\mathbb{R}}^{m_M}\\ z_{b,k}=(z_{1,k},\dots ,z_{N,k},z_{N+1,k},\dots ,z_{N+M,k})\\ \text{hfill}\in {\mathbb{R}}^{n_1}\times \cdots \times {\mathbb{R}}^{n_N}\times {\mathbb{R}}^{m_1}\times \cdots \times {\mathbb{R}}^{m_M}\\ d_{s,k}=({\delta }_{1,k},\dots ,{\delta }_{N,k})\in {\mathbb{R}}^N\\ d_{b,k}=(d_{1,k},\dots ,d_{N,k},d_{N+1,k},\dots ,d_{N+M,k})\\ \in {\mathbb{R}}^{n_1}\times \cdots \times {\mathbb{R}}^{n_N}\times {\mathbb{R}}^{m_1}\times \cdots \times {\mathbb{R}}^{m_M}.\end{array}$$

$$\begin{array}[]{l}\text{for}\;k=0,1,\ldots\\ \left\lfloor\begin{array}[]{l}\mu_{k}\in[\varepsilon,2-\varepsilon]\\ q_{s,k}=x_{s,k}-h_{s,k}\\ q_{b,k}=Ax_{b,k}-{h}_{b,k}\\ s_{k}=x_{s,k}-q_{s,k}/2\\ b_{k}=x_{b,k}-Qq_{b,k}\\ c_{s,k}=\text{\rm proj}_{D}(2s_{k}-x_{s,k})\\ c_{b,k}=W(2b_{k}-x_{b,k})\\ x_{s,k+1}=x_{s,k}+\mu_{k}(c_{s,k}-s_{k})\\ x_{b,k+1}=x_{b,k}+\mu_{k}(c_{b,k}-b_{k})\\ \text{for}\;i=1,\ldots,N\\ \left\lfloor\begin{array}[]{l}z_{i,k}=X_{i}b_{k}\\ (\delta_{i,k},d_{i,k})=(0,y_{i})+\text{\rm prox}_{\gamma\widetilde{\varphi}_{i}}(2\sigma_{i,k}-\eta_{i,k},2z_{i,k}-h_{i,k}-y_{i})\\ \end{array}\right.\\[5.69054pt] \text{for}\;i=1,\ldots,M\\ \left\lfloor\begin{array}[]{l}z_{N+i,k}=L_{i}b_{k}\\ d_{N+i,k}=\text{\rm prox}_{\gamma{\psi}_{i}}(2z_{N+i,k}-h_{N+i,k})\\ \end{array}\right.\\[5.69054pt] h_{s,k+1}=h_{s,k}+\mu_{k}(d_{s,k}-s_{k})\\ h_{b,k+1}=h_{b,k}+\mu_{k}(d_{b,k}-z_{b,k}).\end{array}\right.\\[5.69054pt] \end{array}$$

$$\underset{\begin{array}{c}s\in D,b\in E\\ \Vert s{\Vert }_2+\Vert b{\Vert }_2\to +\infty \end{array}}{\lim }\mathsf{g}(s,Ab)=+\infty$$

$$\big{(}D\times A(E)\big{)}\cap\text{\rm ri}\,\text{\rm dom}\,\mathsf{g}\neq{\varnothing}.$$

$${\lambda }_0=\sqrt{2}{q}_n(r/p),$$

$$\underset{\begin{subarray}{c}{\sigma\in\mathbb{R},\>b\in E}\end{subarray}}{\text{minimize}}\;\;\widetilde{\|\cdot\|_{2}}\big{(}\sigma,Xb-y\big{)}+\lambda\|b\|_{1},$$

$$\underset{\begin{subarray}{c}{s\in D,\>b\in E}\end{subarray}}{\text{minimize}}\;\;\displaystyle\sum_{i=1}^{n}\widetilde{h_{\rho}}\big{(}\sigma_{i},x_{i}b-\eta_{i}\big{)}+\lambda\|b\|_{1},\\ \text{where}\quad D=\big{\{}{(\sigma,\ldots,\sigma)\in\mathbb{R}^{n}}~{}\big{|}~{}{\sigma\in\mathbb{R}}\big{\}}.$$

Full text

1 Introduction

Compositional data sets are ubiquitous in many areas of science, spanning such disparate fields as geology and ecology. In microbiology, compositional data arise from high-throughput sequence-based microbiome profiling techniques, such as targeted amplicon sequencing (TAS) and metagenomic profiling. These methods generate large-scale genomic survey data of microbial community compositions in their natural habitat, ranging from marine ecosystems to host-associated environments. Elaborate bioinformatics processing tools Callahan2015 ; Caporaso2010 ; Edgar2013 ; Lagkouvardos2017 ; Schloss2009 typically summarize TAS-based sequencing reads into sparse compositional counts of operational taxonomic units (OTUs). The quantification of the relative abundances of OTUs in the environment is often accompanied by measurements of other covariates, including physicochemical properties of the underlying habitats, variables related to the health status of the host, or those coming from other high-throughput protocols, such as metabolomics or flow cytometry.

An important step in initial exploratory analysis of such data sets is to infer parsimonious and robust statistical relationships between the microbial compositions and habitat- or host-specific measurements. Standard linear regression modeling cannot be applied in this context because the microbial count data carry only relative or compositional information. One of the most popular approaches to regression modeling with compositional covariates is the log-contrast regression model, originally proposed in Aitchison1984 in the context of experiments with mixtures. The linear log-contrast model expresses the continuous outcome of interest as a linear combination of the log-transformed compositions subject to a zero-sum constraint on the regression vector. This leads to the intuitive interpretation of the response as a linear combination of log-ratios of the original compositions. In a series of papers, the linear log-contrast model has been generalized to the high-dimensional setting via regularization. The sparse linear log-contrast model, introduced in Lin2014 , considers variable selection via $\ell^{1}$ regularization and has been extended (i) to multiple linear constraints for sub-compositional coherence across predefined groups of predictors Shi2016 ; (ii) to sub-composition selection via tree-structured sparsity-inducing penalties Wang2017 ; (iii) to longitudinal data modeling via a constraint group lasso penalty Sun2018 ; and (iv) to outlier detection via a mean shift modeling approach Mishra2019 . A common theme of these statistical approaches to log-contrast modeling is the formulation of the estimators as the solution of a convex optimization problem, and the theoretical analysis of the statistical properties of these estimators under suitable assumptions on the data.

In the present paper, we take a complementary approach and focus on the structure of the optimization problems underlying log-contrast modeling. We propose an general optimization model for linear log-contrast regression which includes previous proposals as special cases and allows for a number of novel formulations that have not yet been explored. A particular feature of our model is the joint estimation of regression vectors and associated scales for log-contrast models, similar to the scaled Lasso approach in high-dimensional linear regression Sun2012 . This is achieved by leveraging recent results on the connection between perspective functions and statistical models Svva18 ; Combettes2018a ; Combettes2018b . We introduce a Douglas-Rachford splitting algorithm that produces an exact solution to the resulting constrained optimization problems with rigorous guarantees on the convergence of the iterates. By contrast, most existing approaches to solve such problems proceed by first approximating it and then employing coordinate descent methods with less demanding convergence guarantees. We illustrate the versatility of our modeling approach by applying novel log-contrast model instances to environmental and gut microbiome data analysis tasks.

2 Linear log-contrast models

We first introduce the statistical log-contrast data formation model under consideration. We then review several prominent estimators for regularized log-contrast regression models.

2.1 Statistical log-contrast data formation model

Let $Z$ be a known $(n\times p)$-dimensional compositional design matrix with rows $(z_{i})_{1\leqslant i\leqslant n}$ in the simplex \big{\{}{(\zeta_{1},\ldots,\zeta_{p})\in\left]0,+\infty\right[^{p}}~{}\big{|}~{}{\sum_{k=1}^{p}\zeta_{k}=1}\big{\}}. In the microbiome context, each row represents a composition of $p$ OTUs or components at a higher taxonomic rank. We apply a log transform $(x_{i})_{1\leqslant i\leqslant n}=(\log{z_{i}})_{1\leqslant i\leqslant n}$ resulting in the design matrix $X\in\mathbb{R}^{n\times p}$. In this context, we introduce the following log-contrast data formation model.

Model 2.1

The vector $y=(\eta_{i})_{1\leqslant i\leqslant n}\in\mathbb{R}^{n}$ of observations is

[TABLE]

where $X\in\mathbb{R}^{n\times p}$ is the aforementioned design matrix with rows $(x_{i})_{1\leqslant i\leqslant n}$, $\overline{b}\in\mathbb{R}^{p}$ is the unknown regression vector (location), $\overline{o}\in\mathbb{R}^{n}$ is the unknown mean shift vector containing outliers, $e\in\mathbb{R}^{n}$ is a vector of realizations of i.i.d. zero mean random variables, $S\in\left[0,+\infty\right[^{n\times n}$ is a diagonal matrix the diagonal of which are the (unknown) standard deviations, and $C\in\mathbb{R}^{p\times K}$ is a matrix expressing $K$ linear constraints on the regression vector.

The linear log-contrast data formation model is similar to the standard (heteroscedastic) linear model with the important difference that there are linear equality constraints on the regression vector. This stems from the fact that the entries in $X\in\mathbb{R}^{n\times p}$ are not independent due to the compositional nature. In the original model Aitchison1984 , the constraint matrix $C\in\mathbb{R}^{p\times K}$ is the $p$-dimensional all-ones vector $\mathbf{1}_{p}$, resulting in a zero-sum constraint on the regression vector. To gain some intuition about the implications of this constraint, consider a two-dimensional example with given estimates $b=(\beta_{1},\beta_{2})$, and denote by $\xi_{i,1}$ and $\xi_{i,2}$ the first and second column entries of $X$. The linear equality constraint enforces $\beta_{2}=-\beta_{1}$, and thus each observation can be expressed as

[TABLE]

Due to the construction of the design matrix as the log transformation of the compositions, this model is equivalent to

[TABLE]

which expresses the response as a linear function of the log-ratios of the original compositional components. This example also shows that the regression coefficients in the log-contrast model bear a different interpretation than in the standard linear model. Combined log-ratio coefficients relate the response to log-fold changes of the corresponding component ratios.

2.2 Statistical estimators for log-contrast models

2.2.1 Sparse log-contrast regression

In the low-dimensional setting, the standard log-contrast model with zero-sum constraints can be estimated by solving a least-squares problem subject to a linear constraint, or alternatively, via standard linear regression applied to isometrically log-ratio transformed compositions Hron2012 . In the high-dimensional setting, we need structural assumptions on the regression vector for consistent estimation. To this end, the sparse log-contrast model was introduced in Lin2014 . It is based on the optimization problem

[TABLE]

where $\|\cdot\|_{1}$ is the $\ell^{1}$ norm and $\lambda\in\left[0,+\infty\right[$ is a tuning parameter that balances model fit and sparsity of the solution. The estimator enjoys several desirable properties, including scale invariance, permutation invariance, and selection invariance. The latter property is intimately related to the principle of sub-compositional coherence Aitchison1986 and means that the estimator is unchanged if one knew in advance the sparsity pattern of the solution and applied the procedure to the sub-compositions formed by the nonzero components. In Lin2014 , model consistency guarantees are derived for the estimator and the underlying optimization problem is approached via penalization. The proposed iterative algorithm alternates between estimating the Lagrange multipliers and solving a convex subproblem with a coordinate descent strategy. Model selection for the regularization parameter $\lambda$ is performed with a generalized information criterion.

2.2.2 Sparse log-contrast regression with side information

In many situations, it is desirable to incorporate side information about the covariates into log-contrast modeling. For instance, for microbial compositions, each component can be associated with taxonomic or phylogenetic information, thus relating the $p$ components through a rooted taxonomic or phylogenetic tree $\mathcal{T}_{p}$. One way to use this hierarchical tree information is to perform log-contrast regression at a higher taxonomic level, effectively reducing the dimensionality of the regression problem. Let $\mathcal{T}_{p}$ be a tree with $1\leqslant i_{h}\leqslant h$ levels and $p$ leaves and assume that, at given level $i_{h}$, the $p$ compositions split into $K$ groups with sizes $(p_{k})_{1\leqslant k\leqslant K}$. Sub-compositional coherence across the groups can be expressed by the linear constraints $C^{\top}b=0$, where $C$ is an orthogonal $(p\times K)$-dimensional binary matrix. The $k^{th}$ column comprises $p_{k}$ ones at the coordinates of the components that belong to the $k^{th}$ group. Sparse log-contrast regression with group-level compositional coherence can thus be achieved by solving the optimization problem

[TABLE]

where $\lambda\in\left[0,+\infty\right[$ is a tuning parameter. In Shi2016 , model consistency guarantees are derived for this estimator as well as a debiasing procedure for the final estimates. This is done by extending results from Javanmard2014 to the log-contrast setting. In Lin2014 , the underlying optimization problem is approached via an augmented Lagrangian approach, while model selection is achieved by scaling a theoretically derived $\lambda_{0}$ with a data-driven heuristic estimate of the standard deviation $\sigma$ Sun2012 , resulting in $\lambda=\lambda_{0}\sigma$.

An alternative way of incorporating tree information has been proposed in Wang2017 . There, the tree structure is encoded in a parameterized matrix $J_{\alpha}\in\mathbb{R}^{{m-1}\times p}$, where $m$ is the number of vertices in the tree. An estimator based on the minimization problem

[TABLE]

is proposed, where $\lambda\in\left[0,+\infty\right[$ is a tuning parameter. The structure of $J_{\alpha}$ promotes tree-guided sparse sub-composition selection and comprises a weighting parameter $\alpha\in[0,1]$. The authors of Wang2017 are unable to solve the optimization in (6) exactly and resort to a heuristic that abandons the linear constraints and solves a generalized Lasso problem instead. The two tuning parameters $\lambda$ and $\alpha$ are selected via an information criterion.

2.2.3 Robust log-contrast regression

The previous estimators assume the response to be outlier-free with respect to the statistical model under consideration. One way to relax this assumption and to guard against outliers in the response is to use a robust data fitting term. In Mishra2019 , the robust log-contrast regression is introduced via mean shift modeling; see, e.g., Antoniadis2007 ; She2011 . One specific instance of this framework considers the estimation problem

[TABLE]

and where nonzero elements in the mean shift vector $o\in\mathbb{R}^{n}$ capture outlier data, and $\lambda_{1}$ and $\lambda_{2}$ are tuning parameters. In Nguyen2013a , the objective function in (7) is approximated in the form of (5) with a single tuning parameter. As shown in Antoniadis2007 for partial linear models and She2011 for outlier detection, an equivalent form of (7) is to use the Huber function Huber1964 as robust data fitting function and the $\ell^{1}$ norm as regularizer. The Huber function is defined as

[TABLE]

where $\rho\in\left]1,+\infty\right[$ is a fixed parameter with default value $\rho=1.345$ that determines the transition from the quadratic to the linear part. The model in (7) can be written as

[TABLE]

After model estimation, each data point in the linear region of the Huber function is considered an outlier. The latter two models thus allow for joint sparse selection of predictors and outliers in a convex framework.

3 Optimization of general log-contrast models

We introduce an optimization model for general log-contrast regression that includes all previous examples as special cases. We assume that the data follow the data formation model outlined in Model 2.1. Our model belongs to the class of perspective M-estimation models Combettes2018b and allows for joint estimation of regression parameters and corresponding scales while preserving the overall convexity of the model. We then present a proximal algorithm that can solve instances of the optimization model with theoretical guarantees on the convergence of the iterates. Finally, we propose two model selection schemes for practical regularization parameter selection that leverage the joint scale estimation capability of our optimization model.

3.1 Convex optimization model

Let us first introduce some notation (see Livre1 ; Rock70 for details). We denote by $\Gamma_{0}(\mathbb{R}^{n})$ the class of lower semicontinuous convex functions $\varphi\colon\mathbb{R}^{n}\to\left]-\infty,+\infty\right]$ such that \text{\rm dom}\,\varphi=\big{\{}{x\in\mathbb{R}^{n}}~{}\big{|}~{}{\varphi(x)<{+\infty}}\big{\}}\neq{\varnothing}. Given $\varphi\in\Gamma_{0}(\mathbb{R}^{n})$ and $x\in\mathbb{R}^{n}$, the unique minimizer of $\varphi+\|x-\cdot\|_{2}^{2}/2$ is denoted by $\text{\rm prox}_{\varphi}x$. Now let $D$ be a convex subset of $\mathbb{R}^{n}$. Then $\iota_{D}$ is the indicator function of $D$ (it takes values [math] on $D$ and ${+\infty}$ on its complement), $\text{\rm ri}\,D$ is the relative interior of $D$ (its interior relative to its affine hull), and, if $D$ is nonempty and closed, $\text{\rm proj}_{D}=\text{\rm prox}_{\iota_{D}}$ is the projection operator onto $D$.

The following general log-contrast optimization model enables the joint estimation of the regression vector $\overline{b}=(\overline{\beta}_{k})_{1\leqslant k\leqslant p}\in\mathbb{R}^{p}$ and of the scale vector $\overline{s}=(\overline{\sigma}_{i})_{1\leqslant i\leqslant N}\in\mathbb{R}^{N}$ in Model 2.1 within a convex optimization setting.

Problem 1

Consider the setting of Model 2.1. Let $N$ and $M$ be strictly positive integers, let $D$ be a vector subspace of $\mathbb{R}^{N}$, let $(n_{i})_{1\leqslant i\leqslant N}$ be strictly positive integers such that $\sum_{i=1}^{N}n_{i}=n$, let $(m_{i})_{1\leqslant i\leqslant M}$ be strictly positive integers, and set $m=\sum_{i=1}^{M}{m_{i}}$. For every $i\in\{1,\ldots,N\}$, let $\varphi_{i}\in\Gamma_{0}(\mathbb{R}^{n_{i}})$, let

[TABLE]

be the perspective of $\varphi_{i}$, let $X_{i}\in\mathbb{R}^{n_{i}\times p}$, and let $y_{i}\in\mathbb{R}^{n_{i}}$ be such that

[TABLE]

Finally, set

[TABLE]

and, for every $i\in\{1,\ldots,M\}$, let $\psi_{i}\in\Gamma_{0}(\mathbb{R}^{m_{i}})$ and $L_{i}\in\mathbb{R}^{m_{i}\times p}$. The objective is to

[TABLE]

Remark 1

Problem 1 comprises four main components which are associated with different aspects of the general log-contrast regression model.

• •

The perspective functions $(\widetilde{\varphi}_{i})_{1\leqslant i\leqslant N}$ play the role of the loss function in statistical estimation and couple the estimation of the scale vector $s$ and the regression vector $b$. Because the functions $({\varphi}_{i})_{1\leqslant i\leqslant N}$ are convex, the overall minimization problem (13) remains a convex one in $(s,b)$.

• •

Problem 1 allows for the partitioning of the design matrix $X$ and response $y$ into $N$ blocks with individual scale parameters $(\sigma_{i})_{1\leqslant i\leqslant N}$. This is beneficial when data from multiple measurement sources are available for the prediction of the response or when heteroscedasticity in the design matrix is expected for different groups of measurements. Introducing multiple scales has also numerical advantages. Indeed, as discussed in Combettes2018b , certain proximity operators of perspective functions are easier to compute in separable form.

• •

The vector subspaces $D$ and $E$ (see (12)) enforce linear constraints on the scale vector $s=(\sigma_{i})_{1\leqslant i\leqslant N}$ and the regression vector $b$, respectively.

• •

Additional properties of the regression vector, such as (structured) sparsity, are promoted through the use of the penalization functions $(\psi_{i})_{1\leqslant i\leqslant M}$ and the matrices $(L_{i})_{1\leqslant i\leqslant M}$. The penalization functions typically contain a free parameter $\lambda$ the setting of which requires a model selection strategy.

Perspective functions are discussed in Livre1 ; Svva18 ; Combettes2018a ; Combettes2018b ; Rock70 . The construction (10) guarantees that $(\forall i\in\{1,\ldots,N\})$ $\widetilde{\varphi}_{i}\in\Gamma_{0}(\mathbb{R}^{n_{i}})$. We provide below two examples of perspective functions that will be used in the numerical investigations of Section 4.

Example 1

Consider the function $\varphi=\|\cdot\|_{2}^{2}+1/2$ defined on the standard Euclidean space $\mathbb{R}^{P}$. Then (10) yields (see Fig. 1)

[TABLE]

Now fix $(\sigma,u)\in\mathbb{R}\times\mathbb{R}^{P}$ and $\gamma\in\left]0,+\infty\right[$. If $4\gamma\sigma+\|u\|_{2}^{2}>2\gamma^{2}$, let $t$ be the unique solution in $\left]0,+\infty\right[$ to the equation

[TABLE]

and set $p=tu/\|u\|_{2}$ if $u\neq 0$, and $p=0$ if $u=0$. Then (Combettes2018b, , Example 2.4) yields

[TABLE]

A prominent estimator where the perspective function (14) is used as a loss function in conjunction with the $\ell^{1}$ norm as penalization function is the scaled Lasso estimator for high-dimensional sparse linear regression Sun2012 .

Example 2

Set $\varphi=h_{1}+1/2$, where $h_{1}$ is the Huber function of (8). Then (10) yields (see Fig. 2)

[TABLE]

Now fix $(\sigma,u)\in\mathbb{R}\times\mathbb{R}$ and $\gamma\in\left]0,+\infty\right[$. Then (Combettes2018b, , Example 2.5) asserts that $\text{\rm prox}_{\gamma\widetilde{\varphi}}(\sigma,u)$ is computed as follows.

1. (i)

Suppose that $|u|\leqslant\gamma\rho$ and $|u|^{2}\leqslant\gamma(\gamma-2\sigma)$. Then $\text{\rm prox}_{\gamma\widetilde{\varphi}}(\sigma,u)=(0,0)$. 2. (ii)

Suppose that $\sigma\leqslant\gamma(1-\rho^{2})/2$ and $|u|>\gamma\rho$. Then

[TABLE] 3. (iii)

Suppose that $\sigma>\gamma(1-\rho^{2})/2$ and $|u|\geqslant\rho\sigma+\gamma\rho(1+\rho^{2})/2$. Then

[TABLE] 4. (iv)

Suppose that $|u|^{2}>\gamma(\gamma-2\sigma)$ and $|u|<\rho\sigma+\gamma\rho(1+\rho^{2})/2$. If $u\neq 0$, let $t$ be the unique solution in $\left]0,+\infty\right[$ to the equation

[TABLE]

Then

[TABLE]

Using the perspective function (17) as a loss function and the $\ell^{1}$ norm as a penalization function recovers a robust version of the scaled Lasso approach Combettes2018b ; Owen2007 .

3.2 Algorithm

Our algorithmic solution method to solve Problem 1 is patterned after that of Combettes2018b , which relies on an application of the Douglas-Rachford splitting algorithm in a higher-dimensional space. To present the algorithm and its convergence properties, it is convenient to introduce the matrices

[TABLE]

together with the function

[TABLE]

and to define, for every iteration index $k\in\mathbb{N}$, the vectors

[TABLE]

Algorithm 3.1

Let $\gamma\in\left]0,+\infty\right[$, $\varepsilon\in\left]0,1\right[$, $x_{s,0}\in\mathbb{R}^{N}$, $x_{b,0}\in\mathbb{R}^{p}$, $h_{s,0}\in\mathbb{R}^{N}$, and $h_{b,0}\in\mathbb{R}^{n+m}$. Iterate

[TABLE]

Proposition 1

Consider the setting of Problem 1. Suppose that

[TABLE]

and that

[TABLE]

Then Problem 1 has at least one solution. Now let $(s_{k})_{k\in\mathbb{N}}$ and $(b_{k})_{k\in\mathbb{N}}$ be sequences generated by Algorithm 3.1. Then $(s_{k})_{k\in\mathbb{N}}$ converges to some $s\in\mathbb{R}^{N}$ and $(b_{k})_{k\in\mathbb{N}}$ converges to some $b\in\mathbb{R}^{p}$ such that $(s,b)$ solves Problem 1.

Proof.

See Appendix A. ∎

In most practical situations, (26) and (27) are typically satisfied. For example the following describes a scenario that will be encountered in Section 4.

Proposition 2

Consider the setting of Problem 1 and suppose that the following additional properties hold:

1. (i)

For every $i\in\{1,\ldots,N\}$, $\varphi_{i}=\theta_{i}+\alpha_{i}$, where $\theta_{i}\colon\mathbb{R}^{n_{i}}\to\left[0,+\infty\right[$ is convex and $\alpha_{i}\in\left]0,+\infty\right[$. 2. (ii)

For every $i\in\{1,\ldots,M\}$, $\psi_{i}\colon\mathbb{R}^{m_{i}}\to\left[0,+\infty\right[$. 3. (iii)

For some $j\in\{1,\ldots,M\}$, $\psi_{j}(L_{j}b)\to{+\infty}$ as $\|b\|_{2}\to{+\infty}$ while $C^{\top}b=0$. 4. (iv)

$D\cap\left]0,+\infty\right[^{N}\neq{\varnothing}$.

Then (26) and (27) are satisfied.

Proof.

See Appendix B. ∎

3.3 Model selection

In the context of log-contrast regression, a number of different model selection strategies have been proposed, including stability selection Lin2014 ; Meinshausen2010 and Generalized Information Criteria Sun2018 . In Shi2016 , a scale-dependent tuning parameter has been derived where the optimal scale has been found via line search. Our joint scale and regression modeling approach makes this line search obsolete, thus yielding a parameter-free model selection scheme. More specifically, we consider two model selection schemes. Firstly, following Shi2016 , we consider

[TABLE]

where $q_{n}(t)=n^{-1/2}\Phi^{-1}(1-t)$, $\Phi^{-1}$ is the quantile function for the standard normal distribution, and $r$ is the solution to the equation $r=q^{4}_{1}(r/p)+2q^{2}_{1}(r/p)$. In practice, this data-independent model selection scheme may lead to inclusion of spurious coefficients. To assess the robustness of the inferred solutions we combine this theoretically derived regularization with stability selection Meinshausen2010 . The original stability selection approach selects, for every subsample, a small set of predictors from the regularization path, e.g., the first $q$ predictors that appear along the path or the $q$ coefficients that are largest in absolute value across the entire path. We here propose to select, for every subsample, the nonzero coefficients present at regularization parameter $\lambda_{0}$. Note that $\lambda_{0}$ is sample-size dependent and hence needs to be adapted to the specific subsample size used in stability selection. As default values, we consider a subsample size of $\lceil n/2\rceil$ and generate 100 subsamples. The key diagnostic in stability selection is the selection frequency profile for each coefficient. To select a stable set of coefficients, a threshold parameter $t_{s}\in[0.6,0.9]$ is recommended Meinshausen2010 , where all the coefficients with selection frequency above $t_{s}$ are included in the final model.

4 Applications to compositional microbiome data

We apply several instances of the general log-contrast model outlined in Problem 1 in the context of microbiome data analysis tasks. We set $M=1$, $m_{1}=m$, $L_{1}=\operatorname{Id}\,$, and employ as a regularization function the $\ell^{1}$ norm $\psi_{1}=\lambda\|\cdot\|_{1}$. We use the functions in Examples 1 and 2 as instances of the perspective loss functions $\widetilde{\varphi}_{1}$. We refer to these instances as Least Squares (LS) and Huber log-contrast model, respectively. Thus, in case of the LS model, (13) becomes

[TABLE]

while in the case of the Huber model it becomes

[TABLE]

Note that the projection of a vector $s\in\mathbb{R}^{n}$ onto $D$, as required in Algorithm 3.1, is given by

[TABLE]

Dependent on the application, we use different zero-sum constraints on $b$ as specified by the matrix $C$. To solve the various instances of Problem 1, we use Algorithm 3.1 and set the parameter $\mu_{k}=1.9$ and $\gamma=1$. We consider that the algorithm has converged when $\|b_{k}-b_{k+1}\|_{2}<10^{-6}$.

4.1 Body mass index prediction from gut microbiome data

We first consider a cross-sectional study that examines the relationship between diet and gut microbiome composition, where additional demographic covariates, including body mass index (BMI) are available, referred to as COMBO data set Wu2011c . After pre-processing and filtering, the data set comprises the log-transformed relative abundances of $p=87$ taxa at the genus level across $n=96$ healthy subjects. Following previous analyses Lin2014 ; Shi2016 , we investigate the relationship between BMI and the microbial compositions in a log-contrast regression framework. We use $C=\mathbf{1}_{p}$ to model the standard zero-sum constraint. In addition to the compositional covariates, two covariate measurements, fat and calorie intake, are also taken into account via joint unpenalized least squares. We investigate the influence of different loss functions, LS and Huber, as well as the sub-compositional constraints on the quality of the estimation, the size of the support set, and the predictive power.

To highlight the ability of the algorithm to jointly estimate regression and scale we solve the two problems over the regularization path with $40$ $\lambda$ values on a log-linear grid in $[0.0069,\ldots,0.6989]$. We also consider the theoretically derived regularization parameter $\lambda_{0}=0.1997$ from (28). Figure 3a and b show the solution path of the regression vector $b$ for the sparse LS log-contrast model and the Huber model, respectively. Figure 3c displays the corresponding joint scale estimates $\sigma$ for the LS and the Huber model. The estimated regression coefficients at $\lambda_{0}$ are highlighted in Figure 3d. Both models agree on a set of six genera, including Clostridium as strongest negative and Acidaminococcus as the strongest positive predictors. This implies that the log-ratio of Acidaminococcus to Clostridium is positively associated with BMI. Other genera include Alistipes, Megamonas, and Coprobacillus with negative coefficients, and Dorea with positive coefficient. In Lin2014 ; Shi2016 , the genera Alistipes, Clostridium, Acidaminococcus, and Allisonella have been identified as key predictors. The solutions of the perspective log-contrast models corroborates these finding for Clostridium and Acidaminococcus, and to a less extent to Alistipes, whereas the genus Allisonella has only a small strictly positive coefficient in both log-contrast models (Figure 3d).

Next, we consider the stability selection scheme introduced in Section 3.3 with default parameters and threshold $t_{s}=0.7$. Figure 4a shows the stability-based frequency profile for the sparse LS and Huber log-contrast model. For both models, only Clostridium and Acidaminococcus are selected. Stability selection thus leads to a simple explanatory log-ratio model formed by the ratio of the relative abundances of Acidaminococcus to Clostridium. However, when considering the final model prediction results, as shown in Figure 4b for the Huber model, this model can only explain normal to overweight participants (BMI 20-30) as $34$ out of 88 participants are considered outliers in the Huber model. The overall refitted $R^{2}$ is $0.19$ under the Huber model but increases to $0.43$ for the 60 inlier participants.

Next, we investigate the influence of sub-compositional constraints on the stability selection frequency for the two estimation procedures. We follow the analysis of Shi2016 and consider a subset of 45 genera that have the highest relative abundances in the data set. These 45 genera belong to $K=4$ distinct phyla: Actinobacteria (two genera), Bacteroides (eight genera), Firmicutes (32 genera), and Proteobacteria (three genera). The constraint matrix $C$ is hence an orthogonal $(45\times 4)$-dimensional binary matrix. Figure 5a and b show stability selection profile for both the LS and the Huber model with and without compositional constraints, respectively. Figure 5c shows the difference in the selection frequency profiles. Although several genera, including Collinsella, Paraprevotella, Parabacteroides, Faecalibacterium, Oscillibacter, and Parasutterela display significant frequency differences, the two genera Clostridium and Acidaminococcus, both belonging to the Firmicutes phylum, demonstrate again the highest stability both with and without sub-compositional constraints.

4.2 Relationship between soil microbiome and pH

concentration

We next consider a dataset put forward in Lauber2009 comprising $n=88$ soil samples from North and South America. Both amplicon sequencing data and environmental covariates, including pH concentrations, are available and have been re-analyzed via a log-ratio approach in Morton2017a . The amplicon data contains $p=116$ OTUs, and we consider $C=\mathbf{1}_{p}$. We perform stability selection with default parameters as outlined in Section 3.3. The selection frequency of the different regression coefficients is shown Figure 6a. At stability threshold $t_{\text{s}}=0.7$, seven taxa are selected in both models, four Acidobacteria, two Proteobacteria, and one Actinobacteria. After re-estimation of the two perspective log-contrast models on the selected subset, the two taxa of order Ellin6513 and one taxon of family Koribacteraceae have negative coefficients whereas two taxa belonging to the genus Balneimonas as well as one Rubrobacter taxon and one taxon of order RB41 have positive associations (Figure 6b). This allows the model to be compactly represented with only four log-ratios of compositions. The coefficients show small variations with respect to the choice of the data fitting term. The Huber model ($R^{2}=0.88$) deems 29 data points to be outliers in the final estimate (Figure 6c). For completeness, we include the mean absolute deviation (MAD) between the model estimates and the data in Figure 6d. The selected taxa cover a wide range of average pH levels (as provided in Morton2017a ), ranging from $4.9$ to $6.75$, implying that the learned model may indeed generalize to other soil types not present in the current data set.

5 Discussion and conclusion

Finding linear relationships between continuous variables and compositional predictors is a common task in many areas of science. We have proposed a general estimation model for high-dimensional log-contrast regression which includes many previous proposals as special cases Lin2014 ; Mishra2019 ; Shi2016 ; Wang2017 . Our model belongs to the class of perspective M-estimation models which allows for scale estimation in the data fitting term while preserving the overall convexity of the underlying model. This is made possible due to recent advances in the theory of perspective functions Svva18 ; Combettes2018a ; Combettes2018b . We have introduced a Douglas-Rachford algorithm that can solve the corresponding constrained nonsmooth convex optimization problem with rigorous guarantees on the convergence of the iterates. We have illustrated the viability of our approach on two microbiome data analysis tasks: body mass index (BMI) prediction from gut microbiome data in the COMBO study and pH prediction from soil microbial abundance data. Our joint regression and scale estimation enabled the use of a universal single tuning parameter $\lambda_{0}$ Shi2016 to control the sparsity of the estimates. We have combined this approach with stability-based model selection Meinshausen2010 to determine sparse stable sets of log-contrast predictors. For the gut microbiome BMI analysis, the robust Huber log-contrast model identified two genera that can predict BMI well for normal to overweight participants while simultaneously identifying outliers with respect to the log-contrast model. In the soil microbiome data set, we derived a parsimonious pH prediction model which requires only four log-ratios of microbial compositions with an overall $R^{2}=0.88$.

Going forward, we believe that the general log-contrast model and the associated optimization and model selection techniques presented here will provide a valuable off-the-shelf tool for log-contrast regression analysis when compositional data such as microbial relative abundance data are used as predictors in exploratory data analysis.

Appendix A – Proof of Proposition 1

Define $\mathsf{g}$ as in (23) and set

[TABLE]

Then $\mathsf{f}\in\Gamma_{0}(\mathbb{R}^{N+p})$ as the indicator of the vector subspace $D\times E$, and

[TABLE]

where the last identity follows from (Livre1, , Proposition 29.17(iii))). On the other hand, it follows from (Svva18, , Proposition 2.3(ii)) and (Livre1, , Proposition 8.6) that $\mathsf{g}\in\Gamma_{0}(\mathbb{R}^{N+n+m})$. Furthermore, we derive from (Livre1, , Propositions 24.11 and 24.8(ii)) that

[TABLE]

In addition, (27) implies that

[TABLE]

Consequently, $\text{\rm dom}\,(\mathsf{f}+\mathsf{g}\circ\mathsf{L})\neq{\varnothing}$. Thus,

[TABLE]

while, using the variable $\mathsf{w}=(s,b)\in\mathbb{R}^{N+p}$, (26) and (32) imply that

[TABLE]

It therefore follows from (Livre1, , Proposition 11.15(i)) that

[TABLE]

Since (13) is equivalent to

[TABLE]

we infer from (38) that Problem 1 admits at least one solution. Note that (39) can be rewritten as

[TABLE]

Now set $\boldsymbol{u}=(\mathsf{w},\mathsf{z})\in\boldsymbol{\mathcal{H}}=\mathbb{R}^{2N+m+n+p}$ and

[TABLE]

Then $\boldsymbol{F}\in\Gamma_{0}(\boldsymbol{\mathcal{H}})$, $\boldsymbol{G}\in\Gamma_{0}(\boldsymbol{\mathcal{H}})$, and (40) is equivalent to

[TABLE]

Moreover, we deduce from (27) that

[TABLE]

Consequently, using standard relative interior calculus (Rock70, , Section 6), (41) yields

[TABLE]

Therefore, given $\gamma\in\left]0,+\infty\right[$, $\varepsilon\in\left]0,1\right[$, $\boldsymbol{v}_{0}\in\boldsymbol{\mathcal{H}}$, and a sequence $(\mu_{k})_{k\in\mathbb{N}}$ in $[\varepsilon,2-\varepsilon]$, the nonsmooth convex minimization problem (42) can be solved using the Douglas-Rachford algorithm

[TABLE]

which produces a sequence $(\boldsymbol{u}_{k})_{k\in\mathbb{N}}$ that converges to a solution to (42) (Livre1, , Section 28.3). Next, it follows from (Livre1, , Proposition 24.11 and Example 29.19(i)) that

[TABLE]

Now define

[TABLE]

Then we derive from (46) that, given $\mathsf{x}_{0}\in\mathbb{R}^{N+p}$ and $\mathsf{h}_{0}\in\mathbb{R}^{N+n+m}$, (45) becomes

[TABLE]

Let us partition the vectors appearing in (48) according to their scale and regression components as

[TABLE]

In terms of these new variables, using the matrix $Q$ of (22), (32) and (47) yield

[TABLE]

and it follows from (22), (24), (33), (34), and (50) that (48) is precisely (25). Altogether, since $(\boldsymbol{u}_{k})_{k\in\mathbb{N}}=(\mathsf{w}_{k},\mathsf{z}_{k})_{k\in\mathbb{N}}$ converges to a solution to (42), $(\mathsf{w}_{k})_{k\in\mathbb{N}}=(s_{k},b_{k})_{k\in\mathbb{N}}$ converges to a solution to Problem 1.

Appendix B – Proof of Proposition 2

• •

If $s\notin\left[0,+\infty\right[^{N}$, then (10) yields $(\forall b\in\mathbb{R}^{p})$ $\mathsf{g}(s,Ab)={+\infty}$. On the other hand if, for some $i\in\{1,\ldots,N\}$, $\sigma_{i}\in\left]0,+\infty\right[$ then we deduce from (i) that $(\forall b\in\mathbb{R}^{p})$ $\widetilde{\varphi}_{i}(\sigma_{i},X_{i}b-y_{i})=\sigma_{i}\theta_{i}((X_{i}b-y_{i})/\sigma_{i})+\alpha_{i}\sigma_{i}\geqslant\alpha_{i}\sigma_{i}\to{+\infty}$ as $\sigma_{i}\to{+\infty}$. Hence, (ii) entails that $(\forall b\in\mathbb{R}^{p})$ $\mathsf{g}(s,Ab)\to{+\infty}$ as $\|s\|_{2}\to{+\infty}$ while $s\in\left[0,+\infty\right[^{N}$. On the other hand, it follows from (iii) that $(\forall s\in\mathbb{R}^{N})(\forall b\in E)$ $\mathsf{g}(s,Ab)\geqslant\psi_{j}(L_{j}b)\to{+\infty}$ as $\|b\|_{2}\to{+\infty}$. Altogether, (26) holds.

• •

It follows from (i) and (10) that $(\forall i\in\{1,\ldots,N\})$ $\text{\rm ri}\,\text{\rm dom}\,\widetilde{\varphi}_{i}=\left]0,+\infty\right[\times\mathbb{R}^{n_{i}}$. Furthermore, (ii) yields $(\forall i\in\{1,\ldots,M\})$ $\text{\rm ri}\,\text{\rm dom}\,{\psi}_{i}=\mathbb{R}^{m_{i}}$. Therefore $\text{\rm ri}\,\text{\rm dom}\,\mathsf{g}=\left]0,+\infty\right[^{N}\times\mathbb{R}^{n}\times\mathbb{R}^{m}$. Since trivially $A(E)\subset\mathbb{R}^{n+m}$, (27) reduces to (iv).