Estimation of causal effects with small data in the presence of trapdoor variables

We consider the problem of estimating causal effects of interventions from observational data when well-known back-door and front-door adjustments are not applicable. We show that when an identifiable causal effect is subject to an implicit functional constraint that is not deducible from conditional independence relations, the estimator of the causal effect can exhibit bias in small samples. This bias is related to variables that we call trapdoor variables. We use simulated data to study different strategies to account for trapdoor variables and suggest how the related trapdoor bias might be minimized. The importance of trapdoor variables in causal effect estimation is illustrated with real data from the Life Course 1971-2002 study. Using this dataset, we estimate the causal effect of education on income in the Finnish context. Bayesian modelling allows us to take the parameter uncertainty into account and to present the estimated causal effects as posterior distributions.


Introduction
Understanding causal relations forms the basis of decision making in society. The role of statistics is to provide tools that allow us to estimate the causal effects of planned interventions. Instead of modelling just associations, a causal model describes the functional relationships present in the system of interest. A core feature of causal models (Pearl, 2009) is the capability to represent the effects of actions on the model via symbolic interventions. These interventions are assumed modular in the sense that they only modify the target of the intervention leaving other causal mechanisms in the model intact. The resulting probability distribution of this post-interventional causal model is defined as the causal effect of the intervention.
Various methods have been suggested for estimating causal effects in different settings. Propensity score matching (Rosenbaum and Rubin, 1983;Imbens, 2000), inverse probability weighting (Rosenbaum and Rubin, 1983;Rosenbaum, 1987), and g-methods (Robins et al., 1992) are some of the most well-known methods. These methods try to mimic a randomized experiment by creating pseudo-samples or by weighting the original sample in various ways, and it is typically assumed that there are no unobserved confounders present. In the classical structural equation modelling (SEM) approach (see, e.g., Kline, 2011), observed variables are assumed to be Gaussian, and unobserved confounders are treated as correlations between observed variables (there are also some extensions, such as LiNGAM for linear non-Gaussian cases (Shimizu, 2014)).
Another approach to causal inference is based on causal graphs, where the notion of identifiability plays a central role (see Section 2). For certain graphs, the well-known back-door and front-door adjustment criteria (Pearl, 1995) allow for identification of interventional distributions from observational data in the presence of unobserved confounders. The back-door criterion tells us whether a set of variables forms an admissible set, so that we only need to condition on these variables when estimating the causal effect, whereas the similar front-door criterion can be applied in the presence of a mediator between the interventional variable and the response variable. More generally, do-calculus (Pearl, 1995) can be used to assess whether the interventional distribution of interest is identifiable given only the known causal graph, without any parametric assumptions about the distributions of the variables or the form of the effects they have on each other. Do-calculus provides an identifying functional, i.e., a nonparametric formula for the interventional distribution consisting of terms that represent observational distributions. Although identifiability does not in general guarantee estimability (Maclaren and Nicholson, 2019), it is usually possible to obtain an estimator for the causal effect by replacing the terms present in the identifying functional by suitable parametric or nonparametric estimators and then combining the results accordingly.
We study the estimation of causal effects with small data in scenarios where standard adjustment criteria are not applicable. By small data, we refer to a case where parameter estimation exhibits non-negligible uncertainty due to the sample size. In addition, we consider the presence of functional equality constraints known as Verma-constraints (Verma and Pearl, 1990;Tian and Pearl, 2002;Robins, 1986). Under certain conditions (see Section 2), these constraints are related to special variables that we call trapdoor variables. These variables can bias the causal effect estimator for finite samples and we refer to this form of bias as trapdoor bias. We demonstrate the practical ramifications of trapdoor variables for the estimation of causal effects via simulations in a number of synthetic scenarios with small sample sizes and compare a variety of estimation strategies in both nonparametric and parametric settings.
As a motivating example, we consider Bayesian estimation of the causal effect of education on yearly income using real data from the Finnish Life Course 1971-2002 study (Kuusinen, 2018). We construct a causal model for this study where we take into account the grade point average (GPA) from primary school, language skills, gender and the socioeconomic status of the parents. In the causal model, we find that GPA is in fact a trapdoor variable due to a functional equality constraint on the causal effect of interest. Bayesian modelling allows us to estimate the full post-interventional distribution of the income on different levels of education, which indicate a clear positive causal effect of education on income. We combine Bayesian estimation with a specialized Monte Carlo approach in order to take the effect of the trapdoor variable into account in a number of scenarios including the Life Course model. All analysis was done in the R environment (R Core Team, 2020), and the codes for the simulation experiments, the Life Course example, and the figures for the simulation results (created with the ggplot2 package (Wickham, 2016)) of this paper are available at https://github.com/helske/trapdoor. The paper is structured as follows. Section 2 introduces the notation, gives a definition for the trapdoor variables and present examples on causal models where such variables are present. Section 3 focuses on various aspects related to the estimation of causal effects in the presence of trapdoor variables including a Bayesian approach and demonstrates how a trapdoor variables manifests under a linear-Gaussian model. Section 4 considers the effect of trapdoor variables and the trapdoor bias of causal effect estimators via simulation in a model with binary variables, a linear-Gaussian model and a synthetic scenario based on the Life Course model. The analysis using the real Life Course data is presented in Section 5. Section 6 provides some concluding remarks.

