Modelling ethnic differences in the distribution of insulin resistance via Bayesian nonparametric processes: an application to the SABRE cohort study

2.1 Prior distributions

The model in Eq. (2) requires the specification of a prior distribution for the vector of regression coefficients β=(β1,…,βp−1), the intercept βi0 and the precision parameter τi2. We adopt a nonparametric prior, the DGDP prior, on the vector (βi0,τi2). As explained below, this choice of prior distribution allows to cluster the observations. Moreover the use of the DGDP prior provides both flexibility and parsimony about the number of parameters that we introduce in the model. We now present a brief review of the main properties of the Dirichlet Process (DP) and its generalisation to the DGDP. The DP (see Ref. [9], [10] for basic properties) is arguably the most widely used nonparametric Bayesian prior, mainly because of computational simplicity and ease of interpretation. In DP based models computational complexity of posterior simulation is in theory dimension independent, allowing specification of possibly high dimensional random effects distributions. A random measure P that is generated by a DP is almost surely discrete. Sethuraman [11] provides a constructive definition of this process, showing that if a random probability measure P is distributed according to a DP, with mass parameter α and base measure G0, then it admits the following representation:

where the atoms θ1,θ2,… are iid realisations from G0, δθk is the Dirac measure that assigns mass probability one in correspondence of the location θk. The weights ψk follow a stick-breaking process (see Ref. [12] for details):

with ϕk∼iidBeta(1,α) and ψ1=ϕ1. By construction we have 0≤ψk≤1 and ∑k=1∞ψk=1. The discreteness of the DP induces clustering of the subjects in the sample based on the unique values of the random effects parameters (in our case θk=(β0k,τk2)), where the number K of clusters is unknown and learned from the data.

In this paper we are interested in modelling the distribution of HOMA IR in each of the three ethnic groups (i) allowing for borrowing information across groups (ii) highlighting differences and similarities (iii) accounting for the effect of covariates. To this end, we employ a generalisation of the DP proposed by Ishwaran and James [12] and Hjort [13]; the Generalised Dirichlet Process (GDP). The GDP employees a richer parametrisation in the stick-breaking construction, allowing greater flexibility in the moments of the random distributions. Consider the stick-breaking in Eq. (5), where the elements ϕk are draws from a Beta (1, α). In the generalisation proposed by Hjort [13], the ϕk still draws from a Beta distribution, but the first hyper-parameter does not need to be fixed to one. In what follows we use an alternative parametrisation of the Beta distribution, where the hyper-parameters are specified in terms of the mean and the concentration parameter. In the GDP the {ϕk} are then independent draws from a Beta(μkυk,(1−μk)υk):

where

and

with υk∈(0,∞), are the expected value and the variance of the Beta random variable respectively. The weights of the GDP admit the same stick-breaking construction as for the DP. Hjort [13] proposes a more parsimonious parametrisation of Eq. (6), setting μk=μ and υk=υ. This simplification does not impose significant restriction in applications. We now explain how we introduce ethnicity information in the distribution of HOMA IR. The final model will contain two main components: one for the clinical covariates and one for the patients effect. The model for the covariates expresses prior information on how covariate influence the clinical outcome, while the nonparametric prior (GDP) is used as random effect distribution to capture inter-patients variability. Moreover, it is desirable to specify a random effect distribution for each ethnicity in a way that the random effect distributions are related (similar or very different), but not necessarily identical.

There is a wealth of literature on how to extend the DP to incorporate covariate information, for example, letting the weights and/or locations of the infinite mixture in Eq. (4) depend on a variable of interest that defines sub-groups in the observations. See the seminal paper of MacEachern [14] on the Dependent Dirichlet Process (DDP). Similarly, also the GDP can be extended in presence of categorical covariates. Barcella et al. [6] introduces the DGDP, where the dependence among random distributions is introduced through the weights ψk of the mixture in four. The parameters ψk are generated from the stick-breaking process, so the dependence is introduced directly on the parameters υ and μ.

Consider G groups defined by a covariate of interest g∈G, where G is the covariate space. We let μ, which represents the mean of the Beta random variables depend on the particular value of g, while we assume the same υ across groups. We denote with μg the mean of the Beta random variable corresponding to group g. In our application groups are defined by the ethnicity. The random measure Pg, i.e. the random distribution associated to group g, is then defined as:

Here θk=(β0k,τk2). In particular, dependence across the μg,g=1,…,G, is obtained by specifying a Beta regression on μg and using a categorical predictor zg, i.e. an indicator variable which denotes to which group the observations are associated to. This strategy allows for group dependent clustering of the observations. See Ref. [6] for details and clustering properties. Finally, the full model for HOMA IR is specified as follows

for g=1,…,G. The parameter μ is linked through a function f (e.g. logit or probit), mapping from (0,1) into (−∞,∞), to the linear predictor zgη where η is a vector of regression coefficients of appropriate dimension to which we assign a standard Normal prior:

where ημ and ηΣ denote the prior mean and covariance matrix respectively. We assume independence a priori between the parameters β0 and τ2, which is reflected in the choice of the base measure G0, defined as the product of a Normal distribution and a Gamma distribution.

We specify a Spike and Slab prior on each of the p−1 regression coefficients βj. This prior specification provides an effective variable selection strategy [7], [15]. We introduce indicator variables ωj:

where p(ωj=1|π)=π is the probability of the slab, i.e. the probability that a covariate is included in the model, while 1−π represents the probability of the spike, i.e. the probability that the regression coefficient corresponding to the jth covariate is equal to 0 and does not affect the response. As before, δ0(βj) is a mass point at zero, representing the spike of the mixture. μβ represents the prior mean (usually set to 0) of the slab component and τβ2 is the prior precision.

