sparx.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This analysis attempts to model the causal relationships between excercise and Parkinsons disease progression metrics from the SPARX clinical trial. The goal of this analysis is to develop a \"better than\" naive model from which further improvements can be made. This analysis relies on the causalnex package which contains methods for causal inference with bayesian networks. Although, the MDS UPDRS ME score is used as the endpoint in this analysis, a number of endpoints could be used in its place. \n",
    "\n",
    "This notebook requires the dbsas.pgsql databse connection (See notes.md).\n",
    "\n",
    "Start by setting up the python environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: psycopg2 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (2.8.4)\n",
      "Requirement already satisfied: causalnex in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (0.4.3)\n",
      "Requirement already satisfied: pgmpy==0.1.6 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (0.1.6)\n",
      "Requirement already satisfied: pandas==0.24.0 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (0.24.0)\n",
      "Requirement already satisfied: matplotlib<4.0,>=3.0.3 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (3.1.1)\n",
      "Requirement already satisfied: scipy<1.3,>=1.2.0 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (1.2.3)\n",
      "Requirement already satisfied: prettytable==0.7.2 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (0.7.2)\n",
      "Requirement already satisfied: numpy<2.0,>=1.14.2 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (1.17.2)\n",
      "Requirement already satisfied: wrapt<1.12,>=1.11.0 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (1.11.2)\n",
      "Requirement already satisfied: scikit-learn==0.20.2 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (0.20.2)\n",
      "Requirement already satisfied: networkx==2.2 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from causalnex) (2.2)\n",
      "Requirement already satisfied: python-dateutil>=2.5.0 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from pandas==0.24.0->causalnex) (2.8.0)\n",
      "Requirement already satisfied: pytz>=2011k in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from pandas==0.24.0->causalnex) (2019.3)\n",
      "Requirement already satisfied: cycler>=0.10 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib<4.0,>=3.0.3->causalnex) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib<4.0,>=3.0.3->causalnex) (1.1.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib<4.0,>=3.0.3->causalnex) (2.4.2)\n",
      "Requirement already satisfied: decorator>=4.3.0 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from networkx==2.2->causalnex) (4.4.0)\n",
      "Requirement already satisfied: six>=1.5 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from python-dateutil>=2.5.0->pandas==0.24.0->causalnex) (1.12.0)\n",
      "Requirement already satisfied: setuptools in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib<4.0,>=3.0.3->causalnex) (41.4.0)\n",
      "Requirement already satisfied: pandas in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (0.24.0)\n",
      "Requirement already satisfied: pytz>=2011k in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from pandas) (2019.3)\n",
      "Requirement already satisfied: numpy>=1.12.0 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from pandas) (1.17.2)\n",
      "Requirement already satisfied: python-dateutil>=2.5.0 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from pandas) (2.8.0)\n",
      "Requirement already satisfied: six>=1.5 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from python-dateutil>=2.5.0->pandas) (1.12.0)\n",
      "Collecting numpy==1.17\n",
      "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/c1/4b/78119133136c20e5ad2e01bf72b0633241defd619939908223cd394a9c32/numpy-1.17.0-cp37-cp37m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (15.0MB)\n",
      "\u001b[K     |████████████████████████████████| 15.0MB 2.2MB/s eta 0:00:01\n",
      "\u001b[?25hInstalling collected packages: numpy\n",
      "  Found existing installation: numpy 1.17.2\n",
      "    Uninstalling numpy-1.17.2:\n",
      "      Successfully uninstalled numpy-1.17.2\n",
      "Successfully installed numpy-1.17.0\n",
      "Requirement already satisfied: matplotlib in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (3.1.1)\n",
      "Requirement already satisfied: cycler>=0.10 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib) (1.1.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib) (2.4.2)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib) (2.8.0)\n",
      "Requirement already satisfied: numpy>=1.11 in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from matplotlib) (1.17.0)\n",
      "Requirement already satisfied: six in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from cycler>=0.10->matplotlib) (1.12.0)\n",
      "Requirement already satisfied: setuptools in /Users/kyle/opt/anaconda3/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib) (41.4.0)\n"
     ]
    }
   ],
   "source": [
    "!pip install psycopg2\n",
    "!pip install causalnex\n",
    "!pip install pandas\n",
    "!pip install numpy==1.17\n",
    "!pip install matplotlib\n",
    "\n",
    "# db module to connect to postgres db with all tables, instructions in README.md\n",
    "from db import connect\n",
    "conn = connect()\n",
    "\n",
    "# imports\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "\n",
    "# from imports\n",
    "from collections import Counter\n",
    "\n",
    "\n",
    "# causalnex\n",
    "from causalnex.plots import plot_structure\n",
    "from causalnex.structure import StructureModel\n",
    "from causalnex.structure.notears import from_pandas\n",
    "from causalnex.network import BayesianNetwork\n",
    "from causalnex.discretiser import Discretiser\n",
    "from causalnex.evaluation import classification_report\n",
    "from causalnex.evaluation import roc_auc\n",
    "from causalnex.inference import InferenceEngine\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.model_selection import train_test_split\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, a quick sanity check to make sure endpoints match with reported endpoints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MDS-UPDRS mean subscores\n",
      "Part 1 - Mean High Intensity 6 months:  1.0263157894736843\n",
      "Part 2 - Mean High Intensity 6 months:  1.5\n",
      "Part 1 - Mean Moderate Intensity 6 months:  -0.11904761904761904\n",
      "Part 2 - Mean Moderate Intensity 6 months:  0.7619047619047619\n",
      "Part 1 - Mean Control 6 months:  0.6216216216216216\n",
      "Part 2 - Mean Control 6 months:  0.8648648648648649\n"
     ]
    }
   ],
   "source": [
    "mdsupdrs = pd.read_sql('SELECT t1.idnum, t1.timepoint, t1.mds_part1, t1.mds_part2, t1.mds_part3, t1.mds_part4, t2.typerandomization1 \\\n",
    "FROM mdsupdrs AS t1 LEFT JOIN participant AS t2 ON t1.idnum = t2.idnum', conn)\n",
    "\n",
    "\n",
    "# baseline (timepoint 4) and 6 months (timepoint 11)\n",
    "data4 = mdsupdrs[(mdsupdrs['timepoint'] == 4)]\n",
    "data11 = mdsupdrs[(mdsupdrs['timepoint'] == 11)]\n",
    "\n",
    "df = data11.merge(data4, on=('idnum'), suffixes=('_6m', '_base'))\n",
    "\n",
    "hi = df[df['typerandomization1_base'] == 1]\n",
    "med = df[df['typerandomization1_base'] == 2]\n",
    "con = df[df['typerandomization1_base'] == 3]\n",
    "\n",
    "print(\"MDS-UPDRS mean subscores\")\n",
    "print(\n",
    "    \"Part 1 - Mean High Intensity 6 months: \",\n",
    "    np.mean(hi.mds_part1_6m - hi.mds_part1_base)\n",
    ")\n",
    "print(\n",
    "    \"Part 2 - Mean High Intensity 6 months: \",\n",
    "    np.mean(hi.mds_part2_6m - hi.mds_part2_base)\n",
    ")\n",
    "print(\n",
    "    \"Part 1 - Mean Moderate Intensity 6 months: \",\n",
    "    np.mean(med.mds_part1_6m - med.mds_part1_base)\n",
    ")\n",
    "print(\n",
    "    \"Part 2 - Mean Moderate Intensity 6 months: \",\n",
    "    np.mean(med.mds_part2_6m - med.mds_part2_base)\n",
    ")\n",
    "print(\n",
    "    \"Part 1 - Mean Control 6 months: \",\n",
    "    np.mean(con.mds_part1_6m - con.mds_part1_base)\n",
    ")\n",
    "print(\n",
    "    \"Part 2 - Mean Control 6 months: \",\n",
    "    np.mean(con.mds_part2_6m - con.mds_part2_base)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MDS-UPDRS mean subscores\n",
    "\n",
    "Part 1 - Mean High Intensity 6 months:  1.0263157894736843\n",
    "\n",
    "Part 2 - Mean High Intensity 6 months:  1.5\n",
    "\n",
    "Part 1 - Mean Moderate Intensity 6 months:  -0.11904761904761904\n",
    "\n",
    "Part 2 - Mean Moderate Intensity 6 months:  0.7619047619047619\n",
    "\n",
    "Part 1 - Mean Control 6 months:  0.6216216216216216\n",
    "\n",
    "Part 2 - Mean Control 6 months:  0.8648648648648649\n",
    "\n",
    "These values are approximately what is reported in the study publication, so it seem sensible to proceed. However, it is a little concerning that these values are not precisely the same and that could be investigated further."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SPARX study design is a randomized controlled trial with 2 interventions;  moderate-intensity and high-intensity endurance exercise with a control group (no excercise). The primary endpoint is taken at 6 months with participants encouraged to continue to excercise an additonal 6 months.\n",
    "\n",
    "Hence, a naive bayesian network to model this study for the first 6 months could take the following structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kyle/opt/anaconda3/lib/python3.7/site-packages/networkx/drawing/nx_pylab.py:563: MatplotlibDeprecationWarning: \n",
      "The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.\n",
      "  if not cb.iterable(width):\n",
      "/Users/kyle/opt/anaconda3/lib/python3.7/site-packages/networkx/drawing/nx_pylab.py:660: MatplotlibDeprecationWarning: \n",
      "The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.\n",
      "  if cb.iterable(node_size):  # many node sizes\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAOHklEQVR4nO3dfbBtdV3H8c/3clEuIhiCD6GDjxMaY06JOo4PlY6ZoymipoKIaT6kMTaY5lgNo6jTAzpmJaQIKmoYoqbppKZMaU6KNZL4lOYDXK5CNxEucEXk1x9rHe/mdu49h8PT9+jrNbOHe9Zev3XW3tz7e++19jr71BgjANDNhpt7BwBgOQIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASwIFQEsCBUBLAgVASxtv7h1YUdXtkjwjyX2S3CbJJUnOTXJaxrj45tw1gHVjHc6lNca4ufdheVWHJXlZkl9PMpJsWrj3yiSV5MNJXpMxPnvT7yDAOrCO59Kegap6XpITk+yV3Z+GvCbJ9iTHZYyTbopdA1g31vlc2u89qB1P6N5Zef82zOudOI+7EXanvllVj7gxtg2wqKqOrKqP3EAbu1Hn0qoaVXWP67eTK+/UmlTVw6vqy1V1RVV9oqoOXuW4R1TVv1fV5VV1flU9eeHOw7LjCb0ulp7Y+13HcQBtjDHeMcZ45PXe0HWcS89Ocqfpj63m0jUFqqoOSHJWkj9Ksn+Sc5KcsYpx907yziQvT7Jfkvsm+dzCKi/LdCi6FnvN428yVdX/IhNgXbiB55N1NZfuylqPoJ6Q5Lwxxt+NMbYnOT7JL1TVISuM+8MkJ48xPjzGuHqMsXWM8fUkOabqlyo5/NRkw52T/EySk5J8NjsuOXnhyo/l0QdU/WlVnV5Vt6yqF1TVl+ZD0Y1JUlVnV9VrquozVfX9qnp/Ve2/tJGqenpVfauqtlbVyxe/QVUdX1Vnztu/dNrtun9VnVNVl1bVd6vqtat9EoE+qupnq+o9VXVxVX2jqo6dl3+oqk5cWO+MqnrL/OdjqupTVfWGeT75clU9fGHd/arqlKraUlWbq+qEqtpjp7Gvq6r/TXL8vOyTC+NHVf1OVf1XVV1WVa+sqrtX1afnOefdVXWLhfUfs2fVF26THP6gZMO5C4/vLkn+PNN8ul+S38z0ptPlma6euDDJPtNtw+bk0Q+ueuT8fS6Z9/8vF7/XKp/T21bVB+Z9/ez8+Hd+fC/Y5QbGGNf5luT1Sd6407IvJDlihXH/neSVSf4zyZYkpyfZf4yRf05enWQ8NxlXJuMfk3HLZDwuGd9NxgXJODAZZydj7Pp2xW8kH03ymSRbk2xLclWmK1c2zvtwdpLNSQ5Ncqsk70ly+nzfvecxD01yyySvTXJ1kkfM9x+f5IdJHp8piJuSfDrJ0+f790nywLU8p25ubjffbf73/Lkkf5zkFknuNs9Xv5bkDkkuSvKrSY6cl996HnfMPEf8XpI9M83731+a15K8L8nJ81xzu3lueu5OY38304/8bJqXfXJhv0aSv0+yb5KfT/KDJP80799+Sb6Y5Bnzur+Y5KJTk7/4YXLFack4OBnb5/nx4GQclozNydiajEOS8cb5vk8k46Cd5tLTk9cleeC8b3dJ8qUkL9pp3+6xwvP6t/Nt73l+PX+Zx/fRXY1f01V8VXVKkovHGH+wsOxTSd40xjhtN+OuyhTqR87/fWuS7WOMIz9fddZ9k8MvSHLQvP5tk/x1pv/jSXJEkockedFu9u0JyY/em+yxsGhkuoxyy8Jmr0py2fz1xiQHzvfvk+kv2feWdjnTX86t85hbZ/rLu3Vh+0vbuzzTlTDA+rNnphM3Fy0s2yfT/HBJptNe+83Lv5fp33wyRWXfJN9dGHdApvngB0lunx1zz9L6e2eaQzZlmlMu2sX9SXLHJP+T6YXx0ra3Z3ohnfl7J8ml8/5dc0ayx5PnS8l/LsnfJHlYpsKckOSoecBL5kEnZXrVflSSC679nLw9Yxy99EVVvSjJw8YYh89fjyT3HGN8LcuYjxS3Jzl0jPGVedkJSX55jPHghW08fIzx8eW2sdZzntuy44lZsm92TPq7cmWSU8cYX5137tVJPpYke8zbu/3CypuW+Xpbdu+A6TnekCk6GzNF4xZJHpDkR0nOTPLBJKctbPZrSR6V5MXztzhhYZP/keTYJP+S5Lgkd821zzbeNcnvZzrq+namo66PrbCbQC+PTfJXmY50lmzIfIYkU8D+LdO/8ccvrPPkTEc9j15YdnKSzyf510xzzc7b/GaSX5nHPm2Z7S0u25zkcfOYZDoie2eSd89fvyTTkdmLk7w9yYN+K9n4nPnOpSOCJXdY+PPeO923s3OSgw6r+mCS+82rb8y1rxlYydIcfP7CsvOXWW+5ZUnWHqjzMv1EcpKkqm6V5O7z8t05N9MRzf/zoynm19te06uNb2f6H3ZCkifNd20ZY1w9H8XtO8a4cN73e2V6dXJekq8nudfCfXtnelW1dYxxYVVdluSKpftnFyb5VFVtyPTe3OlJbjvGuPyGeDzAja+qzk3yjTHGPXdx/59lmiPumuko4l3z8ksyBWLLmE9H1fSJDV/J9OL2B5nmg6uX2eYlSa5anE92XlZVSXLRwtdXJblk4ettSfaZ56evJfn4tuSQ7DhQWt3jX2bZs6ZTcu9I8tQxxmXzEdQTr8NmL850CvNOSb46L7vzMuvt8jTeWi+SeG+SQ6vqiKraK9N523PHGF9eYdypSZ5ZVXebJ/+XZnqFkUuTlcauxpUHT4V/aKYn5vmZfkJ6Z0dV1b3nfXhFkjPHGEtHV4+pqgfPbwa+Iis8R1V1VFUdOMa4JtOpgGQ6UgPWj88kubSqXlpVm6pqj6o6tKoOq6qHJnlmkqPn2xuq6qCFsbdLcmxV7VlVT0pyryQfGmNsSfKRJCdW1b5VtWG+wOFhN9JjeFOS552WfO+a5MrLk/xDVj6tlUxnqrZmevNsduXW6YX7pUm2zRfAPf+67Mw8p56V6eKPvedtHL3CsGtZU6DG9LlNRyR5VabzsQ9I8pRVjHtLkrdlOlT+VqZXF8cmySnTxQrXVx03XSl4Rqajtc8l+cAy67090ym+72Q6t3zsvH/nJXlBpkPoLZke2wXLjF/0qCTnza9kXp/kKWO6shFYJ+bJ9LGZfvTlG5nOxLw503tAb0vywjHG5jHGJ5OckuTUmg9vMs1n95zHvCrJE8cYS+8hHZ3pLYYvZppPzpy3eWM8hnOS/PZzkofsn2y6R3a8j7GSQ5I8NdOVF7dJsjmpO05z4dMyNe5NWcWPEi3jhZneG/tOpnn3XZnm/VXp9VFHVWdlOt+6lnBek+R9GeOI3X+LOjvTVXtvXsP3APixqjomybOX3vRv4yaYS9eiqv4kyR3GGM9YceX0+6ij12S66mMtts/jAX7atZhLq+qQqrpPTe6f5FmZ3iJalZXeX9m2zG2593RWGrOtqh6y4t5Mn6R7XJIrdrfa8/LjHyj78W3PZI9Knr3i9wD4SbfKuXQZV2T6wNhzVjugqs7bxZx/ZKbL6M/KdNn9uzN9/NL7V73tVqf4lqzzT+AFaGGdz6U9A5Vk/rDCl2X6+YJd/Q6TD2X6HSarrj3AT5V1PJf2DdSSqgOz/G+BfGua/hZIgHbW4VzaP1AA/FTqdhUfACQRKACaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWhIoAFoSKABaEigAWvo/pX9q76MknKcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sm = StructureModel()\n",
    "sm.add_edges_from([\n",
    "    ('experimental_group','d_0_6m_updrs')\n",
    "])\n",
    "_ = plot_structure(sm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where d_0_6m_updrs is the difference between the 6 month UPDRS score and the base line UPDRS score and the experimental group is one of high intensity, moderate intensity, or control."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally, we may want to include the 12 month endpoint. Also, including data for the baseline, 3 month, and 6 month endpoint may be more informative. For example, subjects may have different rates of disease progression or slowing based on the baseline. There may also be differences in the individuals response size to the treatment. \n",
    "\n",
    "A model could be made with timepoint -> score and experimental_group -> score, but this would not give us the conditional probalility distributions for individuals with different rates of disease progression/slowing. Therefore, it seem appropriate to model a given enpoint as its difference from its baseline.\n",
    "\n",
    "A network to model this could take the following form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2debyc4/n/358sskiInVhi3wWRtKp8dftSsX9RaqeK2qriS6mWKkVR/aKtpbXXVluJUNVf7VXE1tojRCSxC7Igkev3x3VPzjMzJ2c/53lm5nq/XvPKycw8z3PN5GTec9/XfV+XzIwgCIIgKBq98g4gCIIgCJojBBUEQRAUkhBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSSEFQQBEFQSEJQQRAEQSEJQQVBEASFJAQVBEEQFJIQVBAEQVBIQlBBEARBIQlBBUEQBIUkBBUEQRAUkhBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSSEFQQBEFQSEJQQRAEQSHpk3cAQRB0M9LSwH7AcGAIMB14FrgCs3fzDC0IWkJmlncMQRB0B9Io4ARgG8CAAZlHZwMC7gLOwOzxng8wCFomBBUE9Yh0KHAu0J+Wp/LnAZ8CYzC7qCdCC4K2EoIKgnqjSU4D23HULEJSQcGIRRJB0EVIukLSaW143uuSvtVNQYyi/XIiPf9cpJFdH1QQdIwQVBDkiKRTJF3TBedZV9ITg+DhxWDgt4Dn23+a/njOKggKQQgqCOqDqfvBIR/DF+8BOwB7tP8cvYDRSEt1cWxtQk58JgXziV+GIOggkjaW9KSkTyTdgI9ASo9tJ+lpSdMlPSJpeDPHfxs4Edhd0gxJz6T7D5D0QjrvREmHtBaLmU2/Ar7RC8yA3sCEzOP7A4fhy/kGAV8F3gKOBhYD1gaeSqfCl6Qv6DUfL2lKiu0lSd9M9/eWdKKkV9Nj4yWtmB7bTNLjkj5Kf26WOd99kk6X9DCeB1tV0qKS/ihpWrrWaZJ6t/YeBPVHCCoIOoCkhYDbgKuBxYE/A7ukx0YAlwGHAEsAFwO3S+qXPYeZ3Q38ErjBzAaZ2YbpoXeA7YBFgAOA89I5W2P4EBjQHzgSN1+WG4HTgPeAfsBXgBHp77sCx/jTBuD7pZp7zWsBRwCjzGwwsDXwenr4GOC7wOgU94HALEmLA3cC56f34tfAnZKWyJx6H+BgYDAwCbgSmAusDmwMbAUc1IbXH9QZIagg6BibAn2B35jZHDO7CSjtJfo+cLGZ/cvMvjCzK4HP0jGtYmZ3mtmr5twP3ANs0dIxkvq9C8OmAx8BF+Kf7Fl2BjbBh3k7pz/3xUdbuzN/BAW+mbc5vsDdtq6kvmb2upm9mh47CDjJzF5KcT9jZu8D2wKvmNnVZjbXzK4DXgS2z5z3CjN7zszm4rLfBjjazGaa2TvAeXRoxjKodaKSRBB0jKHAFCvfpzEp/TkM2E/SkZnHFkrHtIqkbYCTgTXxL5EDgX9nHu8HrA+MxJ2zCbDOQ/DpzsDCwKHAUsALwNLpuGUy1xjQzN9nNP11enNxmdkESUcDpwDrSforcIyZTQVWBF5t5rChNL0vJSYBy2f+Pjnz8zBc/NMkle7rVfGcoEGIEVQQdIxpwPLKfIoCK6U/JwOnm9mQzG1gGj1UUrYRMcnnZuAc3CFLAw8DIyRdLGk88CFwOT4iexY4HFhiZzgDrxDBPDyhM6X9r2t2OmezmNm1ZrY5LhIDzsq85tWaOWRqem6WlSpCy74Hk/HR5pKZ924RM1uvfS8jqAdCUEHQMf6J50mOktRH0v8AX0qPXQocKunLaWXawpK2lTS4mfO8DawsqVfKa42iafbtn8DHwJbACvgo6gj8w3u4mR1gZhea2aPA5t+C8XNBH+MJocWAddr/uoTngKofkNaS9I0k0U9xmX2RHv4D8AtJa6TXPDzlmcYBa0raM71PuwPrAmObu4aZTcOnNM+VtEh6X1aTtGX7X0pQ68QUXxB0ADP7PEnpUnztwTjglvTYE5K+j6eC1sA/yB8CHsicoo+kjXEZrQZ8jn9hfB54DPgaLotb0v0TzOzCFkIa8nf4xWDoOwi33N1klhW2nXEtFJDtB5yJe28O8Ai+uAF88UM/XC5L4nmmnc3sTUnbAf8H/B5fXLidmb3XQgz7pus8jy+cmEjTSC1oIKLUURB0M5L6AutRnjNaH//gHZ9uTwDPmNnMTl5sFHAf7a8kwecwZ3u4/x44G3jEzGa0elAQdCMhqCDoQjIyKoloZPr765TL6OlOy2jBQXSoFt+F8I8jfdXdp/jivleBM8zsqm6IMghaJQQVBB0kyWhdmkRUGhlNoklE43EZdXo0IulEqrc3ATxoZttUPLnd1czluaf3KBfbZ3iuDWBdM3ujY9EHQfsJQQVBG5DUB5dRdppuA+ANmkRUktEnecVZhhd+PQHfPLugflDj8H5QT/ghOhdfiNEbX2k33Mw+6smwg6BECCoIKkgyWodyGQ3HZVQ5TVcMGbWE19bb71XY+jn41hZwy2K+QvDKygURkoYBrwH3Ai8DX8cXNbzW43EHDU8IKmhoMjLK5ow2AN6kXEZP1YSMWkDSY/gCv6vMrKV6e6OB+81sZtpsfCLwP2b2zx4KNQiAEFTQQCQZrU25jIbjU1nZnNFTZvZxXnF2B5K+ii8BH4jnnNY2s8oKDws6dlt8Y/BRZnZ990UZBOWEoIK6JFW/LsmoNFW3IV7ZoFJGdZ9jkfQIXh8WfHPtn1oaRTVz/HDgDnxD7mkWHxxBDxAbdYOaJ8loLcpzRhvh5YhKIroVeLIRZFRJKsckfN/VyunPQe05h5k9K2lT4C94ZYiDzOyzro41CLLECCqoKTIyyk7TbYi3N6rMGTVb9LRRkbQYMNHMFuvEOQYCV+F1AndupSJEEHSKEFRQWJKM1qRcRhvh9euy03RPhoxapysElc7TCzgd2A1f4fdiV8QXBJWEoIJCkD70SjIqTdVtBLxLtYw+zCvOWqarBJU53wF4zbw9zezvXXHOIMgSggp6nCSjNSjPGW2MVzHIbnp90sw+yCvOeqOrBZXO+TXgBuAnZvaHrjpvEEAIKuhmMjLKTtOVZJTNGYWMupnuEFQ675p4W/dbgR+b2byuPH/QuISggi4jyWh1qmX0AdXTdO/nFWej0l2CSudeAm8N8gGwd7cVwg0aihBU0CGSjFajPGe0Md4uvFJGsdKrAHSnoNL5FwIuwStxbJ9awQdBhwlBBa2S9tGsRnnOaATwEeU5o/Eho+LS3YJK1xDwY+AwYAcze6q7rhXUPyGooIyMjLLTdCUZVeaMFtR5NSggPSGozLV2A34HHGhmd3T39YL6JCpJNDBJRqtSLaNPaBLRr/CRUcgoaDNm9mdJk4BbJa0O/CbKIwXtJUZQDUKS0SqU54xGADMpzxmNN7N38ooz6D56cgSVueYwvIbfI8CRZjanp64d1D4hqDokyWhlymW0CTCL6pzR2zmFGfQweQgqXXcR4Hp8xuY7UfUjaCshqBqnQkbZ26c0yegJQkYNT16CStfuA5wHfBPYNhogBm0hBFVDJBkNozxnVJJRdgHDeDN7K684g2KSp6AyMRwB/ATYxcweySuOoDYIQRWUJKOVqJ6m+5zqnNG0vOIMaociCCrFMRq4AvihmV2XZyxBsQlBFYAkoxWpltFcqnNGsfkx6BBFEVSKpdQA8TLg1FjhFzRHCKqHSTJagXIRjcS7nFbmjEJGQZdRJEEBSFoWb4D4CnCQmX2ac0hBwQhBdSMZGVXmjIwKGQFT41tk0J0UTVAAkgYAVwJD8QaIsd8umE8IqotIMlqe6mk6qMgZAVNCRkFPU0RBwfy6jr8A9sAbIL6Qc0hBQSh+JQlpaWA/YDgwBC9G+ixwBTl920oyGkq5jEYCoklEl6SfQ0ZB0AKpPcdPJL0C3C9pTzO7N++4gvwp7ghKGgWcAGyDT4kNyDw6G5fBXcAZmD3evaFoKNU5o95UT9O9GTIKikpRR1BZJG0J3AicZGaX5h1PkC/FFJR0KHAu0B/o1cIz5+F7gMZgdlHXXFqlkVF2dNSHin1GwOSQUVBL1IKgACStgTdAvB043sy+yDmkICeKJ6gmOQ1sx1GzqJBUmobbB3jWzJ5u/lJajuqc0UJU54zeCBkFtU6tCApA0uJ4A8TpwF7RALExaZegJL2OLwftnvlhn9a7j/bJqcQsYEvMnpA0EF8ZtAtwkZkdlpa0VuaM+lGxzwiYFDIK6pFaEhTMb4B4MbAh3gBxSs4hBT1MhwUl6RRgdTPbu9NBuFDOGQzfEyy0IfBA+08zD7hNcDzwN3wRw0LAx8AMPIdVmTMKGQUNQ60JCubPhBwPHE40QGw4irKK75KBsPAEmLcE0Ox8XOv0Mhi9FOz0bnneaiDeVmJid8ko/SdSWo0UBEEXkf7PnilpAnCPpO+Z2e15xxX0DC0tQEDSI6kkSeX93wZOBHaXNEPSM+n+AyS9IOkTSRMlHdJaAJLWAnZ4EcYvDdabps1DAPvjvaO3AQYBXwXeAo4GFgPWBkpfqQR2B/we/8b1GF63rk96ncdJmpJie0nSN9P1e0s6UdKr6bHxklZMj20m6XFJH6U/N8vEfZ+k0yU9jE8vrippUUl/lDQtXes0Sb1bew+CIGgZM7sJ2Bb4vaQx6UthUOe0KCh8/vd2Sf2yd5rZ3cAvgRvMbJCZbZgeegfYDlgEOAA4T9KIVq7xZWDSGNhzSRiwAXBzxRNuBE4D3sOTRl/Bh0TvAbsCxzQ9dcCXYREz+5WZfRnfN7UTPtV3BDDKzAYDWwOvp2OOAb4LjE5xHwjMSknaO4HzgSWAXwN3SloiE9o+wMHAYGASnveaC6wObAxsBRzUyusPgqANmNlj+H///YCLJPXNOaSgm2lRUGZ2JfAZsGlbTmZmd5rZq+bcD9wDbNHKYSsA6y8BvacCF+K/fdmt5Dvjo6r+6ef+wL74RqTdaRpBAXwKS0v6jqRrgMnAbcCyuNvWldTXzF43s1fTIQfhey5eSnE/Y2bv49/WXjGzq81sbqq6/CKwfeZyV5jZc2Y2F1gcH+gdbWYzU1fa8/Dd8UEQdAFm9gY+kbICME7SkJxDCrqR1qb4puNVtoe25WSStpH0qKQP0rGjgSVbOWw2MOc38PhCwJbA13GzlVgm8/OAZv4+I/P3P/vo6HpgL3zkMxd4FJ8VPAV4R9L1ab8T6fW9SjVD8VFRlkl4OaMSkzM/DwP6AtMkTU+v/2Jg6WbOHQRBBzGzT4AdgeeBRyStmnNIQTfR2ghqiJkNXEDPlrIFB2ka8GbgHGAZMxsCjMMrPrTEswDz4N+4rDrD7NU8B/UJXh2c9Oe7wBj8F/oUfOrv7PT4ZGC1Zs41FZdOlpWA7FLX7HswGR9tLpnetyFmtoiZrdfxlxMEQXOkWY0fAr8DHpb01bxjCrqe1kZQC0vaVtLgZh5+G1g5FXoEz/P0w2UwV9I2eA6mNR4A3hgGi80BPYxvhNq6ra+gIuTN4GRgTeBJXE7X4lOUl+GCGgF8CV/gMR74EPitpJ0kLSRpeMozjQPWlLSnpD6SdgfWBcY2d+HUNPAe4FxJi0jqJWm1VLolCIJuwMwuxPPGt0raM+94gq6ltWXmE4CHaH5b0p+BvYH3Jb1mZiMkHYWvaeiHNyNrdTmomc2RtOO78IfB0HdV4Cp8dV47mQeMw+xdw1fg4bIqbSreD1gHmAM8AhyFzxaOxBd03IhP0c3Cd7A/gC+gOBYflU3AKy2/10IM+wJn4iIcDEwEzmr/SwmCoK2Y2V2SvgGMTWWSogFinVCsUkedqCTxKcz7MZzwf3BZKxJp4fJaGN+1nq04sQounOwG3+fM7POOXCMI8qIWN+q2h0wDxAnA96IBYu1TLEFBh2rxzYPZh0P/i3wU9QW+VepO4Ki0wq4T4WhhYCPKC8iuAvyH8gKyz5nZnM5cKwi6k3oXFJQ1QFwe2CkaINY2PSIoSSfiG3sredDMtmnmgHZXMxf8N/A/FY+VFl2sm5andgkZaWWLzK5MtbSeD2kFRaERBAXzGyCeCuyJT8s/n3NIQQcp3giqhDQS7wc1mgX3gxqH94N6QtIGwL9wqRmwdU82PZM0iGppDcOllS1IG9IKcqFRBFVC0r74quJogFijFFdQJaSlgP3ugqPWhbeG+WbZZ4ErKzvqSnoCWA84DvgZcIiZ3dLjMTfFM5im6cGSuFbCl9Rnc1rPd3YqMghao9EEBSDpv/AFUCeb2cV5xxO0j+ILKiHpNrxyw20tPGcDYICZPZZKLP0FuAA4uyirepK0NqY8p7UiLt3s9OALIa2gK2lEQcH8Bohj0+24aIBYO9SVoJo5ZgV8uft44LCirrzLSCs7PbgC1dJ6MaQVdJRGFRTMb4B4E76Jfy8zm9HKIUEBqGtBpeMGAX/C9yXtYmYfdkd8XY2kRaiW1vK4tLI5rRfiG2HQFhpZUDC/AeLv8c3625vZmzmHFLRC3QsqHdsb+BVeAHY7M5vQ1fH1BJIWpWl6sCSuocAzlEvrxZBWUEmjCwrm9247DjgS2NHMxuccUtACRWlY2K2kD+sxkl4GHpK0m5k9mHdc7cXMPsI3Mt9Xui9JawQuq9HAT4FlU4+u7EKMl0JaQaOTctFnyRsg3i3p+x350hv0DA0hqBJmdrGkicBNko41s6vzjqmzJGn9I90ASC0IStLaFi/5tIykpynPab0c0goaETO7WdIk4DZJqwPnFmUhVdBEQ0zxNXOuUtHXa4GfNUKr9oy0SlODI/FWIE9RLa26fz8akZjiqyZ1zx6L76E8PPYoFouGFFQ639J4M8PJwP5m1tlWHzVH+sDKSmsTYCngacpzWiGtOiAE1TxpFe11+Cb/3WplIVUj0LCCSufsj7fhWBVPmL7dVeeuVdJy3NL0YElcS9I00iqJ65WQVm0RglowaSHVuXinn+0yHbeDHGloQaXzCq86sT++9PQ/XXn+eiD1xypJqySuJfCeW1lpTQhpFZcQVOtIOgxfaLSrmT2cdzyNTsMLKnP+vYDzgH3N7O7uuEY9kZFWNqe1GE3SKonr1ZBWMQhBtQ1J38bb0v3IzP6UdzyNTAiq/BpfxXeb/8LMftdd16lXJC1JdU5rCD49mM1phbRyIATVdiStj1ehuQo4JVb45UMIqvo6q+K9pP4KjIll2J0jSSs7NbgJsCg+0qqUVm38MtYoIaj2IWkZfCHVa8CB0QCx5wlBNX+tIXhL+8+A75rZJ919zUZCXqG+UlqDqc5pTQxpdR0hqPaTGiBejnch2MnM3sk5pIYiBLXg6/UFLgQ2xRdPdFnDw6CatOx/E8rFNYjyfNZ44LWQVscIQXWM1ADx58De+Aq/53IOqWEIQbV8TQE/Ao4Bdjazx3vq2kGZtLILMQbSJK3SLaTVBkJQnUPSPvhS9L3N7J6842kEQlBtu/YOwB+BQ83s5p6+ftBEygtkpbUJ3m25Mqf1ekirnBBU55G0BT79f4qZXZR3PPVOCKrt1y81QPwtcFZ8+BUHSctSndPqT/X04KRG/ncLQXUNqXbfWOAu4NhYSNV9hKDaF8Py+NLTp/HRVCEbIAYgaTmqc1oLUS2tNxpFWiGoriO9lzcBs/CFVNEAsRsIQbU/jkHANfhS6V3M7IM84wnaTkZa2ZxWX8qnBp8AJtejtEJQXUtaSPV7/HcpGiB2AyGojsXSGzgT2AFf1fNKziEFHUTSUKql1ZvyUdZ46kBaIaiuJy2kOhb4IdEAsctpqH5QXUWac/5fSa8AD0ra3czuzzuuoP2Y2VRgKj51W/rAKUlrE+Ag/FtyL0mV0nqz1qUVdI707392pgHiwWZ2a95x1QshqE5gZpekBog3SjrOzK7MO6agc6QPnCnpdjvMl9byNEnr4PQnSVpZcU0JaTUeZnarpDeAv6RFFOfE70HnCUF1EjO7V9LXgLGS1gR+GnXm6ov0QfNmuv0FyqRVmho8JP1skipzWlPjw6r+MbPxkjbFV/itKemwaIDYOSIH1UWk8j23AtPwiugN1wCx0UnSWoHqnNYXlI+ynkhTiz0dX+SgeoC0kOo6fFP5rtEAseOEoLqQ1ADxj8DqeML0rZxDCnImSWtFqjcXz6ViIUZ3SysE1XOkhVRnA6PxhVQTcg6pJokpvi7EzD6VtDfe8OxRSdub2b/zjivIjzS190a63QrzpbUSTbI6AthE0udU7NMys2l5xB10jrSQ6hhJLwMPSdrNzB7MO65aIwTVxaQPpFPTCr+/S9rPzO7KO66gOKTfkUnpdguUSas0yjoSGCnpMyr2acXIvHYws4vSQqqbJR1jZtfkHVMtEYLqJszsOkmT8F/M083swrxjCopLhbRuhvnSGkbT9OBR+EjrU6qnB0NaBcXM7pH0deAOSWsBJ8dCqrYROahuRtIqeAPEe4FjzGxuziEFNUyS1spU1x6cRbW03q44NnJQOZKq89+GT/ceEAupWicE1QOkBog34onxPczs45xDCuqIJK1VqJbWDMpzWhOAx0JQ+ZEWUl2Of8nYqfJLRFBOCKqHSHW7LgA2w1f1RAPEoNvISCu7cnAksAi+l2t+P63oEtuzpH+bU4B98Rp+/8k3ouISgupB0i/m0Xjtrp3N7LGcQwoaiDTF9zrwfZrENQL4hOrpwXdzCrNhkLQXcB6wj5n9Ne94ikgIKgdSA8Q/AIeZ2U15xxM0Bs3loNKXptUoH2WNAD6iesn7ez0edJ0jaXO8AeKpZvb7vOMpGrGKLwfM7HZJWwG3S1oDODNK4QR5kH7vJqTbDQCSelEurR8DIyRNJ1MNA5fW+3nEXS+Y2UNJUmPTCr8x0QCxiRhB5Uhq9XAH8G/g4GiAGHQnnVnFl5FWNqc1AviQ8n1aIa0OkBZS3QTMBvY0s09yDqkQhKByRtLCeAPExfAGiPGfO+gWunqZeZLW6pSvHNwY+IDqnFY09myFtJDqt8CX8YVUk3MOKXdCUAUg/Uc/E9gJ/8V8OeeQgjqkJ/ZBpd/lNSjPaW0MvEd5TuvJkFY1KSd4DPAjfBn6EzmHlCuRgyoAaVf5caluV6kB4n05hxUE7Sb9Lr+UbtfCfGmtSZO0fgZsLOldynNaTzZ65e+UEzxX0qvAXZIOMbNb8o4rL0JQBcLM/iDpNeAGST82s8vzjikIOkuS1ovp9ieYX+17DZqmBk/GpfUO5TmthpSWmd2WaYC4BvCrRlxIFYIqGGb2d0lb0tQA8SdRtyuoN9JKtZK0roH50iqNtEYC2wMbSXqb8pzWk2Y2PY+4exIzezI1QLwDb4D4g0ZbSBU5qIIiaUm8PcPbeAPEWTmHFNQ4tViLL0lrLcpzWhvhjUErc1of5RVnd5IaIP4JGIwvpGqYEWUIqsBI6odv6F0Lb4AYvYGCDlOLgmqOJK21KZfWhsBUyqX1VL1IK73ms/BR5baN0gAxpvgKjJl9Jmlf4CSaGiA+m3dcQZAnaXrwuXS7CsqkVcpp7QxsKGkq1TmtmivWnF7zsZkGiN8xswfyjqu7CUEVnJQY/UVqgHivpP3NbFzecQVBkaiQ1pUAkvrQNNIaCewCDJc0hfKc1lO1Ii0zuyQtpLpJ0rFmdlXeMXUnMcVXQ0j6Ct7M7gwzuyDveILaol6m+DpDktY6lE8PDgcmUz09WNhqDpLWAcYC1wE/q9eFVCGoGiM1QBwL/AM4OhogBm0lBNU8SVrrUi6tDXBpZacHCyUtSUvhDRCnAPvVYwPEEFQNImlRvAHiPGD3WpmeCPIlBNV2UtmhdSivPbgB3g23cnpwRo5x9gf+iNdJ3LHeGiCGoGqU9K3vfGALvDzSpJxDCgpOCKpzJGmVRlolca0PTKJcWk/3pLRSeaSfAftTZw0QY5FEjWJmcyUdDhwF/FPSzmb2r7zjCoJ6xczmAM+k22UwX1rr0TTK2gtYPy1kyOa0njazmd0UlwE/Twup/p+kfc3s7qonSksD++E5tyHAdOBZ4AoK2qAyRlB1gKTt8f8wh5vZjXnHExSTGEH1DBlpZacH1wdeozyn1eXSkvRVvG3HaWb223TnKOAEYBvAgAGZQ2YDAu4CzsDs8a6Mp7OEoOoESRsBtwMXA79sxLpdQcuEoPJD0kJUS2s9YCLV04OdqhojaVV8IdXf5sCLfeAcoD/Qq4XD5gGfAmMwu6gz1+9KQlB1hKTl8Lpdz+ENED/LOaSgQISgikWS1vqU57TWBV6lXFrPtFdakob8BB45CdbsD73bcegsCiSpEFSdkRogXg0sCewcDRCDEiGoYiDpdeAgM7u3mcf60SStkrjWASZQntN6psVl5dIog/sEAzsQ4ixgSwrQiyoWSdQZZjZT0q7AGXh5pO3M7KW84wqCoHXSrEdJRECZtEqjrAOAdSRNoDynlZXWCfJpvY7QH89Z7dLB47uMEFQdknaVH5/qdj0gaQ8z+0fecQVB0H5akNYGNEnre8Dakl5ZBZ57Gbbv03LOqSV6AaORluro6j5JfbqiiEBHX0BQA5jZH4E9gOslHZh3PEEQzGeUpOclfSjpckn9JS0maaykd9P9YyWtUDpA0v6SJkr6BO+jtZaZXWRm3wcuxDcRr9ILvvGar8xbIIb3lF8aWBRfd17aPDXbH+uzKLws6SNJD0kakGLYQdJzkqZLui+VXCrF97qk4yU9C8yU1EfSUEk3p9f0mqSj2vMmhaDqnDRy+i/gRElnpfbbQRDky17A1ngFiDXxjgW9gMuBYcBKuCsuhPm55fOBbcxsMLAZ8HR6bCfgRLyC++Lfgrf2a2VhxD3AA8DL+GaoG4Al0mPHAk9Bn0fhXmBx4DhgXmqgeh1wNLAUMA64Iy32KPFdYFt8n9U8fNHWM8DywDeBoyVt3dY3KT6sGoCUg9oU+ArwZ0kdSZwGQdB1XGhmk83sA+B04Ltm9r6Z3Wxms1LNv9OBLTPHzMM3AQ8ws2lm9ly6/xC8gPQLQJ+T4ZOn8fIWC6IvUBqGGb4KY7l0gcuA//P7+pnZF2b2SJpm3B2408z+ljYtn4Pvqdosc+rz0+uaDYwCljKzU83sczObCFyKz+q0iRBUg2Bm7wH/DcwA7k9L0oMgyIfJmZ8nAUMlDZR0saRJkj7GBzlDJPVOG3p3Bw4Fpkm6U9JwSV/CO1ffhdEAABvzSURBVAxfJOkLYPbq8FXDK8guiG8ARwCHA8sABwMfA+/hm6FW86dNrzhsKBnvpVz3ZHx01NzrGpZe1/TSDR/pLdPKezOfEFQDkb4F7Y9XQP6XpA3zjSgIGpYVMz+vhHcDHoN3z/6ymS2CT82Dl9vrB3yAt9u5FV8c8RRwCfAFcA0+QzJgJhw/G2ZnhzXNcRS+6uI5fKrvbHxvSn/gFfdUZXPUqbh05geVXkfWhdl9S5OB18xsSOY22MxGtxLafEJQDYY5p+NTzX+TtF3eMQVBA3K4pBUkLY6PKm4ABuN5p1mSvoHnowAeAz7EOxhshi8tvwN42Mw2Ao7Ep9NmmtmnF8ItN7SSg3oc+BcwB1gYl1JvXAgHAsfAQkfA3ZJ6S/pKEuSNwLaSvpnKOY0BPgMeWcBlHgM+TgsnBqRzrS8vvdQmQlANSqrZtwNwiaQfpm9DQRD0DNfiaxUm4lNpU/EptM2Bj/DaeKU9TUfiZZHeAHbEc1Or49N9mNmtwFn4at2Pj4T7roZpeEqpWT4Gvg8shg+JlsC/sQKcDfOWgYm/hb/io7azgF4pl703cAE+G7g9Xj398+aukbocb49PQb6WjvkDvnCwTUQliQZH0sp43a77gR9GA8T6JSpJ5Eta7bYB1VUiXqV6w22n6vGlArH3EZUkglrGzF5PFZBvAMZK2t3MPso7riCoZSrq7JVklK2zNx64gq6QUXOYPY40BjiX9kmqVIsvdzlBCCoAzOyjlIv6DfCIpG3N7PWcwwqCmiBTqTxb9LWyUvlVuIy6pSdUs3HBc32gd79mHmumm2Ihq5mHoALAGyACR0g6EpfU/5jZo3nHFQRFoo29nq6mGxsUthUzexDojzQSr603mgX3gxqH94MqxMipRAgqKMPMLpA0Ebhd0pFmdkPeMQVBHjTTLXdk+vvrNE3TXUMBZNQiLp1dkJai+Y66Vxa1o24IKqjCzO6U9N+4pNbEu3PWxmqaIOgASUbrUi6j9fGNqSUZXYvLqJkZshrAJXRO3mG0hxBU0Cxm9oykTfEuvWtKOigaIAb1gKQ+NMmoNFW3Ab6Mu5Qzuh6X0Sd5xRmEoIIWMLNpkrbEE7z3Sto5lUwKgpogyWgdynNGw/EqB6Wc0Q3AUyGj4hGCClrEzGZJ+g6+ObDUAPHFvOMKgkoyMspO05VkVJqmu5GQUc0QggpaJRWFPCE1QLxf0p5m9ve84woalySjtamW0RSaZHQTLqOP84oz6BwhqKDNmNnlkl4DbpD0EzP7Q94xBfWPpN40yag0VbchXh6olDO6BZdRbDKvI0JQQbsws/skbQHcmVb4/TiNsIKg0yQZrUV5zmgjvLZcKWd0K/BkyKj+CUEF7cbMXk4r/G4BbpK0T6H3gQSFJCOj7DTdhsBbNE3T3YaPjCp7EwUNQAgq6BBm9r6krYCLgQckbW9mU/OOKygmSUZrUj5NtxHwNk3TdLfjI6OQUQCEoIJOYGafSToA+DG+wm8HM3s677iCfJHUC5dRdppuY+AdmmR0By6jD/OKMyg+IaigU6QKE2dImoA3QDzQzO7IO66gZ0gyWoNqGb1HU87oVFxGH+QVZ1CbhKCCLsHM/izpDeBWSasDv4nySPVFRkbZnFFJRqWc0S8IGQVdRAgq6DLM7F+SvoI3QFxT0lFmNifvuIL2k2S0OuU5o43xDqulabrTcBm9n1ecQX0Tggq6FDOblBogXo8vRf9OJL2LTZLRapTLaATwIU0y+iUwPmQU9CQhqKDLMbOPJe0AnEdTA8TX8o4rcNIUbDZnNAL4iKac0Zm4jKLuYpArIaigW0gNEI/MNEDcxcweyTuuRkKSaBoZbQJsCiwK/J2mnNFZ+DRdIfsBBY1NCCroVlIDxFeB2yT90MyuyzumeiTJaFWqp+k+oWma7nxguJkNyyvOIGgPIaig2zGzcZK+BdwhaQ3gF7HCr+MkGa1CuYw2AWbQNE13Dj5N907muMXwlt9BUBOEoIIewcyelfRlyhsgfpp3XEUnyWhlykW0CTCLJhmdS4WMgqAeCEEFPYaZvSXpa8CVNDVAjNxHIiOj7D6jEcBsmnJG5+EyejunMIOgxwhBBT1KaoC4O76h89G0wq/hGiAmGQ2jepruU5pyRr/BZfRWXnEGQZ6EoIIeJ7Xn+ImkV/BCs3ua2b15x9VdJBmtRLWMPqdpmu58XEbT8oozCIpGCCrIDTO7IjVAvFHSSWZ2ad4xdZYkoxUpF9FIYA5NMrqAkFEQtEoIKsgVM7tf0uaUN0D8Iu+42kJGRlkRbQJ8QZOMfovLKFqRBEE7CUEFuWNmr2QaIN4saa+iNUBMMlqB6mm6eTTljH6Xfp4ay+iDoPOEoIJCYGYfZBogPpgaIE7JI5Yko+Upl9FIfA9RaWR0Ufo5ZBQE3UQIKigMZva5pAOB42lqgPhUd14zyWgo1Tkj0SSji4GDgSkhoyDoOUJQQaFIAjgzNUC8R9L3zOz2rjq/pKyMSn/2ommf0aXAocCbIaMgyJcQVFBIzOymTAPENYBft1cYSUaVCxj60JQz+gPwA2ByyCgIikcIKigsZvZYRQPEIxbUAFHSclQvYFiIJhldBhwOvBEyCoLaIAQVFBozeyMtQ78OGCdpN6A/1QsY+tGUM7oCOIKQURDUNCGooNBIWhaX0OPAPsC7wEzgMVxGVwJHAZNCRkFQX4SggsIgaRmqc0YDaVrAcDywHr6I4edm9nBOoQZB0AOEoIJckLQ01TmjQTTljK4BfgS8VjEyuknSY/jiiaPN7NqejTwIgp4iBBV0O5KWojpnVJLReOBa4BiqZdQsZnaXpG8AY9MKv1Njei8I6o8QVNClZGSUnapbhCYZXQccC0zsjFTM7D+pPNJf8BV+34sGiEFQX4Sggg4jaUmqc0aLAk/iMroBOA6X0byuvn5FA8S/S9opGiAGQf0QggrahKQlqJbRYriMngBuxBcxvNodMloQZjZb0h7AqcC/JG1nZs/31PWDIOg+QlBBFZIWp3oBw+LAU7iMbgJOACb0pIwWRIrhJEkvA/fVewPEIGgUQlANTkZG2dHREjRN090CnEhBZNQSZnaVpNfxBognm9nFeccUBEHHCUE1EJIWo3qabkl8ZDQeuA04CXil6DJaEGb2gKQt8BV+awLH1UoDxCAIyglB1SlJRiMol9FSwNM0yeinwMu1KqMFkRogfgW4GbglNUCckXdcQRC0jxBUHSBpCC6jbM5oGVxGTwC3AyfjMmqI0URqgLg18HuaGiC+mXdcQRC0HRV9f6OkvYFdgS8DbwCTgWPM7I1cA8sJSYtSLqORNMmotNfoCRpIRi2RGhIehxeP3dHMnsw5pNxIo+qJZrZY3rEEQVuoBUH9EDgb6JvumgOsaGZv5xdVz5CRUTZvtBzVMnopZNQyknbB27R/38xuyzuePAhBBbVGLQiqPzAV33MzB7jEzI7IN6quR9IiVOeMhgLP0CSi8cCLIaOOIWkknns7jw40QKx1QlBBrVF4QcH8UdSvgbnAMDN7K+eQOkWS0caUy2h54FmaRFSS0dy84qxHJK2IN0D8F3D4ghog1iMhqKDWqBVB9Qc+Bv5mZtvmHU97kDQYl1F2AcOKVMvohZBRz5D+Ta7DGx/uZmYf5hxSjxCCCmqNmhAUgKTRwBNm9k7esSyIjIyyOaOSjLI5o5BRzkjqDZwLbA1sZ2av5hxStxOCCmqN4gvK+wbtBwwHhgDT8Q/8K8ixMKikQVRP060E/JvynNELjTSNVGtIOgzfD7ZrvTdADEEFtUZxBSWNwuu9bQMYMCDz6GxAwF3AGZg93r2haBCwEeUyGgb8h/JpuudDRrWHpG8DVwE/MrM/5R1PdxGCCmqNYgpKOhSffukP9GrhmfOAT4ExmF3UNZfWwriMsjmjVaiW0XMho/pB0vrAHbioTqnHFX4hqKDWKJ6gmuQ0sB1HzaIZSUnaDi9y+mLzl9JAqmW0Ki6jbM4oZNQASFoGb4A4ETiw3hoghqCCWqNFQaXK0Af1WOsCn9a7j/bJqcQsYEvMnpDUFzgfOAQ438yOzsgou4BhVeA5qmX0eWdfSlCbSBoAXI7nE3cq8qKc9hKCCmqNlqbPypB0iqRrOntBSQtJuknS65IsdUQtccKvYMD6wGB8Xu3stp+6P3CCpOVw0eyP56n2kfRv4D1cWusDDwP7AkPMbJSZHWpml5rZUyGnxsbMZgN7An/HGyCul3NIQdCw5FUs9iHgN8Cf59/jq/W2AXQVvmTvVWArfJ32Hq2fs5fB6KXgjXfLX9dg4ADg32b2WRfF3yKS+sQy8tolVXf/aWqA+A9Je5vZPXnHFQSNRptGUGmV04nA7pJmSHom3X+ApBckfSJpoqRDWjuXmX1uZr8xs4eAbMme/QA7Dq/30wdYC9gRH+7MjwX4HbAGbp6f4iL7CrAo9F8GpuAjpQnAZ+lUH2blJGlJSWMlTZf0gaQHJfVKj60o6RZJ70p6X9KF6f5ekk6SNEnSO5KuSrXykLRyGg1+T9IbwP9L928q6ZF0nWcqRotBwTGzq4FdgKvkudEgCHqQNgnKzO4GfgncYGaDzGzD9NA7wHbAIvgo5TxJIzoYy3DKl5JjwINA5RzL3XjC6FHgV8DBwJ/wMufTfa/Uk2a2Bj742g2obLMwJt23FF4J/ETA0ubNscAkYGW8/ND16Zj90+3reO5qEHBhxXm3BNYBtpa0PHAncBreLv1Y4GZJS7Xp3QgKgZk9CGwOHC3pvPQ7EgRBD9DmHFRzmNmdZvaqOfcD9wBbdPB0QyrvOAVfR35Axf3H40ZcD08obYUbY1FgS3gf2EHSr3GH3YQnvLPMwauCDzOzOWb2YFpW/CW8QOv/mtlMM/s0jfQA9sILjE5Mze9OAPaQlJ1OPCUdNxvYGxhnZuPMbJ6Z/Q3PjY1u/1sT5ImZTcAH6RsCt6V9cUEQdDOdEpSkbSQ9mqbJpuMfvkt28HTTs3+5EN+QcifQr+KJy2R+HlDx90/cVTsDR+M/f4H3kcpyNj4FeE+amvxxun9FYNIC8kdD8ZFViUn49GH28pMzPw8DdkvTe9PT+7M5Lsagxkj1+r4NvI03QFwh55CCoO5pj6DK1qNL6oe31D4HWMbMhgDj8DRRR3gWrxDBZcCZ+DKqdn4KzJ7ppYbmZeI1YJqkeyWdKWk3fGrvWDNbFdgeOEbSN3HBrFQxKioxFZdOiZXw6urZvlTZ92gycLWZDcncFjazM9v3koKikFZ4fh+4FnhU0iY5hxQEdU17BPU2sHJpMQGwED64eReYK2kbfLatVST1SxXKARaS1H8GXAnoT3hC6G/48Ked6Cm4F5/WexkfPV2P54V+DczEp94eAz6Q9DfgB3gzxKXS/dOAMyUtLKm/pK+mc18H/EjSKmmKp5STW9BqvWuA7SVtLal3OtfX4pt3bZOms88GjgTulrRz3jEFQb3SHkGVloS/L+lJM/sEOAq4EfgQ3ztyexvP9RI+Wloe+Cswe7Bvzr3rJDyJNApfhTAIaOPyqXnAuA98w+4MvJDrpcC1ZvZ2ygX9wsx2BE5Pz9kC2AeX77n4oo+PgB2At/BR0+7p/JcBVwMPAK/hJZaOXFAwZjYZX4R4Ii7xycD/0slp1aAYmNmt+LaICyT9r6SOzhwEQbAAilXqqIsqSXT88lqOpkoTpWoTC9FUaaJUbeKNeqzVFrSf1ADxDuBx4LAil8SKShJBrVEsQUGX1uLrmnDmS2tk5s++lLfUGE9Iq2FJU77X4b+zuxa1AWIIKqg1ukVQkk7Ep7YqedDMtmnDCTpdzbzTMbQYnoZSPsraBF/RVymtySGtxiDtjzoHn/bbLi1NLxQhqKDWKN4IqoQ0Et9rNJoF94Mah/eD6vC0XleQ8g/NSasX5VOD44E3Q1r1S6o4cQreSv7BnMMpIwQV1BrFFVQJr7zQXEfdK/PsqNsaGWllW3mMxMWaHWU9AUwJadUPkrbCV3EeY2adLrDcVYSgglqj+IKqI5K0lqc6p2VUTw+GtGqYVAV9LC6qk1MB2lwJQQW1RggqZ5K0VqBcWpvg+bVKaU0NadUO8gr9f8GrjhyQSmDlGU8IKqgpQlAFJElrRapzWnOpyGmZ2dS84gxaR94A8TK8+PBOZvZ2y0d0aywhqKCmCEHVCBlpVea05tA0yipJa1pecQbVpH+7U/AmmduZ2XM5xRGCCmqKEFQNkz74VqI6p/UZ5VODIa0CIGlvvOTWPmb21xyuH4IKaooQVJ2RpDWM6pzWp1TktMzsrbzibFQkbY7Xivy5mf2+h68dggpqihBUA5CktTLVOa1ZVOe0csuRNAqSVsNX+P0VGGNmX7RySFddNwQV1BQhqAYlI63KnNZMqnNa7+QUZt2SZHET/iVhz1R8uSeuGYIKaoYQVDCfJK1VqJ4enEF1Tiuk1Ukk9QV+h3dy3i5VwO/O64WggpoiBBW0SJLWqpRLawTwCdU5rcJW9igq6f0dg3eA3sm6sWxXCCqoNUJQQbtJH6qrUT41OALvpVWZ03ovrzhrCUk74f3LDjGzW7rpGiGooKYIQQVdQuq0vCrlU4Ob4LUTK3Na7+cVZ5GRNAKvPHEh8KuurhoSggpqjRBU0G0kaZVGWtnpwQ+pzmmFtABJK+ANEJ8EfmBmn3fhuUNQQU0Rggp6lCSt1SmX1sbAB1TntD7IK848SQ0QrwUGAbt0VQPEEFRQa4SggtxJ0lqD8pzWxsB7lOe0nmwUaaUGiL8CtgO27YoGiCGooNYIQQWFJElrTcrzWSVpZXNaTxa1xXpXIOkQ4OfAd8zsgU6eKwQV1BQhqKBmyEgruxBjY+AdynNadSUtSf8N/Ak41syu6sR5QlBBTRGCCmqaNBVWGmmVxLUR8DblOa0nzWx6XnF2Fknr4uWRrgV+1pEGiCGooNYIQQV1R5LWWpTntDYE3qI6p/VRXnG2l9QA8TZgMrB/exsghqCCWiMEFTQESVprUy2tqZRL66kiS0tSf7wB4qrAju0p7huCCmqNEFTQsGSklc1plaRVmdP6OK84K0mVPE4G9sdr+P2njceFoIKaIgQVBBkk9aFppFUS13BgCuU5rafylpakvYDzgH3N7O42PD8EFdQUIaggaIUkrXUonx4cjueCKqXV7W0zKmL7Kt624zQz+20rzw1BBTVFCCoIOkCS1rqUS2sD4A2qc1ozujmWVYE7gXuAY6oaIPriiv1mwcj7Yedt4AbgWeAKogJ9UGBCUEHQRaT+TutQntMqSSub0+pyaUkago+kPgW+a2afII0CTgC2AQwYkDlkNiDgLuAMzB7vyniCoCsIQQVBN5KkVRpplcS1PjCJ8unBpzsrrXStC4FNX4br14CTgP5ArxYOm4dLbQxmF3Xm+kHQ1YSggqCHSSJZj+qc1lzgFsqlNbOd59YlcONesMtAHyG1lVm0UVKSDFijK+oDBkFLhKCCoABIWohqaa0HTKQ8p/W0mc1q4USjgPuAgW257n3A3sCb/tdZwJa00tU3BBX0FH3yDiAIGh1JfVLfp6fS7Q/p/pK0SlOD+wDrSZpIeU4rK60T8Gm9jtA/Hb9LB4/vEOn1z+3Jawa1QUtz00HQEEgaKulmSe9Kek3SUen+cZLOzTzvBkmXpZ/3l/SwpAskfSTpRUnfzDx3UUl/lDRN0hRJp6WNwdljz5P0AXBKuu+hzPEGHATcCPwaeBfYA199tzIuruHA+cB7kl7vK304BHbeDHo9m3l9KwPnpCcvCuyOJ51m4qsnpuKNpwZBrykwenNpK0n/lDQ9xX9hkmV73tMlJN0h6WNJj6fXX/b6JB0u6RXglXTfZum5H6U/N8s8/3VJ38r8/RRJ16SfV07nO1jS1BTzmPbEGxSTEFTQ0KQK6XcAzwDLA98Ejpa0NXAgsI+kb6RNsaOAH2YO/zI+BbckXtnhFkmLp8euxHNKq+MV17fChVN57NLA6QsI79v4yGlT4DjgEmDPFOfCwENmNhL4OrDIcfDUNPjiEGAH4LPMiW4E7gZeo7S+3E9wFzAUmJFuy4P9wL31o/S6vpLek8Naeh+b4be4A5cF9ku3SnZK78O66X27ExfuEriU75S0RDuu+XW8r9hWwI+zQgtqkxBU0OiMApYys1PN7HMzmwhcCuxhZm8Bh+Ky+T+8YkN2I+47wG/MbI6Z3QC8BGwraRn8Q/5oM5tpZu/gFR/2yBw71cwuMLO5LRR9PcvMPjaz54D/APeY2cRUK/AuXHzgJY9+ezpMGQC99wP6AY9mTnQULqLFge2Bpxf8fgzYC5Yws0dTbK8DFwNbLviQctJIcRfgZDObZWbP4+9hJWeY2Qfp9W8LvGJmV6frXge8mMJtKz9P7/e/gcuB77bj2KCARA4qaHSGAUMlZVtx9AYeTD+PxZduv2RmD1UcO8XKVxlNwj0wDOgLTPOyeYB/GZyceW725wWRLQQ7u5m/L5t5DfsNgr6l/9Cf41N3JZbN/Dyw4rFKnoDlR0lj8dzXQPxzYnwb4i2xVDqmtdebvW8o/v5lmYSPFttK9nyT8D1oQQ0TI6ig0ZkMvGZmQzK3wWY2Oj1+OvACsJykym/kyytjIGAl/LN/Mj7DtmTmnIuY2XqZ53bl8tnJwOkz4PrpwHR8OV5bhg/NrUP/nu/behFfqbcIcOICnrog3sWnN1fI3LdiM8/LvgdTcdFmWQmvgQg+XZhdmbgs1WSvUfq3CGqYEFTQ6DwGfCzpeEkDJPWWtL6kUZL+CzgA2DfdLpCU/Ua/NHCUpL6SdsOrSIwzs2l42aFzJS0iqZek1SS1eZqsnVwKHHoFfDgPZs/EkzltKQq4DPA+kOkvMvt9mAN8DMyQtDbwg/YEk0ot3YIv/hiYzrFvK4eNA9aUtKekPpJ2x0U5Nj3+NLBHeq9HArs2c46fpuuth/+73dCeuIPiEYIKGpr0Ybo93oX3NeA9fJn3csBVwBFmNiVN7/0RuDwzavoXnpR/Dx9p7Wpm76fH9gUWAp4HPsTLEC3XTa/hCeD7B8MWi8OA1fFFEG1hbXyktSowBJgCWg4OxxdjfILLryMf9EfgiwbfAq4GrqN83Ubla3gf2A4YgzvzOLyVyHvpKT8FVsPfy5/jnYUruR+YAPwdOMfM7ulA3EGBiI26QdABJO0PHGRmm+cdSxnSLcCOdOzL5zzgNsy6fB+UpLOAZc2sudV8nT33yviXi76xn6q+iBFUENQXZ+DbnDrCp+n4TiNpbUnD5XwJ+B5wa1ecO2gcQlBBUE94VfIx+DqJ9lCqxddimaMskp6TNKOZ217AYDwPNRPfhnUu8Jd2xhQ0ODHFFwT1iHQoLoWoZh7ULCGoIKhXfLXbCcBoFtwPahzeD6rNI6cg6ClCUEFQ70hL4aWGhuOL9abjFY+ujI66QZEJQQVBEASFJBZJBEEQBIUkBBUEQRAUkhBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSSEFQQBEFQSEJQQRAEQSEJQQVBEASFJAQVBEEQFJIQVBAEQVBIQlBBEARBIQlBBUEQBIUkBBUEQRAUkhBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSSEFQQBEFQSP4/52m1jDWwkhQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sm = StructureModel()\n",
    "sm.add_edges_from([\n",
    "    ('base_score', 'delta_3m_score'),\n",
    "    ('base_score','delta_6m_score'),\n",
    "    ('base_score', 'delta_12m_score'),\n",
    "    ('experimental_group','delta_3m_score'),\n",
    "    ('experimental_group','delta_6m_score'),\n",
    "    ('experimental_group', 'delta_12m_score'),\n",
    "    ('delta_3m_score','delta_6m_score'),\n",
    "    ('delta_6m_score','delta_12m_score')\n",
    "])\n",
    "_ = plot_structure(sm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The base score should not provide information for the experimental group given the randomization, though this could be investigated further.\n",
    "\n",
    "Though this model is still very naive it provides a point of comparison comparison to aid model refinement.\n",
    "\n",
    "Now its time to get the relevant data for this model.\n",
    "\n",
    "In a more detailed analysis, it would be wise to QC the data at this point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# requires connection to pdsas postgresql databse\n",
    "mdsupdrs = pd.read_sql('SELECT t1.idnum, t1.timepoint, t1.mds_part3, t2.typerandomization1, t2.typerandomization2 \\\n",
    "FROM mdsupdrs AS t1 LEFT JOIN participant AS t2 ON t1.idnum = t2.idnum', conn)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C011 [4. 8.]\n",
      "C017 [4. 8.]\n",
      "C019 [4. 8.]\n",
      "C028 [4.]\n",
      "C030 [4.]\n",
      "C032 [4.]\n",
      "C048 [ 4. 11. 17.]\n",
      "C064 [4. 8.]\n",
      "C066 [4.]\n",
      "C080 [4.]\n",
      "D002 [ 4.  8. 11.]\n",
      "D012 [ 4.  8. 11.]\n",
      "D016 [ 4. 11. 17.]\n",
      "D023 [4. 8.]\n",
      "D025 [ 4.  8. 11.]\n",
      "D027 [ 4. 11.  8.]\n",
      "D028 [ 8. 11.  4.]\n",
      "D029 [ 4.  8. 11.]\n",
      "D034 [ 4.  8. 11.]\n",
      "D036 [ 4.  8. 11.]\n",
      "D037 [4.]\n",
      "D039 [ 4.  8. 11.]\n",
      "P030 [4.]\n",
      "P032 [11.  8.  4.]\n"
     ]
    }
   ],
   "source": [
    "# some timepoints dont have values\n",
    "tmp = mdsupdrs[(mdsupdrs['timepoint'] != 0)]\n",
    "for idnum in tmp['idnum'].unique():\n",
    "    if len(tmp[(tmp['idnum']==idnum)].timepoint.unique()) != 4:\n",
    "        print(idnum, tmp[(tmp['idnum']==idnum)].timepoint.unique())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "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>idnum_base</th>\n",
       "      <th>timepoint_base</th>\n",
       "      <th>mds_part3_base</th>\n",
       "      <th>typerandomization1_base</th>\n",
       "      <th>typerandomization2_base</th>\n",
       "      <th>idnum_3m</th>\n",
       "      <th>timepoint_3m</th>\n",
       "      <th>mds_part3_3m</th>\n",
       "      <th>typerandomization1_3m</th>\n",
       "      <th>typerandomization2_3m</th>\n",
       "      <th>idnum_6m</th>\n",
       "      <th>timepoint_6m</th>\n",
       "      <th>mds_part3_6m</th>\n",
       "      <th>typerandomization1_6m</th>\n",
       "      <th>typerandomization2_6m</th>\n",
       "      <th>idnum_12m</th>\n",
       "      <th>timepoint_12m</th>\n",
       "      <th>mds_part3_12m</th>\n",
       "      <th>typerandomization1_12m</th>\n",
       "      <th>typerandomization2_12m</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C002</td>\n",
       "      <td>4.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C002</td>\n",
       "      <td>8.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C002</td>\n",
       "      <td>11.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C002</td>\n",
       "      <td>17.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C003</td>\n",
       "      <td>4.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>C003</td>\n",
       "      <td>8.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>C003</td>\n",
       "      <td>11.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>C003</td>\n",
       "      <td>17.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>C005</td>\n",
       "      <td>4.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C005</td>\n",
       "      <td>8.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C005</td>\n",
       "      <td>11.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C005</td>\n",
       "      <td>17.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C004</td>\n",
       "      <td>4.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C004</td>\n",
       "      <td>8.0</td>\n",
       "      <td>45.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C004</td>\n",
       "      <td>11.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C004</td>\n",
       "      <td>17.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>C007</td>\n",
       "      <td>4.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C007</td>\n",
       "      <td>8.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C007</td>\n",
       "      <td>11.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C007</td>\n",
       "      <td>17.0</td>\n",
       "      <td>28.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  idnum_base  timepoint_base  mds_part3_base  typerandomization1_base  \\\n",
       "0       C002             4.0            33.0                      2.0   \n",
       "1       C003             4.0            37.0                      3.0   \n",
       "2       C005             4.0            17.0                      1.0   \n",
       "3       C004             4.0            31.0                      1.0   \n",
       "4       C007             4.0            33.0                      2.0   \n",
       "\n",
       "   typerandomization2_base idnum_3m  timepoint_3m  mds_part3_3m  \\\n",
       "0                      NaN     C002           8.0          19.0   \n",
       "1                      1.0     C003           8.0          29.0   \n",
       "2                      NaN     C005           8.0          13.0   \n",
       "3                      NaN     C004           8.0          45.0   \n",
       "4                      NaN     C007           8.0          21.0   \n",
       "\n",
       "   typerandomization1_3m  typerandomization2_3m idnum_6m  timepoint_6m  \\\n",
       "0                    2.0                    NaN     C002          11.0   \n",
       "1                    3.0                    1.0     C003          11.0   \n",
       "2                    1.0                    NaN     C005          11.0   \n",
       "3                    1.0                    NaN     C004          11.0   \n",
       "4                    2.0                    NaN     C007          11.0   \n",
       "\n",
       "   mds_part3_6m  typerandomization1_6m  typerandomization2_6m idnum_12m  \\\n",
       "0          24.0                    2.0                    NaN      C002   \n",
       "1          31.0                    3.0                    1.0      C003   \n",
       "2          17.0                    1.0                    NaN      C005   \n",
       "3          33.0                    1.0                    NaN      C004   \n",
       "4          25.0                    2.0                    NaN      C007   \n",
       "\n",
       "   timepoint_12m  mds_part3_12m  typerandomization1_12m  \\\n",
       "0           17.0           13.0                     2.0   \n",
       "1           17.0           30.0                     3.0   \n",
       "2           17.0            8.0                     1.0   \n",
       "3           17.0           37.0                     1.0   \n",
       "4           17.0           28.0                     2.0   \n",
       "\n",
       "   typerandomization2_12m  \n",
       "0                     NaN  \n",
       "1                     1.0  \n",
       "2                     NaN  \n",
       "3                     NaN  \n",
       "4                     NaN  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Need to fill the missing timepoint mds_total with NA via full join\n",
    "# this is easier than pivoting given different index sizes\n",
    "# also this will drop timepoint 0 which is an unknown \"special timepoint\" \n",
    "# 1=Screening visit 1, 2=Screening visit 2,4=Baseline, 8=3 months, 11=6 months, 17=12 months, 0=special time point\n",
    "data4 = mdsupdrs[(mdsupdrs['timepoint'] == 4)]\n",
    "data4 = data4.add_suffix('_base')\n",
    "data8 = mdsupdrs[(mdsupdrs['timepoint'] == 8)]\n",
    "data8 = data8.add_suffix('_3m')\n",
    "data11 = mdsupdrs[(mdsupdrs['timepoint'] == 11)]\n",
    "data11 = data11.add_suffix('_6m')\n",
    "data17 = mdsupdrs[(mdsupdrs['timepoint'] == 17)]\n",
    "data17 = data17.add_suffix('_12m')\n",
    "\n",
    "data = data4.merge(data8, left_on=\"idnum_base\", right_on=\"idnum_3m\", how=\"outer\")\n",
    "data = data.merge(data11, left_on=\"idnum_base\", right_on=\"idnum_6m\", how=\"outer\",)\n",
    "data = data.merge(data17, left_on=\"idnum_base\", right_on=\"idnum_12m\", how=\"outer\",)\n",
    "\n",
    "data.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "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>idnum_base</th>\n",
       "      <th>typerandomization1_base</th>\n",
       "      <th>typerandomization2_base</th>\n",
       "      <th>mds_part3_base</th>\n",
       "      <th>delta_3m</th>\n",
       "      <th>delta_6m</th>\n",
       "      <th>delta_12m</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C002</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>33.0</td>\n",
       "      <td>-14.0</td>\n",
       "      <td>-9.0</td>\n",
       "      <td>-20.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C003</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>-8.0</td>\n",
       "      <td>-6.0</td>\n",
       "      <td>-7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>C005</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>17.0</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C004</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>31.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>C007</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>33.0</td>\n",
       "      <td>-12.0</td>\n",
       "      <td>-8.0</td>\n",
       "      <td>-5.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  idnum_base  typerandomization1_base  typerandomization2_base  \\\n",
       "0       C002                      2.0                      NaN   \n",
       "1       C003                      3.0                      1.0   \n",
       "2       C005                      1.0                      NaN   \n",
       "3       C004                      1.0                      NaN   \n",
       "4       C007                      2.0                      NaN   \n",
       "\n",
       "   mds_part3_base  delta_3m  delta_6m  delta_12m  \n",
       "0            33.0     -14.0      -9.0      -20.0  \n",
       "1            37.0      -8.0      -6.0       -7.0  \n",
       "2            17.0      -4.0       0.0       -9.0  \n",
       "3            31.0      14.0       2.0        6.0  \n",
       "4            33.0     -12.0      -8.0       -5.0  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# make the numeric delt timepoint columns\n",
    "\n",
    "# 3 months - base\n",
    "data['delta_3m'] = data['mds_part3_3m'] - data['mds_part3_base']\n",
    "# 6 months - base\n",
    "data['delta_6m'] = data['mds_part3_6m'] - data['mds_part3_base']\n",
    "# 12 monts - base\n",
    "data['delta_12m'] = data['mds_part3_12m'] - data['mds_part3_base']\n",
    "\n",
    "# Keep relevant columns\n",
    "data = data.filter(items=['idnum_base','typerandomization1_base','typerandomization2_base','mds_part3_base','delta_3m','delta_6m','delta_12m'])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data needs to be discretized appropriately for the model. \n",
    "\n",
    "As reported in the study publication, the within-person minimal clinically important MDS-UPDRS Motor Examination (ME) score is 3.2. Therefore, the delta from the base at each timepoint is discretized as:\n",
    "    < -3.2           = improvement from baseline\n",
    "    > -3.2 and < 3.2 = no change from baseline\n",
    "    > 3.2            = worsening from baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "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>typerandomization1_base</th>\n",
       "      <th>typerandomization2_base</th>\n",
       "      <th>mds_part3_base</th>\n",
       "      <th>delta_3m</th>\n",
       "      <th>delta_6m</th>\n",
       "      <th>delta_12m</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>idnum_base</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C002</th>\n",
       "      <td>2.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C003</th>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C005</th>\n",
       "      <td>1.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C004</th>\n",
       "      <td>1.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C007</th>\n",
       "      <td>2.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            typerandomization1_base  typerandomization2_base  mds_part3_base  \\\n",
       "idnum_base                                                                     \n",
       "C002                            2.0                   8888.0            33.0   \n",
       "C003                            3.0                      1.0            37.0   \n",
       "C005                            1.0                   8888.0            17.0   \n",
       "C004                            1.0                   8888.0            31.0   \n",
       "C007                            2.0                   8888.0            33.0   \n",
       "\n",
       "            delta_3m  delta_6m  delta_12m  \n",
       "idnum_base                                 \n",
       "C002               0         0          0  \n",
       "C003               0         0          0  \n",
       "C005               0         1          0  \n",
       "C004               2         1          2  \n",
       "C007               0         0          0  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discretized_data = data.copy()\n",
    "discretized_data = discretized_data.set_index('idnum_base')\n",
    "\n",
    "# NA with non domain value of 8888\n",
    "discretized_data = discretized_data.fillna(8888)\n",
    "\n",
    "# NA becomes max numeric level\n",
    "discretized_data[\"delta_3m\"] = Discretiser(method=\"fixed\", numeric_split_points=[-3.2, 3.2 , 1000]).transform(discretized_data[\"delta_3m\"].values)\n",
    "discretized_data[\"delta_6m\"] = Discretiser(method=\"fixed\", numeric_split_points=[-3.2, 3.2 , 1000]).transform(discretized_data[\"delta_6m\"].values)\n",
    "discretized_data[\"delta_12m\"] = Discretiser(method=\"fixed\", numeric_split_points=[-3.2, 3.2 , 1000]).transform(discretized_data[\"delta_12m\"].values)\n",
    "discretized_data.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 months:  Counter({1: 56, 2: 37, 0: 27, 3: 9})\n",
      "6 months:  Counter({1: 51, 2: 40, 0: 26, 3: 12})\n",
      "12 months:  Counter({2: 51, 1: 35, 3: 22, 0: 21})\n"
     ]
    }
   ],
   "source": [
    "# see if the delta's are binned ok\n",
    "print(\"3 months: \", Counter(discretized_data.delta_3m ))\n",
    "print(\"6 months: \", Counter(discretized_data.delta_6m ))\n",
    "print(\"12 months: \", Counter(discretized_data.delta_12m ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, descritize the MDS UPDRS Part 3 baseline scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAANpUlEQVR4nO3dX4xc5XnH8e8vmDQpUAHxgizAXRKhFC6KQStKRBURCBEJUSFSIhW11KqQnAuQQKWqXG5Iq1YyUgPtRUXrBIovCBTxJyCgSSyXiCJVtGtwg4mLSKibOLj2IoKAXhAZnl7Mcbtddj3jnZld3uH7kVZzzjvv7Hkeefzz0TtzjlNVSJLa86HVLkCStDwGuCQ1ygCXpEYZ4JLUKANckhplgEtSo9b0m5DkI8BTwC918x+oqluSnAncB5wMPAtcU1W/ONLvWrt2bU1PTw9dtCR9kOzcufPVqppaON43wIG3gUuq6q0kxwJPJ/kH4A+A26vqviR/A1wL3HGkXzQ9Pc3s7OwyypekD64k/7nYeN8llOp5q9s9tvsp4BLggW58G3DVCOqUJA1ooDXwJMck2QUcBLYDPwZer6pD3ZR9wGnjKVGStJiBAryq3qmqDcDpwAXA2YtNW+y1STYlmU0yOzc3t/xKJUn/z1F9C6WqXge+D1wInJjk8Br66cArS7xma1XNVNXM1NR71uAlScvUN8CTTCU5sdv+KPBZYA/wJPDlbtpG4JFxFSlJeq9BvoWyDtiW5Bh6gX9/VT2W5IfAfUn+DHgOuHOMdUqSFugb4FX1A+C8RcZfprceLklaBV6JKUmNMsAlqVGDrIFrBUxvfnzFjrV3yxUrdixJ4+MZuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqO8G+EH0Ere+RC8+6E0Lp6BS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWpU3wBPckaSJ5PsSfJCkhu68a8l+VmSXd3PF8ZfriTpsEHuhXIIuKmqnk1yArAzyfbuudur6i/GV54kaSl9A7yq9gP7u+03k+wBTht3YZKkIzuquxEmmQbOA54BLgKuT/J7wCy9s/SfL/KaTcAmgPXr1w9Z7spZ6Tv2SdLRGvhDzCTHAw8CN1bVG8AdwCeADfTO0L++2OuqamtVzVTVzNTU1AhKliTBgAGe5Fh64X1PVT0EUFUHquqdqnoX+AZwwfjKlCQtNMi3UALcCeypqtvmja+bN+1LwO7RlydJWsoga+AXAdcAzyfZ1Y3dDFydZANQwF7gq2OpUJK0qEG+hfI0kEWeemL05UiSBuWVmJLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSo/oGeJIzkjyZZE+SF5Lc0I2fnGR7kpe6x5PGX64k6bBBzsAPATdV1dnAhcB1Sc4BNgM7quosYEe3L0laIX0DvKr2V9Wz3fabwB7gNOBKYFs3bRtw1biKlCS911GtgSeZBs4DngFOrar90At54JQlXrMpyWyS2bm5ueGqlST9r4EDPMnxwIPAjVX1xqCvq6qtVTVTVTNTU1PLqVGStIiBAjzJsfTC+56qeqgbPpBkXff8OuDgeEqUJC1mkG+hBLgT2FNVt8176lFgY7e9EXhk9OVJkpayZoA5FwHXAM8n2dWN3QxsAe5Pci3wE+Ar4ylRkrSYvgFeVU8DWeLpS0dbjiRpUF6JKUmNGmQJRWrK9ObHV+xYe7dcsWLHkhbyDFySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSovgGe5K4kB5Psnjf2tSQ/S7Kr+/nCeMuUJC00yBn43cDli4zfXlUbup8nRluWJKmfvgFeVU8Br61ALZKkozDMGvj1SX7QLbGcNLKKJEkDWW6A3wF8AtgA7Ae+vtTEJJuSzCaZnZubW+bhJEkLLSvAq+pAVb1TVe8C3wAuOMLcrVU1U1UzU1NTy61TkrTAsgI8ybp5u18Cdi81V5I0Hmv6TUhyL3AxsDbJPuAW4OIkG4AC9gJfHWONkqRF9A3wqrp6keE7x1CLJOkoeCWmJDXKAJekRvVdQpGGNb358dUuQZpInoFLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIa1TfAk9yV5GCS3fPGTk6yPclL3eNJ4y1TkrTQIGfgdwOXLxjbDOyoqrOAHd2+JGkF9Q3wqnoKeG3B8JXAtm57G3DViOuSJPWx3DXwU6tqP0D3eMpSE5NsSjKbZHZubm6Zh5MkLTT2DzGramtVzVTVzNTU1LgPJ0kfGMsN8ANJ1gF0jwdHV5IkaRDLDfBHgY3d9kbgkdGUI0ka1CBfI7wX+Gfgk0n2JbkW2AJcluQl4LJuX5K0gtb0m1BVVy/x1KUjrkWSdBS8ElOSGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1Ki+/6HD+8X05sdXuwRJel/xDFySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSo4a6lD7JXuBN4B3gUFXNjKIoSVJ/o7gXymeq6tUR/B5J0lFwCUWSGjXsGXgB30tSwN9W1daFE5JsAjYBrF+/fsjDSe8vK32XzL1brljR4+n9bdgz8Iuq6nzg88B1ST69cEJVba2qmaqamZqaGvJwkqTDhgrwqnqlezwIPAxcMIqiJEn9LTvAkxyX5ITD28DngN2jKkySdGTDrIGfCjyc5PDv+VZVfWckVUmS+lp2gFfVy8C5I6xFknQU/BqhJDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRw/6v9JJW0PTmx1f0eHu3XLGix9PR8QxckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcqvEUr6QJqEr2R6Bi5JjTLAJalRQwV4ksuTvJjkR0k2j6ooSVJ/yw7wJMcAfw18HjgHuDrJOaMqTJJ0ZMOcgV8A/KiqXq6qXwD3AVeOpixJUj/DBPhpwE/n7e/rxiRJK2CYrxFmkbF6z6RkE7Cp230ryYtDHHMlrQVeXe0ixmSSe4PJ7m9Fe8utK3UkYLL/3MitQ/X3q4sNDhPg+4Az5u2fDryycFJVbQW2DnGcVZFktqpmVruOcZjk3mCy+7O3do2jv2GWUP4VOCvJmUk+DPw28OhoypIk9bPsM/CqOpTkeuC7wDHAXVX1wsgqkyQd0VCX0lfVE8ATI6rl/aa5ZZ+jMMm9wWT3Z2/tGnl/qXrP546SpAZ4Kb0kNcoAB5LcleRgkt3zxk5Osj3JS93jSatZ43IlOSPJk0n2JHkhyQ3dePP9JflIkn9J8m9db3/SjZ+Z5Jmut7/vPmRvUpJjkjyX5LFuf5J625vk+SS7ksx2Y82/LwGSnJjkgST/3v3d+9Q4ejPAe+4GLl8wthnYUVVnATu6/RYdAm6qqrOBC4HrulseTEJ/bwOXVNW5wAbg8iQXArcCt3e9/Ry4dhVrHNYNwJ55+5PUG8BnqmrDvK/XTcL7EuCvgO9U1a8B59L7Mxx9b1XlT+9zgGlg97z9F4F13fY64MXVrnFEfT4CXDZp/QG/DDwL/Aa9iyXWdOOfAr672vUts6fTu7/olwCP0bt4biJ66+rfC6xdMNb8+xL4FeA/6D5jHGdvnoEv7dSq2g/QPZ6yyvUMLck0cB7wDBPSX7fEsAs4CGwHfgy8XlWHuikt3+LhL4E/At7t9j/G5PQGvSu3v5dkZ3fFNkzG+/LjwBzwd93y1zeTHMcYejPAPyCSHA88CNxYVW+sdj2jUlXvVNUGemerFwBnLzZtZasaXpIvAgerauf84UWmNtfbPBdV1fn07mh6XZJPr3ZBI7IGOB+4o6rOA/6bMS0FGeBLO5BkHUD3eHCV61m2JMfSC+97quqhbnhi+gOoqteB79Nb5z8xyeFrHBa9xUMDLgJ+K8leenf6vITeGfkk9AZAVb3SPR4EHqb3D/AkvC/3Afuq6plu/wF6gT7y3gzwpT0KbOy2N9JbO25OkgB3Anuq6rZ5TzXfX5KpJCd22x8FPkvvw6IngS9305rsrar+uKpOr6pperep+Meq+h0moDeAJMclOeHwNvA5YDcT8L6sqv8Cfprkk93QpcAPGUNvXsgDJLkXuJje3dAOALcA3wbuB9YDPwG+UlWvrVaNy5XkN4F/Ap7n/9ZSb6a3Dt50f0l+HdhG71YOHwLur6o/TfJxemetJwPPAb9bVW+vXqXDSXIx8IdV9cVJ6a3r4+Fudw3wrar68yQfo/H3JUCSDcA3gQ8DLwO/T/ceZYS9GeCS1CiXUCSpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmN+h+i40y9Z58YhAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# take a look at the base score distribution\n",
    "_ = plt.hist(discretized_data.mds_part3_base, bins='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15.0\n",
      "20.0\n",
      "28.0\n"
     ]
    }
   ],
   "source": [
    "# qauntiles seem appropriate\n",
    "print(np.quantile(discretized_data.mds_part3_base, .25))\n",
    "print(np.quantile(discretized_data.mds_part3_base, .50))\n",
    "print(np.quantile(discretized_data.mds_part3_base, .75))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "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>typerandomization1_base</th>\n",
       "      <th>typerandomization2_base</th>\n",
       "      <th>mds_part3_base</th>\n",
       "      <th>delta_3m</th>\n",
       "      <th>delta_6m</th>\n",
       "      <th>delta_12m</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>idnum_base</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C002</th>\n",
       "      <td>2.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C003</th>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C005</th>\n",
       "      <td>1.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C004</th>\n",
       "      <td>1.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C007</th>\n",
       "      <td>2.0</td>\n",
       "      <td>8888.0</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            typerandomization1_base  typerandomization2_base  mds_part3_base  \\\n",
       "idnum_base                                                                     \n",
       "C002                            2.0                   8888.0               3   \n",
       "C003                            3.0                      1.0               3   \n",
       "C005                            1.0                   8888.0               1   \n",
       "C004                            1.0                   8888.0               3   \n",
       "C007                            2.0                   8888.0               3   \n",
       "\n",
       "            delta_3m  delta_6m  delta_12m  \n",
       "idnum_base                                 \n",
       "C002               0         0          0  \n",
       "C003               0         0          0  \n",
       "C005               0         1          0  \n",
       "C004               2         1          2  \n",
       "C007               0         0          0  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# discretize on quantiles\n",
    "discretized_data[\"mds_part3_base\"] = Discretiser(method=\"fixed\", numeric_split_points=[15, 20 , 28]).transform(discretized_data[\"mds_part3_base\"].values)\n",
    "discretized_data.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({3: 33, 1: 29, 2: 35, 0: 32})"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check the counts for each base score level\n",
    "Counter(discretized_data.mds_part3_base)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 2 randomization events can be thought of as producing 4 different interventions over the course of the study; only high intensity, only moderate intensity, control and high intensity, and control and moderate intensity. \n",
    "\n",
    "Because there are different timepoints, the differences between the control group and the exercise groups can still be modeled, along with the data for the entire 12 months, without an additional node for randomization event 2. The drawback to this approach is that the control group in the first six months will now be modeled as 2 independent control groups. It would probably be better to begin by structuring the model with the additional node for randomization event 2, as this would be easier to interpret. However, I decided to continue with the structure as previously defined and treat the samples as being randomized into 4 different interventions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "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>mds_part3_base</th>\n",
       "      <th>delta_3m</th>\n",
       "      <th>delta_6m</th>\n",
       "      <th>delta_12m</th>\n",
       "      <th>experimental_group</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>idnum_base</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C002</th>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C003</th>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C005</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C004</th>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C007</th>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            mds_part3_base  delta_3m  delta_6m  delta_12m  experimental_group\n",
       "idnum_base                                                                   \n",
       "C002                     3         0         0          0                 2.0\n",
       "C003                     3         0         0          0                 3.0\n",
       "C005                     1         0         1          0                 1.0\n",
       "C004                     3         2         1          2                 1.0\n",
       "C007                     3         0         0          0                 2.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discretized_data.loc[discretized_data['typerandomization1_base'] == 1,'experimental_group'] = 1\n",
    "discretized_data.loc[discretized_data['typerandomization1_base'] == 2,'experimental_group'] = 2\n",
    "discretized_data.loc[discretized_data['typerandomization2_base'] == 1,'experimental_group'] = 3\n",
    "discretized_data.loc[discretized_data['typerandomization2_base'] == 2,'experimental_group'] = 4\n",
    "discretized_data = discretized_data.drop(['typerandomization1_base', 'typerandomization2_base'], axis=1)\n",
    "\n",
    "discretized_data.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "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>mds_part3_base</th>\n",
       "      <th>delta_3m</th>\n",
       "      <th>delta_6m</th>\n",
       "      <th>delta_12m</th>\n",
       "      <th>experimental_group</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>idnum_base</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C002</th>\n",
       "      <td>above-28</td>\n",
       "      <td>improvement</td>\n",
       "      <td>improvement</td>\n",
       "      <td>improvement</td>\n",
       "      <td>mod-intensity</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C003</th>\n",
       "      <td>above-28</td>\n",
       "      <td>improvement</td>\n",
       "      <td>improvement</td>\n",
       "      <td>improvement</td>\n",
       "      <td>control-and-high</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C005</th>\n",
       "      <td>between-15-and-20</td>\n",
       "      <td>improvement</td>\n",
       "      <td>no-change</td>\n",
       "      <td>improvement</td>\n",
       "      <td>high-intensity</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C004</th>\n",
       "      <td>above-28</td>\n",
       "      <td>worsening</td>\n",
       "      <td>no-change</td>\n",
       "      <td>worsening</td>\n",
       "      <td>high-intensity</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C007</th>\n",
       "      <td>above-28</td>\n",
       "      <td>improvement</td>\n",
       "      <td>improvement</td>\n",
       "      <td>improvement</td>\n",
       "      <td>mod-intensity</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               mds_part3_base     delta_3m     delta_6m    delta_12m  \\\n",
       "idnum_base                                                             \n",
       "C002                 above-28  improvement  improvement  improvement   \n",
       "C003                 above-28  improvement  improvement  improvement   \n",
       "C005        between-15-and-20  improvement    no-change  improvement   \n",
       "C004                 above-28    worsening    no-change    worsening   \n",
       "C007                 above-28  improvement  improvement  improvement   \n",
       "\n",
       "           experimental_group  \n",
       "idnum_base                     \n",
       "C002            mod-intensity  \n",
       "C003         control-and-high  \n",
       "C005           high-intensity  \n",
       "C004           high-intensity  \n",
       "C007            mod-intensity  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# add some readability to the discrete numeric values\n",
    "mds_part3_base_map = {0: \"below-15\", \n",
    "                      1: \"between-15-and-20\", \n",
    "                      2: \"between-20-and-28\", \n",
    "                      3:\"above-28\"}\n",
    "\n",
    "experimental_group_map = {1: \"high-intensity\",\n",
    "                          2: \"mod-intensity\",\n",
    "                          3: \"control-and-high\",\n",
    "                          4: \"control-and-mod\"}\n",
    "\n",
    "delta_map =  {0: \"improvement\",\n",
    "              1: \"no-change\",\n",
    "              2: \"worsening\",\n",
    "              3: \"no-measure\"}\n",
    "\n",
    "discretized_data[\"mds_part3_base\"] = discretized_data[\"mds_part3_base\"].map(mds_part3_base_map)\n",
    "discretized_data[\"experimental_group\"] = discretized_data[\"experimental_group\"].map(experimental_group_map)\n",
    "discretized_data[\"delta_3m\"] = discretized_data[\"delta_3m\"].map(delta_map)\n",
    "discretized_data[\"delta_6m\"] = discretized_data[\"delta_6m\"].map(delta_map)\n",
    "discretized_data[\"delta_12m\"] = discretized_data[\"delta_12m\"].map(delta_map)\n",
    "\n",
    "discretized_data.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'mds_part3_base': array(['above-28', 'between-15-and-20', 'between-20-and-28', 'below-15'],\n",
       "       dtype=object),\n",
       " 'delta_3m': array(['improvement', 'worsening', 'no-change', 'no-measure'],\n",
       "       dtype=object),\n",
       " 'delta_6m': array(['improvement', 'no-change', 'worsening', 'no-measure'],\n",
       "       dtype=object),\n",
       " 'delta_12m': array(['improvement', 'worsening', 'no-measure', 'no-change'],\n",
       "       dtype=object),\n",
       " 'experimental_group': array(['mod-intensity', 'control-and-high', 'high-intensity',\n",
       "        'control-and-mod', nan], dtype=object)}"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check discrete values are ok\n",
    "data_vals = {col: discretized_data[col].unique() for col in discretized_data.columns}\n",
    "data_vals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NaN values in the experimental group column come from individuals who dropout of the study. There is a group of samples that droped out and there is a couple of samples that are control samples that did not drop out and were not assigned an intervention at six months. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "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>mds_part3_base</th>\n",
       "      <th>delta_3m</th>\n",
       "      <th>delta_6m</th>\n",
       "      <th>delta_12m</th>\n",
       "      <th>experimental_group</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>idnum_base</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C028</th>\n",
       "      <td>between-15-and-20</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C032</th>\n",
       "      <td>below-15</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C064</th>\n",
       "      <td>below-15</td>\n",
       "      <td>no-change</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C080</th>\n",
       "      <td>between-20-and-28</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D005</th>\n",
       "      <td>above-28</td>\n",
       "      <td>no-change</td>\n",
       "      <td>no-change</td>\n",
       "      <td>worsening</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D027</th>\n",
       "      <td>between-20-and-28</td>\n",
       "      <td>no-change</td>\n",
       "      <td>worsening</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D028</th>\n",
       "      <td>below-15</td>\n",
       "      <td>no-change</td>\n",
       "      <td>worsening</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D039</th>\n",
       "      <td>above-28</td>\n",
       "      <td>no-change</td>\n",
       "      <td>worsening</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>P013</th>\n",
       "      <td>above-28</td>\n",
       "      <td>no-change</td>\n",
       "      <td>worsening</td>\n",
       "      <td>no-change</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>P032</th>\n",
       "      <td>between-20-and-28</td>\n",
       "      <td>worsening</td>\n",
       "      <td>worsening</td>\n",
       "      <td>no-measure</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               mds_part3_base    delta_3m    delta_6m   delta_12m  \\\n",
       "idnum_base                                                          \n",
       "C028        between-15-and-20  no-measure  no-measure  no-measure   \n",
       "C032                 below-15  no-measure  no-measure  no-measure   \n",
       "C064                 below-15   no-change  no-measure  no-measure   \n",
       "C080        between-20-and-28  no-measure  no-measure  no-measure   \n",
       "D005                 above-28   no-change   no-change   worsening   \n",
       "D027        between-20-and-28   no-change   worsening  no-measure   \n",
       "D028                 below-15   no-change   worsening  no-measure   \n",
       "D039                 above-28   no-change   worsening  no-measure   \n",
       "P013                 above-28   no-change   worsening   no-change   \n",
       "P032        between-20-and-28   worsening   worsening  no-measure   \n",
       "\n",
       "           experimental_group  \n",
       "idnum_base                     \n",
       "C028                      NaN  \n",
       "C032                      NaN  \n",
       "C064                      NaN  \n",
       "C080                      NaN  \n",
       "D005                      NaN  \n",
       "D027                      NaN  \n",
       "D028                      NaN  \n",
       "D039                      NaN  \n",
       "P013                      NaN  \n",
       "P032                      NaN  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# these are sample not in the defined experimental group\n",
    "discretized_data[discretized_data.experimental_group.isna()]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "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>idnum_base</th>\n",
       "      <th>typerandomization1_base</th>\n",
       "      <th>typerandomization2_base</th>\n",
       "      <th>mds_part3_base</th>\n",
       "      <th>delta_3m</th>\n",
       "      <th>delta_6m</th>\n",
       "      <th>delta_12m</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>C028</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>19.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>C032</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>11.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>C064</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>11.0</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>61</th>\n",
       "      <td>C080</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>27.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>D005</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>29.0</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>79</th>\n",
       "      <td>D027</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>-2.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82</th>\n",
       "      <td>D028</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>14.0</td>\n",
       "      <td>-2.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>88</th>\n",
       "      <td>D039</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>28.0</td>\n",
       "      <td>-2.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>P013</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>28.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>128</th>\n",
       "      <td>P032</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>22.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    idnum_base  typerandomization1_base  typerandomization2_base  \\\n",
       "19        C028                      3.0                      NaN   \n",
       "22        C032                      3.0                      NaN   \n",
       "49        C064                      3.0                      NaN   \n",
       "61        C080                      3.0                      NaN   \n",
       "69        D005                      3.0                      NaN   \n",
       "79        D027                      3.0                      NaN   \n",
       "82        D028                      3.0                      NaN   \n",
       "88        D039                      3.0                      NaN   \n",
       "98        P013                      3.0                      NaN   \n",
       "128       P032                      3.0                      NaN   \n",
       "\n",
       "     mds_part3_base  delta_3m  delta_6m  delta_12m  \n",
       "19             19.0       NaN       NaN        NaN  \n",
       "22             11.0       NaN       NaN        NaN  \n",
       "49             11.0      -3.0       NaN        NaN  \n",
       "61             27.0       NaN       NaN        NaN  \n",
       "69             29.0      -3.0       3.0        5.0  \n",
       "79             20.0      -2.0       5.0        NaN  \n",
       "82             14.0      -2.0       9.0        NaN  \n",
       "88             28.0      -2.0       4.0        NaN  \n",
       "98             28.0       1.0       9.0        0.0  \n",
       "128            22.0      30.0      25.0        NaN  "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# same samples in raw data\n",
    "data[data['idnum_base'].isin(discretized_data[discretized_data.experimental_group.isna()].index)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is not enough (n=2) control samples that did not drop out but were not assigned an intervention at six months so these are dropped from the data set.\n",
    "\n",
    "The remaining samples are controls at the initial randomization but dropped out before the second randomization. Thus, the samples can be included by random assignment to experimental groups 3 and 4 and their drop out will be modelled by the \"no-meaure\" level.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Order before shuffle:  ['C028', 'C032', 'C064', 'C080', 'D027', 'D028', 'D039', 'P032']\n",
      "Order after shuffle:  ['C080', 'D027', 'D039', 'P032', 'C064', 'D028', 'C028', 'C032']\n",
      "Group 3:  ['C080', 'D027', 'D039', 'P032']\n",
      "Group 4: ['C064', 'D028', 'C028', 'C032']\n"
     ]
    }
   ],
   "source": [
    "# drop\n",
    "discretized_data = discretized_data.drop(discretized_data[(discretized_data.experimental_group.isna()) & (discretized_data['delta_12m'] != 'no-measure')].index)\n",
    "\n",
    "# random assignment\n",
    "ids = discretized_data[(discretized_data.experimental_group.isna()) & (discretized_data['delta_12m'] == 'no-measure')].index.tolist()\n",
    "print(\"Order before shuffle: \", ids)\n",
    "\n",
    "# seed to reproduce!\n",
    "random.Random(42).shuffle(ids)\n",
    "print(\"Order after shuffle: \", ids)\n",
    "ids_g3 = ids[0:round(len(ids) / 2)]\n",
    "ids_g4 = ids[round(len(ids) / 2)::]\n",
    "print(\"Group 3: \", ids_g3)\n",
    "print(\"Group 4:\", ids_g4)\n",
    "\n",
    "# do the assignment\n",
    "discretized_data.loc[ids_g3,'experimental_group'] = \"control-and-high\"\n",
    "discretized_data.loc[ids_g4,'experimental_group'] = \"control-and-mod\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'mds_part3_base': array(['above-28', 'between-15-and-20', 'between-20-and-28', 'below-15'],\n",
       "       dtype=object),\n",
       " 'delta_3m': array(['improvement', 'worsening', 'no-change', 'no-measure'],\n",
       "       dtype=object),\n",
       " 'delta_6m': array(['improvement', 'no-change', 'worsening', 'no-measure'],\n",
       "       dtype=object),\n",
       " 'delta_12m': array(['improvement', 'worsening', 'no-measure', 'no-change'],\n",
       "       dtype=object),\n",
       " 'experimental_group': array(['mod-intensity', 'control-and-high', 'high-intensity',\n",
       "        'control-and-mod'], dtype=object)}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check discrete values again\n",
    "data_vals = {col: discretized_data[col].unique() for col in discretized_data.columns}\n",
    "data_vals\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The performace of the model can be evaluated with a train test split. One could perform k-fold cross validation to be more precise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Split 90% train and 10% test\n",
    "train, test = train_test_split(discretized_data, train_size=0.9, test_size=0.1, random_state=42)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dd7xcVdW/n29CCaEFQg0lEIrUQAoKKi/6vr6KNMGgIEixIEhTDC+KoDQBCyACKohKEcFABAQkiPhTQKlJCChFQEgICZ0EElIgyfr9sfZkzszc3Dr3njMz6/l85nPvnTvnnDWTm3lm77X3WjIzgiAIgqBo9Ms7gCAIgiBoixBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSSEFQQBEFQSEJQQRAEQSEJQQVBEASFJAQVBEEQFJIQVBAEQVBIQlBBEARBIQlBBUEQBIUkBBUEQRAUkhBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSSEFQQBEFQSEJQQRAEQSEJQQVBEASFZLm8AwiCoJeR1gEOA4YDg4DZwGPAlZi9lmdoQdAeMrO8YwiCoDeQdgJOBj4JGLBS5rfzAQETgHMxe7jvAwyC9glBBUEzIh0FnA8MoP2p/CXAAmAsZpf2RWhB0FlCUEHQbJTlNLALR80jJBUUjFgkEQQFRNKVkr7XicdNlfSxzB070XU5kR5/PtLoLh4XBL1GCCoImgRJp38Cfo9P63WHAcDJkraRNFHSrHS7S9I2dQw1CDpFCCoImoS1YOXBMITu/7/uB+xxHiwE9gfW9NNyC/C7+kQZBJ0nBBUEBUDSCEmTJc2RNI7MKEjSXpKmSJot6T5Jw9s4fvc34YQboP8qwA7p/iuArYFVgWHAZR2HYmNhPzObap6gFrAY2DxzrSsl/UzSBElzJf1D0nqSLkwjrqckjej+qxEETggqCHJG0grAzcBv8FHLDcCY9LuRwK+BI4HBuGNukbRi9hxmdscX4PEDgLnAo+n+dYDbgLdxWZ0ATG4/nJXw/VJImo2v8LsYOKfqcZ8FTsVHWAuB+9Op1wLGAxd0+gUIgmUQggqC/NkZWB640MzeM7PxQGlf0hHAZWb2oJktNrOrcCHsnD2BpBXlm3Ar2BPYDB8G7QZ8HLi343gGAZjZIGB14FjgkarH3GRmk8xsAXATsMDMrjazxcA4IEZQQY+JShJBkD9DgBlWuedjWvo6FDhM0nGZ360AfEjSVvio6mfAhrNA1asjJgBnAE/jG57mAdt3HM/s0jdm9o6kS4HXJG1tZq+mX72Sefz8Nn5epePLBEH7xAgqCPLnJWADScrct3H6OgOfnTsJH5n8J91/ID6Kehefghu8GO5Z4vkiwIdZY4ATcXvMBvbAS0q0w3y8DFKWfvgy9A269rSCoGfECCoI8ud+YBFwgqR7gS/g8tkUz0n1A+4E7sBX060B/NnM5kiaCjxpZvM3ku7aGD62JB3wLi6ptfH/6BPSSbZrPxbtAFMf80UOjwErA98DZgFP1vE5B0GHxAgqCHJA0gpp5d4RwEXAHOA84AG8dt6juE8GA/sA6+KiuBb4XFvnfBF+/RbMGgyMxFfuXYSvZlgjHbhP+2EtAW5/zFNW1wFv4SO2zYHdU74pCPqMKHUUBL2MpOWBbYHRwKh02w54DpiUbhOBR83snR5ebCfgb3S9kgTvwnt7w913wo+A+8xsbo9iCYIeEoIKgjqSkVFJRKPTz1OplNGUHsto2UF0qxbfJfDX43zh3wKgPz56OtfMru6FKIOgQyIHFQTdJMloG8oiKo2MplEW0bW4jPpuNGJ2Kb7eYpnVzM+hcmPTQui/CD6CL7IoLQbcDBdVEORCjKCCoBNIWg6XUXaabnvgBVxEpdHRFDObk1ecFXjh15MpL95rqx/U7Xg/qIl+iM7H9z31x1cQDjezt/oy7CAoEYIKgiqSjLamUkbDcRlVT9MVQ0btIa0NHPYf+MTj8LFd4cY1fOXgVdUddSUNBZ4H7sK3T30U2MvMnu/zuIOWJwQVtDQZGWVzRtsDL1Ipo0caQkbtIOkhYCfgajM7rJ3H7QHcnTbpHgd8G/i0md3fR6EGARCCClqIJKOtqJTRcHwqqySiSbiM3s4rzt5A0ofwZesD8UUQW5nZtPaPWnrsnvhm4ePNLKqaB31GCCpoSiT1pyyj0lTdDsBMamXU9DkWSfcBu6QfFwO/bW8U1cbxw4FbgV8C37N44wj6gFjFFzQ8SUbvozJntCNeQqgkopuAya0go2pSCSXh+642SV+7VCvPzB6TtDPwB2BLSV82s4X1jjUIssQIKmgoMjLKTtPtALxMbc5o9rLO04pIWgN4zszW6ME5BgJX45Ut9jOz1+sVXxBUE4IKCkuS0ZZUymhHvPZpdppucsioY+ohqHSefsDZwGfwFX5P1SO+IKgmBBUUgvSmV5JRaapuR+A1amU0K684G5l6CSpzvi8A3wcOMrO/1OOcQZAlBBX0OUlGW1CZMxoBvE7lptfJZvZmXnE2G/UWVDrnR/A2IKeY2S/rdd4ggBBU0MtkZJSdpivJKJszChn1Mr0hqHTeLYE/4gtRvmVmS+p5/qB1CUEFdSPJaHNqZfQmtdN0b+QVZ6vSW4JK5x4M3Ij/W3++1wrhBi1FCCroFklGm1GZMxqBN26tllGs9CoAvSmodP4VgF/glTj2NrOZvXGdoHUIQQUdkvbRbEZlzmgk3tAumzOaFDIqLr0tqHQNAd8Cjgb2MbNHeutaQfMTggoqyMgoO01XklF1zui1ZZ0nKB59IajMtT4D/Az4opnd2tvXC5qTqCTRwiQZDaNWRnMoi+iH+MgoZBR0GjO7QdI04CZJmwMXRnmkoKvECKpFSDLalMqc0UjgHSpzRpPM7NW84gx6j74cQWWuORSv4XcfcJyZvddX1w4anxBUE5JktAmVMhoFzKM2Z/RKTmEGfUwegkrXXQ34HT5j89mo+hF0lhBUg1Mlo+xtAWUZTSRk1PLkJah07eWAHwP/A+wZDRCDzhCCaiCSjIZSmTMqySi7gGGSmb2cV5xBMclTUJkYjgVOAcaY2X15xRE0BiGogpJktDG103TvUpszeimvOIPGoQiCSnHsAVwJfM3MrsszlqDYhKAKQJLRRtTKaBG1OaPY/Bh0i6IIKsVSaoD4a+DMWOEXtEUIqo9JMtqQShGNxrucVueMQkZB3SiSoAAkrYc3QHwG+LKZLcg5pKBghKB6kYyMqnNGRpWMgJnxKTLoTYomKABJKwFXAUPwBoix3y5YSgiqTiQZbUDtNB1U5YyAGSGjoK8poqBgaV3Hs4AD8QaIT+YcUlAQil9JQloHOAwYDgzCi5E+BlxJTp+2koyGUCmj0YAoi+gX6fuQURC0Q2rPcYqkZ4C7JR1kZnflHVeQP8UdQUk7AScDn8SnxFbK/HY+LoMJwLmYPdy7oWgItTmj/tRO070YMgqKSlFHUFkk7QZcD5xqZpfnHU+QL8UUlHQUcD4wAOjXziOX4HuAxmJ2aX0urdLIKDs6Wo6qfUbA9JBR0Eg0gqAAJG2BN0C8BfimmS3OOaQgJ4onqLKcBnbhqHlUSSpNwx0CPGZmU9q+lNanNme0ArU5oxdCRkGj0yiCApC0Jt4AcTZwcDRAbE3qLihJV+JTXad28Lip+NLSuzJ37gT8ja7JqcQ8YDfMJkoaiK8MGgNcamZHpyWt1TmjFanaZwRMCxkFzUgjCQqWNkC8DNgBb4A4I+eQgj6mEIskJJ0ObG4upgHdPM0A4GRJp+GVk1dN939F0qfwHFZJRlcBxxMyCoLCYmbvSvoi8E3gAUnRALHFKISgAFb20cwnaT/n1B79DPZYEfZb6AsoShjwX/gnx5BREDQQ6f/s9yU9C9wp6UtmdkvecQV9Q3dlsBRJIyRNljRH0jgyIyBJe0maImm2pPtSeZPq43cHvj0fPr0KDNgh3X8FsDU+DBqGj/M74inAvDzQd4GH8Lp1ywH9zMwkXSnpZ5ImSJor6R+S1pN0oaRZkp6SNKInr0cQBPXHzMYDewI/lzQ25ZiDJqdHgkpzxDcDvwHWBG7A8z5IGonX2ToSGIw75hZJK2bPYWZ3AOd8FKbNBR5N968D3Aa8jcvqBGByB/E8DAPW8wZ8a+Jty58BfgBMzTzss8CpwFrAQuD+dOq1gPHABV15DYIg6BvM7CFgF3xf5KWSls85pKCX6ekIamdgebyd83vpU05pT9IRwGVm9qCZLTazq3Ah7NzWiZbz1XNL2RM3jIDdgI8D93YQzIvAC76Zd0PgzvT1m3gPmhI3mdmkVPfrJmCBmV2dlrKOA2IEFQQFxcxeAD6E/9++XdKgnEMKepGeCmoItZUSpqWvQ4GxaXpvtqTZeMXuIW2daJFPxy1lAm6yNXHj3A683kEwK+G7Z/FR3OeA1fEc1PaZh2Wb9s1v4+dVOrhMEAQ5YmZzgE8BTwD3SRqWc0hBL9FTQb0EbFA1H7xx+jodONvMBmVuA5fR/8Xe8v0O88GHWWOAE3F7zAb2wE3THtvBQvPNu3Pw6uCkw86WNAmfHhghaVT1VGMQBI2DmS0ys68BPwP+IelDeccU1J+eCup+fFHC8ZKWk/Rp4P3pd5cDR0n6gJyVJe0padU2zvPKY/De4rT67l1cUmvjKxwm4PN1HbEbmHwE93M8r7QYF9ZOwDHALHxq4Ir0/anAlpK+ImkUPl0ZBEGDYGaXAF8EbpJ0UN7xBPWlR4Iys3eBTwOH42/4B+C7vzGziXge6pL0u2fT49rihnfhvdWh/0h85d5F+GqGNYBrgX06DmfJCnD7YtgbT1ttm657hpk9amYP4Av9bjWz4fiiiGvxRRUfxDt8Tgb6S7pU0hGSRqaFIEEQFBQzmwD8N3COpNNihV/zUKxSRz2oJLEAlnwLTv4J/NrMOkpXLePyWhnftZ6tOLEpPtedLX/0eJJzEDQMjVZJoqtkGiA+C3wpGiA2PsUSFHSrFt8SmH8MDLjUp/MWAy/jxSaPN7NFPQtHKwM7UllAdlPgX1QWkH3czN7rybWCoDdpdkFBRQPEDYB9owFiY1M8QUG71czPSbcsC2HhIngDWB/PYy3GRz07pl4zdQ5vqbSyRWY3oVZaT4S0gqLQCoKCpQ0QzwQOwhsgPpFzSEE3KaagAKTReD+o0gK+tvpB3Y73g5ooaXvgQVxqBnyiL5ueSVqFWmkNxaWVLUgb0gpyoVUEVULSocB5QDRAbFCKK6gS0trAYRPg+G3g5aG+0OEx4KrqjrqSJuKLI07Cyx0daWY39nnM5XhWpTw9WBLXxsA/qcxpPdHTqcgg6IhWExSApP/CGyCeZmadqZgWFIjiCyoh6WbgSjO7uZ3HbA+sZGYPpVJLfwAuBn5UlEKxSVojqMxpbYRLNzs9+GRIK6gnrSgoWNoA8bZ0OykaIDYOTSWoNo7ZELgVf9M/uqgr7zLSyk4PbkittJ4KaQXdpVUFBUsbII7HN/EfbGZzcw4p6ARNLah03CrAb/HtVWPMbFZvxFdvJK1GrbQ2wKWVzWk9GZ8Ig87QyoKCpcWtfw6MxBsgvphzSEEHNL2g0rH9gR/iNWj3MrNn6x1fXyBpdcrTgyVxDcGLwGel9VRIK6im1QUFkDbxngQcB3zKzCblHFLQDoVpWNibpDfrsZKeBv4u6TNm1lFx9MJhZm/hG5n/VrovSWskLqs9gO8A60l6lMqFGP8OaQWtTspF/yA1QLxD0hHd+dAb9A0tIagSZnaZpOeA8ZJONLPf5B1TT0nS+mu6AZBaEJSktSdwGrCupClU5rSeDmkFrYiZ/V7SNOBmSZsD5xdlIVVQpiWm+No41zb4ip5rge/2xmbeopGRVmlqcDTeF/IRaqXV9K9HKxJTfLVI2gh/L3gQOCb2KBaLlhRUOt86eDfg6cDhZja/HudtJNIbVlZao/Ai8lOozGmFtJqAEFTbpFW01+Gb/D/TKAupWoGWFVQ65wC8Lf0wPGH6SgeHND1pOW5perAkrrUoj7RK4nompNVYhKCWTVpIdT7wCXwh1X9yDimgxQWVziu86sTh+NLTf9Xz/M2ApMGUpVUS12C8PUlWWs+GtIpLCKpjJB2NLzTa38z+kXc8rU7LCypz/oOBHwOHmtkdvXGNZiIjrWxOaw3K0iqJ6z8hrWIQguocknYHrgZOMLPf5h1PKxOCqrzGh/Dd5meZ2c966zrNiqS1qM1pDcKnB7M5rZBWDoSgOo+k7fAqNFcDp8cKv3wIQdVeZxjeS+pPwNhYht0zkrSyU4OjgNXxkVa1tBrjj7FBCUF1DUnr4gupnge+GA0Q+54QVNvXGgTcACwEPmdmc3r7mq2EvEJ9tbRWpTan9VxIq36EoLpOaoB4Bd6FYF8zezXnkFqKENSyr7c8cAmwM7544oW+uG6rkpb9j6JSXKtQmc+aBDwf0uoeIajukRogngF8Hl/h93jOIbUMIaj2ryngBOAbwH5m9nBfXTuokFZ2IcZAytIq3UJanSAE1TMkHYIvRf+8md2ZdzytQAiqc9feB/gVcJSZ/b6vrx+USXmBrLRG4d2Wq3NaU0NalYSgeo6kXfHp/9PN7NK842l2QlCdv36pAeJPgR/Em19xkLQetTmtAdROD05r5X+3EFR9SLX7bgMmACfGQqreIwTVtRg2wJeeTsFHU4VsgBiApPWpzWmtQK20XmgVaYWg6kd6LccD8/CFVNEAsRcIQXU9jlWAa/Cl0mPM7M084wk6T0Za2ZzW8lRODU4EpjejtEJQ9SUtpPo5/rcUDRB7gRBU92LpD3wf2Adf1fNMziEF3UTSEGql1Z/KUdYkmkBaIaj6kxZSnQh8jWiAWHdaqh9UvUhzzv8n6RngXkkHmNndeccVdB0zmwnMxKduS284JWmNAr6Mf0ruJ6laWi82urSCnpH+/X+UaYD4FTO7Ke+4moUQVA8ws1+kBojXSzrJzK7KO6agZ6Q3nBnpdgssldYGlKX1lfSVJK2suGaEtFoPM7tJ0gvAH9IiivPi76DnhKB6iJndJekjwG2StgS+E3Xmmov0RvNiuv0BKqRVmho8Mn1vkqpzWjPjzar5MbNJknbGV/htKenoaIDYMyIHVSdS+Z6bgJfwiugt1wCx1UnS2pDanNZiKkdZE9PUYl/HFzmoPiAtpLoO31S+fzRA7D4hqDqSGiD+CtgcT5i+nHNIQc4kaW1E7ebiRVQtxOhtaYWg+o60kOpHwB74Qqpncw6pIYkpvjpiZgskfR5vePaApL3N7J95xxXkR5raeyHdboKl0tqYsqyOBUZJepeqfVpm9lIecQc9Iy2k+oakp4G/S/qMmd2bd1yNRgiqzqQ3pDPTCr+/SDrMzCbkHVdQHNLfyLR0uxEqpFUaZR0HjJa0kKp9WjEybxzM7NK0kOr3kr5hZtfkHVMjEYLqJczsOknT8D/Ms83skrxjCopLlbR+D0ulNZTy9ODx+EhrAbXTgyGtgmJmd0r6KHCrpPcBp8VCqs4ROaheRtKmeAPEu4BvmNminEMKGpgkrU2orT04j1ppvVJ1bOSgciRV578Zn+79Qiyk6pgQVB+QGiBejyfGDzSzt3MOKWgikrQ2pVZac6nMaT0LPBSCyo+0kOoK/EPGvtUfIoJKQlB9RKrbdTHwQXxVTzRADHqNjLSyKwdHA6vhe7mW9tOKLrF9S/q3OR04FK/h9698IyouIag+JP1hfh2v3bWfmT2Uc0hBC5Gm+KYCR1AW10hgDrXTg6/lFGbLIOlg4MfAIWb2p7zjKSIhqBxIDRB/CRxtZuPzjidoDdrKQaUPTZtROcoaCbxF7ZL31/s86CZH0ofxBohnmtnP846naMQqvhwws1skfRy4RdIWwPejFE6QB+nv7tl0GwcgqR+V0voWMFLSbDLVMHBpvZFH3M2Cmf09Seq2tMJvbDRALBMjqBxJrR5uBf4JfCUaIAa9SU9W8WWklc1pjQRmUblPK6TVDdJCqvHAfOAgM5uTc0iFIASVM5JWxhsgroE3QIz/3EGvUO9l5klam1O5cnAE8Ca1Oa1o7NkBaSHVT4EP4AuppuccUu6EoApA+o/+fWBf/A/z6ZxDCpqQvtgHlf6Wt6AypzUCeJ3KnNbkkFYtKSf4DeAEfBn6xJxDypXIQRWAtKv8pFS3q9QA8W85hxUEXSb9Lf873a6FpdLakrK0vguMkPQalTmtya1e+TvlBM+X9B9ggqQjzezGvOPKixBUgTCzX0p6Hhgn6VtmdkXeMQVBT0nSeirdfgtLq31vQXlq8DRcWq9SmdNqSWmZ2c2ZBohbAD9sxYVUIaiCYWZ/kbQb5QaIp0TdrqDZSCvVStK6BpZKqzTSGg3sDewo6RUqc1qTzWx2HnH3JWY2OTVAvBVvgPjVVltIFTmogiJpLbw9wyt4A8R5OYcUNDiNWIsvSet9VOa0dsQbg1bntN7KK87eJDVA/C2wKr6QqmVGlCGoAiNpRXxD7/vwBojRGyjoNo0oqLZI0tqKSmntAMykUlqPNIu00nP+AT6q3LNVGiDGFF+BMbOFkg4FTqXcAPGxvOMKgjxJ04OPp9vVUCGtUk5rP2AHSTOpzWk1XLHm9JxPzDRA/KyZ3ZN3XL1NCKrgpMToWakB4l2SDjez2/OOKwiKRJW0rgKQtBzlkdZoYAwwXNIMKnNajzSKtMzsF2kh1XhJJ5rZ1XnH1JvEFF8DIWkXvJnduWZ2cd7xBI1Fs0zx9YQkra2pnB4cDkyndnqwsNUcJG0N3AZcB3y3WRdShaAajNQA8Tbgr8DXowFi0FlCUG2TpLUNldLaHpdWdnqwUNKStDbeAHEGcFgzNkAMQTUgklbHGyAuAQ5olOmJIF9CUJ0nlR3amsrag9vj3XCrpwfn5hjnAOBXeJ3ETzVbA8QQVIOSPvVdBOyKl0ealnNIQcEJQfWMJK3SSKskru2AaVRKa0pfSiuVR/oucDhN1gAxFkk0KGa2SNIxwPHA/ZL2M7MH844rCJoVM3sPeDTdfg1LpbUt5VHWwcB2aSFDNqc1xcze6aW4DDgjLaT6f5IONbM7ah4orQMchufcBgGzgceAKylog8oYQTUBkvbG/8McY2bX5x1PUExiBNU3ZKSVnR7cDnieypxW3aUl6UN4247vmdlP0507AScDnwQMWClzyHxAwATgXMwermc8PSUE1SRI2hG4BbgMOKcV63YF7ROCyg9JK1ArrW2B56idHuxR1RhJw/CFVH9+D55aDs4DBgD92jlsCbAAGIvZpT25fj0JQTURktbH63Y9jjdAXJhzSEGBCEEVA0lTgS8D9+Ajq2xOaxvgP1RK69GuSkvSoFPgvlNhywHQvwuHzqObkpJ0JfCimZ3a1WOXRXtGDRqMVAppN7xm158lDc45pCAIloGZvWtmk83scjM70sxG441LDwfuw1cNXgK8LukxSVdIOlbSzpJWaufUGGxxFgztjJwEZOomDfwXXLCq9ISkWel2l6Rtuv9Mu08skmgyzOwdSfsD5+LlkfYys3/nHVcQBB2TZj1KOSpgaU3O7SiPsr4AbC3pWSpzWo9m9kKdLJ/WWyaLaFsAG8CK98G04T4F2Q84BvgdvriiT4kRVBNiZkvM7Jt4l957JH0075iCoJmQNFXS/6WRzTuSfiVpXUkTJM1Jo4410mMPkTRN0huSTqk6z/slTZT0tqRXJF1QfS0zW2hmk8zsMuAcYCRwIrA+Xr7pEOBnwBuSnl5Bem0Q7Lc+9DsWyPbnEN5Tfot0+690/w7AKsA4YA3otz18xGCtdMhiYPNOvjRrSfpzeg3uljQ081x/Iml6eq6TJO3a0clCUE2Mmf0KOBD4naQv5h1PEDQZY4D/xXtY7Y2vhPs2/sbeDzg+TY39HJfIEGAwsGHmHD8BfmJmq+GbbTu7CndXYBPgw+m4k4A1gbNOhH/NhMX3A3/BzZXlZuBB4Ak8CQa+bn4ucED5YTbANyUvAC7GxdgZDgbOwl+DKaQGlYmH8VYpa+Ldlm9IG42XSQiqyTGzv+IflL4t6Qfy9ttBEPSci83sFTObAdwLPGhmj6RpupuAEcD+wG1mdk+6/zv4irkS7wGbS1rLzOaa2QOdvPYZZvaOmf0TuAL4nJktMLPfnAMvDoT+mwBHAndXHXgyboh2k1iw0gK4AVgdOBZ4pJNx/THzXE8BdpG0EYCZXWNmb5jZIjM7H1gRbyW0TOLNqgVIOaidgV3wTy0Dcw4pCJqBbFmh+W38vAo+appeujPte3oj87gv4SOwpyQ9LGmvTl57eub7aek6SNpuF9hjPWA1fDj3etWBG3XyAsCgFO+lwNXyjb6djitV03gzE9tYSU9KekvSbFx+a7V3shBUi2Bmr+PTEXOBu9OS9CAIepeXyDghfThcurrWzJ4xs88B6+ANCcdLWrkT591I0gBJ7wf2AYZJmgI8OhgGPgO8jc/LVW8kUudjn52+9gMGAht0Jq6l1/FOwGsCM1O+6ZvAZ4E1zGwQ8FZH4YSgWog07D6cNA0taYd8IwqCpmc8sJekD6fNumeSed+V9HlJa6d2GSUhLK4+iaQV5RUhDk533YePTq7Gp/DvBL4CTF4B7lkZ5j+FJ786Yl18t3CJPwOTYMFC+Kek1YALgFnAk5043R6Z53oWPu05Hd/6sgh4DVhO0nfxQV67hKBaDHPOxlcB/bkLUwpBEHQRM3scX6Z9LT6amgW8mHnI7sDjkubiCyYOBEzSaElHSrpc0uR03C/xBXfgFWNm4zX1TjOz48zsIWDsLTBsNVjpCCoWPSyT0/ECfYPwFRqzgYNgxZXgNHzT8ObA7ma2oBOnuxY/7k3KtQkB/oQvInkan5JcQOU0ZZtEJYkWRtLOwI341MJFUR6puYlKEsWjjT1Oo/A2H89SW01ivqRN8Jp+y7fbC066EfgU3RuELAFuxmxMN46tK7FRt4UxswckfRCv27WlpK9FA8Qg6B3StNf2VDZG3BofpZREdAXdKG3UBucCn8BzR11lQTo+d0JQLY6ZTU0VkMcBt0k6wMzeyjuuIGhkkoxKdfZKMsrW2ZsEXEkbMpJ0sKTL2jjtNGDPTgVg9jDSWOB8uiapUi2+iW39UtLjwNA2fnWkmf22jft7REzxBcDSBogXAh8F9jSzqflGFNSbmOLrHTKVyrNFX/iJjBYAABitSURBVNuqVP5ob/WEaie4o3BJNWQ18xhBBYA3QASOlXQccJ+kT3dh02AQtASd7PX0G3qxQWGXMLsUaSK+P3cPlt0P6na8H1SbI6e8CEEFFZjZxZKeA26RdJyZjcs7piDIA9V2yx2dfp5KeZruGooio2Xh0hmDtDZtd9S9qqgddUNQQQ1m9kdJ/4tLaku8O2djzAUHQTdIMtqGShlth+d9SjK6FpfR3Lzi7BEuofPyDqMrhKCCNjGzR9My9FvwFX5fjgaIQTOQ8q0lGZWm6rbHi6OWcka/w2U0J684gxBU0A5m9pKk3fDd6ndJ2i+VTAqChiDJaGsqc0bD8U2ipZzROOCRkFHxCEEF7WJm8yR9FjibcgPEp/KOKwiqycgoO01XklFpmu56QkYNQwgq6JBUJ+xkSU/jhWYPMrO/5B1X0LokGW1FrYxmUJbReFxGb+cVZ9AzQlBBpzGzKyQ9D4yTdIqZ/TLvmILmR1J/yjIqTdXtAMyknDO6EZdRbDJvIkJQQZcws7+l0vl/TCv8vpVGWEHQY5KM3kdlzmhHvNBqKWd0EzA5ZNT8hKCCLmNmT2cKzY6XdEih94EEhSQjo+w03Q7Ay5Sn6W7GR0azl3WeoHkJQQXdwszekPRxvOz/PZL2NrOZeccVFJMkoy2pnKbbEe9CW5qmuwUfGYWMAiAEFfQAM1so6QvAt/AVfvuY2ZS84wryRVI/XEbZaboRwKuUZXQrLqNZecUZFJ8QVNAjUoWJcyU9izdA/KKZ3Zp3XEHfkGS0BbUyep1yzuhMXEZv5hVn0JiEoIK6YGY3SHoBuEnS5sCFUR6pucjIKJszKsmolDM6i5BRUCdCUEHdMLMHJe1CuQHi8Wb2Xt5xBV0nyWhzKnNGI/BW3qVpuu/hMnojrziD5iYEFdQVM5uWGiD+Dl+K/tlIehebJKPNqJTRSGAWZRmdA0wKGQV9SQgqqDtm9rakfYAf472l9jSz5/OOK3DSFGw2ZzQSeItyzuj7uIyi7mKQKyGooFdIDRCPyzRAHGNm9+UdVyshSZRHRqOAnYHVgb9Qzhn9AJ+mK2Q/oKC1CUEFvUpqgPgf4GZJXzOz6/KOqRlJMhpG7TTdHMrTdBcBw81saF5xBkFXCEEFvY6Z3S7pY8CtkrYAzooVft0nyWhTKmU0CphLeZruPHya7tXMcWvgLb+DoCEIQQV9gpk9JukDVDZAXJB3XEUnyWgTKkU0CphHWUbnUyWjIGgGQlBBn2FmL0v6CHAV5QaIkftIZGSU3Wc0EphPOWf0Y1xGr+QUZhD0GSGooE9JDRAPwDd0PpBW+LVcA8Qko6HUTtMtoJwzuhCX0ct5xRkEeRKCCvqc1J7jFEnP4IVmDzKzu/KOq7dIMtqYWhm9S3ma7iJcRi/lFWcQFI0QVJAbZnZlaoB4vaRTzezyvGPqKUlGG1EpotHAe5RldDEhoyDokBBUkCtmdrekD1PZAHFx3nF1hoyMsiIaBSymLKOf4jKKViRB0EVCUEHumNkzmQaIv5d0cNEaICYZbUjtNN0Syjmjn6XvZ8Yy+iDoOSGooBCY2ZuZBoj3pgaIM/KIJcloAyplNBrfQ1QaGV2avg8ZBUEvEYIKCoOZvSvpi8A3KTdAfKQ3r5lkNITanJEoy+gy4CvAjJBREPQdIaigUCQBfD81QLxT0pfM7JZ6nV9SVkalr/0o7zO6HDgKeDFkFAT5EoIKComZjc80QNwCuKCrwkgyql7AsBzlnNEvga8C00NGQVA8QlBBYTGzh6oaIB67rAaIktandgHDCpRl9GvgGOCFkFEQNAYhqKDQmNkLaRn6dcDtkj4DDKB2AcOKlHNGVwLHEjIKgoYmBBUUGknr4RJ6GDgEeA14B3gIl9FVwPHAtJBREDQXIaigMEhal9qc0UDKCxi+CWyLL2I4w8z+kVOoQRD0ASGoIBckrUNtzmgVyjmja4ATgOerRkbjJT2EL574upld27eRB0HQV4Sggl5H0trU5oxKMpoEXAt8g1oZtYmZTZD038BtaYXfmTG9FwTNRwgqqCsZGWWn6lajLKPrgBOB53oiFTP7VyqP9Ad8hd+XogFiEDQXIaig20hai9qc0erAZFxG44CTcBktqff1qxog/kXSvtEAMQiahxBU0CkkDaZWRmvgMpoIXI8vYvhPb8hoWZjZfEkHAmcCD0ray8ye6KvrB0HQe4SgghokrUntAoY1gUdwGY0HTgae7UsZLYsUw6mSngb+1uwNEIOgVQhBtTgZGWVHR4MpT9PdCHybgsioPczsaklT8QaIp5nZZXnHFARB9wlBtRCS1qB2mm4tfGQ0CbgZOBV4pugyWhZmdo+kXfEVflsCJzVKA8QgCCoJQTUpSUYjqZTR2sAUyjL6DvB0o8poWaQGiLsAvwduTA0Q5+YdVxAEXSME1QRIGoTLKJszWheX0UTgFuA0XEYtMZpIDRA/AfyccgPEF/OOKwiCzqOi72+U9Hlgf+ADwAvAdOAbZvZCroHlhKTVqZTRaMoyKu01mkgLyag9UkPCk/DisZ8ys8k5h5QbaVT9nJmtkXcsQdAZGkFQXwN+BCyf7noP2MjMXskvqr4hI6Ns3mh9amX075BR+0gag7dpP8LMbs47njwIQQWNRiMIagAwE99z8x7wCzM7Nt+o6o+k1ajNGQ0BHqUsoknAUyGj7iFpNJ57+zHdaIDY6ISggkaj8IKCpaOoC4BFwFAzeznnkHpEktEIKmW0AfAYZRGVZLQorzibEUkb4Q0QHwSOWVYDxGYkBBU0Go0iqAHA28CfzWzPvOPpCpJWxWWUXcCwEbUyejJk1Dekf5Pr8MaHnzGzWTmH1CeEoIJGoyEEBSBpD2Cimb2adyzLIiOjbM6oJKNszihklDOS+gPnA58A9jKz/+QcUq8TggoajeILyvsGHQYMBwYBs/E3/CvJsTCopFWonabbGPgnlTmjJ1tpGqnRkHQ0vh9s/2ZvgBiCChqN4gpK2gmv9/ZJwICVMr+dDwiYAJyL2cO9G4pWAXakUkZDgX9ROU33RMio8ZC0O3A1cIKZ/TbveHqLEFTQaBRTUNJR+PTLAKBfO49cAiwAxmJ2aX0urZVxGWVzRptSK6PHQ0bNg6TtgFtxUZ3ejCv8QlBBo1E8QZXlNLALR82jDUlJ2gsvcvpU25fSQGplNAyXUTZnFDJqASStizdAfA74YrM1QAxBBY1GXQQl6UrgRTM7tYPHTQW+vMxWCD6t9ze6JqcS84DdMJsoaXngIuBI4CIz+3pGRtkFDMOAx6mV0bvduH7QBEhaCbgCzyfuW+RFOV0lBBU0Gu1Nn/Uqkk6XdE3V3Sfj03qd5l28DtImMFDwsKRP46I5HM9THSFpAfAO8FfgYOAfwKHAIDPbycyOMrPLzeyRkFNrY2bzgYOAv+ANELfNOaQgaFlyE1QNvlrvk3Qjpg8D1wDrAat7Z9fhlEW3IvAlYOV0//rAHDObHDIK2sLMlpjZd4DvAn+V9PG8YwqCVqRbgpI0QtJkSXMkjSMz6pG0l6QpkmZLuk/S8DaO3x1vgneApLmSHgUO+yX02xpYFZ9760y3uRWAr+OS6gd2lH/yvRB4FliYnuMDZjbPzP6N5xg+lInFJB0t6Zn0fM6StJmk+yW9Lel6SSt053UKGhsz+w0wBrhanhsNgqAP6bKg0pv1zcBv8DbgN+D/iZE0Evg1nvsZjDvmFkkrZs9hZncA5wDjzGwVM9sBGL4+rHAbXjLiCuAEvK1rp2MD7ejnP8HMtsA3yX4GeDHFJ2BXPO+UZXc8L7UzXvn6F/hU4EbAdsDnuhBG0ESY2b3455+vS/px2uAbBEEf0J0R1M54ZfELzew9MxsPlPYhHQFcZmYPmtliM7sKH8Xs3InzDtoT2AxPHO0GfBy4t4vB9YfVJH1M0gXAA8B4POENcDr+nK+oOuwHZva2mT2Or+C708yeM7O38L1WI7oYRtBEmNmzwC7ADsDNaV9cEAS9THcENQSYUbVPZFr6OhQYm6b3ZkuajY9ChnTivLMn4CZbEy8ZcTvweheDe8T7Rv0Jn/kbBiwGXpB0LL4wYk8zW1h1WLZ1x/w2fo43pBYn1evbHf/buFfShjmHFARNT3cE9RKwQZouK1EaoUwHzjazQZnbQDO7ro3zVKxvfwseHwOciL8DzAb2qH5QBxjYYrgLH7VZ+W7eAM4D7gB2kTSsKv4g6JC0qOYI4FrgAUmjcg4pCJqa7gjqfrztxfGSlkvLut+ffnc5cJSkD8hZWdKeqYhqNa8Am0jqB/B7GLcQWBvvQz8BuLOTAS3Ey0kINMFjGAE8jY+eHsT3SH0N7yv1eeAe4HVJf06n2F3SpiGtoCPM+RFwHHCHpP3yjikImpXlunqAmb2bpHQ58D18Ju7G9LuJko4ALgG2wKfH/o4LoZobcFm8Iel5Mxv5nPToZ2GHhcDewD6djOl9lOcYZ/gyc4Ct8Gm+/fBmh+dnDrkGz0eNAj4G7Et5D1Y/YH9Jc/DNu0FQg5ndJGk6npPaHDivGcsjBUGeFKvUUZ0qSXT/8lqfcqWJUrWJFShXmihVm3gh3owCWNoA8VZ8odDRRS6JFZUkgkajWIKCutbiq084S6U1OvN1eSpbakwipNWypFV91+F/s/sXtQFiCCpoNIonKKipZn4Ovmmqml2BCb0op2WHpyFUjrJG4dOl1dKaHtJqDdL+qPPwaih7paXphSIEFTQaxRQUgDQazwuVFvO11Q/qdrwfVLen9epBWlzRlrT6UTk1OAkvqlvQFz3oKanixOl4K/mubuPrVUJQQaNRXEGVkNam7Y66V+XZUbcjMtLKtvIYjYs1O8qaSO2+sqCBSbX7rgG+YWbVBZFzIwQVNBrFF1QTkaS1AbU5LaN2ejCk1cCkKui34aI6zcyW5BxSCCpoOEJQOZOktSGV0hqFdwuultbMkFbjIK/Q/wd8F8QXUiuPPOMJQQUNRQiqgCRpbURtTmsRVTktM5uZV5xBx6QGiL8GNsEbIL7S/hG9GksIKmgoQlANQkZa1Tmt9yiPskrSeimvOINa0r/d6XgtyL1SUeI84ghBBQ1FCKqBSW98G1Ob01pI5dRgSKsASPo8cAFwiJn9KYfrh6CChiIE1WQkaQ2lNqe1gKqclpm9nFecrYqkD+MtYM4ws5/38bVDUEFDEYJqAZK0NqE2pzWP2pxWbjmSVkHSZvgKvz8BY81scR9dNwQVNBQhqBYlI63qnNY71Oa0Xs0pzKYlyWI8/iHhIDOb00fXDEEFDUMIKlhKktam1E4PzqU2pxXS6iGSlgd+hrer2cvMpvfy9UJQQUMRggraJUlrGJXSGgmU2pFkc1qFrexRVNLrOxZvDbOv9WLZrhBU0GiEoIIuk95UN6NyanAk8Ba1Oa3X84qzkZC0L95j7Ugzu7GXrhGCChqKEFRQF1Jn5GFUTg2OwmsnVue03sgrziIjaSReeeIS4If1rhoSggoajRBU0GskaZVGWtnpwVnU5rRCWoCkDfEGiJOBr5rZu3U8dwgqaChCUEGfkqS1OZXSGgG8SW1O68284syT1ADxWmAVYEy9GiCGoIJGIwQV5E6S1hZU5rRGAK9TmdOa3CrSSg0QfwjsBexZjwaIIaig0QhBBYUkSWtLKvNZJWllc1qTi9pivR5IOhI4A/ismd3Tw3OFoIKGIgQVNAwZaWUXYowAXqUyp9VU0pL0v8BvgRPN7OoenCcEFTQUIaigoUlTYaWRVklcOwKvUJnTmmxms/OKs6dI2gYvj3Qt8N3uNEAMQQWNRggqaDqStN5HZU5rB+BlanNab+UVZ1dJDRBvBqYDh3e1AWIIKmg0QlBBS5CktRW10ppJpbQeKbK0JA3AGyAOAz7VleK+Iaig0QhBBS1LRlrZnFZJWtU5rbfzirOaVMnjNOBwvIbfvzp5XAgqaChCUEGQQdJylEdaJXENB2ZQmdN6JG9pSToY+DFwqJnd0YnHh6CChiIEFQQdkKS1NZXTg8PxXFC1tHq9bUZVbB/C23Z8z8x+2sFjQ1BBQxGCCoJukKS1DZXS2h54gdqc1txejmUY8EfgTuAbNQ0QfXHFYfNg9N2w3ydhHPAYcCVRgT4oMCGoIKgTqb/T1lTmtErSyua06i4tSYPwkdQC4HNmNgdpJ+Bk4JOAAStlDpkPCJgAnIvZw/WMJwjqQQgqCHqRJK3SSKskru2AaVROD07pqbTStS4Bdn4afrcFnAoMAPq1c9gSXGpjMbu0J9cPgnoTggqCPiaJZFtqc1qLgBuplNY7XTy3fgHXHwxjBvoIqbPMo5OSkmTAFvWoDxgE7RGCCoICIGkFaqW1LfAclTmtKWY2r50T7QT8DRjYmev+Dfg88KL/OA/YjQ66+oaggr5iubwDCIJWR9Jyqe/TI+n2y3R/SVqlqcFDgG0lPUdlTisrrZPxab3uMCAdP6abx3eL9PwX9eU1g8agvbnpIGgJJA2R9HtJr0l6XtLx6f7bJZ2fedw4Sb9O3x8u6R+SLpb0lqSnJP1P5rGrS/qVpJckzZD0vbQxOHvsjyW9CZye7vt75ngDvgxcD1wAvAYciK++2wQX13DgIuB1SVOXl2YNgv0+CP0eyzy/TYDz0oNXBw7Ak07v4KsnZuKNp1aBfjNgjw9LH5d0v6TZKf5Lkiy78poOlnSrpLclPZyef8Xzk3SMpGeAZ9J9H0yPfSt9/WDm8VMlfSzz8+mSrknfb5LO9xVJM1PMY7sSb1BMQlBBS5MqpN8KPApsAPwP8HVJnwC+CBwi6b/TptidgK9lDv8APgW3Fl7Z4UZJa6bfXYXnlDbHK65/HBdO9bHrAGcvI7zd8ZHTzsBJwC+Ag1KcKwN/N7PRwEeB1U6CR16CxUcC+wALMye6HrgDeJ7S+nI/wQRgCDA33TYA+6p764T0vHZJr8nR7b2ObfBT3IHrAYelWzX7ptdhm/S6/REX7mBcyn+UNLgL1/wo3lfs48C3skILGpMQVNDq7ASsbWZnmtm7ZvYccDlwoJm9DByFy+YneMWG7EbcV4ELzew9MxsH/BvYU9K6+Jv8183sHTN7Fa/4cGDm2JlmdrGZLWqn6OsPzOxtM3sc+Bdwp5k9l2oFTsDFB17y6Kdnw4yVoP9hwIrAA5kTHY+LaE1gb2DKsl+PlQ6GwWb2QIptKnAZsNuyD6kkjRTHAKeZ2TwzewJ/Das518zeTM9/T+AZM/tNuu51wFMp3M5yRnq9/wlcAXyuC8cGBSRyUEGrMxQYIinbiqM/cG/6/jZ86fa/zezvVcfOsMpVRtNwDwwFlgde8rJ5gH8YnJ55bPb7ZZEtBDu/jZ/XyzyHw1aB5Uv/od/Fp+5KrJf5fmDV76qZCBvsJN2G574G4u8TkzoRb4m10zEdPd/sfUPw1y/LNHy02Fmy55uG70ELGpgYQQWtznTgeTMblLmtamZ7pN+fDTwJrC+p+hP5BsoYCNgYf++fjs+wrZU552pmtm3msfVcPjsdOHsu/G42MBtfjteZ4UNb69C/5Pu2nsJX6q0GfHsZD10Wr+HTmxtm7tuojcdlX4OZuGizbIzXQASfLsyuTFyPWrLXKP1bBA1MCCpodR4C3pb0TUkrSeovaTtJO0n6L+ALwKHpdrGk7Cf6dYDjJS0v6TN4FYnbzewlvOzQ+ZJWk9RP0maSOj1N1kUuB466EmYtgfnv4MmczhQFXBd4A8j0F5n/BrwHvA3MlbQV8NWuBJNKLd2IL/4YmM5xaAeH3Q5sKekgSctJOgAX5W3p91OAA9NrPRrYv41zfCddb1v8321cV+IOikcIKmhp0pvp3ngX3ueB1/Fl3usDVwPHmtmMNL33K+CKzKjpQTwp/zo+0trfzN5IvzsUWAF4ApiFlyFav5eew0TgiK/ArmvCSpvjiyA6w1b4SGsYMAiYAVofjsEXY8zB5dedN/pj8UWDLwO/Aa6jct1G9XN4A9gLGIs78yS8lcjr6SHfATbDX8sz8M7C1dwNPAv8BTjPzO7sRtxBgYiNukHQDSQdDnzZzD6cdywVSDcCn6J7Hz6XADdjVvd9UJJ+AKxnZm2t5uvpuTfBP1wsH/upmosYQQVBc3Euvs2pOyxIx/cYSVtJGi7n/cCXgJvqce6gdQhBBUEz4VXJx+LrJLpCqRZfu2WOskh6XNLcNm4HA6vieah38G1Y5wN/6GJMQYsTU3xB0IxIR+FSiGrmQcMSggqCZsVXu50M7MGy+0HdjveD6vTIKQj6ihBUEDQ70tp4qaHh+GK92XjFo6uio25QZEJQQRAEQSGJRRJBEARBIQlBBUEQBIUkBBUEQRAUkhBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSSEFQQBEFQSEJQQRAEQSEJQQVBEASFJAQVBEEQFJIQVBAEQVBIQlBBEARBIQlBBUEQBIUkBBUEQRAUkhBUEARBUEhCUEEQBEEhCUEFQRAEhSQEFQRBEBSS/w9X9DZSuWRw4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Need to change the node labels but same structure as before\n",
    "sm = StructureModel()\n",
    "sm.add_edges_from([\n",
    "    ('mds_part3_base', 'delta_3m'),\n",
    "    ('mds_part3_base','delta_6m'),\n",
    "    ('mds_part3_base', 'delta_12m'),\n",
    "    ('experimental_group','delta_3m'),\n",
    "    ('experimental_group','delta_6m'),\n",
    "    ('experimental_group', 'delta_12m'),\n",
    "    ('delta_3m','delta_6m'),\n",
    "    ('delta_6m','delta_12m')\n",
    "])\n",
    "_ = plot_structure(sm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kyle/opt/anaconda3/lib/python3.7/site-packages/pandas/core/generic.py:5069: FutureWarning: Attribute 'is_copy' is deprecated and will be removed in a future version.\n",
      "  object.__getattribute__(self, name)\n",
      "/Users/kyle/opt/anaconda3/lib/python3.7/site-packages/pandas/core/generic.py:5070: FutureWarning: Attribute 'is_copy' is deprecated and will be removed in a future version.\n",
      "  return object.__setattr__(self, name, value)\n"
     ]
    }
   ],
   "source": [
    "# fit on all data so nodes can see all possible states\n",
    "bn = BayesianNetwork(sm)\n",
    "bn = bn.fit_node_states(discretized_data)\n",
    "\n",
    "# learn conditional probability of each node\n",
    "bn = bn.fit_cpds(train, method=\"BayesianEstimator\", bayes_prior=\"K2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kyle/opt/anaconda3/lib/python3.7/site-packages/sklearn/metrics/classification.py:1145: UndefinedMetricWarning: Recall and F-score are ill-defined and being set to 0.0 in labels with no true samples.\n",
      "  'recall', 'true', average, warn_for)\n"
     ]
    },
    {
     "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>precision</th>\n",
       "      <th>recall</th>\n",
       "      <th>f1-score</th>\n",
       "      <th>support</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>delta_6m_improvement</th>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>delta_6m_no-change</th>\n",
       "      <td>0.571429</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>delta_6m_no-measure</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>delta_6m_worsening</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>macro avg</th>\n",
       "      <td>0.267857</td>\n",
       "      <td>0.325000</td>\n",
       "      <td>0.291667</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>micro avg</th>\n",
       "      <td>0.461538</td>\n",
       "      <td>0.461538</td>\n",
       "      <td>0.461538</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>weighted avg</th>\n",
       "      <td>0.373626</td>\n",
       "      <td>0.461538</td>\n",
       "      <td>0.410256</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      precision    recall  f1-score  support\n",
       "delta_6m_improvement   0.500000  0.500000  0.500000        4\n",
       "delta_6m_no-change     0.571429  0.800000  0.666667        5\n",
       "delta_6m_no-measure    0.000000  0.000000  0.000000        0\n",
       "delta_6m_worsening     0.000000  0.000000  0.000000        4\n",
       "macro avg              0.267857  0.325000  0.291667       13\n",
       "micro avg              0.461538  0.461538  0.461538       13\n",
       "weighted avg           0.373626  0.461538  0.410256       13"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# we can see how well we are predicting with this model \n",
    "classification_report(bn, test, \"delta_6m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Judging from the low scores in the classification report, the model is performing poorly when classifying endpoint levels. The model has too many parameters for the number of training samples (parameters > degress of freedom). There are also only 13 samples in the test set. Since there is a small number of samples in the test set, the endpoint level no-measure at 6 months is not observed. Though the recall of no-change at 6 months (0.8) is relatively high, the precision is still low (0.57). \n",
    "\n",
    "However, supposing the model was performing well, some insights can be gained by examining\n",
    "the conditional probability distributions of the nodes.\n",
    "\n",
    "Normally, to make inferences from the conditional probability distributions, the model should be retrained on the full data set, given the model is performing well. However, I am skipping this step for the sake of brevity.\n",
    "\n",
    "Next, I compare interventions at 6 months, given that there was improvement at 3 months and a base motor score above 28."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "delta_6m\n",
       "improvement    0.555556\n",
       "no-change      0.222222\n",
       "no-measure     0.111111\n",
       "worsening      0.111111\n",
       "Name: (improvement, high-intensity, above-28), dtype: float64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# examining the conditional probabilities at the delta_6m node for the high intensity group \n",
    "tmp = bn.cpds[\"delta_6m\"]\n",
    "idx = pd.IndexSlice\n",
    "tmp.loc[idx[:], idx['improvement','high-intensity','above-28']]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this model, individuals who show some improvement at 3 months and have a base motor score above 28 have a 55% probability of improvement at 6 months with high intensity exercise. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta_6m\n",
      "improvement    0.500\n",
      "no-change      0.250\n",
      "no-measure     0.125\n",
      "worsening      0.125\n",
      "Name: (improvement, mod-intensity, above-28), dtype: float64\n",
      "delta_6m\n",
      "improvement    0.4\n",
      "no-change      0.2\n",
      "no-measure     0.2\n",
      "worsening      0.2\n",
      "Name: (improvement, control-and-high, above-28), dtype: float64\n",
      "delta_6m\n",
      "improvement    0.4\n",
      "no-change      0.2\n",
      "no-measure     0.2\n",
      "worsening      0.2\n",
      "Name: (improvement, control-and-mod, above-28), dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# delta_6m for similar individual with different treatments\n",
    "print(tmp.loc[idx[:], idx['improvement','mod-intensity','above-28']])\n",
    "print(tmp.loc[idx[:], idx['improvement','control-and-high','above-28']])\n",
    "print(tmp.loc[idx[:], idx['improvement','control-and-mod','above-28']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar indivdiuals with the moderate intensity exercise have a 50% probability of improvement. Similar individuals with no prescribed exercise (both control groups) have a 40% probability of improvement. \n",
    "\n",
    "Thus, at 6 months, high intensity exercise gives the best probability of improvement among individuals who have shown improvement after 3 months and have base scores greater than 28. Improvement is defined as less than -3.2 difference from the baseline MDS-UPDRS Motor Examination (ME)- in other words, individuals whose MDS-UPDR ME score decreased more than 3.2 points.\n",
    "\n",
    "Given the model training is sample limited, a refinement to the model could be made be decreasing the number of discrete levels at each endpoint. For example, insted of:\n",
    "    \n",
    "    \"improvement, no-change, no-measure, worsening\"\n",
    "\n",
    "A binary classification into:\n",
    "\n",
    "    \"worsening, not-worsening\"\n",
    "\n",
    "or,\n",
    "    \n",
    "    \"improving, not-improving\"\n",
    "    \n",
    "May yield better model performance. \n",
    "    \n",
    "One of the key factors missing in the model is the degree to which individuals adhere to their prescribed exercise. This could be incorporated via data from the activity monitor log (aml table).\n",
    "\n",
    "Another factor that may be important to model performance is a participants adverse events. This data could be investigate using the data in the adverse events (ae) table.*\n",
    "\n",
    "Additional counterfactuals could be examined by including phenotypic data such as sex, age, and health information (demo, gmh, gmhmed, gmnpdmed, gmpd, mnpd_med tables). One example is usage of PD specific and non-PD specific medication. Maybe, a certain medication is causing the improvement in MDS-UPDRS Motor Examination (ME) scores. \n",
    "\n",
    "Other questions we could begin to answer with this model include:\n",
    "Is high intensity exercise better for individuals with a high or low MDS UPDRS ME base score? What gives a high probability of no-measure at timepoint 12 months given a measure in timepoint 6 months; in other words, what may be causing people to drop out of the study? And conversely, what may be causing people to stick with the study?\n",
    "\n",
    "If the model is expanded to include an adherence-to-excercise node, one could perform an intervention with the causalnex InferenceEngine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# e.g.\n",
    "#ie = InferenceEngine(bn)\n",
    "#ie.do_intervention"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This would allow one to examine the conditional probability of the given endpoint when all individual adhere to their prescribed exercise to the maximum extent.\n",
    "\n",
    "Lastly, instead of looking at the MDS-UPDRS Motor Examination (ME) scores outcomes, a different outcome like depression as measured by the Beck Depression Index (bdi table) could be used. This may be important to practitioners if it was found that moderate excercise could help with depression in patients with Parkinson's disease (PD), since depression is a symptom of PD.\n",
    "\n",
    "\n",
    "*It would be interesting to develop a model for factors causing adverse events. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}