Notation and basic definitions
Our analysis is based on the framework of structural causal models (SCM) and directed graphs, and we assume the reader to be familiar with these concepts and their core probabilistic and graphical properties. For a more detailed discussion on SCMs and graph theoretic concepts, we refer the reader to works such as (Pearl, 2009) and (Koller and Friedman, 2009).
We use capital letters to denote variables (V ) and small letters to denote their values (v). Bold letters are used to denote sets of variables (V) and value assignments (v). The set of all possible value assignments to V is denoted by val(V). Set difference of sets A and B is denoted by A \ B. We use shorthand notation P (Y | x) to denote the probability distribution P θ (Y | X = x) where we typically omit the dependence of (unknown) model parameters θ.
Each SCM M over a set of variables V is associated a joint probability distribution P (V) in a population of interest and a causal graph G over V where directed edges between two observed variables in V correspond to direct causal relationships which are assumed to not form any cycles. Bidirected edges between two observed variables in V are used to denote confounding by an unobserved common cause. In this framework, interventions are represented using the do(·)-operator. An intervention do(X = x) forces the variables in X to take the values specified by x while leaving other mechanisms of the model intact. This intervention induces a submodel M x with the interventional distribution P (V | do(X = x)). A causal effect P (Y | do(X = x)) is said to be identifiable in G if it is uniquely computable from P (V) in any SCM that induces G. A variable Y in the post-intervention model M x is denoted as Y (x). For an identifiable causal effect, an identifying functional is a function f such that f (P (V)) = P (Y | do(X = x)). We assume that P (V = v) > 0 for values v ∈ val(V) making all conditional distributions and interventions well-defined.
Note that identifiability only indicates the existence of an estimator and does not take into account the potential problems stemming from finite data. Therefore, even though determining the identifiability of P (Y | do(X = x)) is an important first step, it does not guarantee that we can estimate the causal effect in practice without additional stability assumptions (Maclaren and Nicholson, 2019). Note also that through this paper we assume that our causal graph is correct, while in practice we are rarely certain of this.
Interventional distributions exhibit conditional independence constraints, which can be characterized by d-separation (Pearl, 1988) in the associated causal graph of the model. However, identifiable causal effects can be subject to Verma-constraints. As an example of such a constraint, we show that in the causal graph of Fig. 1, a causal effect does not depend on the value of a variable W despite it appearing in the identifying functional of the interventional distribution. The causal effect of X on Y is identifiable in this graph which can be verified using do-calculus or by applying an identifiability algorithm such as the one by Huang and Valtorta (2006) or the ID algorithm by Shpitser and Pearl (2006). Application of the ID algorithm implemented in the R package causaleffect (Tikka and Karvanen, 2017) provides the formula where the left-hand side depends on the value of X and Y , but on the right-hand the variable W is also present and it is not subject to summation, unlike the variable Z. However, the right-hand side cannot depend on the value of W , as there is an admissible set Z which gives us an alternative formula i.e., a simple back-door formula where variable W is not present. Y X Z W Figure 1: A causal graph where the identifying functional of P (Y | do(X = x)) obtained by an application of the ID algorithm does not depend on the value of of W and there is an admissible set Z enabling back-door adjustment.
Perhaps surprisingly, Verma-constraints can be expressed as conditional independences in the interventional distribution (Shpitser and Pearl, 2008) and can be used to give an alternative definition for nested Markov models in acyclic directed mixed graphs (Richardson et al., 2017). Verma-constraints have been used for testing edges (Shpitser et al., 2009) and for marginalization via variable elimination (Shpitser et al., 2011).

