Identification of Latent Variables From Graphical Model Residuals

Construction of graphical models (GMs) pursues two related objectives: accurate inference of conditional independencies (causal discovery), and construction of models for outcomes of interest with the purpose of estimating the average causal effect (ACE) with respect to interventions (Pearl 2000; Hernán and Robins 2006) (causal inference). These two distinct goals mandate that certain identifiability conditions for both types of tasks be met. Of particular interest to us is the condition of causal sufficiency ((Spirtes, Glymour, and Scheines 1993)) - namely, that all of the relevant confounders have been observed. When this condition is not met, accurate GMs and identifiability of the causal effects are hard to infer.

Intuitively, the presence of unobserved confounding leads to violations of conditional independence among the affected variables downstream from any latent confounder. Score-based methods for GM construction aim to minimize unexplained variance for all variables in the network by accounting for conditional independencies in the data (Pearl 2000; Friedman and Koller 2013). Unobservability of a causally important variable will induce dependencies among its descendants that are conditionally independent of each other given the latent variable. This gives rise to inferred connectivity that is excessive compared to the true network (Elidan et al. 2001). Further, since none of the observed descendants are perfect correlates of the unobserved ancestor, some of the information from the ancestor will "pollute" model residuals and give rise to correlation patterns in the residual space.

Hitherto, the causal inference literature has largely focused on addressing the subset of problems where the ACE can be estimated reliably in the absence of this guarantee, such as when conditional exchangeability holds (Hernán and Robins 2006). In causal discovery, certain constraint-based methods like the Fast Causal Algorithm (FCI) (Spirtes, Meek, and Richardson 1995) can deal with general latent structures but do not scale well to large problems. While faster variants of FCI like RFCI (Colombo et al. 2012) and FCI+ (Claassen, Mooij, and Heskes 2013) exist, score-based approaches tend to perform better than purely constraint-based methods. This has been demonstrated in problems without latent confounders (Nandy et al. 2018) and has been suggested that this is because constraint-based methods are sensitive to early mistakes which progagate into more errors later on while in score-based approaches the mistakes are localized and do not affect the score of graphs later on (Bernstein et al. 2020). It is likely that these limitations also affect FCI-type algorithms. Furthermore, while FCI can deal with general latent structures, in certain scenarios where a few latent confounders drive many observed variables most observed variables are conditionally dependent which implies dense maximal ancestral graphs where very few edges can be oriented (Frot, Nandy, and Maathuis 2017) (see supplementary Figure LABEL:sup-fig_rfci_ex1 for an illustration). Finally, another disadvantage of constraint based methods is that they generate a single model and do not provide a score summarizing how likely the network is given the data while score based methods can generate this type of scores by using Bayesian or ensemble approaches (Jabbari et al. 2017). Relatively computationally efficient score-based methods have been proposed for inferring latent variables affecting GMs by the means of expectation maximization (EM) (as far back as (Friedman 1997, 1998)). However, for a large enough network, local gradients do not provide a reliable guide, nor do they address the cardinality of the latent space. Methods for using near-cliques for detection of latent variables in directed acyclic graphs (DAGs) with Gaussian graphical models (GGMs) have been proposed that address both problems by analyzing near-cliques in DAGs (Elidan et al. 2001; Silva et al. 2006). A method related to ours has been proposed for calculating latent variables in a greedy fashion in linear and "invertible continuous" networks, and relating such estimates to observed data to speed up structure search and enable discovery of hidden variables (Elidan, Nachman, and Friedman 2007). However, with a clique-driven approach it is impossible to tell whether any cliques have shared parents and, importantly, whether any signal remained to be modeled, resulting in score-based testing rejecting proposed "ideal parents". Additionally, Wang and Blei (2019) introduced the deconfounder approach to detecting and adjusting for latent confounders of causal effects in the presence of multiple causes but in the context of a fixed DAG under relatively strict conditions (Wang and Blei 2019a, b).

