# A Copula Type-Model for Examining the Role of Microbiome as a Potential Tool in Diagnosis

Authors: Enrique Calderín–Ojeda; Guillermo López–Campos; Emilio Gómez–Déniz
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Continuous advancements in biotechnology are generating new knowledge and data sources that might be of interest for the insurance industry. A paradigmatic example of these advancements is genetic information which can reliably notify about future appearance of certain diseases making it an element of great interest for insurers. However, this information is considered by regulators in the highest confidentiality level and protected from disclosure. Recent investigations have shown that the microbiome can be correlated with several health conditions. In this paper, we examine the potential use of microbiome information as a potential tool for cardiovascular diagnosis. By using a recent dataset, we analyze the relation of some variables associated to coronary illnesses and several components of the microbiome in the organism by using a new copula-based multivariate regression model for compositional data in the predictor. Our findings show that the coabundance group associated toRuminococcaceae-Bifidobacteriaceae has a negative impact on the age for nonsedentary individuals. However, one should be cautious with this conclusion since environmental conditions also influence the baseline microbiome.


---

## Body

## 1. Introduction

In recent years, the advances in biomedical sciences and biotechnology have enabled an unprecedented leap forward in the amounts and variety of data and information available for research and other purposes. This has translated in new applications and the development of new approaches such as a personalised or precision medicine where the aim is to use these new data and information sources (mostly related with genetic/genomic information) for diagnostic and therapeutic purposes and tailoring them to individuals or groups of patients (see Ginsburg and Phillips [1]). However, genomic information is considered in the highest level of confidentiality and protected from disclosure. For these reasons, the use of genomic data sparked a debate around the ethics and limits associated with the use of this knowledge and information in different sectors how it could be eventually regulated. In an attempt to overcome this potential limitation, in this paper, we propose to explore new avenues and information sources and the use of microbiome as a potential tool for disease diagnosis.The microbiome is defined as the set of microorganisms that live inside or on the organism and its analysis has attracted a great interest in the biomedical domain, particularly since the development of new technologies that have facilitated and reduced the costs of accessing this information. Microbiome is an extremely dynamic element, and it changes with time and environmental conditions and other external factors, such as diet, geographical location, or physical activity and even with interaction between microbes and microbes and the host. The advances in biotechnology allow researchers to measure dynamic behaviors of the microbiota at a large scale (see [2]). Recent studies have shown that differences in the microbiome composition are correlated with an increasing number of conditions ranging from cardiovascular diseases, autoimmune diseases, metabolic diseases, or neurological disorders and mental health aspects [3–7]. Another interesting reference linking a well-established cohort in the biomedical domain (the Framingham Study) with changes in the microbiome for multiple relevant health parameters such as cardiovascular risk, metabolic syndrome, and diabetes is Walker et al. [8]. BMI and physical activity have been also studied in the context of the microbiome finding relationships with different microbiome compositions [9–12]. Quite often, these relationships have been studied under the umbrella of either different ages or other health conditions. For a recent review highlighting that physical activity has an impact in gut microbiota and that physical exercise could be used to control obesity and health (see [13]). Both BMI and physical activity are important factors considered in insurance underwriting. For example, a higher BMI was predominantly related to blood pressure and lipids, which is consistent with results found in the literature (see [14] or [15]). There also exists an association between obesity and higher BMI with all-cause mortality (see [16]). Besides, BMI is connected to increased cancer risk as was recently described by Bhaskaran et al. [17] in a recent paper. On the other hand, microbiome data present a singular challenge due to its inherently high-dimensional and sparse structure. To handle the high dimensionality and compositional nature of the data, Wang et al. [18] proposed a sparse microbial causal mediation model specifically; also, Zhang et al. [19] used an isometric log-ratio transformation of the relative abundances as the mediator variables between treatment and outcome. A statistical approach that enables the inclusion of all daily activity behaviors, based on the principles of compositional data analysis was described by Dumuid et al. [20].In our cross-sectional analysis, using a sample of eligible individuals with unique microorganisms across the Indian microbiome population, due to the large number of operational taxonomic units available in the gut microbiome across the sample, a clustering analysis to reduce the dimensionality of the dataset was initially carried out. Then, the resulting proportions of each of the groups of bacterial coabundance are combined in compositional data predictors that will be used to jointly explain the relationship of age, BMI, and level of physical activity by using a new bivariate regression model for compositional data in the predictor. In this paper, the margins are a beta regression and mixture of logistic regression models. As the age in years of the individuals is restricted to the interval 18–65, a beta regression model indexed by mean and dispersion parameters is considered. This regression family is useful in situations where the dependent variable is continuous and restricted to a bounded interval. On the other hand, regardless of the gender and physical activity level, the empirical distribution of BMI in humans is bimodal. Therefore, choosing suitable parametric models that can capture this feature is crucial; for that reason, a mixture of logistic regression model has been chosen due to its flexibility and simplicity. These margins are linked via at-copula. This is an elliptical copula that is particularly well suited for this purpose as they not only allow for separate modeling of the univariate marginal distributions from the dependency structure but also for covariate adjustment in the margins and uncertainty quantification of their dependence estimates. In addition, they can specify different levels of correlation between the marginals. The compositional data included in the linear predictor are rewritten as logarithms of ratios. Then, we perform estimation via inference for margins method to explain the age, BMI, and level of physical activity using as margins a beta regression and mixture of regression models. Although copula models have been widely applied to model the joint distributions with mixed margins, copula models with the margins proposed in this work with compositional data in the predictor have not been extensively studied in the literature.The rest of the paper is structured as follows: in Section2, an examination of a human microbiome dataset is carried out. Here, a cluster analysis to classify the proportion of the most significant bacterial coabundance groups in the sample is completed. Furthermore, an approach to deal with the implementation of microbiome data as compositional data in the predictor is presented. The relationship of age, body mass index, and physical activity level with microbiome is analyzed in Section 3. Here, we firstly consider the marginal relation of age given a level of physical activity with microbiome through a beta regression model. Next, we examine the connection of BMI with the microbiome given the level of physical activeness by using a mixture of logistics regression model. Later, the joint relationship of these variables is examined by using a t-copula. Finally, discussion and extensions conclude the paper.


## 2. Analysis of a Human Microbiome Dataset

