This chapter introduces a novel bi-level approach for uncovering latent variables in Bayesian networks (BNs). The proposed methodology addresses the dual challenges of data completion by discovering latent variables and BN optimization. At the upper level, an estimation of the values of the hidden variable is performed using the expectation-maximization algorithm. Simultaneously, the lower level employs an evolutionary process to discover the optimal BN structure that represents the underlying relationships within the generated data by using the minimum description length principle. The upper-level and lower-level tasks, forming a hierarchical framework, interact through a leader-follower Stackelberg game. Our approach provides a comprehensive solution for latent variable identification within BNs. Through empirical evaluations, we demonstrate the efficacy of our proposed method in extracting latent variables in many datasets from real-world scenarios.