We show that there exist circumstances when causal sufficiency can be asymptotically achieved and exchangeability ensured even when the causal drivers of outcome are confounded by a number of latent variables. This can be done when confounding is pleiotropic, that is, when the latent variable affects a "large enough" number of variables, some driving an outcome of interest and others not (we tie this to the expansion property defined in (Anandkumar et al. 2013)). Notably, this objective cannot be achieved when confounding affects only the variables of interest and their causal parents ((D’Amour 2019)). Our approach is provably optimal given a proposed graph. This leads to the iterative EM-like optimization approach over structure and latent space proposed below.

The outline of this paper is as follows: In the first two sections, we discuss the process of diagnosting latent confounding in the special case of GGMs where closed-form local solutions exist. In the third section, we derive the improvement cap for predictive performance of the graph thus learned. In the fourth section, we discuss implementation, and in the fifth demonstrate the approach on simulated and real data. Finally, we propose some extensions to the approach for more complex graphical model structures and functional forms.

Consider a factorized joint probability distribution over a set of observed and latent variables Table 1 summarizes the notation to be used to distinguish variables, their relevant partitions , and parameters.

Assume that the joint distribution $D_{u}$ over the full data or $D$ over the observed data is factorized as a directed acyclic graph, $G_{u}$ or $G$. We will consider individual conditional probabilities describing nodes and their parents, $P(V|Pa^{V},\theta)$, where $\theta$ refers to the parameters linking the parents of $V$ to $V$. $\hat{\theta}$ will refer to an estimate of these parameters. We will further assume that $G$ is constructed subject to regularization and using unbiased estimators for $P(V|Pa^{V},\hat{\theta})$. We will further assume that $D_{u}$ plus any given constraints are sufficient to infer the true graph up to Markov equivalence. For convenience, we’ll focus on the actual true graph’s parameters, so that, using unbiased estimators, $E[\hat{\theta_{m}}|D_{u}]=\theta_{m},\forall m$.

Mirroring $D$ (or $D_{u}$), we will define a matrix $R$ of the same dimensions ($s\times(v+o)$) that captures the residuals of modeling every variable $N\in\{V,O,(U)\}$ via $G$ (or $G_{u}$). Note that $R_{u}$ - the residuals matrix that contains residuals of $U$ - doesn’t make sense and isn’t defined. In the linear case, these would be regular linear model residuals, but more generally we will consider probability scale residuals (PSR, (Shepherd, Li, and Liu 2016)). That is, we define $R[i,j]=PSR(P(V_{j}|Pa^{V_{j}},\hat{\theta}_{j})|D[i,j])$, the residuals of $V_{j}$ given its graph parents. The use of probability-scale residuals allows us to define $R$ for all ordinal variable types, up to rank-equivalence.

Here, we consider the special case of GGMs where our approach for determining the existence of latent confounders can be written down in a closed form.

Recall that, for some $V_{j}\in\{V,U,O\}$, $P^{V_{j}}_{k}$ denotes the $k$th parent of $V_{j}$. For GGMs, we can write down a fragment of any DAG G as a linear equation where a child node is parameterized as a linear function of its parents:

For example, consider $O_{1}$ in Figure 2. We can write:

For any variable N that has parents in $G_{u}$, we can group variables in $P^{N}$ into three subsets: $X_{U}\in\{P^{U},C^{U}\}$, $X_{\cancel{U}}\notin\{P^{U},C^{U}\}$, and the set $U$ itself, and write down the following general form using matrix notation:

Explicit dependence of $N$ on $U$ happens when $\beta_{U}\neq 0$.

Note that if we deleted $U$ and its edges from $G_{u}$ without relearning the graph, Equation 3 from $G_{u}$ would read:

where the residual term $R_{N}$ is simply equal to the direct contribution of $U$ to $N$.

