Computer Science
Computer Science-2021.gz:143324
Computer Science
7e9a7d611587950c7aed7fd31f92d547b93136bc
Simultaneous Missing Value Imputation and Structure Learning with Groups

Learning structures between groups of variables from data with missing values is an important task in the real world, yet difficult to solve. One typical scenario is discovering the structure among topics in the education domain to identify learning pathways. Here, the observations are student performances for questions under each topic which contain missing values. However, most existing methods focus on learning structures between a few individual variables from the complete data. In this work, we propose VISL, a novel scalable structure learning approach that can simultaneously infer structures between groups of variables under missing data and perform missing value imputations with deep learning. Particularly, we propose a generative model with a structured latent space and a graph neural network-based architecture, scaling to a large number of variables. Empirically, we conduct extensive experiments on synthetic, semi-synthetic, and real-world education data sets. We show improved performances on both imputation and structure learning accuracy compared to popular and recent approaches.


Introduction
Understanding the structural relationships among different variables provides critical insights in many real-world applications, such as medicine, economics and education (Sachs et al., 2005;Zhang et al., 2013). However, it is commonly impossible to perform randomized controlled trials for many real-world applications due to ethical or cost considerations. Thus, learning graphs from observed data, known as structure learning, has recently made remarkable progress (Fatemi et al., 2021;Yu et al., 2019;Zheng et al., 2018;2020 For many applications, variables in the data can be gathered into semantically meaningful groups, where useful insights are at group level. For example, in finance, one may be interested in how a financial situation influences different industries (i.e. groups) instead of individual companies (i.e. variables). Similarly, in education, the data can contain student responses to thousands of individual questions (i.e. variables), where each question belongs to a broader topic (i.e. groups). Again, it is insightful to find relationships between topics instead of individual questions. Moreover, real-world data such as educational data is inherently sparse since it is not feasible to ask every question to every student; the dimensions of the data in terms of the number of variables and the number of observations are very high, posing a scalability challenge. Despite the progress in structure learning, no existing method can discover group-wise relationships given large-scale partially observed data.
In this work, we present VISL (missing value imputation with structural learning), a novel approach to simultaneously tackle group-wise structure learning and missing value imputations driven by the real-world topic relationship discovery in an education setting. This is accomplished by combining variational inference with a generative model that leverages a structured latent space and a decoder based on message-passing Graph Neural Networks (GNN) (Gilmer et al., 2017). Namely, the structured latent space endows each group of variables with its latent subspace, and the interactions between the subspaces are regulated by a GNN whose behavior depends on the inferred graph from variational inference, see Fig. 1(a). VISL satisfies all the desired properties: it leverages continuous optimization of the structure learning to achieve scalability (Zheng et al., 2018;2020); the VISL formulation naturally handles missing values, and it can discover relations at different levels of granularity with pre-defined groups. Empirically, we evaluate VISL on one synthetic and two real-world problems including the aforementioned education scenario. VISL shows improved performance in both missing data imputation and structure learning accuracy compared to popular and recent approaches for each task. We worked closely with an education domain expert to evaluate the learned topic relationships, and our model has provided insightful results as recognized by the domain experts. VISL is a variational auto-encoder based framework. Observations from each group are encoded into low dimensional latent variables. The structure is treated as a global latent variable. A GNN based decoder is used to decode the latent variables to observations. (b) Probabilistic graphical model for VISL, where the partial observation x is generated from its local latent variable z and the global latent variable G which characterizes the structures.

Model Description
In the following, we present the formulation of VISL for scalable group-wise structure learning with partial observations using a novel deep generative model based framework.

Problem setting
Assume a training data set X = {x n } N n=1 with x n ∈ R D . The observed and missing values are denoted as X O and X U , respectively, where we assume the data are missing completely at random (MCAR) or missing at random (MAR). In Appx. A, we explain how to handle MAR. In particular, variables can be gathered into M pre-defined groups, where each can be denoted as χ n,m = [x n,i ] i∈Im . I m containing the variable indices belonging to group m (e.g., I 2 = [4, 5, 6] indicates group 2 includes the 4 th , 5 th and 6 th variables). One should note that each I m may have varying sizes for different m (i.e. varying group sizes). The goal of VISL is to (i) perform missing value imputation for test samples and (ii) infer structures between groups of variables. We use the adjacency matrix G ∈ [0, 1] M ×M to represent a graph, where G ij = 1 or 0 indicates whether there is a directed edge from i−th to j−th group or not. In the context of the education domain, the above formulation can be rephrased as follows: variable x n containing the student's responses to a set of questions. x i,j = 1 represents student i has answered question j correctly. Groups can be defined as the topic associated with each question. I m contains the question IDs that belong to the same topic, and χ m represents a group of responses related to that topic. Clearly, not all students can answer every question. Thus, X O , X U represent the existing responses and un-answered questions, respectively. The goal of VISL is to (i) predict students' responses to un-answered questions, which by itself is important in the education domain (Wang et al., 2020;2021), and (ii) discover the relationships between topics, which can help education experts optimize the learning experience and the curriculum. For structure learning, we adopt a Bayesian approach for graphs (Heckerman et al., 2006). Namely, we seek to maximize the posterior probability of G given partially observed training data X O within the space of all DAGs: G = arg max G∈DAGs p(X O |G)p(G). (1) To optimize over the structure with the DAG constraint in Eq. 1, we resort to recent continuous optimization techniques (Kyono et al., 2020;Zheng et al., 2018;2020), In the following two sections, we present our detailed formulation, training and imputation algorithms of VISL, that allows the model to infer the latent structure G and impute missing valuesx U in a test samplex ∈ R D based on the observedx O .

Generative model and variational inference
Algorithm 1 Generative process Gij ∼ Bernoulli(pij) for n ∈ {1, 2, · · · , N } do Zn ∼ N (0, σ 2 z I) xn ∼ N (f θ (Zn, G), σ 2 x I) end for For the generation of observation X, we adopt the latent variable model of Fig. 1. Particularly, given an inferred graph G and latent Z, the generative path from Z to X is provided in Fig. 1, where we use a graph neural network (GNN) decoder that respects the learned graph structure G and the provided grouping structure. Then the joint model likelihood is Amortized variational inference. The true posterior distribution over Z and G in Eq. 3 is intractable since we use a complex deep learning architecture. Therefore, we resort to an efficient amortized variational inference as in Kingma & Welling (2013); Kingma et al. (2019). Here, we consider a fully factorized variational distribution q(Z, G) = q φ (G) N n=1 q φ (Z n |x n ), where q φ (Z n |x n ) is a Gaussian whose mean and (diagonal) covariance matrix are given by an encoder. For q(G), we consider the product of independent Bernoulli distributions over the edges; that is, the presence of each edge from i to j is associated with a probability p ij to be estimated. With the above formulation, the evidence lower bound (ELBO) is Next, we explain our choice of the generator (decoder), which uses a GNN over a learned graph G to model the interactions between latent variables, representing the information about each group. Then, we focus on the inference network (encoder), representing the mapping from the group of observed variables to its corresponding latent representation.
Generator. The generator (i.e., decoder) takes Z n and G as inputs and outputs the reconstructedx n = f θ (Z n , G), where θ are the decoder parameters. In order to respect the pre-defined group structure, as shown in Fig. 1, Z n is partitioned into M parts, where z n,m represents the latent variable for the group of observations χ n,m . This defines a group-wise structured latent space. We adopt a two-step process for the generative path Z n to X n : (i) GNN message passing with respect to the learned graph G between latent z n,m ; (ii) final read-out layer to generate X n .
GNN message passing in the generator. In message passing, the information flows between nodes in T consecutive node-to-edge (n2e) and edge-to-node (e2n) operations (Gilmer et al., 2017). At the t-th step, we compute an embedding h f i→j for each edge i → j, called forward embedding, which summarizes the information sent from node i to j. Specifically, the n2e/e2n operations in VISL are Here, t refers to the t-th iteration of message passing (that is, Z (0) = Z n , notice that we omit subindex n for clarity).
Finally, MLP f , and MLP e2n are MLPs to be trained.
Interestingly, the message passing updates indicate that the information flows between latent nodes if a directed edge is specified in graph G. Hence, the inferred structure G directly defines relations for latent space Z which contains the information of pre-defined groups. We show that under certain conditions, the inferred graph G also represents the group-wise structure in observational space, and the corresponding model can be reformulated to a general structural equation model (SEM) (Peters et al., 2017) (see Appx. B).
Read-out layer in the generator. After T iterations of GNN message passing, we have Z (T ) . We then apply a final function that maps Z (T ) to the reconstructedx, which also respects the pre-defined group structure. Since the observation x = [χ 1 , . . . , χ M ] may contain χ m with different dimensions, we adopt M different MLPs, one for each group as the final read-out layer, to respect the group structure. Namely,x = (g 1 (z T 1 ), . . . , g M (z T M )) , where g m represents the MLP for group m. Thus, the decoder parameters θ include the parameters of the following neural networks: MLP f , MLP e2n and g m for m = 1, . . . , M .
Inference network. As in standard VAEs, the encoder maps a sample x n to its latent representation Z n . As discussed before, Z n is partitioned into M parts, where each z n,m contains the information of the observation in group m. Similar to the read-out layer, we utilize the M MLPs to map groups of observations to the mean/variance of the latent variables: Here, µ m φµ m and σ m φσ m are neural networks for group m. When missing values are present, we replace them with a constant as in (Nazabal et al., 2020). A graphic representation of how the encoder respects the structure of the latent space is shown in the appendix, Fig. 6(b).

Training VISL
Given the model described above, we propose the training objective to minimize w.r.t. θ, φ and G: where ELBO is given by Eq. 4 and the DAG regulariser R(G) was introduced in Eq. 2 to favor the DAG-ness of learned graph G.
Evaluating the training loss L VISL . VISL can work with any type of data.
The log-likelihood term (log p θ (x n |Z n , G) in Eq. 4) is defined according to the data type. We use a Gaussian likelihood for continuous variables and a Bernoulli likelihood for binary ones. For the inference of Z and G, the standard reparametrization trick is used to sample Z n from the Gaussian q φ (Z n |x n ) (Kingma & Welling, 2013;Kingma et al., 2019). To backpropagate the gradients through the discrete variable G, we resort to the Algorithm 2 Training VISL.
Input :Training dataset X, possibly with missing values. for each batch of samples {xn}n∈B do Drop a percentage of the data for each sample xn. Encode xn through the reparametrization trick to sample Zn ∼ N (µ φ (xn), σ 2 φ (xn)) using Eq.7. Use the Gumbel-softmax to sample G from q(G). Use decoder to reconstructxn = f θ (Zn, G). Calculate the training loss LVISL (Eq. 8). Gradient step w.r.t. φ (encoder parameters), θ (decoder parameters) and G (posterior edge probabilities). Output :Encoder parameters φ, decoder parameters θ, and posterior probabilities over the edges G.
Gumbel-softmax trick to sample from q(G) (Jang et al., 2017;Maddison et al., 2017). The KL[q φ (Z n |x n )||p(Z n )] and KL[q(G)||p(G))] terms can be obtained in closedform since they are Gaussian distributions and independent Bernoulli distributions over the edges, respectively. This formulation brings additional advantages in real-life applications since one can easily incorporate domain knowledge and prior information into the VISL framework. For example, if the existence of a specific edge is known a priori, the edge probability can be set to 0/1 in the prior distribution. Finally, the DAG-loss regulariser in Eq. 8 can be computed by evaluating the function R on a Gumbel-softmax sample from q(G). To adapt the model to different missing levels in the training data X, we adopt the masking strategy (Ma et al., 2019;Gong et al., 2019), which drops a random percent of the observed values during training. The entire training procedure for VISL is summarised in Algorithm 2.
Two-step training. After training, we obtain the posterior of the graph G, which respects the underlying structure of the groups as shown in Appx. B. With the trained network, we can impute missing values in the groups where their ancestors contain some observations but if a group has no ancestors no information can be propagated during imputation. After learning the graph structure and to facilitate the imputation task, we introduce a backwards edge: for an edge from j to i we denote the backwards edge information as h b i→j which codifies the information that the i → j edge lets flow from the j-th to the ith node. It is defined in the same way as Eq. 5, i.e.,: h where MLP b is the backward MLP; and the e2n update (Eq.6) is modified to z In summary, we propose a two-stage training process, where the first stage -described in previous sections -focuses on discovering the edge directions between nodes without the MLP b (i.e., we do not train the MLP b ). In the second stage, we fix the graph structure G and continue to train the model with the backward MLP. This two-stage training process allows VISL to leverage the backward MLP for the imputation task without updating the graph structure.
Revisiting the learning objectives. The optimal graph of relationships, denoted as G in Eq. 2, is given by the estimated posterior probabilities of graph G. In addition, the regularizer R(G) provides a way to evaluate if the resulting graph is a DAG. By tuning the regularizer strength λ, one can ensure that the resulting G * represents a proper DAG.
Therefore, the distribution over x U (missing values) is obtained by applying the encoder and decoder with x as input. One important distinction of VISL compared to Ma et al. (2019);Nazabal et al. (2020) is that it incorporates the learned structure G into the imputation, which helps the model avoid over-fitting due to spurious correlations (Kyono et al., 2020).
Special case: variable-wise relations. In the above formulation, we have defined VISL for group-wise structure learning. Variable-wise relations can be regarded as a special case. In particular, we can set M = D and I m = {m} (see Fig. 5 (a) in the appendix), i.e. each group only contains a single variable. Through this modification, we can further simplify the encoder and read-out layer. Instead of using M different MLPs, a single MLP can be shared across all variables since each group has dimension of 1. The mean function for the encoder is then defined as µ n = (µ φ (x n,1 ), . . . , µ φ (x n,D )) .
One can define encoder variance σ ( Fig. 6 (a) in the appendix) and the read-out layer g analogously.

