A Bayesian Nonparametric Approach for Inferring Drug Combination Effects on Mental Health in People with HIV

04/11/2020 ∙ by Wei Jin, et al. ∙ 0

Although combination antiretroviral therapy (ART) is highly effective in suppressing viral load for people with HIV (PWH), many ART agents may exacerbate central nervous system (CNS)-related adverse effects including depression. Therefore, understanding the effects of ART drugs on the CNS function, especially mental health, can help clinicians personalize medicine with less adverse effects for PWH and prevent them from discontinuing their ART to avoid undesirable health outcomes and increased likelihood of HIV transmission. The emergence of electronic health records offers researchers unprecedented access to HIV data including individuals' mental health records, drug prescriptions, and clinical information over time. However, modeling such data is very challenging due to high-dimensionality of the drug combination space, the individual heterogeneity, and sparseness of the observed drug combinations. We develop a Bayesian nonparametric approach to learn drug combination effect on mental health in PWH adjusting for socio-demographic, behavioral, and clinical factors. The proposed method is built upon the subset-tree kernel method that represents drug combinations in a way that synthesizes known regimen structure into a single mathematical representation. It also utilizes a distance-dependent Chinese restaurant process to cluster heterogeneous population while taking into account individuals' treatment histories. We evaluate the proposed approach through simulation studies, and apply the method to a dataset from the Women's Interagency HIV Study, yielding interpretable and promising results. Our method has clinical utility in guiding clinicians to prescribe more informed and effective personalized treatment based on individuals' treatment histories and clinical characteristics.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 5

page 15

page 23

page 39

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Early initiation and adherence to antiretroviral therapy (ART) regimens optimizes health outcomes in people with HIV (PWH) and prevents further HIV transmission (Bangsberg et al., 2001; de Olalla Garcia et al., 2002; Yun et al., 2005; Saag et al., 2018). However, viral rebound is possible due to the high viral evolutionary dynamics and the occurrence of drug-resistant mutations, ultimately resulting in treatment failure. ART agents fall into several classes including nucleotide reverse transcriptase inhibitor (NRTI), non-nucleotide reverse transcriptase inhibitor (NNRTI), protease inhibitor (PI), integrase inhibitor (INSTI), and entry inhibitor (EI). Different drug classes target HIV via different mechanisms. While each individual ART drug is susceptible to certain resistant mutations, a combination of drugs from different drug classes can successfully suppress the virus. Therefore, modern ART regimens typically combine three or more drugs of different classes.

Despite the remarkable success of effective ART reducing disease-related morbidity and mortality (Trickey et al., 2017), numerous observational studies have reported ART-related adverse effects on central nervous system (CNS) function including depression, anxiety, sleep disturbances, suicidal ideation, developmental disorders, and neurological toxicities (Rekha, 2018; Zash et al., 2018). For example, significant neuropsychiatric side effects such as nightmares, hallucinations, and depression (Gaida et al., 2016), have been reported for EFV (efavirenz), an NNRTI. Peripheral neuropathy has been reported for D4T (stavudine), an NRTI, especially when used in combination with other NRTIs (Arenas-Pinto et al., 2016; Saag et al., 2018). These side effects may result in ART discontinuation with downstream consequences such as difficulties in work performance and functioning, HIV disease progression, and increased likelihood of HIV transmission (Fazeli et al., 2015; Watkins and Treisman, 2015). Since ART is recommended for PWH indefinitely, it is critical to understand and quantify drug effects, especially the drug combination effect, on CNS function to facilitate the design and effectiveness of ART regimens.

In this paper, we focus on ART-related effect on depressive symptoms. Depression is one of the leading mental health comorbidities in PWH, affecting from 20% to 60% of those with the virus (Bengtson et al., 2016). Depression is associated with numerous adverse consequences including poor ART adherence (Chattopadhyay et al., 2017), rapid disease progression (Ironson et al., 2017), and increased risk-taking behaviors (Brickman et al., 2017). The high prevalence and the harmful effect of depression among PWH highlights the need for effective clinical management and adequate treatment for depression. To date, few studies are dedicated to investigating the effects of ART regimens on depression, many of which present inconsistent findings. For instance, Pearson et al. (2009) and Mollan et al. (2014) reported increased ART-related depression, whereas Okeke and Wagner (2013) and Jelsma et al. (2005) reported the opposite. One possible explanation is that the effects of ART regimens are heterogeneous and may be confounded by numerous factors such as socioeconomic status, behavioral factors, and clinical performance. ART may alleviate depressive symptoms for some individuals through viral suppression and physical health improvement. For others, ART-related neurotoxicities may aggravate depression and lead to treatment failure. Therefore, investigating the heterogeneous effects of ART on depression among PWH while accounting for major confounders can help identify individual factors that drive depression with ART exposure, thereby facilitating precision medicine for PWH.

Large-scale HIV datasets, such as the Women’s Interagency HIV Study (WIHS), provide an opportunity and challenge to study the effects of ART regimens on depressive symptoms. The WIHS is a prospective, observational, multicenter study which includes women with HIV and women at-risk for HIV infection in the United States (Bacon et al., 2005). Sociodemographics, medication use, clinical diagnoses, and laboratory test results are collected longitudinally with the goal of investigating the impact of HIV infection on multimorbidity. For example, Figure 1(a, b) presents two individuals’ ART medication data versus their clinic visits denoted by calendar dates. They were followed for different time periods with distinct visit dates and drug uses. Their corresponding four depression scores were also recorded at each visit measuring somatic symptoms, negative affect, lack of positive affect, and interpersonal symptoms. The depression scores were summarized from a self-report questionnaire using the Center for Epidemiological Studies Depression Scale where a higher score reflects worse symptoms (CES-D, Radloff 1977; Carleton et al. 2013), as shown in Figure 1(c, d).

(a) ART use for individual 1 (b) ART use for individual 2
(c) Depression score for individual 1 (d) Depression score for individual 2
Figure 1: ART medication data and depression scores of two individuals versus their clinical visits denoted by calendar dates.

