UW STAT 535: Statistical Machine Learning (2026 Autumn)
Lecture 14: Causal inference

Yen-Chi Chen (University of Washington)

Potential outcome model

The causal inference is a very popular topic in modern statistics and machine learning research. To draw causal conclusion in a statistical models, we need to setup a framework for making causal statements. There are two popular frameworks: the potential outcome model and the graphical model. We first introduce the potential outcome model.

Let \(Y\in\mathbb{R}\) be the response variable/variable of interest and \(A\in\{0,1\}\) be the binary treatment. \(A=1\) refers to the case where the individual receives a treatment (treatment group) and \(A=0\) refers to the case where the individual receives a placebo (control group). You can think of \(Y\) as a measure of health condition (such as blood pressure) and the binary treatment \(A\) refers to whether this individual receives a certain treatment or not. The goal is to study the causal effect of \(A\) on \(Y\).

Under this scenario, our data consists of pairs \[(Y_1,A_1),\cdots, (Y_n,A_n),\] where \(Y_i\) is the outcome of the \(i\)-th individual and \(A_i\) is the treatment indicator of the \(i\)-th individual.

If the treatment \(A\) indeed has a causal effect on \(Y\), then we should think of two versions of \(Y\), denoted as \(Y(0)\) and \(Y(1)\). \(Y(0)\) is the outcome variable if the individual does not receive any treatment (\(A=0\)). \(Y(1)\) is the outcome variable in the case where the individual receives a treatment (\(A=1\)). The above model is called the potential outcome model. Here is a key concept of the potential outcome model: \[``Y|A=0" = ``Y(0)|A=0",\quad ``Y|A=1" = ``Y(1)|A=1".\] Namely, given \(A=0\), we can replace \(Y\) by \(Y(0)\) and given \(A=1\), we can replace \(Y\) by \(Y(1)\).

In the potential outcome model, every individual has two outcome variables. One is the observed outcome that we can observe and the other is the potential outcome that we do not observe. For instance, suppose that \(A_i=0\) (no treatment), then \(Y_i = Y_i(0)\) is the observed outcome. The other outcome \(Y_i(1)\) is the potential outcome of the response \(Y_i\) (if the individual receives the treatment). Thus, we only observe one of the two outcomes. Note that the observed response \(Y_i = Y_i(A_i)\).

The causal effect can be viewed as the distributional difference between random variables \(Y(1)\) and \(Y(0)\). A simple summary of the difference is the mean difference, which is also known as the average treatment effect (ATE): \[\tau = \mathbb{E}(Y(1)) - \mathbb{E}(Y(0)).\] So we will think of methods for estimating the ATE.

One may think of using the difference in conditional mean to estimate the ATE. Namely, we use the estimator \[\widehat\tau_{\sf naive} = \frac{\sum_{i=1}^n Y_i I(A_i=1)}{\sum_{i=1}^n I(A_i=1)} - \frac{\sum_{i=1}^n Y_i I(A_i=0)}{\sum_{i=1}^n I(A_i=0)} = \frac{\sum_{i=1}^n Y_i A_i}{\sum_{i=1}^n A_i} - \frac{\sum_{i=1}^n Y_i (1-A_i)}{\sum_{i=1}^n (1-A_i)}.\] However, this estimator may be biased if the responses \(Y_i(0),Y_i(1)\) and \(A\) are dependent. For a concrete example, suppose that a doctor always gives a treatment to patients who look very sick and the treatment only has a small effect. Then even if the treatment is effective, the average outcome of those who received the treatment will still be lower than the average outcome of those without the treatment.

Thus, a common requirement to ensure that the naive estimator converges to the ATE is \[\begin{equation} (Y(0),Y(1))\perp A. \label{eq::ATE} \end{equation}\] The independence of the two versions of outcomes and the treatment assignment. Under this assumption (and some other mild conditions such as the absolute mean exists and there is a positive probability of an individual receiving/not receiving a treatment), we have \[\widehat\tau_{\sf naive}\overset{P}{\rightarrow} \tau.\] In a clinical trial, one scenario to ensure \((Y(0),Y(1))\perp A\) is the randomized controlled trial: every individual is randomly assigned to treatment or control without using any additional information about this individual.