In our analyses we use a dataset available in Dubey et al.’s [21] LogMPIE study. This dataset is freely accessible, and it may be downloaded from the European Nucleotide Archive (ENA) portal of the European Bioinformatics Institute (https://www.ebi.ac.uk/ena/data/view/PRJEB25642). In this study, as it was portrayed in the original description of the dataset, they identify and map the Indian gut microbiome. It was carried out in fourteen geographical locations. Individuals were uniformly selected across geographical regions and some variables associated with changes in the structure of microbiome such as BMI, age in years of the individual, restricted to the interval 18–65 and level of physical activity (sedentary-nonsedentary) and gender (male-female) were also considered in their study design. In addition, a subject is classified as an obese if his/her BMI is greater than 30. This study recorded data from 1004 eligible individuals and reported 993 unique microorganisms across the Indian microbiome population. Unfortunately, in this dataset neither a longitudinal analysis across time of individuals nor changes in the composition of microbiome in old subjects are available.


### 2.1. Cluster Analysis

In general, microbiome empirical distribution includes a high proportion of zero observation and a truncation point mass to account for high values that are too sparse to model; for that reason, models that gives an accurate estimates of the true proportion of zeros have been considered in the literature (see [22] and [23]). In addition, given the dynamic character of the microbiome other techniques such as functional response regression on correlated longitudinal microbiome sequencing data has been recently considered in the literature [24]. In this work, in order to facilitate further analyses and reduce the dimensionality of our data set, we started carrying out a cluster analysis. A main task of exploratory data mining, to group a set of bacterial coabundance collections in such a way that objects in the same group or cluster are more similar to each other than to those in other clusters. We performed our clustering based on the coabundance of genus-like groups at a taxonomic level of species within a sample of 1004 subjects. A total of 993 bacterial genera were identified. The core microbiota analysis was completed by using the Hierarchical Ordered Partitioning and Collapsing Hybrid (HOPACH) package in R that can be downloaded from the Bioconductor website http://www.bioconductor.org/. This package includes the HOPACH clustering algorithm that assembles a hierarchical tree of clusters by recursively portioning the whole dataset while ordering and collapsing clusters at each level. In our analysis, we have discarded redgenus that contain at least a minimum relative abundance of 30%, i.e., 70% of zeros in the sample of 1004 individuals. The algorithm uses the MSS (Mean/Median Split Silhouette) criteria to identify the level of the tree with maximally homogeneous clusters. The correlation distance (cor) was the metric selected for clustering the microbiome species by calculating dissimilarities between variables. We have also used a nonparametric bootstrap to estimate the probability that each species belongs to each cluster and to better understand the variability of each cluster. For that reason, we employed the “boothopach” function by taking 1000 bootstrap resample datasets to obtain a suitable balance between precision and speed. As a result of this, we were able to group the microbiome in five groups containing different numbers of genera (see supplementary tables in Table 1). The five different group of bacteria (classes) identified from the cluster analyses could be associated with different taxonomic groups according to the most abundant or representative genus for each of the identified clusters. Groups 1 and 4 are the two largest groups in terms of number of taxonomic elements. Also, as in Group 1, a majority of members comes from the Bacteroidales-Bacteroidaceae group that represents almost 2/3 of the species contained in this cluster (17 out of 27 members), it could be related to Bacteroidales-Bacteroidaceae cluster. Group 4 is associated with Lachnospiraceae which represent almost 1/3 of the total in this group (7 out of 23 members). The other three groups (2, 3, and 5) were assigned to the Ruminococcaceae-Bifidobacteriaceae group (5 out of 19 members), Negativicutes group (4 out of 19 members), and Pasteurellaceae group (3 out of 15 members), respectively. The results and relationships between the different elements on each of the clusters are presented in Figure 1. Here, species close to each other in the tree are shown in a similar way. The ordered distance matrix shows the clustering structure. Similar clusters appear as blocks on the diagonal of this heatmap. Darker colours represent small distances whereas the lighter colours represent large distances. The identified clusters have different sizes and compositions, with two large coabundance clusters, grouping the majority of the genus analyzed. It is important to note that we have combined under the name Group 0 all the discarded operational taxonomic units, that is, all the species with more that 70% of zeroes in the sample.Table 1
Composition of the five clusters identified in the analysis of the Indian dataset (Dubey et al., [21]).
Group 1Bacteroidetes-Bacteroidia-Bacteroidales-Rikenellaceae-Alistipes-onderdonkiiBacteroidetes-Bacteroidia-Bacteroidales-Rikenellaceae-Alistipes-putredinisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-coprocolaBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-coprophilusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-doreiBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-fragilisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-intestinalisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-plebeiusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-thetaiotaomicronBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-uniformisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-vulgatusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-xylanisolvensBacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Barnesiella-intestinihominisFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-biformeFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-ventriosumFirmicutes-Clostridia-Clostridiales-Clostridiales-Flavonifractor-plautiiProteobacteria-Gammaproteobacteria-Enterobacteriales-Enterobacteriaceae-Klebsiella-variicolaFirmicutes-Bacilli-Lactobacillales-Lactobacillaceae-Lactobacillus-rogosaeFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Megasphaera-sp.Bacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Odoribacter-splanchnicusBacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Parabacteroides-distasonisBacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Parabacteroides-merdaeBacteroidetes-Bacteroidia-Bacteroidales-Prevotellaceae-Prevotella-copriFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-sp.Firmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-bromiiFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-callidusProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Sutterella-wadsworthensisGroup 2Actinobacteria-Actinomycetales-Actinomycineae-Actinomycetaceae-Actinomyces-odontolyticusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-massiliensisActinobacteria-Actinobacteria-Bifidobacteriales-Bifidobacteriaceae-Bifidobacterium-adolescentisActinobacteria-Actinobacteria-Bifidobacteriales-Bifidobacteriaceae-Bifidobacterium-bifidumFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-lutiFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-wexleraeFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Butyrivibrio-crossotusActinobacteria-Coriobacteriia-Coriobacteriales-Coriobacteriaceae-Collinsella-aerofaciensFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-eligensFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Faecalibacterium-prausnitziiFirmicutes-Clostridia-Clostridiales-Clostridiales-Howardella-ureilyticaProteobacteria-Gammaproteobacteria-Enterobacteriales-Enterobacteriaceae-Klebsiella-pneumoniaeFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Megasphaera-micronuciformisProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Parasutterella-excrementihominisFirmicutes-Clostridia-Clostridiales-Peptostreptococcaceae-Peptostreptococcus-stomatisFirmicutes-Negativicutes-Acidaminococcales-Acidaminococcaceae-Phascolarctobacterium-faeciumSpirochaetes-Spirochaetes-Spirochaetales-Spirochaetaceae-Treponema-succinifaciensFirmicutes-Erysipelotrichia-Erysipelotrichales-Erysipelotrichaceae-Turicibacter-sanguinisLentisphaerae-Lentisphaeria-Victivallales-Victivallaceae-Victivallis-vadensisGroup 3Proteobacteria-Alphaproteobacteria-Rhodospirillales-Acetobacteraceae-Acidiphilium-sp.Verrucomicrobia-Verrucomicrobiae-Verrucomicrobiales-Akkermansiaceae-Akkermansia-muciniphilaProteobacteria-Gammaproteobacteria-Aeromonadales-Succinivibrionaceae-Anaerobiospirillum-succiniciproducensActinobacteria-Actinobacteria-Bifidobacteriales-Bifidobacteriaceae-Bifidobacterium-longumProteobacteria-Deltaproteobacteria-Desulfovibrionales-Desulfovibrionaceae-Bilophila-wadsworthiaFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-obeumFirmicutes-Erysipelotrichia-Erysipelotrichales-Erysipelotrichaceae-Bulleidia-p-1630-c5Firmicutes-Clostridia-Clostridiales-Lachnospiraceae-Coprococcus-catusProteobacteria-Deltaproteobacteria-Desulfovibrionales-Desulfovibrionaceae-Desulfovibrio-pigerFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Dialister-succinatiphilusFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-siraeumFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Megasphaera-elsdeniiFirmicutes-Negativicutes-Selenomonadales-Selenomonadaceae-Mitsuokella-jalaludiniiFirmicutes-Negativicutes-Selenomonadales-Selenomonadaceae-Mitsuokella-multacidaBacteroidetes-Bacteroidia-Bacteroidales-Prevotellaceae-Prevotella-stercoreaFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Romboutsia-ilealisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminiclostridium-siraeumActinobacteria-Coriobacteriia-Eggerthellales-Eggerthellaceae-Slackia-isoflavoniconvertensProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Sutterella-sp.Group 4Bacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-caccaeBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-ovatusFirmicutes-Erysipelotrichia-Erysipelotrichales-Erysipelotrichaceae-Catenibacterium-mitsuokaiFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-bartlettiiFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Coprococcus-comesFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Coprococcus-eutactusFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Dorea-formicigeneransFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Dorea-longicatenaFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-hadrumFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-halliiFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-ramulusProteobacteria-Alphaproteobacteria-Rhizobiales-Hyphomicrobiaceae-Gemmiger-formicilisFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Lachnoclostridium-clostridioformeFirmicutes-Clostridia-Clostridiales-Oscillospiraceae-Oscillibacter-sp.Firmicutes-Negativicutes-Acidaminococcales-Acidaminococcaceae-Phascolarctobacterium-succinatutensFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-faecisFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-inulinivoransFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-faecisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-gnavusFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-sp.Firmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-torquesFirmicutes-Negativicutes-Selenomonadales-Selenomonadaceae-Selenomonas-bovisProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Sutterella-stercoricanisGroup 5Proteobacteria-Gammaproteobacteria-Pasteurellales-Pasteurellaceae-Actinobacillus-minorFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-faecisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-productaFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-disporicumFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-perfringensFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-sp.Proteobacteria-Deltaproteobacteria-Desulfovibrionales-Desulfovibrionaceae-Desulfovibrio-D168Proteobacteria-Gammaproteobacteria-Enterobacteriales-Enterobacteriaceae-Escherichia-coliProteobacteria-Gammaproteobacteria-Pasteurellales-Pasteurellaceae-Haemophilus-parainfluenzaeProteobacteria-Gammaproteobacteria-Pasteurellales-Pasteurellaceae-Haemophilus-pittmaniaeFirmicutes-Bacilli-Lactobacillales-Lactobacillaceae-Lactobacillus-ruminisProteobacteria-Gammaproteobacteria-Pseudomonadales-Pseudomonadaceae-Pseudomonas-liniFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-intestinalisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-gauvreauiiFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Veillonella-disparFigure 1
Heatmap of coabundance groups. This figure represents the distance among the abundances of the different genus of bacteria characterised in the sequencing analyses. The dashed lines depict the boundaries between the five different clusters. Darker colours represent closer distances (coabundance) between genus whereas lighter colours displays larger distances.Table2 displays the mean, median, and standard deviation for each one of the coabundance groups. It is noticeable that the proportion of bacteria that belongs to Group 1 is higher in average than the proportion in the other groups. The variability is also larger for the first coabundance group.Table 2
Mean, median, and standard deviation for each coabundance group.
Coabundance groupMeanMedianStandard deviationGroup 10.55480.56690.1542Group 20.18990.16700.1024Group 30.08160.06240.0688Group 40.08990.08130.0517Group 50.06030.02390.0859Group 00.02350.01200.0632In Figure2, some ternary plots for different combinations of the bacterial groups are displayed. In particular we have compared the coabundance Group 1, with Group 2 (top left), Group 3 (top right), Group 4 (bottom left), and Group 5 (bottom right). In order to ensure that the total sum is one, we have combined the coabundance proportion for the rest of the groups in each graph as Others. Group 1 is always located at the top of each triangle. The proportion of coabundance of Group 1 is measured in terms of the horizontal lines, i.e., 0% of coabundance is measured in terms of base of the triangle (farthest from the vertex Group 1). In the lower left apex of each triangle is represented the groups compared to Group 1. The right side of the triangle now becomes the baseline for the percentage of the groups located in this vertex. Finally, the combined groups are located at the lower right apex of the triangle.Figure 2
Ternary plots associated to different coabundance groups. (a) Group1/Group 2/Other groups. (b) Group 1/Group 3/Other groups. (c) Group 1/Group 4/Other groups. (d) Group 1/Group 5/Other groups.
(a)(b)(c)(d)The rate of coabundance for the combined groups is calculated from the left side of the triangle (0% abundance) to the lower right corner (100% abundance). It is observable that the data lie from a high amount of coabundance of Group 1 and Group 2 with a low coabundance of third, fourth, and fifth groups (top left graph). From the rest of the graphs, it can be inferred that when Group 1 is compared to the other groups, the coabundance of these groups is lower than in the former graph. Also, as Group 2 has been included in the lotOthers, the corresponding coabundance of the combined group is higher than in the top left graph.


#### 2.1.1. Compositional Data Predictor

Compositional data can be defined as arrays of strictly positive numbers for which ratios between them are important without any further requirement [25]. Microbiome data are compositional, that is, the distance between component values is only meaningful proportionally (see [26]). The elements of the composition are non-negative and sum to unity. An important issue in microbiome data is the large presence of zeros; however, the issue of zero values in some components is not addressed in most papers and especially in the task of regression. In general, in compositional research problems, most of the basic statistical analysis tools are incorrect unless the variables are rewritten in terms of logarithms of ratios as proposed in the log-ratio methodology for compositional data. After computing these log-ratios, standard regression methods can be used since the relative character of the information is considered when analyzing the results, as one group or variable can only increase in relative terms if some other group or groups reduce. In this work we focus on the case of compositional data being included in the predictor variables. The effect of increasing one of the variables in relative terms in the predictor therefore depends on which other variables are decreased when this occur. In log-ratio parlance, the effect of increasing one log-ratio is interpreted while keeping all other log-ratios constant as the same log-ratio may have different meaning depending on the way that the other log-ratios in the model are assembled. Thus, the interpretation of log-ratios as explanatory variables is usually different from other approaches. Several different approaches of building and interpreting the log-ratios have been considered in the literature, often leading to the same predictions and residuals [27]. Among the different parametrizations, in this work, we have chosen centred log-ratios [28]. In our analysis, we consider a vector of 6-dimensional real space that carries information on the relative importance of its components,(1)xi=xi0,xi1,…,xi5∈ℝ+6,withxij>0,j=0,1,2,…,5,∑j=05xij=1,where i=1,…,1004. Note that the explanatory variables xij represent the proportion of the bacterial coabundance proportion of Group j in individual i. Centred log-ratios are calculated by using a quotient between each variable and the geometric mean of all components (see [28]),(2)log2xj∏j=05xj6,withj=0,1,2,…,5.The fact that we are using logarithms to base 2 means that a unit increase in this logarithm leads to a double increase in the original magnitude. In order to avoid perfect collinearity one centred log-ratio must be deleted from the regression equation. Since all six centred log-ratios add-up to zero, by increasing a fixed centred log-ratio while keeping the other four remaining log-ratios in the regression equation (with regressorsβj with j=0,1,2,…,5) constant implies increasing the given centred log-ratio whilst reducing the omitted centred log-ratio by the same amount. In this regard, a positive statistically significant regression βj coefficient indicates an increasing value of the covariate xj at the expense of decreasing the amount of the omitted component has a significant positive effect on expected value of the response variable. This is equivalent to say that in terms of logarithm of base 2, βj is interpreted as the expected change in the response variable when the ratio between xj and the omitted explanatory variable is multiplied by six. Finally, in order to obtain the estimates and their corresponding p values for all possible pair combinations, the model needs to be repeated six different times by ignoring each time a different centred log-ratio.


## 2.1. Cluster Analysis

In general, microbiome empirical distribution includes a high proportion of zero observation and a truncation point mass to account for high values that are too sparse to model; for that reason, models that gives an accurate estimates of the true proportion of zeros have been considered in the literature (see [22] and [23]). In addition, given the dynamic character of the microbiome other techniques such as functional response regression on correlated longitudinal microbiome sequencing data has been recently considered in the literature [24]. In this work, in order to facilitate further analyses and reduce the dimensionality of our data set, we started carrying out a cluster analysis. A main task of exploratory data mining, to group a set of bacterial coabundance collections in such a way that objects in the same group or cluster are more similar to each other than to those in other clusters. We performed our clustering based on the coabundance of genus-like groups at a taxonomic level of species within a sample of 1004 subjects. A total of 993 bacterial genera were identified. The core microbiota analysis was completed by using the Hierarchical Ordered Partitioning and Collapsing Hybrid (HOPACH) package in R that can be downloaded from the Bioconductor website http://www.bioconductor.org/. This package includes the HOPACH clustering algorithm that assembles a hierarchical tree of clusters by recursively portioning the whole dataset while ordering and collapsing clusters at each level. In our analysis, we have discarded redgenus that contain at least a minimum relative abundance of 30%, i.e., 70% of zeros in the sample of 1004 individuals. The algorithm uses the MSS (Mean/Median Split Silhouette) criteria to identify the level of the tree with maximally homogeneous clusters. The correlation distance (cor) was the metric selected for clustering the microbiome species by calculating dissimilarities between variables. We have also used a nonparametric bootstrap to estimate the probability that each species belongs to each cluster and to better understand the variability of each cluster. For that reason, we employed the “boothopach” function by taking 1000 bootstrap resample datasets to obtain a suitable balance between precision and speed. As a result of this, we were able to group the microbiome in five groups containing different numbers of genera (see supplementary tables in Table 1). The five different group of bacteria (classes) identified from the cluster analyses could be associated with different taxonomic groups according to the most abundant or representative genus for each of the identified clusters. Groups 1 and 4 are the two largest groups in terms of number of taxonomic elements. Also, as in Group 1, a majority of members comes from the Bacteroidales-Bacteroidaceae group that represents almost 2/3 of the species contained in this cluster (17 out of 27 members), it could be related to Bacteroidales-Bacteroidaceae cluster. Group 4 is associated with Lachnospiraceae which represent almost 1/3 of the total in this group (7 out of 23 members). The other three groups (2, 3, and 5) were assigned to the Ruminococcaceae-Bifidobacteriaceae group (5 out of 19 members), Negativicutes group (4 out of 19 members), and Pasteurellaceae group (3 out of 15 members), respectively. The results and relationships between the different elements on each of the clusters are presented in Figure 1. Here, species close to each other in the tree are shown in a similar way. The ordered distance matrix shows the clustering structure. Similar clusters appear as blocks on the diagonal of this heatmap. Darker colours represent small distances whereas the lighter colours represent large distances. The identified clusters have different sizes and compositions, with two large coabundance clusters, grouping the majority of the genus analyzed. It is important to note that we have combined under the name Group 0 all the discarded operational taxonomic units, that is, all the species with more that 70% of zeroes in the sample.Table 1
Composition of the five clusters identified in the analysis of the Indian dataset (Dubey et al., [21]).
Group 1Bacteroidetes-Bacteroidia-Bacteroidales-Rikenellaceae-Alistipes-onderdonkiiBacteroidetes-Bacteroidia-Bacteroidales-Rikenellaceae-Alistipes-putredinisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-coprocolaBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-coprophilusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-doreiBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-fragilisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-intestinalisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-plebeiusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-thetaiotaomicronBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-uniformisBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-vulgatusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-xylanisolvensBacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Barnesiella-intestinihominisFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-biformeFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-ventriosumFirmicutes-Clostridia-Clostridiales-Clostridiales-Flavonifractor-plautiiProteobacteria-Gammaproteobacteria-Enterobacteriales-Enterobacteriaceae-Klebsiella-variicolaFirmicutes-Bacilli-Lactobacillales-Lactobacillaceae-Lactobacillus-rogosaeFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Megasphaera-sp.Bacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Odoribacter-splanchnicusBacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Parabacteroides-distasonisBacteroidetes-Bacteroidia-Bacteroidales-Porphyromonadaceae-Parabacteroides-merdaeBacteroidetes-Bacteroidia-Bacteroidales-Prevotellaceae-Prevotella-copriFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-sp.Firmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-bromiiFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-callidusProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Sutterella-wadsworthensisGroup 2Actinobacteria-Actinomycetales-Actinomycineae-Actinomycetaceae-Actinomyces-odontolyticusBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-massiliensisActinobacteria-Actinobacteria-Bifidobacteriales-Bifidobacteriaceae-Bifidobacterium-adolescentisActinobacteria-Actinobacteria-Bifidobacteriales-Bifidobacteriaceae-Bifidobacterium-bifidumFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-lutiFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-wexleraeFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Butyrivibrio-crossotusActinobacteria-Coriobacteriia-Coriobacteriales-Coriobacteriaceae-Collinsella-aerofaciensFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-eligensFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Faecalibacterium-prausnitziiFirmicutes-Clostridia-Clostridiales-Clostridiales-Howardella-ureilyticaProteobacteria-Gammaproteobacteria-Enterobacteriales-Enterobacteriaceae-Klebsiella-pneumoniaeFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Megasphaera-micronuciformisProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Parasutterella-excrementihominisFirmicutes-Clostridia-Clostridiales-Peptostreptococcaceae-Peptostreptococcus-stomatisFirmicutes-Negativicutes-Acidaminococcales-Acidaminococcaceae-Phascolarctobacterium-faeciumSpirochaetes-Spirochaetes-Spirochaetales-Spirochaetaceae-Treponema-succinifaciensFirmicutes-Erysipelotrichia-Erysipelotrichales-Erysipelotrichaceae-Turicibacter-sanguinisLentisphaerae-Lentisphaeria-Victivallales-Victivallaceae-Victivallis-vadensisGroup 3Proteobacteria-Alphaproteobacteria-Rhodospirillales-Acetobacteraceae-Acidiphilium-sp.Verrucomicrobia-Verrucomicrobiae-Verrucomicrobiales-Akkermansiaceae-Akkermansia-muciniphilaProteobacteria-Gammaproteobacteria-Aeromonadales-Succinivibrionaceae-Anaerobiospirillum-succiniciproducensActinobacteria-Actinobacteria-Bifidobacteriales-Bifidobacteriaceae-Bifidobacterium-longumProteobacteria-Deltaproteobacteria-Desulfovibrionales-Desulfovibrionaceae-Bilophila-wadsworthiaFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-obeumFirmicutes-Erysipelotrichia-Erysipelotrichales-Erysipelotrichaceae-Bulleidia-p-1630-c5Firmicutes-Clostridia-Clostridiales-Lachnospiraceae-Coprococcus-catusProteobacteria-Deltaproteobacteria-Desulfovibrionales-Desulfovibrionaceae-Desulfovibrio-pigerFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Dialister-succinatiphilusFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-siraeumFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Megasphaera-elsdeniiFirmicutes-Negativicutes-Selenomonadales-Selenomonadaceae-Mitsuokella-jalaludiniiFirmicutes-Negativicutes-Selenomonadales-Selenomonadaceae-Mitsuokella-multacidaBacteroidetes-Bacteroidia-Bacteroidales-Prevotellaceae-Prevotella-stercoreaFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Romboutsia-ilealisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminiclostridium-siraeumActinobacteria-Coriobacteriia-Eggerthellales-Eggerthellaceae-Slackia-isoflavoniconvertensProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Sutterella-sp.Group 4Bacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-caccaeBacteroidetes-Bacteroidia-Bacteroidales-Bacteroidaceae-Bacteroides-ovatusFirmicutes-Erysipelotrichia-Erysipelotrichales-Erysipelotrichaceae-Catenibacterium-mitsuokaiFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-bartlettiiFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Coprococcus-comesFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Coprococcus-eutactusFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Dorea-formicigeneransFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Dorea-longicatenaFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-hadrumFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-halliiFirmicutes-Clostridia-Clostridiales-Eubacteriaceae-Eubacterium-ramulusProteobacteria-Alphaproteobacteria-Rhizobiales-Hyphomicrobiaceae-Gemmiger-formicilisFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Lachnoclostridium-clostridioformeFirmicutes-Clostridia-Clostridiales-Oscillospiraceae-Oscillibacter-sp.Firmicutes-Negativicutes-Acidaminococcales-Acidaminococcaceae-Phascolarctobacterium-succinatutensFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-faecisFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-inulinivoransFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-faecisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-gnavusFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-sp.Firmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-torquesFirmicutes-Negativicutes-Selenomonadales-Selenomonadaceae-Selenomonas-bovisProteobacteria-Betaproteobacteria-Burkholderiales-Sutterellaceae-Sutterella-stercoricanisGroup 5Proteobacteria-Gammaproteobacteria-Pasteurellales-Pasteurellaceae-Actinobacillus-minorFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-faecisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Blautia-productaFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-disporicumFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-perfringensFirmicutes-Clostridia-Clostridiales-Clostridiaceae-Clostridium-sp.Proteobacteria-Deltaproteobacteria-Desulfovibrionales-Desulfovibrionaceae-Desulfovibrio-D168Proteobacteria-Gammaproteobacteria-Enterobacteriales-Enterobacteriaceae-Escherichia-coliProteobacteria-Gammaproteobacteria-Pasteurellales-Pasteurellaceae-Haemophilus-parainfluenzaeProteobacteria-Gammaproteobacteria-Pasteurellales-Pasteurellaceae-Haemophilus-pittmaniaeFirmicutes-Bacilli-Lactobacillales-Lactobacillaceae-Lactobacillus-ruminisProteobacteria-Gammaproteobacteria-Pseudomonadales-Pseudomonadaceae-Pseudomonas-liniFirmicutes-Clostridia-Clostridiales-Lachnospiraceae-Roseburia-intestinalisFirmicutes-Clostridia-Clostridiales-Ruminococcaceae-Ruminococcus-gauvreauiiFirmicutes-Negativicutes-Veillonellales-Veillonellaceae-Veillonella-disparFigure 1
Heatmap of coabundance groups. This figure represents the distance among the abundances of the different genus of bacteria characterised in the sequencing analyses. The dashed lines depict the boundaries between the five different clusters. Darker colours represent closer distances (coabundance) between genus whereas lighter colours displays larger distances.Table2 displays the mean, median, and standard deviation for each one of the coabundance groups. It is noticeable that the proportion of bacteria that belongs to Group 1 is higher in average than the proportion in the other groups. The variability is also larger for the first coabundance group.Table 2
Mean, median, and standard deviation for each coabundance group.
Coabundance groupMeanMedianStandard deviationGroup 10.55480.56690.1542Group 20.18990.16700.1024Group 30.08160.06240.0688Group 40.08990.08130.0517Group 50.06030.02390.0859Group 00.02350.01200.0632In Figure2, some ternary plots for different combinations of the bacterial groups are displayed. In particular we have compared the coabundance Group 1, with Group 2 (top left), Group 3 (top right), Group 4 (bottom left), and Group 5 (bottom right). In order to ensure that the total sum is one, we have combined the coabundance proportion for the rest of the groups in each graph as Others. Group 1 is always located at the top of each triangle. The proportion of coabundance of Group 1 is measured in terms of the horizontal lines, i.e., 0% of coabundance is measured in terms of base of the triangle (farthest from the vertex Group 1). In the lower left apex of each triangle is represented the groups compared to Group 1. The right side of the triangle now becomes the baseline for the percentage of the groups located in this vertex. Finally, the combined groups are located at the lower right apex of the triangle.Figure 2
Ternary plots associated to different coabundance groups. (a) Group1/Group 2/Other groups. (b) Group 1/Group 3/Other groups. (c) Group 1/Group 4/Other groups. (d) Group 1/Group 5/Other groups.
(a)(b)(c)(d)The rate of coabundance for the combined groups is calculated from the left side of the triangle (0% abundance) to the lower right corner (100% abundance). It is observable that the data lie from a high amount of coabundance of Group 1 and Group 2 with a low coabundance of third, fourth, and fifth groups (top left graph). From the rest of the graphs, it can be inferred that when Group 1 is compared to the other groups, the coabundance of these groups is lower than in the former graph. Also, as Group 2 has been included in the lotOthers, the corresponding coabundance of the combined group is higher than in the top left graph.


### 2.1.1. Compositional Data Predictor

Compositional data can be defined as arrays of strictly positive numbers for which ratios between them are important without any further requirement [25]. Microbiome data are compositional, that is, the distance between component values is only meaningful proportionally (see [26]). The elements of the composition are non-negative and sum to unity. An important issue in microbiome data is the large presence of zeros; however, the issue of zero values in some components is not addressed in most papers and especially in the task of regression. In general, in compositional research problems, most of the basic statistical analysis tools are incorrect unless the variables are rewritten in terms of logarithms of ratios as proposed in the log-ratio methodology for compositional data. After computing these log-ratios, standard regression methods can be used since the relative character of the information is considered when analyzing the results, as one group or variable can only increase in relative terms if some other group or groups reduce. In this work we focus on the case of compositional data being included in the predictor variables. The effect of increasing one of the variables in relative terms in the predictor therefore depends on which other variables are decreased when this occur. In log-ratio parlance, the effect of increasing one log-ratio is interpreted while keeping all other log-ratios constant as the same log-ratio may have different meaning depending on the way that the other log-ratios in the model are assembled. Thus, the interpretation of log-ratios as explanatory variables is usually different from other approaches. Several different approaches of building and interpreting the log-ratios have been considered in the literature, often leading to the same predictions and residuals [27]. Among the different parametrizations, in this work, we have chosen centred log-ratios [28]. In our analysis, we consider a vector of 6-dimensional real space that carries information on the relative importance of its components,(1)xi=xi0,xi1,…,xi5∈ℝ+6,withxij>0,j=0,1,2,…,5,∑j=05xij=1,where i=1,…,1004. Note that the explanatory variables xij represent the proportion of the bacterial coabundance proportion of Group j in individual i. Centred log-ratios are calculated by using a quotient between each variable and the geometric mean of all components (see [28]),(2)log2xj∏j=05xj6,withj=0,1,2,…,5.The fact that we are using logarithms to base 2 means that a unit increase in this logarithm leads to a double increase in the original magnitude. In order to avoid perfect collinearity one centred log-ratio must be deleted from the regression equation. Since all six centred log-ratios add-up to zero, by increasing a fixed centred log-ratio while keeping the other four remaining log-ratios in the regression equation (with regressorsβj with j=0,1,2,…,5) constant implies increasing the given centred log-ratio whilst reducing the omitted centred log-ratio by the same amount. In this regard, a positive statistically significant regression βj coefficient indicates an increasing value of the covariate xj at the expense of decreasing the amount of the omitted component has a significant positive effect on expected value of the response variable. This is equivalent to say that in terms of logarithm of base 2, βj is interpreted as the expected change in the response variable when the ratio between xj and the omitted explanatory variable is multiplied by six. Finally, in order to obtain the estimates and their corresponding p values for all possible pair combinations, the model needs to be repeated six different times by ignoring each time a different centred log-ratio.


## 2.1.1. Compositional Data Predictor

Compositional data can be defined as arrays of strictly positive numbers for which ratios between them are important without any further requirement [25]. Microbiome data are compositional, that is, the distance between component values is only meaningful proportionally (see [26]). The elements of the composition are non-negative and sum to unity. An important issue in microbiome data is the large presence of zeros; however, the issue of zero values in some components is not addressed in most papers and especially in the task of regression. In general, in compositional research problems, most of the basic statistical analysis tools are incorrect unless the variables are rewritten in terms of logarithms of ratios as proposed in the log-ratio methodology for compositional data. After computing these log-ratios, standard regression methods can be used since the relative character of the information is considered when analyzing the results, as one group or variable can only increase in relative terms if some other group or groups reduce. In this work we focus on the case of compositional data being included in the predictor variables. The effect of increasing one of the variables in relative terms in the predictor therefore depends on which other variables are decreased when this occur. In log-ratio parlance, the effect of increasing one log-ratio is interpreted while keeping all other log-ratios constant as the same log-ratio may have different meaning depending on the way that the other log-ratios in the model are assembled. Thus, the interpretation of log-ratios as explanatory variables is usually different from other approaches. Several different approaches of building and interpreting the log-ratios have been considered in the literature, often leading to the same predictions and residuals [27]. Among the different parametrizations, in this work, we have chosen centred log-ratios [28]. In our analysis, we consider a vector of 6-dimensional real space that carries information on the relative importance of its components,(1)xi=xi0,xi1,…,xi5∈ℝ+6,withxij>0,j=0,1,2,…,5,∑j=05xij=1,where i=1,…,1004. Note that the explanatory variables xij represent the proportion of the bacterial coabundance proportion of Group j in individual i. Centred log-ratios are calculated by using a quotient between each variable and the geometric mean of all components (see [28]),(2)log2xj∏j=05xj6,withj=0,1,2,…,5.The fact that we are using logarithms to base 2 means that a unit increase in this logarithm leads to a double increase in the original magnitude. In order to avoid perfect collinearity one centred log-ratio must be deleted from the regression equation. Since all six centred log-ratios add-up to zero, by increasing a fixed centred log-ratio while keeping the other four remaining log-ratios in the regression equation (with regressorsβj with j=0,1,2,…,5) constant implies increasing the given centred log-ratio whilst reducing the omitted centred log-ratio by the same amount. In this regard, a positive statistically significant regression βj coefficient indicates an increasing value of the covariate xj at the expense of decreasing the amount of the omitted component has a significant positive effect on expected value of the response variable. This is equivalent to say that in terms of logarithm of base 2, βj is interpreted as the expected change in the response variable when the ratio between xj and the omitted explanatory variable is multiplied by six. Finally, in order to obtain the estimates and their corresponding p values for all possible pair combinations, the model needs to be repeated six different times by ignoring each time a different centred log-ratio.


## 3. Relation of Age, BMI, and Level of Physical Activity with Microbiome

In this section, we firstly consider the marginal relation of age given a level of physical activity with microbiome through a beta regression model. Next, we examine the connection of BMI with the microbiome given the level of physical activeness by using a mixture of logistics regression model. Finally, the joint relationship of these variables is examined via at-copula.


### 3.1. Relation of Age and Level of Physical Activity with Microbiome

It is our interest to model the relationship between the age of the subject and the proportion of each coabundance genus-like groups at taxonomic level of species via a beta regression model (see Ferrari and Cribari-Neto [29]). This model assumes that the response variable is beta distributed using a parametrization of the beta law that is indexed by mean and dispersion parameters. This regression family is useful for modeling rates and proportions, that is, in situations where the dependent variable of interest is continuous and restricted to a bounded interval a,b where a and b are known scalars with a<b. This model is related to other variables through a regression structure. Our goal is to explain a continuous response variable Y1 with a<y1<b. The density of Y1 is defined as follows:(3)f1y1i|ωi,ϕ=Γϕb−a1−ϕΓωiϕΓ1−ωiϕy1i−ab−aωiϕ−1b−y1ib−a1−ωiϕ−1,with 0<ωi<1 and ϕ>0 with i=1,…,n. This parametrization allows us to obtain a regression structure for the mean of the response along with a dispersion parameter ϕ. Here, n is the sample size and EY1i=ωib−a+a. The variance of the response variable can be easily explained in terms of its mean by the following expression VarY1i=b−a2ωi1−ωi/1+ϕ. The variance decreases with the value of the dispersion parameter.Let us now consider that a random variableY1i denoting age of the individual i in the sample is related to a compositional data predictor related to each one the coabundance groups, xi=1,ui1,…,ui5⊤ where ui1,…,ui5 are chosen among all combinations without repetition from the vector xi0/∏j=05xj6,…,xi5/∏j=05xj6⊤ taking 5 components at a time. Then, by using the logit link (i.e., hωi=logω1/1−ωi), we have that ωi=expxi⊤β¯/1+expxi⊤β¯, where β¯=β0,β1,…,β5⊤ is a vector of regressors. Other choices for the link function link functions for the response model are feasible.We have fitted this beta regression model to this dataset to explain the response variableAge by considering two levels of physical activity: sedentary and nonsedentary by assuming a=17.5 and b=65.5. Below in Table 3, the estimates and p values associated with the six predictors for individuals classified as sedentary for each microbiome coabundance group’s proportion obtained under the regression model (1). Similarly, in Table 4 estimates and p values results for each predictor are shown for nonsedentary subjects. From these tables, it is discernible that for the first predictor the regressor associated to Group 5 for sedentary individuals is statistically significant at the 5% level whereas it is not for nonsedentary subjects. Its value, −0.1260, is interpreted as the decrease in the covariate x5 (i.e., Pasteurellaceae) at the expense of increasing the amount of x0 has a significant negative effect on hωi while keeping the remaining four log-ratios in the equation constant. In a similar fashion for the third predictor for the nonsedentary subjects the explanatory variables, x1, x3, x4, and x5 are significant at the same level while they are not for the sedentary individuals. For all these covariates, the sign of their regressors is positive; therefore, an increase in these regressors at the expense of decreasing the value of x2 has a significant effect on the transformation of the expected value of the mean of the model. For the fourth predictor, the regression coefficient associated to the second group is only significant for the nonsedentary party. Similarly for the fifth predictor, x5 is significant for the sedentary individuals whereas the regressor associated to Group 2, i.e., Ruminococcaceae-Bifidobacteriaceae is significant for the nonsedentary individuals. Similar situation is also verified for the sixth predictor. On the other hand, for the sedentary subjects, the variable x0 is a positive significant variable.Table 3
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as sedentary. The response variable is age.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 1−0.0341—0.13690.0047−0.00500.09210.6026—0.13880.95930.94730.1262Group 2−0.1708−0.1368—−0.1313−0.1500−0.04480.01170.1388—0.08610.06860.4267Group 3−0.0404−0.00640.1305—−0.00760.08570.494040.94390.0880—0.79570.0956Group 4−0.0448−0.01120.1253−0.0046—0.08090.49920.90170.19230.9515—0.2157Group 5−0.1260−0.09210.0446−0.0863−0.1023—0.00440.12610.42860.09310.0284—Group 0—0.03410.17080.03970.01450.1261—0.60250.01170.50100.63340.0044Intercept−0.1209−0.1210−0.1214−0.1178−0.1088−0.12100.51750.51730.51570.52850.56590.5173ϕ2.97972.97952.97962.97892.98072.9795<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC3481.073481.073481.073481.073480.863481.07Table 4
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as nonsedentary. The response variable is age.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 1−0.0465—0.1209−0.0105−0.07450.01590.4022—0.02830.88510.35830.7459Group 2−0.1673−0.1209—−0.1328−0.1954−0.10520.00630.1151—0.04610.03440.0459Group 3−0.03580.01080.1318—−0.06370.02650.49570.88140.0477—0.40140.5497Group 40.02770.07450.19540.0641—0.08990.66230.35830.03440.3983—0.1771Group 5−0.0623−0.01570.1052−0.0263−0.0902—0.11960.74840.04580.55220.1753—Group 0—0.04650.16740.0349−0.02800.0622—0.40210.00630.50690.65890.1206Intercept−0.1641−0.1642−0.1641−0.1667−0.1643−0.16490.28900.28880.28920.28200.28870.2869ϕ2.90022.89992.90002.89812.89992.8999<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC4155.084155.084155.084155.084155.084155.08In Figure3, we have plotted the histograms of the empirical distribution of the response variable Age for the nonsedentary (top left panel) and sedentary (bottom left panel) subjects. For both histograms, we have superimposed the probability density function of the beta distribution. It is observable that this distribution provides a better fit to empirical data for the nonsedentary group than for the sedentary party. Furthermore, we have performed a diagnostic analysis to check the goodness-of-fit of the estimated model by providing a global measure of explained variation and graphical tools based on QQ-plots, to detect departures from the given model and influential observations. Residuals are used to check the appropriateness of a chosen model and to identify outliers. For that reason, randomized quantile residuals (Dunn and Smyth [30]) are used since other type of residuals, i.e., Pearson’s and deviance residuals are far from normality when the parameters of the model are known and they fail to provide useful information of the inadequacy of the model. The ith randomized quantile residuals for a discrete response variable is defined as rq,i=Φ−1Fy1i;ω^i,ϕ^ where Φ−1⋅ is the quantile function of the standard normal distribution and Fy1i;ω^i,ϕ^ is the cumulative distribution function associated to the beta regression model evaluated at the estimated parameters for i=1,…,n. In the right panels of Figure 3, the QQ-plots of the randomized quantile residuals of the beta regression models when the predictor 1 is considered for the nonsedentary (top right) and sedentary (bottom right) subjects. Each dot on the plots represents an empirical residual. A perfect alignment with the 45° line implies that the residuals are normally distributed. In general, it is observable that the residuals for the nonsedentary group adhere closer to the line in the whole distribution.Figure 3
Histograms of theAge distributions and the beta density superimposed and QQplots of the randomized quantile residuals (RQRs) for each regression model by using Predictor 1. (a) Histogram of age and nonsedentary. (b) RQR of age and nonsedentary against Predictor 1 regression model. (c) Histogram of age and sedentary. (d) RQR of age and sedentary against Predictor 1 regression model.
(a)(b)(c)(d)


### 3.2. Relation of BMI and Level of Physical Activity with Microbiome

Regardless of the gender, the empirical distribution of BMI in humans is bimodal. Then, finding appropriate statistical models that have the capacity to explain bimodal datasets is an issue of vital importance. In this work, we use a mixture of two logistic distributions with different locations and scale parameters. We have chosen this family for its flexibility and simplicity. It is now our interest to explain the BMI in the population in terms of a random variableY2∈ℝ. The probability density function of this random variable is(4)f2y2i|wi,μ1,σ1,μ2,σ2=wi4σ1sech2y2i−μ12σ1,(5)+1−wi4σ2sech2y2i−μ22σ2,where μ1, μ2∈ℝ are location parameters and σ1, σ2>0 are scale parameters. In this model, it will also be assumed that the weight parameter for each individual in the sample is again expressed as a function of the same group of covariates, wi=expxi⊤α¯/1+expxi⊤α¯, where α¯=α0,α1,…,α5⊤ is a vector of regressors. Other choices for the link function link functions for the response model are also possible. This parametrization enable us to obtain a regression structure for the mean of the response in the following way, EY2i=wiμ1+1−wiμ2. The variance of the response variable is written in terms of the following linear combination of the scale parameters:(6)VarY2i=π23wiσ12+1−wiσ22.We have now fitted the mixture of logistics regression model given by (5) to this dataset to explain the dependent variable BMI by again considering two levels of physical activity: sedentary and nonsedentary. In Table 5, the estimates and p values associated with the six predictors for individuals classified as sedentary for each microbiome coabundance groups proportion obtained under this mixture of logistics regression model (5). In a similar way, in Table 6, estimates and p values results for each predictor are shown for nonsedentary subjects. From these tables, it is apparent that for the first predictor the regressor associated to Group 1 for sedentary individuals is statistically significant at the 5% level whereas it is not for non-sedentary subjects. The estimated value is −0.3631, that it is interpreted as the decrease in the covariate x1 (i.e., Bacteroidales-Bacteroidaceae) at the expense of increasing the amount of x0 has a significant negative effect on hωi while keeping the remaining four log-ratios in the equation of the predictor constant. Similarly, for the second predictor and the sedentary subjects, the explanatory variables, x0 is statistically significant at the 10% significance level while it is not for the nonsedentary individuals. The sign of this regression coefficient is positive; therefore, an increase in this regressor at the expense of decreasing the value of x1, while keeping the other four log-ratios in the equation constant has a significant effect on the transformation of the expected value of the mean of the model. Finally, for predictor 3, the regression coefficient associated to the first group is only significant at the 10% significance level for the sedentary individuals. The sign of this regressor is negative.Table 5
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as sedentary. The response variable is BMI.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 1−0.3631—−0.5055−0.2915−0.1853−0.18290.0492—0.05340.23930.47110.2635Group 20.22000.3914—0.21540.32150.14210.24290.1217—0.31170.21730.3455Group 3−0.02670.2575−0.2111—0.1062−0.00710.87160.29660.3209—0.61990.9604Group 4−0.17770.0996−0.3113−0.1064—−0.11450.34500.69390.22810.6195—0.5254Group 5−0.05270.2175−0.2280−0.01810.0881—0.66480.16410.13790.90040.6246—Group 0—0.3143−0.12920.08490.19110.1057—0.08190.48380.60610.30490.3826Intercept0.31010.38940.36260.36330.36340.33740.56880.48730.49830.49890.49890.5406μ124.492324.7000024.425024.453724.453224.6451<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.347131.468331.350531.361631.361331.4405<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s12.12692.17612.09052.09842.09842.1679<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.66601.61231.66641.66431.66441.6371<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC2717.152717.252716.832716.832716.832717.54Table 6
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as nonsedentary. The response variable is BMI.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 10.0419—0.14530.02640.21680.05260.6145—0.45000.88470.26960.6742Group 2−0.1095−0.1466—−0.11850.0662−0.12770.46240.4459—0.46610.76540.3291Group 30.0164−0.02820.1186—0.17930.02020.89800.87690.4660—0.32580.8540Group 4−0.1438−0.1902−0.0425−0.1618—−0.16690.34640.33220.84810.3746—0.2921Group 50.0084−0.02520.12040.00220.1827—0.92950.83660.35580.98430.2488—Group 0—−0.03280.1126−0.00560.1725−0.0027—0.80680.44970.96500.25940.9771Intercept0.20750.22820.23080.23070.19020.17310.61450.57980.57540.57540.64380.6438μ123.906523.892323.895323.892123.890623.8869<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.260831.258231.260631.257631.256331.2567<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s11.94781.94441.94341.94171.94081.9418<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.47691.47851.47791.47881.47901.4803<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC3264.493264.453264.453264.453264.473264.58In Figure4, we have plotted the histograms of the empirical distribution of the response variable BMI for the nonsedentary (top left panel) and sedentary (bottom left panel) individuals. For both histograms, we have superimposed the density function of the mixture of logistics distributions. It can be seen that this distribution is able to reproduce the two modes of the empirical distribution for both cohorts. Note that for the group of sedentary individuals, the second modal value located around the BMI value of 31 is clearly more predominant. Once again, we have plotted the QQ-plots of the randomized quantile residuals of this mixture of logistics regression when the first predictor 1 is considered for the nonsedentary (top right) and sedentary (bottom right) individuals. In general, it is observable that the residuals for the nonsedentary group adhere closer to the line in the whole distribution but it underestimates the top part of the distribution of residuals.Figure 4
Histograms of the BMI distributions and the mixture of logistic density superimposed and QQplots of the randomized quantile residuals (RQR) for each regression model by using Predictor 1. (a) Histogram of BMI and nonsedentary. (b) RQR of BMI and nonsedentary against Predictor 1 regression model. (c) Histogram of BMI and sedentary. (d) RQR of BMI and sedentary against Predictor 1 regression model.
(a)(b)(c)(d)


### 3.3. Joint Relation of Age, BMI, and Level of Physical Activity with Microbiome

The degree of association between the two variables age and BMI in the sample for the different levels of physical activity can be summarized in terms of some measures of correlation for bivariate data. In Table7, Pearson’s, Spearman’s, and Kendall’s measures of correlation for these continuous random variables are displayed. It is noticeable that there exists weak positive correlation between these two variables. The degree of association is less intense for the nonsedentary individuals.Table 7
Pearson’s, Spearman’s, and Kendall’s measures of correlation for the variables age and BMI and different levels of physical activity.
Measure of correlationPhysical activeness levelSedentaryNonsedentaryPearson’s0.22380.1738Spearman’s0.20600.1828Kendall’s0.14210.1242We model the joint dependence of age and BMI for different level of physical activity and their relationship with the proportion of each coabundance genus-like groups at taxonomic level of species via at-copula with degrees of freedom (df) parameter ν>1 with marginal distributions given by the beta regression model given in (1) and the mixture of logistics distributions provided in (5). The density of this of this multivariate distribution is defined as(7)gy¯i|Θ1,Θ2,Θ3=Γν+1/2Γν/21+z¯Σ−1z¯⊤/ν−ν+2/2Γν+1/22Σ1/2∏j=121+zj2/ν−ν+1/2×f1yi1|Θ1×f2yi2|Θ2,where y¯i:=yi1,yi2⊤, Θ1=β¯,ϕ, Θ2=α¯,μ1,σ1,μ2,σ2, and Θ3=ν,Σ. Here, Σ:=ρij1≤i,j≤2 is a symmetric and positive definite scatter matrix with dimension 2×2 with unit diagonal entries and −1<ρij<1, ⋅ denotes the determinant of a matrix, and Γ⋅ is the complete gamma function. Also, z¯:=z1,z2⊤ with zj=tν−1Fjyj|Θj with j=1,2, where tν−1⋅ is the quantile function of univariate t-distribution with ν df and Fjyij|Θj is the cdf associated to the regression models presented above with i=1,…,n.The corresponding log-likelihood function, given a sampley:=y¯1,…,y¯n is provided by(8)ℓΘ1,Θ2,Θ3|y=ℓCΘ1,Θ2,Θ3|y+ℓMΘ1,Θ2|y=nlogΓν+12+nlogΓν2−2nlogΓν+12−n2logΣ−ν+22∑i=1nlog1+z¯iΣ−1z¯i⊤ν+ν+12∑i=1n∑j=12log1+zij2ν+nlogΓϕ+n1−ϕlogb−a−∑i=1nlogΓωiϕ+logΓ1−ωiϕ+∑i=1nωiϕ−1logyi−ab−a+∑i=1n1−ωiϕ−1logb−yib−a+∑i=1nlogwi4σ1sech2y2i−μ12σ1+1−wi4σ2sech2y2i−μ22σ2,where ℓC⋅ and ℓM⋅ are the log-likelihood functions of the copula and marginal model, respectively. Maximum likelihood estimation can be used to estimate the parameters of expression (8) via an adaptive maximization by parts (MBP) algorithm as described in [31], by using initial estimates Θ^10,Θ^20,Θ^30 generated by inference for margins algorithm. In the step k of this algorithm for k=1,2,…, we find,(9)Θ^1k,Θ^2k=argmaxℓMΘ1,Θ2|y+ℓCΘ1,Θ2|y,Θ^3k−1,Θ3k=argmaxℓCΘ3|y,Θ^1k,Θ^2k.The algorithm stops when a terminating condition between two consecutive iterations is reached, i.e.,(10)ℓΘ^1k,Θ^2k,Θ^3k|y−ℓΘ^1k−1,Θ^2k−1,Θ^3k−1|y1<10−3.Finally, we have fitted the bivariate distribution given in (5) to the bivariate data set. Once again, the two levels of physical activeness have been considered. Results are shown in Table 8 for the sedentary case and Table 9 for the nonsedentary situation. When using the first predictor, i.e., the omitted covariate is x0, the variable x5 and x1 are statistically significant at the 5% level of significance for the variables age and BMI, respectively, for sedentary individuals while they are not for the nonsedentary group. Also, the regressor associated with the covariates x2 and x0 for the BMI are significant at the same level of significance for the sedentary individuals whereas they are not for the nonsedentary group. The sign of these regression coefficients is negative. Conversely, for the third predictor and the variable Age, the regressors for the variables x3 and x5 are positive significant only for the nonsedentary subjects. In addition, the variable x3 is positive significant for the variable BMI for the same level of physical activeness. Regression coefficients associated to x2 in the fourth and fifth predictors (omitted variables x3 and x4, respectively) are negative significant for the response variable age. Finally, when the explanatory variable x5 is deleted, the regressor linked to x0 is positive significant at the 5% level for the sedentary subjects and response variable Age. Similarly, for this sixth predictor, the regressor associated to the variable x2 for the same response variable is negative significant for individuals classified as nonsedentary.Table 8
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as sedentary (bivariate regression model).
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6AgeGroup 1−0.0051—0.0603−0.00000.03460.08570.9316—0.50810.99970.69750.1475Group 2−0.1580−0.1524—−0.1233−0.0983−0.04370.01780.0944—0.10110.29850.4303Group 3−0.0462−0.02930.1203—0.01670.08240.42650.74360.1103—0.82320.1031Group 4−0.0430−0.03710.0880−0.0065—0.07650.51010.67610.01160.0745—0.2339Group 5−0.1125−0.10450.0161−0.0825−0.0603—0.00990.07790.77110.10270.3485—Group 0—0.01260.14630.04400.06450.1257—0.84480.02820.44800.32380.0039Intercept−0.2183−0.1924−0.0095−0.1094−0.1543−0.11740.23610.29580.95870.55170.40140.5226ϕ3.12413.12613.11433.13543.12583.1411<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001BMIGroup 1−0.3887—−0.4573−0.2960−0.1711−0.25300.0292—0.06420.21540.48500.1132Group 20.10770.4876—0.18270.31420.22270.54290.0500—0.36800.21200.1314Group 3−0.06470.3060−0.1777—0.13550.04210.68350.20180.3813—0.50900.7622Group 4−0.20220.1727−0.3104−0.1363—−0.09200.25980.48130.21690.5057—0.5961Group 5−0.11690.2616−0.2148−0.04170.0902—0.31820.10230.14440.76400.6033—Group 0—0.3713−0.11110.06400.19790.1079—0.03620.52860.68570.26850.3546Intercept0.45460.39990.34520.37000.39620.37640.37940.43650.49950.46920.44030.4631μ124.537024.533724.552124.530924.534524.5355<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.451131.447731.455731.450931.448631.4552<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s11.97931.97641.97861.97081.97521.9791<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.55621.55761.55301.55261.55721.5599<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ρ120.25540.25500.25710.25570.25580.2557<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ν17.115117.104216.929516.814317.065316.99710.00500.00500.00460.00430.00490.0047Table 9
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as nonsedentary (bivariate regression model).
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6AgeGroup 1−0.0256—0.1122−0.0160−0.07830.01230.6443—0.14190.82520.33300.8009Group 2−0.1632−0.1279—−0.1313−0.1980−0.10320.00760.0952—0.04780.03170.0493Group 3−0.03240.00750.1344—−0.06500.02830.53600.91740.0425—0.39060.5205Group 40.03530.07490.19170.0617—0.09130.57670.35500.03710.4146—0.1690Group 5−0.0546−0.02050.1022−0.0299−0.0950—0.17140.67480.05130.49800.1527—Group 0—0.04200.16400.0323−0.03280.0622—0.44800.00720.53830.60370.1197Intercept−0.2207−0.1789−0.1543−0.1667−0.1726−0.15950.15320.24710.31700.28020.26360.3015ϕ2.92312.92152.93572.92692.92342.9294<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001BMIGroup 10.0183—0.11220.01820.18860.01390.8906—0.14190.91980.33380.9090Group 2−0.1063−0.1419—−0.11200.0545−0.11030.47300.4582—0.48830.80520.3944Group 30.0070−0.02860.1344—0.16820.00030.95610.87460.0426—0.35370.9933Group 4−0.1635−0.1953−0.0549−0.1703—−0.16440.28250.31680.80360.3474—0.2959Group 50.0037−0.02810.1133−0.00230.1680—0.96870.81680.38150.98310.2856—Group 0—−0.03440.1074−0.00670.1623−0.0054—0.79630.46820.95760.28600.9546Intercept0.25180.20810.22630.23520.23270.24710.53800.61000.57900.56410.56890.5455μ123.880923.874223.881823.886623.879323.8795<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.266231.259331.264731.271631.266731.2640<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s11.92261.92030.91881.91691.92241.9238<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.46531.46561.46241.46341.46651.4655<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ρ120.19700.19680.19650.19680.19700.1967<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ν114.9900110.699098.4523103.0250110.6760106.21800.59780.58400.53930.58790.58160.5676


## 3.1. Relation of Age and Level of Physical Activity with Microbiome

It is our interest to model the relationship between the age of the subject and the proportion of each coabundance genus-like groups at taxonomic level of species via a beta regression model (see Ferrari and Cribari-Neto [29]). This model assumes that the response variable is beta distributed using a parametrization of the beta law that is indexed by mean and dispersion parameters. This regression family is useful for modeling rates and proportions, that is, in situations where the dependent variable of interest is continuous and restricted to a bounded interval a,b where a and b are known scalars with a<b. This model is related to other variables through a regression structure. Our goal is to explain a continuous response variable Y1 with a<y1<b. The density of Y1 is defined as follows:(3)f1y1i|ωi,ϕ=Γϕb−a1−ϕΓωiϕΓ1−ωiϕy1i−ab−aωiϕ−1b−y1ib−a1−ωiϕ−1,with 0<ωi<1 and ϕ>0 with i=1,…,n. This parametrization allows us to obtain a regression structure for the mean of the response along with a dispersion parameter ϕ. Here, n is the sample size and EY1i=ωib−a+a. The variance of the response variable can be easily explained in terms of its mean by the following expression VarY1i=b−a2ωi1−ωi/1+ϕ. The variance decreases with the value of the dispersion parameter.Let us now consider that a random variableY1i denoting age of the individual i in the sample is related to a compositional data predictor related to each one the coabundance groups, xi=1,ui1,…,ui5⊤ where ui1,…,ui5 are chosen among all combinations without repetition from the vector xi0/∏j=05xj6,…,xi5/∏j=05xj6⊤ taking 5 components at a time. Then, by using the logit link (i.e., hωi=logω1/1−ωi), we have that ωi=expxi⊤β¯/1+expxi⊤β¯, where β¯=β0,β1,…,β5⊤ is a vector of regressors. Other choices for the link function link functions for the response model are feasible.We have fitted this beta regression model to this dataset to explain the response variableAge by considering two levels of physical activity: sedentary and nonsedentary by assuming a=17.5 and b=65.5. Below in Table 3, the estimates and p values associated with the six predictors for individuals classified as sedentary for each microbiome coabundance group’s proportion obtained under the regression model (1). Similarly, in Table 4 estimates and p values results for each predictor are shown for nonsedentary subjects. From these tables, it is discernible that for the first predictor the regressor associated to Group 5 for sedentary individuals is statistically significant at the 5% level whereas it is not for nonsedentary subjects. Its value, −0.1260, is interpreted as the decrease in the covariate x5 (i.e., Pasteurellaceae) at the expense of increasing the amount of x0 has a significant negative effect on hωi while keeping the remaining four log-ratios in the equation constant. In a similar fashion for the third predictor for the nonsedentary subjects the explanatory variables, x1, x3, x4, and x5 are significant at the same level while they are not for the sedentary individuals. For all these covariates, the sign of their regressors is positive; therefore, an increase in these regressors at the expense of decreasing the value of x2 has a significant effect on the transformation of the expected value of the mean of the model. For the fourth predictor, the regression coefficient associated to the second group is only significant for the nonsedentary party. Similarly for the fifth predictor, x5 is significant for the sedentary individuals whereas the regressor associated to Group 2, i.e., Ruminococcaceae-Bifidobacteriaceae is significant for the nonsedentary individuals. Similar situation is also verified for the sixth predictor. On the other hand, for the sedentary subjects, the variable x0 is a positive significant variable.Table 3
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as sedentary. The response variable is age.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 1−0.0341—0.13690.0047−0.00500.09210.6026—0.13880.95930.94730.1262Group 2−0.1708−0.1368—−0.1313−0.1500−0.04480.01170.1388—0.08610.06860.4267Group 3−0.0404−0.00640.1305—−0.00760.08570.494040.94390.0880—0.79570.0956Group 4−0.0448−0.01120.1253−0.0046—0.08090.49920.90170.19230.9515—0.2157Group 5−0.1260−0.09210.0446−0.0863−0.1023—0.00440.12610.42860.09310.0284—Group 0—0.03410.17080.03970.01450.1261—0.60250.01170.50100.63340.0044Intercept−0.1209−0.1210−0.1214−0.1178−0.1088−0.12100.51750.51730.51570.52850.56590.5173ϕ2.97972.97952.97962.97892.98072.9795<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC3481.073481.073481.073481.073480.863481.07Table 4
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as nonsedentary. The response variable is age.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 1−0.0465—0.1209−0.0105−0.07450.01590.4022—0.02830.88510.35830.7459Group 2−0.1673−0.1209—−0.1328−0.1954−0.10520.00630.1151—0.04610.03440.0459Group 3−0.03580.01080.1318—−0.06370.02650.49570.88140.0477—0.40140.5497Group 40.02770.07450.19540.0641—0.08990.66230.35830.03440.3983—0.1771Group 5−0.0623−0.01570.1052−0.0263−0.0902—0.11960.74840.04580.55220.1753—Group 0—0.04650.16740.0349−0.02800.0622—0.40210.00630.50690.65890.1206Intercept−0.1641−0.1642−0.1641−0.1667−0.1643−0.16490.28900.28880.28920.28200.28870.2869ϕ2.90022.89992.90002.89812.89992.8999<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC4155.084155.084155.084155.084155.084155.08In Figure3, we have plotted the histograms of the empirical distribution of the response variable Age for the nonsedentary (top left panel) and sedentary (bottom left panel) subjects. For both histograms, we have superimposed the probability density function of the beta distribution. It is observable that this distribution provides a better fit to empirical data for the nonsedentary group than for the sedentary party. Furthermore, we have performed a diagnostic analysis to check the goodness-of-fit of the estimated model by providing a global measure of explained variation and graphical tools based on QQ-plots, to detect departures from the given model and influential observations. Residuals are used to check the appropriateness of a chosen model and to identify outliers. For that reason, randomized quantile residuals (Dunn and Smyth [30]) are used since other type of residuals, i.e., Pearson’s and deviance residuals are far from normality when the parameters of the model are known and they fail to provide useful information of the inadequacy of the model. The ith randomized quantile residuals for a discrete response variable is defined as rq,i=Φ−1Fy1i;ω^i,ϕ^ where Φ−1⋅ is the quantile function of the standard normal distribution and Fy1i;ω^i,ϕ^ is the cumulative distribution function associated to the beta regression model evaluated at the estimated parameters for i=1,…,n. In the right panels of Figure 3, the QQ-plots of the randomized quantile residuals of the beta regression models when the predictor 1 is considered for the nonsedentary (top right) and sedentary (bottom right) subjects. Each dot on the plots represents an empirical residual. A perfect alignment with the 45° line implies that the residuals are normally distributed. In general, it is observable that the residuals for the nonsedentary group adhere closer to the line in the whole distribution.Figure 3
Histograms of theAge distributions and the beta density superimposed and QQplots of the randomized quantile residuals (RQRs) for each regression model by using Predictor 1. (a) Histogram of age and nonsedentary. (b) RQR of age and nonsedentary against Predictor 1 regression model. (c) Histogram of age and sedentary. (d) RQR of age and sedentary against Predictor 1 regression model.
(a)(b)(c)(d)


## 3.2. Relation of BMI and Level of Physical Activity with Microbiome

Regardless of the gender, the empirical distribution of BMI in humans is bimodal. Then, finding appropriate statistical models that have the capacity to explain bimodal datasets is an issue of vital importance. In this work, we use a mixture of two logistic distributions with different locations and scale parameters. We have chosen this family for its flexibility and simplicity. It is now our interest to explain the BMI in the population in terms of a random variableY2∈ℝ. The probability density function of this random variable is(4)f2y2i|wi,μ1,σ1,μ2,σ2=wi4σ1sech2y2i−μ12σ1,(5)+1−wi4σ2sech2y2i−μ22σ2,where μ1, μ2∈ℝ are location parameters and σ1, σ2>0 are scale parameters. In this model, it will also be assumed that the weight parameter for each individual in the sample is again expressed as a function of the same group of covariates, wi=expxi⊤α¯/1+expxi⊤α¯, where α¯=α0,α1,…,α5⊤ is a vector of regressors. Other choices for the link function link functions for the response model are also possible. This parametrization enable us to obtain a regression structure for the mean of the response in the following way, EY2i=wiμ1+1−wiμ2. The variance of the response variable is written in terms of the following linear combination of the scale parameters:(6)VarY2i=π23wiσ12+1−wiσ22.We have now fitted the mixture of logistics regression model given by (5) to this dataset to explain the dependent variable BMI by again considering two levels of physical activity: sedentary and nonsedentary. In Table 5, the estimates and p values associated with the six predictors for individuals classified as sedentary for each microbiome coabundance groups proportion obtained under this mixture of logistics regression model (5). In a similar way, in Table 6, estimates and p values results for each predictor are shown for nonsedentary subjects. From these tables, it is apparent that for the first predictor the regressor associated to Group 1 for sedentary individuals is statistically significant at the 5% level whereas it is not for non-sedentary subjects. The estimated value is −0.3631, that it is interpreted as the decrease in the covariate x1 (i.e., Bacteroidales-Bacteroidaceae) at the expense of increasing the amount of x0 has a significant negative effect on hωi while keeping the remaining four log-ratios in the equation of the predictor constant. Similarly, for the second predictor and the sedentary subjects, the explanatory variables, x0 is statistically significant at the 10% significance level while it is not for the nonsedentary individuals. The sign of this regression coefficient is positive; therefore, an increase in this regressor at the expense of decreasing the value of x1, while keeping the other four log-ratios in the equation constant has a significant effect on the transformation of the expected value of the mean of the model. Finally, for predictor 3, the regression coefficient associated to the first group is only significant at the 10% significance level for the sedentary individuals. The sign of this regressor is negative.Table 5
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as sedentary. The response variable is BMI.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 1−0.3631—−0.5055−0.2915−0.1853−0.18290.0492—0.05340.23930.47110.2635Group 20.22000.3914—0.21540.32150.14210.24290.1217—0.31170.21730.3455Group 3−0.02670.2575−0.2111—0.1062−0.00710.87160.29660.3209—0.61990.9604Group 4−0.17770.0996−0.3113−0.1064—−0.11450.34500.69390.22810.6195—0.5254Group 5−0.05270.2175−0.2280−0.01810.0881—0.66480.16410.13790.90040.6246—Group 0—0.3143−0.12920.08490.19110.1057—0.08190.48380.60610.30490.3826Intercept0.31010.38940.36260.36330.36340.33740.56880.48730.49830.49890.49890.5406μ124.492324.7000024.425024.453724.453224.6451<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.347131.468331.350531.361631.361331.4405<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s12.12692.17612.09052.09842.09842.1679<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.66601.61231.66641.66431.66441.6371<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC2717.152717.252716.832716.832716.832717.54Table 6
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as nonsedentary. The response variable is BMI.
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6Group 10.0419—0.14530.02640.21680.05260.6145—0.45000.88470.26960.6742Group 2−0.1095−0.1466—−0.11850.0662−0.12770.46240.4459—0.46610.76540.3291Group 30.0164−0.02820.1186—0.17930.02020.89800.87690.4660—0.32580.8540Group 4−0.1438−0.1902−0.0425−0.1618—−0.16690.34640.33220.84810.3746—0.2921Group 50.0084−0.02520.12040.00220.1827—0.92950.83660.35580.98430.2488—Group 0—−0.03280.1126−0.00560.1725−0.0027—0.80680.44970.96500.25940.9771Intercept0.20750.22820.23080.23070.19020.17310.61450.57980.57540.57540.64380.6438μ123.906523.892323.895323.892123.890623.8869<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.260831.258231.260631.257631.256331.2567<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s11.94781.94441.94341.94171.94081.9418<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.47691.47851.47791.47881.47901.4803<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001AIC3264.493264.453264.453264.453264.473264.58In Figure4, we have plotted the histograms of the empirical distribution of the response variable BMI for the nonsedentary (top left panel) and sedentary (bottom left panel) individuals. For both histograms, we have superimposed the density function of the mixture of logistics distributions. It can be seen that this distribution is able to reproduce the two modes of the empirical distribution for both cohorts. Note that for the group of sedentary individuals, the second modal value located around the BMI value of 31 is clearly more predominant. Once again, we have plotted the QQ-plots of the randomized quantile residuals of this mixture of logistics regression when the first predictor 1 is considered for the nonsedentary (top right) and sedentary (bottom right) individuals. In general, it is observable that the residuals for the nonsedentary group adhere closer to the line in the whole distribution but it underestimates the top part of the distribution of residuals.Figure 4
Histograms of the BMI distributions and the mixture of logistic density superimposed and QQplots of the randomized quantile residuals (RQR) for each regression model by using Predictor 1. (a) Histogram of BMI and nonsedentary. (b) RQR of BMI and nonsedentary against Predictor 1 regression model. (c) Histogram of BMI and sedentary. (d) RQR of BMI and sedentary against Predictor 1 regression model.
(a)(b)(c)(d)


## 3.3. Joint Relation of Age, BMI, and Level of Physical Activity with Microbiome

The degree of association between the two variables age and BMI in the sample for the different levels of physical activity can be summarized in terms of some measures of correlation for bivariate data. In Table7, Pearson’s, Spearman’s, and Kendall’s measures of correlation for these continuous random variables are displayed. It is noticeable that there exists weak positive correlation between these two variables. The degree of association is less intense for the nonsedentary individuals.Table 7
Pearson’s, Spearman’s, and Kendall’s measures of correlation for the variables age and BMI and different levels of physical activity.
Measure of correlationPhysical activeness levelSedentaryNonsedentaryPearson’s0.22380.1738Spearman’s0.20600.1828Kendall’s0.14210.1242We model the joint dependence of age and BMI for different level of physical activity and their relationship with the proportion of each coabundance genus-like groups at taxonomic level of species via at-copula with degrees of freedom (df) parameter ν>1 with marginal distributions given by the beta regression model given in (1) and the mixture of logistics distributions provided in (5). The density of this of this multivariate distribution is defined as(7)gy¯i|Θ1,Θ2,Θ3=Γν+1/2Γν/21+z¯Σ−1z¯⊤/ν−ν+2/2Γν+1/22Σ1/2∏j=121+zj2/ν−ν+1/2×f1yi1|Θ1×f2yi2|Θ2,where y¯i:=yi1,yi2⊤, Θ1=β¯,ϕ, Θ2=α¯,μ1,σ1,μ2,σ2, and Θ3=ν,Σ. Here, Σ:=ρij1≤i,j≤2 is a symmetric and positive definite scatter matrix with dimension 2×2 with unit diagonal entries and −1<ρij<1, ⋅ denotes the determinant of a matrix, and Γ⋅ is the complete gamma function. Also, z¯:=z1,z2⊤ with zj=tν−1Fjyj|Θj with j=1,2, where tν−1⋅ is the quantile function of univariate t-distribution with ν df and Fjyij|Θj is the cdf associated to the regression models presented above with i=1,…,n.The corresponding log-likelihood function, given a sampley:=y¯1,…,y¯n is provided by(8)ℓΘ1,Θ2,Θ3|y=ℓCΘ1,Θ2,Θ3|y+ℓMΘ1,Θ2|y=nlogΓν+12+nlogΓν2−2nlogΓν+12−n2logΣ−ν+22∑i=1nlog1+z¯iΣ−1z¯i⊤ν+ν+12∑i=1n∑j=12log1+zij2ν+nlogΓϕ+n1−ϕlogb−a−∑i=1nlogΓωiϕ+logΓ1−ωiϕ+∑i=1nωiϕ−1logyi−ab−a+∑i=1n1−ωiϕ−1logb−yib−a+∑i=1nlogwi4σ1sech2y2i−μ12σ1+1−wi4σ2sech2y2i−μ22σ2,where ℓC⋅ and ℓM⋅ are the log-likelihood functions of the copula and marginal model, respectively. Maximum likelihood estimation can be used to estimate the parameters of expression (8) via an adaptive maximization by parts (MBP) algorithm as described in [31], by using initial estimates Θ^10,Θ^20,Θ^30 generated by inference for margins algorithm. In the step k of this algorithm for k=1,2,…, we find,(9)Θ^1k,Θ^2k=argmaxℓMΘ1,Θ2|y+ℓCΘ1,Θ2|y,Θ^3k−1,Θ3k=argmaxℓCΘ3|y,Θ^1k,Θ^2k.The algorithm stops when a terminating condition between two consecutive iterations is reached, i.e.,(10)ℓΘ^1k,Θ^2k,Θ^3k|y−ℓΘ^1k−1,Θ^2k−1,Θ^3k−1|y1<10−3.Finally, we have fitted the bivariate distribution given in (5) to the bivariate data set. Once again, the two levels of physical activeness have been considered. Results are shown in Table 8 for the sedentary case and Table 9 for the nonsedentary situation. When using the first predictor, i.e., the omitted covariate is x0, the variable x5 and x1 are statistically significant at the 5% level of significance for the variables age and BMI, respectively, for sedentary individuals while they are not for the nonsedentary group. Also, the regressor associated with the covariates x2 and x0 for the BMI are significant at the same level of significance for the sedentary individuals whereas they are not for the nonsedentary group. The sign of these regression coefficients is negative. Conversely, for the third predictor and the variable Age, the regressors for the variables x3 and x5 are positive significant only for the nonsedentary subjects. In addition, the variable x3 is positive significant for the variable BMI for the same level of physical activeness. Regression coefficients associated to x2 in the fourth and fifth predictors (omitted variables x3 and x4, respectively) are negative significant for the response variable age. Finally, when the explanatory variable x5 is deleted, the regressor linked to x0 is positive significant at the 5% level for the sedentary subjects and response variable Age. Similarly, for this sixth predictor, the regressor associated to the variable x2 for the same response variable is negative significant for individuals classified as nonsedentary.Table 8
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as sedentary (bivariate regression model).
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6AgeGroup 1−0.0051—0.0603−0.00000.03460.08570.9316—0.50810.99970.69750.1475Group 2−0.1580−0.1524—−0.1233−0.0983−0.04370.01780.0944—0.10110.29850.4303Group 3−0.0462−0.02930.1203—0.01670.08240.42650.74360.1103—0.82320.1031Group 4−0.0430−0.03710.0880−0.0065—0.07650.51010.67610.01160.0745—0.2339Group 5−0.1125−0.10450.0161−0.0825−0.0603—0.00990.07790.77110.10270.3485—Group 0—0.01260.14630.04400.06450.1257—0.84480.02820.44800.32380.0039Intercept−0.2183−0.1924−0.0095−0.1094−0.1543−0.11740.23610.29580.95870.55170.40140.5226ϕ3.12413.12613.11433.13543.12583.1411<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001BMIGroup 1−0.3887—−0.4573−0.2960−0.1711−0.25300.0292—0.06420.21540.48500.1132Group 20.10770.4876—0.18270.31420.22270.54290.0500—0.36800.21200.1314Group 3−0.06470.3060−0.1777—0.13550.04210.68350.20180.3813—0.50900.7622Group 4−0.20220.1727−0.3104−0.1363—−0.09200.25980.48130.21690.5057—0.5961Group 5−0.11690.2616−0.2148−0.04170.0902—0.31820.10230.14440.76400.6033—Group 0—0.3713−0.11110.06400.19790.1079—0.03620.52860.68570.26850.3546Intercept0.45460.39990.34520.37000.39620.37640.37940.43650.49950.46920.44030.4631μ124.537024.533724.552124.530924.534524.5355<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.451131.447731.455731.450931.448631.4552<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s11.97931.97641.97861.97081.97521.9791<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.55621.55761.55301.55261.55721.5599<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ρ120.25540.25500.25710.25570.25580.2557<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ν17.115117.104216.929516.814317.065316.99710.00500.00500.00460.00430.00490.0047Table 9
Parameter estimates (first row) andp values (second row) for the regressors associated with the six predictors for individuals classified as nonsedentary (bivariate regression model).
Predictor 1Predictor 2Predictor 3Predictor 4Predictor 5Predictor 6AgeGroup 1−0.0256—0.1122−0.0160−0.07830.01230.6443—0.14190.82520.33300.8009Group 2−0.1632−0.1279—−0.1313−0.1980−0.10320.00760.0952—0.04780.03170.0493Group 3−0.03240.00750.1344—−0.06500.02830.53600.91740.0425—0.39060.5205Group 40.03530.07490.19170.0617—0.09130.57670.35500.03710.4146—0.1690Group 5−0.0546−0.02050.1022−0.0299−0.0950—0.17140.67480.05130.49800.1527—Group 0—0.04200.16400.0323−0.03280.0622—0.44800.00720.53830.60370.1197Intercept−0.2207−0.1789−0.1543−0.1667−0.1726−0.15950.15320.24710.31700.28020.26360.3015ϕ2.92312.92152.93572.92692.92342.9294<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001BMIGroup 10.0183—0.11220.01820.18860.01390.8906—0.14190.91980.33380.9090Group 2−0.1063−0.1419—−0.11200.0545−0.11030.47300.4582—0.48830.80520.3944Group 30.0070−0.02860.1344—0.16820.00030.95610.87460.0426—0.35370.9933Group 4−0.1635−0.1953−0.0549−0.1703—−0.16440.28250.31680.80360.3474—0.2959Group 50.0037−0.02810.1133−0.00230.1680—0.96870.81680.38150.98310.2856—Group 0—−0.03440.1074−0.00670.1623−0.0054—0.79630.46820.95760.28600.9546Intercept0.25180.20810.22630.23520.23270.24710.53800.61000.57900.56410.56890.5455μ123.880923.874223.881823.886623.879323.8795<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001μ231.266231.259331.264731.271631.266731.2640<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s11.92261.92030.91881.91691.92241.9238<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001s21.46531.46561.46241.46341.46651.4655<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ρ120.19700.19680.19650.19680.19700.1967<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001ν114.9900110.699098.4523103.0250110.6760106.21800.59780.58400.53930.58790.58160.5676


## 4. Discussion and Extensions

Although genetic information can reliably inform about the future appearance of certain diseases and it is an element of great interest for different stakeholders, this genomic information is considered in the highest level of confidentiality and protected from disclosure. In this sense, in the insurance industry, a particularity of genomic information in this context is that it does not only provide information and knowledge about the individual taking the insurance but also in respect to their ancestors and descendants. As a consequence of these limitations, international regulators, e.g., the Council of Europe encourages insurers to update their actuarial bases according to relevant and new scientific knowledge and this may open the gates to explore new avenues and data types and information sources. As part of this new vision, we have examined the potential use of microbiome information in some variables associated with the insurance underwriting. Recent investigations have shown that changes in the gut microbiome are associated to certain risk of pathologies could be a potential proximal predictor of disease onset.Recently, in an unpublished work by using text mining techniques in life insurance literature and microbiome research, a significant overlap between certain diseases and health conditions and other elements that are considered in insurance underwriting. One of these elements is the body mass index (BMI). This is one of the variables considered in the standard health declaration. Traditionally, this declaration is the first step in the risk assessment in health insurance underwriting practice. Certainly, depending on the level of insured capital and age of the policyholder, extra medical examination will be the obligatory required guarantee regardless of the outcome of the standard health declaration. However, medical examinations are expensive, disturbing for the applicant and time-consuming in the underwriting process (see [32]). The importance of BMI is linked to obesity that will lead to large number of chronic diseases, and consequently increase health expenditures and claims costs. Therefore, an early detection of obesity is crucial to safeguard the financial structure of the health insurance provider. Similar conclusions can be drawn about the early detection of cardiovascular, mental metabolic or immune diseases. Then, it is extremely important for the private health insurers to monitor their policyholders’ health status in order to reduce future claims costs.The main findings of our analysis show that the second bacterial coabundance group associated toRuminococcaceae-Bifidobacteriaceae has a significant negative effect on the expected value of the response variable Age for the nonsedentary individuals. This issue is verified not only for the marginal model associated with the response variable age but also for the joint regression model. Concerning the fifth coabundance group related to Pasteurellaceae, it was observed a positive impact on the expected value of age for the sedentary individuals in the marginal model; in contrast, for some predictors, a negative impact in the mean of the response variable age is noticed under the bivariate regression model. This fact is somehow consistent with the recent work of Jollet et al. [33], where it is described a randomized clinical trial that shows in two of the bacterial groups of the Pasteurellaceae group are related to the level of physical activity. However, the degree that this conclusion is valid is arguable since, in general it appears that the majority of people have a unique baseline microbiome that is influenced by environmental conditions. The standard microbiome composition is affected not only by the level of physical activity but also for other factors such as age, diet, and medication. In addition, we have only analyzed a single observation of the gut microbiome per individual limited to the age interval 18–65. It should be highly recommended to perform a longitudinal analysis to relate microbiome information to risk mortality factors and probability of developing certain pathologies. For example, to deal with the high proportion of zeroes in the operational taxonomic units observed, a two part zero-inflated regression model with random effects could be used [26]. In this regard, the American Gut Project (a citizen science project containing more than 10,000 samples not only from the USA but also from several other countries around the world including Australia, UK, or Spain) represents an opportunity to access microbiome data for a variety of age group. This source of information together with the Human Mortality Database available in https://www.mortality.org could be used to analyze how changes in the gut microbiome are related to human longevity in different countries.



---

*Source: 8033806-2022-06-06.xml*