In 2014, the United Nations Joint Programme on HIV and AIDS (UNAIDS) set the 90-90-90 targets as the primary strategy to end the HIV/AIDS epidemic (Joint United Nations Programme on HIV/AIDS 2014), which includes identifying 90% of people living with HIV through expanded testing, placing 90% of positively identified individuals on antiretroviral therapy, and ensuring that 90% of those on therapy can achieve undetectable viral loads by 2020. Considerable progress has been made toward Tri90 (Gaolathe et al. 2016;Gisslen et al. 2017; Joint United Nations Programme on HIV/AIDS 2016, 2017; Marukutira et al. 2018). Yet, as of 2019, there was a significant gap; instead of the targets of , it was 81-67-59 globally (Joint United Nations Programme on HIV/AIDS 2020a); and now the COVID-19 crisis has the potential to undermine existing efforts toward the HIV/AIDS epidemic (Joint United Nations Programme on HIV/AIDS 2020b).
To end HIV/AIDS epidemic, we shall learn from the past efforts toward Tri90 and identify important features that could guide more targeted and effective health policies. Novel datasets and sophisticated modeling tools are needed to enhance our understanding of Tri90 achievements. The Population-based HIV Impact Assessment (PHIA) survey is a nationally representative HIV-focused survey that started data collection in 2015. It is designed to measure the reach and impact of HIV programs.
We believe it offers a unique opportunity to better understand the key factors that drive the HIV epidemics in the most affected countries in sub-Saharan Africa. In this article, we propose a novel causal structural learning algorithm to discover important covariates and potential causal pathways for the Tri90. Causal structural learning aims to build a directed acyclic graph (DAG) that shows direct causal relations among variables of interest in a given domain. The resulting DAG helps us to understand the mechanism behind data.
Many classical structural learning algorithms are constraintbased, such as the PC algorithm (Spirtes and Glymour 1991), the PC-stable algorithm (Colombo and Maathuis 2014), and the MMPC algorithm (Tsamardinos, Aliferis, and Statnikov 2003b). The constraint-based algorithms learn graphical structures by d-separation set searching. D-separation is an important graphical concept for causal structural learning (Geiger and Pearl 1990;Geiger, Verma, and Pearl 1990), where the "d" stands for dependence or directed. Roughly speaking, in a directed acyclic graph (DAG), two vertices X and Y are d-separated by a set of vertices Z = {Z 1 , . . . , Z d } if and only if all the paths/information between X and Y are blocked by vertices in Z (see Section 2.1 for a rigorous definition of d-separation). The d-separation relationship in a DAG can be related to the conditional independence through the Markov condition and the faithfulness assumption (Spirtes 2010) (see definitions in Section 2.1). That is to say, vertices X and Y are d-separated by a vertex set Z = {Z 1 , . . . , Z d } is equivalent to the conditional independence of the corresponding variables X and Y given the corresponding set of variables Z = {Z 1 , . . . , Z d } under the Markov condition and the faithfulness assumption.
However, existing constraint-based structural learning algorithms can be quite aggressive in edge removal: if a Type II error is made such that two connected covariates X and Y are thought to be conditionally independent given some Z, the edge between X and Y is removed mistakenly, and useful information about important features and possible causal pathways get lost during this edge removal process. Especially in the case of categorical variables and relatively small sample sizes, the conditional independence tests used by structural learning algorithms can have high Type II error rates that lead to many false edge-removals and a severe information loss. As illustrated in the PHIA data analysis, very few features are connected by using the constraintbased structural learning algorithms. Results and discussions in detail can be found later in Section 3.
In literature dealing with the unreliable conditional independence tests, most literature focuses on their negative effects in the orientation procedure but much less on those in the skeleton learning procedure, with a few exceptions. Bromberg and Margaritis (2009) proposed a method to resolve the inconsistencies in the conditional independence tests in the skeleton learning procedure. Their idea is to deduce a "preference" score on the conditional independence test results from a certain set of axioms, such as Pearl's axioms (Dawid 1979;Pearl 1988), which all true conditional independence relationships should follow. These axioms can be seen as integrity constraints that can avoid certain inconsistent test outcomes. However, the computational complexity of the algorithm proposed by Bromberg and Margaritis (2009) significantly increases as the number of vertices increases. Thus, the algorithm cannot be employed on the MPHIA dataset for which the graph is not sparse enough.
In this article, we propose a new causal structural learning algorithm that aims to preserve more information on important features and potential causal pathways. We apply the proposed algorithm on the Malawi PHIA (Ministry of Health, Malawi et al. 2017) dataset, which we refer to as MPHIA, and obtain interesting findings related to Tri90 pathways. We further compare and validate the proposed algorithm with some classical structural learning algorithms using information criteria and simulations.
The remaining part of the article is organized as follows. In Section 2.1, we provide definitions for important concepts in causal structural learning, such as d-separation, and also a basic summary for the graphical notations used in the article. In Section 2.2, we use a simple example to illustrate the problem of aggressive edge removal of existing structural learning algorithms. In Sections 2.3-2.6, we propose a new causal structural learning algorithm to overcome the aggressive edge removal issue. In Section 3, we apply the proposed algorithm on the MPHIA dataset, compare the results obtained from the proposed algorithm with those of the existing algorithms in Section 3.1, and discuss the discovered Tri90 pathways in detail in Section 3.2. In Section 4, we compare the numerical performance of the proposed algorithm with classical structural learning algorithms in simulation studies. Section 5 provides a summary and discussion for the article. To save space, technical details, additional numerical results, and the relevant codebooks are provided in the supplementary materials.
In this section, we propose our causal structural learning algorithm. We first provide definitions for important concepts in causal structural learning, such as d-separation, and also a basic summary for the graphical notations used in the article in Section 2.1. We then discuss the problem of aggressive edge removal of existing structural learning algorithms with a simplified example in Section 2.2. To deal with this problem, we propose our new algorithm and provide a high-level overview of the algorithm in Section 2.3. The proposed algorithm consists of two main steps, the forward step and the maximization step, which are explained in detail in Sections 2.4 and 2.5, respectively. Section 2.6 provides the orientation procedure of our proposed algorithm and also summarizes the proposed method.
Suppose G is a directed acyclic graphical (DAG) model which represents the joint probability distribution over the vertex set V with directed edges and no directed loop. Each vertex in the graph represents a variable. We use variables X, Y, etc., to refer to the variables corresponding to the vertices X, Y, etc., and use the edge X → Y or Y ← X to refer to the directed edge from X to Y, where X is a parent vertex of Y, and Y is a child of X. X -Y denotes an undirected edge that could be either X → Y or Y → X. We use a path to refer to an acyclic sequence of adjacent vertices and a causal path from X to Y to refer to a path that all arrows are pointing away from X and into Y. If there is a causal path from X to Y, we say that X is an ancestor of Y and that Y is a descendant of X.
Next, we provide the formal definition of d-separation (Geiger and Pearl 1990;Geiger, Verma, and Pearl 1990). A collider is a vertex on a path with two incoming arrows. More specifically, a vertex Z is a collider (v-structure) on a path U if and only if the path U contains a subpath X → Z ← Y. For vertices X, Y and a vertex set Z which does not contain X and Y, X is d-connected to Y given Z if and only if there is an acyclic path U between X and Y such that every collider on U is either a member of Z or an ancestor of a member of Z, and no non-
Moreover, a set of variables V is causally sufficient if and only if no variable outside V is a direct cause of more than one variable in V. For a causally sufficient set of variables V with probability distribution P(V), the Markov condition assumes that the d-separation in the DAG G implies conditional independence in P(V), that is, if X is d-separated from Y by Z in G, then X is independent of Y conditional on Z in P(V); and the faithfulness assumption assumes that every conditional independence relationship in P(V) is entailed by the d-separation relationship for the causal DAG G, that is, if X is independent of Y conditional on Z, then X is d-separated from Y by Z. Therefore, the d-separation in the graph is equivalent to conditional independence in the distribution under the Markov condition and the faithfulness assumption. See Spirtes (2010) for a more detailed introduction and discussion. In this article, we always assume the causal sufficiency, the Markov condition, and the faithfulness assumption.
In the MPHIA dataset, classical constraint-based structural learning algorithms discover very few important features for Tri90, as shown later in Section 3. In this section, we explain the main cause of the problem and illustrate the motivation of our method with a simplified example from the MPHIA data. A complete analysis of the MPHIA data can be seen in Section 3.
Example 1. Suppose our graph contains only three vertices X, Y, and Z, and X is Tri90Aware, the indicator variable for awareness of the HIV positive status, which is one of the variables in which we are mainly interested in our Tri90 goal study. We want to check whether Y (AlcoholFrequency, an ordinal variable for alcohol drinking frequency) and Z (WealthQuintile, an ordinal variable for the wealthiness) are neighbors of X (Tri90Aware). That is to say, we are interested in whether the wealthiness and alcohol drinking frequency have direct causal relationships with the awareness of the HIV positive status. We have the following four (conditional) independence test results among X, Y, and Z for a confidence level α = 0.05 in the MPHIA dataset of males who are included in the Tri90 study.
X ⊥ Y is rejected, and X ⊥ Z is rejected, X ⊥ Y|Z is not rejected, and X ⊥ Z|Y is not rejected.
(1)
Since both Y (AlcoholFrequency) and Z (WealthQuintile) are not independent with X (Tri90Aware) marginally from the testing results, they should be connected to X (Tri90Aware) either directly or indirectly, according to the Markov condition. However, based on the faithfulness assumption, neither Y (AlcoholFrequency) nor Z (WealthQuintile) should be a neighbor of X (Tri90Aware) since X (Tri90Aware) is conditionally independent with Y (AlcoholFrequency) given Z (WealthQuintile) and X (Tri90Aware) is conditionally independent with Z (WealthQuintile) given Y (AlcoholFrequency). Those four testing results lead to contradicted conclusions under the Markov condition and the faithfulness assumption. Existing structural learning algorithms, such as the PC-stable algorithm, remove both the edge between Tri90Aware and AlcoholFrequency (X -Y) and the edge between Tri90Aware and WealthQuintile (X -Z) from the conditional independence testing results and conclude that there are no neighbors of X (Tri90Aware), which might be too strict in detecting edges.
Many covariates in the MPHIA data are categorical. When used as the conditional set, those categorical covariates lead to relatively high Type II error rates for the conditional independence tests. In addition, triples (X, Y, Z) with contradictory/inconsistent testing results such as Equation ( 1) are quite common in the MPHIA data. Such contradiction/inconsistency also leads to aggressive edge-removal for the constraint-based causal structural learning algorithms. To solve the false edgeremoval issue, we propose a new graphical structural learning algorithm, and we illustrate how the new algorithm successfully finds more edges later in Example 2.
In a directed acyclic graph (DAG) G with vertex set V and X ∈ V. Suppose N = (N 1 , N 2 , . . . , N q ) is the parents and children set of X, then under the faithfulness assumption and the Markov condition, respectively, we have
(2)
where M = {M 1 , M 2 , . . . , M r } = V\(N ∪ {X}), which is the set of vertices not connected to X. Furthermore, under the Markov condition and the faithfulness assumption, suppose a variable set N satisfies the conditions (2) and (3), then it is easy to verify that N is the set of parents and children of X.
For each vertex X, we want to get the best subset of V that fulfills the conditions (2) and (3). Our procedure consists of two steps: a forward step and a maximization step. The forward step finds all sets N that satisfy (2) for vertex X and will be illustrated in Algorithm 1 in Section 2.4. The maximization step picks the best set that fulfills (3) among those sets found by Algorithm 1 and will be provided in Section 2.5.
Before deriving the details of the algorithms, let us come back to Example 1 and see how algorithms based on (2) and (3) can solve the problem of aggressive edge removal that presents in the existing classical structural learning algorithms.
Example 2. Same as Example 1, suppose our graph only contains three vertices X (Tri90Aware), Y (AlcoholFrequency), and Z (WealthQuintile). Further, assume that we have the same (conditional) independence test results as in Example 1. Although the four testing results in (1) are incompatible under the Markov condition and the faithfulness assumption, it is possible to have a valid graphical structure with only three of them. For example, (4) and
(5)
Note that the neighborhood of X (Tri90Aware) satisfies (2) and (3) for X. In the situations of (4), the neighbor of X (Tri90Aware) is Z (WealthQuintile); and in the situations of (5), the neighbor of X (Tri90Aware) is Y (AlcoholFrequency). However, {Y, Z} does not satisfy (2) and is not the neighborhood of X. Algorithm 1 finds both {Z} and {Y} as potential neighborhoods for X (Tri90Aware) since both of them satisfy (2) for X. Algorithm 2 compares between {Z} and {Y} and chooses the set which fulfills (3) better as the neighborhood of X. In sum, the new algorithm chooses either {Z} or {Y} as the neighborhood of X, and concludes that either AlcoholFrequency (Y) or WealthQuintile (Z) is the neighbor of Tri90Aware (X) but not both of them. So, the proposed algorithm is less aggressive in edge removal than the existing classical structural learning algorithms such as the PC-stable algorithm as discussed in Example 1.
Before stating Algorithm 1, we first provide some useful definitions and their properties. Let p be the total number of vertices. For each vertex X, let T = V\{X} = {T i , i = 1, 2, . . . , p -1} and vertices are considered/added sequentially in the order of T 1 , T 2 , . . . , T p-1 to form a candidate neighborhood of X. That is to say, for any already formed non-empty candidate neighborhood
we consider whether an additional vertex from L X (S) := {T s q +1 , T s q +2 , . . . , T p-1 } can be added into the variable set, S. We also define L X (∅) = T, which means that we need to consider all the vertices in T when S starts from an empty set. Furthermore, let C X (S) be the vertices in L X (S) that can be added into S while still satisfying Equation (2). That is to say,
(6) Proposition 1 establishes some useful properties of C X (S), which are used in Algorithm 1 to facilitate the calculation of C X (S).
Proposition 1. C X (S) has the following properties.
2. If S does not satisfy (2), then
and
The properties of C X (S) in Proposition 1 can be shown by using its definition, and the proof can be found in Section S.1, supplementary materials. Proposition 1 establishes a recursive structure for C X (S). Equation ( 7) shows that C X (S) consists of all vertices that are marginally correlated with X when S is an empty set. Equation (10) states that C * X (S) is an upper-bound (concerning the partial order of inclusion) for C X (S) when S is non-empty. Furthermore, (11) tells us that we only need to examine the conditional independence relationship with conditional set of size |S| to get C X (S) from C * X (S) ∩ L X (S) for non-empty S. Equation ( 11) is used in Algorithm 1 to reduce the number of conditional independence tests and speed up the finding process. We can show that using (11), no conditional independence test is repeated in Algorithm 1. Indeed, let s max = max{s 0 , s 1 , . . . , s q }, S max = T s max , s -max = {s 0 , s 1 , . . . , s q }\{s max }, and S -max = {T s |s ∈ s -max }, then the conditional independence test of X and T s 0 given S = {T s 1 , T s 2 . . . , T s q } can only happen when we check whether S max can be added into candidate set S -max .
Step: Find all potential parent and children sets of X that satisfy (2). The algorithm sequentially adds the candidate set satisfying (2) into the preliminary result set, R. Require: (i) a vertex set V, (ii) a vertex X in V, and (iii) size α for conditional independence tests. Ensure: The set N X including all possible N that (approximately) satisfies Equation ( 2).
Let S = {∅} and R = ∅. while S = ∅ do Calculate C X (S) for each S ∈ S of cardinal from low to high as follows:
and N is not a proper subset of any other sets in R}.
Remark 1. Note that the final output of Algorithm 1, N X , only contains the candidate sets which are not proper subsets of other candidate sets for the following reasons. If N is a candidate set that satisfies (2) for X, then all proper subsets of N satisfy (2) for X, but no proper subset of N will satisfy (3). This is because for any N 1 N, take M ∈ N\N 1 , and we have X ⊥ M|S for any S ⊆ N 1 N since N satisfies (2). Hence, N 1 does not satisfy (3).
Remark 2. In the calculation of C X (S) and C * X (S), the size of conditional sets in the conditional independence test is not restricted. However, conditional independence tests with large conditional sets are not accurate, and the number of all possible conditional independence tests also grows exponentially with the size of conditional sets. Hence, we set the upper-bound of sizes of conditional sets in Algorithm 1, which is also discussed by Tsamardinos, Brown, and Aliferis (2006) and is commonly implemented in causal structural learning software such as Scutari (2010). In addition, as shown by ( 11), if we set the upperbound of sizes of conditional sets to be M CI , then we can approximate C X (S) by
Remark 3. M CI is a tuning parameter in Algorithm 1. Theoretically, one should choose M CI that is no less than the largest degree of vertices in the graph. However, there are concerns about using a large M CI , and the true degrees are also unknown in practice. First, a large M CI corresponds to big conditional sets, which will increase the computing cost. Second, when the sample size is not large enough, a large M CI leads to less reliable results due to the increased Type II errors. Generally, it is recommended to set M CI = 3 when the true graph is expected to be sparse or moderately sparse (Yan and Zhou 2020). We also carry on a simulation study on different choices of M CI , and more details and discussions can be seen in Section S.5.1, supplementary materials.
To choose the best neighbor set from the set of candidate sets N X , we check how well each candidate N ∈ N X satisfies Equation (3). Define
where CI(X, Y|S) is the p-value of some chosen conditional independence test for X and Y given S. To see why we use max in (13), note that conditional independence relationship implies big p-values, and S N (X, Y) measures whether any subset of N makes X and Y conditional independent. The maximum in ( 13) is from the fact that (3) only requires one S that makes X and Y conditionally independent given S as suggested by the large p-value. The idea of using conditional independence set of the largest p-value in the causal structural learning algorithm can also be found in Ramsey (2016). Define
which measures how well it is for subsets of N to "separate" X from any vertices not in N ∪ {X}. Here we are concerned about whether there is any M i ∈ V\(N ∪ {X}) violates the conditional independence between M i and X given N, and hence the minimum of p-values across M i 's is used. Equations ( 13) and ( 14) together can be seen as a minimax procedure. If N satisfies Equation (3), then for any M i not in the neighborhood of X, there should exist
) is large and it is likely that M i ⊥ X|S. So N with larger Q X (N) is likely to be the true neighbor set of X. Furthermore, we have the following useful properties of S N (X, Y) in Proposition 2, which are used in Algorithm 2 to facilitate the calculation of S N (X, Y).
Proposition 2. S N (X, Y) has the following properties.
1. Let ∅ be the empty set. Then
The properties of S N (X, Y) in Proposition 2 can be shown by using its definition, and the proof can be found in Section S.2, supplementary materials. Similar to Propositions 1 and 2 shows a recursive structure in S N (X, Y). (15) shows how to calculate S N (X, Y) for an empty set. ( 17) states that S * N (X, Y) is a lowerbound for S N (X, Y) for a nonempty set N. Furthermore, (18) tells us that we only need to calculate the p-value of conditional independence test CI(X, Y|N) to get S N (X, Y) from S * N (X, Y) when N is non-empty. ( 18) is used in Algorithm 2 to reduce the number of conditional independence tests and speed up the searching process. It is easy to prove that using (18), no conditional independence test is repeated in Algorithm 2 because the conditional independence test of X and Y given N only happens when we calculate S N (X, Y).
Remark. Similar to Algorithm 1, we propose to set upper-bound for sizes of conditional sets to be M CI , then S N (X, Y) can be approximated by
Furthermore, from (18), we can approximate S N (X, Y) by
Based on ( 18) and (20), we have Algorithm 2 which selects the N with the largest Q X (N) to be the neighbor set of X.
With Algorithms 1 and 2, we can learn the neighbor set of every vertex X in the DAG G and hence the skeleton of the DAG. We can further orient edges according to the conditional independence relationship. The overall structural learning algorithm is summarized in Algorithm 3. The orientation procedure is quite similar to those of classic constraint-based causal structural learning algorithms, such as the PC algorithm (Spirtes and Glymour 1991). We also incorporate some prior knowledge into Algorithm 2 Maximization Step: Find N ∈ N X with the largest Q X (N). Require: (i) a vertex set V, (ii) a vertex X in V, and (iii) candidate neighbor set N X of X. Ensure: N ∈ N X with the largest Q X (N).
Let r = 0 and N
. end if end for return N X as the output. the orientations. For example, we know covariates such as Age-Group cannot be affected by other covariates such as Education, so if there is an edge between AgeGroup and Education, then we orient the edge as AgeGroup to Education. If there are multiple d-separation sets S(X, Y) for a nonadjacent pair (X, Y), we shall use the d-separation set S(X, Y) with the largest p-value to make the orientation results stable. For Z adjacent to both X and Y, we call X -Z-Y a v-structure and make the orientation X → Z ← Y if Z ∈ S(X, Y). Finally, note that there could be conflicting vstructures. For example, if there is a chain X -Z -Y -W, and both
In Algorithm 3, we resolve the conflict between v-structures by comparing the p-value of the v-structure. Going back to the previous example, if S(X, Y) has a larger p-value than S(Z, W), then we orient Z -Y as Z ← Y following v-structure X -Z -Y and vice-versa.
Remark 1. Note that multiple orientations may satisfy the inferred d-separation (conditional independence) structure. Hence, in Algorithm 3, we first use prior knowledge in the orientation process to establish orientations and to enhance the interpretability of the orientation result. For instance, we presume that age may lead to education status but not the other way around. The prior knowledge is provided in the sample code of the supplementary materials.
Remark 2.
Form G by recursive orientation according to the following two rules: 1. If X-Y and there is a directed path from X to Y, then orient X -Y as X → Y; 2. If X and Y are not adjacent and there is a Z such that X → Z and Z -Y, then orient Z -Y as Z → Y.
the largest p-value to stabilize the orientation results. It agrees with Equation (13) used in the maximization step for skeleton learning. The approach of using the conditional independence set with the largest p-value in the orientation process has also been used in Ramsey (2016). Sections 2.5 and 2.6 both use the largest p-value among the conditional independence tests but for different purposes. Section 2.5 is about skeleton learning, and Section 2.6 is about orientation. Ramsey (2016) only concerned the orientation phase but not the skeleton learning.
Remark 3. To make the graph learned by the proposed algorithm more interpretable, we calculate p-values for the significance of the (undirectional) connections for all edges in the graph learned by Algorithm 3. More specifically, for edge X -Y between the covariates X and Y, define P(X, Y) to measure the significance of X -Y as follows
where N X and N Y are the neighbors of X and Y, respectively. Note that the measure P(X, Y) is undirectional, that is, the orientation of the edge X -Y has no effect on P(X, Y) from its definition, and P(X, Y) also has no information on the edge orientation.
The MPHIA survey is a new HIV-focused, cross-sectional, household-based, nationally representative survey of adults and adolescents aged 15 years and older as well as children aged 0-14 years. In addition to HIV testing results, the survey contains demographic questions, such as age group, gender (preference to the gender collected in the adult or adolescent questionnaire), ethnic group, and HIV-related questions, such as access to preventive care and treatment services. There were 26,871 survey participants and 1407 covariates in total. We will focus on adults and adolescents aged 15 years and older with complete Tri90 related information. There are 2217 such individuals in the MPHIA survey data, including 712 males and 1505 females. That is to say, these 2217 individuals are HIV-positive with known status of HIV awareness (Aware), antiretroviral therapy (ART), and viral load suppression (VLS).
In PHIA surveys, the 1407 covariates include 34 continuous variables such as age and time, 32 discrete variables such as the number of partners, and 1341 categorical variables (nominal ones such as gender and ethnic group, ordinal ones such as alcohol frequency). Some covariates are applicable to males or females only, that is, the questions related to pregnancy are only applicable to females. Therefore, the numbers of vertices are quite different between the female graph and the male graph. After a data-preprocessing procedure (illustrated later), we will apply the causal structural learning algorithm to six datasets for each combination of gender and 90-90-90 goal separately. Overview of the data-preprocessing procedure: We first drop all the covariates with a dominant level (one level has > 99% samples), because the MPHIA sample size is not large enough to reject any null hypothesis of the conditional independence involving those variables. Also, we "merge" all the closely related covariates into a single covariate in the MPHIA dataset. For example, some multi-option questions are dummy-coded by many two-option ones in the MPHIA dataset, and we combine them together to create multi-level categorical covariates. By "merging" these kinds of covariates, we can reduce the number of covariates and improve the interpretability of our results. Also, some covariates appear in the MPHIA dataset multiple times with the same meaning but different names, such as Eth-nicGroup and EthnicCode, and we keep one covariate and drop the others in such situations. We further remove some covariates that are direct indicators of the Tri90 goals from the MPHIA data. For example, there are questions like whether the subject takes a certain ART medicine or not in the MPHIA survey. These covariates are strongly correlated to the ART status, but not helpful for the purpose of building the causal pathways because they may block the connection between the ART status and other meaningful covariates. Furthermore, some categorical variables have many levels which can complicate the analysis. To reduce the number of levels while keeping the main information of each covariate, we keep the biggest levels of each covariate which cover at least 95% of the individuals, and combine the remaining levels into the "Others" category. We subsample the whole dataset by gender and by Tri90 goal to create six datasets, and within each of the six datasets, we further drop the covariates with a dominant level (one level has > 99% samples). The sample sizes and the numbers of vertices in each dataset are summarized in Table 1. There is one covariate included in the female awareness dataset but not in the female ART and VLS datasets: LiveHere (whether the individual lives here or not).
The covariate is categorical and has a dominating level, which is below the threshold of 99% in the female awareness dataset but is above the threshold in the female ART and VLS datasets. More details about the six datasets are provided in Section S.3, supplementary materials.
In the next section, we compare the results obtained from the proposed algorithm with those of the existing algorithms in Section 3.1. We further discuss the potential Tri90 pathways discovered by the proposed algorithm in detail in Section 3.2.
We also apply the following existing structural learning algorithms to each dataset: PC-stable algorithm proposed by Colombo and Maathuis (2014), which is a stable/orderindependent variant of the original PC (initials of the first names) algorithm proposed by Spirtes and Glymour (1991); MMPC (Max Min Parents and Children) algorithm proposed by Tsamardinos, Aliferis, and Statnikov (2003b); IAMB (Incremental Association Markov Blanket) algorithm proposed by Tsamardinos, Aliferis, and Statnikov (2003a); GS (Grow-Shrink) algorithm proposed by Margaritis (2003). Notice that we use the PC-stable algorithm instead of the PC algorithm since the PC-stable algorithm is order-independent. Order independence means that a random reordering of the variables does not affect the graphical learning result, which enhances the stability of the algorithm and also makes the results more interpretable.
We use the existing graphical learning algorithms implemented in the R package bnlearn (Scutari 2010) in Sections 3 and 4. The default conditional independence test statistic in the bnlearn package is the mutual information for categorical variables and the linear correlation for continuous variables. We set the upper bound of sizes of conditional sets M CI = 3 and the size of the conditional independence test α = 0.05 by default for all the causal structural learning algorithms used in the comparison.
We first summarize the structural learning results by different graphical learning methods. The results in Table 2 shows that the proposed algorithm makes new discoveries regarding Tri90 goals. More specifically, Table 2 show that the numbers of edges (NE) and the number of directed edges (NDE) are both much smaller than the number of vertices (NV) in the graphs learned by the existing algorithms including the PC-stable, the MMPC, the IAMB, and the GS algorithms. These existing algorithms lean toward fractured graphs and do not have much conditional independence information for the orientation. On the contrary, the proposed algorithm produces larger numbers of edges and well-connected graphs and has more potential to infer the directions of edges.
Remark. It is important to note that the goal of the causal graphical algorithms is not to produce as many edges and directed edges as possible. Later in this section, we will use the Bayesian Information Criterion (BIC) to show that these discoveries made by the proposed algorithm provide useful information about the MPHIA data. Furthermore, in Section 3.2, we discuss the discovered pathways in detail, which are reasonable, also confirmed in other HIV Tri90 literature, and can provide useful insight for the three Tri90 goals.
Let d be the distance (defined as the length of shortest path regardless of direction) from a particular 90-90-90 goal (awareness of HIV, ART, or VLS) to a covariate and N(d ≤ k) be the number of covariates whose distances to a 90-90-90 goal are smaller than or equal to k. Very few covariates are close to 90-90-90 goals in the graphs learned by the existing algorithms with only two covariates whose distance to the 90-90-90 goals are smaller or equal to three. On the contrary, the proposed algorithm discovers many covariates that are of a small distance to the 90-90-90 goals including several direct neighbors (d ≤ 1), which includes the ones discovered by the existing algorithms. See details in Table S.1, supplementary materials.
The Bayesian information criterion (BIC, Schwarz 1978) is a classical statistical tool for model selection. We compare our proposed graphical learning algorithm with the aforementioned classical PC-stable, MMPC, IAMB, and GS algorithms using BIC criterion and summarize the results in Table 3. Let D k , k = 1, 2, . . . , 6, be the datasets corresponding to the three 90-90-90 goals of each gender, respectively. For k = 1, 2, . . . , 6, we use each of the graphical learning algorithm A i , for i = 1, 2, . . . , 5, to learn a DAG G i,k on D k , where A i , i = 1, 2, . . . , 5, stand for the PC-stable, MMPC, IAMB, GS, and the proposed algorithm, respectively.
The log-likelihood of a DAG G can be decomposed as follows:
where θ j 's are parameters of the model, n is the sample size, p is the number of covariates,
is the set of parents of X j in the DAG G, and π ij (G) is the ith observation of π j (G).
It shows that the log-likelihood of G can be decomposed as the sum of the log-likelihood of local structures of a covariate given its parents. For the estimation of the local structures, existing literature on causal graphical models often assumes a linear model of a covariate on its parents (Spirtes 2010;Valente et al. 2010;Bolla, Abdelkhalek, and Baranyi 2019). Note that in the MPHIA dataset, many covariates are categorical, so we assume generalized linear models (GLM) of a covariate given its parents instead of linear models, and we fit the GLM of local structures by MLE. The degree of freedom (DF) and the log-likelihood of a DAG are the sums of DFs and log-likelihoods of the local structures, respectively. Note that all the Log-likelihood* in Table 3 are positive since they are the differences between the log-likelihood and the log-likelihood of the null model (model with intercept only). Furthermore, in Table 3, BIC score is calculated by -2Log-likelihood* + DF log(n), where n is the sample size. Hence, lower BIC scores correspond to better models.
Our proposed algorithm learns a much larger number of edges in all the six graphs compared with the existing algorithms in Table 2, and thus it has much larger degrees of freedom for the log-likelihood defined in Equation ( 22). Table 3 shows that the proposed algorithm has the largest degree of freedom and the largest log-likelihood. The much larger log-likelihoods imply that the proposed algorithm discovers a lot more useful information in MPHIA, which is further confirmed by the best (smallest) BIC scores in Table 3.
Table 4 lists the neighbors of the 90-90-90 goals discovered by the proposed causal structural learning algorithm. Details of the graphs learned by the proposed algorithm are provided in Figures 1 andS.1-S.5, which render parts of the graphs surrounding each 90-90-90 goal for each gender. The partial graphs we present in Figures 1 andS.1-S.5 are uniquely determined using our proposed algorithm with prior knowledge in the orientation process. Hence, the problem of multiple orientations does not affect the interpretation of the results of the MPHIA data. See Remark 1 of Algorithm 3. Also note that the covariates are renamed to have more intuitive meaning than their original names in MPHIA data codebooks. The MPHIA data codebook and the meaning of important covariates are provided in Section S.7, supplementary materials. Those covariates discovered by the proposed algorithm are consistent with findings in the existing literature, and we discuss them by Tri90 goals and by genders. However, we should be cautious about the interpretation of the results of the proposed algorithm. Notably, the proposed algorithm assumes causal sufficiency, implying no hidden confounding (or common causes) covariates. Although the MPHIA survey has included 1407 related covariates, and most of the confounding covariates are likely observed, there may still be some unobserved confounding factors that can affect the interpretation of the proposed algorithm results.
For HIV awareness among females, we use Figure 1 to visually illustrate the covariates that are closely connected to HIV awareness (distances ≤ 2). The learned graphs for other Tri90 goals and genders are provided in Figures S.1-S.5 in the Section S.4, supplementary materials. We find that HIV-positive females with more pregnant times are more likely to know their HIV status. We think this is because pregnant women are more likely to be tested for HIV during antenatal clinic visits or laboring. Peltzer et al. (2009) reported that age, condom usage, and education are associated with HIV status awareness. The proposed algorithm finds that AgeGroup and Education are potential "reasons" for female's HIV awareness, and that EasyCondom and PLWHSupport-Group may be results of female's HIV awareness: (i) older individuals are more likely to know their HIV status; (ii) individuals who are unaware of their HIV status often do not know whether it is easy to get condoms or have difficulties in getting condoms; (iii) females unaware of their status never answer the PLWH-SupportGroup question; (iv) individuals with higher education are less likely to be aware of their HIV positive status, which is opposite of the marginal association reported in Peltzer et al. (2009). It is because (iv) is the conditional association given AgeGroup (which affects both Education and HIV awareness): while females with higher education levels in the age group 35-44 are more likely to be aware of HIV positive status, females with higher education levels in other age groups are less likely to be aware of the HIV status. Further studies will be needed to establish the relationship between education and HIV status awareness.
Figure 1 also presents other covariates that indirectly connect to the HIV awareness among female HIV patients. Among those, employment status (WorkLast12Mo) and marital status are impacted by both urban residence and education; and education is impacted by ethnic group, urban residence, and age group. AgeGroup is a potential "reason" for SellSexEver and SyphilisTestInPreg: younger generations are less likely to have SyphilisTestInPreg; older generations are less likely to have sold sex ever. CircumcisedHIVRisk is a score measuring whether the individual agrees that men who are circumcised are not at risk of HIV at all, do not need to use condoms, and can have multiple sexual partners without the risk of HIV. The graph suggests females that are more cautious toward HIV risk (lower CircumcisedHIVRisk score) have smaller numbers of pregnancies on average. Finally, more sexual partners (PartNumber12Mo), being a partner to the head of house (RelationToHeadOfHouse), and lower tolerance toward violence (ViolenceOK) lead to easiness in getting condoms.
We find males whose partners are older are more likely to be aware of their HIV status, both marginally, or conditioning on their age groups. Comparing the male awareness pathway and the female awareness pathway, we suspect that older females contribute more to HIV awareness because someone's age group is associated with his/her partner's age.
In Figure S.1, the proposed algorithm finds that Partner-Number12Mo (sexual partners in the last 12 months) and WifeNumLiveElsewhere (number of wive/partners that live elsewhere) to be important covariates for male HIV awareness. These discoveries are supported by Peltzer et al. (2009) which showed that HIV awareness is significantly associated with the number of sexual partners. WifeNumLiveElsewhere and PartnerNum12Mo are strongly correlated and have closely related meanings. Instead of studying their impacts on HIV awareness separately, we study the relationship between HIV awareness and WealthQuintile which is the parent of both WifeNumLiveElsewhere and PartnerNum12Mo. We find that "poor" HIV positive males (with wealth in lower 20%) are more likely to have no partner in the last 12 months, more likely to have wives/partners live elsewhere, and less likely to be aware of their HIV status than the wealthier (in the upper 80%) HIV positive males. This finding is consistent with Dokubo et al. (2014).
We also find that ViolenceOK and AlcoholFrequency may be possible "reasons" of HIV awareness in males, where Vio-lenceOK is a score to measure whether the individual believes it is right for a man to beat his wife/partner under various scenarios. Further investigation shows that among males who are HIV positive: (i) those with a higher violence score are more likely to be unaware of their HIV status; (ii) those who never drink are most likely to be aware of their HIV status while those with a high frequency of alcohol use (more than four drinks a week) are least likely to be aware of their HIV positive status, and such an effect is more significant in males who have not worked in the last 12 months. Alcohol usage is an important factor associated with sexual risk behavior (Kalichman et al. 2007;Hahn, Woolf-King, and Muyindike 2011), and violence score can also be an important factor for sexual behavior. These factors deserve more attention when advocating HIV awareness.
The proposed causal structural learning algorithm also finds that WorkLast12Mo has a small distance to HIV awareness which is supported by Peltzer et al. (2009). PLWHSupportGroup is a result of HIV awareness because males unaware of their status never answer the PLWHSupportGroup question.
Next, we discuss the ART pathways for females and males in one section because they are relatively simpler compared with the HIV awareness pathways. For the same reason, we discuss the VLS pathways for females and males in the same section. Hodgson et al. (2014), Tomori et al. (2014), and Kebaabetswe et al. (2019) showed the connection between TravelTime to the receiving of ART. The proposed algorithm indicates that among people living with HIV and being aware of their positive HIV status, those who do not know the travel time to HIV care are less likely to receive ART treatment than the individuals who know the travel time. Friedman et al. (2015), Hatcher et al. (2015), and Sullivan, Messer, and Quinlivan (2015) showed that violence and other kinds of abuse are significantly associated with ART initiation. We find females with smaller violence scores are more likely to be on ART treatment.
We find that wealthier males are more likely to initiate the ART treatment, which is in agreement with Gebru, Lentiro, and Jemal (2018). On the other hand, the males with "poor" wealth are likely to have no partners in the last 12 months, and more likely to have wives or partners live elsewhere. The proposed algorithm also discovers that females with known status of SyphilisTestInPreg are more likely to be on ART treatment. The learned female ART pathway reveals that those received ART treatments had easy access to the HIV care (short travel time) and easy access to antenatal care where they could be offered the syphilis test during pregnancy. They were both indicators of receiving good health care services which increased the chance of receiving ART treatment.
Finally, we find AbnormPenisDischarge to be a potential "reason" for male ART initiation. Further investigation shows that males with the missing response on abnormal penis discharge problems are more likely to be on ART treatment among the HIV-positive males who are aware of their HIV status. Unfortunately, it was unclear why people did not respond AbnormPenisDischarge question.
We find females who attend the support groups more frequently are more likely to have their viral load suppressed. It is supported by Roberts (2000), Tomori et al. (2014), and Rangarajan et al. (2016). Friedman et al. (2015), Hatcher et al. (2015), and Sullivan, Messer, and Quinlivan (2015) showed that violence and other kinds of abuse are significantly associated with viral load suppression (VLS), especially among females. It supports the connection between VLS and ForceSexTimes found in the female VLS pathway.
We also find SeekMedicalHelp and WifeNum to be potential "reasons" of VLS among males who are on ART. Males who seek help from doctors or nurses because of health issues such as abnormal penis discharge and painful urination are more likely to have viral load unsuppressed; and males with more wives are more likely to have viral load suppressed.
In Section 3, we see that the proposed causal structural learning algorithm discovers many more new edges than the existing ones, and we validate the results of the proposed algorithm through the Bayesian information criterion (BIC). However, since the true causal graph for MPHIA data is unknown, we cannot validate the edge discoveries directly. To get more insights on the Type I and Type II error rates and their tradeoff for the proposed algorithm against existing ones, we carry out simulation studies with settings mimicking the MPHIA data in this section.
We check the true positive and negative rates of the proposed algorithm against the existing algorithms on synthetic datasets that mimic the MPHIA data. More specifically, we use the DAG learned from MPHIA data as the truth to generate the simulation data. Here we choose the graphs learned by the proposed algorithm as the truth for simulation purposes since Table 3 shows that the graphs learned by the proposed algorithm are better fits for the MPHIA data in the BIC criterion than those learned by the other algorithms. That is to say, let G k be the DAG learned by the proposed algorithm on the 90-90-90 MPHIA dataset D k for k = 1, 2, . . . , 6. Then we fit the data distribution P k based on G k on the data D k . We further randomly generate simulated datasets based on the distribution P k with sample size n. The distribution estimation and the random datasets generation use the bn.fit (Bayesian network fitting) and rbn (random Bayesian network) functions in the R package bnlearn (Scutari 2010). We then carry out the proposed algorithm together with the aforementioned PC-stable, GS, MMPC, and IAMB algorithms on the generated datasets. Finally, we calculate true positive rates and true negative rates of edges disregarding the orientation for each algorithm.
Here we set the sample size n = 250, 500, 1000 for a sample size similar with our real data and to check the performance of the proposed algorithm with different sample sizes. We repeat the Monte Carlo simulation 500 times for each setting and summarize the results in Table 5. The left and right panels of Table 5 summarize the empirical true positive and negative rates of the proposed algorithm as well as those of existing algorithms, respectively. From the right panel of Table 5, we can see that the proposed algorithm has similar true negative rates with existing algorithms. Furthermore, from the left panel of Table 5, we can see that the proposed algorithm has better true positive rates than existing algorithms.
UNAIDS 90-90-90 goals are important milestones to end AIDS.
To understand the progress on the 90-90-90 goals better, we analyze the Malawi PHIA (MPHIA) dataset to discover important covariates and potential causal pathways for the 90-90-90 goals through causal structural learning in the paper. Existing classical constraint-based causal structural learning algorithms are quite aggressive in edge removal and can lead to information losses while building directed graphical models, especially in the case of categorical variables and relatively small sample sizes. To deal with the problem, we propose a new causal structural learning algorithm. The proposed algorithm can preserve more information about important features and potential causal pathways as shown by various numerical studies when many covariates in the domain are categorical. More specifically, the proposed algorithm improves true positive rates over the existing classical algorithms while having a comparable true negative rate. It shows that our proposed algorithm has a great potential to discover important features and potential causal pathways, especially in a domain with many categorical variables.
Carrying out the causal graphical analysis on the MPHIA dataset, we obtain interesting results on important covariates and possible causal pathways related to the UNAIDS 90-90-90 goals. For example, the proposed algorithm discovers age and condom usage to be important covariates for female HIV awareness and number of sexual partners to be important for male HIV awareness, which agrees with literature, such as Peltzer et al. (2009). The proposed algorithm also discovers similarities as well as differences between female and male pathways. For example, travel time is discovered to be an important covariate for both female and male ART. However, there are also different important covariates for female and male ART, such as attitude toward violence for female ART and partner numbers for male ART.
It is also important to pay attention to the assumptions of the proposed causal structural learning algorithm when we interpret the results. One important assumption behind the proposed algorithm is causal sufficiency which is critical for our algorithm and many other constraint-based causal structural learning algorithms. Although the MPHIA survey provides many related covariates, the causal sufficiency assumption may still not hold perfectly.
In the article, we stratify the MPHIA dataset by sex and learn the graphical models for each sex and each Tri90 goal. It is also possible to learn the DAG with further stratified data by other covariates such as age, but the sample size would be too small for each stratum to make reliable inferences.
The discoveries on causality are important extensions for existing literature where only correlation instead of causation is established. More studies can be carried out to further validate the potential causal discoveries, and other statistical inference tools such as mediation analysis can be applied to further the understanding of the causal relationship. The discoveries on causality can help develop better HIV response strategies and related policies.
The supplementary materials contain the proofs of Propositions 1 and 2. The supplementary materials also present additional details for the six Tri90 datasets. The supplementary materials also present some analysis results of the MPHIA data, more simulation studies, as well as the code to carry out the analysis in the paper. A summary for the notations in the proposed algorithms and parts of MPHIA codebook are also provided in the supplementary materials.