<jupyter_start><jupyter_text>Causal Discovery in PythonOver the last decade, causal inference gained a lot of traction in academia and in the industry. Causal models can be immensely helpful in various areas – from marketing to medicine and from finance to cybersecurity. To make these models work, we need not only data as in traditional machine learning, but also a causal structure. Traditional way to obtain the latter is through well-designed experiments. Unfortunately, experiments can be tricky – difficult to design, expensive or unethical. Causal discovery (also known as structure learning) is an umbrella term that describes several families of methods aiming at discovering causal structure from observational data. During the talk, we will review the basics of causal inference and introduce the concept of causal discovery. Next, we will discuss differences between various approaches to causal discovery. Finally, we will see a series of practical examples of causal discovery using Python.  Installing the environment* Using **Conda**:    `conda env create --file econml-dowhy-py38.yml`* Installing `gcastle` only:    `pip install gcastle==1.0.3rc3`<jupyter_code>def get_n_undirected(g):
    
    total = 0
    
    for i in range(g.shape[0]):
        for j in range(g.shape[0]):
            if (g[i, j] == 1) and (g[i, j] == g[j, i]):
                total += .5
    
    return total<jupyter_output><empty_output><jupyter_text>PC algorithm**PC algorithm** starts with a **fully connected** graph and then performs a series of steps to remove edges, based on graph independence structure. Finally, it tries to orient as many edges as possible.Figure 1 presents a visual representatrion of these steps.Figure 1. Original graph and PC algorithm steps. (Gylmour et al., 2019)Interested in more details? [Gylmour et al. - Review of Causal Discovery Methods Based on Graphical Models (2019)](https://www.frontiersin.org/articles/10.3389/fgene.2019.00524/full)<jupyter_code># Let's implement this structure
x = np.random.randn(1000)
y = np.random.randn(1000)

z = x + y + .1 * np.random.randn(1000)
w = .7 * z + .1 * np.random.randn(1000)
# To matrix
pc_dataset = np.vstack([x, y, z, w]).T
# Sanity check
pc_dataset, pc_dataset.shape
# Build the model
pc = PC()
pc.learn(pc_dataset)
pc.causal_matrix
# Get learned graph
learned_graph = nx.DiGraph(pc.causal_matrix)

# Relabel the nodes
MAPPING = {k: v for k, v in zip(range(4), ['X', 'Y', 'Z', 'W'])}
learned_graph = nx.relabel_nodes(learned_graph, MAPPING, copy=True)

# Plot the graph
nx.draw(
    learned_graph, 
    with_labels=True,
    node_size=1800,
    font_size=18,
    font_color='white'
)<jupyter_output><empty_output><jupyter_text>Let's do some more discovery!  Generate datasetsWe'll use a [scale-free](https://en.wikipedia.org/wiki/Scale-free_network) model to generate graphs.Then we'll use three different causal models on this graph:* linear Gaussian* linear exp* non-linear quadratic<jupyter_code># Data simulation, simulate true causal dag and train_data.
true_dag = DAG.scale_free(n_nodes=10, n_edges=15, seed=18)


DATA_PARAMS = {
    'linearity': ['linear', 'nonlinear'], 
    'distribution': {
        'linear': ['gauss', 'exp'],
        'nonlinear': ['quadratic']
    }
}

datasets = {}

for linearity in DATA_PARAMS['linearity']:
    for distr in DATA_PARAMS['distribution'][linearity]:
        
        datasets[f'{linearity}_{distr}'] = IIDSimulation(
            W=true_dag, 
            n=2000, 
            method=linearity, 
            sem_type=distr)

# Sanity check
datasets
plt.figure(figsize=(16, 8))
for i, dataset in enumerate(datasets):
    X = datasets[dataset].X
    
    plt.subplot(4, 2, i + 1)
    plt.hist(X[:, 0], bins=100)
    plt.title(dataset)
    plt.axis('off')
    
    plt.subplot(4, 2, i + 5)
    plt.scatter(X[:, 8], X[:, 4], alpha=.3)
    plt.title(dataset)
    plt.axis('off')
    
plt.subplots_adjust(hspace=.7)
plt.show()<jupyter_output><empty_output><jupyter_text>Visualize the true graph<jupyter_code>nx.draw(
    nx.DiGraph(true_dag), 
    node_size=1800,
    alpha=.7,
    pos=nx.circular_layout(nx.DiGraph(true_dag))
)
GraphDAG(true_dag)
plt.show()<jupyter_output><empty_output><jupyter_text>Method comparison<jupyter_code>methods = OrderedDict({
    'PC': PC,
    'GES': ges,
    'LiNGAM': ICALiNGAM,
    'GOLEM': GOLEM
})
%%time

results = {}

for k, dataset in datasets.items():
    print(f'************* Current dataset: {k}\n')
    X = dataset.X
    
    results[dataset] = {}
    
    for method in methods:
        
        if method not in ['GES', 'CORL']:
            print(f'Method: {method}')
            
            # Fit the model
            if method == 'GOLEM':
                model = methods[method](num_iter=2.5e4)
            else:
                model = methods[method]()
            
            model.learn(X)
            
            pred_dag = model.causal_matrix

        elif method == 'GES':
            print(f'Method: {method}')
            
            # Fit the model
            pred_dag, _ = methods[method].fit_bic(X)
              
        # Get n undir edges
        n_undir = get_n_undirected(pred_dag)

        # Plot results
        GraphDAG(pred_dag, true_dag, 'result')

        mt = MetricsDAG(pred_dag, true_dag)
        print(f'FDR: {mt.metrics["fdr"]}')
        print(f'Recall: {mt.metrics["recall"]}')
        print(f'Precision: {mt.metrics["precision"]}')
        print(f'F1 score: {mt.metrics["F1"]}')
        print(f'No. of undir. edges: {n_undir}\n')
        print('-' * 50, '\n')

        results[dataset][method] = pred_dag      
            
    print('\n')<jupyter_output>************* Current dataset: linear_gauss

Method: PC
{'fraction_non_alphanumeric': 0.12206654991243432, 'fraction_numerical': 0.01768826619964974, 'mean_word_length': 3.1595047341587765, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 2, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '17588', 'license': 'MIT', 'n_tokens_mistral': 1952, 'n_tokens_neox': 1811, 'n_words': 534, 'path': '2022-03--pydata-hamburg/02 - Causal discovery.ipynb', 'repo_name': 'AlxndrMlk/talks-and-workshops', 'size': '5710'}