Related Work
Since VISL simultaneously tackles missing value imputation and structure learning, we review both fields. Moreover, we review recent works that utilize structure learning to improve the performance of another deep learning task, similar to VISL. Finally, as one of the focused applications of this work is in the education domain, we review recent advances of AI in education.
Structure learning. Structure learning aims to infer the underlying structures associated with some observations. There are mainly three types of methods: constrained-based, score-based, and hybrid. Constraint-based ones exploit (conditional) independence tests to find the underlying structure, such as PC (Spirtes & Glymour, 1991) and Fast Causal Inference (FCI) (Spirtes et al., 2000). They have recently been extended to handle partially observed data through testwise deletion and adjustments (Strobl et al., 2018;Tu et al., 2019a). Score-based methods find the structure by optimizing a proper scoring function. The core difficulty lies in the number of possible graphs growing super-exponentially with the number of nodes (Chickering et al., 2004 Fatemi et al., 2021), the classification objective is supplemented with a self-supervised task that learns a graph of interactions between variables through a GNN. However, these works focused on the supervised classification task, and they did not advance the performance of the structure learning itself.
Missing values imputation. The relevance of missing data in real-world problems has motivated a long history of research (Dempster et al., 1977;Rubin, 1976 to perform missing value imputation and personalization. However, partial VAE does not consider the structural relations between questions/topics and cannot perform structure learning. With the additional insights from structure learning, VISL can provide more information to teachers to help curriculum design than just imputations.

Experiments
We evaluate the performance of VISL in three different problems: a synthetic experiment where the data generation process is controlled, a semi-synthetic problem (simulated data from a real-world problem) with many more variables (Neuropathic Pain), and the real-world problem that motivated the development of the group-level structure learning (Eedi). To compare the model performance with related work, the first two datasets are on the variable level, which means that each group only has one variable. In the education setting, we focus on the real-world usage of the method and have worked closely with the domain expert to evaluate the results. Additional experiments are presented in the appendix.
Baselines. We consider five baselines for the structure discovery task at the variable level. PC (Spirtes et al., 2000) and GES (Chickering, 2002) are the most popular methods in constraint-based and score-based approaches, respectively. We also consider three recent algorithms based on contin-   Figure 3. Probability of edges obtained by VISL in the synthetic experiment. By using a 0.5 threshold, we get the predicted graph in Fig. 2    . Unlike VISL, these baselines cannot deal with missing values in the training data. Therefore, we work with fully observed training data in the first two sections when using these baselines. In contrast, the real-world data in the last section comes with partially observed training data, and the goal is to discover group-wise relationships. These baselines are not applicable. For the missing data imputation task, we also consider five baselines. Metrics. Imputation performance is evaluated with standard metrics such as RMSE (continuous variables) and accuracy (binary variables). For binary variables , we also provide the area under the ROC and the Precision-Recall curves (AU-ROC and AUPR, respectively), which are especially useful for imbalanced data (such as Neuropathic Pain). We follow common practice (Glymour et al., 2019; Tu et al., 2019a) regarding structure discovery performance, and consider metrics on the adjacency and orientation. While the former does not take into account the direction of the edges, the latter does. For both adjacency and orientation, we compute recall, precision and F 1 -score. We also provide causal accuracy, a popular structure discovery metric that considers edge orientation (Claassen & Heskes, 2012).

Synthetic experiment
We simulate fifteen synthetic datasets. For each simulated dataset, we first sample the true structure G; see Fig. 2(a) for an example. We obtain the samples of the datasets by computing each variable based on its parents using a nonlinear mapping based on the sin function. The appendix provides further details, including a visualisation of the generated data in Fig. 7. For each dataset, we simulate 5000 training and 1000 test samples.
Imputation performance. VISL outperforms the baselines in terms of imputation across all synthetic datasets (  Table 1).The results grouped by the number of variables are presented by Table 8 in the appendix. Therefore, VISL exploits the learned graph to improve imputation by avoiding spurious correlations.
Structure discovery performance. VISL obtains better performance than the baselines, see Table 2. The results split by the number of variables are shown in the appendix, Table 10. Notice that NOTEARS (NL) is slightly better in terms of orientation precision. However, this is at the expense of a significantly lower capacity to detect true edges; see the recall and the trade-off between both (F 1 -score). In this small synthetic experiment, it is possible to visually inspect the predicted graph. Fig. 3 shows the posterior probability of each edge (i.e. the estimated matrix G) for the simulated dataset that uses the true graph in Fig. 2(a). Using a threshold of 0.5, we obtain the predicted graph in Fig. 2(b). We observe that all the true edges are captured by VISL, with some additional edges due to finite data and non-convex optimization.
Finally, VISL can scale to large data both in terms of data  Table 3. Structure discovery results for neuropathic pain data (mean and std error over five runs).
points (which benefits naturally from its SGD-based optimization) and dimensionality (thanks to the continuous optimization over the graph space). We demonstrate the computational efficiency with synthetic data ranging from 4 nodes to 512 nodes in the appendix,  Table 4. Imputation results for neuropathic pain data (mean and std error over five runs).
We evaluate our method using a machine learning benchmark in healthcare applications (Tu et al., 2019b). The dataset contains records of patients regarding the symptoms associated with neuropathic pain. There are 222 variables in this dataset. Unlike the previous experiment with continuous data, this dataset has binary variables indicating the symptoms. The train and test sets have 1000 and 500 patients, respectively.
Imputation performance. VISL shows competitive or superior performance when compared to the baselines, see Table 4. Notice that AUROC and AUPR allow for an appropriate threshold-free assessment in this imbalanced scenario. Indeed, as expected from medical data, the minority of values are 1 (symptoms); here, the prevalence of symptoms is around 8% in the test set. Interestingly, it is precisely in AUPR where the differences between VISL and the rest of the baselines are larger except MICE, whose performance is very similar to that of VISL in this dataset.
Structure discovery results. As in the synthetic experiment, VISL outperforms the causality-based baselines; see Table 3. Notice that NOTEARS (NL) is slightly better in terms of adjacency-precision, i.e. the edges that it predicts are slightly more reliable. However, this is at the expense of a significantly lower capacity to detect true edges, see the recall and the trade-off between both (F 1 -score).

Eedi topics dataset
Finally, we evaluate our method on an even more challeng-  Table 5. Imputation results for Eedi topics dataset (mean and standard error over five runs).
ing real-world dataset in education requiring group-wise structure discovery. This is an important real-world problem in the field of AI-powered educational systems (Wang et al., 2021;2020). In this setting, we are interested in relationships between topics while the observations are questionanswer pairs under these topics. The dataset is very sparse, with 74.1% of the values missing. The dataset contains the responses by 6147 students to 948 mathematics questions. The 948 variables are binary (1 if the student provided the correct answer and 0 otherwise). These 948 questions target very specific mathematical concepts and are grouped within a meaningful hierarchy of topics; see Fig. 4. Here we apply our proposed model to find the relationships among the topics using the third level of the topic hierarchy (Fig. 4), resulting in 57 group-level nodes.  Table 6. Average expert evaluation of the topic relationships. Cohen's κ inter-annotator agreement is 0.72 for adjacency and 0.76 for orientation (substantial agreement).
Imputation results. VISL achieves competitive or superior performance when compared to the baselines (Table 5).
Although the dataset is relatively balanced (54% of the values are 1), we provide AUROC and AUPR for completeness. Notice that this setting is more challenging than the previous ones since we learn relationships between groups of variables (topics). Indeed, whereas the group extension allows for more meaningful relationships, the information flow happens at a less granular level. Interestingly, even in this case, VISL obtains similar or improved imputation results compared to the baselines.
Structure discovery results between groups. Most of the baselines used so far cannot be applied here because i)  Figure 4. Hierarchy of topics in Eedi data. All the questions are related to maths (level 0 topic). The number of topics at levels 1, 2 and 3 are 3, 25 and 57. Each question is associated with only one topic at level 3 (thus, to only one topic at any higher level).  To cope with arbitrary groups, we apply the group-specific mappings described for VISL. Finally, as an additional reference, we also compare with randomly generated relationships, which we refer to as Random.
Moreover, as there are no ground truth relationships in this real-world application, we ask two experts (teachers) to assess the validity of the relationships found by VISL, DAG-GNN, and Random. For each relationship, they evaluate the adjacency (whether it is sensible to connect the two topics) and the orientation (whether the first one is a prerequisite for the second one). They provide an integer value from 1 (strongly disagree) to 5 (strongly agree), i.e. the higher, the better. The complete list of relationships and expert evaluations for VISL, DAG-GNN, and Random can be found in the appendix; see Table 15, Table 16, and Table 17, respectively. In summary, Table 6 shows here the average evaluations: we see that the relationships discovered by VISL score much more highly across both metrics than the baseline models.
Another interesting aspect is how the relationships between level-3 topics are distributed across higher-level topics (recall Fig. 4). Intuitively, it is expected that most of the relationships happen inside higher-level topics (e.g. Numberrelated concepts are more probably related to each other than to Geometry-related ones). Table 7 shows such a distribution for the compared methods. Indeed, notice that the percentage of inside-topic relationships is higher for VISL (82%) and DAG-GNN (42%) than for Random (34%). An analogous analysis for the 25 level-2 topics is provided in the appendix; see Education Impact. Lastly, to make a real-world impact, we have been provided with an additional education dataset in the same format as Eedi by an education organization to help provide insight for math curriculum building. The final structure among all topics found by VISL is presented by figure Fig. 8 in the appendix. The predicted relationships allowed insights into which topics are foundational and need to be covered earlier (topics with many originating edges), as well as which topics are more complex and should be covered later (topics with many incoming edges). This allowed us to re-evaluate the order of topics in a nationwide used secondary curriculum. Specifically, topics such as "arithmetic" or "properties of shapes" were moved earlier in the curriculum, while topics such as "negative numbers" or "proportion and similarity" were moved to a later stage in the curriculum. Another interesting example found by the domain expert is the Venn diagram, which was originally taught in year 9/10 and is now suggested to move to year 7. Experts found that the Venn diagram has been a useful tool in teaching other topics which are currently taught before year 10. Moving this earlier will help students to learn other topics better. This emphasises the real-world impact our model, VISL, can have in planning curricula.