Trapdoor variables
Here we give a broad definition that captures the notion that a set of variables Z may appear in an identifying functional of a causal effect, but the value of the causal effect is not dependent on the value of Z.
Definition 1 (Functional independence) Let X, Y, Z ⊂ V be disjoint sets and let P (Y | do(X = x)) be an identifiable causal effect with an identifying functional g(v) = f (P (v)). If the domain of g is val(X) × val(Y) × val(Z), and g(x, y, z 1 ) = g(x, y, z 2 ) for all x ∈ val(X), y ∈ val(Y) and z 1 , z 2 ∈ val(Z), then g is functionally independent from Z.
The constraint defined above is specific to the given identifying functional. In some instances we may be able to find identifying functionals that do not exhibit functional equality constraints for any subset Z of V. Our interest lies in the opposite direction, where every identifying functional exhibits a specific type of functional independence. Before characterizing this property of interest, we must first define an operation known as the latent projection of a causal graph (Pearl and Verma, 1991).
Definition 2 (Latent projection) Let G be a causal graph over a set of vertices V ∪L. The latent projection L(G, V) is a causal graph over V where for every pair of distinct vertices Z, W ∈ V it holds that (a) L(G, V) contains an edge Z −→ W if there exists a directed path Z −→ · · · −→ W in G on which every vertex except Z and W is in L.
(b) L(G, V) contains an edge Z ←→ W if there exists a path from Z to W in G that does not contain the pattern Z −→ M ←− W (a collider), on which every vertex except Z and W is in L, the first edge of the path has as arrowhead into W , and the last edge has an arrowhead into Z.
Latent projections can be used to derive identifying functionals for causal effects such that they do not contain a specific variable, i.e., the variable is considered latent, and the causal effect of interest is identified in the corresponding latent projection (Tikka and Karvanen, 2018). For this reason the presence of a functional constraint in some identifying functional does not rule out the possibility of obtaining another identifying functional that is not subject to the same constraint (Recall the example on the identifying functional of P (Y | do(X = x)) in the graph of Fig. 1). We rule out this possibility of finding alternative identifying functionals by restricting our attention to settings where a trapdoor variable is present.
Finding trapdoor variables is straightforward as outlined by their definition. Given a causal effect of interest, we first determine its identifiability and whether the identifying functional is subject to functional equality constraints. If this is the case, we proceed to verify whether the causal effect might be identifiable when the set Z is considered latent using a suitable latent projection. The algorithm by Tian and Pearl (2002) can be used to systematically enumerate functional equality constraints implied by a causal model. Evans (2018) showed that this constraint-finding algorithm is complete for categorical variables, but it is not known whether the algorithm is complete in general, i.e., whether all such constraints can be found by the algorithm.
If there exists a trapdoor variable for an identifying functional of a causal effect of interest, then the estimate of the causal effect may depend on the value of the trapdoor variable. This dependency can introduce bias in the estimate.
Definition 4 (Trapdoor bias) Letĝ(v) be an estimator of an identifying functional g(v) of a causal effect P (Y | do(X = x)) and let Z be a set of trapdoor variables with respect to g(v). Let B(ĝ(v)) denote the bias of this estimator. If there exists z 1 , z 2 ∈ val(Z) so that B(ĝ(x, y, z 1 )) = B(ĝ(x, y, z 2 )), thenĝ(v) exhibits trapdoor bias with respect to Z.
For a consistent estimator, the effect of the trapdoor bias becomes negligible as the sample size grows, but for small samples the choice of how to take the trapdoor variables into account may be significant. Figure 2: Three causal graphs of increasing complexity, where the interest is in the causal effect of X on Y .

