<jupyter_start><jupyter_text>The Deconfounder in Action In this notebook, we are going to see **the deconfounder in action**.We will perform **causal inference** with the deconfounder on a **breast cancer** dataset.**Goal:**To convince all of us that the deconfounder is **easy** to use! The **deconfounder** operates in three steps:1.   **Fit** a factor model to the assigned causes; it leads to a candidate substitute confounder.2.   **Check** the factor model with a predictive check.3.   **Correct** for the substitute confounder in a causal inference.Let's get started!  Getting ready to work!<jupyter_code>!pip install tensorflow_probability
import tensorflow as tf
import numpy as np
import numpy.random as npr
import pandas as pd
import tensorflow as tf
import tensorflow_probability as tfp
import statsmodels.api as sm

from tensorflow_probability import edward2 as ed
from sklearn.datasets import load_breast_cancer
from pandas.plotting import scatter_matrix
from scipy import sparse, stats
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression 
from sklearn.metrics import confusion_matrix, classification_report, roc_auc_score, roc_curve

import matplotlib
matplotlib.rcParams.update({'font.sans-serif' : 'Helvetica',
                            'axes.labelsize': 10,
                            'xtick.labelsize' : 6,
                            'ytick.labelsize' : 6,
                            'axes.titlesize' : 10})
import matplotlib.pyplot as plt

import seaborn as sns
color_names = ["windows blue",
               "amber",
               "crimson",
               "faded green",
               "dusty purple",
               "greyish"]
colors = sns.xkcd_palette(color_names)
sns.set(style="white", palette=sns.xkcd_palette(color_names), color_codes = False)
!pip show tensorflow
!pip show tensorflow_probability
# set random seed so everyone gets the same number
import random
randseed = 123
print("random seed: ", randseed)
random.seed(randseed)
np.random.seed(randseed)
tf.set_random_seed(randseed)<jupyter_output>random seed:  123<jupyter_text>The scikit-learn breast cancer dataset *   It is a data set about **breast cancer**. *   We are interested in how tumor properties **affect** cancer diagnosis. *   The **(multiple) causes** are tumor properties, e.g. sizes, compactness, symmetry, texture. *   The **outcome** is tumor diagnosis, whether the breast cancer is diagnosed as malignant or benign.<jupyter_code>data = load_breast_cancer()
print(data['DESCR'])<jupyter_output>Breast Cancer Wisconsin (Diagnostic) Database
=============================================

Notes
-----
Data Set Characteristics:
    :Number of Instances: 569

    :Number of Attributes: 30 numeric, predictive attributes and the class

    :Attribute Information:
        - radius (mean of distances from center to points on the perimeter)
        - texture (standard deviation of gray-scale values)
        - perimeter
        - area
        - smoothness (local variation in radius lengths)
        - compactness (perimeter^2 / area - 1.0)
        - concavity (severity of concave portions of the contour)
        - concave points (number of concave portions of the contour)
        - symmetry 
        - fractal dimension ("coastline approximation" - 1)

        The mean, standard error, and "worst" or largest (mean of the three
        largest values) of these features were computed for each image,
        resulting in 30 features.  For instance, field 3 is Mean Radius, field
        13 is [...]<jupyter_text>***For simplicity, we will work with the first 10 features, i.e. the mean radius/texture/perimeter***/....<jupyter_code>num_fea = 10