The complexity of longitudinal observations, a heterogeneous population, and dynamic and mixed ART assignments present four analytical and modeling challenges. The first challenge is high-dimensionality. With more than 20 ART drugs on the market, there are nearly a half million possible drug combinations, making the estimation of drug combination effect a high-dimensional problem. The second issue is unbalancedness. Some ART combinations are frequent whereas others are rare. For example, D4T+LAM+NFV (two NRTIs + one PI) was recorded 993 times in the WIHS, while a similar ART regimen D4T+LAM+ATZ was only recorded 12 times. The third issue is sparseness. Only a tiny portion (hundreds of drug combinations) of the high-dimensional ART regimen space were observed in the WIHS, while the inference on the entire space is desired. Finally, there is the issue of non-stationarity. Since PWH are given different ART regimens during the course of treatment, the effect of the current regimen on depression is likely dependent on prior ART regimen use.

Statistical learning methods, such as logistic regression, tree-based models, and neural networks have been used to study the effect of ART regimens on survivals and to predict virological responses

(Altmann et al., 2007; Larder et al., 2007; Altmann et al., 2009; Caniglia et al., 2017)

. However, the representations of ART regimens in these models were simplistic, either using a binary variable to indicate whether an individual is on ART, or lumping ART regimens together into a few coarse types. For example,

Lundgren et al. (2002) dichotomizes ART regimens into those with or without a PI. Bogojeska et al. (2010) proposed to predict binary virological responses to ART regimens by fitting a separate logistic regression model for each regimen and borrowing information from similar regimens. The similarity between regimens was defined by a linear additive function. Although computationally efficient, the linear approach treats drugs as exchangeable and does not account for drug classes. Moreover, it cannot be easily extended for modeling ordinal or continuous outcomes, adjusting for covariates, and considering treatment histories in a longitudinal setup.

To address the aforementioned challenges, we develop a novel Bayesian nonparametric model using and extending subset-tree kernels (Collins and Duffy, 2002) and the distance-dependent Chinese restaurant process (Blei and Frazier, 2011; Dahl et al., 2017) to estimate the effects of ART regimens on depressive symptoms after adjusting for relevant socio-demographic, behavioral, and clinical factors. The subset-tree kernel method represents drug combinations in a way that synthesizes known regimen structure and the corresponding drugs into a single mathematical representation to induce an appropriate similarity among different ART regimens. This formulation enables us to efficiently borrow information across ART regimens and develop inferences and predictions on unobserved ART regimens. Furthermore, we extend the subset-tree kernel among ART regimens to sequences of ART regimens in a longitudinal setup, and use the distance-dependent Chinese restaurant process as a prior to capture heterogeneity among individuals by considering individuals’ treatment histories. The R code implementing our model can be found at https://drive.google.com/open?id=1FB8o0cHx0lVq-PdEGZciVoknB8nCUXCI.

The rest of the paper proceeds as follows. In Section 2, we present the proposed Bayesian nonparametric model along with the posterior inference. We evaluate the performance of the proposed approach through simulation studies and sensitivity analyses in Section 3. In Section 4, we apply the proposed model to a large-scale HIV clinical dataset to study the effects of ART regimens on depressive symptoms. Finally, we conclude with a discussion in Section 5.

2 Model and Inference

2.1 Probability Model

Denote to be the score of depression item for individual at visit , where , , and . Let denote the ART regimen used by individual at visit . For example, D4T+LAM+NFV if individual takes a combination of drugs D4T, LAM, and NFV at visit . Let be an

-dimensional vector including an intercept, time-invariant covariates (e.g., race), and time-varying covariates (e.g., BMI, CD4 count) for individual

at visit . We construct a sampling model for the depression score as follows,

(2.1)

where is a dimensional matrix, is a -dimensional vector-valued function, is a

-dimensional vector following a multivariate normal distribution

that models the dependency among different depression items, and is an independent normal error. For identifiability, we assume to be a correlation matrix. The first term in (2.1) captures the dependence of the outcome on the covariate . The second term is the key component of our model characterizing the combination effect of ART regimen , the details of which are given below.

Combination effect . We construct with two desired properties: 1) sharing-of-information - encouraging similar effects for similar ART regimens; and 2) parsimony - reducing the high dimensional ART regimen space to a manageable size. Specifically, we first pick a number of representative ART regimens, denoted by , which are similar to the notion of knots in splines. Then we model the combination effect of ART regimen using a kernel smoother,

(2.2)

where ’s are -dimensional vectors and the kernel weights are defined by an ART regimen similarity function . How to choose these representative regimens will be discussed later in the simulation studies and applications. The concept of similarity between different regimens has been introduced in HIV studies. Bogojeska et al. (2010) proposed a linear kernel method to compute the similarity between regimens based on the proportion of common drugs that two regimens share, , where is a binary vector indicating the drugs comprising the regimen and is a vector of 1’s. Although conceptually simple, the linear kernel treats all ART drugs as exchangeable and does not account for drug classes. For example, the regimen D4T (NRTI) + LAM (NRTI) + NFV (PI) should be more similar to D4T + LAM + ATZ (PI) than D4T + LAM + EFV (NNRTI) since NFV and ATZ belong to the same drug class PI whereas EFV belongs to another class NNRTI. However, the linear kernel approach gives rise to the same similarity score for these two pairs.

We propose to use a subset-tree (ST) kernel method, which was originally developed in natural language processing to represent sentence structure

(Collins and Duffy, 2002). We represent each ART regimen as a rooted tree which encodes the knowledge about the regimen structure such as drug classes and the number of distinct drug classes under each regimen. Figure 2 illustrates this idea using three regimens (A, B, C) as an example. Regimen A contains D4T (NRTI) + LAM (NRTI) + EFV (NNRTI); regimen B contains D4T (NRTI) + LAM (NRTI) + IDV (PI); and regimen C contains FTC (NRTI) + TDF (NRTI) + ATZ (PI) + RTV (PI). The linear kernel used in Bogojeska et al. (2010) would assign 0 similarity score to regimens A and C since they share no common drugs. But the ST kernel will be able to capture the similarity on the drug class level (highlighted by the yellow boxes). In fact, ST kernel calculates the similarity score between regimens across all levels of the tree representation. This feature will later be further exploited to compute the similarity between (longitudinal) sequences of regimens of possibly different lengths.

(a) Regimen A (b) Regimen B
(c) Regimen C (d) Similarity score matrix
Figure 2: Tree representations for ART regimens with their similarity matrix.

The main idea of the ST kernel is to compute the number of common substructures between two trees and . Let denote the set of nodes for any tree and let denote the set of children nodes of the node (i.e., nodes immediately below ). The similarity score, between two regimen trees and , is calculated by

(2.3)

where is defined for each pair of nodes as follows. (i) If and are terminal nodes , then . (ii) If and have different sets of children nodes , then . (iii) If and have the same nonempty set of children nodes, then , where is the cardinality of a set and is a child of for