The parameter π is assigned a Beta prior:

where πa and πb are the hyper-parameters of the Beta distribution (e.g. setting πa=πb=1 gives a uniform distribution). Appropriate choices of these hyper-parameter allows us to impose sparsity in the variable selection.

2.2 Posterior inference

Posterior inference is performed through Markov Chain Monte Carlo (MCMC) methods. A detailed description of the algorithm is provided in Supplementary Material. We run the MCMC for 30,000 iterations, discarding a burn-in period of 10,000, thinning every five iterations. We specify the following hyper-parameters: the truncation level of the stick-breaking is set to L=30. This threshold works well in our application since the number of non-empty clusters is always below 20. The base measure parameters are set to m0=0, κ02=0.1, τa=0.5, τb=0.5. The DGDP concentration parameter υ has a Gamma prior with aυ=2 and aυ=1, with expected value aυ/bυ=2. The regression coefficient η is parametrised with prior mean ημ=(0,0,0) and ηΣ=I3. The slab of the regression coefficient βj is a Normal distribution with prior mean μβ=0 and prior precision τβ2=0.1. The prior inclusion probability π has a Beta prior with parameters πa=πb=1.

3.1 HOMA IR: cluster analysis

Posterior inference for HOMA IR shows evidence of 10 clusters. We use the Binder loss function [16] available in the R package mcclust, to estimate the number of clusters in the sample and the clustering allocation based on the MCMC output. In Figure 2 we show the empirical distribution of the outcome HOMA IR in each of the 10 estimated clusters. The overlap between some of the curves is due to the fact that the clusters are estimated conditionally to the covariates and metabolite levels included in the regression model. Table 1 summarises the ethnic composition of each cluster, while Table 2 provides some basic information in terms of age, smoking habits, percentage of females and percentage of first generations (i.e. foreign-born) migrants in each cluster. It is worth noting that clusters 8, 9 and 10, the most insulin resistant clusters, are mostly composed of first generation migrants. In Table 3 we report the ethnic composition of each cluster based on the sub-region of origin. The majority of South-Asians come from the Punjabi-Sikh minority, which represents the major South-Asian component in each cluster, with the exception of cluster 9, where there is a higher percentage of South-Asians of Muslim origin. To understand which covariates are the most important determinant of the response, we examine the posterior probability of each regression coefficient to be different from zero, p(ωj=1|rest,data). 10 predictors have the respective p(ωj=1|rest,data)>0.5 and are considered for further analysis (Figures 3–5).

Figure 2:

Empirical distribution of HOMA IR in each cluster. Black lines denote clusters with a higher proportion of Europeans, while red lines denote a higher proportion of South-Asians. A description of cluster main characteristics is given in Table 1. The numbers above each distribution denote the cluster.

Figure 3:

From top to bottom, boxplots of Acetoacetate, Alanine and Glycine. Black boxplots indicate clusters with a majority of Europeans, while red boxplots indicate clusters with a majority of South-Asians.

Figure 4:

From top to bottom, boxplots of Histidine, Isoleucine and Phospholipids in large HDL. Black boxplots indicate clusters with a majority of Europeans, while red boxplots indicate clusters with a majority of South-Asians.

Cluster 1 is the largest and least insulin resistant group (n = 1249, 57% of participants). Its ethnic composition is: 71% of Europeans, 39% of South-Asians and 70% of Africans-Caribbean. The second largest group is cluster 6, comprising 243 participants (11% of the total, of which, 8% of Europeans, 13% of South-Asians and 24% of Africans-Caribbean). It is evident the net distinction between clusters with a South-Asian majority, compared with Europeans, which are all characterised by higher levels of HOMA IR, with the exception of cluster 4. Cluster 4 is entirely composed of South-Asians. In particular the cluster is characterized by almost only first-generation migrants, a higher proportion of females (0.33), relative to the other clusters and a lower proportion of current smokers (0.16), compared with the other cluster containing South-Asians. Cluster 5 (entirely Europeans) compared with cluster 1, presents both measures of adiposity(subscapular skinfold and sagittal diameter) are modestly higher, while levels of the amino acids tyrosine and isoleucine are significantly higher. Moreover, acetoacetate levels are lower compared to cluster 1, while the levels of alanine aminotransferase (ALT) in cluster 5 are higher than in cluster 1 (Figure 5), the latter suggesting that raised HOMA IR levels may be characterised in this cluster by increased insulin levels with reduced clearance of insulin by the liver [17], implying relatively intact pancreatic beta cell function. The metabolite patterns for cluster 5 also indicate associations with both central and subcutaneous adiposity and amino acid perturbations.

Figure 5:

From top to bottom, boxplots of Tyrosine, sagittal diameter and Alanine Aminotransferase. Black boxplots indicate clusters with a majority of Europeans, while red boxplots indicate clusters with a majority of South-Asians.

Each of the 10 clusters has a distinctive metabolic and phenotypic profile, consistent with suggestions that there are different pathways to type 2 diabetes [18] and that some pathways may be more strongly associated with a particular ethnic group. For example clusters 4 and 9 are entirely composed of South-Asians, while clusters 2, 3 and 5 are entirely Europeans. Of these clusters, 8, 9 and 10 are among the most insulin resistant with high levels of tyrosine, alanine, ALT and subcutaneous adiposity. Some of the clusters identified are very small and will need replication in larger studies together with formal pathway analysis. However, these methods have generated intriguing, novel and persuasive clusters, which highlight the complexity and potential multiplicity of mechanisms underlying the development of insulin resistance and type 2 diabetes.