Conclusions
We introduced VISL, a novel approach that simultaneously performs group-wise structure discovery and learns to impute missing values. Both tasks are performed jointly: imputation is informed by the discovered relationships and vice-versa, leading to improved performance for both tasks. Moreover, motivated by a real-world problem, VISL shows its impact in the real-world education domain to aid domain experts in setting up curriculum. Fatemi, B., Asri, L. E., and Kazemi, S. M. Slaps: Selfsupervision improves structure learning for graph neural networks. arXiv preprint arXiv:2102.05034, 2021.

A. VISL for missing-at-random scenario
In this section, we briefly explain why VISL can handle MAR problem by leveraging the results from Rubin (1976). Let's denote r i as the missing mask, where r i,d = 1 indicates x i,d is observed. For the random variable R, we use p λ (r|x) as the missing mechanism with parameters λ. To explicitly state the dependence of VISL and its model parameter θ, we use p θ (x) to denote the corresponding model density. We can now formally define the concept of MAR. Definition A.1 (Missing at Random (Rubin, 1976)). The missing data are missing at random if for each value of λ, p λ (r|x) takes teh same value for all x u . Namely, p λ (r|x) = p λ (r|x o ).
Recall that our VISL is trained by maximizing the ELBO (Eq.4) based on the observed values x o . However, this formulation ignores the missing mechanisms p φ (r|x). In order to perform missing value imputation, one need to ensure that the inference for θ should be correct. In the following, we show under MAR, ignoring missing mechanism does not affect the correctness of inferring θ under ELBO. The following proof is an adaptation of Theorem 7.1 in Rubin (1976).
When explicitly modelling the missing mechanism, the joint likelihood can be written as where the third equality is from the definition of MAR and the last inequality is from the standard ELBO derivation. The above equation explicitly lower bounds the joint likelihood by two separate terms regarding θ and φ. Thus, when performing inference over θ, one can safely ignore the missing mechanism involving φ, resulting in the same optimization objective as Eq.4.