To see why randomization \((Y_i(0),Y_i(1))\perp A\) makes the naive estimator work, note that \(\mathbb{E}(Y|A=a) = \mathbb{E}(Y(a)|A=a)\) for \(a=0,1\). Then under randomization \[\begin{align*} Y(1)\perp A &\Rightarrow \mathbb{E}(Y|A=1) = \mathbb{E}(Y(1)|A=1) = \mathbb{E}(Y(1)),\\ Y(0)\perp A &\Rightarrow \mathbb{E}(Y|A=0) = \mathbb{E}(Y(0)|A=0) = \mathbb{E}(Y(0)). \end{align*}\] Thus, \[\tau = \mathbb{E}(Y(1)) - \mathbb{E}(Y(0)) = \mathbb{E}(Y|A=1)- \mathbb{E}(Y|A=0)\] and the right-hand side is what \(\widehat\tau_{\sf naive}\) is consistently estimating.

Relaxing randomization

In practice, the total randomization on the treatment may be very challenging or even unethical (this basically requires that a doctor has to choose not to treat someone who is very sick when the randomized decision is \(A=0\)). And randomization is often note the case in an observational study. So we would like to think of relaxing the condition \((Y_i(0),Y_i(1))\perp A\).

One possible approach is to use the confounder (confounding variable) \(X\). Namely, in our data, we not only observe the outcome \(Y\) and the treatment \(A\) but also some additional information about each individual, denoted as \(X\) (could be univariate or multivariate). So our data is \[(Y_1,A_1,X_1),\cdots, (Y_n,A_n,X_n).\] In a medical study, \(X\) is often the demographic variables (gender, educational level, ...etc) but it could also be a clinical variable of other diseases or health conditions.

We allow the outcomes \((Y_i(0),Y_i(1))\) and the treatment \(A\) to be dependent, but they are conditionally independent given the observed confounding variable \(X\). Namely, \[\begin{equation} (Y(0),Y(1))\perp A|X. \label{eq::PSM} \end{equation}\] Under this assumption, we have \[\begin{align*} Y(1)\perp A |X&\Rightarrow \mathbb{E}(Y|A=1,X) = \mathbb{E}(Y(1)|A=1,X) = \mathbb{E}(Y(1)|X),\\ Y(0)\perp A |X&\Rightarrow \mathbb{E}(Y|A=0,X) = \mathbb{E}(Y(0)|A=0,X) = \mathbb{E}(Y(0)|X). \end{align*}\]

R-Learner and U-Learner

The R- and U-learners are built on a slightly different generative model. We start with a regression model on the potential outcomes: \[Y(1) = m_1(X) +\epsilon_1(X),\qquad Y(0) = m_0(X) +\epsilon_0(X).\] Also, using the fact that \(Y = A Y(1) + (1-A)Y(0)\), we have \[m(x) = \mathbb{E}(Y|X=x) = \pi_1(x) m_1(x) + (1-\pi_1(x)) m_0(x).\] As a result, we obtain \[\begin{align*} Y - m(X) &= A Y(1) + (1-A)Y(0) -\pi_1(X) m_1(X) - (1-\pi_1(X)) m_0(X) \\ & = Y(0)- m_0(X) + A(Y(1)-Y(0)) - \pi_1(X) (m_1(X) - m_0(X)) \\ & = \epsilon_0(X) + (A-\pi_1(X)) \tau(X) + A(\epsilon_1(X) +\epsilon_0(X))\\ & = (A-\pi_1(X)) \tau(X) + \epsilon(X,A), \end{align*}\] where \(\epsilon(X,A) = \epsilon_0(X) + A(\epsilon_1(X) +\epsilon_0(X))\) is a conditional mean 0 quantity.