df = pd.DataFrame(data["data"][:,:num_fea], columns=data["feature_names"][:num_fea])
df.shape
df.head()
dfy = data["target"]
dfy.shape, dfy[:100] # binary outcomes<jupyter_output><empty_output><jupyter_text>Preparing the dataset for the deconfounder  Only one step of preprocessing needed! We need to get rid of the highly correlated causes **Why** do we need to get rid of highly correlated causes?If two causes are **highly correlated**, a valid substitute confounder will largely **inflate the variance** of causal estimates downstream.This phenomenon is **closely related to** the variance inflation phenomenon in linear regression. ***A more technical explanation (ignorable)***Think of the extreme case where two causes are perfectly collinear $A_1 = 5A_2$. The only random variable Z that $$A_1 \perp A_2 | Z,$$$(A_1, A_2)$ **must** be a **deterministic function** of Z. For example, $Z = A_1$ or $Z = A_2$.Such a substitute confounder Z **breaks one of the conditions** the deconfounder requires. See "***A note on overlap***" in the theory section of the paper. **How** do we get rid of highly correlated causes?*   We first make a **scatter plot** of **all pairs** of the causes.*   It reveals which causes are **highly correlated**.*   We will **exclude** these highly correlated causes by hand.<jupyter_code>sns.pairplot(df, size=1.5)
# perimeter and area are highly correlated with radius
fea_cols = df.columns[[(not df.columns[i].endswith("perimeter")) \
                     and (not df.columns[i].endswith("area")) \
                     for i in range(df.shape[1])]]
dfX = pd.DataFrame(df[fea_cols])

print(dfX.shape, dfy.shape)<jupyter_output>(569, 8) (569,)<jupyter_text>How does the dataset look like after preprocessing?<jupyter_code># The causes
dfX.head()
# The outcome
dfy[:25]<jupyter_output><empty_output><jupyter_text>The dataset is ready. Let's do causal inference with the deconfounder!  Step 1: Fit a factor model to the assigned causes; it leads to a substitute confounder.  We start with trying out a random factor model. How about a probabilistic PCA model? The matrix of assigned causes $X$*   It has N=569 rows and D=8 columns. *   N is the number of subjects/data points.*   D is the number of causes/data dimension.  Step 1.1: Some chores first...  Standardize the dataThis step is optional to the deconfounder. It only makes finding a good probabilistic PCA model easier.<jupyter_code># dfX.std()
# standardize the data for PPCA
X = np.array((dfX - dfX.mean())/dfX.std())<jupyter_output><empty_output><jupyter_text>Then holdout some data! We will later need to check the factor model with some heldout data.So let's holdout some now.<jupyter_code># randomly holdout some entries of X
num_datapoints, data_dim = X.shape

holdout_portion = 0.2
n_holdout = int(holdout_portion * num_datapoints * data_dim)

holdout_row = np.random.randint(num_datapoints, size=n_holdout)
holdout_col = np.random.randint(data_dim, size=n_holdout)
holdout_mask = (sparse.coo_matrix((np.ones(n_holdout), \
                            (holdout_row, holdout_col)), \
                            shape = X.shape)).toarray()

holdout_subjects = np.unique(holdout_row)

x_train = np.multiply(1-holdout_mask, X)
x_vad = np.multiply(holdout_mask, X)<jupyter_output><empty_output><jupyter_text>Step 1.2: We are ready to fit a probabilistic PCA model to x_train. This step of "**fitting** a factor model" involves **inferring latent variables** in probability models. We will rely on **Tensorflow Probability**, a library for probabilistic reasoning and statistical analysis in TensorFlow.There are many **other probabilistic programming toolboxes** for fitting factor models, e.g. Pyro, Stan. Some of the latent variable models can also be fit with **scikit-learn**. We are free to use any of these with the deconfounder! **What does a probabilistic PCA model look like?***  Probabilistic PCA is a dimensionality reduction technique. It models data with a lower dimensional latent space.*  We consider the assigned causes of the $n$th subject. We write it as $\mathbf{x}_n$, which is a $D=8$ dimensional vector.*  The probabilistic PCA assumes the following data generating process for each $\mathbf{x}_n$, $n = 1, ..., N$:\begin{equation*}\mathbf{z}_{n} \stackrel{iid}{\sim} N(\mathbf{0}, \mathbf{I}_K),\end{equation*}\begin{equation*}\mathbf{x}_n \mid \mathbf{z}_n\sim N(\mathbf{z}_n\mathbf{W}, \sigma^2\mathbf{I}_D).\end{equation*}*  We construct a $K$-dimensional substitute confounder $\mathbf{z}_{n}$ for each subject $n$, $n = 1, ..., N$. *  Each $\mathbf{z}_{n}$ is a $K$-dimensional latent vector, $n = 1, ..., N$.<jupyter_code>def ppca_model(data_dim, latent_dim, num_datapoints, stddv_datapoints):
    w = ed.Normal(loc=tf.zeros([latent_dim, data_dim]),
                scale=tf.ones([latent_dim, data_dim]),
                name="w")  # parameter
    z = ed.Normal(loc=tf.zeros([num_datapoints, latent_dim]),
                scale=tf.ones([num_datapoints, latent_dim]), 
                name="z")  # local latent variable / substitute confounder
    x = ed.Normal(loc=tf.multiply(tf.matmul(z, w), 1-holdout_mask),
                scale=stddv_datapoints * tf.ones([num_datapoints, data_dim]),
                name="x")  # (modeled) data
    return x, (w, z)

