<jupyter_start><jupyter_text>Causal Inference In Statistics - A Primer        1.5 Structural Causal Models        <NAME>        www.data4sci.com            @bgoncalves, @data4sci<jupyter_code>from collections import Counter
from pprint import pprint

import pandas as pd
import numpy as np

import matplotlib
import matplotlib.pyplot as plt 

import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf

from CausalModel import CausalModel

import watermark

%load_ext watermark
%matplotlib inline<jupyter_output><empty_output><jupyter_text>We start by print out the versions of the libraries we're using for future reference<jupyter_code>%watermark -n -v -m -g -iv<jupyter_output>numpy           1.18.1
watermark       2.0.2
statsmodels     0.11.1
matplotlib      3.1.3
autopep8        1.5
json            2.0.9
statsmodels.api 0.11.1
pandas          1.0.1
Sat Jul 25 2020 

CPython 3.7.3
IPython 6.2.1

compiler   : Clang 4.0.1 (tags/RELEASE_401/final)
system     : Darwin
release    : 19.5.0
machine    : x86_64
processor  : i386
CPU cores  : 8
interpreter: 64bit
Git hash   : 8e58bf1ca04167a82c7ab09e68e7b88b5ae50adb<jupyter_text>Load default figure style<jupyter_code>plt.style.use('./d4sci.mplstyle')
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']<jupyter_output><empty_output><jupyter_text>Modeling Causal Assumptions We can load the structure of this model directly from disk:<jupyter_code>G = CausalModel('dags/Primer.Fig.1.9.dot')<jupyter_output><empty_output><jupyter_text>And plot it quickly<jupyter_code>fig, ax = plt.subplots(1, figsize=(3, 2.2))
G.plot(ax=ax, legend=False)<jupyter_output><empty_output><jupyter_text>We generate our fake data<jupyter_code>x = np.linspace(20, 30)
y = np.linspace(100, 300, 50)
X, Y = np.meshgrid(x, y)
X = X.flatten()
Y = Y.flatten()

Z = 2*X+3*Y

model = pd.DataFrame({'X': X, 'Y': Y, 'Z': Z})<jupyter_output><empty_output><jupyter_text>And inspect it:<jupyter_code>model<jupyter_output><empty_output><jupyter_text>And to go the other way, we simply perform an OLS Linear Regression fit to estimate the dependency of Z as a function of X and Y<jupyter_code>ols = smf.ols('Z ~ X + Y', model).fit()<jupyter_output><empty_output><jupyter_text>From where we easily extract the coefficients of our model for $f_Z$<jupyter_code>ols.summary()<jupyter_output><empty_output><jupyter_text>Let us now build the model for SCM 1.5.3<jupyter_code>G2 = CausalModel()
G2.add_causation(r'$U_X$', 'X')
G2.add_causation(r'$U_Y$', 'Y')
G2.add_causation(r'$U_Z$', 'Z')
G2.add_causation('X', 'Y')
G2.add_causation('Y', 'Z')
G2.pos = {r'$U_X$': (-1, 1), r'$U_Y$':(0,1), r'$U_Z$': (1,1), 'X': (-1, 0), 'Y':(0,0), 'Z':(1, 0)}<jupyter_output><empty_output><jupyter_text>Visualize it<jupyter_code>fig, ax = plt.subplots(1, figsize=(2.2, 2.2))
G2.plot(ax=ax)<jupyter_output><empty_output><jupyter_text>And save it as a dot file so that we can reuse it later on.<jupyter_code>G2.save_model('dags/Primer.SCM.1.5.3.dot')<jupyter_output><empty_output><jupyter_text>Since this is a stochastic model, we have to use random numbers. We'll assume that the exogenous variables are normally distributed and generate 100,000 values for each of them<jupyter_code>N = 100_000
Ux = np.random.normal(size=N) 
Uy = np.random.normal(size=N)
Uz = np.random.normal(size=N)<jupyter_output><empty_output><jupyter_text>And define the 3 functions that determine the values of X, Y and Z<jupyter_code>fX = lambda ux: ux
fY = lambda x, uy: 1/3*x+uy
fZ = lambda y, uz: 1/16*y+uz<jupyter_output><empty_output><jupyter_text>With these two ingredients at hand, we can easily build up the DataFrame that describes our model<jupyter_code>X = fX(Ux)
Y = fY(X, Uy)
Z = fZ(Y, Uz)