Example on trapdoor variables
Consider the three causal graphs in Fig. 2. We are interested in estimating the causal effect of X on Y . For example, in our application in Section 5, X will be the education level and Y is the yearly income. Furthermore, Z corresponds to the GPA from primary school, and W to the socioeconomic status of the parents. In all three graphs we have an arrow from X to Y , meaning that we assume that there is a direct causal effect of education on income. In addition, we assume that there may be some unobserved confounders between socioeconomic status of the parents and income. In the graphs of Fig. 2b and Fig. 2c, we also assume that there is confounding between socioeconomic status of the parents and education level of the participant. In Fig. 2c, we assume that the effect of socioeconomic status on the education level is mediated by the GPA. We will further extended the third graph with additional variables, such as gender, in Section 5. Now consider the estimation of the interventional distribution P (Y | do(X = x)), i.e., the distribution of Y when we intervene on X by setting it to x, which differs from a simple conditional distribution of P (Y | X = x) in these graphs. In order to estimate P (Y | do(X = x)), we need to find a formula for it in terms of the observed variables only, i.e., Y , X, W , and Z in our example. In the graph of Fig. 2a we obtain the so called back-door adjustment formula By conditioning on W we block all back-door paths from X to Y i.e., those paths between X and Y which have arrows into X. Thus by estimating the parameters of the terms P (Y | w, x) and P (W ) we can estimate the interventional distribution of interest, P (Y | do(X = x)). For example, assuming that all variables are Gaussian and their relationships are linear, we have , with b ij denoting the estimated regression coefficient of variable j on variable i, a i denoting the intercept term and s 2 i corresponding to the estimated residual variance. Now consider the graph of Fig. 2b. In this case, conditioning on W blocks the path Y ←→ W −→ X but opens the path Y ←→ W ←→ X meaning that we cannot apply the back-door adjustment here. In fact, adding the unobserved confounder between W and X renders the causal effect nonidentifiable.
In the graph of Fig. 2c, we assume that we have obtained data on variable Z which lies on the directed path from W to X. Given this additional information, the causal effect P (Y | do(X = x)) is again identifiable, and we have . (1) However, there is no term for the distribution of Z in equation (1). Thus after estimating the distributions P (Y | x, z, w), P (X | z, w), and P (W ), Z is essentially reduced to a fixed but unknown parameter in the context of P (Y | do(X = x)). This graph is also considered by Tian and Pearl (2002), Pearl and Mackenzie (2018), and Jung et al. (2020). Tian and Pearl (2002) show that this graph contains a functional equality constraint of an interventional distribution which cannot be expressed as a conditional independence constraint using the observed variables. The constraint states that the formula for P (Y | do(X = x)) given in equation (1) is functionally independent of Z, meaning that its value is independent on the choice of z, just as we would intuitively expect. However, when estimating equation (1) from the data we clearly must choose some value for Z, even though the constraint states that the actual value should not matter. In fact, Z is a trapdoor variable with respect to the identifying functional of equation (1) in this graph. This follows from the functional equality constraint and the fact that the causal effect is not identifiable in the causal graph of Fig. 2b which is the latent projection of Fig. 2c when Z is considered latent. Trapdoor variables are in no way limited to the scenario of Fig. 2c. Fig. 3a presents a generalization of Fig. 2c with an additional set of variables B = {B 1 , . . . , B k }, that are possibly connected to each other with directed arrows or through unobserved confounders. In addition, any member of B is possibly connected to W via a bidirected edge or to Z, X or Y via a directed edge. Fig. 4 depicts three additional examples of graphs where trapdoor variables are present. In the following sections we investigate the challenges that trapdoor variables impose on the estimation of causal effects. Figure 3: A generalization of the setting presented in Fig. 2c with an additional set of variables B is shown in (a). The square node for B denotes an arbitrary causal graph over B. Edges to and from B mean that such an edge may exist between any member of B and the other endpoint of the edge. An instance of (a) is shown in (

Plug-in estimation
A straightforward strategy for causal effect estimation given the assumed causal graph is to construct a parametric model for the conditional distributions that appear in the formula of an identifying functional of a causal effect. As an example, consider the graph of Fig. 2c, and assume that we wish to estimate E(Y | do(X = x)). We need to estimate the unknown model parameters θ corresponding to the causal effect of equation (1) and we let Pθ(·) denote the density where the unknown parameters θ are assigned some fixed values such as their maximum likelihood (ML) estimates. The plug-in estimator for the expected value of the interventional distribution takes the following form . (2) Note that this estimate may depend on the value of z when using a small sample estimate for θ. This dependency and the corresponding trapdoor bias vanishes when the models are correctly specified and the true parameter values are used (by the definition of functional independence). However, even if we assume an unbiased estimatorθ, the estimator (2) can still exhibit bias since it is a nonlinear function ofθ that depends on the choice of z.
For large enough datasets, the choice of z in equation (2) should not have a large impact, but we will show that for small samples, the value of the trapdoor variable z can have a crucial role in estimating the correct causal effect of X on Y .
An alternative strategy for defining the plug-in estimator is to eliminate the functional constraint by averaging the identifying functional over a (conditional) distribution of the trapdoor variable. If Z is a trapdoor variable and T ⊆ (V \ Z), then since f (P (v)) is functionally independent of Z. This leads to the following alternative estimator for the expected value of the interventional distribution in the causal graph of Fig. 2c Eθ . ( Note that while this strategy eliminates the problem of having to choose a specific value for z, it also introduces a new problem of selecting the set T and estimating the model parameters for the distribution P θ (Z | t). As before, this estimator is a nonlinear function of the parameters that can result in bias. Note that in a frequentist framework, there is no way to divide the total bias of the final estimate into trapdoor bias, plug-in bias, or other possible sources of bias.

Bayesian estimation of causal effects
We advocate the use of Bayesian methods for jointly estimating the parameters of the distributions of identifying functionals of interest. By drawing samples from the joint posterior of the model parameters (using any generic Bayesian inference engine, typically some type of Markov chain Monte Carlo (MCMC) algorithm), we propagate the parameter estimation uncertainty to the final causal effect estimates and avoid the plug-in bias due to the nonlinear formula of the identifying functional with respect to model parameters. This allows us to focus on the effects of the trapdoor variable and the trapdoor bias of the estimators. We can also obtain samples from the full posterior of P (Y | do(X = x)) which can be used for straightforward evaluation of any properties of interest (such as mean and variance) of this posterior.
In parametric estimation of causal effects, we typically do not have an analytical formula for the types of plug-in estimators that are shown in equations (2) and (3). In these cases, we can use a Monte Carlo approach to draw samples from Pθ(Y | do(X = x)), and estimate the desired quantity using these samples. The specific Monte Carlo algorithm depends on the causal graph, the identifying functional and the corresponding conditional distributions. In simple cases, we simulate variables from their (conditional) distributions with fixed x, whereas for example in the case of equation (1) where P (X | ·) is present, we need additional weighting of the samples. Consider our extended example graph in Fig. 3a, which leads to the following formula for P (Y | do(X = x)): Algorithm 1 describes the Monte Carlo algorithm for equation (4) based on the second approach discussed in Section 3.1, given the parameter vectorθ (the same algorithm is Algorithm 1 Monte Carlo algorithm for sampling from Pθ(Y | do(X = x)) defined by equation (4) 3: Set or sample the value of the trapdoor variable, for example as For j = 1, . . . , M :

5:
Sample w ij ∼ Pθ(W ) 6: Compute the normalized weights: also suitable for (1) by omission of the references to the variables in B). Algorithm 1 draws samples from the marginal and conditional distributions defined in (4) and gives us N × M weighted replications from Pθ(Y | do(X = x)) which allows us to compute, e.g., where the weightsγ ij are defined in the last step of Algorithm 1.
Here the trapdoor variable has a more crucial role than in the case of analytical formulas. With poor choices of z, our weightsγ ij can become degenerate, i.e., most of the weights are near zero. In order to avoid this, it is natural to condition the trapdoor variable on x and perhaps other variables such as members of b i in Algorithm 1 which should make the weights well behaved.
The suitable number of Monte Carlo samples N × M can determined by computing the Monte Carlo standard error (MCSE) of our causal effect estimate. The MCSE measures the additional uncertainty in our estimate due to the finite Monte Carlo sample size. For example, the MCSE for estimator (5) can be computed as When combining Algorithm 1 with Bayesian parameter estimation, the algorithm is used at each MCMC iteration given the current values of the model parameters θ.
With the MCMC approach, we can compute functions of Y (x) at each iteration, giving us samples from its posterior distribution, or we can store all N × M weighted Monte Carlo samples leading to posterior distribution of Y (x), i.e., the variable Y in the post-interventional distribution, which can be further used to evaluate functions of interest. In the latter case, it can be practical to resample (using the corresponding weights) and store only one replication y ij at each iteration if memory constraints limit the storing of all samples.

Analytical causal effect for a linear-Gaussian model
We will now consider the estimation of causal effects in a common linear-Gaussian case for the causal graph of Fig. 2c. We assume that the underlying model is defined as Here all the parameters µ · , α · , β ·· , and σ · are unknown, and U and V are unobserved confounders. Our observational model needed for estimating P (Y | do(X = x)) is where parameters a · , b ·· , and s · are unknown and have to be estimated from the data. Now using equations of model (7) to equation (1) yields and While the variance V ar(Y | do(X = x)) does not contain the trapdoor variable Z, the expected value E(Y | do(X = x)) does, which is not surprising since P (Y | do(X = x)) = P (Y | do(X = x, Z = z)) in the graph of Fig. 2c. Also, if the interest is only in the difference E(Y | do(X = x + 1)) − E(Y | do(X = x)) then the effect of trapdoor variable cancels out in this linear case. As our model is linear-Gaussian, we can apply path analysis (Wright, 1934) to our causal graph (see, e.g., Pearl (2013) for examples) to find the marginal covariance matrix Σ for (Y, X, Z, W ) (shown in Appendix). From Σ, using the properties of multivariate normal distribution, we can obtain the theoretical formulas for the parameters (a · , b ·· , s · ) of model (7) in terms of the true parameters (µ · , α · , β ·· , σ · ) of the causal graph (see Appendix for details). Plugging these into equation (8), we obtain E(Y | do(X = x)) = α y + β yu µ u + β yx x, and V ar(Y | do(X = x)) = β 2 yu σ 2 u − 2β wu β zw β xz β yx β yu σ 2 u + σ 2 y . As expected, given the true causal model and its known parameters, the effect of the trapdoor variable cancels out. Nevertheless, with a finite dataset, the estimate of expected value (8) depends on z, unless b yz − bywbxws 2 w b 2 xw s 2 w +s 2 x b xz happens to estimate to zero.

Binary model
We will now illustrate different choices for the trapdoor variable with nonparametric estimation of P (Y | do(X = x)) in a case where all variables are binary. Consider Bernoulli variables defined in accordance with Fig. 2c as where V and U correspond to the bidirected edges between X and W and between W and Y , respectively. By solving equation (1) analytically we obtain P (Y = y | do(X = x)) = (0.2 + 0.4x) y (0.8 − 0.4x) 1−y , y ∈ {0, 1}, which does not depend on z. However, in practice when U and V are unobserved and we need to estimate distributions P (W ), P (X | z, w), and P (Y | x, z, w) from finite data, we will show how the estimate of P (Y | do(X = x)) can depend on the chosen value z.
As an example, we simulated 100 000 datasets according to model (9) with sample sizes 100, 300, and 500, and estimated P (Y | do(X = x)) nonparametrically with various methods for dealing with the trapdoor variable. Besides fixing z to zero or one and using the estimator (2), we also computed a weighted average of these estimates, where the weights were based on the estimated marginal distribution P (Z), or the conditional distribution P (Z | x), corresponding to the estimator (3). With sample sizes 100 and 300, there were some cases where fixing z to zero or one lead to undefined probabilities (for example, there were no observations for which x = 0, z = 0 and w = 0). In these cases (about 10% of cases with sample size 100, and less than 0.1% with sample size 300), the entire replication was discarded. Fig. 5 shows the results of the simulation. We see that using the fixed values with z = 0 or z = 1 leads to some bias. The weighted average estimators perform better, and the one based on P (Z | x) outperforms the estimator based on P (Z). Overall, we are not far off from the ground truth in this simple setting.

Linear-Gaussian model
Now consider a model based on equation (6) with We will compare various approaches to account for the trapdoor variable Z. Instead of the nonparametric approach used in Section 4.1, we now switch to parametric Bayesian modelling. For comparative purposes, we use uniform priors for all unknown model parameters (a · , b ·· , and s · in equation (7)). As stated in Section 3, the Bayesian approach takes into account the uncertainty in P (Y | do(X = x)) due to parameter estimation by integrating over the posterior distribution of the parameters and avoids the plug-in bias due to the nonlinearity of equation (8).
For model estimation, we wrote a Stan model (Carpenter et al., 2017;Stan Development Team, 2019) which simultaneously estimates all unknown model parameters and E(Y | do(X = x)) using MCMC. We estimate the true expected causal effect E(Y | do(X = x)) using the causal graph with observed confounders, with x ∈ {0, 3, 6, 9}, and compare it to the estimates obtained by six different methods: (a) Fix the trapdoor variable Z to the constant z = 0 in equation (8).
(b) Marginalize over P (Z) i.e., in addition of estimating model (7) also estimate P (Z), so that the estimator (8) uses z = E(Z) at each MCMC iteration.
(c) As above, but estimate P (Z | x) and use z = E(Z | x).
so that the contribution of z is fixed to zero.
(e) Fit a structural equation model (SEM), a common approach for linear-Gaussian causal modelling (Kline, 2011). We used the R package lavaan (Rosseel, 2012) for this purpose.
(f) Use the composition of weighting operators (CWO) (Jung et al., 2020), which applies weighted regression to estimate functions of interventional distributions for arbitrary identifying functionals.
Based on model (10) we have E(X) = 6, E(Z) = 4, and E(Z | X = 0) = 0.25, E(Z | X = 3) = 2.1255, E(Z | X = 6) = 4 and E(Z | X = 9) = 5.875. We sampled 1 000 replications of varying sample size from model (10) and for each replication estimated E(Y | do(X = x)) using a single MCMC chain with 10 000 postwarmup iterations and averaging over the iterations (thus taking account the uncertainty from parameter estimation), except for the SEM and CWO approaches which are based on the maximum likelihood estimates. Fig. 6 shows how, despite theoretical equivalence, results depend heavily on the chosen strategy. Perhaps the most natural choices of z = 0 and z = E(Z) are prone to bias which depends on the value x, but the case where z is adapted based on x performs well. The SEM approach, which explicitly models the covariance structures between variables, also performs well due to the structure of the graph (see Shpitser et al. (2018) for more details), but is not applicable in more general graphs with non-Gaussian or nonlinear equations. The CWO method shows considerable bias when x = E(X).

Analytical vs. Monte Carlo approach for the linear-Gaussian model
Consider again model (10) but now with i.e., smaller measurement error in X, leading also to a narrower density of P (X | z, w). Because of this, a poor choice for the value of the trapdoor variable can cause most of the normalized weightsγ in Algorithm 1 to be close to zero, leading to inefficient Monte Carlo sampling. As an example, we simulated one replication of size n = 100 from this model, and estimated the posterior distribution of E(Y | do(X = 0)), both using the analytical formula and Monte Carlo with N = 500. For MCMC, we ran 4 chains with a total of 100 000 post-warmup iterations. In addition, at each MCMC iteration we sampled one replication from the posterior predictive density P (Y | do(X = 0)) from the weighted Monte Carlo sample. Fig. 7 shows the posterior distributions. We can see that with suitable z = E(Z | X = 0) both analytical and Monte Carlo methods give approximately equal results. A poor choice of z = E(Z) biases results compared to the analytical solution (which is also biased) as can be seen from the discrepancy between

Non-Gaussian model
In Section 5 we study the causal effect of education level on income using real data.
Here we perform a simulation experiment where we emulate the real-data case using the assumed graph of the real-data case and define the data generating process so that it reflects the nature of the variables in Section 5. We use graph of Fig. 3b where we obtain the following formula for P (Y | do(X = x)): s,g P (g)P (s) w P (Y | x, z, w, s, g)P (x | z, w, s, g)P (s | w)P (w) w P (x | z, w, s, g)P (s | w)P (w) .
Compared to the linear-Gaussian experiment of Section 4.2, here additional difficulties arise due to the fact that in addition to the unknown parameters θ, our distributional assumptions for the terms in equation (11) are not correct. For example, while (Y | X = x, U 1 = u 1 , S = s, G = g) is Gamma distributed by definition, (Y | X = x, S = s, G = g, W = w, Z = z) might not be. This can naturally bias our causal estimates, but the question remains whether different choices for the trapdoor variable affect the bias or precision of the causal estimates.
For the terms in (11), we assume a Gamma distribution for Y , and model its expected value given other variables via log-link, assuming a monotonic effect (Bürkner and Charpentier, 2020) of X and W . We treat X and W as ordinal variables and model them with ordered logistic regression and a sequential model with a logit-link (Tutz, 1990;Bürkner and Vuorre, 2019) respectively. We use a normal distribution for S, and a Bernoulli distribution for G.
We simulated 1 000 replications of sample size 500 from this model, and estimated the posterior mean of P (Y | do(X = x)) using the four different strategies for the trapdoor variable: We use a Beta distribution for Z in all cases by modelling the expected value via a logit-link and precision via a log-link. Based on the data generating process, the true causal effects are 19 690, 24 049, and 32 463 for x = 0, 1, 2 respectively. The root mean square error (RMSE) and bias of our causal estimates compared to ground truth are shown in Table 1. We see that RMSE increases with respect to the intervention variable X and the differences between trapdoor strategies are relatively small except when x = 2, where conditioning on the intervention variable results in approximately 20% smaller RMSE than when using the marginal distribution of the trapdoor variable or when conditioning only with covariates S and G. The differences in bias estimates are within the Monte Carlo error in the case of X = 0. With X = 1, the trapdoor strategies without conditioning on X results in positive bias whereas conditioning on X leads to a bias of similar magnitude but in a different direction. For X = 2, the bias is substantially larger with trapdoor strategies not involving X. In this simulation experiment, conditioning on the covariates B does not improve the estimates. Table 1: RMSEs and biases of the estimates for the causal effect of X on Y in the non-Gaussian simulation experiment using four different strategies for sampling the trapdoor variable Z, and the true data generating process with observed confounders. MCSEs (computed with simhelpers R package (Joshi and Pustejovsky, 2020)) were between 0.4 and 2 for the RMSE estimates, and between 14 and 62 for the bias estimates, increasing with respect to X in both cases.

Causal effect of education on income
As an example with real data, we analyse Life Course 1971-2002 dataset from Finnish Social Science Data Archive (Kuusinen, 2018). This longitudinal study consists of life courses of 634 Finnish children born in 1964-1968 in Jyväskylä, Finland. In the early 1970s, when the children were aged between three and nine, the Illinois Test of Psycholinguistic Abilities (ITPA) was used to test their verbal intelligence. Participants were then followed up, with further information on their life events gathered in 1984, 1991, and 2002. While the dataset has been used in various studies regarding the intelligence and school achievement (see, e.g., Kuusinen and Leskinen, 1988), we use this data to study the causal effect of education level (secondary or less, lowest/lower tertiary, or higher tertiary) on yearly income (euros, in 2000). Note that due to the limited geographical scope of the data, the participants do not necessarily form a fully representative sample of the corresponding birth cohorts in Finland. In addition to earned income Y and education level X, we have information on the GPA from primary school Z, socioeconomic status (SES) of the parents W (low, middle, high), gender G, and the ITPA score S (the General Language composite variable, a combination of the subtests). Out of the 634 participants in the full dataset, we use 509 participants with fully observed aforementioned variables and avoid considerations related to missing data. The assumed causal graph is shown in Fig. 3b, and as before, all analysis are based on the assumption that this graph is correct. We use same distributional assumptions as in our simulation experiment in Section 4.4, after transforming the original scale of GPA from 4-10 to 0-1).
It could be argued that there should be direct arrows from parental SES to child's education level and income. However, in Finland it is likely that the effect of SES to education is strongly mediated by the child's school achievement (variable Z in our graph) (Acacio-Claro et al., 2018). Also, the socioeconomic status in our data was coded based on the occupations of the parents, which in turn depends on their education. We thus assume that the unobserved confounder between W and X contains, among other things, the education levels of the parents. Similarly, there are studies suggesting that in Finland, after accounting for a child's education, the effect of family income on children's income is low (Österbacka, 2001). Finally, the intergenerational mobility of education and income are among the highest in Finland (and Nordic countries in general) compared to other countries (Pfeffer, 2008;Björklund et al., 2002). This suggests that even if there is a direct arrow from W to X or Y , these direct effects are likely negligible. Nevertheless, it is, of course, possible (and perhaps likely) that our causal graph is an oversimplification of the complex causal mechanisms related to education and income. We used four chains with total of 100 000 post-warmup iterations and N = M = 250, leading to MCSE around 50-80, depending on x with both trapdoor strategies. Table 2 shows the posterior summary statistics for median(Y | do(X = x)) and E(Y | do(X = x)) for all education groups. We see a clear discrepancy between the estimates based on the conditional trapdoor variable Z ∼ P (Z | x, s, g) and the marginal trapdoor variable Z ∼ P (Z). Finally, Fig. 8 shows the full posterior distributions P (Y | do(X = x)). The two strategies for the trapdoor variable give somewhat different results: When conditioning the trapdoor variable grade on the education level the interventional distributions between educational levels differ more than with sampling the trapdoor variable from P (Z). For example, we estimate the mean annual income for the highest education level as 29 500 euros using the conditional trapdoor variable compared to 26 600 euros when using the marginal trapdoor variable. On the basis of the simulation results of Section 4.4, we prefer the estimates with Z ∼ P (Z | x, s, g). Unsurprisingly, obtaining a higher education level has a positive causal effect on income. Our causal estimates also differ somewhat from the observed incomes in different education groups, with the observed median incomes being (from lowest to highest) 17 700, 19 300, and 26 800 euros, and the mean incomes 18 700, 21 800, and 28 700 euros.

Conclusion
We have shown how it is possible to estimate causal effects when the back-door and front-door adjustments are not applicable and highlighted the potential issues related to the application of theoretical identifying formulas to finite datasets. Our real-data example on the effect of education on income illustrates how Bayesian causal inference can be applied to complex causal graphs with multiple types of variables, and how trapdoor variables can have a substantial effect on the estimated causal effects. The basic structure that leads to implicit functional constraints and trapdoor variables was presented in the graph of Fig. 2c. This graph has been considered in the literature earlier but according to our knowledge, it has not been fully analysed from the viewpoint of estimation. This basic structure can be extended in several directions so that the causal effect is still identifiable while retaining the implicit functional constraint. The graph for the Life Course 1971-2002 study (Fig. 3b) is an example of such an extension. There are other interesting graphs such as those shown in Fig. 4 with similar or even more complex identifying functionals which exhibit the problems related to trapdoor variables.
Determining the optimal method to account for the presence of a trapdoor variable remains a challenging problem. The trapdoor bias exhibited for small sample sizes depends on the (possibly parametric) assumptions about the causal model as well as the properties of the estimators used for the conditional distributions appearing in the identifying functional. Even in the simple case of the linear model with a single trapdoor variable, the optimal method is not apparent. It may also be the case that minimizing the trapdoor bias might result in a large variance of the causal effect estimator necessitating further considerations about the estimation problem at hand.
Our simulation experiments in Bernoulli and linear-Gaussian cases illustrated how choosing the value of the trapdoor variable can have substantial effect in estimating the interventional distribution P (Y | do(X = x)). Our results suggest that a good default for the trapdoor variable should be based on its conditional distribution given the interventional variable X. On the other hand, for nonlinear and non-Gaussian models the effect of the trapdoor variable can have nonlinear effects on P (Y | do(X = x)) which can further bias the results. Therefore as a general sensitivity check, we recommend computing the causal effect using various strategies to account for the trapdoor variables and reporting how sensitive the results are with respect to these strategies.

Acknowledgements
This work belongs to the thematic research area "Decision analytics utilizing causal models and multiobjective optimization" (DEMO) supported by Academy of Finland (grant number 311877). We thank Yonghan Jung for providing example codes for the CWO method, and Satu Helske for providing insights on the intergenerational mobility.