Response Function Variables

This yields a functional causal model \(\{F_V:pa(V)\times U_V\to V\mid V\in\mathcal{V}\}\), although not in a computationally convenient form. However, by the assumption of finite-valued measured variables \(V\in\mathcal{V}\), we can in fact partition the ranges of the unmeasured ones \(U_V\in\mathcal{U}\) into finitely many sets; they simply enumerate the finitely many distinct functions \(pa(V)\to V\) for each \(V\in\mathcal{V}\). We call such a canonical partition \(R_V\) of \(U_V\) the response function variable of \(V\).

Original Example: Response Function Variables.

(The dashed line signifies \(R_X\not\!\perp\!\!\!\perp R_Y\).) In the context of our specific example, where e.g. \(X\) has a single parent only so \(R_X\) takes exactly \(4\) distinct values, say \(\{0,1,2,3\}\) where \(\forall z\in\{0,1\}\), we have

\[\begin{matrix} \textbf{Response Pattern}&\textbf{Function}&\textbf{Interpretation}\\ r=0&f_X(z,0)=0&\text{never takes $X$, regardless of $Z$}\\ r=1&f_X(z,1)=z&\text{full compliance with assignment $Z$}\\ r=2&f_X(z,2)=1-z&\text{total defiance of assignment $Z$}\\ r=3&f_X(z,3)=1&\text{always takes $X$, regardless of $Z$}\\ \end{matrix}\]

Parameters and their Relationships

We can observe the outcomes of the variables in \(\mathcal{V}\), so we can estimate their joint distribution, but in fact we will only need their conditional distribution given \(\mathcal{L}\), so let’s establish some notation for this. For our specific example we define \[p_{(y,x;z)}:=P(y,x|z),\quad\forall x,y,x\in\{0,1\}\] Although unknown, we need notation for the distribution of the response function variables as well. Our specific example requires only \[q_{(r_X,r_Y)}:=P(r_X,r_Y),\quad\forall r_X,r_Y\in\{0,1,2,3\}\] By the Markovian property, our example yields \[p_{(y,x;z)}=\sum_{r_X,r_Y}P(y|x,r_Y)P(x|z,r_X)q_{(r_X,r_Y)}\] Note that all coefficients are deterministic binary, since e.g. \(P(y|x,r_Y)=\begin{cases}1&\text{, if }y=f_Y(x,r_Y)\\0&\text{, otherwise}\end{cases}\), and by recursion the same applies to more complex examples. Thus, in our particular example, we have \[p=Rq\text{ for some matrix }R\in\{0,1\}^{8\times16}\] Further, we can express potential outcomes of \(Y\) under intervention on \(X\) in terms of the parameters \(q\) by the adjustment formula; \[P(y|do(x))=\sum_{r_X,r_Y}P(y|x,r_Y)q_{(r_X,r_Y)}\] Hence, we have \[\begin{align*} \text{ACE}_q(X\rightarrow Y) &=P(Y=1|do(X=1))-P(Y=1|do(X=0))\\ &=c^Tq\text{ for some vector }c\in\{0,1,-1\}^{16} \end{align*}\]

Linear Programming