. Here we include a hyperparameter

, which is a decay factor to control the relative influence from nodes near the root to alleviate the peakiness of the ST kernel when the depth of tree fragments is considerably large (Beck et al., 2015). Figure 2(d) presents the similarity score matrix among A, B, and C when . Note that the self-similarity of regimen C is higher than those of regimens A or B because regimen C consists of more drugs and therefore has a higher similarity score due to the additive definition.

To estimate ’s in (2.2), another challenge arises from the potential high-dimensionality and multicollinearity of the kernel weights calculated from the similarity function. Following the idea from principal component regression (Kendall et al., 1965)

, we consider a principal component analysis on the design matrix that consists of kernel weights. Specifically, let

be a -dimensional vector whose -th element is the kernel weight , and let be the design matrix () for the kernel regression in (2.2). We perform the principal component analysis on and retain the first

principal components that explain at least 99.9% of the total variance. The resulting

matrix is denoted by . Then the combination effect can be approximated by , where is the matrix that needs to be estimated.

2.2 Priors

To capture heterogeneity among individuals and individual treatment histories, we use the distance-dependent Chinese restaurant process (ddCRP, Blei and Frazier 2011; Dahl et al. 2017) to induce clustering of individuals depending on the similarity between their treatment histories. Let . We assume that with a mass parameter , a base measure , and a similarity function .

We give a brief introduction to the ddCRP below. Due to the discrete nature of ddCRP, are likely to have ties. Let denote the unique values of where . Let denote a partition of such that . The ties among ’s naturally give rise to a partition, i.e., if individual belongs to cluster , . Following Dahl et al. (2017), can be written as

(2.4)

Denote a permutation of , a partition of , and the number of subsets in ,

. The probability mass function of the partition

is defined as the product of increasing conditional probabilities (Dahl et al., 2017),

(2.5)

where

(2.6)

and for by convention.

Similarity function . The similarity function depends on individuals’ treatment histories. Let and denote two sequences of treatment regimens for individuals and

, respectively. The proposed ddCRP prior assumes that the prior probability that they belong to the same cluster is proportional to the similarity score

. Therefore, individuals with similar ART regimen histories are a priori more likely to be clustered together. To measure the similarity between two regimen sequences, we extend the ST kernel in by combining multiple regimen trees into a single tree under the common root “ART.” Figure 3 shows an example of the tree structure for one ART regimen sequence with three distinct ART regimens. Then we apply the ST kernel in to calculate the similarity score between regimen sequences in the same fashion as before.

Figure 3: Tree representation (bottom) for a sequence of ART regimens (top).

We complete the model by assigning hyperpriors. We use a conjugate gamma prior

on , a conjugate inverse-gamma prior on

, and a uniform distribution on the permutation, i.e.,

for all , for ease of posterior computation. In addition, we use a conjugate normal prior as the base measure . Specifically, let and be the -th row of and , respectively. We assume that and for , and , where , , , and . For the correlation matrix , we assume that following Lewandowski et al. (2009), where denotes the determinant of a matrix.

2.3 Posterior Inference

We carry out posterior inference with the Markov chain Monte Carlo (MCMC) algorithm. The posterior sampling for

’s, , and is straightforward through standard Gibbs sampler and Metropolis-Hastings sampler. To draw posterior samples for the parameters related to the ddCRP prior, the key step is to compute the full conditional distribution of the partition based on the probability mass function in and . Suppose at the current state, the partition is and let , , denote these subsets without individual . Let be the partition obtained by moving from its current subset to the subset . Here we let denote the index of a new empty subset , and let denote the partition after moving to a new subset. Then the full conditional distribution for the allocation of individual is given by,

for , where the new parameters are drawn from the base measure . Note that is calculated by evaluating and at the partition . More details of the MCMC can be found in the Supplementary Material Section A.

3 Simulation Study

In this section, we conducted simulation studies to evaluate the performance of the proposed model by comparing the posterior inference to the simulation truth. To demonstrate the advantages of using the ddCRP prior for taking into account individuals’ heterogeneity and treatment histories, and the ST kernel for inducing an appropriate ART regimen similarity, we compared the proposed model to two alternative methods. The first alternative replaces the ddCRP prior on with independent conjugate multivariate normal priors on ’s and ’s that do not take into account individuals’ heterogeneity and treatment histories, and replaces the ST kernel with a linear kernel (Bogojeska et al., 2010) based on the proportion of common drugs that two regimens share, , where is a binary vector indicating the drugs comprising the regimen and is a vector of 1’s. The linear kernel only considers the number of same drugs in each pair of regimens, but ignores the drug class information. We call this method Normal+Linear. The second alternative, called DP+Linear, replaces the ddCRP prior on with a Dirichlet process prior that does not account for individuals’ treatment histories and replaces the ST kernel with the linear kernel as in the first alternative method. Furthermore, the robustness of the decay factor in the ST kernel was demonstrated by sensitivity analyses.

3.1 Simulation setup

Assume that there were individuals with depression items and covariates with one intercept, one time-invariant covariate, and one time-varying covariate, i.e., , where ’s and ’s were generated from independent standard normal distributions, , . Individuals’ treatment histories were randomly sampled from the WIHS dataset without replacement, resulting in the number of visits per individual to range from 2 to 38. We set the simulated true decay factor , and then computed the similarity scores among individuals’ treatment histories. Based on the similarity scores among different individuals, we randomly generated one realization from the ddCRP prior, yielding the simulated true number of clusters to be and the number of individuals in each clusters is 67, 61, and 72, respectively. Figure 4(a) presents the simulated true clustering scheme. Conditional on the clustering memberships, we generated the simulated true from a standard multivariate normal distribution, the values of which are shown in Supplementary Table S1.

(a) Simulated true clustering scheme

(b) Posterior probabilities of individual co-clustering

Figure 4: Simulated true clustering scheme and posterior probabilities of individual co-clustering averaged over 100 repeated simulations.

We selected representative drug regimens if a regimen has been used in more than 10 visits among all the 200 individuals, yielding . We generated the simulated true in the kernel regression from a standard multivariate normal distribution, and computed the kernel weight matrix by applying the ST kernel in (2.3) to individuals’ treatment regimens and the selected representative drug regimens. We performed the principal component analysis on the design matrix , and chose the first principal components for that explains 99.9% variation of the original matrix. We set the non-diagonal elements of the correlation matrix to be , and . Lastly, we generated the depression scores from .

