examples/CausalBasedSimulator_estimating_functional_relationships.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "32812f50-ddcb-4e18-8fdb-8b2b4c11bf39",
   "metadata": {},
   "source": [
    "# Estimating functional relationship\n",
    "\n",
    "This notebook details the method train(), which estimates the functional relationships between variables."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a10e8787-1ede-470d-b52f-de25706b857d",
   "metadata": {},
   "source": [
    "## Import and settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5196757c-0ab0-4f54-ae3e-2a28a7fa2702",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from scipy.special import expit\n",
    "\n",
    "from IPython.display import display, Markdown\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "\n",
    "from lingam.causal_based_simulator import CausalBasedSimulator\n",
    "from lingam.utils import make_dot\n",
    "\n",
    "np.set_printoptions(precision=3, suppress=True)\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a275a04d-400d-4601-b369-715700c50c63",
   "metadata": {},
   "source": [
    "## Creating test data\n",
    "Create test data with variables from x0 to x5, of which x2 is a discrete variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "766331ef-fd25-4401-b6e5-ec4812016296",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "### Causal graph"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"263pt\" height=\"218pt\"\n",
       " viewBox=\"0.00 0.00 263.00 218.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 214)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-214 259,-214 259,4 -4,4\"/>\n",
       "<!-- x0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>x0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"154\" cy=\"-105\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"154\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0</text>\n",
       "</g>\n",
       "<!-- x1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>x1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"137\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"137\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x1</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x1 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M154.32,-86.69C154.16,-76.97 153.39,-64.69 151,-54 150.32,-50.97 149.41,-47.86 148.38,-44.8\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"151.63,-43.5 144.81,-35.39 145.09,-45.98 151.63,-43.5\"/>\n",
       "<text text-anchor=\"middle\" x=\"169\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">3.00</text>\n",
       "</g>\n",
       "<!-- x4 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>x4</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x4</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x4 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x4</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M138.58,-89.96C124.52,-77.3 104.67,-59.81 96,-54 84.29,-46.15 70.59,-38.86 58.49,-32.98\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"59.74,-29.7 49.2,-28.59 56.75,-36.03 59.74,-29.7\"/>\n",
       "<text text-anchor=\"middle\" x=\"130\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">8.00</text>\n",
       "</g>\n",
       "<!-- x5 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>x5</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"228\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"228\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x5</text>\n",
       "</g>\n",
       "<!-- x0&#45;&gt;x5 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>x0&#45;&gt;x5</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M169.06,-89.74C175.4,-83.54 182.74,-76.08 189,-69 196.33,-60.71 203.91,-51.22 210.42,-42.76\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"213.27,-44.8 216.53,-34.71 207.69,-40.56 213.27,-44.8\"/>\n",
       "<text text-anchor=\"middle\" x=\"217\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">4.00</text>\n",
       "</g>\n",
       "<!-- x2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>x2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63\" cy=\"-105\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-101.3\" font-family=\"Times,serif\" font-size=\"14.00\">x2</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>x2&#45;&gt;x1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M59.11,-86.95C57.73,-76.63 57.8,-63.64 64,-54 72.53,-40.75 87.59,-32.39 101.83,-27.18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"103.03,-30.47 111.46,-24.06 100.87,-23.81 103.03,-30.47\"/>\n",
       "<text text-anchor=\"middle\" x=\"80\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">1.00</text>\n",
       "</g>\n",
       "<!-- x2&#45;&gt;x4 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>x2&#45;&gt;x4</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M40.01,-95.39C28.49,-89.76 15.62,-81.17 9,-69 4.65,-60.99 6.34,-51.71 10.05,-43.34\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"13.21,-44.85 14.89,-34.39 7.06,-41.52 13.21,-44.85\"/>\n",
       "<text text-anchor=\"middle\" x=\"25\" y=\"-57.8\" font-family=\"Times,serif\" font-size=\"14.00\">1.00</text>\n",
       "</g>\n",
       "<!-- x3 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>x3</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"98\" cy=\"-192\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"98\" y=\"-188.3\" font-family=\"Times,serif\" font-size=\"14.00\">x3</text>\n",
       "</g>\n",
       "<!-- x3&#45;&gt;x0 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>x3&#45;&gt;x0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M108.53,-175.01C116.82,-162.43 128.52,-144.68 137.97,-130.33\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"140.94,-132.19 143.51,-121.92 135.09,-128.34 140.94,-132.19\"/>\n",
       "<text text-anchor=\"middle\" x=\"145\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">3.00</text>\n",
       "</g>\n",
       "<!-- x3&#45;&gt;x2 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>x3&#45;&gt;x2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M86.85,-175.54C82.99,-169.64 78.9,-162.72 76,-156 72.87,-148.75 70.36,-140.58 68.4,-132.98\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"71.78,-132.03 66.09,-123.1 64.96,-133.63 71.78,-132.03\"/>\n",
       "<text text-anchor=\"middle\" x=\"92\" y=\"-144.8\" font-family=\"Times,serif\" font-size=\"14.00\">1.50</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f992c29cbe0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "### Test data"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x0</th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "      <th>x4</th>\n",
       "      <th>x5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.725033</td>\n",
       "      <td>-0.853352</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.298043</td>\n",
       "      <td>-5.518573</td>\n",
       "      <td>-3.361619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.092471</td>\n",
       "      <td>5.580223</td>\n",
       "      <td>1</td>\n",
       "      <td>0.449011</td>\n",
       "      <td>17.970192</td>\n",
       "      <td>7.668833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.251115</td>\n",
       "      <td>10.669606</td>\n",
       "      <td>1</td>\n",
       "      <td>0.965044</td>\n",
       "      <td>25.752890</td>\n",
       "      <td>11.510047</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.808976</td>\n",
       "      <td>13.150129</td>\n",
       "      <td>1</td>\n",
       "      <td>1.217832</td>\n",
       "      <td>30.951193</td>\n",
       "      <td>14.725043</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.023790</td>\n",
       "      <td>8.491654</td>\n",
       "      <td>1</td>\n",
       "      <td>1.096086</td>\n",
       "      <td>23.642615</td>\n",
       "      <td>11.167459</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         x0         x1  x2        x3         x4         x5\n",
       "0 -0.725033  -0.853352   1 -0.298043  -5.518573  -3.361619\n",
       "1  2.092471   5.580223   1  0.449011  17.970192   7.668833\n",
       "2  3.251115  10.669606   1  0.965044  25.752890  11.510047\n",
       "3  3.808976  13.150129   1  1.217832  30.951193  14.725043\n",
       "4  3.023790   8.491654   1  1.096086  23.642615  11.167459"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "display(Markdown(\"### Causal graph\"))\n",
    "\n",
    "G = np.array([[0.0, 0.0, 0.0, 3.0, 0.0, 0.0],\n",
    "              [3.0, 0.0, 1.0, 0.0, 0.0, 0.0],\n",
    "              [0.0, 0.0, 0.0, 1.5, 0.0, 0.0],\n",
    "              [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "              [8.0, 0.0, 1.0, 0.0, 0.0, 0.0],\n",
    "              [4.0, 0.0, 0.0, 0.0, 0.0, 0.0]])\n",
    "\n",
    "causal_order = [3, 2, 0, 4, 1, 5]\n",
    "\n",
    "display(make_dot(G))\n",
    "\n",
    "display(Markdown(\"### Test data\"))\n",
    "\n",
    "N = 1000\n",
    "n_features = len(G)\n",
    "discrete_index = 2\n",
    "\n",
    "e = np.random.uniform(-np.sqrt(3), np.sqrt(3), size=(n_features, N))\n",
    "X = np.empty((n_features, N))\n",
    "for i in causal_order:\n",
    "    X[i] = G[i] @ X + e[i]\n",
    "    if i == discrete_index:\n",
    "        X[i] = (expit(X[i]) >= np.random.uniform(size=N)).astype(int)\n",
    "X = pd.DataFrame(X.T, columns=[f\"x{i}\" for i in range(n_features)])\n",
    "\n",
    "X[f\"x{discrete_index}\"] = X[f\"x{discrete_index}\"].astype(int)\n",
    "\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f72f6d67-e0af-4cd5-8a52-573272416a53",
   "metadata": {},
   "source": [
    "## Estimate functional relationships between variables\n",
    "\n",
    "By executing train() given test data and a causal graph, the functional relationships between each variable in the test data can be estimated. The model representing the functional relationships of each variable can be specified in the models argument. If not specified, a linear or logistic regression model is used by default.\n",
    "\n",
    "The explanatory variables used for the model representing the functional relationships are determined by analyzing the causal graph in train(). To train each model, a column of explanatory and objective variables is taken from the given *X* and used. The most upstream variables do not have models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "af6b5e19-49f4-4295-a4dc-23ac5290abfb",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<lingam.causal_based_simulator.CausalBasedSimulator at 0x7f98d7234e50>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim = CausalBasedSimulator()\n",
    "sim.train(X, G, is_discrete=[False, False, True, False, False, False])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9672527c-67af-4de9-b8c1-b2e60083421e",
   "metadata": {},
   "source": [
    "You can check *sim.train_result_* for the estimated model and a list of variables that are inputs to the estimated model. *train_reuslt_* is a dictionary whose keys are variable names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "65ec4b64-7092-4aa5-b528-6891b9fde564",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['x0', 'x1', 'x2', 'x3', 'x4', 'x5'])\n"
     ]
    }
   ],
   "source": [
    "print(sim.train_result_.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b856a93-7936-485d-856d-1942ddcfbde5",
   "metadata": {},
   "source": [
    "### Contents of estimated results for x1\n",
    "\n",
    "`sim.train_result_[\"x1\"]` is a dictionary. This dictionary contains four items."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f528cd72-764e-4c6c-b726-12260e9500e9",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['model', 'parent_names', 'predicted', 'residual'])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.train_result_[\"x1\"].keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a996b71f-a17d-4b2d-9501-17efbd4ad174",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "- model"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LinearRegression()\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "- parent_names"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['x0', 'x2']\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "- predicted"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.181  7.272 10.749 12.423 10.067]\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "- residual\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.328 -1.692 -0.079  0.728 -1.575]\n"
     ]
    }
   ],
   "source": [
    "display(Markdown(\"- model\"))\n",
    "print(sim.train_result_[\"x1\"][\"model\"])\n",
    "\n",
    "display(Markdown(\"- parent_names\"))\n",
    "print(sim.train_result_[\"x1\"][\"parent_names\"])\n",
    "\n",
    "display(Markdown(\"- predicted\"))\n",
    "print(sim.train_result_[\"x1\"][\"predicted\"][:5])\n",
    "\n",
    "display(Markdown(\"- residual\\n\"))\n",
    "print(sim.train_result_[\"x1\"][\"residual\"][:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80e07a2a-bf08-4b8e-86df-2e5fad8a8096",
   "metadata": {},
   "source": [
    "## Estimate functional relationships between variables by specifying models\n",
    "\n",
    "The functional relationships between variables is estimated by default using linear regression when the objective variable is a continuous variable and logistic regression when it is a discrete variable, but another machine learning model can be used by using the *models* option. Here we set up *x1* and *x2* to apply the random forest model as the models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cf7133b2-2cb3-47cf-85c7-3d0150f9e9fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x1 RandomForestRegressor()\n",
      "x2 RandomForestClassifier()\n"
     ]
    }
   ],
   "source": [
    "models = {\n",
    "    \"x1\": RandomForestRegressor(),\n",
    "    \"x2\": RandomForestClassifier(),\n",
    "}\n",
    "\n",
    "sim2 = CausalBasedSimulator()\n",
    "sim2.train(X, G, models=models, is_discrete=X.dtypes==int)\n",
    "\n",
    "print(\"x1\", sim2.train_result_[\"x1\"][\"model\"])\n",
    "print(\"x2\", sim2.train_result_[\"x2\"][\"model\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}