Using DoWhy for categorical network causal effect estimation

[NianzuMa]

Question
I want to use DoWhy to learn parameter for a predefined network structure B -> A, A -> T, B -> T given a simulated dataset.
After this network is learned with parameter. Then I want to do manipulation of variable A (do operation, hard intervention) and see how the posterior of T change.

I found that this is very straightforward in pgmpy and SMILE, but I am not sure how to do it in DoWhy.

My sample code is as follows:

# %%
import pandas as pd
import numpy as np
import os
import scipy
import math
from functools import partial
import dowhy.gcm as gcm
import networkx as nx


# %%
output_folder = "./simulated_data"
os.makedirs(output_folder, exist_ok=True)

simulated_dataset_path = "./simulated_data/simulated_discrete_dataset.csv"

# %%
def simulate_bayesian_network():
    # Set random seed for reproducibility
    np.random.seed(1)

    N = 50000

    # Step 1: B ~ Bernoulli(0.5)
    B = np.random.binomial(1, 0.5, size=(N,))

    # Step 2: A depends on B
    A = np.zeros(N, dtype=int)
    A[B == 1] = np.random.binomial(1, 0.9, size=(np.sum(B == 1),))
    A[B == 0] = np.random.binomial(1, 0.1, size=(np.sum(B == 0),))

    # Step 3: T depends on A and B
    T = np.zeros(N, dtype=int)
    T[(A == 1) & (B == 1)] = np.random.binomial(1, 0.9, size=np.sum((A == 1) & (B == 1)))
    T[(A == 1) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 1) & (B == 0)))
    T[(A == 0) & (B == 1)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 1)))
    T[(A == 0) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 0)))
    
    data = {}
    data['A'] = A
    data['B'] = B
    data['T'] = T

    return pd.DataFrame(data)

# Simulate and preview the dataset
data = simulate_bayesian_network()
print(data.head())

# %%
# gcm stands for "Graphical Causal Models"

causal_model = gcm.StructuralCausalModel(nx.DiGraph([('B', 'A'), ('B', 'T'), ('A', 'T')]))


# %%

auto_assignment_summary = gcm.auto.assign_causal_mechanisms(causal_model, data)


# %%
print(auto_assignment_summary)

# %% [markdown]
# ## Comments
# 
# Auto detected causal mechanism is not right. 
# 
# This is the code for all possible causal mechanism: 
# https://github.com/py-why/dowhy/blob/main/dowhy/gcm/causal_mechanisms.py
# 
# 

# %%
# causal_model.set_causal_mechanism('B', gcm.EmpiricalDistribution())
# causal_model.set_causal_mechanism('A', gcm.classification.LogisticRegressionClassifier())
# causal_model.set_causal_mechanism('T', gcm.classification.LogisticRegressionClassifier())

# %%

gcm.fit(causal_model, data)

# %%
gcm.average_causal_effect(causal_model,
                         'T',
                         interventions_alternative={'A': lambda x: 1},
                         interventions_reference={'A': lambda x: 0},
                         num_samples_to_draw=100000)

# %%
def print_out_posterior(df, variable_n):
    posterior_counts = df[variable_n].value_counts(normalize=True)
    print(posterior_counts)
    

# %%
samples = gcm.interventional_samples(causal_model,
                                     {'A': lambda x: 1},
                                     num_samples_to_draw=1000)
samples.head()

print_out_posterior(samples, "T")

# %%
mean_T = samples['T'].mean()
print(f"The mean of column 'T' is: {mean_T}")

# %%
samples_2 = gcm.interventional_samples(causal_model,
                                     {'A': lambda x: 0},
                                     num_samples_to_draw=1000)
samples_2.head()

# %%
mean_T_2 = samples_2['T'].mean()
print(f"The mean of column 'T' is: {mean_T_2}")

# %%
print_out_posterior(samples_2, "T")

In pgmpy and SMILE, there is very explicit way to defined a categorical network, pass network structure to it. Then the packages will learn the conditional probability table from the data for the network.

(1)
However, in DoWhy I found that auto_assignment_summary = gcm.auto.assign_causal_mechanisms(causal_model, data) is not work as expected and also I cannot find a good way using causal_model.set_causal_mechanism to define a categorical network.

In the auto_assignment_summary, Node A is assigned "Discrete AdditiveNoiseModel using LinearRegression", Node T is assigned "Discrete AdditiveNoiseModel using Pipeline".

