Identification and Estimation of the Bi-Directional MR with Some Invalid Instruments

Authors: feng  xie
zhen  yao
lin  xie
yan  zeng
zhi  geng
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 04:23:36

Abstract: We consider the challenging problem of estimating causal effects from purely observational data in the bi-directional Mendelian randomization (MR), where some invalid instruments, as well as unmeasured confounding, usually exist. 
To address this problem, most existing methods attempt to find proper valid instrumental variables (IVs) for the target causal effect by expert knowledge or by assuming that the causal model is a one-directional MR model. 
As such, in this paper, we first theoretically investigate the identification of the bi-directional MR from observational data. In particular, we provide necessary and sufficient conditions under which valid IV sets are correctly identified such that the bi-directional MR model is identifiable, including the causal directions of a pair of phenotypes (i.e., the treatment and outcome).
Moreover, based on the identification theory, we develop a cluster fusion-like method to discover valid IV sets and estimate the causal effects of interest.
We theoretically demonstrate the correctness of the proposed algorithm.
Experimental results show the effectiveness of our method for estimating causal effects in both one-directional and bi-directional MR models.

Main paper: 
 introduction: 
 Mendelian randomization (MR) is a powerful method for estimating the causal effect of a potentially modifiable risk factor (treatment) on disease (outcome) in observational studies, which has been used across a wide variety of contexts, such as clinic [Mokry et al., 2015 , Carreras-Torres et al., 2017] , socioeconomics [Cuellar-Partida et al., 2016 , Brumpton et al., 2020] , and drug targets [Ference et al., 2015 , Li et al., 2020] . Most existing MR methods assume a one-directional causal relationship. However, bi-directional relationships are ubiquitous in real-life scenarios, for example, obesity and vitamin D status [Vimaleswaran et al., 2013] , body mass index (BMI) and type 2 diabetes (T2D), etc. Understanding and analyzing bi-directional relationships can provide deeper insights into complex systems, leading to more effective interventions and solutions.Within a causal inference framework, MR can be implemented as a form of instrumental variables (IVs) analysis in which genetic variants 2 serve as IVs for the modifiable risk factor. Specifically, as shown in Figure 1 , suppose obesity (X) and Vitamin D status (Y ) are the treatment and outcome of interest, respectively, and lifestyle factors (U) are the unmeasured confounders between them. One can use genetic variants G as IVs to explore the causal effect of X on Y in observational data [Vimaleswaran et al., 2013] if G are valid IVs, i.e., G satisfy the following three assumptions [Wright, 1928 , Goldberger, 1972 , Bowden and Turkington, 1990 , Hernán and Robins, 2006 ]:A1.[Relevance] The genetic variants G are associated with the exposure X; A2. [Exclusion Restriction] The genetic variants G have no direct pathway to the outcome Y ; Figure 1 illustrates these three assumptions for a genetic variant G that is a valid IV relative to the causal relationship X → Y . Given a valid IV G in the linear model, one can estimate the causal effect of risk factor X on the outcome Y of interest consistently by using Two-Stage Least Squares (TSLS) [Burgess et al., 2017] . However, in most MR studies, the genetic variants may include some invalid IVs due to horizontal pleiotropy-where some genetic variants may have pleiotropic effects on both the treatment X and the outcome Y via different biological pathways [Davey Smith and Ebrahim, 2003, Lawlor et al., 2008] . How to select valid IVs from such genetic variants becomes the major issue in the MR study. In general, valid IVs need to be chosen based on domain knowledge or expert experience [Richmond et al., 2017 , Martens et al., 2006 , Zhao et al., 2023] . However, it is usually difficult to select appropriate valid IVs without precise prior knowledge of genetic variants, and an invalid IV may cause biased estimation of the effect of X on Y [Sanderson et al., 2022, Richmond and Smith, 2022] . Therefore, it is desirable to investigate ways of selecting genetic variants as valid IVs solely from observed variables.Many contributions have been made to handle MR studies with some invalid IVs under specific assumptions. Existing methods can be roughly divided into four strategies. The first strategy assumes the InSIDE assumption, which states that the genetic variants' pleiotropic effects on the outcome Y are uncorrelated with their effects on the exposure X, and requires that the number of IVs increases with the sample size [Bowden et al., 2015] . The second strategy is based on the Majority rule assumption, which assumes that over 50% of the candidate IVs are valid [Kang et al., 2016 , Bowden et al., 2016 , Windmeijer et al., 2019 , Hartford et al., 2021] . The third strategy assumes the Plurality rule assumption, which means that the number of valid IVs is larger than any number of invalid IVs with the same ratio estimator limit [Guo et al., 2018 , Windmeijer et al., 2021] . The last strategy assumes the Two-Valid IVs assumption, which requires that the number of valid IVs is greater than or equal to 2, and the Rank-faithfulness assumption, which states that the probability distribution P is rank-faithful to a directed acyclic graph G if every rank constraint on a sub-covariance matrix that holds in P is entailed by every linear structural model with respect to G [Silva and Shimizu, 2017 , Xie et al., 2022 , Cheng et al., 2023] . Although these methods have been used in various fields, they assumed that the one-directional causal relationship between X and Y is known and may fail to handle scenarios where bi-directional relationships are present.Recently, Darrous et al. [2021] proposed a latent heritable confounder MR method to estimate bidirectional causal effects with a hierarchical prior on the true parameters. However, their method is sensitive to prior misspecification. Later, Li and Ye [2022] introduced a focusing framework to test bi-directional causal effects under the noise distribution of the effect estimates assumption, which is exciting. However, this assumption is strictly parameterized (See Condition 1 in [Li and Ye, 2022] for more details). Moreover, they focus on the two-sample MR model that assumes homogeneity across samples, whereas we consider the one-sample MR model. See Appendix A for more details.measured genetic variants, which may include invalid IVs. Analogous to Li and Ye [2022] , for each sample, indexed by i, the phenotype X i is taken as a linear function of the phenotype Y i , genetic variants G i , and an error term ε Xi , while phenotype Y i is modeled as a linear function of the phenotype X i , genetic variants G i , and an error term ε Yi . Without loss of generality, we assume that all variables have a zero mean (otherwise can be centered). Specifically, the generating process of the data is given below:EQUATIONwhere β Y →X is the causal effect of Y on X and β X→Y is that of X on Y . γ X = (γ X,1 , ..., γ X,g ) ⊺ and γ Y = (γ Y,1 , ..., γ Y,g ) ⊺ are the direct effects of genetic variants G on X and Y , respectively. Note that ε X and ε Y are error terms that are correlated due to unmeasured confounders U between X and Y . Besides, there may exist a single genetic variant G i ∈ G with both γ X,i ̸ = 0 and γ Y,i ̸ = 0, which implies that G i also has a direct pathway to the outcome, violating Assumption A2 [Exclusion Restriction] . Like Li and Ye [2022] , we here restrict the independence between any two genetic variants in G (i.e., genetic variants are randomized). In Section 5, we discuss the situation where this restriction is violated. Unlike Li and Ye [2022] , we allow genetic variants to cause unmeasured confounders (partially violating Assumption A3 [Randomness]).Following Hausman [1983] , we assume that 1) can be reorganized to avoid recursive formulations as follows:β X→Y β Y →X ̸ = 1. Then let ∆ = 1/(1 -β X→Y β Y →X ), Eq.(X = (G ⊺ γ X + G ⊺ γ Y β Y →X + ε X + ε Y β Y →X )∆, Y = (G ⊺ γ X β X→Y + G ⊺ γ Y + ε X β X→Y + ε Y )∆.(2)If β X→Y β Y →X = 1, the causal effects β X→Y and β Y →X are not identifiable, even in the presence of a valid IV.Goal. In this paper, we aim to address whether the bi-directional MR model in Eq.( 1) is identifiable; that is, whether we can estimate the causal effects β X→Y and β Y →X given infinite data, even without prior knowledge about which candidate IVs are valid or invalid. 
iv estimator: 
 In this section, we briefly describe the classical IV estimator, the two-stage least squares (TSLS) estimator, which is capable of consistently estimating the causal effects β X→Y and β Y →X of interest in the above bi-directional MR causal model of Eq.(1) when the valid IVs are known [Wooldridge, 2010] . We denote by G X→Y Proposition 1 (Two Stage Least Square (TSLS) Estimator). Assume the system is a linear bidirectional causal model 1. For a given causal relationship X → Y in the system, the causal effect of X on Y can be identified byEQUATION)where P = (G X→Y V ) ⊺ G X→Y V (G X→Y V ) ⊺ -1 G X→Y V is the projection matrix.Remark 1. Note that, when utilizing the TSLS estimator, the causal effect is biased if the candidate set G includes invalid IVs. For instance, consider the causal relationship X → Y , given an invalid IV setG X→Y I , the causal effect of X on Y is given by βX→Y = X ⊺ PX -1 X ⊺ PY = β X→Y + X ⊺ PX -1 X ⊺ P(G ⊺ γ Y + ε Y ) β bias , where P = (G X→Y I ) ⊺ G X→Y I (G X→Y I ) ⊺ -1 G X→Y I .In the remaining of the paper, we denote such a procedure by βX→Y = TSLS(X, Y, G X→Y ).Similarly, for causal relationship Y → X, one can also have the corresponding causal effect of Y on X given the set of IVs, i.e., βY →X = TSLS(Y, X, G Y →X ). For details about the proof of Proposition 1 and Remark 1, please refer to Appendix G.1 and G.2. 
