{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Structure Learning in Bayesian Networks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we show a few examples of Causal Discovery or Structure Learning in pgmpy. pgmpy currently has the following algorithm for causal discovery:\n", "\n", "1. **PC**: Has $3$ variants original, stable, and parallel. PC is a constraint-based algorithm that utilizes Conditional Independence tests to construct the model.\n", "2. **Hill-Climb Search**: Hill-Climb Search is a greedy optimization-based algorithm that makes iterative local changes to the model structure such that it improves the overall score of the model.\n", "3. **Greedy Equivalence Search (GES)**: Another score-based method that makes greedy modifications to the model to improve its score iteratively.\n", "4. **ExpertInLoop**: An iterative algorithm that combines Conditional Independence testing with expert knowledge. The user or an LLM can act as the expert.\n", "5. **Exhaustive Search**: Exhaustive search iterates over all possible network structures on the given variables to find the most optimal one. As it tries to enumerate all possible network structures, it is intractable when the number of variables in the data is large.\n", "\n", "The following Conditional Independence Tests are available to use with PC algorithm.\n", "1. **Discrete Data**: When all variables are discrete/categorical.\n", " 1. **Chi-square test**: `ci_test=\"chi_square\"`\n", " 2. **G-squared**: `ci_test=\"g_sq\"`\n", " 3. **Log-likelihood**: Is equivalent to G-squared test. `ci_test=\"log_likelihood`\n", "2. **Continuous Data**: When all variables are continuous/numerical.\n", " 1. **Partial Correlation**: `ci_test=\"pearsonr\"`\n", "3. **Mixed Data**: When there is a mix of categorical and continuous variables.\n", " 1. **Pillai**: `ci_test=\"pillai\"`\n", "\n", "For Hill-Climb, Exhausitive Search, and GES the following scoring methods can be used:\n", "1. **Discrete Data**: When all variables are discrete/categorical. \n", " 1. **BIC Score**: `scoring_method=\"bic-d\"`\n", " 2. **AIC Score**: `scoring_method=\"aic-d\"`\n", " 3. **K2 Score**: `scoring_method=\"k2\"`\n", " 4. **BDeU Score**: `scoring_method=\"bdeu\"`\n", " 5. **BDs Score**: `scoring_method=\"bds\"`\n", "2. **Continuous Data**: When all variables are continuous/numerical.\n", " 1. **Log-Likelihood**: `scoring_method=\"ll-g\"`\n", " 2. **AIC**: `scoring_method=\"aic-g\"`\n", " 3. **BIC**: `scoring_method=\"bic-g\"`\n", "3. **Mixed Data**: When there is a mix of discrete and continuous variables.\n", " 1. **AIC**: `scoring_method=\"aic-cg\"`\n", " 2. **BIC**: `scoring_method=\"bic-cg\"`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 0. Simulate some sample datasets" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from itertools import combinations\n", "\n", "import networkx as nx\n", "import numpy as np\n", "from sklearn.metrics import f1_score\n", "\n", "from pgmpy.estimators import PC, HillClimbSearch, GES\n", "from pgmpy.utils import get_example_model\n", "from pgmpy.metrics import SHD" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "aa35b1a7897646d7ac141e27e04a6165", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/37 [00:00\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
b1191cspGeutGfixCcspAyecOyedEsucAcchByceP...dnaKfolKycgXlacZnuoMdnaGb1583mopByaeMftsJ
01.0606412.0444770.3412161.448399-0.3517162.189750-1.689554-0.2284562.871002-0.433597...1.5228171.6459831.5955022.465247-0.5329871.126289-0.3025890.7734833.8848570.729557
10.6321510.9643210.8302290.6965980.6392040.058108-0.7361890.7120951.4674980.320727...1.2226021.7907271.7635902.945772-2.5324641.4606992.7325950.0979822.5660640.652853
20.5857662.8624370.9222910.3700000.7239322.487161-1.916624-0.3003592.050980-0.064301...0.5993051.3020910.5097173.090268-1.7456130.1680430.851346-0.6404725.8007120.031888
31.8028662.2770380.6085592.1802830.1164532.539035-1.656839-1.4205402.605192-0.160302...2.0156632.8235882.1016253.521299-1.2123911.4853692.4492962.0321163.8889173.036650
41.8685482.4809991.0793642.4138620.9617432.280195-1.740610-0.4726923.338063-0.509144...1.1421391.251832-0.1457892.057217-2.862230-0.0897211.0208320.0405894.9258480.575977
\n", "

5 rows × 46 columns

\n", "" ], "text/plain": [ " b1191 cspG eutG fixC cspA yecO yedE \\\n", "0 1.060641 2.044477 0.341216 1.448399 -0.351716 2.189750 -1.689554 \n", "1 0.632151 0.964321 0.830229 0.696598 0.639204 0.058108 -0.736189 \n", "2 0.585766 2.862437 0.922291 0.370000 0.723932 2.487161 -1.916624 \n", "3 1.802866 2.277038 0.608559 2.180283 0.116453 2.539035 -1.656839 \n", "4 1.868548 2.480999 1.079364 2.413862 0.961743 2.280195 -1.740610 \n", "\n", " sucA cchB yceP ... dnaK folK ycgX lacZ \\\n", "0 -0.228456 2.871002 -0.433597 ... 1.522817 1.645983 1.595502 2.465247 \n", "1 0.712095 1.467498 0.320727 ... 1.222602 1.790727 1.763590 2.945772 \n", "2 -0.300359 2.050980 -0.064301 ... 0.599305 1.302091 0.509717 3.090268 \n", "3 -1.420540 2.605192 -0.160302 ... 2.015663 2.823588 2.101625 3.521299 \n", "4 -0.472692 3.338063 -0.509144 ... 1.142139 1.251832 -0.145789 2.057217 \n", "\n", " nuoM dnaG b1583 mopB yaeM ftsJ \n", "0 -0.532987 1.126289 -0.302589 0.773483 3.884857 0.729557 \n", "1 -2.532464 1.460699 2.732595 0.097982 2.566064 0.652853 \n", "2 -1.745613 0.168043 0.851346 -0.640472 5.800712 0.031888 \n", "3 -1.212391 1.485369 2.449296 2.032116 3.888917 3.036650 \n", "4 -2.862230 -0.089721 1.020832 0.040589 4.925848 0.575977 \n", "\n", "[5 rows x 46 columns]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Discrete variable dataset\n", "alarm_model = get_example_model(\"alarm\")\n", "alarm_samples = alarm_model.simulate(int(1e3))\n", "alarm_samples.head()\n", "\n", "# Continuous variable dataset\n", "ecoli_model = get_example_model(\"ecoli70\")\n", "ecoli_samples = ecoli_model.simulate(int(1e3))\n", "ecoli_samples.head()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Function to evaluate the learned model structures.\n", "def get_f1_score(estimated_model, true_model):\n", " nodes = estimated_model.nodes()\n", " est_adj = nx.to_numpy_array(\n", " estimated_model.to_undirected(), nodelist=nodes, weight=None\n", " )\n", " true_adj = nx.to_numpy_array(\n", " true_model.to_undirected(), nodelist=nodes, weight=None\n", " )\n", "\n", " f1 = f1_score(np.ravel(true_adj), np.ravel(est_adj))\n", " print(\"F1-score for the model skeleton: \", f1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. PC algorithm" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:pgmpy: Datatype (N=numerical, C=Categorical Unordered, O=Categorical Ordered) inferred from data: \n", " {'INTUBATION': 'C', 'ANAPHYLAXIS': 'C', 'VENTLUNG': 'C', 'HISTORY': 'C', 'LVEDVOLUME': 'C', 'VENTMACH': 'C', 'TPR': 'C', 'MINVOL': 'C', 'VENTTUBE': 'C', 'HRBP': 'C', 'BP': 'C', 'PAP': 'C', 'PULMEMBOLUS': 'C', 'DISCONNECT': 'C', 'ERRLOWOUTPUT': 'C', 'HYPOVOLEMIA': 'C', 'CVP': 'C', 'HREKG': 'C', 'CO': 'C', 'PRESS': 'C', 'HRSAT': 'C', 'CATECHOL': 'C', 'EXPCO2': 'C', 'PVSAT': 'C', 'ARTCO2': 'C', 'LVFAILURE': 'C', 'KINKEDTUBE': 'C', 'VENTALV': 'C', 'FIO2': 'C', 'SHUNT': 'C', 'PCWP': 'C', 'SAO2': 'C', 'HR': 'C', 'MINVOLSET': 'C', 'STROKEVOLUME': 'C', 'INSUFFANESTH': 'C', 'ERRCAUTER': 'C'}\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c277a6a8093a4592a2a3bbf953cb5f27", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/5 [00:00