log_joint = ed.make_log_joint_fn(ppca_model)<jupyter_output><empty_output><jupyter_text>**Let's fit a probabilistic PCA model with K = 2 latent dimensions.**<jupyter_code>latent_dim = 2
stddv_datapoints = 0.1

model = ppca_model(data_dim=data_dim,
                   latent_dim=latent_dim,
                   num_datapoints=num_datapoints,
                   stddv_datapoints=stddv_datapoints)<jupyter_output><empty_output><jupyter_text>The cell below implements **variational inference** for probabilistic PCA  in tensorflow probability.You are free to fit the probabilistic PCA in your favourite ways with your favourite package. Note: approximate inference is perfectly fine!It is orthogonal to our discussion around the deconfounder.Let's **ignore** that for now (and forever).<jupyter_code>def variational_model(qw_mean, qw_stddv, qz_mean, qz_stddv):
    qw = ed.Normal(loc=qw_mean, scale=qw_stddv, name="qw")
    qz = ed.Normal(loc=qz_mean, scale=qz_stddv, name="qz")
    return qw, qz


log_q = ed.make_log_joint_fn(variational_model)

def target(w, z):
    """Unnormalized target density as a function of the parameters."""
    return log_joint(data_dim=data_dim,
                   latent_dim=latent_dim,
                   num_datapoints=num_datapoints,
                   stddv_datapoints=stddv_datapoints,
                   w=w, z=z, x=x_train)

def target_q(qw, qz):
    return log_q(qw_mean=qw_mean, qw_stddv=qw_stddv,
               qz_mean=qz_mean, qz_stddv=qz_stddv,
               qw=qw, qz=qz)


qw_mean = tf.Variable(np.ones([latent_dim, data_dim]), dtype=tf.float32)
qz_mean = tf.Variable(np.ones([num_datapoints, latent_dim]), dtype=tf.float32)
qw_stddv = tf.nn.softplus(tf.Variable(-4 * np.ones([latent_dim, data_dim]), dtype=tf.float32))
qz_stddv = tf.nn.softplus(tf.Variable(-4 * np.ones([num_datapoints, latent_dim]), dtype=tf.float32))

qw, qz = variational_model(qw_mean=qw_mean, qw_stddv=qw_stddv,
                           qz_mean=qz_mean, qz_stddv=qz_stddv)


energy = target(qw, qz)
entropy = -target_q(qw, qz)

elbo = energy + entropy


optimizer = tf.train.AdamOptimizer(learning_rate = 0.05)
train = optimizer.minimize(-elbo)

init = tf.global_variables_initializer()

t = []

num_epochs = 500

with tf.Session() as sess:
    sess.run(init)

    for i in range(num_epochs):
        sess.run(train)
        if i % 5 == 0:
            t.append(sess.run([elbo]))

        w_mean_inferred = sess.run(qw_mean)
        w_stddv_inferred = sess.run(qw_stddv)
        z_mean_inferred = sess.run(qz_mean)
        z_stddv_inferred = sess.run(qz_stddv)
        