identifiability of bi-directional mr model: 
 In this section, we first show by a simple example that the valid and invalid IV sets can impose different constraints. Then, we formulate these constraints and present the necessary and sufficient conditions that render the model identifiable.U G1 G2 G3 G4 G5 X YFigure 2: An illustrative example where valid and invalid IV sets induce distinct constraints, whereG X→Y V = (G 1 , G 3 ) ⊺ is a valid IV set, while G X→Y I = (G 2 , G 4 , G 5 ) ⊺ is invalid due to pathways G 2 → Y , G 4 → Y and G 5 → Y .Figure 2 gives an illustration of our basic idea. Suppose there are 2 valid IVs, i.e., G 1 and G 3 , and 3 invalid IVs, i.e., G 2 , G 4 and G 5 , relative to the causal relationship X → Y . For the inverse Y → X, identical conclusions can be derived. Assume the generating process of X and Y follows Eq.( 1), where specifically we set all variables to have a zero mean and unit variance andX = Y β Y →X +G 1 γ X,1 +G 2 γ X,2 +G 3 γ X,3 +G 4 γ X,4 + G 5 γ X,5 + ε X , and Y = Xβ X→Y + G 2 γ Y,2 + G 4 γ Y,4 + G 5 γ Y,5 + ε Y , where β X→Y = β Y →X = 0.6, γ X,1 = γ X,2 = γ X,3 = 1, γ X,4 = 1.8, γ X,5 = 1.2, γ Y,2 = 0.5, γ Y,4 = 0.6, γ Y,5 = 0.3.Observation 1. Given a valid IV set where there are at least two IVs , e.g.,G X→Y V = (G 1 , G 3 ) ⊺ , we have: corr(Y -Xω {G3} , G 1 ) = 0, corr(Y -Xω {G1} , G 3 ) = 0,where corr(•) denotes the Pearson's correlation coefficient between two random variables, andω {Gi} = TSLS(X, Y, {G i }) with i ∈ {1, 3}.Observation 2. However, given an invalid IV set, e,g.,G X→Y I = (G 2 , G 4 , G 5 ) ⊺ , we have: corr(Y -Xω {G4,G5} , G 2 ) ̸ = 0, corr(Y -Xω {G2,G5} , G 4 ) ̸ = 0, corr(Y -Xω {G2,G4} , G 5 ) ̸ = 0,where ω {Gi,Gj } = TSLS(X, Y, {G i , G j }) with i ̸ = j and i, j ∈ {2, 4, 5}.Based on these two observations, we can see that the valid and invalid IV sets entail clearly different constraints regarding the correlation coefficients3 , which inspires us to further derive our identification analysis. For details about the observations, please refer to Appendix C.We begin with the following definition and the main assumptions. Definition 1. (Pseudo-Residual) Let G be a subset of candidate genetic variants. A pseudo-residual of {X, Y } relative to G is defined asEQUATION)where ω G = TSLS(X, Y, G).This "pseudo-residual" is applied to measure the uncorrelated relationship with some genetic variants. Note that concepts to "pseudo-residual" have been developed to address different tasks [Drton and Richardson, 2004 , Chen et al., 2017 , Cai et al., 2019 , Xie et al., 2020] , but our formalization is different from theirs in terms of the parameter ω G . To the best of our knowledge, it has not been realized that the uncorrelated property involving such pseudo-residuals reflects the validity of the IVs in the bi-directional MR causal model (see Propositions 2 ∼ 3).As analyzed by Chu et al. [2001] , within the framework of linear models, a variable being a valid IV imposes no constraints on the joint marginal distribution of the observed variables. In other words, there is no available test to determine whether a variable is a valid IV without making further assumptions. Therefore, we introduce the following assumption.Assumption 1 (Valid IV Set). For a given causal relationship (if the relationship exists), there exists a valid IV set that consists of at least two valid IVs. For example, for the causal relationshipX → Y , |G X→Y V | ≥ 2.Assumption 1 is the same as in [Silva and Shimizu, 2017, Cheng et al., 2023] if the model is a one-directional MR model. Notice that this assumption is much milder than the Majority rule assumption [Kang et al., 2016 , Bowden et al., 2016 , Windmeijer et al., 2019 , Hartford et al., 2021] and Plurality rule assumption [Guo et al., 2018 , Windmeijer et al., 2021] : Assumption 1 does not rely on the total number of genetic variants, while Majority rule and Plurality rule assumptions do.We now show that under Assumption 1, one can identify the class of invalid IV sets.Proposition 2 (Identifying Invalid IV Sets). Let G = {G 1 , . . . , G p }, p ≥ 2 be a subset of candidate genetic variants G. Suppose that Assumption 1 and the size of participants n → ∞ hold. If there exists a G j ∈ G such that,EQUATION)then G is an invalid IV set for any one of the causal relationships.Proposition 2 is valuable when given a set of genetic variants but lack prior knowledge of the validity of IV sets. By testing specific correlations in Eq.( 5), we can conclude that the current IV set is invalid.Example 1. Consider the example in Figure 2 .Suppose G = {G 1 , G 2 }, according to Proposition 2 we have corr(PR (X,Y | {G2}) , G 1 ) ̸ = 0. Thus G is an invalid IV set. U G2 G1 X Y (a) U G2 G1 X Y (b)Figure 3 : An illustrative example that valid and invalid IV sets may induce the same constraints in Eq.( 5) of Proposition 2, where {G 1 , G 2 } is the set of valid IVs in (a), while invalid in (b). It implies Assumption 1 is not sufficient to find valid IV sets.One may naturally raise the following question: Can all the invalid IV sets be identified by Proposition 2? In other words, can we identify all the valid IV sets? Unfortunately, Assumption 1 is an insufficient condition for identifying the sets of valid IVs. See an illustrative example below.Example 2 (Counterexample). Consider the example in Figure 3 . Suppose G = {G 1 , G 2 }. In subgraph (a), for G 1 ∈ G, according to Proposition 2 we obtaincorr(PR (X,Y | {G2}) , G 1 ) = 0 (6)Similarly, for G 2 ∈ G, we havecorr(PR (X,Y | {G1}) , G 2 ) = 0 (7)In subgraph (b), for G 1 ∈ G, we check the condition of Proposition 2 and obtainEQUATION)Similarly, for G 2 ∈ G, we haveEQUATION)We know that {G 1 , G 2 } is an invalid IV set in this case. However, if the proportion of effects of G 1 on X and Y is equal to the proportion of the effects of G 2 on X andY , i.e., γ X,1 γ Y,2 = γ Y,1 γ X,2 , we have corr(PR (X,Y | {G2}) , G 1 ) = 0, corr(PR (X,Y | {G1}) , G 2 ) = 0,which is consistent with the conclusion of subgraph (a). These observations imply that one can not identify valid IV sets under only Assumption 1.To give sufficient conditions for identifying valid IV sets, we first introduce the following assumption.Assumption 2 (Generic Identifiability). For a given MR causal model, parameters γ and β live in a set of Lebesgue measure non-zero.In other words, Assumption 2 specifically implies that given a general subset of candidate genetic variants G, there exists aG j ∈ G that satisfies | Gi∈G/Gj βi cov(G i , G i )[(γ X,i + γ X,U γ U,i )(γ Y,j + γ Y,U γ U,j ) -(γ Y,i + γ Y,U γ U,i )(γ X,j + γ X,U γ U,j )]| ̸ = 0, where βi ̸ = 0 is the estimated coefficient between G i and X. For instance, in Example 2, such an assumption simplifies to γ X,1 γ Y,2 ̸ = γ Y,1 γ X,2 . Whenever γ X,1 γ Y,2 = γ Y,1 γ X,EQUATION)then G is a proper valid IV set for one of the causal relationships, i.e.,G ⊆ G X→Y V or G ⊆ G Y →X V .Proposition 3 indicates that, under Assumptions 1 and 2, we are able to identify a subset that consists of valid IVs for one of the relationships, from the candidate genetic variants G. Here, it should be noted that in a bi-directional MR model (with both β X→Y ̸ = 0 and β Y →X ̸ = 0), one cannot determine whether the identified IV set G is related to the causal relationship X → Y or Y → X.To address this, we here introduce the following conditions that allow the estimated valid IV sets to determine the direction.Assumption 3. In a bi-directional MR model, for a given causal relationship X → Y , β 2 X→Y • Var(X) < Var(Y ). Similarly, for a given causal relationship Y → X, β 2 Y→X • Var(Y ) ≤ Var(X).Assumption 3 is a very natural condition that one expects to hold for the unique identifiability of valid IVs for the bi-directional model. This assumption is often reasonable in MR models as discussed in Xue and Pan [2020] . One may refer to Pages 4∼5 in Xue and Pan [2020] for more details. Proposition 4 (Identifying Direction of Causal Influences). Suppose Assumption 3 holds. Given a valid IV set denoted asG, if each G j ∈ G satisfies the condition |corr(G j , Y )/corr(G j , X)| < 1, then G is a valid IV set for the causal relationship X → Y , i.e., G ⊆ G X→Y V . Conversely, if this condition is not satisfied, G is a valid IV set for the causal relationship Y → X, i.e., G ⊆ G Y →X V .Proposition 4 indicates that, under Assumption 3, it is possible to ascertain the directional orientation of the identified IV sets. Thus, with Propositions 2∼4, we have sufficient and necessary conditions under which the bi-directional MR model is fully identifiable, formalized as follows: Theorem 1 (Identifiability of Bi-directional MR model). Suppose that Assumptions 1 ∼ 3, and the size of participants n → ∞ hold, the Bi-directional MR model is fully identifiable, which includes the identification of valid IV sets, the causal relationships the identified sets relate to, the causal directions between X and Y , and the causal effects of X on Y as well as of Y on X. 
