Printing the structural causal equations used for generating a sample

[Stan201]

Hi all,

For a simulation study I have sampled data from a SCM by first using gcm.auto.assign_causal_machanisms(causal_model, data) and also used causal_model.set_causal_mechanism(node, gcm.AdditiveNoiseModel(gcm.ml.create_linear_regressor()). For my thesis, I want to print the linear equation which is used to generate data from the causal model, but cant find how I can do that. I want to know it in the form: Node1 = 0.8*Node2 + 0.05Node3 + N(0,20) for example. Is this possible for linear equations? Is this also possible for non-linear equations?

Thanks in advance!
Stan

[bloebp]

Hi Stan,

there is no direct way to print it, but you can access the model parameters. These, however, depend on the used models. Generally, you can access the fitted sklearn models and their parameters. For non-linear models, these might be difficult to interpret (depending on the model type). The easiest is with the linear models as you did by setting them explicitly. Note, however, that the AdditiveNoiseModel (and auto assignment) typically uses an EmpiricalDistribution for the noise which is parameter free. But you can also explicitly assign, e.g., a Gaussian distribution. I prepared a code snippet that prints something like your example, but it currently only works with linear models:

import networkx as nx
import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import LinearRegression

import dowhy.gcm as gcm
from dowhy.gcm import EmpiricalDistribution, ScipyDistribution
from dowhy.gcm.causal_mechanisms import StochasticModel
from dowhy.graph import is_root_node, get_ordered_predecessors

X = np.random.normal(loc=0, scale=1, size=1000)
Y = 2 * X + np.random.normal(loc=0, scale=1, size=1000)
Z = 3 * Y + np.random.normal(loc=0, scale=1, size=1000)
data = pd.DataFrame(data=dict(X=X, Y=Y, Z=Z))

causal_model = gcm.StructuralCausalModel(nx.DiGraph([("X", "Y"), ("Y", "Z")]))
causal_model.set_causal_mechanism("X", EmpiricalDistribution())
causal_model.set_causal_mechanism("Y", gcm.AdditiveNoiseModel(gcm.ml.create_linear_regressor(),
                                                              ScipyDistribution(stats.norm)))
causal_model.set_causal_mechanism("Z", gcm.AdditiveNoiseModel(gcm.ml.create_linear_regressor()))
gcm.fit(causal_model, data)

print(causal_model.causal_mechanism("Y").prediction_model.sklearn_model.coef_)
print(causal_model.causal_mechanism("Y").noise_model.parameters)
print(causal_model.causal_mechanism("Z").prediction_model.sklearn_model.coef_)

print("Put together:")


def get_distribution_string(distribution_mechanism: StochasticModel):
    if isinstance(distribution_mechanism, ScipyDistribution):
        return f"{str(distribution_mechanism)} with parameters {distribution_mechanism.parameters}"
    else:
        return str(distribution_mechanism)


for n in causal_model.graph.nodes:
    mechanism = causal_model.causal_mechanism(n)
    if is_root_node(causal_model.graph, n):
        print(f"{n} ~ {get_distribution_string(mechanism)}")
    else:
        if isinstance(mechanism, gcm.AdditiveNoiseModel):
            if isinstance(mechanism.prediction_model.sklearn_model, LinearRegression):
                coefficients = mechanism.prediction_model.sklearn_model.coef_
                ordered_parents = get_ordered_predecessors(causal_model.graph, n)

                equation = n + " :="
                for i, parent in enumerate(ordered_parents):
                    if i > 0:
                        if coefficients[i] >= 0:
                            equation += " +"
                        else:
                            equation += " -"

                    coeff = coefficients[i] if coefficients[i] >= 0 else abs(coefficients[i])
                    equation += f" {coeff} * {parent}"

                    equation += f" + Noise[{get_distribution_string(mechanism.noise_model)}]"
                    print(equation)
            else:
                raise NotImplementedError("Not implemented")
        else:
            raise NotImplementedError("Not implemented")

The example should run by itself and also with your causal model if all models are linear.

[Stan201]

Sweet this is exactly what I was looking for!