Could you please help me to explicitly define a categorical network?

I want explore DoWhy because my real world project contains (1) categorical variable depend on continuous variable, (2) continuous variable depend on categorical variable, (3) continuous variable depend on continuous variable. (4) categorical variable depend on categorical variable. DoWhy looks like can deal with all of these cases (other package has limitations) and can do hard intervention of a variable and see the effect estimation of the target variable.

(2)
Is it OK to use gcm.interventional_samples to do hard intervention and sample data, then estimate posterior probability of target variable T? In my sample code, am I doing in the right way?
In current setting, I found that

T 
0 0.482 
1 0.479 
2 0.021 
-1 0.018

It is surprising that T contains value 2 and -1. But in my input data, it should only have two value 0 and 1, which is a categorical variable.
Why it behaves in this way and what I can do to make it right?

Expected behavior
The posterior of T should only contains value 0 and 1.

Version information:

• DoWhy version [0.13]

[bloebp]

Hi, thanks for detailed question!

(1)
Based on the data generation process and the information that the auto assignment chose "Discrete Additive Noise Model using Linear Regression" indicates that the dataframe is numerical, i.e., the auto assignment assumes discrete ordinal data in that case (not categorical). If you want to explicitly model this as categorical data and not as discrete data, then the best is to explicitly convert it into strings. For instance, return pd.DataFrame(data).astype(str).

Or you just give the categories specific string names.

(2)

Is it OK to use gcm.interventional_samples to do hard intervention and sample data, then estimate posterior probability of target variable T?

Yes that would be possible.

Why it behaves in this way and what I can do to make it right?

This is most likely because the "Discrete Additive Noise Model" is, by default, not setting a limit on the support. That is, if the noise of the data contains 0 and 1, then what can happen is +/-1 is added to "T", which leads to a range of [-1, 2]. That being said, I believe once these are correctly modeled as categorical (i.e. using strings instead of numbers), it should be binary again.

[NianzuMa]

Thanks for your suggestion.

I have modified according to your suggestion as follows:

# %%
import pandas as pd
import numpy as np
import os
import scipy
import math
from functools import partial
import dowhy.gcm as gcm
import networkx as nx


# %%
import dowhy
print(dowhy.__version__)

# %%
output_folder = "./simulated_data"
os.makedirs(output_folder, exist_ok=True)

simulated_dataset_path = "./simulated_data/simulated_discrete_dataset.csv"

# %%
def simulate_bayesian_network():
    # Set random seed for reproducibility
    np.random.seed(1)

    N = 50000

    # Step 1: B ~ Bernoulli(0.5)
    B = np.random.binomial(1, 0.5, size=(N,))

    # Step 2: A depends on B
    A = np.zeros(N, dtype=int)
    A[B == 1] = np.random.binomial(1, 0.9, size=(np.sum(B == 1),))
    A[B == 0] = np.random.binomial(1, 0.1, size=(np.sum(B == 0),))

    # Step 3: T depends on A and B
    T = np.zeros(N, dtype=int)
    T[(A == 1) & (B == 1)] = np.random.binomial(1, 0.9, size=np.sum((A == 1) & (B == 1)))
    T[(A == 1) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 1) & (B == 0)))
    T[(A == 0) & (B == 1)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 1)))
    T[(A == 0) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 0)))
    
    data = {}
    data['A'] = A
    data['B'] = B
    data['T'] = T

    #return pd.DataFrame(data)
    return pd.DataFrame(data).astype(str)

# Simulate and preview the dataset
data = simulate_bayesian_network()
print(data.head())

# %%
# gcm stands for "Graphical Causal Models"

causal_model = gcm.StructuralCausalModel(nx.DiGraph([('B', 'A'), ('B', 'T'), ('A', 'T')]))


# %%

auto_assignment_summary = gcm.auto.assign_causal_mechanisms(causal_model, data)


# %%
print(auto_assignment_summary)

# %% [markdown]
# ## Comments
# 
# Auto detected causal mechanism is not right. 
# 
# This is the code for all possible causal mechanism: 
# https://github.com/py-why/dowhy/blob/main/dowhy/gcm/causal_mechanisms.py
# 
# 