B. Does VISL respect the graph G in observational space
From the formulation of the decoder in VISL (Eq.5 and 6), the inferred graph G seems to define whether the information flow between nodes is allowed or not. Namely, when G ij = 1, the information is allowed to pass from z i to z j at each iteration t. Thus, G directly defines a structure for latent space Z, and indirectly defines a structure in observation X through the GNN updates and the final read-out layer. A natural question to ask is whether the resulting observations x from VISL also respect the graph G. In the following, we show that when GNN is in equilibrium and the read-out layer is invertible without additional observational noise (σ x = 0), the VISL is in fact a SEM for observation x, which respects the graph G.
In the following, for the clarity of notations, we consider the structure learning between variables. For group-wise relations, it is trivial to generalize, since the going from variable-wise to group-wise only changes the read-out layer, where we use M different MLPs instead of one.
First, let's clarify what do we mean by "respect a graph G". Definition B.1 (Respect a graph G). For a given VISL model p(x, z; G) with a specific graph G, we say the model p(x; G) = p(x, z; G)dz respects the graph G if it can be factorized where P A(d) is a set of parents of node i specified by graph G.

B.1. GNN at steady state
From the GNN message passing equations, we can re-organize Eq.5 and 6 into one equation: where P A(i) t−1 is a set of parents' value for node i at iteration t − 1, z t i is the value for node i at iteration t and F (·) represents the GNN message passing updates.
The above equation resembles a fixed-point iteration procedure for function F . Indeed, under the context of GNN, this has been considered as a standard procedure to search for the equilibrium state due to its exponential convergence (Dai et al. 2018, Eq.1;Gu et al. 2020, Eq.2(b)). Thus, we assume that the GNN updates F has unique equilibrium states given the initial conditions AN (i) 0 ∪ z 0 i for each i, where AN (i) 0 represents the initial values of ancestors for node i. For a sufficient condition of its existence, one can refer to Gu et al. (2020, Theorem 4.1). We note that this is only a sufficient condition, meaning that the GNN without the conditions in Gu et al. (2020) can still have equilibrium states. Since discussing a necessary and sufficient conditions for the existence of the equilibrium state is out of the scope of this paper, we simply made an assumption that function F has steady states. The reason we consider the initial ancestor values rather than just parent values is due to the message passing nature, where the value P A(i) t contains the information from the nodes that is at most t-hops away.
Since graph G represents a DAG, one can always find a permutation π of the original index i = 1, . . . , D based on the topological order. For concise notations, we assume the identity permutation. When the GNN is in equilibrium, we can rewrite Eq.11 as where the superscript ∞ represents that steady state of the node. From the assumption, since the steady state for each z ∞ i depends on the initial values AN (i) 0 ∪ z 0 i , it is trivial to see that the steady state P A ∞ (i) depends on AN 0 (i). Therefore, the steady state z ∞ i is uniquely determined by P A ∞ (i) and z 0 i . Namely, for i = 1, . . . , D, where H i is a mapping from P A(i) ∞ and z 0 i to the steady state of node i. This is exactly the general form of an structural equation model (SEM) defined by graph G. If we further assume that the read-out layer g(·) is invertible, we can obtain which is also an SEM based on G for observation x. Thus, it is trivial that VISL respects the graph G based on the above assumptions.
In practice, due to the exponential convergence of fixed point iteration, we found out that one does not need to use large iteration t. To balance the performance and computational cost, we found that 3 iterations of GNN message passing is enough to obtain reasonable performances.