print("Inferred axes:")
print(w_mean_inferred)
print("Standard Deviation:")
print(w_stddv_inferred)

plt.plot(range(1, num_epochs, 5), t)
plt.show()

def replace_latents(w, z):

    def interceptor(rv_constructor, *rv_args, **rv_kwargs):
        """Replaces the priors with actual values to generate samples from."""
        name = rv_kwargs.pop("name")
        if name == "w":
            rv_kwargs["value"] = w
        elif name == "z":
            rv_kwargs["value"] = z
        return rv_constructor(*rv_args, **rv_kwargs)

    return interceptor<jupyter_output>Inferred axes:
[[ 1.0657768   0.7117974   0.06351132  0.5483043   0.7629369   0.8615003
   0.10290438 -0.40677893]
 [-0.20814618 -0.34985837  0.8877295   0.7185301   0.48308808  0.34881124
   0.8212925   1.0651304 ]]
Standard Deviation:
[[0.00943145 0.0087002  0.00819345 0.00831227 0.00870339 0.01042197
  0.01169398 0.01013054]
 [0.01030399 0.01017774 0.00799988 0.00906306 0.00936898 0.00806279
  0.0075188  0.00970751]]<jupyter_text>So we just played some **magic** to **fit the probabilistic PCA** to the matrix of assigned causes $\mathbf{X}$.**The only important thing here is: **We have **inferred** the latent variables $\mathbf{z}_n, n=1, ..., N$ and the parameters $\mathbf{W}$.Specifically, we have obtained from this step```w_mean_inferred,w_stddv_inferred,z_mean_inferred,z_stddv_inferred.```  Step 2: Check the factor model with a predictive check. Now we are ready to **check** the probabilistic PCA model.The checking step is **very important** to the deconfounder. Pleeeeeze **always** check the factor model! **How** do we perform the predictive check?1.  We will **generate** some replicated datasets for the heldout entries.2.  And then **compare** the replicated datasets with the original dataset on the heldout entries.3.  If they **look similar**, then we are good to go.  Step 2.1: We generate some replicated datasets first. *  We will start with generating some **replicated datasets** from the predictive distribution of the assigned causes $X$:\begin{align}  p(\mathbf{X^{rep}_{n,heldout}} \,|\, \mathbf{X_{n, obs}}) =  \int p(\mathbf{X_{n, heldout}} \,|\, \mathbf{z}_n) p(\mathbf{z_n} \,|\, \mathbf{X}_{n, obs}) \mathrm{d} \mathbf{z_n}.\end{align}*  That is, we generate these datasets from a probabilistic PCA model given the **inferred** latent variables $\hat{p}(\mathbf{z}_n)$ and $\hat{p}(\mathbf{W})$:\begin{equation*}\mathbf{z}_{n} \sim \hat{p}(\mathbf{z}_n),\end{equation*}\begin{equation*}\mathbf{W} \sim \hat{p}(\mathbf{W}),\end{equation*}\begin{equation*}\mathbf{x}_n \mid \mathbf{z}_n\sim N(\mathbf{z}_n\mathbf{W}, \sigma^2\mathbf{I}_D).\end{equation*}*  These replicated datasets tell us what the assigned causes $X$ **should look like** if it is indeed generated by the fitted probabilistic PCA model.<jupyter_code>n_rep = 100 # number of replicated datasets we generate
holdout_gen = np.zeros((n_rep,*(x_train.shape)))