# %%
# causal_model.set_causal_mechanism('B', gcm.EmpiricalDistribution())
# causal_model.set_causal_mechanism('A', gcm.classification.LogisticRegressionClassifier())
# causal_model.set_causal_mechanism('T', gcm.classification.LogisticRegressionClassifier())

# %%

gcm.fit(causal_model, data)

# %%
gcm.average_causal_effect(causal_model,
                         'T',
                         interventions_alternative={'A': lambda x: "1"},
                         interventions_reference={'A': lambda x: "0"},
                         num_samples_to_draw=100000)

# %%
def print_out_posterior(df, variable_n):
    posterior_counts = df[variable_n].value_counts(normalize=True)
    print(posterior_counts)
    

# %%
samples = gcm.interventional_samples(causal_model,
                                     {'A': lambda x: "1"},
                                     num_samples_to_draw=10000)
samples.head()

print_out_posterior(samples, "T")

# %%
samples_2 = gcm.interventional_samples(causal_model,
                                     {'A': lambda x: "0"},
                                     num_samples_to_draw=10000)
samples_2.head()

# %%
print_out_posterior(samples_2, "T")

The auto_assignment_summary now shows:

--- Node: B
Node B is a root node. Therefore, assigning 'Empirical Distribution' to the node representing the marginal distribution.

--- Node: A
Node A is a non-root node with categorical data. Assigning 'Classifier FCM based on LogisticRegression(max_iter=10000)' to the node.
This represents the causal relationship as A := f(B,N).
For the model selection, the following models were evaluated on the (negative) F1 metric:
LogisticRegression(max_iter=10000): -0.8985975589105193
HistGradientBoostingClassifier: -0.8985975589105193
Pipeline(steps=[('polynomialfeatures',
                 PolynomialFeatures(degree=3, include_bias=False)),
                ('logisticregression', LogisticRegression(max_iter=10000))]): -0.8985975589105193

--- Node: T
Node T is a non-root node with categorical data. Assigning 'Classifier FCM based on LogisticRegression(max_iter=10000)' to the node.
This represents the causal relationship as T := f(A,B,N).
For the model selection, the following models were evaluated on the (negative) F1 metric:
LogisticRegression(max_iter=10000): -0.8997770692506857
HistGradientBoostingClassifier: -0.8997770692506857
Pipeline(steps=[('polynomialfeatures',
                 PolynomialFeatures(degree=3, include_bias=False)),
                ('logisticregression', LogisticRegression(max_iter=10000))]): -0.8997770692506857

It looks reasonable now.

Question A
My questions is, it is not like pgmpy and SMILE that deal with everything using conditional probability table for all variables. Is there anyway that DoWhy doing similar thing like SMILE? Why DoWhy use Logistic Regression for categorical network?

To give a clear impression, the learned network file of SMILE looks like below:

<nodes>
		<cpt id="B">
			<state id="B___0" />
			<state id="B___1" />
			<probabilities>0.4988 0.5012000000000001</probabilities>
		</cpt>
		<cpt id="A">
			<state id="A___0" />
			<state id="A___1" />
			<parents>B</parents>
			<probabilities>0.8965918203688853 0.1034081796311147 0.09940143655227453 0.9005985634477254</probabilities>
		</cpt>
		<cpt id="T">
			<state id="T___0" />
			<state id="T___1" />
			<parents>B A</parents>
			<probabilities>0.9036268503197531 0.09637314968024685 0.8941450174486235 0.1058549825513765 0.8988358089120836 0.1011641910879165 0.1015109220612344 0.8984890779387655</probabilities>
		</cpt>
	</nodes>

I found that the printed posterior of T
(1) when do(A="1")
T
1 0.6353
0 0.3647

(2) when do(A="0")
T
0 0.7696
1 0.2304

The ground truth is can be generated by this python script:

import numpy as np

# Example 1: With confounding
print('example 1: with confounding')
np.random.seed(1)
N = 50000
B = np.random.binomial(1, 0.5, size=(N,))

# Values of A depend on B
A = np.zeros(N, dtype=int)
A[B == 1] = np.random.binomial(1, 0.9, size=np.sum(B == 1))
A[B == 0] = np.random.binomial(1, 0.1, size=np.sum(B == 0))