We applied the proposed model to the simulated dataset with 100 repeated simulations. The hyperparameters were set to be , , , , , , , , , and . For each analysis, we ran 10,000 MCMC iterations with an initial burn-in of 5,000 iterations and a thinning factor of 10. Convergence diagnostic assessed using R package coda, including autocorrelation plots and trace plots (Figures S2 and S3 in the Supplementary Material) of the post-burn-in MCMC samples for some randomly selected parameters, showed no issues of non-convergence.

3.2 Simulation results

We first report on the performance in terms of recovering the individual clustering. Our model successfully identified as it only overestimated the true number of clusters by 1 in 2.24% of the post-burn-in MCMC posterior samples among all the 100 repeated simulations. We further calculated the posterior probabilities of individual co-clustering based on the empirical proportions of individuals being clustered in the same cluster over the post-burn-in MCMC samples. The co-clustering probability matrix averaged over 100 repeated simulations is shown in Figure 4(b), indicating that the proposed method assigns individuals to their simulated true clusters with high probabilities.

Next, we examine whether we can recover the drug combination effect . We randomly selected one simulated dataset from 100 repeated simulations, and plotted the histogram of the true drug combination effects overlaid with the empirical density of the posterior expected combination effects in Figure S4 of the Supplementary Material. As for individual-specific drug combination effects across visits, Figure 5 compares the simulation truths and the estimated combination effects for two randomly selected individuals in this dataset. Both Figure S4 and Figure 5 show that our model can well recover the drug combination effects.

(a) Combination effects for individual 1 (b) Combination effects for individual 2
Figure 5: Combination effects for two randomly selected individuals from one randomly selected simulated dataset. The horizontal axis is the index of visit, and the vertical axis is the combination effect. The black lines represent the simulated truths of combination effects, the green lines represent the estimations under the Normal+Linear method, the blue lines represent the estimations under the DP+Linear method, and the red lines represent the estimations under the proposed method (ddCRP+ST). The shaded area represents the posterior 95% credible bands under the proposed method.

For parameter estimation, Figure S6 in the Supplementary Material plots the 95% estimated credible intervals (CI) for

’s using the same simulated dataset, where the triangles represent the simulation truths. As shown in Figure S6, all the 95% CI are centered around the simulated true values. As another metric of performance, we computed, for each simulated dataset, the mean squared error (MSE) taken as the averaged squared errors between the post-burn-in MCMC posterior samples and the simulated truth. Table S5 in the Supplementary Material summarizes the mean and standard deviation of MSE across 100 simulated datasets for

’s. Both Table S5 and Figure S6 show that the proposed method performs well in terms of estimating the parameter values.

In addition, we compared the proposed method to two alternatives: the Normal+Linear and the DP+Linear methods. Figure 5 compares the estimated combination effects under the proposed model to those under the two alternative methods. The proposed method with the ddCRP prior and the ST kernel well recovered the ground truth, while both the Normal+Linear and DP+Linear methods had larger bias in estimating the drug combination effects.

Lastly, to explore the sensitivity of the posterior inference with respect to the decay factor , we conducted inference under several values of = 0.1, 0.3, 0.5, 0.8, 1 for one randomly selected simulated dataset. The decay factor was originally introduced in natural language processing to alleviate the peakiness of the ST kernel when the depth of the tree fragments is considerably large, in which case self similarities are disproportionately larger than similarities between two different trees. Therefore, the decay factor , which down-weights the contribution of large tree fragments to the kernel exponentially with their sizes, could have significant influence on the inference if the tree structure is deep. However, this is not the case in our application with relatively shallow trees. Figure 6 compares the parameter estimations under different values of , showing that there is no significant difference among all these experiments.

(a) Cluster 1 (b) Cluster 2 (c) Cluster 3
Figure 6: credible intervals for randomly selected in the sensitivity analyses, where the triangles represent the simulated true values, and the colors represent different values of the decay factor .

4 Application: WIHS Data Analysis

The Women’s Interagency HIV Study (WIHS) is a multisite, longitudinal cohort study of women living with HIV and women at-risk for HIV in the United States (Barkan et al., 1998; Adimora et al., 2018). Full details of the study design and prospective data collection are described at https://statepi.jhsph.edu/wihs/wordpress. Participants provide biological specimens, complete physical examinations, and undergo extensive assessment of demographic, clinical, and behavioral data via interviews at each visit. Included in this assessment was the Center for Epidemiological Studies Depression Scale (CES-D, Radloff 1977), which is a self-report assessment of depressive symptoms spanning somatic (e.g., sleep and appetite difficulties), negative affect (e.g., loneliness and sadness), lack of positive affect (e.g., hopelessness), and interpersonal symptoms (e.g., people are unfriendly). For the present analysis, we included all women from the Washington, D.C. site in the WIHS with at least five visits and complete CES-D data, which yielded individuals. We also extracted the following sociodemographic, behavioral, and clinical risk factors for depressive symptoms: age, race, smoking status, substance use (e.g., marijuana, cocaine, and heroin), body mass index (BMI), hypertension, CD4 count, and viral load. We selected representative ART regimens in (2.2) using the same criterion as in the simulation study. In particular, these representative ART regimens are combinations of 24 ART agents in five drug classes: NRTI, NNRTI, PI, INSTI, and EI.

We applied the proposed model to the WIHS dataset using the same hyperparameters as in the simulation study and set the decay factor to be . We performed the principal component analysis on the kernel weight matrix based on these 87 representative ART regimens, and selected the first principal components that explain 99.9% variation of the original matrix. We used 5,000 post burn-in samples after 5,000 iterations with a thinning factor of 10 for posterior inference. The proposed model identified three clusters, with the number of women in each cluster being 132, 84, and 43 respectively. Table S7 in the Supplementary Material summarizes the demographic, clinical, and behavioral characteristics of women in the three clusters at their initial visits, and Table S8 reports the frequency of the 24 ART agents and their corresponding drug classes used by women in the three clusters, respectively.

Figure 7 summarizes the posterior means and the corresponding credible intervals of the estimated coefficients with respect to age, CD4 count, viral load, and substance use on four depression items in each cluster. As seen in Figure 7, the effects of covariates on depressive symptoms were distinct among the three clusters. Panel (a) shows that younger people had higher depressive symptoms in cluster 2, but lower depressive symptoms in clusters 1 and 3. Panels (b) and (c) indicate that higher CD4 and lower viral load are associated with lower depressive symptoms. Panel (d) shows a positive relationship between substance use and depressive symptoms. These findings are consistent with the literature (Berg et al., 2007; Springer et al., 2009; Grov et al., 2010; Taniguchi et al., 2014).