prebim algorithm: 
 In this section, we leverage the above theoretical results and propose a data-driven algorithm, PRe-BiM (Pseudo-Residual-based Bi-directional MR), for estimating the causal effects of interest. It contains two key steps, finding valid IV sets (Step I), and inferring causal direction and estimating the causal effects given the identified valid IV sets (Step II). The algorithm is designed with the following rules, which are proved correct: (1) valid IV sets can be correctly discovered given the measured genetic variants (Section 4.1), and (2) the causal direction and causal effects can be uniquely identified (Section 4.2). The entire process is summarized in Algorithm 1.Input: A dataset of measured genetic variants G = (G 1 , ..., G g ) ⊺ , two phenotypes X and Y , significance level α, and parameter W , maximum number of IVs to consider.1: Valid IV sets V ← FindValidIVSets(G, X, Y , α, W ); ▷ Step I 2: ( βX→Y , βY →X ) ← InferCausalDirectionEffects(X, Y , V); ▷Step II Output: βX→Y and βY →X , the causal effects of X on Y and Y on X, respectively.Propositions 2 and 3 have paved the way to distinguish the invalid and valid IV sets, respectively. This enables us to use a cluster fusion-like method to discover the valid IV sets4 . To improve the algorithm's effectiveness with a finite sample size, especially when the pool of candidate genetic variants G is substantial, we introduce a parameter W , representing the maximum number of instrumental variables (IVs) to consider. This approach aims to circumvent the necessity of identifying the largest IV set within the collection, thereby minimizing the need for IV validity testing.The basic process is provided in Algorithm 2. Due to limited space, the specific details of algorithm execution are provided in Appendix E. In practice, we check for the valid IV test, i.e., Eq.( 10), with the Pearson Correlation Test using the Fisher transformation [Anderson, 2003] , denoted by PCT. Input: A dataset of measured genetic variants G = (G 1 , ..., G g ) ⊺ ,V ′ new . If V is empty, add V ′ new to V.If not, and the combined set {V ∪ V ′ new } passes the PCT, merge them; otherwise, discard V ′ new . Stop the iteration when V contains two sets of IVs or when there is one or fewer genetic variants left in G.Output: V, a set that collects the valid IV sets.In this section, we infer the causal direction and estimate the causal effects given the identified valid IV sets. The basic process is provided in Algorithm 3. Due to limited space, the specific details of algorithm execution are provided in Appendix F. In this section, we show that, in the large sample limit, for the causal relationships of interest X ↔ Y , the causal effects βX→Y and βY →X obtained through the proposed method are unbiased. Theorem 2 (Correctness). Assume that the data {X, Y } and G strictly follow Eq.( 1) and Assumptions 1 ∼ 3 hold. Given infinite samples, the PReBiM algorithm outputs the true causal effects β X→Y and β Y →X correctly.We give the computational complexity analysis as follows. First, the complexity of finding a valid IV set V new of length 2 is O 2 g 2 , where g is the number of IVs. Then, the complexity of expandingV new to the length W is O (2g-W -2)(W -1), where W is the maximum length of the IV set. For finding additional valid IV sets, the number of loops, t, must satisfy t ≤ g-4 W . The complexity of finding the second valid IV set of length W is given by Ot k=1 2 g-kW 2 + (2g-(2k+1)W -2)(W -1) 2 + 2W (t -1) + 1 . Therefore, the total computational complexity of the PReBiM algorithm is O t k=0 2 g-kW 2 + (2g-(2k+1)W -2)(W -1) 2 + 2W (t -1) + 1 . 5 Discussion U G2 G1 X Y (a) U G2 G1 X Y (b)Figure 4 : Simple examples, where{G 1 , G 2 } is a valid IV set in (a), whereas invalid in (b).The preceding sections demonstrate that a valid IV set can be identified when any two genetic variants in G are independent. In this section, we explore whether the main results of Propositions 2 and 3 remain applicable when some genetic variants in G are dependent. As shown in Figure 4 (b), we evaluate the condition of Proposition 2 and find that corr(PR(X,Y,|G2) , G 1 ) ̸ = 0 and corr(PR (X,Y,|G1) , G 2 ) ̸ = 0. This indicates that {G 1 , G 2 } is an invalid IV set. Conversely, in Fig- ure 4 (a), when assessing the condition of Proposition 3, we observe that corr(PR (X,Y,|G2) , G 1 ) = 0 and corr(PR (X,Y,|G1) , G 2 ) = 0, suggesting that {G 1 , G 2 } is a valid IV set.Consequently, it is interesting to note that even with such dependence between genetic variants, our main results may still be effective in identifying valid IV sets. In Appendix H.1 ∼ H.2, we further demonstrate the practicality of our method for this example. As for the complex scenarios where there is confounding among genetic variants or between the genetic variants and some phenotype, or where there exist genetic variants that are associated with another phenotype while this phenotype causes X and Y , we leave them into our future work. 