model = pd.DataFrame({'X': X, 'Y': Y, 'Z': Z})<jupyter_output><empty_output><jupyter_text>And inspect it<jupyter_code>model<jupyter_output><empty_output><jupyter_text>If we observe that Y=3, then the expected (mean) value of Z is:<jupyter_code>Z_Y3 = np.mean(fZ(3, Uz))
print(Z_Y3)<jupyter_output>0.18488105416791895<jupyter_text>Which is close to the analytically expected value of 3/16=0.1875 Another way of calculating this value would be through a direct fit using the numerical values in our DataFrame<jupyter_code>ols = smf.ols('Z ~ Y', model).fit()
ols.summary()
ols.predict({'Y': 3})<jupyter_output><empty_output><jupyter_text>With the added advantage that this approach takes into account also the values of the exogenous (error) variables. If, on the other hand, we observe that $X=3$ then our calculation requires 2 steps, as this information must *flow* from X to Y and finally to Z<jupyter_code>Y_X3 = fY(3, Uy)
Z_X3 = np.mean(fZ(Y_X3, Uz))
print(Z_X3)<jupyter_output>0.05970934593248185<jupyter_text>Which again is close to the expected analytical value of 1/16=0.0625 Or a simple fit<jupyter_code>ols = smf.ols('Z ~ X', model).fit()
ols.summary()
ols.predict({'X': 3})<jupyter_output><empty_output><jupyter_text>We can also run the computation in the opposite direction (contrary to the direction of causality). If we want to know the expected value of $X$ given that $Y=2$. In this case we must fit a linear regression model to obtain a the dependency of X as a function of Y<jupyter_code>ols = smf.ols('X ~ Y', model).fit()
ols.summary()
ols.predict({'Y': 2})<jupyter_output><empty_output><jupyter_text>We can also use a similar approach to estimate the values of variables in the middle of a causal chain when we know the values before and after. For example, if we want to estimate the value of Y given that we observe specific values for X and Z, all we have to do is fit a linear model for Y as a function of X and Z<jupyter_code>ols = smf.ols('Y ~ X + Z', model).fit()
ols.summary()<jupyter_output><empty_output><jupyter_text>And now that we have a description of Y as a function X and Z, we can simply plugin the numbers<jupyter_code>ols.predict({'X':1, 'Z':3})<jupyter_output><empty_output><jupyter_text>Blocking variables. An interesting point to make clear (and that we will delve into more deeply in Chapt 2) is that to determine the value of a given variable, we only need the know the values of it's parents. Let's say we observe that Y=3 and are interested in finding out the value of Z. We can simply do:<jupyter_code>ols = smf.ols('Z ~ Y', model).fit()
ols.predict({'Y': 3})<jupyter_output><empty_output><jupyter_text>On the other hand, let's say that we also observe that X=1. In this case, we would do:<jupyter_code>ols = smf.ols('Z ~ Y + X', model).fit()
ols.summary()
ols.predict({'Y': 3, 'X':1})<jupyter_output><empty_output><jupyter_text>which returns the same result (up to numerical precision). While counterintuitive, this result is simple to understand. In our causal model, the value of Z is determined only by the value of Y. Once Y is determined, the factors that might have led to its value (a statistical fluctuation in Uy, or the value of X) are irrelevant. All that Z "sees" is the current value of Y and Y *shields* Z from everything that happened before in its causal chain In this simple observation lies much of the power of this approach!  Product Decomposition<jupyter_code>G = CausalModel('dags/Primer.Fig.1.10.dot')
fig, ax = plt.subplots(1, figsize=(2.2, 2.2))
G.plot(ax=ax)
PZ = pd.DataFrame({'Z': [0, 1], 'P(Z)': ['1-r', 'r']})
PZ
MedicationProb = pd.DataFrame({
    'X': [1, 1, 0, 0],
    'Z': [0, 1, 0, 1],
    ' P(X | Z)': [r'$q_{0}$', r'$q_{1}$', r'$1-q_0$', r'$1-q_{1}$']
})
MedicationProb
DeathProb = pd.DataFrame({
    'Y':[1, 1, 1, 1, 0, 0, 0, 0],
    'X':[1, 0, 0, 1, 1, 0, 0, 1],
    'Z':[0, 0, 1, 1, 0, 0, 1, 1],
    'P(Y | X,Z)':[r'$p_{10}$', r'$p_{00}$', r'$p_{01}$', r'$p_{11}$', r'$1-p_{10}$', r'$1-p_{00}$', r'$1-p_{01}$', r'$1-p_{11}$']})
DeathProb<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13297255508310785, 'fraction_numerical': 0.038654812524159254, 'mean_word_length': 5.2445695897023334, 'pattern_counts': {'":': 0, '<': 103, '<?xml version=': 0, '>': 103, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '322409', 'license': 'MIT', 'n_tokens_mistral': 3105, 'n_tokens_neox': 2913, 'n_words': 881, 'path': '1.5 - Structural Causal Models.ipynb', 'repo_name': 'm5l14i11/Causality', 'size': '7770'}