# Values of T depend on A and B
T = np.zeros(N, dtype=int)
T[(A == 1) & (B == 1)] = np.random.binomial(1, 0.9, size=np.sum((A == 1) & (B == 1)))
T[(A == 1) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 1) & (B == 0)))
T[(A == 0) & (B == 1)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 1)))
T[(A == 0) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 0)))

# 1.1 Observing A=1
obs_idx = (A == 1)
print(f'average value of B when observing A==1: {np.mean(B[obs_idx])}')
print(f'average value of T when observing A==1: {np.mean(T[obs_idx])}')

# 1.2 Manipulating A=1 (do(A=1))
A = np.ones(N, dtype=int)  # intervene A
T = np.zeros(N, dtype=int)
T[(A == 1) & (B == 1)] = np.random.binomial(1, 0.9, size=np.sum((A == 1) & (B == 1)))
T[(A == 1) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 1) & (B == 0)))
T[(A == 0) & (B == 1)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 1)))
T[(A == 0) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 0)))

print("\n1.2 Manipulating A=1 (do(A=1))")
print(f'average value of B when manipulating A==1: {np.mean(B)}')
print(f'average value of T when manipulating A==1: {np.mean(T)}')
print('the difference in T when observing vs. manipulating A is a result of confounding by B.')
print('The Value of B when observing vs. manipulating A is different, when B is a cause of A')
print('This difference propagates to the difference in T, since B also causes T.')

# 1.3 Manipulating A=0 (do(A=0))
A = np.zeros(N, dtype=int)  # intervene A
T = np.zeros(N, dtype=int)
T[(A == 1) & (B == 1)] = np.random.binomial(1, 0.9, size=np.sum((A == 1) & (B == 1)))
T[(A == 1) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 1) & (B == 0)))
T[(A == 0) & (B == 1)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 1)))
T[(A == 0) & (B == 0)] = np.random.binomial(1, 0.1, size=np.sum((A == 0) & (B == 0)))

print("\n1.3 Manipulating A=0 (do(A=0))")
print(f'average value of B when manipulating A==0: {np.mean(B)}')
print(f'average value of T when manipulating A==0: {np.mean(T)}')
print('the difference in T when observing vs. manipulating A is a result of confounding by B.')
print('The Value of B when observing vs. manipulating A is different, when B is a cause of A')
print('This difference propagates to the difference in T, since B also causes T.')

The ground truth should be
do(A="1")
P(T=1) = 0.49728

do(A="0")
P(T=1) = 0.09906

I used SMILE package, it gives
(1)do(A="1")
P(T=T___0)=0.49687680884046415
P(T=T___1)=0.503123191159536

(2)do(A="0")
P(T=T___0)=0.901225580366229
P(T=T___1)=0.09877441963377087

Is result from SMILE is very similar to ground truth. However, DoWhy gives very different result.

Question B
Am I missing something?
Thank you very much for your answer.

[bloebp]

Thanks for flagging this. It is indeed surprising that the logistic regression model is so badly calibrated.

I took a closer look and the issue is with the model selection criteria that only focuses on the F1 score, but that might not necessarily ensure a 'well' calibrated model. I opened a PR to address this (#1349).

For now, you can override the selected model manually. For instance, replace

auto_assignment_summary = gcm.auto.assign_causal_mechanisms(causal_model, data)

with

causal_model.set_causal_mechanism('B', gcm.EmpiricalDistribution())
causal_model.set_causal_mechanism('A', gcm.ClassifierFCM(gcm.ml.create_hist_gradient_boost_classifier()))
causal_model.set_causal_mechanism('T', gcm.ClassifierFCM(gcm.ml.create_hist_gradient_boost_classifier()))

or if you want a model closer to mimic a conditional table, I guess you can also modify the parameters of a scikit-learn DecisionTreeClassifier:

causal_model.set_causal_mechanism('B', gcm.EmpiricalDistribution())
causal_model.set_causal_mechanism('A', gcm.ClassifierFCM(gcm.ml.SklearnClassificationModel(DecisionTreeClassifier())))
causal_model.set_causal_mechanism('T', gcm.ClassifierFCM(gcm.ml.SklearnClassificationModel(DecisionTreeClassifier())))

The reason why we use these models is to have more flexibility for the input data (it supports an arbitrary mix of continuous, discrete and categorical input variables).

[github-actions]

This issue is stale because it has been open for 30 days with no activity.

[github-actions]

This issue was closed because it has been inactive for 14 days since being marked as stale.