experimental results: 
 In this section, we conduct various simulation studies to evaluate the performance of the proposed PReBiM method. The methods we compare against are: 1) NAIVE, the least-squares regression method; 2) MR-Egger [Bowden et al., 2015] ; 3) sisVIVE algorithm [Kang et al., 2016] 5 ; 4) IV-TETRAD algorithm [Silva and Shimizu, 2017] 6 ; 5) TSHT algorithm [Guo et al., 2018] 7 ; and 6) PReBiM, our method. The second method, MR-Egger, is specifically for MR studies with invalid IVs. The last four methods are capable of handling the one-sample MR model in the presence of unmeasured confounders. Our source code can be found in the Supplementary Materials.Experimental setup. To compare the performance of these methods in a realistic setting, analogous to Slob and Burgess [2020] , the genetic variants are modeled as Single Nucleotide Polymorphisms (SNPs), with a varying minor allele frequency maf j , and take values 0, 1, or 2. The minor allele frequencies are drawn from a uniform distribution. Specifically, the data generation process for the bi-directional model is as follows:EQUATIONwhere the error terms ε 1 , ε 2 , ε 3 each follow an independent normal distribution with mean 0 and unit variance. The causal effects β Y →X and β X→Y are generated from a uniform distribution between [-1, -0.5] ∪ [0.5, 1]. Note that by setting β Y →X = 0 in Eq.( 12), the data generation process for the one directional model can be obtained.We here conducted experiments with the following tasks.• T1: Sensitivity to Sample Size. We evaluated the impact of different sample sizes: n = 2k, 5k, and 10k, where k equals 1, 000.• T2: Sensitivity to the Number of Valid IVs. We explored three distinct IV combination scenarios: S(2, 2, 6), S(3, 3, 8), and S(4, 4, 10). The first two elements in the parentheses represent the number of valid instrumental variables for causal relationships X → Y and Y → X, respectively. The last element in the parentheses represents the total number of genetic variants.For a comprehensive analysis, we examined the data from the Bi-directional MR model and the One-directional MR model, as detailed in Section 6.1 and 6.2, respectively. For each experimental setting, we repeated the simulation 500 times and averaged the results.Metrics: We adopted the following evaluation metrics to evaluate the results of all methods.• (i) Correct-Selecting Rate (CSR): the number of correctly identified valid IVs divided by the total number of genetic variants, for the direction X → Y or Y → X. • (ii) Mean Squared Error (MSE): The squared difference between the estimated causal effect and the true causal effect.To further validate the effectiveness and practicality of our algorithm, we have expanded our experimental settings to include additional scenarios such as small sample sizes, different numbers of valid IVs, large-scale datasets, and cases where all IVs are valid. Within these settings, we compare our algorithm with the five aforementioned methods to assess algorithm performance. Detailed experimental results are provided in Appendix H.3 ∼ H.6.In this section, we evaluate the performance of six methods on bi-directional MR data, i.e., β X→Y ̸ = 0 and β Y →X ̸ = 0. Table 1 summarizes the results of metrics CSR and MSE in a bidirectional MR model with various sample sizes. It is noted that sisVIVE, IV-TETRAD, and TSHT are not specifically designed for the bi-directional model, so we run them in both forward and reverse directions to obtain the corresponding results. As expected, our proposed PReBiM method gives the best results in all three scenarios, with all sample sizes, indicating that it can identify valid IVs and derive consistent causal effects from purely observational data, even without knowledge of the causal direction between X and Y . Furthermore, the performance of the other three comparison methods of selecting IVs was poor, indicating that they are not suitable for bi-directional MR models.The MSE of the NAIVE method does not tend toward 0, indicating that if unobserved confounding is not considered, the estimated causal effect is certainly biased. Similarly, the MSE of MR-Egger is poor if estimated with all instrumental variables.Table 1 : Performance comparison of NAIVE, MR-Egger, sisVIVE, IV-TETRAD, TSHT, and PRe-BiM in estimating bi-directional MR models across various sample sizes and three scenarios. Note: The symbol "-" means methods have no output. "↑" means a higher value is better, and vice versa.S(2, 2, 6) S(3, 3, 8) S(4, 4, 10) X → Y Y → X X → Y Y → X X → Y Y → X Size AlgorithmIn this section, we evaluate the performance of six methods on one-directional MR data, e.g., β X→Y ̸ = 0 and β Y →X = 0. For fairness, we assume that the causal direction X → Y is known for all methods. Figure 5 summarizes the results of metrics CSR (solid lines) and MSE (dashed lines) for all methods in a one-directional MR model with different scenarios. As expected, the metrics CSR and MSE of our method approach the expected value as the sample size increases. The overall results are comparable to the IV-TETRAD algorithm and superior to the sisVIVE algorithm and TSHT algorithm. This is because our method, similar to the IV-TETRAD method, requires fewer valid IVs in the system (at least 2), whereas the other two methods require the number of valid IVs that are related to the overall genetic variants, namely the majority rule assumption (over 50% of the candidate IVs are valid) and plurality rule assumption (the number of valid IVs is larger than any number of invalid IVs). It has been noted that under the IV-TETRAD method for scenario S(2, 0, 6), the metrics CSR and MSE outperform our method. This is because the IV-TETRAD method requires a predetermined number of IV sets to be specified in advance, and in this case, we have predetermined the number of valid IVs for it. However, our method does not have this prior information. We first apply our method to analyze the bi-directional causal relationships between obesity and vitamin D status using the GWAS data from Vimaleswaran et al. [2013] . Note that the original data only consists of summary data, and thus we generate individual-level data based on the data generation mechanism described in Eq.( 1). The data we used comprises 16 related SNPs as candidate IVs, including seven SNPs associated with obesity (BMI): rs9939609 (FTO), rs2867125 (TMEM18), rs4074134 (BDNF), rs7647305 (ETV5), rs7138803 (FAIM2), rs3101336 (NEGR1), and rs10938397 (GNPDA2); four SNPs related to 25-hydroxyvitamin D (25(OH)D): rs12785878 (DHCR7), rs10741657 (CYP2R1), rs2282679 (GC), and rs6013897 (CYP24A1); and five SNPs associated with BMI and 25(OH)D: rs4498364, rs1381660, rs4549702, rs12760109, and rs77913838. We observed that FTO and FAIM2 serve as valid IVs for the X → Y , while DHCR7 and CYP24A1 are valid IVs for the Y → X, with estimated causal effects of -1.15 and -0.05, respectively. These findings are consistent with those findings in [Vimaleswaran et al., 2013] .Next, we apply our method to analyze the causal effect of institutions on economic development using the Colonial Origins dataset from [Acemoglu et al., 2001] . The data includes eight candidate IVs: absolute latitude (lat_abst), the density of European settlements in 1900 (euro1900), the level of democracy in 1900 (democ1), constitutional constraints in 1900 (cons1), the timing of industrialization (indtime), democracy in 2000 (democ00a), constitutional constraints in 2000 (cons00a), and log-transformed European settler mortality rates (logem4). Our method identifies euro1900 and logem4 as valid IVs, with an estimated causal effect of 0.861, which is consistent with the findings reported in Acemoglu et al. [2001] .Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: We placed the detailed experimental settings of the experiments in Section 6 of the paper and the complete details with the code in the supplementary materials. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusion: 
 We investigated the identifiability conditions for the bi-directional MR model using purely observational data. First, we proposed a testable constraint for identifying valid IV sets, and we demonstrated that, under appropriate conditions, these valid IV sets can be partially or even fully identified by checking the proposed constraints. Experimental results on simulation data further verified the effectiveness of our algorithm. Currently, we assume that genetic variants are randomized. Allowing for general dependent genetic variants is a direction for future work. Additionally, while we focus on a linear bi-directional MR model, in future research, we will generalize to address nonlinear or even more complex data.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: Our experiments included mean square error (MSE) for the evaluation metrics, which expresses the statistical significance of the experiment. We show it in all the tables in the paper. It is obtained by taking the mean by calculating the square of the difference between the estimated causal effect and the true value.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: The basic goal of epidemiologic studies is to assess the effect of changes in exposure on outcomes.The causal effect of exposure on outcome is often different from the observed correlation due to confounding factors. Correlations between exposure and outcome cannot be used as reliable evidence for inferring causality. In contrast, Mendelian Randomization (MR) studies use genetic variants to infer the causal effect of exposure on outcome. The idea behind MR studies is to find a genetic variant (or multiple genetic variants) that is associated with exposure but not with other confounders and that does not directly affect the outcome. Such a genetic variant in principle fulfills the basic assumptions of instrumental variables (IVs). It is then used to assess the causal effect of the exposure on the outcome [Thomas and Conti, 2004] . Since the assumptions of IVs cannot be fully tested and may be violated during IVs' selection, a number of methods have emerged to identify valid IVs.MR studies can be conducted using a single sample, i.e., one-sample Mendelian Randomization, in which individuals are tested for genetic variation, exposure, and outcome in the same population.On the contrary, in two-sample Mendelian Randomization, the association between genetic variation and exposure is estimated in one dataset while that between genetic variation and outcome is estimated in another dataset. Compared with the Two-sample MR, one-sample MR provides more direct evidence by eliminating potential biases that can arise from individual-level data, which turns out to be more meaningful yet challenging.Additionally, most existing MR methods assume a one-directional causal relationship between exposure and outcome, whereas bi-directional relationships are ubiquitous in real-life scenarios. For instance, there exist bi-directional relationships between obesity and vitamin D status [Vimaleswaran et al., 2013] , body mass index (BMI) and type 2 diabetes (T2D), diastolic blood pressure and stroke [Xue and Pan, 2022] , insomnia and five major psychiatric disorders [Gao et al., 2019] , smoking and BMI [Carreras-Torres et al., 2018] , etc. It is thus desirable to take further research on bi-directional MR. Bi-directional Mendelian Randomization assesses not only the causal effect of the exposure on the outcome but also the effect of the outcome on the exposure. To achieve this, valid IVs for both the exposure and the outcome are needed. However, this is tough because genetic variants associated with the exposure are often also associated with the outcome, and vice versa. Our goal here is to identify and differentiate valid IVs for each direction within a one-sample MR framework and infer bi-directional causal effects.Regarding our bi-directional MR causal model, we give further relative descriptions in the following.(i) Similar to Li and Ye [2022] , our model in Eq.( 1) describes bi-directional causality in equilibrium, and in econometrics, it is known as the simultaneous equation model. According to Hausman [1983] , to avoid the recursive formalization, Eq.( 1) equivalently is transformed into Eq.( 2) withβ X→Y β Y →X ̸ = 1. (ii)For such a bi-directional (cyclic) structure, Mooij et al. [2013] take an Ordinary Differential Equation (ODEs) perspective and show under which conditions the equilibrium state of a system of the first-order ODEs can be described by a deterministic structural causal model (SCM). This is a very interesting work. Please note that our setting is different from those in Mooij et al. [2013] : our model is motivated by bi-directional MR studies, other than a dynamical system from the control field. Further, due to the specificity of MR studies, we can only obtain the snapshot data or equilibrium data, other than dynamical ones for the treatment, the outcome as well as genetic variants. G X→Y V , G Y →X V A set consisting of all valid IVs in G relative to X → Y , and Y → X, respectively G X→Y I , G Y →X I A set consisting at least one invalid IV in G relative to X → Y or Y → X, respectively TSLS(•)The procedure of estimating causal effects using two-stage least square methodω G Estimated causal effect derived from G PR (X,Y | G\G j )Pseudo-residual of {X, Y } relative to G\Gj P CT Pearson Correlation Test C More Details on Observations 1 and 2C.1 Observation 1 Since G X→Y V = (G 1 , G 3) ⊺ is a valid IV set, employing G 1 or G 3 in TSLS(•) entails unbiased causal effects, i.e., ω {Gi} = TSLS(X, Y, {G i }) = β X→Y with i ∈ {1, 3}. Thus, we have,Y -Xω {G1} = Y -Xβ X→Y = G 2 γ Y,2 + G 4 γ Y,4 + G 5 γ Y,5 + ε Y , Y -Xω {G3} = Y -Xβ X→Y = G 2 γ Y,2 + G 4 γ Y,4 + G 5 γ Y,5 + ε Y .Because both G 1 and G 3 are valid IVs, they are independent of ε Y ; and any two genetic variants are assumed to be independent of each other, we get,corr(Y -Xω {G3} , G 1 ) = 0, corr(Y -Xω {G1} , G 3 ) = 0.Without loss of generality, we here give details about how corr(Y -Xω {G4,G5} , G 2 ) ̸ = 0 is derived. The other two inequalities could be obtained in the same way. Following the procedure of TSLS(•), we first get the estimate of X with {G 4 , G 5 }, denoted by X, and obtain ω {G4,G5} using X and Y .X = G 4 β4 + G 5 β5 = G 4 (γ X,4 + γ Y,4 β Y →X )∆ + G 5 (γ X,5 + γ Y,5 β Y →X )∆ = (G 4 γ X,4 + G 4 γ Y,4 β Y →X + G 5 γ X,5 + G 5 γ Y,5 β Y →X )∆,where since G 4 and G 5 are neither correlated with the latent confounder U nor the other IVs as shown in Figure 2 , the second equality holds. Then, we could get the estimate of ω {G4,G5} , i.e.,ω {G4,G5} = cov(G 4 , G 4 )(γ X,4 + γ Y,4 β Y →X )(γ Y,4 + γ X,4 β X→Y ) cov(G 4 , G 4 )(γ X,4 + γ Y,4 β Y →X ) 2 + cov(G 5 , G 5 )(γ X,5 + γ Y,5 β Y →X ) 2 + cov(G 5 , G 5 )(γ X,5 + γ Y,5 β Y →X )(γ Y,5 + γ X,5 β X→Y ) cov(G 4 , G 4 )(γ X,4 + γ Y,4 β Y →X ) 2 + cov(G 5 , G 5 )(γ X,5 + γ Y,5 β Y →X ) 2 . Then we get, corr(Y -Xω {G4,G5} , G 2 ) = (γ Y,2 + γ X,2 β X→Y )∆ -(γ X,2 + γ Y,2 β Y →X )∆ω {G4,G5} = [γ X,2 (β X→Y -ω {G4,G5} ) + γ Y,2 (1 -β Y →X ω {G4,G5} )]∆. Since, β X→Y -ω {G4,G5} = β X→Y (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + β X→Y (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) - (γ X,4 + γ Y,4 β Y →X )(γ Y,4 + γ X,4 β X→Y )cov(G 4 , G 4 ) (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) - (γ X,5 + γ Y,5 β Y →X )(γ Y,5 + γ X,5 β X→Y )cov(G 5 , G 5 ) (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) = (β X→Y β Y →X -1)[γ Y,4 (γ X,4 + γ Y,4 β Y →X ) + γ Y,5 (γ X,5 + γ Y,5 β Y →X )] (γ X,4 + γ Y,4 β Y →X ) 2 + (γ X,5 + γ Y,5 β Y →X ) 2 = (β X→Y β Y →X -1)(C 44 Γ X,4 γ Y,4 + C 55 Γ X,5 γ Y,5 ) C 44 Γ 2 X,4 + C 55 Γ 2 X,5, and1 -β Y →X ω {G4,G5} = (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) - β Y →X (γ X,4 + γ Y,4 β Y →X )(γ Y,4 + γ X,4 β X→Y )cov(G 4 , G 4 ) (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) - β Y →X (γ X,5 + γ Y,5 β Y →X )(γ Y,5 + γ X,5 β X→Y )cov(G 5 , G 5 ) (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) = [γ X,4 (γ X,4 + γ Y,4 β Y →X )cov(G 4 , G 4 ) + γ X,5 (γ X,5 + γ Y,5 β Y →X )cov(G 5 , G 5 )]∆ (γ X,4 + γ Y,4 β Y →X ) 2 cov(G 4 , G 4 ) + (γ X,5 + γ Y,5 β Y →X ) 2 cov(G 5 , G 5 ) = ∆(C 44 Γ X,4 γ X,4 + C 55 Γ X,5 γ X,5 ) C 44 Γ 2 X,4 + C 44 Γ 2 X,5 , Where we define Γ X,i = γ X,i + γ Y,i β Y →X , Γ Y,j = γ Y,j + γ X,j β X→Y , C ij = cov(G i , G j ) and ∆ ≜ 1 -β Y →X β X→Y . we can finally obtain corr(Y -Xω {G4,G5} , G 2 ), that is, corr(Y -Xω {G4,G5} , G 2 ) = γ X,2 (β X→Y β Y →X -1)(C 44 Γ X,4 γ Y,4 + C 55 Γ X,5 γ Y,5 )∆ C 44 Γ 2 X,4 + C 55 Γ 2 X,5 + γ Y,2 (1 -β Y →X β X→Y )(C 44 Γ X,4 γ X,4 + C 55 Γ X,5 γ X,5 )∆ C 44 Γ 2 X,4 + C 55 Γ 2 X,5 = C 44 Γ X,4 (γ X,2 γ Y,4 -γ X,4 γ Y,2 ) + C 55 Γ X,5 (γ X,2 γ Y,5 -γ X,5 γ Y,2 ) C 44 Γ 2 X,4 + C 55 Γ 2 X,5 ̸ = 0.Figure 6 gives an illustration of our idea, demonstrating that even when IVs are correlated, we can still identify them as valid IVs using our proposed method. In subgraph (a), there are two valid IVs: G 1 and G 2 . However, in subgraph (b), these two IVs become invalid as G 1 is correlated with Y , rendering them invalid.U G2 G1 X Y (a) U G2 G1 X Y (b)Figure 6 : A simple example, where the set {G 1 , G 2 } is a valid IV set in subgraph (a), whereas it is not valid in subgraph (b).Assume the generating process of X and Y in subgraph (a) follows:EQUATION)where ε X and ε Y are correlated, and ε 2 represents the independent error term.For G 1 ∈ G, we check the condition of Proposition 3:ω {G1} = TSLS(X, Y, {G 1 }) = cov(Y, G 1 ) cov(X, G 1 ) = β X→Y + cov(ε Y , G 1 ) cov(X, G 1 ) = β X→Y PR (X,Y | {G1}) = Y -ω {G1} X = ε Y corr(PR (X,Y | {G1}) , G 2 ) = 0Similarly, for G 2 ∈ G, we check the condition of Proposition 3:ω {G2} = TSLS(X, Y, {G 2 }) = cov(Y, G 2 ) cov(X, G 2 ) = β X→Y + cov(ε Y , G 2 ) cov(X, G 2 ) = β X→Y PR (X,Y | {G2}) = Y -ω {G2} X = ε Y corr(PR (X,Y | {G2}) , G 1 ) = 0Assume the generating process of X and Y in subgraph (b) follows Equation 16:EQUATION)where ε X and ε Y are correlated, and ε 2 represents the independent error term. Then, Equation 16 can be reorganized to avoid recursive formulations as follows:EQUATIONwhereΓ X,1 = γ X,1 + γ 1,2 γ X,2 + γ Y,1 β Y →X and Γ Y,1 = γ Y,1 + β X→Y (γ X,1 + γ X,2 γ 1,2 ). It is noted that Γ X,1 ̸ = 0 and Γ Y,1 ̸ = 0.In subgraph (b), for G 2 ∈ G, we check the condition of Proposition 3: ε2,ε2) cov(G1,G1) ̸ = 0. Based on these two observations, we can conclude that the two IVs in subgraph (a) can be deemed valid. However, in subgraph (b), they can be deemed invalid.ω {G2} = TSLS(X, Y, {G 2 }) = cov(Y, G 2 ) cov(X, G 2 ) = β X→Y + γ Y,1 cov(G 1 , G 2 ) cov(X, G 2 ) = β X→Y + γ Y,1 γ 1,2 ∆ cov(G 1 , G 1 ) cov(G 1 , G 1 )(γ 1,2 Γ X,1 ) + cov(ε 2 , ε 2 )γ X,2 = β X→Y + γ Y,1 γ 1,2 ∆ 1 γ 1,2 Γ X,1 + γ X,2 M , PR (X,Y | {G2}) = Y -ω {G2} X = G 1 γ Y,1 + ε Y - γ Y,1 γ 1,2 ∆ 1 γ 1,2 Γ X,1 + γ X,2 M X = G 1 γ Y,1 + ε Y - γ Y,1 γ 1,2 (G 1 Γ X,1 + ε 2 γ X,2 + ε X + ε Y β Y →X ) γ 1,2 Γ X,1 + γ X,2 M , corr(PR (X,Y | {G2}) , G 1 ) = γ Y,1 - γ Y,1 γ 1,2 Γ X,1 γ 1,2 Γ X,1 + γ X,2 M = γ Y,1 γ 1,2 Γ X,1 + γ Y,1 γ X,2 M -γ Y,1 γ 1,2 Γ X,1 γ 1,2 Γ X,1 + γ X,2 M = γ Y,1 γ X,2 M γ 1,2 Γ X,1 + γ X,2 M (18) where M = cov(It is noteworthy that, to mitigate the unreliability of statistical tests in finite samples, in practice, we choose to take the smallest correlation coefficient to select a valid IV set, as in Line 4 or 10 of Algorithm 4. V new ← arg min G∈G 2 i=1 corr(PR (X,Y | G\Gi) , G i ). 5: repeat 6: G′ ← G 7:For each variable G k ∈ G′ , remove all variables from G′ corresponding to those that set {V, G k } fail to the PCT.if G′ ̸ = ∅ then 9:Compute the correlation between pseudo-residual PR (X,Y | V) and G k .10:G ′ ← arg min G k ∈ G′ corr(PR (X,Y | V) , G k ).11:V new ← V new ∪ {G k }.12:else 13:Break the loop of line 5;14:end if 15:until the length of the set V new has reached our predefined threshold W .16:if Set V = ∅ then 17:Add set V new into V. Delete set V new from G.end ifend if 25: until V contains two sets or | G| ≤ 1. Output: V, a set that collects the valid IV sets. if |corr(G j , Y )/corr(G j , X)| < 1 then 4: Add G j into G X→Y V ; 5: else 6: Add G j into G Y →X V ; 7: end if 8: end for 9: βX→Y ← TSLS(X, Y, G X→Y V ). 10: βY →X ← TSLS(Y, X, G Y →X V ). Output: βX→Y and βY →X G Proofs G.1 Proof of Proposition 1Proof. From procedure "two-stage least squares", we obtained the fitted values X from the regression in the first stage:X = (G X→Y V ) ⊺ α. According to OLS regression, X = (G X→Y V ) ⊺ (G X→Y V (G X→Y V ) ⊺ ) -1 G X→Y V X = (G X→Y V ) ⊺ (G X→Y V (G X→Y V ) ⊺ ) -1 G X→Y V X = PX.Then run the OLS regression for Y on X in second stage: Y = β X. where:β = ( X⊺ X) -1 ( X⊺ Y ) = ((PX) ⊺ PX) -1 ((PX) ⊺ Y ) = (X ⊺ PX) -1 (X ⊺ PY ) = (X ⊺ PX) -1 (X ⊺ P) Xβ X→Y + (G X→Y V ) ⊺ γ Y + ε Y = β X→Y + (X ⊺ PX) -1 (X ⊺ P(G X→Y V ) ⊺ γ Y ) + (X ⊺ PX) -1 (X ⊺ Pε Y ) Since G X→Y V is a valid IV set, for G j ∈ G X→Y V , both γ Y,j = 0 (A2) and cov(G j , ε Y ) = 0 (A3).Hence, we obtained(X ⊺ PX) -1 (X ⊺ P(G X→Y V ) ⊺ γ Y ) = 0 (X ⊺ PX) -1 (X ⊺ Pε Y ) = (X ⊺ PX) -1 (X ⊺ (G X→Y V ) ⊺ (G X→Y V (G X→Y V ) ⊺ ) -1 G X→Y V ε Y = 0, which implies that β = β X→Y .Proof. Similar to G.1, from procedure "two-stage least squares", we obtained the fitted values X from the regression in the first stage:X = (G X→Y I ) ⊺ α. According to OLS regression, X = (G X→Y I ) ⊺ (G X→Y I (G X→Y I ) ⊺ ) -1 G X→Y I X = (G X→Y I ) ⊺ (G X→Y I (G X→Y I ) ⊺ ) -1 G X→Y I X = PX.Then run the OLS regression for Y on X in second stage: Y = β X. where:  and A3 ). Hence, we obtainedβ = ( X⊺ X) -1 ( X⊺ Y ) = (( PX) ⊺ PX) -1 (( PX) ⊺ Y ) = (X ⊺ PX) -1 (X ⊺ PY ) = (X ⊺ PX) -1 (X ⊺ P) Xβ X→Y + (G X→Y I ) ⊺ γ Y + ε Y = β X→Y + X ⊺ PX -1 X ⊺ P(G X→Y I ) ⊺ γ Y + ε Y ) β bias = β X→Y + (X ⊺ PX) -1 (X ⊺ P(G X→Y I ) ⊺ γ Y ) + (X ⊺ PX) -1 (X ⊺ Pε Y ) Since G X→Y I is a valid IV set, for each G j ∈ G X→Y I , either γ Y,j ̸ = 0 (A2), or cov(G j , ε Y ) ̸ = 0 (A3), or both γ Y,j ̸ = 0 and cov(G j , ε Y ) ̸ = 0 (A2(X ⊺ PX) -1 (X ⊺ P(G X→Y I ) ⊺ γ Y ) ̸ = 0 (X ⊺ PX) -1 (X ⊺ Pε Y ) = (X ⊺ PX) -1 (X ⊺ (G X→Y I ) ⊺ (G X→Y I (G X→Y I ) ⊺ ) -1 G X→Y I ε Y ̸ = 0,which implies that β bias ̸ = 0. That is to say the estimated causal effect β = β X→Y + β bias is biased.Proof. We prove it by contradiction, i.e., if G is a valid IV set for one of the causal relationships, then for each G j ∈ G, corr(PR (X,Y | G\Gj ) , G j ) = 0. Without loss of generality, we consider the case for X → Y . The proof with the inverse causal relationship can be derived accordingly. i) assume that G is the valid IV set relative to the causal relationship X → Y , i.e., G ⊆ G X→Y V . Let G j be a genetic variant in G. We can see that G\G j is also a valid IV set. Because the size of participants n → ∞, and according to the consistent estimator of TSLS, we haveω {G\Gj } = TSLS(X, Y, G\G j ) = β X→Y .Based on Eq.( 1), we havePR (X,Y | G\Gj ) = Y -Xω {G\Gj } = G ⊺ γ Y + ε Y Because G is a valid IV set relative to X → Y , for G k ∈ G, k ̸ = j, either γ Y,k = 0 (if G k ∈ G) or corr(G k , G j ) = 0 (if G k / ∈ G).It is noted that G k can not be correlated with G j , otherwise, the genetic variant G j will dissatisfy A2 or A3. Hence, we obtain corr(G ⊺ γ Y + ε Y , G j ) = 0, which implies that corr(PR (X,Y | G\Gj ) , G j ) = 0. ii) assume that G is the valid IV set relative to Y → X, i.e., G ⊆ G Y →X V. Let G j be a genetic variant in G. We can see that G\G j is also a valid IV set relative to Y → X. It is noted that we here still need to verify Eq.( 5). According to the TSLS estimator and Eq.( 1), we haveω {G\Gj } = TSLS(X, Y, G\G j ) = X ⊺ PX -1 X ⊺ PY = X ⊺ PX -1 X ⊺ P(1/β Y →X )(X -G ⊺ γ X -ε X ) = (1/β Y →X ) 1 -X ⊺ PX -1 X ⊺ P(G ⊺ γ X + ε X ) , = 1 -X ⊺ PX -1 X ⊺ P(G ⊺ γ X + ε X ) β -1 Y →X , where P = {G\G j } [{G\G j } ⊺ {G\G j }] -1 {G\G j } ⊺ , X, Y and ε X are n × 1 vectors, G is n × p matrix and γ X is a g × 1 vector.Hence, we havePR (X,Y | G\Gj ) = Y -Xω {G\Gj } = Y -X 1 -X ⊺ PX -1 X ⊺ P(G ⊺ γ X + ε X ) β -1 Y →X = Y -Xβ -1 Y →X + X X ⊺ PX -1 X ⊺ P(G ⊺ γ X + ε X )β -1 Y →X = -(G ⊺ γ X + ε X )β -1 Y →X + X X ⊺ PX -1 X ⊺ P(G ⊺ γ X + ε X )β -1 Y →X = -I + X X ⊺ PX -1 X ⊺ P (G ⊺ γ X + ε X ) β -1 Y →X ,whereI is a n × n identity matrix. Because G is a valid IV set relative to Y → X, for G k ∈ G, k ̸ = j, either γ X,k = 0 (if G k ∈ G) or corr(G k , G j ) = 0 (if G k / ∈ G). Hence, corr(G ⊺ γ X + ε X ), G j ) = 0Therefore, we obtaincorr((X/β Y →X ) X ⊺ PX -1 X ⊺ P(G ⊺ γ X + ε X ), G j ) = 0 which implies that corr(PR (X,Y | G\Gj ) , G j ) = 0.Proof. We prove it by contradiction, i.e., if G is an invalid IV set for one of the causal relationships, then there exists at least oneG j ∈ G such that corr(PR (X,Y | G\Gj ) , G j ) ̸ = 0.Without loss of generality, we consider the case for X → Y . The proof with the inverse causal relationship Y → X can be derived accordingly. There are three cases to be considered. i) Suppose G j is an invalid IV in G and G\G j is a valid IV set; ii) Suppose G j is an invalid IV in G and G\G j is an invalid IV set; iii) Suppose G j is a valid IV in G and G\G j is an invalid IV set.i) Suppose G j is invalid and all IVs in G\G j are valid. Due to the validity of all IVs in G\G j , according to the consistent estimator of TSLS, we have,ω {G\Gj } = TSLS(X, Y, G\G j ) = β X→Y .In order to consider the case of all violating assumptions for IVs, we extend Eq.(1) to Eq.( 6), we havePR (X,Y | G\Gj ) = Y -Xω {G\Gj } = G ⊺ γ Y + ε Y = G ⊺ γ Y + U γ Y,U + ε 3 ,where we make ε Y = U γ Y,U + ε 3 to tell the influential difference between the latent confounder U and the external noise ε 2 on Y . Here G j ⊥ ⊥ ε 3 but G j may be correlated with U by γ U,j . Denote by G i the other IV with i ̸ = j, and we haveG j ⊥ ⊥ G i . Hence, corr(PR (X,Y | G\Gj , G j ) = γ Y,j corr(G j , G j ) + γ Y,U corr(U, G j ) = γ Y,j + γ Y,U γ U,j .Due to the faithfulness assumption, corr(PR (X,Y | G\Gj , G j ) ̸ = 0.ii) Suppose G j is invalid and G\G j is an invalid IV set, where there are m genetics in G\G j = {G 1 , ..., G m }. Following the procedure of TSLS, we first get the estimate of X with G\G j , denoted by X, and obtain ω {G\Gj } using X and Y .X = m k=1 G k βk = m k=1 G k ∆[(γ X,k + γ Y,k β Y →X ) + (γ X,U + γ Y,U β Y →X )γ U,k ].ω {G\Gj } = m k=1 [ βk (γ Y,k + γ X,k β X→Y ) + βk γ U,k (γ Y,U + γ X,U β X→Y )]C kk ∆ m k=1 ( βk ) 2 = m k=1 [ βk (γ Y,k + γ X,k β X→Y ) + βk γ U,k (γ Y,U + γ X,U β X→Y )]C kk m k=1 ( βm ) 2 , = m k=1 βk [(γ Y,k + γ X,k β X→Y ) + γ U,k (γ Y,U + γ X,U β X→Y )]C kk m k=1 ( βk ) 2 = m k=1 (A X,k + A X,U γ U,k )(A Y,k + A Y,U γ U,k )C kk m k=1 ( βk ) 2,where βk = βk /∆ = (γ X,k + γ Y,k β Y →X ) + γ U,k (γ X,U + γ Y,U β Y →X ) = A X,k + A X,U γ U,k , A X,k = γ X,k + γ Y,k β Y →X and A X,U = γ X,U + γ Y,U β Y →X . A Y,k = γ Y,k + γ X,k β X→Y and A Y,U = γ Y,U + γ X,U β X→Y and C kk = cov(G k , G k ).Note that such notations A • enjoy some great properties shown as below,A X,k -A Y,k β Y →X = γ X,k (1 -β X→Y β Y →X ), A X,k β X→Y -A Y,k = γ Y,k (β X→Y β Y →X -1), A X,1 A Y,2 -A X,2 A Y,1 = (1 -β X→Y β Y →X )(γ X,1 γ Y,2 -γ X,2 γ Y,1 ). (19) Hence, corr(PR (X,Y | G\Gj ), G j ) = (γ Y,j + γ X,j β X→Y )∆ -(γ X,j + γ Y,j β Y →X )∆ω {G\Gj } + γ U,j (γ Y,U + γ X,U β X→Y )∆ -γ U,j (γ X,U + γ Y,U β Y →X )∆ω {G\Gj } = ∆{[γ Y,j (1 -β Y →X ω {G\Gj } ) + γ X,j (β X→Y -ω {G\Gj } )] + γ U,j [γ Y,U (1 -β Y →X ω {G\Gj } ) + γ X,U (β X→Y -ω {G\Gj } )]} = ∆[(1 -β Y →X ω {G\Gj } )(γ Y,j + γ Y,U γ U,j ) + (β X→Y -ω {G\Gj } )(γ X,j + γ X,U γ U,j )].Employing the properties in Eqs.( 19), we get1 -β Y →X ω {G\Gj } = m k=1 (A X,k + A X,U γ U,k ) 2 C kk -β Y →X (A X,k + A X,U γ U,k )(A Y,k + A Y,U γ U,k )C kk m k=1 ( βk ) 2 C kk = (1 -β X→Y β Y →X ) m k=1 (A X,k + A X,U γ U,k )(γ X,k + γ X,U γ U,k )C kk m k=1 ( βk ) 2 C kk ,andβ X→Y -ω {G\Gj } = β X→Y m k=1 (A X,k + A X,U γ U,k ) 2 C kk - m k=1 (A X,k + A X,U γ U,k )(A Y,k + A Y,U γ U,k )C kk m k=1 ( βk ) 2 C kk = (β X→Y β Y →X -1) m k=1 (A X,k + A X,U γ U,k )(γ Y,k + γ Y,U γ U,k )C kk m k=1 ( βk ) 2 C kk . Thus, corr(PR (X,Y | G\Gj ), G j ) = m k=1 C kk (A X,k + A X,U γ U,k )(γ X,k + γ X,U γ U,k ) m k=1 ( βk ) 2 C kk (γ Y,j + γ Y,U γ U,j ) - m k=1 C kk (A X,k + A X,U γ U,k )(γ Y,k + γ Y,U γ U,k ) m k=1 ( βk ) 2 C kk (γ X,j + γ X,U γ U,j ) = m k=1 βk C kk [(γ X,k + γ X,U γ U,k )(γ Y,j + γ Y,U γ U,j ) -(γ Y,k + γ Y,U γ U,k )(γ X,j + γ X,U γ U,j )] m k=1 ( βk ) 2 C kk = ∑ m k=1 βk C kk [(γ X,k γ Y,j -γ Y,k γ X,j )+γ Y,U (γ X,k γ U,j -γ U,k γ X,j )+γ X,U (γ U,k γ Y,j -γ Y,k γ U,j )] m k=1 ( βk ) 2 C kk ̸ = 0.Because of the faithfulness assumption, βk ̸ = 0. Due to Assumption 2, for G k ∈ G\G j and the confounder U ,m k=1 βk C kk [(γ X,k γ Y,j -γ Y,k γ X,j )+γ Y,U (γ X,k γ U,j -γ U,k γ X,j )+γ X,U (γ U,k γ Y,j - γ Y,k γ U,j )] ̸ = 0.iii) Suppose G j is a valid IV in G and G\G j is an invalid IV set. Similar to ii), we get ω {G\Gj } :ω {G\Gj } = m k=1 (A X,k + A X,U γ U,k )(A Y,k + A Y,U γ U,k )C kk m k=1 ( βk ) 2 C kk .Since the valid IV G j is uncorrelated with neither other invalid IVs nor the confounder, we then derive,corr(PR (X,Y | G\Gj ), G j ) = (γ Y,j + γ X,j β X→Y )∆ -(γ X,j + γ Y,j β Y →X )∆ω {G\Gj } = ∆γ X,j (β X→Y -ω {G\Gj } ).where:β X→Y -ω {G\Gj } = β X→Y m k=1 C kk [(A X,k + A X,U γ U,k ) 2 -(A X,k + A X,U γ U,k )(A Y,k + A Y,U γ U,k )] m k=1 ( βk ) 2 C kk = (β X→Y β Y →X -1) m k=1 C kk (A X,k + A X,U γ U,k )(γ Y,k + γ Y,U γ U,k ) m k=1 ( βk ) 2 C kk .Thus, we getcorr(PR (X,Y | G\Gj ), G j ) = ∆γ X,j (β X→Y -ω {G\Gj } ) = - m k=1 C kk (A X,k + A X,U γ U,k )(γ Y,k + γ Y,U γ U,k )γ X,j m k=1 ( βk ) 2 C kk = - m k=1 βk C kk (γ Y,k + γ Y,U γ U,k )γ X,j m k=1 ( βk ) 2 C kk .Since G j is a valid IV (γ Y,j = γ U,j = 0), and from Assumption 2, we havem k=1 βk C kk (γ Y,k + γ Y,U γ U,k )γ X,j ̸ = 0. Proposition 3 is proved.Proof. We prove this proposition by contradiction, i.e., i) if G is the valid IV set relative to the causal relationship Y → X, then ∃ G j ∈ G satisfies | corr(Gj ,Y ) corr(Gj ,X) | ≥ 1; and ii) if G is the valid IV set relative to the causal relationship X → Y , then ∃ G j ∈ G satisfies | corr(Gj ,Y ) corr(Gj ,X) | < 1. i) Suppose G is the valid IV set relative to the causal relationship Y → X. According to Assumption 3, we have |β Y →X | ≤ Var(X)Var(Y ) . Hence,corr(G j , Y ) corr(G j , X) = Cov(G j , Y ) Var(G j )Var(Y ) Var(G j )Var(X) Cov(G j , X) = γ Y,j ∆ Var(Y ) Var(X) γ Y,j ∆β Y →X = 1 β Y →X Var(X) Var(Y ) ≥ 1, which proves case i).ii) Suppose G is the valid IV set relative to the causal relationship X → Y . According to Assumption 3, we have|β X→Y | < Var(Y ) Var(X) . Hence, corr(G j , Y ) corr(G j , X) = Cov(G j , Y ) Var(G j )Var(Y ) Var(G j )Var(X) Cov(G j , X) = γ X,j ∆β X→Y Var(Y ) Var(X) γ X,j ∆ = β X→Y Var(X) Var(Y ) < 1,which proves case ii) and concludes the proposition.Proof. Firstly, given Assumptions 1 and 2, in light of Propositions 2 and 3, we can identify all sets of valid instrumental variables (IVs). Next, given Assumption 3, by Proposition 4, we can further identify what the causal relationships corresponding to these IV sets. It is important to note that if all the IV sets correspond to a single causal relationship, such as X → Y , we can conclude that the model is a one-directional MR model. Finally, with the corresponding IV set identified for each causal influence, we can employ the classical IV estimator, namely the two-stage least squares (TSLS) estimator, to accurately estimate the causal effects.Based on the above analysis, we can conclude that the bi-directional MR model is fully identifiable.Proof. The correctness of PReBiM originates from the following observations:• Firstly, in Step 1 (Algorithm 4), by Propositions 2 and 3, valid IV sets in G have been exactly discovered in V (Lines 2 ∼ 25 of Algorithm 4).• Secondly, in Step 2 (Algorithm 5), by the condition of Proposition 4, the direction of causal influence can be determined (Line 2 ∼ 8 of Algorithm 5).• Lastly, by Eq.(1) of Proposition 1, unbiased causal effects can be estimated exactly given the identified IV sets (Lines 9 ∼ 10 of Algorithm 5).Based on the above analysis, we can identify all valid IV sets and obtain the causal effects correctly.H.1 Evaluating the Performance on Bi-directional MR Model with Dependent Valid IVs Table 3 summarizes the CSR and MSE metrics' results across six methods within a bi-directional MR model with varying sample sizes. As expected, our proposed method exhibits superior performance, even in scenarios where the included valid instrumental variables (IVs) are correlated. This highlights its capability to accurately identify effective IVs and provide consistent causal effect estimates solely from observational data, regardless of whether the relevant IVs are correlated. Importantly, this effectiveness persists without prior knowledge of the causal direction between X and Y . Conversely, the performance of the other three methods remains suboptimal in such contexts, indicating their limited applicability in cases with correlated valid IVs. Note: This scenario represents a bi-directional model with correlated valid IVs. The symbol "-" means that methods do not output information. "↑" means a higher value is better, and vice versa.In this section, we assess the efficacy of four methods within a one-directional Mendelian Randomization (MR) framework. For a fair comparison, we operate under the assumption that the causal pathway from X to Y is established for all methods. Figure 7 displays the CSR and MSE metrics for these methods, considering diverse sample sizes and numbers of valid instrumental variables (IVs).As expected, our approach maintains strong performance, even amidst correlated valid IVs. The outcomes are comparable with those of the TSHT algorithm, with our method notably surpassing both the sisVIVE and IV-TETRAD algorithms under various conditions. In this section, we will evaluate the performance of these six methods in a bi-directional Mendelian randomization framework with small samples. Note: The symbol "-" means that methods do not output this information. "↑" means a higher value is better, and vice versa.X → Y Y → X X → Y Y → X X → Y Y → X Size AlgorithmIn this section, we evaluate the performance of these methods in a bi-directional MR framework with large-scale datasets. Table 6 show results of the settings S(10, 10, 30) and S (15, 15, 40) . We observe that the overall performances of all methods decrease with large-scale datasets, whereas our method still performs more superior than other baselines.In this section, we evaluate the performance of these methods in a bi-directional MR framework with all IVs being valid. That is, only valid IVs are included in the model, either to estimate the direction X → Y or the direction Y → X. The results in Table 7 below show the settings S(2, 2, 4).To ensure the thoroughness of the experiment, we've added another setting, S (3, 3, 6) . As expected, our method gives the best results across all sample sizes and directions. As the sample size increases, the metric CSR approaches 1, indicating that our method can identify valid IVs and also validating our theoretical guarantees. d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 5. Open access to data and code Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: We anonymized and packaged our code and data, and put them in supplementary materials.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.Please note that these constraints can be tested easily by observed data.A cluster is considered a valid IV set, and the term "fusion-like" refers to the specific process of identifying and merging these clusters.For sisVIVE algorithm, we used the implementations in the R sisVIVE package, which can be downloaded at https://cran.r-project.org/web/packages/sisVIVE/.For IV-TETRAD method, we used the implementations in the R package, which can be downloaded at https://www.homepages.ucl.ac.uk/~ucgtrbd/code/iv_discovery/.For the TSHT algorithm, we used the implementations in the R RobustIV package, which can be downloaded at https://cran.r-project.org/web/packages/RobustIV/.Note: The symbol "-" means that methods do not output this information. "↑" means a higher value is better, and vice versa. 


