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 drugresistant mutations, ultimately resulting in treatment failure. ART agents fall into several classes including nucleotide reverse transcriptase inhibitor (NRTI), nonnucleotide 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 diseaserelated morbidity and mortality (Trickey et al., 2017), numerous observational studies have reported ARTrelated 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 (ArenasPinto 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 ARTrelated 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 risktaking 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 ARTrelated 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, ARTrelated 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.
Largescale 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 atrisk 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 selfreport questionnaire using the Center for Epidemiological Studies Depression Scale where a higher score reflects worse symptoms (CESD, 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 
The complexity of longitudinal observations, a heterogeneous population, and dynamic and mixed ART assignments present four analytical and modeling challenges. The first challenge is highdimensionality. 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 highdimensional 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 highdimensional ART regimen space were observed in the WIHS, while the inference on the entire space is desired. Finally, there is the issue of nonstationarity. 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, treebased 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 subsettree kernels (Collins and Duffy, 2002) and the distancedependent 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 sociodemographic, behavioral, and clinical factors. The subsettree 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 subsettree kernel among ART regimens to sequences of ART regimens in a longitudinal setup, and use the distancedependent 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=1FB8o0cHx0lVqPdEGZciVoknB8nCUXCI.
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 largescale 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, timeinvariant covariates (e.g., race), and timevarying 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 vectorvalued 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) sharingofinformation  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 subsettree (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 
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 selfsimilarity 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 highdimensionality 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 firstprincipal 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 distancedependent 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.We complete the model by assigning hyperpriors. We use a conjugate gamma prior
on , a conjugate inversegamma 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 MetropolisHastings 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 timeinvariant covariate, and one timevarying 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 coclustering 
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 nondiagonal 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 burnin 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 postburnin MCMC samples for some randomly selected parameters, showed no issues of nonconvergence.
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 postburnin MCMC posterior samples among all the 100 repeated simulations. We further calculated the posterior probabilities of individual coclustering based on the empirical proportions of individuals being clustered in the same cluster over the postburnin MCMC samples. The coclustering 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 individualspecific 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 
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 postburnin 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 downweights 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 
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 atrisk 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 (CESD, Radloff 1977), which is a selfreport 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 CESD 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 burnin 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 
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 firstline therapies (Günthard et al., 2014), and previous clinical studies also reported that RAL was welltolerated 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 
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 
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 
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’ sociodemographic, 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 longterm 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 geno2phenotheo on a large clinical database. The Journal of infectious diseases, 199(7):999–1006.
 ArenasPinto et al. (2016) ArenasPinto, 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 subsaharan africa failing firstline therapy and the response to secondline 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). Nonadherence to highly active antiretroviral therapy predicts progression to aids. Aids, 15(9):1181–1183.
 Barkan et al. (1998) Barkan, S. E., Melnick, S. L., PrestonMartin, 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 hivtransmission 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ándezDíaz, S., Meyer, L., Seng, R., et al. (2017). Comparison of dynamic monitoring strategies based on CD4 cell counts in virally suppressed, HIVpositive individuals on combination antiretroviral therapy in highincome 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 hivinfected patients: a prospective randomized controlled trial from eastern india. HIV & AIDS Review. International Journal of HIVRelated 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ópezColomés, J. L., and Caylà, J. A. (2002). Impact of adherence and highly active antiretroviral therapy on survival in hivinfected 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 efavirenzinduced neuropsychiatric side effects and the implications. Expert review of antiinfective 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 hiv1 integrase inhibitor raltegravir (mk0518) in treatmentexperienced patients with multidrugresistant 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 hivrelated stigma explain depression among older hivpositive 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 17year 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 healthrelated 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érezElí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., MoralesRamirez, J. O., Crumpacker, C. S., Isaacs, R. D., et al. (2007). Rapid and durable antiretroviral effect of the hiv1 integrase inhibitor raltegravir as part of combination therapy in treatmentnaive patients with hiv1 infection: results of a 48week 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 HIV1 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 hivpositive 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 zidovudinelamivudine and lopinavirritonavir for hiv infection. Clinical infectious diseases, 40(2):303–305.
 Radloff (1977) Radloff, L. S. (1977). The CESD scale: A selfreport 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., AndradeVillanueva, 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 antiretroviralnaive hiv1 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 HIVinfected 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 HIVpositive 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 hivinfected 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). Neuraltube 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 hyperparameters

Update ,
where and .

Update ,
where and .

Update ,
where and .

Update ,
where and .
A2.6: Update , the normal correlation term
where