(a) Age (b) CD4 count
(c) Viral load (d) Substance use
Figure 7: Posterior means and CIs for the estimated coefficients corresponding to age, CD4 count, viral load, and substance use in the real data analysis. The dots represent the posterior means and the colors indicate different depressive symptoms.

Next, we report the effects of ART regimens, i.e., drug combinations, on depressive symptoms in each cluster. Figure 8 plots the association between ART regimens and depressive symptoms with respect to the first two principal components in each cluster. To explore the patterns and interpret the estimated drug combination effects, we further list the top five positively and negatively related ART regimens for each principal component in terms of the coefficients of the loading matrix in Table 1. As shown in Figure 8, the first principal component was negatively associated with all the depressive symptoms in cluster 1 and 3, but had little effects in cluster 2. In addition, the first principal component was positively associated with ART regimens consisting of two NRTI drugs FTC + TDF, an NNRTI drug EFV, RPV or NVP, and an additional INSTI drug RAL (Table 1), which indicates a beneficial or protective effect for these ART regimens on depressive symptoms. In fact, combining two NRTI drugs as backbone with an additional NNRTI drug was recommended as one of the first-line therapies (Günthard et al., 2014), and previous clinical studies also reported that RAL was well-tolerated and provided desirable viral suppression when used with certain NRTIs such as TDF (Grinsztejn et al., 2007; Markowitz et al., 2007). Conversely, negative relationships were observed between the first principle component and ART regimens consisting of two NRTI drugs AZT + LAM and a PI drug such as LPV, revealing worse depressive symptoms for women using these drug combinations. Rabaud et al. (2005) reported that a large proportion of individuals receiving AZT + LAM + LPV experienced serious adverse effects, especially gastrointestinal side effects such as nausea and vomiting, leading to poor tolerability of this regimen and treatment discontinuation. Furthermore, the second principal component was positively associated with depressive symptoms in clusters 1 and 3. ART regimens consisting of two NRTI drugs AZT + LAM and an NNRTI drug such as EFV were positively related to the second principal component, while regimens consisting of two NRTI drugs FTC + TDF and two PI drugs such as ATZ + RTV were negatively related to the second principal component. Therefore, a combination of AZT, LAM, and EFV was estimated to have adverse effects on depressive symptoms whereas a combination of FTC, TDF, ATZ and RTV was estimated to have beneficial effects. Indeed, Gallant et al. (2006) reported more frequent adverse effects and treatment discontinuation when individuals were on EFV combined with AZT and LAM instead of FTC and TDF. Conversely, adding PI drugs ATZ and RTV to NRTI drugs FTC and TDF yields both significant antiviral efficacy and safety (Soriano et al., 2011).

(a) Principal component (b) Principal component
Figure 8: Posterior means and CIs for the estimated combination effects on four depressive symptoms with respect to the first two principal components in the WIHS data analysis. The dots represent the posterior means and the colors indicate different depressive symptoms.
ART Regimens Loading Coefficients
lightgray Principal Components
FTC + TDF + EFV 0.182
FTC + TDF + RPV 0.182
FTC + TDF + NVP 0.181
FTC + TDF + EFV + RAL 0.174
DDI + TDF + EFV 0.162
AZT + LAM + LPV -0.171
AZT + LAM + SQV -0.165
DDI + LAM + LPV -0.163
ABC + LAM + LPV -0.162
AZT + LAM + IDV -0.162
lightgray Principal Components
AZT + LAM + EFV 0.179
AZT + LAM + NVP 0.178
ABC + AZT + LAM + EFV 0.169
ABC + AZT + EFV 0.159
AZT + DDI + NVP 0.156
FTC + TDF + DRV + RTV -0.191
FTC + TDF + FPV + RTV -0.190
FTC + TDF + ATZ + RTV -0.189
DDI + TDF + ATZ + RTV -0.182
FTC + TDF + ATZ -0.175
Table 1: Top five positively and negatively related ART regimens for the first two principal components in terms of the coefficients of the loading matrix.

The U.S. Department of Health and Human Services provides general guidelines on ART treatments; however, these guidelines do not take into account individual heterogeneity and treatment histories. To make clinical decisions tailored to each person (precision medicine), understanding the individualized adverse effect of each possible drug combination will be one of the key contributors. The proposed method can accurately predict individuals’ adverse effects of ART based on their clinical profiles, which can help guide clinicians to prescribe ART regimens. For illustration, we randomly selected an individual from the WIHS dataset with seven visits in total, who started AZT (NRTI) at the first visit, added LAM (NRTI) at the second visit, and used the drug combination AZT + LAM + SQV (PI) from her fourth to sixth visits. Then we considered two hypothetical scenarios. In the first scenario, we assumed that the individual kept using the similar NRTI + PI drug combination as before but only replaced the PI drug SQV with a different PI drug LPV. In the second scenario, this individual was switched to a distinct NRTI + NNRTI drug combination FTC (NRTI) + TDF (NRTI) + EFV (NNRTI). Figure 9

plots the posterior predictive depression scores for this individual at the last visit based on the information from her previous six visits under the two hypothetical scenarios. As shown in Figure

9(c), there were no significant differences between using different PI drugs when combined with NRTI drugs AZT and LAM as the backbone treatment. However, using NRTI drugs FTC and TDF as the backbone with NNRTI drug EFV demonstrated superior performance on alleviating depressive symptoms. As a result, we would recommend the clinician to select the ART regimen FTC + TDF + EFV instead of AZT + LAM + LPV for this particular individual. This example demonstrates that the proposed method has the potential to guide more informed and effective personalized medicine in HIV clinical practice.

(a) ART use for scenario 1 (b) ART use for scenario 2
(c) Predictive depression score for scenario 1 (d) Predictive depression score for scenario 2
Figure 9: Predictive depression scores for an individual in the WIHS dataset with two different hypothetical scenarios of ART medication use. The dashed lines represent the predictive 95% credible bands with respect to each depressive symptom.