Decision: The paper addresses the challenging problem of causal effect estimation from observational data in Mendelian randomization (MR). It studies the identification problem in bi-directional MR and presents complete theoretical results with both necessary and sufficient conditions. Additionally, the paper introduces an algorithm and reports experimental results on both synthetic and real-world datasets.

All reviewers agree that the work tackles a challenging and practical problem. The theoretical results are sound and complete, the algorithm is well-founded, and the evaluations are comprehensive. However, there was some concern about the initial lack of experiments on real-world datasets. In response, the authors added two new experiments involving real-world data in the discussion phase. Reviewers were also interested in understanding the impact of the strong assumptions on the feasibility of the algorithm. The authors provided explanations for these assumptions, which were satisfactory to the reviewers.

Reviews: review: paper_summary: paper_summary: This paper studies the identifiability problem of the bi-directional Mendelian randomization (MR) model, where $X$ and $Y$ are a pair of phenotypes of interest and causes of each other, and $\textbf{G}$ is the set of measured genetic variants, which may include invalid instrumental variables (IVs). Under some assumptions, the paper has identified and proved correct the sufficient and necessary conditions for identifying valid IV sets from $\textbf{G}$  based on observational data, without requiring prior knowledge about which candidate IVs in $\textbf{G}$ are valid or invalid. Supported by the theoretical result, an algorithm is proposed for finding the valid IV sets from the set of measured genetic variants  $\textbf{G}$ using observation data and estimating the bi-directional causal effects using the found valid IVs. Experiments are conducted with synthetic data to show the effectiveness of the proposed algorithm.
strength_weakness: strength_weakness: 1. The experimental evaluation is done with synthetic data only. Although the presented experiments with synthetic data are comprehensive and the identification conditions have been theoretically proved, as the theoretical result relies on several assumptions, it would be necessary to conduct some case studies with real world data to evaluate how the method works in practice (where domain knowledge or literature can be used to justify the correctness of the found IV sets)