for i in range(n_rep):
    w_sample = npr.normal(w_mean_inferred, w_stddv_inferred)
    z_sample = npr.normal(z_mean_inferred, z_stddv_inferred)

    with ed.interception(replace_latents(w_sample, z_sample)):
        generate = ppca_model(
            data_dim=data_dim, latent_dim=latent_dim,
            num_datapoints=num_datapoints, stddv_datapoints=stddv_datapoints)

    with tf.Session() as sess:
        x_generated, _ = sess.run(generate)

    # look only at the heldout entries
    holdout_gen[i] = np.multiply(x_generated, holdout_mask)<jupyter_output><empty_output><jupyter_text>Step 2.2: Then we compute the test statistic on both the original and the replicated dataset. *   We use the **test statistic** of **expected heldout log likelihood**:\begin{align}  t(\mathbf{X_{n,heldout}}) = \mathbb{E}_{\mathbf{Z}, \mathbf{W}}[{\log p(\mathbf{X_{n,heldout}} \,|\, \mathbf{Z}, \mathbf{W}) \,|\,  \mathbf{X_{n,obs}}}].\end{align}*   We calculate this test statistic **for each $n$** and for **both** the **original** dataset $\mathbf{X_{n,heldout}}$ and the **replicated** dataset $\mathbf{X^{rep}_{n,heldout}}$.<jupyter_code>n_eval = 100 # we draw samples from the inferred Z and W
obs_ll = []
rep_ll = []
for j in range(n_eval):
    w_sample = npr.normal(w_mean_inferred, w_stddv_inferred)
    z_sample = npr.normal(z_mean_inferred, z_stddv_inferred)
    
    holdoutmean_sample = np.multiply(z_sample.dot(w_sample), holdout_mask)
    obs_ll.append(np.mean(stats.norm(holdoutmean_sample, \
                        stddv_datapoints).logpdf(x_vad), axis=1))

    rep_ll.append(np.mean(stats.norm(holdoutmean_sample, \
                        stddv_datapoints).logpdf(holdout_gen),axis=2))
    
obs_ll_per_zi, rep_ll_per_zi = np.mean(np.array(obs_ll), axis=0), np.mean(np.array(rep_ll), axis=0)<jupyter_output><empty_output><jupyter_text>Step 2.3: Finally we compare the test statistic of the original and the replicated dataset. *   We compare the test statistics via the $p$-values.\begin{equation*}  \text{$p$-value} = p\left(t(\mathbf{X_{n,heldout}^{rep}}) < t(\mathbf{X_{n, heldout}})\right).\end{equation*}*   The **smaller** the $p$-value is, the **more different** the original dataset is from the replicated dataset.*   We **fail** the check if the $p$-value is **small**.*   Note this goes in the opposite direction to the conventional usage of $p$-values. We compute a $p$-value for each $n$ and output the average $p$-values.<jupyter_code>pvals = np.array([np.mean(rep_ll_per_zi[:,i] < obs_ll_per_zi[i]) for i in range(num_datapoints)])
holdout_subjects = np.unique(holdout_row)
overall_pval = np.mean(pvals[holdout_subjects])
print("Predictive check p-values", overall_pval)<jupyter_output>Predictive check p-values 0.6842410714285714<jupyter_text>**We passed the check!**The substitute confounder $\mathbf{z}_n$ constructed in Step 1 is valid. We are ready to move on!  An optional stepWe can also peak at **the predictive check of individual subjects**.This step is just for fun. It is how we generate Figure 2 of the paper.*   We randomly choose a subject.*   Plot the kernel density estimate of the test statistic on the replicated datasets.*   Plot the test statistic on the original dataset (the dashed vertical line).<jupyter_code>subject_no = npr.choice(holdout_subjects) 
sns.kdeplot(rep_ll_per_zi[:,subject_no]).set_title("Predictive check for subject "+str(subject_no))
plt.axvline(x=obs_ll_per_zi[subject_no], linestyle='--')<jupyter_output>/usr/local/lib/python3.6/dist-packages/statsmodels/nonparametric/kde.py:475: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.
  grid,delta = np.linspace(a,b,gridsize,retstep=True)<jupyter_text>Step 3: Correct for the substitute confounder in a causal inference. **How** to estimate causal effects?*  For simplicity, we fit a logistic regression as an outcome model here. *  The target is the observed outcome $y_n$, $n=1,\ldots, N$.*  The regressor is the multiple causes $\mathbf{X}_n$, $n=1,\ldots, N$.**How** to correct for the substitute confounder?*  We include the substitute confounder $\mathbf{Z}_n$, $n=1,\ldots, N$, into the regressors.<jupyter_code># approximate the (random variable) substitute confounders with their inferred mean.