To facilitate the implementation of the proposed method in the decision process of HIV clinicians, and for broad application in personalized medicine, we have created an interactive web application to illustrate this example using R package shiny (Chang et al., 2019), available at https://wjin.shinyapps.io/Rshiny/. The web user interface interactively displays the predictive depression scores of an individual in response to the user’s choice of the individual’s clinical characteristics and ART medication use. Figure S9 in the Supplementary Material shows a screenshot of the web application.

5 Conclusion

To facilitate a precision medicine approach, we proposed a novel Bayesian nonparametric approach to estimate the effects of ART regimens on depressive symptoms. The method is built upon the ST kernel method that quantifies similarities among ART regimens and the ddCRP that accounts for individuals’ heterogeneity in both treatment histories and clinical characteristics. Through simulation studies and analysis of the WIHS dataset, we have demonstrated that the proposed model can accurately estimate the drug combination effects and yield meaningful and interpretable results.

There are several potential extensions. First, the current similarity score is parameterized by a hyperparameter . We could impose a prior on and estimate it from the posterior inference. It will require us to develop more efficient posterior samplers because in each iteration of MCMC the similarity matrix needs to be recalculated at the current value of . Second, the similarity between ART regimens may also depend on the individuals’ socio-demographic, behavioral, and clinical characteristics. We could extend the model to account for these factors by modifying the parameter in (2.2) as a function of these variables. Finally, combination therapies are needed for many complex diseases beyond HIV such as cancer and chronic diseases. Each chronic condition requires long-term medication use. The proposed method can be applied to such electronic health records datasets (Gill et al., 2010) to examine the side effects of combination therapies, potentially yielding better therapy management for the elderly population and which has the potential to reduce public healthcare costs.

Acknowledgment

This work was supported by the Johns Hopkins University Center for AIDS Research NIH/NIAID fund (P30AI094189) 2019 faculty development award to Dr. Xu, NSF 1940107 to Dr. Xu, NSF DMS1918854 to Drs. Xu and Rubin, and NSF DMS1918851 to Dr. Ni.