Now consider $G$, the graph learned over the variables $\{V,O\}$ excluding the latent space $U$. The network $G$ would have to adjust to the missingness of $U$ (e.g., Figure 2 vs Figure 2). As a result, $R_{N}$ will be partially substituted by other variables in $\{P^{U},C^{U}\}$. Still, unless $U$ is completely explained by $\{P^{U},C^{U}\}$ (as described in (D’Amour 2019)) and in the absence of regularization (when a high enough number of covariates may lead to such collinearity), $R_{N}$ will not fully disappear in $G$. Hence, even after partially explaining the contribution of $U$ to $N$ by some of the parents of $N$ in $G$,

Therefore, the columns in the residuals table corresponding to $G$ (Table 2) that represent the parents and children of $U$ will contain residuals collinear with $U$:

Rearranging and combining,

Equation 7 tells us that, for graphical Gaussian models, components of $U$ are obtainable from linear combinations of residuals, or principal components (PCs) of the residual table $R$. In other words, $U$ is identifiable by principal component analysis (PCA). Whether the residuals needed for this identification exist depends on the expansion property as defined in (Anandkumar et al. 2013).

Theorem 1 does require G to be known, which will give rise to an iterative algorithm later on. However, from it we see that focusing on residuals of a graph permits identification of the latent space in some circumstances. Supplementary theorem LABEL:sup-thm:deconfounding relaxes the need for orthogonality when latent space’s structure is irrelevant but controlling for it is important.

Note that this approach is similar, in the linear case, to that in (Elidan, Nachman, and Friedman 2007) except insofar as we show the principal components of the entire residual matrix to be optimal for the discovery of the whole latent space of a given DAG, and thus couple structure inference and latent space discovery via EM (1).

The most important utility of causal modeling is to identify drivers as well as drivers of drivers of outcomes of interest. These drivers of drivers may have practical importance. For example, in the development of a drug, direct predictors of an outcome (e.g. $V_{6}$ pointing to $O_{1}$ in Figure 2) may not be viable targets but upstream variables (e.g. $V_{5}$) may in fact be be promising targets and the direct parents of the outcome mediate the effect of these potential targets. In the presence of latent confounding, the identifiability of direct causal effects of outcomes is also jeopardized (Hernán and Robins 2006). For example, in Figure 2 misisng $U$ induces apparent correlation among $V_{3}$, $V_{4}$, $V_{5}$, and $V_{6}$ even though some of these variables are conditionally independent given $U$.

The extent of the problem of unmeasured confounding can be quantified in more detail. Suppose we model $O_{1}$ without controlling for $U$ (2):

Setting the coefficient of determination for the model

equal to $\rho_{3}^{2}$. Then the estimated variance of $\beta_{3}$ in the presence of collinearity can be related to the variance when collinearity is absent via the following formula (Rawlings, Pantula, and Dickey 1998):

Formula 8 describes the variance inflation factor (VIF) of $\beta_{3}$. Note that $\lim_{\rho\to 1}\frac{1}{1-\rho^{2}}=\infty$, so even mild collinearity induced by latent variables can severely distort coefficient values and signs, and thus estimation of ACE. Our approach reduces the VIFs of coefficients related to outcomes and thus make all causal statements relating to outcomes, such as calculation of ACE, more reliable, by controlling for $\bar{U}$ - the estimate of $U$ - in the network,

Consider an output $O_{j}$. While it is difficult to describe the limit of error on the coefficients of the drivers of $O_{j}$, it is straightforward to put a ceiling on the improvement in the likelihood obtainable from modeling $U$ and approximating $G_{u}$ with $G_{\bar{u}}$. Suppose we eventually model $U$ as a linear combination of a set of variables $X\subset V$, and denote by $X\setminus W$ the set difference: members of $X$ not in $W$. Then for any outcome $O_{i}$ predicted by a set of variables $W$ in the graph $G$ and in truth predicted by the set $Z+U$, we can contrast three expressions (from $G$ and $G_{\bar{U}}$ respectively):