Z_hat = z_mean_inferred 
# augment the regressors to be both the assigned causes X and the substitute confounder Z
X_aug = np.column_stack([X, Z_hat])
# holdout some data from prediction later
X_train, X_test, y_train, y_test = train_test_split(X_aug, dfy, test_size=0.2, random_state=0)
dcfX_train = sm.add_constant(X_train)
dcflogit_model = sm.Logit(y_train, dcfX_train)
dcfresult = dcflogit_model.fit_regularized(maxiter=5000)
print(dcfresult.summary())
res = pd.DataFrame({"causal_mean": dcfresult.params[:data_dim+1], \
                  "causal_std": dcfresult.bse[:data_dim+1], \
                  "causal_025": dcfresult.conf_int()[:data_dim+1,0], \
                  "causal_975": dcfresult.conf_int()[:data_dim+1,1], \
                   "causal_pval": dcfresult.pvalues[:data_dim+1]})
res["causal_sig"] = (res["causal_pval"] < 0.05)
res = res.T
res.columns = np.concatenate([["intercept"], np.array(dfX.columns)])
res = res.T
res<jupyter_output><empty_output><jupyter_text>We check the predictions to see if the logistic outcome model is a good outcome model.<jupyter_code># make predictions with the causal model 
dcfX_test = X_test
dcfy_predprob = dcfresult.predict(sm.add_constant(dcfX_test))
dcfy_pred = (dcfy_predprob > 0.5)
print(classification_report(y_test, dcfy_pred))<jupyter_output>precision    recall  f1-score   support

          0       0.96      0.91      0.93        47
          1       0.94      0.97      0.96        67