References

  • Adimora et al. (2018) Adimora, A. A., Ramirez, C., Benning, L., Greenblatt, R. M., Kempf, M.-C., Tien, P. C., Kassaye, S. G., Anastos, K., Cohen, M., Minkoff, H., et al. (2018). Cohort profile: the women’s interagency HIV study (WIHS). International journal of epidemiology, 47(2):393–394i.
  • Altmann et al. (2007) Altmann, A., Beerenwinkel, N., Sing, T., Savenkov, I., Däumer, M., Kaiser, R., Rhee, S.-Y., Fessel, W. J., Shafer, R. W., and Lengauer, T. (2007). Improved prediction of response to antiretroviral combination therapy using the genetic barrier to drug resistance. Antiviral therapy, 12(2):169.
  • Altmann et al. (2009) Altmann, A., Däumer, M., Beerenwinkel, N., Peres, Y., Schülter, E., Büch, J., Rhee, S.-Y., Sönnerborg, A., Fessel, W. J., Shafer, R. W., et al. (2009). Predicting the response to combination antiretroviral therapy: retrospective validation of geno2pheno-theo on a large clinical database. The Journal of infectious diseases, 199(7):999–1006.
  • Arenas-Pinto et al. (2016) Arenas-Pinto, A., Thompson, J., Musoro, G., Musana, H., Lugemwa, A., Kambugu, A., Mweemba, A., Atwongyeire, D., Thomason, M. J., Walker, A. S., and Paton, N. I. (2016). Peripheral neuropathy in hiv patients in sub-saharan africa failing first-line therapy and the response to second-line art in the earnest trial. Journal of Neurovirology, 22(1):104–113.
  • Bacon et al. (2005) Bacon, M. C., Von Wyl, V., Alden, C., Sharp, G., Robison, E., Hessol, N., Gange, S., Barranday, Y., Holman, S., Weber, K., and Young, M. A. (2005). The women’s interagency hiv study: an observational cohort brings clinical sciences to the bench. Clinical and diagnostic laboratory immunology, 12(9):1013–1019.
  • Bangsberg et al. (2001) Bangsberg, D. R., Perry, S., Charlebois, E. D., Clark, R. A., Roberston, M., Zolopa, A. R., and Moss, A. (2001). Non-adherence to highly active antiretroviral therapy predicts progression to aids. Aids, 15(9):1181–1183.
  • Barkan et al. (1998) Barkan, S. E., Melnick, S. L., Preston-Martin, S., Weber, K., Kalish, L. A., Miotti, P., Young, M., Greenblatt, R., Sacks, H., and Feldman, J. (1998). The women’s interagency hiv study. Epidemiology, pages 117–125.
  • Beck et al. (2015) Beck, D., Cohn, T., Hardmeier, C., and Specia, L. (2015). Learning structural kernels for natural language processing. Transactions of the Association for Computational Linguistics, 3:461–473.
  • Bengtson et al. (2016) Bengtson, A. M., Pence, B. W., Crane, H. M., Christopoulos, K., Fredericksen, R. J., Gaynes, B. N., Heine, A., Mathews, W. C., Moore, R., Napravnik, S., Safren, S., and Mugavero, M. J. (2016). Disparities in depressive symptoms and antidepressant treatment by gender and race/ethnicity among people living with hiv in the united states. PloS one, 11(8):e0160738.
  • Berg et al. (2007) Berg, C. J., Michelson, S. E., and Safren, S. A. (2007). Behavioral aspects of hiv care: adherence, depression, substance use, and hiv-transmission behaviors. Infectious Disease Clinics of North America, 21(1):181–200.
  • Blei and Frazier (2011) Blei, D. M. and Frazier, P. I. (2011). Distance dependent Chinese restaurant processes.

    Journal of Machine Learning Research

    , 12(Aug):2461–2488.
  • Bogojeska et al. (2010) Bogojeska, J., Bickel, S., Altmann, A., and Lengauer, T. (2010). Dealing with sparse data in predicting outcomes of HIV combination therapies. Bioinformatics, 26(17):2085–2092.
  • Brickman et al. (2017) Brickman, C., Propert, K. J., Voytek, C., Metzger, D., and Gross, R. (2017). Association between depression and condom use differs by sexual behavior group in patients with hiv. AIDS and Behavior, 21(6):1676–1683.
  • Caniglia et al. (2017) Caniglia, E. C., Cain, L. E., Sabin, C. A., Robins, J. M., Logan, R., Abgrall, S., Mugavero, M. J., Hernández-Díaz, S., Meyer, L., Seng, R., et al. (2017). Comparison of dynamic monitoring strategies based on CD4 cell counts in virally suppressed, HIV-positive individuals on combination antiretroviral therapy in high-income countries: A prospective, observational study. The Lancet HIV, 4(6):e251–e259.
  • Carleton et al. (2013) Carleton, R. N., Thibodeau, M. A., Teale, M. J., Welch, P. G., Abrams, M. P., Robinson, T., and Asmundson, G. J. (2013). The center for epidemiologic studies depression scale: a review with a theoretical and empirical examination of item content and factor structure. PloS one, 8(3):e58067.
  • Chang et al. (2019) Chang, W., Cheng, J., Allaire, J., Xie, Y., and Mcpherson, J. (2019). shiny: web application framework for r. r package version 1.4.0.
  • Chattopadhyay et al. (2017) Chattopadhyay, S., Ball, S., Kargupta, A., Talukdar, P., Roy, K., Talukdar, A., and Guha, P. (2017). Cognitive behavioral therapy improves adherence to antiretroviral therapy in hiv-infected patients: a prospective randomized controlled trial from eastern india. HIV & AIDS Review. International Journal of HIV-Related Problems, 16(2):89–95.
  • Collins and Duffy (2002) Collins, M. and Duffy, N. (2002). Convolution kernels for natural language. In Advances in neural information processing systems, pages 625–632.
  • Dahl et al. (2017) Dahl, D. B., Day, R., and Tsai, J. W. (2017). Random partition distribution indexed by pairwise information. Journal of the American Statistical Association, 112(518):721–732.
  • de Olalla Garcia et al. (2002) de Olalla Garcia, P., Knobel, H., Carmona, A., Guelar, A., López-Colomés, J. L., and Caylà, J. A. (2002). Impact of adherence and highly active antiretroviral therapy on survival in hiv-infected patients. Journal of acquired immune deficiency syndromes (1999), 30(1):105–110.
  • Escobar and West (1995) Escobar, M. D. and West, M. (1995). Bayesian density estimation and inference using mixtures. Journal of the american statistical association, 90(430):577–588.
  • Fazeli et al. (2015) Fazeli, P. L., Marquine, M. J., Dufour, C., Henry, B. L., Montoya, J., Gouaux, B., Moore, R. C., Letendre, S. L., Woods, S. P., Grant, I., Jeste, D. V., and Moore, D. J. (2015). Physical activity is associated with better neurocognitive and everyday functioning among older adults with HIV disease. AIDS and Behavior, 19(8):1470–1477.
  • Gaida et al. (2016) Gaida, R., Truter, I., Grobler, C., Kotze, T., and Godman, B. (2016). A review of trials investigating efavirenz-induced neuropsychiatric side effects and the implications. Expert review of anti-infective therapy, 14(4):377–388.
  • Gallant et al. (2006) Gallant, J. E., DeJesus, E., Arribas, J. R., Pozniak, A. L., Gazzard, B., Campo, R. E., Lu, B., McColl, D., Chuck, S., Enejosa, J., et al. (2006). Tenofovir df, emtricitabine, and efavirenz vs. zidovudine, lamivudine, and efavirenz for hiv. New England Journal of Medicine, 354(3):251–260.
  • Gill et al. (2010) Gill, J. M., Klinkman, M. S., and Chen, Y. X. (2010). Antidepressant medication use for primary care patients with and without medical comorbidities: a national electronic health record (ehr) network study. The Journal of the American Board of Family Medicine, 23(4):499–508.
  • Grinsztejn et al. (2007) Grinsztejn, B., Nguyen, B.-Y., Katlama, C., Gatell, J. M., Lazzarin, A., Vittecoq, D., Gonzalez, C. J., Chen, J., Harvey, C. M., Isaacs, R. D., et al. (2007). Safety and efficacy of the hiv-1 integrase inhibitor raltegravir (mk-0518) in treatment-experienced patients with multidrug-resistant virus: a phase ii randomised controlled trial. The Lancet, 369(9569):1261–1269.
  • Grov et al. (2010) Grov, C., Golub, S. A., Parsons, J. T., Brennan, M., and Karpiak, S. E. (2010). Loneliness and hiv-related stigma explain depression among older hiv-positive adults. AIDS care, 22(5):630–639.
  • Günthard et al. (2014) Günthard, H. F., Aberg, J. A., Eron, J. J., Hoy, J. F., Telenti, A., Benson, C. A., Burger, D. M., Cahn, P., Gallant, J. E., Glesby, M. J., et al. (2014). Antiretroviral treatment of adult HIV infection: 2014 recommendations of the international antiviral society–USA panel. JAMA, 312(4):410–425.
  • Ironson et al. (2017) Ironson, G., Fitch, C., and Stuetzle, R. (2017).

    Depression and survival in a 17-year longitudinal study of people with hiv: Moderating effects of race and education.

    Psychosomatic medicine, 79(7):749–756.
  • Jelsma et al. (2005) Jelsma, J., MacLean, E., Hughes, J., Tinise, X., and Darder, M. (2005). An investigation into the health-related quality of life of individuals living with hiv who are receiving haart. AIDS Care, 17(5):579–588. PMID: 16036244.
  • Kendall et al. (1965) Kendall, M. G. et al. (1965).

    course in multivariate analysis.

  • Larder et al. (2007) Larder, B., Wang, D., Revell, A., Montaner, J., Harrigan, R., De Wolf, F., Lange, J., Wegner, S., Ruiz, L., Pérez-Elías, M. J., et al. (2007). The development of artificial neural networks to predict virological response to combination HIV therapy. Antiviral therapy, 12(1):15.
  • Lewandowski et al. (2009) Lewandowski, D., Kurowicka, D., and Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of multivariate analysis, 100(9):1989–2001.
  • Lundgren et al. (2002) Lundgren, J. D., Mocroft, A., Gatell, J. M., Ledergerber, B., Monforte, A. D., Hermans, P., Goebel, F.-D., Blaxhult, A., Kirk, O., and Phillips, A. N. (2002). A clinically prognostic scoring system for patients receiving highly active antiretroviral therapy: results from the eurosida study. The Journal of infectious diseases, 185(2):178–187.
  • Markowitz et al. (2007) Markowitz, M., Nguyen, B.-Y., Gotuzzo, E., Mendo, F., Ratanasuwan, W., Kovacs, C., Prada, G., Morales-Ramirez, J. O., Crumpacker, C. S., Isaacs, R. D., et al. (2007). Rapid and durable antiretroviral effect of the hiv-1 integrase inhibitor raltegravir as part of combination therapy in treatment-naive patients with hiv-1 infection: results of a 48-week controlled study. JAIDS Journal of Acquired Immune Deficiency Syndromes, 46(2):125–133.
  • Mollan et al. (2014) Mollan, K. R., Smurzynski, M., Eron, J. J., Daar, E. S., Campbell, T. B., Sax, P. E., Gulick, R. M., Na, L., O’keefe, L., Robertson, K. R., and Tierney, C. (2014). Association between efavirenz as initial therapy for HIV-1 infection and increased risk for suicidal ideation or attempted or completed suicide: an analysis of trial data. Annals of internal medicine, 161(1):1–10.
  • Okeke and Wagner (2013) Okeke, E. N. and Wagner, G. J. (2013). Aids treatment and mental health: evidence from uganda. Social science & medicine, 92:27–34.
  • Pearson et al. (2009) Pearson, C., Micek, M., Pfeiffer, J., Montoya, P., Matediane, E., Jonasse, T., Cunguara, A., Rao, D., and Gloyd, S. (2009). One year after art initiation: psychosocial factors associated with stigma among hiv-positive mozambicans. AIDS and Behavior, 13(6):1189.
  • Rabaud et al. (2005) Rabaud, C., Burty, C., Grandidier, M., Christian, B., Penalba, C., Béguinot, I., Jeanmaire, H., and May, T. (2005). Tolerability of postexposure prophylaxis with the combination of zidovudine-lamivudine and lopinavir-ritonavir for hiv infection. Clinical infectious diseases, 40(2):303–305.
  • Radloff (1977) Radloff, L. S. (1977). The CES-D scale: A self-report depression scale for research in the general population. Applied psychological measurement, 1(3):385–401.
  • Rekha (2018) Rekha, M. M. (2018). A case study & report on drug related adverse effect observed and addressed during optimal therapy by doctor of pharmacy, in an antiretroviral therapy ward of a tertiary care teaching hospital. Research & Reviews: Journal of Medicine, 8(1):1–5.
  • Saag et al. (2018) Saag, M. S., Benson, C. A., Gandhi, R. T., Hoy, J. F., Landovitz, R. J., Mugavero, M. J., Sax, P. E., Smith, D. M., Thompson, M. A., Buchbinder, S. P., del Rio, C., Eron, J. J., Fätkenheuer, G., Günthard, H. F., Molina, J.-M., Jacobsen, D. M., and Volberding, P. A. (2018). Antiretroviral drugs for treatment and prevention of hiv infection in adults: 2018 recommendations of the international antiviral society–usa panel. JAMA, 320(4):379–396.
  • Soriano et al. (2011) Soriano, V., Arastéh, K., Migrone, H., Lutz, T., Opravil, M., Andrade-Villanueva, J., Antunes, F., Di Perri, G., Podzamczer, D., Taylor, S., et al. (2011). Nevirapine versus atazanavir/ritonavir, each combined with tenofovir disoproxil fumarate/emtricitabine, in antiretroviral-naive hiv-1 patients: the arten trial. Antiviral therapy, 16(3):339.
  • Springer et al. (2009) Springer, S. A., Chen, S., and Altice, F. (2009). Depression and symptomatic response among HIV-infected drug users enrolled in a randomized controlled trial of directly administered antiretroviral therapy. AIDS care, 21(8):976–983.
  • Taniguchi et al. (2014) Taniguchi, T., Shacham, E., Önen, N. F., Grubb, J. R., and Overton, E. T. (2014). Depression severity is associated with increased risk behaviors and decreased CD4 cell counts. AIDS care, 26(8):1004–1012.
  • Trickey et al. (2017) Trickey, A., May, M. T., Vehreschild, J.-J., Obel, N., Gill, M. J., Crane, H. M., Boesecke, C., Patterson, S., Grabar, S., Cazanave, C., Cavassini, M., Shepherd, L., d’Arminio Monforte, A., Sighem, A., Saag, M., Lampe, F., Hernando, V., Montero, M., Zangerle, R., Justice, A. C., Sterling, T., MIngle, S., and Sterne, J. A. (2017). Survival of HIV-positive patients starting antiretroviral therapy between 1996 and 2013: a collaborative analysis of cohort studies. The Lancet HIV, 4(8):e349–e356.
  • Watkins and Treisman (2015) Watkins, C. C. and Treisman, G. J. (2015). Cognitive impairment in patients with aids–prevalence and severity. HIV/AIDS (Auckland, NZ), 7:35.
  • Yun et al. (2005) Yun, L. W., Maravi, M., Kobayashi, J. S., Barton, P. L., and Davidson, A. J. (2005). Antidepressant treatment improves adherence to antiretroviral therapy among depressed hiv-infected patients. JAIDS Journal of Acquired Immune Deficiency Syndromes, 38(4):432–438.
  • Zash et al. (2018) Zash, R., Makhema, J., and Shapiro, R. L. (2018). Neural-tube defects with dolutegravir treatment from the time of conception. New England Journal of Medicine, 379(10):979–981.