2. It would be very helpful if the assumptions and their feasibility (and consequences/limitations) in practice can be illustrated and justified with real world examples.
questions: questions: 1. Line 106: Does the assumption regarding the multiplication of the two effects have any practical meaning/implication?
2. Could you explain what "cluster fusion" means exactly in the paper and why the proposed algorithm is said to be "cluster fusion-like"?
3. Section 6.2 - how the one-directional data used in this section generated?
4. The work is based on the assumed structure in Figure 2 (plus some invalid IVs as illustrated in the other figures) , but in practice there would be more complicated situations than those, e.g. the vertical pleiotropy effect in biology where the IVs (genetic variants) are associated with another phenotype (or biological pathway) and this in turn causes the two phenotypes of interest ($X$ and $Y$).
limitations: limitations: Some limitations of the paper have been discussed briefly, but as mentioned above, the consequence and limitations due to the assumptions should be discussed a bit more.
paper_summary: paper_summary: The paper addresses the problem of estimating causal effects in bi-directional Mendelian randomization (MR) models with some invalid instrumental variables (IVs) and unmeasured confounding. It proposes a framework for identifying valid IV sets under the assumption that the IV set consists of genetic variants that are independent of each other and that at least two of them are valid IVs.  The authors introduce a cluster fusion-like algorithm based on this framework and demonstrate its effectiveness through theoretical proofs and experimental results.
strength_weakness: strength_weakness: While the paper discusses various assumptions (such as the independence of genetic variants and existence of two valid IVs), it would benefit from a more in-depth exploration of the limitations and potential pitfalls of these assumptions in real-world data. Addressing how violations of these assumptions impact the results could strengthen the paper.