avg / total       0.95      0.95      0.95       114<jupyter_text>We are done! We have computed the average causal effect of raising the causes by one unit (see the "causal mean" column above).  Is the deconfounder worth the effort? We finally compare the **causal** estimation (with the deconfounder) with the **noncausal** estimation (with vanilla regression).  The classical logistic regression! Note it is noncausal :-(<jupyter_code># regress the outcome against the causes only (no substitute confounders)
nodcfX_train = sm.add_constant(X_train[:,:X.shape[1]])
nodcflogit_model = sm.Logit(y_train, nodcfX_train)
nodcfresult = nodcflogit_model.fit_regularized(maxiter=5000)
print(nodcfresult.summary())
res["noncausal_mean"] = np.array(nodcfresult.params)
res["noncausal_std"] = np.array(nodcfresult.bse)
res["noncausal_025"] = np.array(nodcfresult.conf_int()[:,0])
res["noncausal_975"] = np.array(nodcfresult.conf_int()[:,1])
res["noncausal_pval"] = np.array(nodcfresult.pvalues)
res["noncausal_sig"] = (res["noncausal_pval"] < 0.05)
res["diff"] = res["causal_mean"] - res["noncausal_mean"]
res["pval_diff"] = res["causal_pval"] - res["noncausal_pval"]
nodcfX_test = sm.add_constant(X_test[:,:X.shape[1]])
nodcfy_predprob = nodcfresult.predict(nodcfX_test)
nodcfy_pred = (nodcfy_predprob > 0.5)<jupyter_output><empty_output><jupyter_text>**Causal models do not hurt predictions here!**<jupyter_code>dcflogit_roc_auc = roc_auc_score(y_test, dcfy_pred)
dcffpr, dcftpr, dcfthresholds = roc_curve(y_test, dcfy_predprob)
nodcflogit_roc_auc = roc_auc_score(y_test, nodcfy_pred)
nodcffpr, nodcftpr, nodcfthresholds = roc_curve(y_test, nodcfy_predprob)
plt.figure()
plt.plot(nodcffpr, nodcftpr, label='Noncausal Logistic Regression (area = %0.9f)' % nodcflogit_roc_auc)
plt.plot(dcffpr, dcftpr, label='Causal Logistic Regression (area = %0.9f)' % dcflogit_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic')
plt.legend(loc="lower right")
plt.savefig('Log_ROC')
plt.show()<jupyter_output><empty_output><jupyter_text>**But causal models do change the regression coefficients and which features are significant.** *   The mean smoothness is a feature **significantly correlated** with the cancer diagnosis.*   But it does **not significantly** **causally affect** the cancer diagnosis.*   The effect of all features are **over-estimated** with the noncausal model except the "mean compactness".<jupyter_code>res.sort_values("pval_diff", ascending=True)[["pval_diff", "causal_pval", "noncausal_pval", "causal_sig", "noncausal_sig", "causal_mean", "noncausal_mean"]]<jupyter_output><empty_output><jupyter_text>*   We include causes into the regression **one-by-one**.*   The deconfounder coefficients **does not** flip signs.*   But classical logistic regression coefficients **does** flip signs.*   This suggests that **the deconfounder is causal**. *   It is because **causal** coefficients **do not change** as we include more variables into the system; causal estimation already controls for confounders so that it is causal.*   However, **correlation** coefficients **can change** as we include more variables into the system; if the added variable is a confounder, than the regression coefficients change to account for the confounding effects.<jupyter_code># The deconfounder with causes added one-by-one
# The first i coefficient is the causal coefficient of the first i causes.
# i = 1, ..., 8.
for i in range(X.shape[1]):
    print(i, "causes included")
    # augment the regressors to be both the assigned causes X and the substitute confounder Z
    X_aug = np.column_stack([X[:,:i], Z_hat])
    # holdout some data from prediction later
    X_train, X_test, y_train, y_test = train_test_split(X_aug, dfy, test_size=0.2, random_state=0)
    dcfX_train = sm.add_constant(X_train)
    dcflogit_model = sm.Logit(y_train, dcfX_train)
    dcfresult = dcflogit_model.fit_regularized(maxiter=5000)
    print(dcfresult.summary())    
# Logistic regression with causes added one-by-one
# The first i coefficient is the causal coefficient of the first i causes.
# i = 1, ..., 8.
for i in range(X.shape[1]):
    print(i, "causes included")
    # augment the regressors to be both the assigned causes X and the substitute confounder Z
    X_aug = np.column_stack([X[:,:i]])
    # holdout some data from prediction later
    X_train, X_test, y_train, y_test = train_test_split(X_aug, dfy, test_size=0.2, random_state=0)
    dcfX_train = sm.add_constant(X_train)
    dcflogit_model = sm.Logit(y_train, dcfX_train)
    dcfresult = dcflogit_model.fit_regularized(maxiter=5000)
    print(dcfresult.summary())<jupyter_output>0 causes included
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.6549205599225008
            Iterations: 6
            Function evaluations: 6
            Gradient evaluations: 6
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                  455
Model:                          Logit   Df Residuals:                      454
Method:                           MLE   Df Model:                            0
Date:                Wed, 26 Sep 2018   Pseudo R-squ.:              -6.217e-15
Time:                        03:50:48   Log-Likelihood:                -297.99
converged:                       True   LL-Null:                       -297.99
                                        LLR p-value:                       nan
===============================================================[...]
{'fraction_non_alphanumeric': 0.12086347545020779, 'fraction_numerical': 0.023279975373249193, 'mean_word_length': 4.053276297880615, 'pattern_counts': {'":': 7, '<': 78, '<?xml version=': 0, '>': 76, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '426403', 'license': 'MIT', 'n_tokens_mistral': 9015, 'n_tokens_neox': 8330, 'n_words': 2427, 'path': 'deconfounder_tutorial.ipynb', 'repo_name': 'tl2624/deconfounder_tutorial', 'size': '25988'}