A: Details of MCMC

A1: Summary of Model

where is the drug combination effect that will be estimated by a principal component regression, the covariates of which are derived from the first principal components of the kernel weight matrix. Let be the corresponded coefficients for the principal component regression, and ’s be the unique values of induced by ddCRP, where .

We complete the model by assuming , , and independently for , and , where , , , , , and .

A2: Posterior Computation

A2.1: Update , the partition induced by ddCRP

for and , where is a new and independent draw from .

A2.2: Update , the mass parameter of ddCRP

Following the idea of Escobar and West (1995)

, we assign a Gamma distribution prior on

and introduce an auxiliary variable . Then a closed form full conditional posterior will be available for the mass parameter , which is a mixture of Gamma distributions, i.e.,

A2.3: Update , the permutation of subjects

We use Metropolis–Hastings algorithm to update . We first propose a new permutation by shuffling some randomly chosen integers in the current permutation and leaving the rest in their current positions. As a symmetric proposal distribution, the proposed new is accepted with probability

since we are assuming a uniform distribution prior on the permutation.

A2.4: Update , the cluster specific parameters

  • Update , ,

    where , , and if .

  • Update , ,

    where , , and if .

A2.5: Update , , , , the hyper-parameters

  • Update ,

    where and .

  • Update ,

    where and .

  • Update ,

    where and .

  • Update ,

    where and .

A2.6: Update , the normal correlation term