The equation \[\begin{equation} Y - m(X) = (A-\pi_1(X)) \tau(X) + \epsilon(X,A) \label{eq::RU} \end{equation}\] is the key to the R- and U-learners. Both learners will first use the whole data to estimate \(\widehat m(x)\) and \(\widehat\pi_1(x)\), which can be done by any conventional approaches. Equation [eq::RU] first appeared in

Robinson, P. M. (1988). Root-N-consistent semiparametric regression. Econometrica: Journal of the Econometric Society, 931-954.

R-learner. Because both \(Y-m(X)\) and \(A-\pi_1(X)\) can be estimated from the data. The R-learner tries to find \(\tau(X)\) by minimizing the following R-learning loss function: \[L_n(f) = \frac{1}{n}\sum_{i=1}^n\left[(Y_i-\widehat m(X_i)) - (A_i - \widehat\pi_1(X_i)) f(X_i)\right]^2.\] We can also add a penalization to the above loss function to regularize the estimator as well. See the following paper for more details:

Nie, X., & Wager, S. (2021). Quasi-oracle estimation of heterogeneous treatment effects. Biometrika, 108(2), 299-319.

As a concrete example, we may model \(f\) as a linear function and use \(L_1\) penalty, which leads to \[\widehat\beta_\lambda = {\sf argmin}_\beta \frac{1}{n}\sum_{i=1}^n\left[(Y_i-\widehat m(X_i)) - (A_i - \widehat\pi_1(X_i)) \beta^TX_i\right]^2 + \lambda \|\beta\|_1.\] The CATE estimator is \(\widehat\beta_\lambda^Tx\). Note that here we include the constant \(1\) in the covariate (intercept term of the linear model).

U-learner. The U-learner is a method proposed in [KSBY2019]. Instead of viewing the problem as a minimization problem, the U-learner further defines \[U= \frac{Y-m(X)}{A-\pi_1(X)}\] and uses the fact that \(\mathbb{E}(U|X) = \tau(X)\). Thus, based on the nuisance estimators \(\widehat m(x)\) and \(\widehat\pi_1(x)\), we then compute \[\widehat U_i = \frac{Y_i-\widehat m(X_i)}{A_i-\widehat\pi_1(X_i)}\] and then regress \(\widehat U_i\) with \(X_i\) to obtain \(\widehat\tau(x)\).

Graphical models

Some useful references:

Now we discuss another framework for causal inference: the graphical model approach. The graphical model approach is particularly useful when the treatment variable is non-binary and even continuous. In the graphical model, the causal effect is defined through the do operator.

Suppose that we have several variables \(V_1,\cdots, V_d\) of interest and we use the DAG to specify the underlying generating model (Bayesian network). To simplify the problem, suppose that we are interested in estimating the causal effect \(V_1\) on the variable \(V_2\). We often relabel the variables as \(V_1=X\), \(V_2=Y\), and make the rest of them as \(Z_1,\cdots, Z_m\), where \(m=d-2\). With this notation, the parameter of interest is the causal effect from \(X\) on \(Y\). Let \(p(x,y,z_1,\cdots, z_m)\) be the joint PDF (it can be generalized to PMF as well) and \(G = (V, E)\) where \(V = (V_1,\cdots, V_d)= (X,Y,Z_1,\cdots, Z_m)\) and \(E_{ij} = 1\) if there is a directed arrow from node \(i\) to node \(j\). The DAG implies \[\begin{equation} p(x,y,z_1,\cdots, z_m) = p(x|{\sf PA}_x)p(y|{\sf PA}_y)\prod_{j=1}^m p(z_j|{\sf PA}_{z_j}), \label{eq::DAG1} \end{equation}\] where \({\sf PA}_v\) denotes the set of parent nodes of variable \(v\).

Defining the causal effect from \(X\) on \(Y\) is not easy because they may be interacting with variables \(Z_1,\cdots, Z_m\). The do operator provides a solution to this. The do operator defines the causal effect using the conditional PDF \[\begin{equation} p(y|{\bf do}(x)) \equiv p(y|{\bf do}(X=x)). \label{eq::do1} \end{equation}\] We often define \(\tau(x) = \frac{\partial}{\partial x}\mathbb{E}(Y|{\bf do}(X)=x) = \frac{\partial}{\partial x} \int y p(y|{\bf do}(x))dy\) as the causal effect on \(Y\) from \(X\). Note that in general, \[p(y|{\bf do}(x)) \neq p(y|x)\] except for the simple case where there is only an arrow \(X\rightarrow Y\) and no other arrows toward \(Y\).