C. Experimental details
Here we specify the complete experimental details for full reproducibility. We first provide all the details for the synthetic experiment. Then we explain the differences for the neuropathic pain and the Eedi topics experiments.

C.1. Synthetic experiment
Data generation process. To understand how the number of variables affects VISL, we use D = 5, 7, 9 variables (five datasets for each value of D). We first sample the underlying true structure. An edge from variable i to variable j is sampled with probability 0.5 if i < j, and probability 0 if i ≥ j (this ensures that the true structure is a DAG, which is just a standard scenario, and not a requirement for any of the compared algorithms). Then, we generate the data points. Root nodes (i.e. nodes with no parents, like variables 1 and 2 in Fig. 2(a) in the paper) are sampled from N (0, 1). Any other node v i is obtained from its parents Pa(i) as v i = j∈Pa(i) sin(3v j ) + ε, where ε → N (0, 0.01) is a Gaussian noise. We use the sin function to induce non-linear relationships between variables. Notice that the 3-times factor inside the sin encourages that the whole period of the sin function is used (to favor non-linearity). To evaluate the imputation methods, 30% of the test values are dropped. As an example of the data generation process, Fig. 7 below shows the pair plot for the dataset generated from the graph in Fig. 2(a) in the paper.
Model parameters. We start by specifying the parameters associated to the generative process. We use a prior probability p ij = 0.05 in p(G) for all the edges. This favours sparse graphs, and can be adjusted depending on the problem at hand.
The prior p(Z) is a standard Gaussian distribution, i.e. σ 2 z = 1. This provides a standard regularisation for the latent space. The output noise is set to σ 2 x = 0.02, which favours the accurate reconstruction of samples. As for the decoder, we perform T = 3 iterations of GNN message passing. All the MLPs in the decoder (i.e. MLP f , MLP b , MLP e2n and g) have two linear layers with ReLU non-linearity. The dimensionality of the hidden layer, which is the dimensionality of each latent subspace, is 256. Regarding the encoder, it is given by a multi-head neural network that defines the mean and standard deviation of the latent representation. The neural network is a MLP with two standard linear layers with ReLu non-linearity. The dimension of the hidden layer is also 256. When using groups, there are as many such MLPs as groups. Finally, recall that the variational posterior q(G) is the product of independent Bernoulli distributions over the edges, with a probability G ij to be estimated for each edge. These values are all initialised to G ij = 0.5.
Training hyperparameters. We use Adam optimizer with learning rate 0.001. We train during 300 epochs with a batch size of 100 samples. Each one of the two stages described in the two-step training takes half of the epochs. The percentage of data dropped during training for each instance is sampled from a uniform distribution. When doing the reparametrization trick (i.e. when sampling from Z n ), we obtain one sample during training (100 samples in test time). For the Gumbel-softmax sample, we use a temperature τ = 0.5. The rest of hyperparameters are the standard ones in torch.nn.functional.gumbel softmax, in particular we use soft samples. To compute the DAG regulariser R(G), we use the exponential matrix implementation in torch.matrix exp. This is in contrast to previous approaches, which resort to approximations (Zheng et al., 2018;Yu et al., 2019). When applying the encoder, missing values in the training data are replaced with the value 0 (continuous variables).
Baselines details. Regarding the structure learning baselines, we ran both PC and GES with the Causal Command tool offered by the Center for Causal Discovery https://www.ccd.pitt.edu/tools/. We used the default parameters in each case (i.e. disc-bic-score for GES and cg-lr-test for PC). NOTEARS (L), NOTEARS (NL) and DAG-GNN were run with the code provided by the authors in GitHub: https://github.com/xunzheng/notears (NOTEARS (L) and NOTEARS (NL)) and https://github.com/fishmoon1234/DAG-GNN (DAG-GNN). In all cases, we used the default parameters proposed by the authors. Regarding the imputation baselines, Majority Vote and Mean Imputing were implemented in Python. MICE and Missforest were used from Scikit-learn library with default parameters https://scikit-learn.org/stable/modules/generated/ sklearn.impute.IterativeImputer.html#sklearn.impute.IterativeImputer. For PVAE, we use the authors implementation with their proposed parameters, see https://github.com/microsoft/EDDI.
Other experimental details. VISL is implemented in PyTorch. The code is available in the supplementary material. The experiments were run using a local Tesla K80 GPU and a compute cluster provided by Azure Machine Learning platform with NVIDIA Tesla V100 GPU.