Now we have our constraints on \(q\) as well as our effect of interest in terms of \(q\) and we are ready to optimize! By adding the probabilistic constraints on \(q\) we have arrived at e.g. the following LP giving a tight lower bound on \(\text{ACE}(X\rightarrow Y)\): \[\begin{matrix} \min&c^Tq\\ st&\Sigma q=1\\ \&&Rq=p\\ \&&q\geq0 \end{matrix}\] By the Strong Duality Theorem³ of convex optimization, the optimal value of this primal problem equals that of its dual. Furthermore, its constraint space is a convex polytope and by the Fundamental Theorem of Linear Programming⁴, this optimum is attained at one of its vertices. Consequently, we have \[\begin{align*} \min_q ACE_q(X\rightarrow Y) &=\min\{c^Tq\mid q\in\mathbb{R}^{16},q\geq0_{16\times1},1_{1\times16} q=1,Rq=p\}\\ &=\max\{\begin{pmatrix}1&p^T\end{pmatrix}y\mid y\in\mathbb{R}^{9},y\geq0_{9\times1},\begin{pmatrix}1_{16\times1}&R^T\end{pmatrix}\leq c\}\\ &=\max\{\begin{pmatrix}1&p^T\end{pmatrix}\bar{y}\mid\bar{y}\text{ is a vertex of }\{y\in\mathbb{R}^{9}\mid y\geq0_{9\times1},\begin{pmatrix}1_{16\times1}&R^T\end{pmatrix}\leq c\}\} \end{align*}\] Since we allow the user to provide additional linear inequality constraints (e.g. it may be quite reasonable to assume the proportion of “defiers” in the study population of our example to be quite low), the actual primal and dual LP’s look slightly more complicated, but this small example still captures the essentials.
In general, with \(A_e:=\begin{pmatrix}1_{n\times1}^T\\R\end{pmatrix}\), \(b_e:=\begin{pmatrix}1\\p\end{pmatrix}\) and user provided matrix inequality \(A_{\ell}q\leq b_{\ell}\), we have \[\begin{align*}\min_q(ACE(A\to Y))&=\min(\{c^Tq\mid q\in\mathbb{R}^n,A_{\ell}q\leq b_{\ell},A_eq=b_e,q\geq0_{n\times1}\})\\ &=\max(\{\begin{pmatrix}b_{\ell}^T&b_e^T\end{pmatrix}y\mid y\in\mathbb{R}^{m_{\ell}+m_e},\begin{pmatrix}A_{\ell}^T&A_e^T\end{pmatrix}y\leq c,\begin{pmatrix}I_{m_l\times m_l}&0_{m_l\times m_e}\end{pmatrix}y\leq0_{m_l\times1}\}\\ &=\max(\{\begin{pmatrix}b_{\ell}^T&b_e^T\end{pmatrix}\bar{y}\mid\bar{y}\text{ is a vertex of }\{y\in\mathbb{R}^{m_{\ell}+m_e}\mid\begin{pmatrix}A_{\ell}^T&A_e^T\\I_{m_l\times m_l}&0_{m_l\times m_e}\end{pmatrix}y\leq\begin{pmatrix}c\\0_{m_l\times1}\end{pmatrix}\}\})\end{align*}\] The coefficient matrix and right hand side vector of the dual polytope are constructed in the few lines of code below (where \(c_0:=c\)).

a1 <- rbind(cbind(t(A_l), t(A_e)),
            cbind(diag(x = 1, nrow = m_l, ncol = m_l), matrix(data = 0, nrow = m_l, ncol = m_e)))
b1 <- rbind(c0,
            matrix(data = 0, nrow = m_l, ncol = 1))
if (opt == "max") {
  a1 <- -a1
  b1 <- -b1
}

Vertex Enumeration

The last step of vertex enumeration has previously been the major computational bottleneck in causaloptim. It has now been replaced by cddlib⁵, which has an implementation of the Double Description Method (dd). Any convex polytope can be dually described as either an intersection of half-planes (which is the form we get our dual constraint space in) or as a minimal set of vertices of which it is the convex hull (which is the form we want it in) and the dd algorithm efficiently converts between these two descriptions. cddlib also uses exact rational arithmetic, so we don’t have to worry about any numerical instability issues, e.g. Also, it already interfaces with R⁶. The following lines of code extract the vertices of the dual polytope and store them as rows of a matrix.

library(rcdd)
hrep <- makeH(a1 = a1, b1 = b1)
vrep <- scdd(input = hrep)
matrix_of_vrep <- vrep$output
indices_of_vertices <- matrix_of_vrep[ , 1] == 0 & matrix_of_vrep[ , 2] == 1
vertices <- matrix_of_vrep[indices_of_vertices, -c(1, 2), drop = FALSE]

The rest is simply a matter of plugging them into the dual objective function, evaluating the expression and presenting the results!
The first part of this is done by the following line of code, where the utility function evaluate_objective pretty much does what you expect it to do (here (c1_num,p)\(=(\begin{pmatrix}b_{\ell}^T&1\end{pmatrix},p)\) separates the dual objective gradient into its numeric and symbolic parts).

expressions <- apply(vertices, 1, function(y) evaluate_objective(c1_num = c1_num, p = p, y = y))

There are still a few places in causaloptim that are worthwhile to optimize, but with this major one resolved we looking forward to researchers taking advantage of the generalized class of problems that it can handle!

Happy bounding!