The experiments are performed on synthetic datasets. While this is a good starting point, additional validation on real-world datasets would provide more robust evidence of the method’s practical utility.
questions: questions: 1) How does your method perform when the assumptions (e.g., independence of genetic variants) are violated in practice? Are there any robust techniques or adjustments to handle such cases?

2) Regarding the construction of the IV set, do you find that a larger IV set generally leads to more robust causal estimates, or does it introduce more complexity and potential for bias with invalid instruments? 

3) Is the process of constructing the IV set dependent on the order in which instruments are considered? Specifically, does sequentially adding IVs versus a simultaneous assessment of all potential IVs impact the validity and effectiveness of the identified set?

4) Can the proposed algorithm handle large-scale datasets efficiently? What are the computational complexities and potential bottlenecks?
limitations: limitations: See weaknesses.
paper_summary: paper_summary: The authors take up a _very useful_ topic, of trying to identify instruments in models where bidirectional adjacencies exist, at least for the Mendelian randomization application.
strength_weakness: strength_weakness: From _my_ perspective, there were no glaring weaknesses to this paper. Perhaps other reviewers have issues to mention.

The only possible weakness I saw was the strong reliance on the assumption of linearity, though in the discussion this was mentioned as an assumption that could possibly be relaxed in future work.
questions: questions: No particular questions.
limitations: limitations: I did not see a discussion of societal impact.
paper_summary: paper_summary: The paper addresses the challenge of estimating causal effects in bi-directional Mendelian randomization (MR) studies using observational data, where invalid instruments and unmeasured confounding are common. It investigates theoretical conditions for identifying valid instrumental variable (IV) sets and proposes a cluster fusion-like algorithm to discover these IV sets and estimate causal effects accurately. Experimental results demonstrate the effectiveness of the method in handling bi-directional causal relationships, providing insights crucial for improving causal inference in complex systems.
strength_weakness: strength_weakness: The paper has no major weaknesses. However, the setting is restricted with causal relations limited to being linear and assuming genetic variants are randomized, which limits the practical applicability of the proposed approach. Additionally, the experimental results provided are mainly synthetic in nature.
questions: questions: I have few Questions/Suggestions for the Authors:

* In line 106, the authors mention, "Following Hausman [1983], we assume that $\beta_{X ->Y} \beta_{X->X} \neq 1$." It would be useful to discuss in the main paper why this assumption is necessary and what happens when it is violated.

* Similarly, regarding Assumption 3, the authors mention it as a very natural condition that one expects to hold for the unique identifiability of valid IVs. It would be useful to explain briefly in the main paper why this assumption is necessary for the identifiability of IVs.

* The authors in Section 5 claim that with dependence between genetic variants, main results may still be effective in identifying valid IV sets. Does this claim still hold when there is confounding among the genetic variants or between the genetic variants and some phenotype? Or does the dependence just mean direct causal effect here?. 

* At the moment, the proposed solution is restricted to linear causal relationships. Can one apply the proposed method using some linearization technique for scenarios where causal relationships are not necessarily linear?
limitations: limitations: The authors clearly state all the assumptions. The paper could benefit from adding some more discussion on the necessity of these assumptions in the main paper. I don't think the paper has any potential negative societal impacts.

Comments: title: Thank all reviewers again for greatly improving our paper!
Thank you for the responses.
Thanks.
Thanks for your responses
Re.
comment: Dear all reviewers,

We sincerely appreciate all your positive comments! We are grateful for your valuable and inspiring suggestions, which are of great help in improving our paper!

Best wishes,

The authors
Thank you for your detailed responses. I appreciate the clarification and will maintain my positive score.
Thanks for the rebuttal. I will stick to my original assessment.
Thanks the authors for your detailed responses. The extra experiments and discussions will be very helpful. I am happy to keep my positive rating.
Thanks for responding to my questions. The suggested changes by the authors, including additional experiments and clarifications, would be very beneficial for the paper. I will keep my decision and score for the paper.