C.2. Neuropathic pain experiment
Data generation process. We use the Neuropathic Pain Diagnosis Simulator in https://github.com/TURuibo/ Neuropathic-Pain-Diagnosis-Simulator. We simulate five datasets with 1500 samples, and split each one randomly in 1000 training and 500 test samples. To evaluate the imputation methods, 30% of the test values are dropped. These five datasets are used for the five independent runs reported in experimental results.
Model and training hyperparameters. Most of the hyperparameters are identical to the synthetic experiment. However, in this case we have to deal with 222 variables, many more than before. In particular, the number of possible edges is 49062. Therefore, we reduce the dimensionality of each latent subspace to 32, the batch size to 25, and the amount of test samples for Z n to 10 (in training we still use one as before). Moreover, we reduce the initial posterior probability for each edge to 0.2. The reason is that, for 0.5 initialization, the DAG regulariser R(G) evaluates to extremely high and unstable values for the 222 × 222 matrix. Since this is a more complex problem (no synthetic generation), we run the algorithm for 1000 epochs. When applying the encoder, missing values in the training data are replaced with the value 0.5 (binary variables).
Data generation process. The real-world Eedi topics dataset contains 6147 samples, and can be downloaded from the website https://eedi.com/projects/neurips-education-challenge (task3 4 folder). The mapping from each question to its topics (also called "subjects") is given by the file "data/metadata/question metadata task 3 4.csv". For those questions that have more than one topic associated at the same level, randomly sample one of them. The hierarcy of topics (recall Fig. 4 in the paper) is given by the file "data/metadata/subject metadata.csv". We use a random 80%-10%-10% train-validation-test split. The validation set is used to perform Bayesian Optimization (BO) as described below. The five runs reported in the experimental section come from different (random) initializations for the model parameters.
Baselines details. As explained in the paper, in this experiment DAG-GNN is adapted to deal with missing values and groups of arbitrary size. For the former, we adapt the DAG-GNN code to replace missing values with 0.5 constant value, as in VISL. For the latter, we also follow VISL and use as many different neural networks as groups (as described in the paper), all of them with the same architecture as the one used in the original code (https://github.com/fishmoon1234/DAG-GNN).
Other experimental details. The list of relationships found by VISL (Table 15) and DAG-GNN (Table 16) aggregates the relationships obtained in the five independent runs. This is done by setting a threshold of 0.35 on the posterior probability of edge (which is initialized to 0.2) and considering the union for the different runs. This resulted in 50 relationships for VISL and 57 for DAG-GNN. For Random, we simulated 50 random relationships. Also, the probability reported in the first column of Table 15 is the average of the probabilities obtained for that relationship in the five different runs.     D = 5, 7, 9). The values are the mean and standard error over five different simulations.  Table 9. Mapping between indexes for row/column names in Table 12 and Table 14 D = 5, 7, 9). The values are the mean and standard error over five different simulations.  Table 11. Running times (in minutes) for different structure learning approaches in an extended synthetic experiment. For each number of variables, three datasets were simulated following the same data generation process described above, and the results show the mean and standard error. Notice that we have considered three additional baselines (MMHC, Tabu, HillClimb). We observe three different types of methods. MMHC scales poorly (NA means that the training took more than 24 hours), probably due to its hybrid nature that combines constraint-based and score-based approaches. VISL and the other deep learning based methods (DAG-GNN, NOTEARS) can scale to large numbers of variables. Of course, simpler methods such as PC, GES, Tabu, and HillClimb are significantly faster than VISL (note also that these baselines are from highly optimized libraries that leverage e.g. dynamic programming and parallelization), but their performance is worse. Indeed, the structure learning performance for this experiment is shown in Table 18.   Table 9 for a mapping between indexes shown here in row/column names and the actual level-2 topic names.  Table 13. How the 57 relationships found by DAG-GNN are distributed across level 2 topics. The item (i, j) refers to edges in the direction i → j. There are 8 relationships inside level 2 topics (14%). See Table 9 for a mapping between indexes shown here in row/column names and the actual level-2 topic names.  Table 14. How the 50 relationships found by Random are distributed across level 2 topics. The item (i, j) refers to edges in the direction i → j. There are 3 relationships inside level 2 topics (6%). See Table 9 for a mapping between indexes shown here in row/column names and the actual level-2 topic names.  Table 15. Full list of relationships found by VISL in the Eedi topics dataset. Each row refers to one relationship (one edge). From left to right, the columns are the posterior probability of the edge, the sending node (topic), the receiving node (topic), and the adjacency and orientation evaluations from each expert. For each topic, the brackets contain its parent level 2 and level 1 topics.  Figure 8. The structures discovered by VISL based on an alternative education dataset. In particular, the blue node represents the curriculum that should be taught later. On the other hand, the green node represents the fundamental topics that should be taught earlier.