To see how do-calculus is associated with the potential outcome model, consider the DAG on the left of Figure 1. The conditional mean of \(Y\) given \({\bf do}(x)\) is \[\mathbb{E}(Y|{\bf do} (x)) = \int y p(y|{\bf do}(x)) dy = \int y p(y|x,z) p(z)dzdy = \int \left[\int yp(y|x,z)dy \right] p(z)dz = \int m(x,z) p(z)dz.\] When \(X\) is binary (a binary treatment), the potential outcome model with confounder being \(Z\) implies \(\mathbb{E}(Y(x)) = \int m(x,z) p(z)dz\); see equation [eq::relax]. Therefore, we conclude that \[\mathbb{E}(Y|{\bf do} (x)) = \mathbb{E}(Y(x)).\] Thus, do-calculus and the potential outcome model lead to the same conclusion.

The conditional PDF \(p(y|{\bf do}(X=x))\) is interpreted as: we change the system in a way that the variable \(X\) is set to \(x\), which leads to a density function of \(Y\) and this density function is \(p(y|{\bf do}(x))\).

Given a DAG \(G= (V,E)\) where \(V = (X,Y,Z_1,\cdots, Z_m)\), the do operation defines a new DAG \(G' = (V,E) = G({\bf do}(x)) = (V,E({\sf do}(x)))\) such that all directed arrows to \(X\) are removed. This leads to a new factorization of the joint PDF: \[\begin{equation} p({\bf do} (x), y,z_1,\cdots, z_m) = p ({\bf do} (x)) p(y|{\sf PA}_y)\prod_{j=1}^m p(z_j|{\sf PA}_{z_j}) \label{eq::DAG2} \end{equation}\] or the corresponding conditional density \[\begin{equation} p(y,z_1,\cdots, z_m|{\bf do} (x)) = p(y|{\sf PA}_y)\prod_{j=1}^m p(z_j|{\sf PA}_{z_j}). \label{eq::do2} \end{equation}\] Equation [eq::do2] is known as g-formula (by J. Robins), or truncated factorization formula (by J. Pearl).

If we use the DAG to interpret the result, the new DAG \(G'\) preserves all causal effects except for the ones that are affecting \(X\). This is exactly how we (commonly) think about the causal effect due to \(X\)–we keep the entire system as is except we add an intervention at variable \(X\) that sets it to be \(x\).

\(\qquad\qquad\qquad\)

Left: original DAG \(G\). Right: the DAG after the \({\bf do}(x)\) operation.

The power of equation [eq::do2] is that the left-hand-side \(p(y,z_1,\cdots, z_m|{\bf do} (x))\) is the conditional density due to the do operation \({\bf do}(x)\), which is a theoretical entity, and the right-hand-side is what we can identify using the original DAG. With equation [eq::do2], we can identify equation [eq::do1] using \[p(y|{\bf do}(x)) = \int p(y,z_1,\cdots, z_m|{\bf do} (x)) dz_1\cdots dz_m = \int p(y|{\sf PA}_y)\prod_{j=1}^m p(z_j|{\sf PA}_{z_j})dz_1\cdots dz_m\] and \(\tau(x)\) accordingly.

Note that in this case, we often assume that the DAG is known so we can estimate all the conditional densities \(p(z_j|{\sf PA}_{z_j})\) (and \(p(y|{\sf PA}_{y})\)) using the data. Equation [eq::do2] shows that we can identify the causal effect from the data.

Example 1. In the example of Figure 1, we have an original DAG in the left and a new DAG due to \({\bf do}(x)\). The original DAG implies a factorization of the joint PDF \[p(x,y,z) = p(x|z)p(y|x,z)p(z).\] All the three conditionals can be estimated/identified from the data if we know this DAG in advance. Using the g-formula (equation [eq::do2]), the \({\bf do}(x)\) operation leads to a conditional density \[p(y,z|{\bf do}(x)) = p(y|x,z)p(z),\] which is still identifiable from the three conditionals provided in the original DAG.

\(\qquad\qquad\qquad\)

A more complicated DAG. Left: original DAG \(G\). Right: the DAG after the \({\bf do}(x)\) operation.

Example 2. In Figure 2, we provide a more complicated example where there are \(6\) variables. The left panel displays the original DAG and the right panel displays the DAG after a \({\bf do}(x)\) operation. The original DAG implies the following factorization \[p(x,y,z_1,z_2,z_3,z_4) = p(z_1)p(z_2|z_1)p(x|z_1,z_2)p(z_3|x)p(z_4|z_1,z_3)p(y|x,z_2,z_3).\] All these conditionals are identifiable from the data. After the \({\bf do}(x)\) operation, the conditional density is \[p(y,z_1,z_2,z_3, z_4|{\bf do}(x)) = p(z_1)p(z_2|z_1)p(z_3|x)p(z_4|z_1,z_3)p(y|x,z_2,z_3).\] Each element in the right-hand side is identifiable so we can identify the entire conditional density. Note that if we are only interested in \(p(y|{\bf do}(x))\), we can write it as \[\begin{align*} p(y|{\bf do}(x)) &= \int p(y,z_1,z_2,z_3, z_4|{\bf do}(x))dz_1dz_2dz_3dz_4\\ &=\int p(y|x,z_2,z_3)p(z_1)p(z_2|z_1)p(z_3|x)\left(\int p(z_4|z_1,z_3) dz_4\right) dz_1dz_2dz_3\\ &=\int p(y|x,z_2,z_3)p(z_1)p(z_2|z_1)p(z_3|x) dz_1dz_2dz_3. \end{align*}\] So we only need to estimate these 4 conditionals. In a sense, we do not need to consider estimating any effect of \(Z_4\) (since it does not have a causal effect onto \(Y\)). Using the DAG induced by the do operator, we have \[Z_3\perp Z_1,Z_2|X,\quad Y\perp X,Z_1|Z_2,Z_3, \quad Z_2\perp X|Z_1,\] we can further write the above equality as \[\begin{align*} p(y|{\bf do}(x))&=\int p(y|x,z_2,z_3)p(z_1)p(z_2|z_1)\underbrace{p(z_3|x)}_{=p(z_3|z_2,x)} dz_1dz_2dz_3\\ & = \int p(y,z_3|x,z_2)p(z_1)p(z_2|z_1)dz_1dz_2dz_3\\ & = \int \underbrace{p(y|x,z_2)}_{=p(y|x,z_1,z_2)}p(z_1)\underbrace{p(z_2|z_1)}_{=p(z_2|z_1,x)}dz_1dz_2\\ & = \int p(y,z_2|x,z_1)p(z_1)dz_2dz_1\\ & = \int p(y|x,z_1)p(z_1)dz_1. \end{align*}\] The conditional density due to the do operator is essentially the conditional density after adjusting \(p(z_1)\), so variable \(Z_1\) is called the adjustment set; see Definition 3.6 of the following paper:

Maathuis, M. H., & Colombo, D. (2015). A generalized back-door criterion. The Annals of Statistics, 43(3), 1060-1088.

The adjustment set offers an elegant way to further simplify the g-formula – the conditional density of \(Y\) given the do operator \({\bf do}(x)\) is the same as we adjust the conditional density of \(Y\) given \(X\) and the variables in the adjustment sets. In a sense, the adjustment set represents the possible sources of interaction from other variables onto the causal effect from \(X\) onto \(Y\). So we have to adjust for these variables to obtain the desired causal effect.

Formally, speaking the collection of random variables \(W \subset\{Z_1,\cdots, Z_n\}\) is called an adjustment set relative to \((X,Y)\) in the DAG if \[p(y|{\bf do}(x)) = \begin{cases} p(y|x),\qquad&\mbox{if } W=\emptyset\\ \int p(y|x, w)p(w)dw,\qquad&\mbox{if } W\neq\emptyset \end{cases}\] In a sense, the adjustment set behaves like confounders when we are investigating the causal effect from \(X\) to \(Y\).

Note that all the above analysis relies on the fact that we know the DAG in advance. This is possible if we have additional scientific knowledge about each variable. However, in a general observational study, all we can estimate (using the data) is conditional independence, which is an undirected graph. We may have some partial knowledge about each edge, leading to a mixed graph (a graph with some directed and some undirected edges). Here are some papers related to finding the adjustment sets for different types of graphs (different types of graphs representing situations where we have different prior knowledge about the relations among variables):

1. Perković, E., Textor, J., Kalisch, M., & Maathuis, M. H. (2018). Complete graphical characterization and construction of adjustment sets in Markov equivalence classes of ancestral graphs. The Journal of Machine Learning Research, 18(1), 8132-8193.
2. Perković, E., Textor, J., Kalisch, M., & Maathuis, M. H. (2015). A complete generalized adjustment criterion. In Uncertainty in Artificial Intelligence (pp. 682-691). AUAI Press.

Remark: linear structural equation modeling (Linear SEM). A popular method that uses the DAG to make inference is the structural equation modeling (SEM). In the simplest form, the SEM assumes a linear effect for every arrow in the DAG. Suppose \(X\rightarrow Y\) and \(U\rightarrow X\) and \(U\rightarrow Y\), then an SEM will be \[Y = \alpha_Y + \beta X+\gamma U+\epsilon_Y,\quad X = \alpha_X + \eta U+\epsilon_X\] and \(\epsilon_X,\epsilon_Y\sim N(0,\sigma^2)\) with \(U\sim p(u)\). If we are interested in the causal effect from \(X\) onto \(Y\), \(\beta\) will be the parameter of interest. So the question is: how do we properly apply the regression to obtain a consistent estimate of \(\beta\). In general, we need to observe \(X,Y,U\) to properly estimate \(\beta\) (using multiple linear regression). Note that \(U\) is called the confounder for the causal effect from \(X\) to \(Y\).

Remark: instrumental variables and linear SEM. Under the setup of the previous remark, if \(U\) is unobserved (unobserved confounder problem), then we cannot identify the causal effect \(\beta\). However, IV (instrumental variable; Section [sec::IV]) offers a solution to this problem. In the DAG, the IV is written in a different form. Suppose that we do not observe \(U\) (so we cannot identify the causal effect \(\beta\)) but we observe another variable \(Z\) such that

there is an arrow \(Z\rightarrow X\) and no other arrows related to \(Z\).

This variable \(Z\) is called a valid IV. The above condition on the graphical structure of \(Z\) can be viewed as a generalization of the IV conditions in Section [sec::IV] to a graphical model framework. Under the linear SEM, \[X = \alpha_X + \eta U + \xi Z+\epsilon_X.\] Putting this into \(Y\), we obtain \[\begin{align*} Y &= \alpha_Y + \beta (\alpha_X + \eta U + \xi Z+\epsilon_X)+\gamma U+\epsilon_Y\\ & = \alpha' + \eta' U + \beta \xi Z + \epsilon'. \end{align*}\] Thus, regressing \(Y\) with \(Z\) leads to a slope \(\beta\xi\) (note that \(Z\perp U\) since there are no mutual arrows). Regressing \(X\) with \(Z\) yields the slope \(\xi\) (again using the fact that \(Z\perp U\)). So we can estimate \(\beta\) by the ratio of the two regression coefficient even if we do not observe the confounder \(U\).

Continuous treatment

When the treatment variable \(A\) is continuous or contains infinite amount of possible values, the causal inference becomes very challenging. In particular, the potential outcomes \(Y(a)\) will be hard to characterize since there could be many of them,

Let \(Y\in \mathbb{R}\) be the outcome of interest, \(A\in\mathbb{R}\) be the treatment variable, and \(X\in \mathbb{R}^d\) be the confounders. Similar to the conventional setup, our data is the collection of \[(Y_1,A_1,X_1),\cdots, (Y_n,A_n,X_n).\] In the continuous treatment problem, the parameter of interest is the dose-response curve \[m(a) = \mathbb{E}(Y(a))\] and we often work with the following conditions:

These two conditions are essentially the same conditions as in the binary treatment case.

Under (C1-2), we can rewrite the dose-response curve as the following form: \[\begin{equation} \begin{aligned} m(a) & = \mathbb{E}(\mathbb{E}(Y|A=a, X))\\ & = \mathbb{E}(\mu(a,X)), \end{aligned} \label{eq::DC} \end{equation}\] where \[\mu(a,x) = \mathbb{E}(Y|A=a,X=x)\] is the regression function. To estimate the dose-response curve, there are three popular approaches.

Direct nonparametric estimation

Regression adjustment. Equation \(\ref{eq::DC}\) suggests a plug-in approach to estimate \(m(a)\): \[\widehat m(a) = \frac{1}{n}\sum_{i=1}^n \widehat\mu(a,X_i),\] where \(\widehat\mu(a,x)\) is a regression estimator, which is the nuisance parameter in this case. This is essentially a regression adjustment method.

Inverse probability weighting. While it is not easy to think about the IPW approach for continuous treatment because the conventional IPW uses the indicator function \(I(A=a)\), which will almost always be \(0\), it is possible to use a kernelized method to construct an IPW estimator. The IPW (kernel smoothing) method is \[\widehat m_h(a) = \frac{1}{n}\sum_{i=1}^n \frac{1}{h}K\left(\frac{A_i-a}{h}\right)\frac{Y_i}{\widehat p(a|X_i)},\] where \(\widehat p(a|x)\) is the estimated conditional PDF of \(A=a\) given \(X=x\) and \(K\) is a smoothing kernel and \(h>0\) is a smoothing bandwidth. Essentially, we replace the indicator function \(I(A=a)\) by its kernelized version \(\frac{1}{h}K\left(\frac{A_i-a}{h}\right)\). The probability weighting (propensity score) is replaced by the conditional PDF \(p(a|x)\). More details about this idea can be found in

Huber, M., Hsu, Y. C., Lee, Y. Y.,& Lettry, L. (2020). Direct and indirect effects of continuous treatments based on generalized propensity score weighting. Journal of Applied Econometrics, 35(7), 814-840

You can also combine both estimators to form a doubly-robust estimator as well.

Positivity conditions. In both regression adjustment and IPW, the positivity conditions become very crucial. The two nuisances \(\widehat\mu(a,x)\) and \(\widehat p(a|x)\) will be evaluated at every \(X_i\) with \(A=a\). Thus, both estimators rely heavily on the uniform consistency of the two nuisances on the region \(\{a\} \times {\sf supp}(X)\). It is possible to bypass the positivity condition using a method called the integral estimator; see the following paper:

Zhang, Y., Chen, Y. C., & Giessing, A. (2024). Nonparametric Inference on Dose-Response Curves Without the Positivity Condition. arXiv preprint arXiv:2405.09003.

Pseudo-outcome approach

Another approach to estimate the dose-response curve is the pseudo-outcome approach introduced in the following paper:

Kennedy, E. H., Ma, Z., McHugh, M. D., & Small, D. S. (2017). Non-parametric methods for doubly robust estimation of continuous treatment effects. Journal of the Royal Statistical Society Series B: Statistical Methodology, 79(4), 1229-1245.

A challenge of estimating the dose-response curve is that the parameter of interest \(m(a)\) is not pathwise differentiable, making the analysis on its efficient influence function very challenging. However, it is still possible to construct an estimator with double robustness property by the use of pseudo-outcomes.

Let \(\widehat\mu(a,x), \widehat p(a|x)\) be the estimated regression and conditional PDF as in the previous sections. The pseudo-outcome of \((Y,A,X)\) is \[\widehat Y = \frac{Y - \widehat\mu(A,X)}{\widehat p(A|X)} \cdot \frac{1}{n}\sum_{j=1}^n \widehat p(A|X_j) + \frac{1}{n}\sum_{j=1}^n \widehat\mu(A, X_j).\] Applying this to every observation, we obtain \[\widehat Y_1,\cdots, \widehat Y_n.\] Then we estimate \(m(a)\) by applying a nonparametric regression of \(\widehat Y\) versus \(A\).

The pseudo-outcome approach enjoys the double robustness property – we just need either \(\widehat\mu\) or \(\widehat p\) to be correctly specified to obtain consistency. However, this method still heavily relies on the positivity condition.

Marginal structural modeling

The marginal structural modeling is a popular approach to deal with problems involving multiple or even infinite potential outcomes. The idea is very simple: we put a parametric model on the dose-response curve, i.e., \[m(a) \equiv \theta(a; \gamma),\] where \(\gamma\) is the underlying parameter.

Examples:

  1. Linear Model. \(\mathbb{E}(Y(a)) = \theta(a; \gamma) =\gamma_0+ \gamma_1^T a .\)

  2. Log-linear model. \(\log \mathbb{E}(Y(a)) = \gamma_0+\gamma_1^Ta\) or equivalently, \[\theta(a;\gamma) = \exp(\gamma_0+\gamma_1^Ta).\]

  3. Logistic model. For \(Y\in[0,1]\), or a binary \(Y\), we may use \(\log \left\{\frac{\mathbb{E}(Y(a))}{1-\mathbb{E}(Y(a))}\right\} = \gamma_0+\gamma_1^Ta\) or equivalently, \[\theta(a;\gamma) = \frac{\exp(\gamma_0+\gamma_1^Ta)}{1-\exp(\gamma_0+\gamma_1^Ta)}.\]

To estimate the underlying parameter \(\gamma\), there are two possible approaches–inverse probability weighting and regression adjustment.

Inverse probability weighting. The IPW approach estimates \(\widehat\gamma\) via the following procedure:

  1. Estimate \(\widehat p(a), \widehat p(a|x)\) from \((A_1,X_1),\cdots, (A_n,X_n)\).

  2. Construct the estimating equation \[\begin{equation} \widehat\Psi(\gamma; Y,A,X) = \frac{\widehat p(A)}{\widehat p(A|X)} (Y-\theta(A;\gamma)) s(A;\gamma),\qquad s(a;\gamma) = \frac{\partial}{\partial \gamma} \theta(a;\gamma). \label{eq::MSM1} \end{equation}\]

  3. Find \(\widehat\gamma_n\) by solving \[0 = \frac{1}{n}\sum_{i=1}^n\widehat\Psi(\widehat\gamma_n; Y_i,A_i,X_i).\]

  4. Obtain the estimator \(\widehat\theta(a) = \theta(a;\widehat\gamma_n)\).

Essentially, the IPW is derived from the usual M-estimator when we assume a parametric model with inverse probability weighting. If we ignore \(\widehat p(A)\) in equation \(\ref{eq::MSM1}\), the whole estimating equation behaves just like the score equation. The multiplier \(\widehat p(A)\) is to stabilize the estimation procedure, see the following note for more discussion: https://faculty.washington.edu/yenchic/short_note/note_msmc.pdf.

Regression adjustment (g-computation). Alternatively, we may put a parametric model of \(m(a)\) indirectly from a model on \(\mu(a,x)\) and estimate it accordingly. Specifically, we place a model \(\mu(a,x;\beta)\) and then estimate \(\beta\) by a least squares approach, leading to \(\widehat\beta\). We can then obtain an estimator of \(m(a)\) accordingly.

Here is what we will do in practice. We first compute the estimator \[\widehat\beta_n = {\sf argmin}_\beta \frac{1}{n}\sum_{i=1}^n (Y_i- \mu(A_i,X_i;\beta))^2.\] Then we construct the estimator of the MSM via \[\widehat\theta(a) = \frac{1}{n}\sum_{i=1}^n \mu(a,X_i;\widehat\beta_n).\]