Model (a) is the model that was actually accepted, subject to regularization, in G. Model (b) is the "complete" model of $O_{i}$ that controls for $U$ non-parsimoniously by controlling for all variables affected by $U$ and not originally in the model. The third model, (c), is the ideal parsimonious model when U is known. We can compare the quality of these models via the Bayesian Score, and the full score, which can approximated by the Bayesian Information Criterion (BIC) in large samples (Koller and Friedman 2009). We assume that model (c) would have the lowest BIC (being the best model), and model (a) would be slightly better, since we know that the set of variables $X\setminus W$ didn’t make it into the first equation subject to regularization by BIC. Assuming $n$ samples,

The "complete model" - model (b) - includes all of the true predictors of $O_{j}$. Therefore its score will be the same as that of the true model - model (c) - plus the BIC penalty, $log(n)$, for each extra term, minus the cost of having $U$ in the true model (that is, the cardinality of the relevant part of the latent space). We know that the extra information carried by this model was not big enough to improve upon model (a), that is $b_{a}<b_{c}+k\log(n)$ for some $k$. Rearranging:

Any improvement in $G_{\bar{U}}$ owing to modeling of $\bar{U}$ cannot, therefore, exceed $k\log(n)$ logs, where $k=|X\setminus W|-|U|$: the information contained in the "complete" model is smaller than its cost.

Although the available improvement in predictive power is also capped in some way, it is still important to aim for that limit. The reason is, correct inference of causality, especially in the presence of latent variables, is the only way to ensure transportability of models in real-world (heterogeneous-data) applications (see, e.g., (Bareinboim and Pearl 2016)).

This approach is a generalization of work presented in (Anandkumar et al. 2013), where the authors show that, under some assumptions. the latent space can be learned exactly, which is also related to the deconfounder approach described by (Wang and Blei 2019a). However, our approach does not require that the observables be conditionally independent given the latent space and instead generate such independence by the use of causal network’s residuals, which are conditionally independent of each other given the graph and the latent space. However, since the network among the observables is undefined in the beginning, the structure of the observable network must be learned at the same time as the structure of the latent space, which leads us to the iterative/variational bayes approach presented in Algorithm 1. Lastly, the use of the entirety of the residual space is different from the work described in (Elidan, Nachman, and Friedman 2007), where local residuals are pursued with the goal to accelerate structure learning while simultaneously discovering the latent space.

Algorithm 1 below describes our approach to learning the latent space and can be viewed as a type of an expectation-maximization algorithm.

How do we learn $\bar{U}=f(R_{\bar{U}})$? In the linear case, we can use PCA, as described above, and in the non-linear case, we can use non-linear PCA, autoencoders, or other methods, as alluded to above. However, the linear case provides a useful constraint on dimensionality, and this constraint can be derived quickly. A useful notion of the ceiling constraint on the linear latent space dimensionality can be found in (Gavish and Donoho 2014). From a practical standpoint, the dimensionality can be set even tighter, and we use a previously described heuristic approach(Buja and Eyuboglu 1992).

In this section, we discuss extensions to this approach for GGMs with interactions, nonlinear functional forms, and categorical variables.

In this work we present a method for describing the latent variable space which is locally optimal under linearity, up to rank-linearity. The method does not place a priori constraints on the number of latent variables, and will infer the upper bound on this dimensionality automatically.

In real-data situations, some assumptions of our linear approach may not hold - PCA and related methods may not suffice. For such cases, we provide a heuristic autoencoder implementation. Autoencoders have shown utility when the structure of the causal model is already known (Louizos et al. 2017). Learning structure over observed data as well as unconstrained latent space via autoencoders, as in our work, results in a hybrid "deep causal discovery" generalization of that approach.

Many domains can benefit from joint causal discovery and inference of latent space. For instance, in epidemiology neither the model structure nor the latent space are typically well-understood beyond simple examples. Multi-omic biological datasets are another area of applicability, since it’s impossible to collect all biological data modalities in practice. In clinical trials, epidemiological features of multi-omic datasets may be hard to fully account for.

Combining causal discovery and causal inference improves statements about ACE, and helps assess data set limitations. Further, we provide a relatively fast implementation suitable to real problems. Therefore, we hope that our work will open new practical applications of both paradigms.