Functional connectivity (FC) has become a primary means of understanding brain function by identifying interactions and, ultimately, how those interactions produce cognitions. Various methods estimate the association network to represent the brain FC, which could be problematic since they only provide spatial connections but not the real causal interactions. Hence, it is necessary to switch the focus from associations to causations. Directed acyclic graph (DAG) models, also known as Bayesian networks, are commonly used to model causal relationships in complex systems. Existing methods have focused on estimating a single directed graphical model. However, in many brain studies, we have data from related classes, such as different developmental stages and different disease states. Regarding statistical models, this corresponds to jointly estimating multiple DAGs under distinct but related conditions. Therefore, we propose a Bayesian incorporated linear Non-Gaussian Acyclic Model (BiLiNGAM) and apply it to the fMRI images from the Philadelphia Neurodevelopmental Cohort (PNC), which include 855 individuals aged 8–22 years who were divided into five adolescence-related stages.