Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
65,277 characters · 15 sections · 92 citation commands
Bayesian shrinkage in mixture of experts models: Identifying robust determinants of class membership
MSC (2010): 62F15, 62J07, 62H30, 90-08\\ Keywords: mixture of experts, classification, shrinkage, bayesian inference, normal gamma prior
Modeling heterogeneity in datasets is a common problem in applied statistics. The task is to find underlying clusters of similar observations that can be used to describe the data. A widespread and known method to accomplish this is finite mixture modeling, where the main idea is to model a single probability distribution as the weighted sum of a finite number of mixture densities. This technique can be used for model based clustering as well as density estimation. Finite mixtures are widely used in different research fields -- a rather common application in marketing research is discussed in lenk2000bayesian, who employ mixture modeling techniques to find clusters of customers with similar behaviour. Earlier references discussing marketing applications are allenby1995 and rossi1996. lubrano2016 use mixtures to find homogenous groups in a study of the income distribution of the United Kingdom. However, the model family also extends to time series analysis naturally as shown in frohwirth2008model. An applied example is the Markov mixture model that fruhwirth2012labor use to model the earning dynamics in the Austrian labour market. For a comprehensive overview of mixture models and estimation strategies, see fruhwirth2006.\\ The main contribution of this article lies, however, in a popular extension of the standard mixture framework. In the most basic Bayesian mixture models, prior class membership is modeled using the component weights, that is the relative size of the mixture clusters. Essentially, this means that the highest prior membership probability is assigned to the largest group in the population. This assumption implicitly claims that each observation has the same prior probability of belonging to a specific group, neglecting other observable characteristics of the data point. To make use of additional information, it is also possible to model mixture parameters as a function of external covariates. Such a specification usually allows for a richer interpretation of the model output and might permit a more holistic use of datasets. This modeling technique is usually referred to as a mixture-of-experts (MOE). Despite the name originating in the machine learning literature\footnote{jacobs1991adaptive calls the mixture components experts and considers the mixture weights as gating networks resulting in the now widely used nomenclature.}, mixture-of-experts models have a wide range of applications, similar to standard mixture models. gormley2008mixture develop a MOE model for rank data and gormley2010mixture use the framework to model network data. MOE models also apply to time series (huerta2003time; fruhwirth2012labor) and longitudinal data (tang2016mixture). Related models have been discussed under different labels for quite some time now, for instance switching regression models (quandt1972new) or concomitant variable latent-class models (dayton1988concomitant). For a comprehensive overview of mixture-of-experts models, refer to gor-fru:mod.\\ This article focuses on a specific problem arising when dealing with mixture-of-experts models where covariates are included to model class membership\footnote{More general mixture-of-experts models also allow for variables to be included to model the the mixture component parameters, see gor-fru:mod.}. There is severe model uncertainty regarding the relevant covariates to include to model prior class membership (as pointed out by anderson2016). Both estimated coefficients and class membership estimates might be sensitive to the particular choice of explanatory variables included in the cluster membership part of the model. One way to resolve this issue is to rerun the model using cross validation as a crude sensitivity analysis. However, the process of choosing which variables to include remains arbitrary. Thus, various approaches for variable selection in mixture frameworks have been discussed in literature, both tackling the questions of which variables to include in the class membership part of the model and which covariates enter the component parameter part of the model. See for instance the generalized smooth finite mixture model from VILLANI2012121, linear cluster-weighted models (ingrassia2014model; generalized in ingrassia2015generalized) or the models for high dimensional mixture regressions in gupta2007variable and devijver2015finite.\\ We propose the use of a continuous shrinkage prior in latent class mixture modeling to isolate robust determinants of class membership to overcome this problem. More specifically, we specify a Normal-Gamma prior (griffin2010inference) and use the P\'{o}lya-Gamma sampler from polson2013bayesian for computations to isolate important predictors of class membership. This results in more efficient shrinkage and improved performance when compared to related methods like standard stochastic search variable priors (george1993variable; introduced to latent class models in ghosh2011), especially when the number of possible predictors is large and/or the sample size is small. We illustrate our approach through simulation studies and an empirical example using Demographics and Health Survey (DHS) data from Mozambique.\\ The remainder of this work is organized as follows: Section 2 introduces the modeling framework. In Section 3, a simulation study is conducted to evaluate the performance of the proposed prior setup. Section 4 illustrates the framework in an application to HIV information source data from Mozambique. Section 5 concludes.
Let $y_i$ denote an observation of data point $i=1,\ldots,N$. This dependent variable can be univariate or multivariate, discrete or continuous, or of a more general structure such as time series or network data. Let $x_i$ be a set of $P$ (for $p=1,\ldots,P$) covariates of $y_i$. Assume $K$ (for $k=1,\ldots,K$) clusters exist in the population that follow known probability density functions $f_k(\cdot|\Xi_k)$ with component specific parameters $\Xi_k$. Denote the component weights as $\eta_k(x_i)$ where $\eta_k(x_i) \geqslant 0$ and $\sum_{k=1}^K \eta_k(x_i) = 1$. Then $y_i$ follows the mixture distribution
We assume the component weights $\eta_k(x_i)$ to be a function of the concomitant variables $x_i$. These covariates influence the distribution of $y_i$ indirectly via the individual prior class membership probabilities $Pr(S_i = k | x_i) = \eta_k(x_i)$. $S_i$ is the latent class membership indicator of individual $i$, where $S_i = k$ if $y_i$ belongs to cluster $k$. This relationship is typically modeled via a multinomial logit link with
where we set $\beta_K = 0$ to achieve identification of the model. This directly results in the interpretation of the coefficients in terms of a change in the log odds relative to the baseline category $K$. Other possibilities to model this ”gating function” are discussed in yuksel2012twenty.
There are several ways to model Bayesian multinomial logistic regression. We choose the method proposed by polson2013bayesian for simplicity and efficiency reasons. The Bayesian framework requires the specification of a prior on $\beta_k$. As we are interested in implicit variable selection (i.e. shrinking coefficients of unpromising explanatory variables to zero), we implement a modified version of the normal gamma prior, a global local shrinkage prior introduced in griffin2010inference:
where $\tau^2_{k,p}$ denotes the local shrinkage parameter of coefficient $p$ in regression $k$. As opposed to griffin2010inference, who apply this prior to a standard regression model, we have to deal with $K-1$ separate sets of coefficients in the multinomial logit framework. Thus, we do not use a single global shrinkage parameter $\lambda$, but introduce global shrinkage parameters $\lambda_k$ per equation. This allows for more flexibility and allows for conducting variable selection for each group of the multinomial logit separately. This might be sensible, taking into consideration that the relevant variables responsible for accurately describing class membership might well alter between classes. Similar prior structures have been implemented into Bayesian time series analysis (see for example huber2017adaptive; bitto2016achieving or kastner2016sparse) and high dimensional spatial models (pfarrhofer2019) recently. To complete the prior setup, we specify the hierarchical structure for $\lambda_k$ and $\tau^2_{k,p}$ to be
The priors on the component parameters $\Xi_k$ are application specific. The choice of values for the hyperparameters $c_0$, $c_1$ and $\theta$ is discussed in appendix (ref).
We implement a Gibbs Sampler to sample the parameters from the full conditional posterior distributions using Markov Chain Monte Carlo (MCMC) methods (robert2013monte). The posterior of the latent class membership indicator $S_i$ is drawn from a multinomial distribution $M(1; p_{i,1},\ldots,p_{i,K})$ with success probabilities $(p_{i,1},\ldots,p_{i,K})$ where
where $f_k$ denotes the probability density function of the components of the mixture distribution. A posteriori, the regression coefficients are normally distributed with
The parameters $m_k$ and $V_k$ of this normal distribution can be derived using the following identities:
where $\mathds{1}(\cdot)$ denotes the indicator function.\\ Let $\omega_{i,k}$ be a latent auxiliary variable that is conditionally P\'{o}lya-Gamma\footnote{For further details on the P\'{o}lya-Gamma sampler as well as standard hyperparameter values and the multinomial setup used here, see polson2013bayesian and their technical supplement.} distributed with
Using this auxiliary variable, the posterior parameters for sampling $\beta_k$ can be derived\footnote{For details on the derivation see the technical supplement of polson2013bayesian as well as for instance koop2003bayesian.} as
where $\Omega_k = \textrm{diag}(\omega_{i,k})$.\\ Finally, the posterior distributions of $\lambda_k$ and $\tau^2_{k,p}$ are both of well-known form and can be derived as
where $P$ is the number of covariates entering the model and $GIG$ denotes the Generalized Inverse Gaussian distribution. hormann2014generating provide an efficient adaptive rejection sampling algorithm that makes it possible to easily draw from the $GIG$. This algorithm is implemented in the R package GIGrvg (leydold2015gigrvg) which we use in our computations. This completes the simulation setup.
Selecting the number of mixture components still remains a challenging issue. Proposed solutions are the use of reversible jump MCMC algorithms (green1995reversible) or shrinkage on the component weights (malsiner2016model). A further commonly used approach is to estimate the marginal likelihoods of models with different number of components and use these likelihoods to decide how many components are suitable (modelselection).\\ Estimating the marginal likelihood is a non-trivial integration problem that involves a number of possible numerical and computational issues. Starting with the purely statistical problem, we need to compute the marginal likelihood given by
where $M_G$ denotes the model with $G$ components\footnote{The notation differentiates between $G$ and $K$ in this subsection to make it clear that $K$ refers to the number of clusters in the data generating process.} and $\Theta_G = (\Xi_1,\ldots,\Xi_G,\beta_1,\ldots,\beta_{G-1})$ denotes the set of all unknown model parameters in a model with $G$ components\footnote{Note that we assume $\beta_G=0$ to achieve identification in the multinomial logistic framework. Thus, only $G-1$ $\beta$ parameters have to be estimated.}. In the overwhelming majority of cases, this integral does not have a closed form solution. However, several methods may be employed to estimate the value of this integral. We use random permutation bridge sampling to estimate the marginal likelihood for model selection purposes. Bridge sampling was first introduced by meng1996simulating and has been thoroughly described for Markov switching and mixture models by fruhwirth2004estimating, who concludes that the bridge sampling estimator is the preferable estimator for the marginal likelihood of this model class and superior to related approaches like importance sampling (geweke1989bayesian) or the harmonic mean estimator (newton1994approximate).\\ To get an estimate of the marginal likelihood, we first need to construct an importance density $q(\Theta_G)$ and generate $L$ i.i.d. draws from this density, denoted by $\Theta^{(l)}_G$ with $l = 1,\ldots,L$. This importance density should have the same domain as the posterior distribution and closely resemble the posterior distribution (GRONAU201780). As shown by fruhwirth2006, the bridge sampling estimator can then be derived as
where $\Theta^{(m)}_G$ with $m=1,\ldots,M$ are the $M$ posterior draws from the Gibbs sampler output using $G$ components and $p^{\star}(\cdot)$ denotes the non-normalized posterior distribution. The choice of $\alpha(\Theta_G)$ is arbitrary, however, meng1996simulating discuss an asymptotically optimal choice which minimizes the expected relative error of the estimator. It is given by
The bridge sampling estimate of the marginal likelihood $\hat{p}_{BS}$ can be obtained using the following algorithm:
In general, both the construction and the evaluation of the importance density for mixture-of-experts models are challenging due to the multimodal nature of the likelihood function. We follow the approach described in fruhwirth2004estimating, who states that the importance density for mixture models can be constructed in a fully automatic manner by saving the posterior distribution moments of $S$ randomly selected MCMC draws. In a random permutation sampler, this results in a multimodal importance density that approximates the modes of the posterior distribution. An i.i.d. sample from this importance density can then be generated by drawing from a uniform mixture of the $S$ saved densities. Additional details on the construction of an importance density and the implementation of a bridge sampler for the proposed model are provided in appendices (ref) and (ref).\\ To compute a marginal likelihood estimate using this procedure, it is necessary to choose a starting value for the bridge sampler. Reasonable choices include alternative estimates of the marginal likelihood. fruhwirth2006 suggests using the importance sampling estimator or the reciprocal importance sampling estimator of the marginal likelihood. Both estimators can be derived from the same functional values that are needed to compute the bridge sampling estimate\footnote{However, other starting values are possible. GRONAU201780 choose $0$ as starting value, stating that ”usually the exact choice of the initial value does not seem to influence the convergence of the bridge sampler much.”.}.\\ Note that in order to evaluate the non-normalized posterior distribution, it is necessary to use the marginal prior densities of the parameters that are specified using a hierarchical prior setup. The marginal prior of $\beta_{k,p}$ is available in closed form and can be derived as
where $K_x(\cdot)$ is the modified Bessel function of the second kind with index $x$ and $\Gamma(\cdot)$ is the gamma function (see for instance bitto2016achieving).\\ A thorough discussion of the bridge sampling technique is out of scope of this article. However, so far literature has been rather sparse on the practical computation of bridge sampling estimates in the context of mixture models and especially mixture-of-experts models. An exception is the recent review by gor-fru:mod who give details on the procedure for mixture-of-experts models.
Parameter estimation in this model family poses various difficulties, especially in a Bayesian framework. Label switching is a known issue when estimating mixture models (hurn2003estimating; jasra2005markov). It is the result of the multimodal likelihood function being invariant to relabeling the components as pointed out by redner1984mixture. This can be problematic as switching labels during MCMC sampling might result in heavily distorted, multimodal posterior distributions that are difficult to summarize. Deriving point estimates such as posterior means then becomes inappropriate (stephens2000dealing).\\ Early approaches deal with label switching by introducing simple restrictions on the mixture parameters such as $\eta_1 < \ldots < \eta_K$ (see for instance lenk2000bayesian). However, identifying simple restrictions in high-dimensional models might be cumbersome or infeasible. In addition, if the restriction does not result in the MCMC sampler visiting all modes of the multimodal likelihood evenly, estimates of the marginal likelihood of the model might be biased according to fruhwirth2004estimating.\\ Early references for other relabeling algorithms include celeux1996stochastic. However, their suggestions require known true parameter values, which makes them difficult to apply in real data settings. stephens2000bayesian suggests to relabel the draws such that the marginal parameter posterior distributions are as unimodal as possible. stephens2000dealing provides a literature review as well as a decision theoretic framework to deal with label switching.\\ We employ the approach described in fruhwirth2006 and identify the posterior draws using a postprocessing procedure via $k$-means clustering. In addition, to force the sampler to explore the full mixture posterior distribution, random permutation sampling is introduced (fruhwirth2001markov). That is, every MCMC iteration is concluded by a random permutation step before storing the parameter draws to achieve balanced label switching.\\ The identification algorithm employed is based on the idea of clustering the parameter draws using distance based measures in the point process representation of the MCMC output. After $M$ saved unconstrained MCMC iterations, $k$-means clustering is applied to all $MK$ posterior draws within a suitable parameter subset. The idea is that draws belonging to the same mixture component will be sorted into the same group by the clustering algorithm. The permutation sequence that results from this $k$-means procedure can then be used to reorder the posterior draws and obtain unique identification for further parameter inference. More formally, we use the following two block algorithm:
Step 2(d) is implemented to ensure that we only use draws where a unique labeling can be found. By removing draws where $\rho^{(t)}$ is not a permutation of $(1,\ldots,K)$, we remove draws where clusters are overlapping in the point process representation. When two or more clusters overlap, no unique labeling for each of the $K$ parameter draws in MCMC draw $m$ is achievable through $k$-means centroid analysis. The ratio of removed draws to the number of saved MCMC draws can be used as an indicator for how well the model is able to separate the mixture clusters. A high rate of non-permutations usually points in the direction of an over-fitting model. For further information on this identification algorithm, refer for instance to malsiner2016model. For further and more specific information on the identifiability of mixture-of-experts models, see for instance jiang1999identifiability or the excellent discussion with many examples in gor-fru:mod.
To illustrate the performance of the proposed prior structure, we conduct a number of simulation studies to compare our approach to other possible model setups. The normal gamma shrinkage prior is compared to the standard prior setup suggested in polson2013bayesian and the stochastic search variable selection prior (SSVS; george1993variable). The basic concept of the SSVS prior is similar to ours in terms of model structure and computational approach. Therefore, a simulation based comparison of the two models seems advisable. The SSVS prior relies on the idea of specifying a mixture of two normal densities as prior for each multinomial logit coefficient. Both normal densities are centered at 0. One has a large variance (”slab”) while the other one has a small variance (”spike”). Using standard mixture modeling techniques, it is possible to estimate whether a particular coefficient will be drawn from the slab or from the spike component of the mixture. Formally, we specify
where $\zeta_1^2 << \zeta_2^2$ and $\delta_{k,p}$ is the binary inclusion indicator of covariate $p$ in group $k$. For details, see george1993variable. Following ghosh2011, we set $\gamma_2^2 = 1$. The normal spike component is specified with variance $\gamma_1^2 = 0.01$.\\ A variety of simulation exercises is conducted. The first study evaluates the performance of the prior only. That is, the relative performance of the NG prior is explored in a multinomial logistic regression setup. In a second step, the three priors are compared in various classification problems where they are employed to cluster observations arising from bivariate normal distributions. Overall, the simulation studies suggest a rather similar performance of the SSVS prior and the NG prior when it comes to estimating the coefficients in the class membership part of the model. However, the NG prior usually shows slight benefits, especially in shrinking unnecessary coefficients to zero, in high sparsity settings and when estimating marginal likelihoods. Details are discussed below.
In this subsection, the prior of polson2013bayesian applying no shrinkage (hereafter ”Standard Prior”\footnote{We set the prior variance of the standard prior to 10 as proposed in polson2013bayesian.}) and the SSVS prior are compared to the NG prior in a multinomial logistic regression setup. This preliminary analysis allows us to evaluate the shrinkage performance independently of the mixture setup.\\ Using the data generating process in Eq. (ref), we simulate four groups with 750 observations and 20 explanatory variables each. The true parameter vectors are chosen to be sparse, thus creating the need for considerable shrinkage within the estimation of the multinomial logistic regression. The true coefficient values are $\beta_1 = (0.8, 1, 2, 0.5, 0, \ldots)'$, $\beta_2 = (0.3, 0, 0, 0, -1, 1.7, -2, 0, \ldots)'$ and $\beta_3 = (0.3,1,-2,0.8,0.9, 0, \ldots)'$\footnote{This setup corresponds to the simulation study conducted in ghosh2011}. All explanatory variables are drawn from a standard normal distribution. Note that this setup implies that we need to deal with group specific relevant membership predictors. Thus, to obtain good estimates, group specific shrinkage is necessary. As this simulation uses a large number of observations, a quite informative likelihood results. Hence, we extend the setup described above by two scenarios using 300 and 100 observations, respectively. This should enable us to evaluate the prior performance in an environment with comparatively uninformative data. We implement a Gibbs sampler using 25000 draws after a burn-in period of 5000 draws. The mean estimates of 25 simulation runs are then compared.\\
Table (ref) reports the root mean squared error (RMSE) with respect to the coefficients that are truly zero, the coefficients that are truly different from zero, all coefficients and the predicted probabilities (P.P., defined in Eq. (ref)) resulting from the estimation. This enables us to separately examine how well the priors are able to shrink unimportant coefficients to 0, how precise the point estimates are and whether they are able to give useful estimates of the predicted probabilities. These predicted probabilities are of utter importance in the mixture-of-experts framework, as they will directly influence class membership and therefore all estimated model parameters.\\ The results suggest that the first simulation using 3000 observations is not a very competitive environment. The likelihood is quite informative, resulting in precise estimates even for the standard setup without introducing shrinkage. Figure (ref) plots the true values against the posterior mean estimates of the respective models. The uncertainty surrounding the posterior means is given by an interval of $\pm 1.96 * SD$ where $SD$ is the posterior standard deviation. Scatterplots suggest that all three models are able to revocer the true coefficient values well. Nevertheless, the NG setup performs particularly well and even outperforms the SSVS setup in terms of precision. However, it comes at the cost of a slightly prolonged computation time.\\ Using just 10% of the observations, estimation becomes more difficult as the data becomes less informative as seen in Figure (ref). The point estimates become considerably worse. The standard prior has problems to recover the true values, as the enlarged RMSEs indicate. The NG prior shows slight advantages in terms of shrinkage and in predicting cluster membership probabilities. However, the SSVS setup is able to provide more accurate point estimates and therefore has a slightly lower RMSE with respect to the true non-zero coefficients and regarding the overall coefficient RMSE. Further reducing the number of observations to $N=100$ leads to inflated coefficient estimates and increasing uncertainty when applying no shrinkage, as depicted in Figure (ref). This results in enlarged RMSEs. The performance of the SSVS and NG prior remains rather similar to the case with $N=300$, however, the shrinkage priors also show larger uncertainty surrounding the posterior means. The SSVS prior produces better point estimates, but is not as efficient as the NG prior when it comes to shrinking unnecessary coefficients to zero. The NG prior shows a slightly better performance when predicting the class membership probabilities. We would like to note that both the NG prior and the SSVS prior perform rather well in absolute terms, producing small RMSEs in general\footnote{It should also be noted that the performance of the analyzed priors depends on the imposed prior variances. However, simulations using different variances did not change the results qualitatively. See sfswagner for a thorough discussion and comparison of various shrinkage priors.}. All coefficient estimates of these simulations are provided in Tables 4-6 in Appendix (ref). The performance of the three priors in a full mixture-of-experts classification setting are discussed in the next subsection.\\
\\
To further examine the performance of the priors, four simulation studies in a full mixture-of-experts framework are conducted. These simulations differ from each other with respect to the degree of overlapping of the clusters, the number of regressors and the complexity of the sparsity structure in the true coefficient vectors.\\ Datasets with 300 observations and four clusters arising from bivariate normal distributions are generated. We differentiate between a ”well separated” case and an ”overlapping” case using $\mu_{1} = (-1.5,-0.5)'$, $\mu_{2}=(0,1.3)'$, $\mu_{3}=(1,-1)'$, $\mu_{4}=(3,-2)'$ for the ”well separated” case and $\mu_{1} = (-1.5,0)'$, $\mu_{2}=(0.5,0.5)'$, $\mu_{3}=(1,-0.5)'$, $\mu_{4}=(3,-0.5)'$ for the ”overlapping” case, respectively. The variance covariance matrices $\Sigma_k$ are chosen to be $0.25I$ for the separated case and $0.2I$ for the overlapping case for all $k=1,\ldots,4$. Figure (ref) depicts two example datasets, representing the the two cases. In addition, we differentiate between a medium number of regressors (corresponding to the same true coefficient vectors as in section (ref)) and a ”high sparsity” case where 60 covariates that are not part of the data generating process are added to the covariate dataset, resulting in a total of 80 predictors in the model. Finally, we look at a case with a more complex sparsity structure as compared to section (ref). In this scenario, the first three predictors are only relevant to the first group, the second three predictors are only relevant to the second group and the third set of four predictors is only relevant to the third group. The last predictor is relevant for all groups. This setup requires the shrinkage priors to flexibly vary the amount of shrinkage by group. An overview of the four different simulation setups is given in Table (ref).\\
For every setup, various summary statistics are computed. As before, RMSEs with respect to zero and non-zero coefficients as well as overall RMSEs and RMSEs with respect to the predicted probabilities are reported for the models with $K=4$. In addition, the misclassification rate of each model is computed. To assess the ability of the models to recover the true number of clusters, we run each simulation study for $K=(2,\ldots,6)$ and report plots of the average log Bayes factors\footnote{Where all models are assigned equal probabilities a priori.} with respect to the model with the true number of clusters $K=4$. Gibbs samplers using 2000 draws after a burn-in period of 5000 draws are implemented. Again, simulations are repeated 25 times and the resulting means of computed statistics across simulations are reported.
Table (ref) reports the main simulation study results. On average, we find that the performance of the SSVS and NG prior is very similar in all simulation exercises. However, a slightly better performance of the NG prior can be found in some cases, especially in the high sparsity setting. Nevertheless, as mentioned before, RMSEs are small for SSVS and NG in absolute terms, suggesting both priors are in principal useful when selecting covariates in a mixture-of-experts framework. We do not report the estimation results for $\mu_k$ and $\Sigma_k$ as all three models perform extremely well in this regard, showing very similar results. The only issue that stands out is that model with the standard prior has a tendency to overestimate the variances. As pointed out by one of the reviewers, it might be illuminating to look at RMSEs that are separately computed by groups. This does not lead to any significant qualitative or quantitative variation in the results. Thus, RMSEs by group are not reported for brevity reasons. More detailed results are available from the author upon request.\\
Figures (ref) - (ref) plot the average log Bayes factors relative to the true model with $K=4$ for the three priors. Positive values suggest that the respective model scored a higher marginal likelihood than the true model with $K=4$ and vice versa. In most cases, this model selection criterion suggests to choose $K=4$, with two exceptions: First, SSVS seems to have a slight tendency to favor models with a smaller number of clusters as compared to the other models. This leads to positive log Bayes factors for the models with $K=2$ and $K=3$ in the ”Overlapping” simulation setup. Second, in the ”High Sparsity” scenario, all log Bayes factors lean towards models with $K=2$ and $K=3$, suggesting influence of the number of predictors on the bridge sampling estimates. However, an in-depth examination of this issue is out of scope of this article and thus left for future research.\\ All in all, we conclude that the NG prior is a very useful alternative to SSVS in mixture-of-experts frameworks. Generally speaking, the results of the SSVS prior and the NG prior will be very similar, although there are some performance gains of the NG prior visible in terms of shrinkage as well as with respect to model selection issues.
Mozambique is a country in Southeastern Africa that is considered one of the poorest and most underdeveloped countries in the world, scoring low in both economic and human development rankings. In the year 2008, Mozambique had the 8th highest HIV prevalence in the world with 1,600,000 people infected (11.6% of the population) of whom around 990,000 were women and children. According to the Joint United Nations Programme on HIV/AIDS, there are around 590,000 HIV orphans living in Mozambique, 180,000 of whom are infected with the virus themselves, a large part due to mother-child transmission. 75% of the infected population between the age of 15 and 19 is female. Moreover, a large gender disparity regarding the level of information on the disease can be observed. While around half of the male adolescent population has comprehensive knowledge on HIV, only 27.4% of adolescent women have enough information to adjust their behaviour to protect themselves and their children according to the United Nations Children’s Fund. This disparity is suspected to be largely due to socioeconomic and sociocultural reasons, with the main drivers being traditional gender roles and religious involvement (AGADJANIAN20051529).\\ Consequently, it is crucial to isolate channels that can be used by the government and non-governmental organizations to disseminate vital information on HIV, especially to the female population. Informing females about HIV has proven not only to decrease the infection rate but also increase the economic and social independence of women (audet2010sociocultural). Our empirical example contributes to this relevant and important issue by clustering women in Mozambique into groups that are relatively homogenous with respect to their information sources on HIV, similar to dias2010modeling. In addition, we use a large dataset of potential geographic and socioeconomic explanatory variables and isolate the most important factors that determine membership in those information clusters. The results may be used to derive for instance information campaign strategies for respective subgroups.
We use a set of binary variables that indicates whether a particular woman uses a specific source to gather information on HIV or not. A convenient choice of mixture distribution is the Bernoulli distribution, which proves useful when clustering binary vectors (see for example the vast literature on market segmentation; wedel2012market).\\ Let $y_i = (y_{i,1},\ldots,y_{i,J})$ be a $J$-dimensional vector of 0s and 1s that describe the HIV information sources used by woman $i = 1,\ldots,N$. Assume that this vector is the realization of a binary multivariate random variable $Y = (Y_1,\ldots,Y_J)$. Now suppose there exist $K$ groups in the population that cause differences in occurence probabilities $\gamma_{k,j} = \text{Pr}(Y_j = 1 |S_i = k)$ in $K$ different classes for $J$ different binary variables. $S_i$ is the latent class indicator of woman $i$. We can rewrite Eq. (ref) where $y_i$ follows the mixture distribution
The $K$ components correspond to the latent classes in the population. This model is widely used in various research fields, starting as early as lazarsfeld1959latent. For details, see fruhwirth2006. We assume that all probabilities $\gamma_{k,j}$ are a priori independent and specificy a beta prior of the form
and derive the posterior distribution conditional on the latent class indicators $S_i$, given by
where
We apply the proposed model to data compiled from the Demographics and Health survey (DHS) for Mozambique from 2003. The DHS is a nationally representative household survey on a wide range of topics, including HIV information sources and various socioeconomic, geographic and health related variables.\\ The dataset includes information on 11,922 women. Ten different information sources are used to cluster these women into groups and a set of around 40 external covariates enters the model to explain class membership. These variables cover socioeconomic characteristics like age and education, region of residence, relationship status and sexual behavior as well as poverty related measures and dwelling characteristics. Table (ref) provides a detailed overview of the candidate explanatory variables.
We estimate the model with a NG shrinkage prior for different values of $K$ and compare the resulting models using the marginal likelihood estimates obtained via bridge sampling.\footnote{The model has been implemented in R (R). Computational time for $K=4$ is around 50 minutes for 5000 draws after a burn in period of 1000 draws on an Intel i7 @ 2.4 GHZ.} We choose the model that maximizes the marginal likelihood. The bridge sampling estimates of the log marginal likelihood for $K=2,\ldots,6$ are provided in Figure (ref). The model with $K=4$ scores highest and is therefore discussed below.\footnote{This is of course not the only way to proceed here. Especially in a development context, other, more informal model selection criteria that take into account long term campaigning strategies or financial constraints may be employed. For example, the groups "Modern & Educated" and "TV/Radio" could be merged as they are both have a distinctive dependence on TV. However, in this paper the statistical possibilities of the proposed model are emphasized and hence we make use of the purely statistical approach.}
The estimates for $\gamma_{k,j}$ are presented in Figure (ref). The uncertainty surrounding these estimates is usually extremely small. At first glance we find that the radio as well as friends and relatives seem to be important information sources for all groups. For the purpose of further interpretation of the model results, we name the groups with respect to their most distinctive HIV information source as described below.
Around 8% of the female population use modern information sources such as television, newspapers and posters. In addition, this group obtains a relatively high amount of information from schools. Thus, we label this group as ”Modern & Educated”. A somewhat larger group (around 13% of the population) relies mostly on TV but is highly unlikely to inform themselves in schools. Hence, we name this group ”TV/Radio”. The third group, ”School”, which comprises around 6% of the female population of Mozambique, relies heavily on schools for obtaining information on HIV. The baseline group ”Community” (73%) has an above average dependence on friends and relatives, community meetings and local churches in terms of information on the disease.\\ Figure (ref) provides a plot of the point estimates of the logit coefficients.
When interpreting the multinomial logit coefficients, one has to keep in mind that the effects are always interpreted with respect to a baseline group. For convenience in the estimation process, we choose the largest group as baseline group (”Community”). In terms of the other categorical variables, the ”baseline woman” is residing in Maputo City, has no religion, is married and is a member of the richest wealth group.\\ Strong effects of the wealth distribution on the probability of being a member of the ”Modern & Educated” and ”TV/Radio” group are observable. These groups also share an above average probability of having access to a flush toilet and electricity. In addition, it is relatively unlikely that a woman lives in the countryside and is a member of one of those groups. These findings are in line with what theory suggests in a poverty plagued country like Mozambique.\\ Unmarried women with above average education are also likely to rely on schools as HIV information sources. This seems puzzling, as female education is a primary development issue in various African countries. However, one should keep in mind that this group is extremely small and comprises just above 6% of the female population. Interestingly, wealth variables seem to be not strongly correlated with the probability of being a member of the ”School” group as compared to the other groups. However, we see that geographic variables determine prior class membership for this group, implying that we can find spatially clustered communities with above average female educational attainment throughout specific provinces.\\ These results are particularly relevant to policy makers. Important insights that can be derived are, for instance, that HIV information campaigns that are targeted on disseminating educational materials via schools are likely to be more effective in Maputo City as compared to Sofala. It also might be a good idea to target folders that are distributed in schools towards single women as opposed to married women. However, these are mere examples. A detailed discussion of the policy implications of the results is out of scope of this article.\\
Finite mixture models are a commonly used tool for model based clustering and density estimation. They can be extended to mixture-of-experts models, allowing to use information from several covariates when clustering dependent variables of arbitrary form. We propose the usage of continous shrinkage priors to find robust predictors of class membership in this context. This enables us to simultaneously identify underyling groups in a population, cluster said observations into these groups and find the important predictors of being a group member. In particular, we suggest a combination of the normal gamma prior (griffin2010inference) and the P\'{o}lya-Gamma sampler (polson2013bayesian) for implicit variable selection in a multinomial logistic regression that is used to model prior class membership.\\ This setup solves the issue of model uncertainty that arises in this context and reduces the sensitivity of the model with respect to included variables. The proposed framework slightly outperforms related approaches and makes more precise clustering in setups with a large number of predictor variables possible.\\ We illustrate the model in a real data application where we apply the model with a mixture of Bernoulli distributions to HIV information sources of women in Mozambique. Model selection is based on the bridge sampling estimate of the marginal likelihood. We find four clusters of women who are relatively homogenous with respect to their HIV information sources. Somewhat unsurprising, we find that wealth plays an important role in the access to information on HIV. Moreover, geographical patterns of information seeking behavior seem to be prevalent.\\ Further research may be pointed into the direction of comparing the performance of different shrinkage priors in this context in a more detailed way as seen in sfswagner. One promising candidate is for example the Dirichlet-Laplace prior from bhattacharya2015dirichlet. It might also be possible to extend various other Bayesian variable selection methods to mixture-of-experts frameworks, for example Bayesian compression (guhaniyogi2015bayesian) or its extension using targeted random projections (mukhopadhyay2017targeted). Another interesting problem is how to apply the idea of shrinkage introduced through the prior class membership weights (e.g. malsiner2016model) for model selection purposes into a mixture-of-experts framework. Also, the evaluation of the forecasting performance of the model was not attempted in this article and is left for further research.