{ "cells": [ { "cell_type": "markdown", "id": "726afc96", "metadata": {}, "source": [ "# Bayesian Network Example with `pgmpy`, `pomegranate`, and `bnlearn`" ] }, { "cell_type": "markdown", "id": "74fa9e85", "metadata": {}, "source": [ "Starting with `pgmpy` (probabilistic graphical models in python), we'll do some simple Bayesian Networks. \n", "\n", "This demo uses the `TabularCPD` object to create tables of Conditional Probability Distributions (CPD). Everything should be self-explanatory and well-documented in the help, but here's some that I didn't understand at first glance:\n", "* `variable_card` = variable cardinality, i.e. the number of states the variable can take\n", "* `evidence` = list of variable *names*\n", "* `evidence_card` = cardinality of the evidence, should be a list" ] }, { "cell_type": "markdown", "id": "a42f03fa", "metadata": {}, "source": [ ">In pgmpy we define the network structure and the CPDs (conditional probability distributions) separately and then associate them with the structure. Here’s an example for defining the above model:" ] }, { "cell_type": "code", "execution_count": 1, "id": "e13176a0", "metadata": {}, "outputs": [], "source": [ "import pandas as pd" ] }, { "cell_type": "code", "execution_count": 2, "id": "7005e698", "metadata": {}, "outputs": [], "source": [ "from pgmpy.models import BayesianNetwork\n", "from pgmpy.factors.discrete import TabularCPD\n", "from pgmpy.inference import VariableElimination" ] }, { "cell_type": "markdown", "id": "ad5f6ae8", "metadata": {}, "source": [ "## Example 1\n", "\n", "This example comes from the excellent and short course by [Phillip Loick on Bayesian Statistics (Udemy)](https://www.udemy.com/course/bayesian-statistics/):" ] }, { "cell_type": "markdown", "id": "43cdda72", "metadata": {}, "source": [ "![fig2.png](fig2.png)" ] }, { "cell_type": "markdown", "id": "47a85b25", "metadata": {}, "source": [ "This is a fun model that models whether or not two people (John and Kate) will run based on the temperature and whether or not these two will meet. " ] }, { "cell_type": "code", "execution_count": 3, "id": "9aa4db63", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 1. Instantiate the model with node edges as a list\n", "model = BayesianNetwork([('T', 'J'), ('T', 'K'), ('J', 'M'), ('K', 'M')])\n", "\n", "# 2. Define the distributions\n", "cpd_t = TabularCPD(variable='T', variable_card=2, values=[[0.4], [0.6]], state_names={'T':['low', 'high']})\n", "cpd_j = TabularCPD(variable='J', variable_card=2, \n", " values=[[0.5, 0.7], \n", " [0.5, 0.3]], \n", " evidence=['T'],\n", " evidence_card=[2],\n", " state_names={'J':['yes', 'no'],\n", " 'T':['low', 'high']})\n", "cpd_k = TabularCPD(variable='K', variable_card=2, \n", " values=[[0.4, 0.75], \n", " [0.6, 0.25]], \n", " evidence=['T'],\n", " evidence_card=[2],\n", " state_names={'K':['yes', 'no'],\n", " 'T':['low', 'high']})\n", "cpd_m = TabularCPD(variable='M', variable_card = 2,\n", " values=[[.5, 0, 0, 0],\n", " [.5, 1, 1, 1]],\n", " evidence=['J', 'K'],\n", " evidence_card=[2, 2],\n", " state_names={'M': ['yes', 'no'],\n", " 'J': ['yes', 'no'],\n", " 'K': ['yes', 'no']})\n", "\n", "# 3. Add CPDs to the model\n", "model.add_cpds(cpd_t, cpd_j, cpd_k, cpd_m)\n", "# 4. Check the model validity (i.e. probabilities all sum to 1)\n", "model.check_model()" ] }, { "cell_type": "markdown", "id": "37dc4ae4", "metadata": {}, "source": [ "Let's doubleclick into the 'M' node because the evidence table looks **different** than the conditional probability table that was in the first figure:" ] }, { "cell_type": "code", "execution_count": 4, "id": "1a3b3f58", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+--------+--------+--------+--------+-------+\n", "| J | J(yes) | J(yes) | J(no) | J(no) |\n", "+--------+--------+--------+--------+-------+\n", "| K | K(yes) | K(no) | K(yes) | K(no) |\n", "+--------+--------+--------+--------+-------+\n", "| M(yes) | 0.5 | 0.0 | 0.0 | 0.0 |\n", "+--------+--------+--------+--------+-------+\n", "| M(no) | 0.5 | 1.0 | 1.0 | 1.0 |\n", "+--------+--------+--------+--------+-------+\n" ] } ], "source": [ "# Printing a CPD with it's state names defined.\n", "print(model.get_cpds('M'))" ] }, { "cell_type": "markdown", "id": "b6f7f89f", "metadata": {}, "source": [ "This is because `pgmpy` expects you to have the `variable` states on rows and then like a multi-index of `evidence` on columns. Recall the `evidence` and `evidence_card` items we called out earlier as `['J', 'K']` in the `TabularCPD` call for variable `M`, so the columns show up in the order of the `evidence` and the `state_names` for each." ] }, { "cell_type": "markdown", "id": "38661aa8", "metadata": {}, "source": [ "## Inference with Variable Elimination" ] }, { "cell_type": "markdown", "id": "35661af2", "metadata": {}, "source": [ "Let’s take an example of inference using Variable Elimination in `pgmpy`. Here we'll use our model to compute the probability distribution that the two meet, and `pgmpy` will make a pretty table: \n", "\n", "In other words, what is the probability that they meet $P(M)=?$" ] }, { "cell_type": "code", "execution_count": 5, "id": "a0efc9ac", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "53a743b5b6f3456abe871c889db86218", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/3 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "T.plot()" ] }, { "cell_type": "markdown", "id": "cec0f78e", "metadata": {}, "source": [ "A few other useful methods:\n", "\n", "* `sample(n)`: draw `n` samples from the distribution\n", "* `probability(X)`: predict the probability of X under this distribution" ] }, { "cell_type": "code", "execution_count": 18, "id": "210e7417", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['low', 'high', 'high'], dtype='" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# 1. Instantiate the BayesianNetwork model\n", "model = pm.BayesianNetwork(\"Calculating Posterior\")\n", "\n", "# 2. Add more distributions\n", "K = pm.ConditionalProbabilityTable(\n", " table=[['low', 'yes', 0.4],\n", " ['low', 'no', 0.6],\n", " ['high', 'yes', 0.75],\n", " ['high', 'no', 0.25]], parents=[T])\n", "M = pm.ConditionalProbabilityTable(\n", " table=[['yes', 'yes', 'yes', 0.5],\n", " ['yes', 'yes', 'no', 0.5],\n", " ['yes', 'no', 'yes', 0],\n", " ['yes', 'no', 'no', 1],\n", " ['no', 'yes', 'yes', 0],\n", " ['no', 'yes', 'no', 1],\n", " ['no', 'no', 'yes', 0],\n", " ['no', 'no', 'no', 1]], parents=[J, K])\n", "\n", "# 3. Define nodes in our network that follow these distributions\n", "n0 = pm.Node(T, name='Temperature')\n", "n1 = pm.Node(J, name='John')\n", "n2 = pm.Node(K, name='Kate')\n", "n3 = pm.Node(M, name='Meet')\n", "model.add_states(n0, n1, n2, n3)\n", "\n", "# 4. Define the Edges for each Node in our model\n", "model.add_edge(n0, n1)\n", "model.add_edge(n0, n2)\n", "model.add_edge(n1, n3)\n", "model.add_edge(n2, n3)\n", "\n", "# 5. Bake the model\n", "model.bake()\n", "\n", "# Optional: Plot\n", "model.plot()" ] }, { "cell_type": "markdown", "id": "949e4c19", "metadata": {}, "source": [ "I like the ability to plot (requires `matplotlib` and `pygraphviz`). Great feature!" ] }, { "cell_type": "markdown", "id": "a8214e43", "metadata": {}, "source": [ "## Inference\n", "\n", "We would like to do a simple inference - what's the probability that they meet?\n", "\n", "In other words, what is the probability that they meet $P(M=\\text{'yes'})$?\n", "\n", "In `pomegranate` you can do this by specifying all of the probabilities like such:" ] }, { "cell_type": "code", "execution_count": 22, "id": "048ed4f5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.19749999999999995" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.probability([['high', 'yes', 'yes', 'yes'],\n", " ['high', 'yes', 'no', 'yes'],\n", " ['high', 'no', 'yes', 'yes'], \n", " ['high', 'no', 'no', 'yes'],\n", " ['low', 'yes', 'yes', 'yes'],\n", " ['low', 'yes', 'no', 'yes'],\n", " ['low', 'no', 'yes', 'yes'], \n", " ['low', 'no', 'no', 'yes'], \n", " ]).sum()" ] }, { "cell_type": "markdown", "id": "38f5a580", "metadata": {}, "source": [ "There is also a `marginal` method to get the marginal distribution of each of the variables:" ] }, { "cell_type": "code", "execution_count": 23, "id": "dce7c67a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(4,)" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.marginal().shape" ] }, { "cell_type": "code", "execution_count": 24, "id": "740191e3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{\n", " \"class\" : \"Distribution\",\n", " \"dtype\" : \"str\",\n", " \"name\" : \"DiscreteDistribution\",\n", " \"parameters\" : [\n", " {\n", " \"yes\" : 0.18910000000000005,\n", " \"no\" : 0.8109000000000001\n", " }\n", " ],\n", " \"frozen\" : false\n", "}" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.marginal()[3]" ] }, { "cell_type": "markdown", "id": "844815e1", "metadata": {}, "source": [ "Interesting that the marginal distribution is different than what I had calculated manually! There must be a bug in my model somewhere..." ] }, { "cell_type": "markdown", "id": "4bf8fc22", "metadata": {}, "source": [ "## Impressions\n", "\n", "It looks like it's still too early for `pomegranate` when it comes to Bayesian Networks when compared to `pgmpy`. It felt still very much WIP to use `pomegranate` with docstrings and help/error messages seem incomplete, while `pgmpy` has extensive documentation. " ] }, { "cell_type": "markdown", "id": "f7f7e927", "metadata": {}, "source": [ "# Classic Disease Model with `pgmpy` and `bnlearn`" ] }, { "cell_type": "markdown", "id": "107cbde0", "metadata": {}, "source": [ "For this next example, we'll use the classic Disease model: \n", "\n", "Let's say that there's a disease that affects 2% of the population, and you have a 95% chance of testing positive (correctly) if you have the disease, and a 10% change of testing positive if you don't have the disease. What's the probability that you have the disease, given that you tested positive?\n", "\n", "The Conditional Probability Table would look like:\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", "\t \n", "\t \n", " \n", " \n", " \n", " \n", "\t \n", " \n", " \n", " \n", "\n", "
Test Result
PositiveNegative
DiseaseYes0.95.05
No0.900.1
" ] }, { "cell_type": "markdown", "id": "0b2a1b7c", "metadata": {}, "source": [ "An alternative to `pgmpy` is to use `bnlearn` which extends some of the features from `pgmpy` and just...does more. It can even use `pgmpy` objects like `TabularCPD`, but it looks like there are lot more convenience functions included in it to help learn from data (including the ability to handle dataframes)." ] }, { "cell_type": "code", "execution_count": 25, "id": "5c06ec3c", "metadata": {}, "outputs": [], "source": [ "# Import the library\n", "import bnlearn as bn" ] }, { "cell_type": "markdown", "id": "22f93dc1", "metadata": {}, "source": [ "`bnlearn.make_DAG` takes in an argument `DAG` which should be your list of edges, and a `CPD` argument for your conditional probability distributions that were created with `pgmpy.TabularCPD`. \n", "\n", "For fun, we'll create this model with two tests ($T1$, $T2$) so we can look at statistics if you get multiple tests." ] }, { "cell_type": "code", "execution_count": 26, "id": "505599fd", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >bayes DAG created.\n", "[bnlearn] >Add CPD: D\n", "[bnlearn] >Add CPD: T1\n", "[bnlearn] >Add CPD: T2\n", "[bnlearn] >Checking CPDs..\n", "[bnlearn] >Check for DAG structure. Correct: True\n" ] } ], "source": [ "# Define the network structure\n", "edges = [('D', 'T1'), ('D', 'T2')]\n", "\n", "d = TabularCPD(variable='D', variable_card=2, values=[[0.02], [0.98]],\n", " state_names={'D':['yes', 'no']})\n", "t1 = TabularCPD(variable='T1', variable_card=2, \n", " values=[[0.95, 0.1],\n", " [0.05, 0.9]],\n", " evidence=['D'],\n", " evidence_card=[2],\n", " state_names={'D':['yes', 'no'],\n", " 'T1':['positive', 'negative']})\n", "t2 = TabularCPD(variable='T2', variable_card=2, \n", " values=[[0.95, 0.1],\n", " [0.05, 0.9]],\n", " evidence=['D'],\n", " evidence_card=[2],\n", " state_names={'D':['yes', 'no'],\n", " 'T2':['positive', 'negative']})\n", "\n", "# Make the actual Bayesian DAG with the previously defined CPD's \n", "model = bn.make_DAG(DAG=edges, CPD=[d, t1, t2])" ] }, { "cell_type": "code", "execution_count": 27, "id": "85761655", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn]> Set node properties.\n", "[bnlearn]> Set edge properties.\n", "[bnlearn] >Plot based on Bayesian model\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAAHBCAYAAACMieH9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA83klEQVR4nO3deXxcd33v//fMaJdlWV5kybJsyZYlWcvoHAKJQ0pWh5CFOIZAVpI4TYmdBhI3ISS/XAKh3N7LbX+3hcfvUvooF7gthV6WAoXbwu1NetOWFGjInNFiLZblRbYkb5K1jaTZzu+PxKdR4kUejXQ0M6/n45HHgyjSzFvBcfTO53w/X49t27YAAAAAAJfE63YAAAAAAEhFlCkAAAAASABlCgAAAAASQJkCAAAAgARQpgAAAAAgAZQpAAAAAEgAZQoAAAAAEkCZAgAAAIAEUKYAAAAAIAGUKQAAAABIAGUKAAAAABJAmQIAAACABFCmAAAAACABlCkAAAAASABlCgAAAAASQJkCAAAAgARQpgAAAAAgAZQpAAAAAEgAZQoAAAAAEkCZAgAAAIAEUKYAAAAAIAGUKQAAAABIAGUKAAAAABJAmQIAAACABFCmAAAAACABlCkAAAAASABlCgAAAAASQJkCAAAAgARQpgAAAAAgAVluBwAAAACw8EYnQvrF6z2yug6rrbdfx0+NKhKNKcvn1aoVy9RYs15mfZXea9aqdOVyt+OmBI9t27bbIQAAAAAsjIPHTuivfvoLvfTLDsVtW7ZtKzc7WznZPnk8Htm2rUg0pplwVLZseTweXdlSo/tv+y01bal0O/6SRpkCAAAA0lAkGtO3fvLP+h8//mfF4raKCvLk8138lE88HtdEaFq2pB3XXaY9d21XQX7uwgdOQZQpAAAAIM2cPjOhp/7wW+o7ekKF+bnK8vku+TVib5aq1SuK9CfPfkwbylcvQNLURpkCAAAA0sipkXHt+f2v68TwqIoK8uXxeOb1ehOhaeXlZutPP/OwqirWJCllemCbHwAAAJAmwpGo9v6Xb+nE8JiWFxbMu0hJ0rKCPE3PRPTJ//Q/dGZsMgkp0wdlCgAAAEgT3/jhKzo8cFLLC/OT+rrLCvJ0ZnxKf/TN/yUebPt3lCkAAAAgDfQcGtS3/+5VLcvPW5DXLyrM0z/9pku/CPQsyOunIsoUAAAAkAa+8aNXJNue08a+RHg9HmX7fPrqd19iOvUmLu0FAAAAUtzJkTG9au3XsoKLT6Wsl/5aR7t/I0nyeL3Kzs1XUUmZyjc3a0PDNnkvsPkvPy9H/UOn1dk3oIbNFUnLn6ooUwAAAECKe/lXHbJtW17v3KZSq9dvkbH9btlxW+HpCZ0+2qvuf/vfOtrzurbd/qiysnPO+XUej0d23NZPX3mdMiUe8wMAAABS3mvtB+WbY5GSJK8vS3kFy5W/rFjFqyu0ybhGV+7Yo9GTx3Qg8I8X/Nrc3GwFug7PN3JaoEwBAAAAKa7z4IByc+b30NnyVWUq3VCnob62C35eTnaWjh0f1nQ4Mq/3SweUKQAAACCFzYSjGh2fVNYFzjrN1bKStZocG77g53g9HmX5vDp+anTe75fqKFMAAABACgtHovJ4vEm5oFeyNZdX8Xg8CkeiSXi/1EaZAgAAAFJYls+btFXl48PHVbB81UU/z7a1YCvYUwl/BwAAAIAUlpuTpZwsn2KxuGzb1tTUlE6fPq3+/n6NjIzM+XXGTg/pZH+3yjc3X/DzbNtWNBbTquJl842e8liNDgAAAKSoUCikYDAobyykQ0fGFAlPKScnRwUFBSpeUazCwsJzfl08FtV0aEyybYWnJnXq6H7tf/1lFa9Zr03GtRd8z2g0puJlBSouKliA7yi1UKYAAACAFDE8PCzLshQIBBQIBHT48GFt3bpV5SWVOjNlq3RV5Zzumjp1dL/+zzd/Xx6PV1m5eSpaWaba97xfGxuukNd34YowHY6opX5jsr6llOaxk/WAJQAAAICksW1bx44dc4qTZVkaGRmR3++XaZoyTVNbt25VTk6OOvuOaffnv66igrwkLaI4v/HJaX32sZ267vLGBX2fVMBkCgAAAFgC4vG4ent7Z02ePB6PTNOUYRi65557tHnz5nNOnuqr12lD2SoNnBxRQV7ugmUMR6LKy83WVWbdgr1HKqFMAQAAAC4Ih8Pat2+fU5xaW1u1atUqmaapq666So8//rjWrVs3p0mTx+PRAzvep9//6g9l2/aCTaemZsLadcfVysmmRkg85gcAAAAsiomJCbW2tjqP7HV1damqqsqZPBmGoZUrVyb8+vF4XE/8579UsPuIipflX/BzQ6GQCgoubYHERGhapSuX6y/+02PKzaFMSZQpAAAAYEGcPn3amToFAgH19/eroaHBKU9+v/+SC83FHD89qvs//RVJtvJyc875OXE7rp7uHtXV1c15ghWJRDUVjuirLzysrZsqkpg4tVGmAAAAgHmybVv9/f2zzjuNjY2ppaXFKU9bt25Vdnb2gmd5raNPT//Rt5Xt8ykv953vFwqFdPzEcVVXVc/p9SLRqELTYT3727frlquNJKdNbZQpAAAA4BLF43Ht379/1qY9n8/nbNkzTVPV1dVzWlO+EH7V2qv/50v/U7F4XMvyZ2/4O3XqlGKxmNauXXvR15mcmlEsHtenHrpVt137roWMnJIoUwAAAMBFhMNhtbe3O5On1tZWrVmzxilOhmGovLx8wdeSX4qDx07oxa/8jfqOnlB+bo6zNOJI/xGVrChRUVHReb82Go1pcnpGa0qW67OPfUj+2g2LFTulUKYAAACAtxkfH1cwGHTKU09Pj6qrq53y1NLSopKSErdjXlQ0FtN3f/ZL/cXf/oumZ8KSRzrWf0Q1NZuV9bbLeWOxuKbDEcXicWX5fPrwje/RrjuuUX7euc9egTIFAAAA6OTJk87jeoFAQMeOHVNjY6MMw5Bpmmpqakr6sojFFI5E9YtAj/7yx/+oV19rV1n5Onm9Xp0dpMXitrwej7ZUlemD17xL11/eoIL8hbuvKl1QpgAAAJBRbNvWkSNHZm3am5iYcIqTaZqqr69XVlb6rf/+3ve+p3379mn3409q8OQZRaJRZfl8Wl1SpIrSEtfOeKWq9PsVAgAAALxFLBZTT0+PM3myLEs5OTlOcXrwwQdVVVWVEUUiGAzq8ssv19pVxVq7qtjtOCmPyRQAAADSyszMjNrb252pU1tbm8rKymZNnsrKytyO6YrbbrtNX/nKV7RhAwslkoHJFAAAAFLa2NiYM3GyLEs9PT2qqamRaZq666679Ad/8AcqLmYKMzQ0pJmZGVVWVrodJW1QpgAAAJBSTpw4Meu80+DgoJqammQYhvbs2aOmpibl5+e7HXPJCQaDMgxjSa1vT3WUKQAAACxZtm3r0KFDs847hUIh55G922+/XbW1tWm5LCLZzpYpJA+/6gAAALBkxGIxdXd3O1Mny7JUUFDglKddu3apqqqK6UoCLMvSLbfc4naMtMICCgAAALhmenpabW1tTnFqb29XeXm5syjCNE2Vlpa6HTPlTUxM6Oabb9bLL7+s7Oxst+OkDSZTAAAAWDSjo6MKBoPO5Km3t1e1tbUyDEP33HOPWlpatHz5crdjpp3W1lY1NDRQpJKMMgUAAIAFMzQ0NOuRvaGhITU3N8s0TX3iE59QY2Oj8vLy3I6Z9izL4rzUAqBMAQAAICni8bgOHToky7KcAjUzMyPTNGUYhnbu3Kna2lr5fD63o2acYDCohx56yO0YaYcyBQAAgIREo1F1dXU5xSkYDGrZsmUyDEOXXXaZHnnkEW3YsIFlES6LRCLat2+fmpub3Y6SdihTAAAAmJNQKKT29nbnkb2Ojg5VVFTINE194AMf0LPPPsuyiCWou7tblZWVWrZsmdtR0g5lCgAAAOd05syZWY/sHTx4ULW1tTJNU/fdd59aWlpUVFTkdkxcBOelFg5lCgAAALJtW4ODg7PK08mTJ+X3+2Wapp588kk1NjYqNzfX7ai4RMFgUDfeeKPbMdIS90wBAABkoHg8roMHDzrFKRAIKBqNzrrfacuWLfJ6vW5HxTzYtq33v//9+ta3vqW1a9e6HSftMJkCAADIAJFIRJ2dnc7kKRgMavny5TJNU1dccYUeffRRVVZWsiwizfT39ys3N5citUAoUwAAAGkoFAqptbXVKU+dnZ1av369TNPUrbfequeff16rV692OyYWGOelFhZlCgAAIA0MDw/LsiynPB06dEj19fUyDEMPPvig/H4/29wyEGVqYVGmAAAAUoxt2xoYGJh13ml4eNhZFvHUU0+poaFBOTk5bkeFyyzL0r333ut2jLRFmQIAAFji4vG4Dhw44NzvFAgEZNu2syjirrvuUk1NDcsiMMvw8LBGRka0adMmt6OkLcoUAADAEhMOh9XZ2elMnVpbW1VSUiLDMHTllVfqscceU0VFBcsicEHBYFDNzc2U7AVEmQIAAHDZ5OSkgsGgc+aps7NTGzdulGma2rFjh1544QWtWrXK7ZhIMZyXWniUKQAAgEU2PDw867zTkSNHtHXrVhmGoV27dsnv96uwsNDtmEhxwWBQn/zkJ92Okda4tBcAAGAB2bato0ePOuedLMvSyMiIWlpanDNP9fX1LItAUk1PT2v79u166aWXlJub63actMVkCgAAIIni8bh6e3udqZNlWfJ6vTIMQ4Zh6J577tHmzZs5x4IF1dHRoS1btlCkFhhlCgAAYB7C4bA6Ojqc4tTa2qpVq1bJNE29733v0yc/+UmVl5ezLAKLivNSi4MyBQAAcAkmJibU2trqTJ66u7tVXV0twzC0c+dOfe5zn9PKlSvdjokMZ1mW7rzzTrdjpD3OTAEAAFzAqVOnZj2y19/fr4aGBue8U3NzswoKCtyOCTji8biuv/56/fCHP1RJSYnbcdIakykAAIA32bat/v5+52LcQCCgsbEx57zTc889p/r6emVnZ7sdFTiv3t5erVq1iiK1CChTAAAgY8XjcfX09MyaPGVlZTlTp/vvv1/V1dUsi0BK4bzU4qFMAQCAjDEzM/OOZRGlpaUyTVPXXnut9u7dq/LycrdjAvNiWZauvPJKt2NkBM5MAQCAtDU+Pq5gMOhMnvbv369NmzbJNE3n0b0VK1a4HRNIqltuuUV/9md/psrKSrejpD0mUwAAIG2cOHFi1nmngYEBNTY2yjRN7d69W83NzcrPz3c7JrBghoaGFI1GtX79erejZATKFAAASEm2bevIkSNOcQoEApqcnJRhGDJNU7fddpvq6uqUlcWPO8gcZ89Lca/Z4uB3FwAAkBJisZi6u7udyZNlWcrNzXWWRTz44IOqqqpiWQQyGssnFhdlCgAALEnT09Nqb293ylNbW5vKyspkmqZuuOEGPfXUUyorK3M7JrCkWJal22+/3e0YGYMFFAAAYEkYGxuTZVlOedq/f7+2bNniPLbX0tKi4uJit2MCS9b4+LhuvfVWvfzyyzzeukj4uwwAAFxx/PjxWeedhoaG1NTUJNM09bu/+7tqampSXl6e2zGBlNHW1qaGhgaK1CLi7zQAAFhwtm3r0KFDzlmnQCCgqakp57zTjh07VFdXJ5/P53ZUIGVxXmrxUaYAAEDSRaNRdXd3O1Mny7JUWFjoPLL38MMPa+PGjWwcA5LIsiw9/PDDbsfIKJyZAgAA8zY1NaW2tjbnzFN7e7vWrVvnTJ4Mw1BpaanbMYG0FYlEdP311+tnP/uZCgsL3Y6TMZhMAQCASzY6OjrrctwDBw6otrZWhmHonnvuUUtLi5YvX+52TCBjdHV1qbKykiK1yChTAADgogYHB53H9SzL0vHjx9Xc3CzTNPXEE0+osbFRubm5bscEMhbnpdxBmQIAALPE43FnWcTZPyKRiAzDkGEY2rlzp2pra1kWASwhlmXppptucjtGxuHMFAAAGS4Siairq2vW5KmoqMg572SapiorK1kWASxRtm3rxhtv1Le//W3OJi4yJlMAAGSYUCik9vZ2Z+q0b98+rV+/XoZh6Oabb9Zzzz2nNWvWuB0TwBwdOXJEBQUFFCkXUKYAAEhzIyMjzrIIy7LU19enuro6maapj33sY/L7/SoqKnI7JoAEWZallpYWt2NkJMoUAABpxLbtWcsiAoGATp48qZaWFhmGob1796qxsVE5OTluRwWQJCyfcA9lCgCAFBaPx9XX1zfrctxoNOqcdbrzzju1ZcsWeb1et6MCWCCWZen+++93O0ZGokwBAJBCIpGIOjs7neIUDAZVXFws0zS1bds27dmzR+vXr2dZBJAhhoeHNTo6qurqarejZCTKFAAAS1goFFJra6szeers7NSGDRtkmqZuvfVWPf/881q9erXbMQG4xLIsNTc3M312CWUKAIAlZHh42DnrFAgEdPjwYdXX18s0TT300EPy+/1atmyZ2zEBLBGcl3IXZQoAAJfYtq2BgYFZl+OOjIzI7/fLNE09/fTTamhoYFkEgPOyLEt79+51O0bGokwBALBI4vG4ent7Z60pt23bWRZx1113qaamhsd1AMzJ1NSU+vr61NjY6HaUjEWZAgBggYTDYe3bt88pT62trSopKZFpmrrqqqv0+OOPa926dSyLAJCQjo4O1dbWMr12EWUKAIAkmZycVDAYdMpTV1eXqqqqZJqmduzYoc9+9rNauXKl2zEBpAnOS7mPMgUAQIJOnz4967xTf3+/tm7dKtM09du//dvy+/0qKChwOyaANGVZlj760Y+6HSOjeWzbtt0OAQDAUmfbto4ePeqcdQoEAhodHVVLS4sMw5Bpmqqvr+dxGwCLIhaL6frrr9ePf/xjrVixwu04GYvJFAAA5xCPx7V//35n6mRZlnw+n1Oc7r33Xm3atIllEQBc0dvbq9LSUoqUyyhTAADojWUR7e3tsixLlmWptbVVq1evlmmauuaaa/TEE0+ovLycZREAloRgMKiWlha3Y2Q8yhQAICONj4+rtbXVmTz19PSourpahmFo586devHFF1VSUuJ2TAA4J8uy9N73vtftGBmPM1MAgIxw8uRJ53E9y7J09OhRNTQ0OHc8NTU1sSwCQEqwbVu33HKL/vzP/1zr1693O05GYzIFAEg7tm2rv79/1qa98fFxGYYhwzD03HPPqb6+XtnZ2W5HBYBLNjg4qHg8roqKCrejZDzKFAAg5cViMfX09MyaPGVnZztTpwceeEBVVVUsiwCQFoLBoAzD4AznEkCZAgCknJmZGXV0dDhTp7a2Nq1du1aGYei6667T3r17VV5e7nZMAFgQlmWxfGKJoEwBAJa8sbExBYNBZ/LU09OjzZs3yzRNfeQjH9EXvvAF1gMDyBiWZWnHjh1ux4AoUwCAJejEiROzLscdGBhQU1OTDMPQnj171NTUpPz8fLdjAsCiGxsb0+DgoOrq6tyOAlGmAAAus21bhw8fnnU57uTkpHM57m233aa6ujplZfGvLABobW1VU1OTfD6f21EgyhQAYJHFYjF1d3fPWhaRl5cn0zRlGIZ27dqlqqoqDlYDwDlwWe/SQpkCACyo6elptbe3O5On9vZ2lZeXyzRN3XDDDXr66ae1du1at2MCQEqwLEuPPPKI2zHwJsoUACCpxsbGnLNOgUBAvb292rJli0zT1N133y3DMLR8+XK3YwJAygmHw+rq6lJzc7PbUfAmyhQAYF6GhoZmLYsYGhpSc3OzTNPU448/rqamJuXl5bkdEwBSXldXlzZu3KiCggK3o+BNlCkAwJzZtq2DBw/OmjzNzMw4yyJ27Nihuro6DkYDwAKwLEuGYbgdA29BmQIAnFc0GlVXV5dTnILBoAoLC2Wapi677DI98sgj2rBhA8siAGARWJalm2++2e0YeAuPbdu22yEAAEvD1NSU2tranMlTR0eHKioqnE17hmGotLTU7ZgAkHHi8bhuvPFGfec73+H34SWEyRQAZLAzZ87MemSvr69PtbW1Mk1T9913n/x+P8siAGAJOHz4sAoLCylSSwxlCgAyhG3bGhwcdO52CgQCOnHihPx+vwzD0JNPPqnGxkbl5ua6HRUA8Dacl1qaKFMAkKbi8bgOHjzoTJ0CgYCi0aizLOJDH/qQtmzZwrIIAEgBwWCQMrUEUaYAIE1EIhF1dnY6kyfLsrR8+XKZpqkrrrhCjz76qCorK1kWAQApyLIsPfDAA27HwNtQpgAgRYVCIbW1tTlTp87OTq1fv16GYejmm2/Wc889pzVr1rgdEwAwT6dOndLY2JiqqqrcjoK3oUwBQIoYHh6eNXU6ePCg6urqZJqmHnzwQfn9fi1btsztmACAJAsGg2ppaZHX63U7Ct6GMgUAS9DZZRFvPe90+vRpZ1nE3r171djYqJycHLejAgAWGOelli7KFAAsAfF4XAcOHFAgEHA27cXjcZmmKdM09dGPflQ1NTX8V0kAyECWZen3fu/33I6Bc6BMAYALwuGwOjs7nfIUDAZVUlIiwzB05ZVXas+ePVq/fj3LIgAgw4VCIfX19amhocHtKDgHyhQALILJyUm1trY65amzs1MbNmyQaZr64Ac/qM985jNatWqV2zEBAEtMe3u76urqeKx7iaJMAcACGB4ennXe6ciRI9q6dasMw9CuXbvk9/tVWFjodkwAwBLHZb1LG2UKAObJtm0dO3bMKU6WZWlkZER+v1+maeqZZ57R1q1b+a+KAIBLFgwGddddd7kdA+fhsW3bdjsEAKSSeDyu3t7eWeXJ4/HINE0ZhiHTNLV582aWRQAA5iUWi+m6667TT37yExUXF7sdB+fAZAoALiIcDqujo8PZstfa2qpVq1bJNE391m/9lj7xiU9o3bp1LIsAACTV/v37VVZWRpFawihTAPA2ExMTzrKIQCCgrq4uVVdXyzRN3XHHHfrc5z6nlStXuh0TAJDmOC+19FGmAGS8U6dOzXpkr7+/Xw0NDTJNU4888oj8fr8KCgrcjgkAyDDBYFDve9/73I6BC+DMFICMYtu2+vv7nUf2AoGAxsbG1NLS4px52rp1q7Kzs92OCgDIYLZt65ZbbtHXvvY1VVRUuB0H58FkCkBai8fj6unpmTV5ysrKkmmaMk1T999/v6qrq1kWAQBYUgYHB2XbttatW+d2FFwAZQpAWpmZmXnHsojS0lKZpqlrr71WTz75pMrLy1kWAQBY0s6el+LfV0sbZQpAShsfH1cwGHQmTz09Pdq0aZNM09SHP/xhff7zn1dJSYnbMQEAuCQsn0gNlCkAKeXEiRPO1MmyLB07dkyNjY0yDEO7d+9WU1MTyyIAACnPsizt3LnT7Ri4CMoUgCXLtm0dOXLEmToFAgFNTEw4F+Peeuutqq+vV1YWv5UBANLH2NiYhoaGVFtb63YUXAQ/gQBYMmKxmLq7u2dNnnJzc51lEQ8++KCqqqpYFgEASGvBYFBNTU3y+XxuR8FFUKYAuGZ6elodHR3O1KmtrU1lZWUyDEPXX3+9nnrqKZWVlbkdEwCARcV5qdTBPVMAFs3Y2Jgsy3ImT/v371dNTY0zeWppaVFxcbHbMQEAcNUjjzyij3/847r88svdjoKLYDIFYMEcP37ceVwvEAhocHBQTU1NMgxDjz32mJqampSfn+92TAAAloxwOKzu7m41NTW5HQVzQJkCkBS2bevQoUOzytPU1JSzLOL2229XbW0tyyIAALiAzs5OVVVVsZk2RfBTDYCERKNRdXd3O+edLMtSYWGhU5527dqlqqoqLhsEAOAScF4qtVCmAMzJ1NSU2tvbneLU3t6u8vJymaap97///fr0pz+t0tJSt2MCAJDSLMvSbbfd5nYMzBELKACc0+joqPO4XiAQ0IEDB1RbW+tMnlpaWrR8+XK3YwIAkDbi8bi2b9+u7373u1q9erXbcTAHTKYASJIGBwdn3e80NDSk5uZmmaapT37yk2psbFReXp7bMQEASFuHDh1SUVERRSqFUKaADBSPx51lEWf/CIfDMk1ThmFo586dqq2t5bJAAAAWEeelUg9lCsgAkUhEXV1dztTJsiwVFRXJNE29+93v1u/8zu9ow4YNLIsAAMBFlmXJNE23Y+AScGYKSEOhUMhZFhEIBLRv3z5VVFQ4kyfDMFgWAQDAEnP77bfrS1/6kqqrq92OgjliMgWkgZGRkVnnnfr6+lRXVyfTNHX//ferpaVFRUVFbscEAADncfLkSU1MTGjjxo1uR8EloEwBKca27VnLIgKBgE6ePCm/3y/TNPXkk0+qsbFRubm5bkcFAABzFAwGZRiGvF6v21FwCShTwBIXj8fV19c363LcaDQq0zRlmqbuvPNObdmyhd98AQBIYSyfSE2UKWCJiUQi6uzsdIpTMBhUcXGxTNPUtm3btHv3blVWVrIsAgCANGJZlp555hm3Y+ASsYACcFkoFFJra6szeers7FRlZaUzeTIMg/smAABIY6FQSDfddJNeeukl5eTkuB0Hl4DJFLDIhoeHZ513OnTokLZu3SrDMPTQQw/J7/dr2bJlbscEAACLpL29XXV1dRSpFESZAhaQbdsaGBiYdTnu8PCwsyzi6aefVkNDA795AgCQwTgvlbooU0ASxeNx9fb2zlpTbtu288jeXXfdpZqaGpZFAAAAh2VZuvfee92OgQRwZgqYh3A4rH379jnlqbW1VSUlJTIMwylQFRUVLIsAAADnFIvFdN111+mnP/2pli9f7nYcXCImU8AlmJycVDAYdMpTV1eXNm7cKNM0tWPHDn32s5/VypUr3Y4JAABSRE9Pj8rKyihSKYoyBVzA6dOnZ5136u/vd5ZFPPzww/L7/SosLHQ7JgAASFGcl0ptlCngTbZtq7+/X5ZlOZOnM2fOqKWlRaZp6tlnn1V9fT3LIgAAQNJYlqVrrrnG7RhIEGemkLHi8bj279/vTJ0sy5LP53POOxmGoc2bN7MsAgAALAjbtvWBD3xA3/jGN7Ru3Tq34yABTKaQMcLhsNrb253JU2trq1atWiXTNHX11VfriSeeUHl5OcsiAADAohgYGJDX61V5ebnbUZAgyhTS1vj4uFpbW53JU3d3tzZt2iTDMLRz5069+OKLKikpcTsmAADIUGfPS/EfclMXZQpp4+TJk87jepZlqb+/Xw0NDTJNUx//+MfV3NysgoICt2MCAABIYvlEOqBMISWdXRbx1k174+PjMgxDhmHoueeeU319vbKzs92OCgAAcE6WZenOO+90OwbmgTKFlBCLxdTT0zNr8pSdne0sivjYxz6m6upqlkUAAICUMDo6qhMnTqimpsbtKJgHyhSWpJmZGWdZRCAQUFtbm0pLS2Wapq699lrt3buXw5oAACBlBYNBNTU1yefzuR0F80CZwpIwNjamYDDoTJ56enq0efNmmaapj3zkI/rCF76gFStWuB0TAAAgKTgvlR4oU3DFiRMnnOIUCAQ0MDCgxsZGmaap3bt3q7m5Wfn5+W7HBAAAWBCWZWnPnj1ux8A8Uaaw4Gzb1uHDh2ddjjs5Oelcjnvbbbeprq5OWVn8cgQAAOkvHA6rp6dHTU1NbkfBPPHTK5IuFoupu7t7VnnKz893lkU89NBDqqqqYlkEAADISPv27dOmTZt4CicNUKYwb9PT02pvb3fKU3t7u8rKymSaprZv365PfepTWrt2rdsxAQAAlgTOS6UPyhQu2djYmHPWKRAIqLe3V1u2bJFhGLr77rvV0tKi4uJit2MCAAAsSZZl6fbbb3c7BpLAY9u27XaIVGHbtk6fmVBv/3EdGTilqZmwfF6vipbla/P6Um2uXKv8vBy3Yybd0NDQrGURQ0NDampqkmmaMk1TTU1NysvLczsmAADAkhePx7V9+3Z973vf06pVq9yOg3liMjUHZ8Ym9Xf/EtR3f/ZLjYxNyuf1KhyJ6mwP9fm8yvL5FI3F1FizXvfc/F5daWxRdlbq3Rtg27YOHjw4a/I0PT3tFKcdO3aorq6OOxEAAAAScPDgQRUXF1Ok0gSTqQuYCUf1zR+9or/++39VPG4rJydLudlZ8ng85/x827Y1OTUjSSoqzNOzj+zQe40tixn5kkWjUXV1dTnFKRgMqrCw0Nm0Z5qmNm7ceN7vGQAAAHP3N3/zN2ptbdXnPvc5t6MgCShT59FzaFAv/H/f18DJES3Lz5PPd2mb56ZmwgpHorrhikb93oO3qqjw4o/BHTt2TOvWrVvQ4jI1NaW2tjZn8tTR0aF169Y5xckwDJWWli7Y+wMAAGSyF154Qe9617t0xx13uB0FSUCZOoffdBzUp/7rt2XHbS0rSPwsUDxua2xyStXr1+jLzz6gFcsLz/l5tm3ra1/7mr74xS/qb//2b9XQ0JDwe77dmTNnZj2y19fXp9raWmfy1NLSouXLlyft/QAAAHB+t99+u7785S+rqqrK7ShIAs5MvU37/n596v/9K3m9XuUX5M7rtbxej4qX5evQsZN64j//hf7bf9j1jnIWDof1/PPP6wc/+IGi0ahee+21hMuUbdsaHByUZVlOgTp+/Lj8fr9M09STTz6pxsZG5ebO7/sCAADApTtx4oQmJye1ceNGt6MgSShTbzE+Oa3n/uR/yuPxKD83OVv5PB6Plhfmq+/oSX3pL3+m5x+9w/lrw8PD+vjHP67XX39dK1as0MTEhF555RU98MADc3rteDyugwcPOlOnQCCgaDTqTJ127typ2tpalkUAAAAsAcFgUC0tLZxFTyOUqbf40rd+ptGJkIqXFST1dd8oVHn6+aut2n5lk67w1+jAgQN64IEHNDQ0pJKSkjcKXH6+fv3rX8u27XP+QxaJRNTZ2elMnizLUlFRkUzT1OWXX65HH31UlZWV/AMKAACwBHFZb/rhzNSbOvuOaffnv65l+bnyei9t2cRcTU2HVbQsX09++DL97mOPKRwOv+O80tjYmH7+85+rurpaoVBIbW1tztRp3759qqysdCZPhmFozZo1C5IVAAAAyXXffffp05/+tPx+v9tRkCRMpt70vZ//SrZtzypSP/3Kpy74NevrLlPte96v/a/9H50eOKDpyTHlFS5X+eYW1b7nRvmysmd9fn5ejg4fHdTDjz2t3Hj8HUUqGo0qFArpxRdfVCwW08GDB1VXVyfTNPXAAw/I7/erqKgoed80AAAAFsXk5KSOHDmi+vp6t6MgiShTeuNS3n/89T4ty5+9HGL7Q59x/veJQ51q/b/fn/Uxny9bI8cPy7bjarr6QyosXq2JkRNqfeX7isyE5L/2TudzbdvWwMCAhs+MKbukWgVTBxSJRBQKhZw/otGofD6fBgYG9MUvflGNjY3KyUnO2S0AAAC4p62tTVu3buVnuzRDmZLUtr9fkt5xl1Rewb9PjrJy89/xMUkq3VCv0g3//l8YCotXacvEDer+1c+dMhWLxXT48GFNTEzI4/Foxluo/b29km2roKBA+fn5KikpUV5enmZmZjQ6OirTNBfkewUAAMDiO7t8AumFMiWp6+CgItFY0l4vGp5W9pvlKxaLqaOjQ5FoRF6PV16fV15vltZX1yvfF33H1+bk5Oj06dM6fvy41q5dm7RMAAAAcI9lWbr//vvdjoEkW5hNCymmreeIcrKT0ytD4yM6YL2iqqYrJUler1elpaVavXq18vLz3jiT5fFoMuLR1NSU3r7/w+PxyOv1yrKspOQBAACAu6LRqDo6OtTc3Ox2FCQZkylJZyZC8iVhg99MaFy//snXtGZ9rapbrpb0RjkqKytzPicSjejU8JgaGq/XUM+/6cCBA8rOztbMzIxyc3OVn5+vSCSiV199VTfddNO8MwEAAMBdPT09Wrdu3TuWjyH1Uaakd0yHEjEdGtMvf/xnKlpVJmP73ee96yk7K1v5+XnaseNmffjG/6jJyUm1tbXptdde0yuvvKLW1lbZtq1f/vKX884EAAAA91mWxXmpNEWZklSQl6u4HU/466cnx/SvP/6qilaulXnjvfJ6fRf8fK/Xq7zcN9amFxYWatu2bdq2bZsef/xxRaNRdXd3a3x8POE8AAAAWDosy9J1113ndgwsAM5MSdq6qULh8DuXQczF9OSo/vVHf6rcgiI1XrVD4elJTYfGNB0akx0/d0Hzeb3aWL76nH8tKytLjY2N2rZtW0J5AAAAsHTYti3LsmQYhttRsACYTElqrKnQj16+8DTpfE7292hy9JQmR0/ppb/8j7P+2vX3P6eC5Stnfcy2bUVjMW2uZFMfAABAujt69KiysrJmnaFH+qBMSWrcvF62bcu27fOedVq32a91j/3hOz5eWf8eVda/Z87vNTUdVnXFGuXncWEbAABAujs7lTrfz5hIbTzmJ2ldaYn8tZWaCE0v+HvF4nHdffOVC/4+AAAAcF8wGOQRvzRGmXrTPTe/V7aSs9nvfCKRqLKzs3TtexoW7D0AAACwdHBeKr1Rpt60raVG9VXlGl+g6ZRt2wrNhPXxO6/nET8AAIAMMDIyopMnT6qmpsbtKFgglKk3eb1e/YfdO+XzehSJJLbZ70LGJqe1dVOFPnzj3M9XAQAAIHW1trbK7/fL6+VH7nTF/7NvUVm2SnsfuEWT02FFY7Gkve7k1LQK83P0wp6d/MMEAACQIYLBIJf1pjl+sn+bD177Lu3+6A2aDM0onIQJ1fjklHKys/Xl5x5URenKi38BAAAA0gLnpdIfZeoc7rvtKn36kQ8qEo1pbHIqoaUUsVhcoxMhrV29Ql994WFt2cjdAgAAAJliZmZGPT09amxsdDsKFhD3TJ3HrVebMuo36gtf/ZH29R2TbGlZQZ683gvfERCJxhSanpHX69V9t16lh+64Rrk5/G0GAADIJB0dHdq8ebPy8/PdjoIF5LEXchd4GojH43q985D++u//Vf/W1iev16NINKbsLJ9TrKKxuOJxW9lZPvl8Xu247jLdccO7tX4tj/UBAABkoq9//esaHR3V3r173Y6CBcTI5CK8Xq/e3bhJ727cpFMj4+o6OKCugwPad+CYQtMz8nm9KikulFG3UTUbylRfXa68XFafAwAAZLJgMKgdO3a4HQMLjMkUAAAAkETxeFw33HCDfvCDH2jlSp5USmcsoAAAAACSqK+vTyUlJRSpDECZAgAAAJKIleiZgzIFAAAAJBGX9WYOyhQAAACQREymMgdlCgAAAEiS48ePa2pqShs2bHA7ChYBZQoAAABIkmAwKMMw5PF43I6CRUCZAgAAAJKER/wyC2UKAAAASBLLslg+kUEoUwAAAEASTExMqL+/X/X19W5HwSKhTAEAAABJ0NbWpq1btyo7O9vtKFgklCkAAAAgCTgvlXkoUwAAAEASnN3kh8xBmQIAAADmKRqNat++fWpubnY7ChYRZQoAAACYp+7ublVUVKioqMjtKFhElCkAAABgnjgvlZkoUwAAAMA8UaYyE2UKAAAAmAfbthUMBrmsNwNRpgAAAIB56O/vV3Z2tsrKytyOgkVGmQIAAADmgUf8MhdlCgAAAJgHylTmokwBAAAA88B5qcxFmQIAAAASNDIyotOnT6umpsbtKHABZQoAAABIUDAYlN/vl9fLj9WZiP/XAQAAgARxXiqzUaYAAACABAWDQcpUBqNMAQAAAAmYnp7W/v371dDQ4HYUuIQyBQAAACRg3759qqmpUV5enttR4BLKFAAAAJAAzkuBMgUAAAAkgDIFyhQAAABwieLxuFpbW+X3+92OAhdRpgAAAIBLdODAAa1cuVIrV650OwpcRJkCAAAALhGP+EGiTAEAAACXjDIFiTIFAAAAXDIu64VEmQIAAAAuydDQkGZmZlRZWel2FLiMMgUAAABcgrNTKY/H43YUuIwyBQAAAFwCzkvhLMoUAAAAcAkoUziLMgUAAADM0cTEhI4dO6a6ujq3o2AJoEwBAAAAc9Ta2qqtW7cqKyvL7ShYAihTAAAAwBzxiB/eijIFAAAAzBFlCm9FmQIAAADmIBKJqKurS36/3+0oWCIoUwAAAMAcdHd3a/369SosLHQ7CpYIyhQAAAAwBzzih7ejTAEAAABzQJnC21GmAAAAgIuwbVvBYJAyhVkoUwAAAMBF9Pf3Kzc3V6WlpW5HwRJCmQIAAAAugkf8cC6UKQAAAOAiKFM4F8oUAAAAcBGUKZwLZQoAAAC4gOHhYY2MjGjTpk1uR8ESQ5kCAAAALiAYDKq5uVleLz86YzZ+RQAAAAAXwCN+OB/KFAAAAHABlCmcD2UKAAAAOI/p6WkdOHBAjY2NbkfBEkSZAgAAAM6jo6NDW7ZsUW5urttRsARRpgAAAIDz4BE/XAhlCgAAADgPyhQuhDIFAAAAnEM8Hldra6v8fr/bUbBEUaYAAACAc+jt7dXq1atVUlLidhQsUZQpAAAA4Bx4xA8XQ5kCAAAAzoEyhYuhTAEAAADnQJnCxVCmAAAAgLcZGhpSNBrV+vXr3Y6CJYwyBQAAALzN2amUx+NxOwqWMMoUAAAA8DY84oe5oEwBAAAAb0OZwlxQpgAAAIC3GB8f18DAgGpra92OgiWOMgUAAAC8RVtbmxoaGpSVleV2FCxxlCkAAADgLXjED3NFmQIAAADegjKFuaJMAQAAAG+KRCLq7OxUc3Oz21GQAihTAAAAwJu6urq0YcMGFRYWuh0FKYAyBQAAALzJsiy1tLS4HQMpgjIFAAAAvInzUrgUlCkAAABAkm3blClcEsoUAAAAIOnw4cMqLCxUaWmp21GQIihTAAAAgKRgMMh5KVwSyhQAAAAgzkvh0lGmAAAAAFGmcOkoUwAAAMh4w8PDGh0dVXV1tdtRkEIoUwAAAMh4lmXJ7/fL6+XHY8wdv1oAAACQ8bisF4mgTAEAACDjcV4KiaBMAQAAIKNNTU2pr69PjY2NbkdBiqFMAQAAIKO1t7ertrZWOTk5bkdBiqFMAQAAIKMFg0Ee8UNCKFMAAADIaJyXQqIoUwAAAMhYsVhMbW1t8vv9bkdBCqJMAQAAIGP19vaqtLRUK1ascDsKUhBlCgAAABmL81KYD8oUAAAAMhaX9WI+KFMAAADISLZtKxAIMJlCwihTAAAAyEiDg4OKx+OqqKhwOwpSFGUKAAAAGenseSmPx+N2FKQoyhQAAAAyEuelMF+UKQAAAGQkLuvFfFGmAAAAkHHGxsY0ODiouro6t6MghVGmAAAAkHFaW1vV1NQkn8/ndhSkMMoUAAAAMg6X9SIZKFMAAADIOCyfQDJQpgAAAJBRwuGwurq61Nzc7HYUpDjKFAAAADJKV1eXNm7cqIKCArejIMVRpgAAAJBRWImOZKFMAQAAIKNwXgrJQpkCAABAxojH4woGg5QpJAVlCgAAABnj8OHDKiwsVGlpqdtRkAYoUwAAAMgYnJdCMlGmAAAAkDG4rBfJRJkCAABAxmAyhWSiTAEAACAjnDp1SmNjY6qqqnI7CtIEZQoAAAAZ4ewWP6+XH4GRHPxKAgAAQEbgvBSSjTIFAACAjMBlvUg2yhQAAADSXigUUl9fnxoaGtyOgjRCmQIAAEDaa29vV11dnXJyctyOgjRCmQIAAEDaYyU6FgJlCgAAAGnv7CY/IJkoUwAAAEhrsVhMbW1tlCkkHWUKAAAAaW3//v0qKytTcXGx21GQZihTAAAASGucl8JCoUwBAAAgrXFZLxYKZQoAAABpy7ZtLuvFgqFMAQAAIG0NDg7Ktm2tW7fO7ShIQ5QpAAAApK2z56U8Ho/bUZCGKFMAAABIWyyfwEKiTAEAACBtcV4KC4kyBQAAgLQ0NjamoaEh1dbWuh0FaYoyBQAAgLQUDAbV1NQkn8/ndhSkKcoUAAAA0hLnpbDQKFMAAABIS1zWi4VGmQIAAEDaCYfD6u7uVlNTk9tRkMYoUwAAAEg7nZ2dqqqqUkFBgdtRkMYoUwAAAEg7nJfCYqBMAQAAIO1wXgqLgTIFAACAtBKPx7msF4uCMgUAAIC0cujQIRUVFWn16tVuR0Gao0wBAAAgrXBeCouFMgUAAIC0QpnCYqFMAQAAIK2wfAKLhTIFAACAtHHy5EmNj49r48aNbkdBBqBMAQAAIG2cnUp5vfyYi4XHrzIAAACkDc5LYTFRpgAAAJA2OC+FxUSZAgAAQFoIhUI6dOiQ6uvr3Y6CDEGZAgAAQFpob29XXV2dcnJy3I6CDEGZAgAAQFrgvBQWG2UKAAAAaYEyhcVGmQIAAEDKi8Viam9vl9/vdzsKMghlCgAAACmvp6dHZWVlWr58udtRkEEoUwAAAEh5POIHN1CmAAAAkPIoU3ADZQoAAAApzbZtLuuFKyhTAAAASGkDAwPyeDwqLy93OwoyDGUKAAAAKe3sI34ej8ftKMgwWW4HAAAAAObimWee0auvvqqrrrpKV111lQzDUGVlJeel4BrKFAAAAFJCRUWFDh48qMHBQf3whz+UJOXn52tsbEy7d+9We3u76urqlJ2d7XJSZAqPbdu22yEAAACAi3n11Ve1a9cuFRUVOR+bnp7WgQMHtG7dOmVlZcm2bfn9fl199dW69tprucQXC4oyBQAAgJQwOTkpv9+v4uJieb1vHP2fmJjQ8PCwNmzYIEmKxWKamprS8PCwbrjhBn3nO99xMzLSHAsoAAAAkBIKCwu1efNmTU9POx8LhUIqKChw/tzn80mS1q9frz/+4z9e9IzILJQpAAAApIxrrrlGU1NTzp+HQiHl5+c7fz4zMyPbtvXNb35TZWVlbkREBqFMAQAAIGVs27ZNOTk5kt64rHdmZsYpU9FoVJOTk/qjP/ojzkphUVCmAAAAkDJM01Q0GpVt25qamlJOTo68Xq9s29bY2Jj27NmjD37wg27HRIZgNToAAACWrEg0psMDJ9V75Lj6h4Y1NR1WYaWh8clRheMxFRQWybZtnTlzRtu3b9dTTz3ldmRkELb5AQAAYMnpPXJcf/MPv9bPftEqSYrH44rG4vJ4PDp9+rRCoUnFYzF5fVnyzIyoYllUf/+Db6moaJnLyZFJKFMAAABYMk6OjOmL//0n+rf2Ptm2rcL8XGW9uaHvrDOjZ3S0/6gikYiysrPl9WVrXUWFylav0At7PiR/7QaX0iPTUKYAAACwJPzDv7bpv3z9p5oJR7S8MF8ej+ecnxcOh9Xd3a1oNKrs7Gxt2rRJBQUFmpyaUSwW187t79Fjd29XTjYnWrCwWEABAAAA133rp/+i3//TH8ojqXhZwXmLlCRnm588b9wndfaeqcL8XBXm5+oH//BrPfNfv6PpmfAiJEcmo0wBAADAVd//37/Sn333ZRUW5Co3J3tOX5Obm6uSFSVasWLFrI/7fF4VL8vXbzoO6vkvf1fRWGwBEgNv4DE/AAAAuGb/4SH9zmf/XHm5OcrO8l38C95ky5ZH559e2batsYkp7b7rBt1761XJiAq8A5MpAAAAuCIciepzX/mBPB5dUpGSdMEiJUkej0cF+bn68+//ow4dOzmfmMB5UaYAAADgir/7J0tHhk6rMD9vQV4/O8unuG3rS9/62YK8PsCKEwAAACy6eDyuv/pfv1BudtasZRM//cqnLvh16+suk3HD3dr/2ks6caRTY6cGFItGdNtjf3jOz19WkKfXOw/p2IlhVZSuTOr3AFCmAAAAsOhae/p1cmRcy/JzZ318+0Ofcf73iUOdav2/35/1MZ/vjQUV8XhUZZuatWrdZvW+/vJ538fr8ci2bf345d/osbtvTPJ3gUzHY34AAABYdK919Ckajb1jBXpewXLnj6zc/Hd8LPvNj9VdfpM2G9do+ZqKi75Xbna2/uX17uR/E8h4lCkAAAAsOqvr8KJdqpuTk6Wjx4c1E44uyvshc1CmAAAAsOj2HxlSbs7ilCmvx6Msn09HBk8tyvshc1CmAAAAsOhCU2H5vIv3o6jHI02Ephft/ZAZKFMAAABYfLa96G8Zd+E9kd4oUwAAAFh0eXk5iscXsdzYUkFezuK9HzICZQoAAACLrrpijWYikYS/fmp8RKOnjmlqbFiSNHrqmEZPHVM0MvOOz7VtW+FoTFXr1iT8fsC5cM8UAAAAFp1Zv1H7DhxTQV7uxT/5HLp//XMd7f6N8+f//N0/kSRt27Fbqys2z/rccCSqslXFymcyhSTz2DYPjwIAAGBx/bK1V8/+8V+rqCBvwd9rdCKk26+9TE/vunXB3wuZhcf8AAAAsOje07hJRYV5mgkn/qjfXNi2LY/HoztuePeCvg8yE2UKAAAAi87n8+quD2zT9MzClqmJ0IxqN5apZsPaBX0fZCbKFAAAAFyx84b3qKS4UKHpdy6NSIZYLC7btvXEx25ekNcHKFMAAABwRWF+rl7Y8yFFojHF4vGkvrZt25oITesjN12hppr1SX1t4CzKFAAAAFxzWUO17vrANk2EphVPUqGybVvjk9Oq2bBWj3z42qS8JnAulCkAAAC4as9d2/XBa9+lsclpRWOxeb2Wbdsam5zSxnWr9cef/pjyclmHjoXDanQAAAC4zrZt/cXf/rO+8aN/kmxbywry5PF4Luk1ZsIRTYcj2uav0Wd2f0hFhQu/dh2ZjTIFAACAJaP3yHG9+Kc/UP/QsGTbKizIk/cCpcq2bU2HI4pEY8rLzdYzu27T9Vc0XnIRAxJBmQIAAMCSEovF9au2A/rrv39Vwe4j8nm9CkeiyvL5dLYjRaIxZWX5JNvWmpXLde+t79X2bc1Mo7CoKFMAAABYsk6NjKvn0KC6Dw2q98iQpsNRZWf5VL5mhRo2V6hmw1ptLF8tr5dVAFh8lCkAAAAASAAVHgAAAAASQJkCAAAAgARQpgAAAAAgAZQpAAAAAEgAZQoAAAAAEkCZAgAAAIAEUKYAAAAAIAGUKQAAAABIAGUKAAAAABJAmQIAAACABFCmAAAAACABlCkAAAAASABlCgAAAAASQJkCAAAAgARQpgAAAAAgAZQpAAAAAEgAZQoAAAAAEkCZAgAAAIAEUKYAAAAAIAGUKQAAAABIAGUKAAAAABJAmQIAAACABFCmAAAAACABlCkAAAAASABlCgAAAAASQJkCAAAAgARQpgAAAAAgAZQpAAAAAEgAZQoAAAAAEvD/AxsT6ob6CwYkAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Set interactive=False for web publishing\n", "# If you have `pyviz` installed you can get an interactive graph!\n", "bn.plot(model, title='Disease Test', interactive=False);" ] }, { "cell_type": "markdown", "id": "b3caf7b0", "metadata": {}, "source": [ "There's also a handy `print_CPD` method that lets you pretty print all the CPD's." ] }, { "cell_type": "code", "execution_count": 28, "id": "fcdc123d", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPD of D:\n", "+--------+------+\n", "| D(yes) | 0.02 |\n", "+--------+------+\n", "| D(no) | 0.98 |\n", "+--------+------+\n", "CPD of T1:\n", "+--------------+--------+-------+\n", "| D | D(yes) | D(no) |\n", "+--------------+--------+-------+\n", "| T1(positive) | 0.95 | 0.1 |\n", "+--------------+--------+-------+\n", "| T1(negative) | 0.05 | 0.9 |\n", "+--------------+--------+-------+\n", "CPD of T2:\n", "+--------------+--------+-------+\n", "| D | D(yes) | D(no) |\n", "+--------------+--------+-------+\n", "| T2(positive) | 0.95 | 0.1 |\n", "+--------------+--------+-------+\n", "| T2(negative) | 0.05 | 0.9 |\n", "+--------------+--------+-------+\n", "[bnlearn] >Independencies:\n", "(T2 ⟂ T1 | D)\n", "(T1 ⟂ T2 | D)\n", "[bnlearn] >Nodes: ['D', 'T1', 'T2']\n", "[bnlearn] >Edges: [('D', 'T1'), ('D', 'T2')]\n" ] } ], "source": [ "bn.print_CPD(model)" ] }, { "cell_type": "markdown", "id": "7604b7cc", "metadata": {}, "source": [ "### Inference\n", "\n", "You can also use `bnlearn.inference` to perform inference using Variable Elimination, similar to what we saw in `pgmpy`:" ] }, { "cell_type": "code", "execution_count": 29, "id": "2ee90ac5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >Variable Elimination..\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f71bd228c62148268e82b961eb89c365", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "12c596bcd25c49da86da1815db84d766", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "+----+-----+----------+\n", "| | D | p |\n", "+====+=====+==========+\n", "| 0 | 0 | 0.162393 |\n", "+----+-----+----------+\n", "| 1 | 1 | 0.837607 |\n", "+----+-----+----------+\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bn.inference.fit(model, variables=['D'], evidence={'T1':'positive'})" ] }, { "cell_type": "markdown", "id": "aa7ba68c", "metadata": {}, "source": [ "*Note: Interesting that it doesn't show the 'state' names here, so the above result is 'yes' for 0 and 'no' for 1 because that's the order in which we defined the states.*\n", "\n", "Surprisingly, with the given assumptions, with a single positive test there would only be a 16% probability that you'd have the disease! What's it look like if you had two positive tests?" ] }, { "cell_type": "code", "execution_count": 30, "id": "65c03a34", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >Variable Elimination..\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c2b1a32bf04547119ef10418ea6dcb8d", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "67023c90e90c4d3381cc2d1003fe9064", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "+----+-----+----------+\n", "| | D | p |\n", "+====+=====+==========+\n", "| 0 | 0 | 0.648115 |\n", "+----+-----+----------+\n", "| 1 | 1 | 0.351885 |\n", "+----+-----+----------+\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bn.inference.fit(model, variables=['D'], evidence={'T1':'positive', 'T2':'positive'})" ] }, { "cell_type": "markdown", "id": "cac16449", "metadata": {}, "source": [ "So we see here that there's a 65% probability of having this disease after testing positive *twice*. A pretty dramatic improvement from the 16% earlier. But still, this is really low!" ] }, { "cell_type": "markdown", "id": "554877c9", "metadata": {}, "source": [ "## What about a real-world example with Covid testing?\n", "\n", "The above example uses numbers to help surprise the learner about how counterintuitive statistics can be. But should we be skeptical of *any* test, especially in 2022?\n", "\n", "If we change the prevalence and sensitivity/specificity to, say, real-world numbers for Covid-19, the numbers become more reassuring.\n", "\n", "We can check real world numbers thanks to resources like:\n", "* [FDA Emergency Use Authorization statistics for Covid Tests](https://www.fda.gov/medical-devices/coronavirus-disease-2019-covid-19-emergency-use-authorizations-medical-devices/eua-authorized-serology-test-performance)\n", "* [CDC Covid Data Tracker for disease prevalence](https://covid.cdc.gov/covid-data-tracker/#datatracker-home)" ] }, { "cell_type": "markdown", "id": "7221e59c", "metadata": {}, "source": [ "### Estimating disease prevalence\n", "\n", "You can use the [CDC Covid Tracker](https://covid.cdc.gov/covid-data-tracker/#datatracker-home) to get a sense of positive cases per 100K people. Let's pull the data for Multnomah County, host of the fine city of Portland, Oregon:\n", "\n", "|![covid_prevalence](covid_prevalence.png)|\n", "|:---:|\n", "|Source: [CDC Covid Tracker](https://covid.cdc.gov/covid-data-tracker/#datatracker-home)|\n" ] }, { "cell_type": "markdown", "id": "03b32cab", "metadata": {}, "source": [ "Let's say that this published rate only captures 1/4 of actual cases, due to folks not reporting (i.e. people testing at home, people who are asymptomatic and don't get tested, etc), then the true prevalence would be about 2.6%:" ] }, { "cell_type": "code", "execution_count": 31, "id": "c182253d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.025874400000000002" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(646.86 * 4)/100000" ] }, { "cell_type": "markdown", "id": "3ffa4c08", "metadata": {}, "source": [ "### FDA data on covid test performance\n", "\n", "You can also get sensitivity and specificity data on the covid tests that the FDA approved for emergency use in the early days of the pandemic: \n", "\n", "Here we'll use data on the Abbott AdviseDx SARS-CoV-2 IgG II (Alinity), which has a sensitivity of 98.1% (51/52), specificity of 99.6% (2000/2008). Because of the way that we've formulated this problem, we can substitute the following:\n", "\n", "$$\n", "\\begin{aligned}\n", "\\text{Sensitivity} &= P(D=\\text{yes}|T=\\text{'positive'}) = 51/52 = 98.0769\\%\\\\\n", "\\text{Specificity} &= P(D=\\text{no}|T=\\text{'negative'}) = 2000/2008 = 99.6016\\%\\\\\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "id": "e7f4fafd", "metadata": {}, "source": [ "Or as a conditional probability table:\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", "\t \n", "\t \n", " \n", " \n", " \n", " \n", "\t \n", " \n", " \n", " \n", "\n", "
Test Result
PositiveNegative
DiseaseYes0.9810.019
No0.0040.996
" ] }, { "cell_type": "markdown", "id": "81efc987", "metadata": {}, "source": [ "For convenience, from here on out we'll redefine this as a function:" ] }, { "cell_type": "code", "execution_count": 32, "id": "d18b7c0e", "metadata": {}, "outputs": [], "source": [ "def disease_test(prevalence, sens, spec):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " prevalence : float\n", " Estimated percent of population that has the disease\n", " sens : float\n", " Sensitivity, or true positive rate\n", " spec : float\n", " Specificity, or true negative rate\n", " \"\"\"\n", " edges = [('D', 'T1'), ('D', 'T2')]\n", "\n", " d = TabularCPD(variable='D', variable_card=2, values=[[prevalence], [1-prevalence]],\n", " state_names={'D':['yes', 'no']})\n", " t1 = TabularCPD(variable='T1', variable_card=2, \n", " values=[[sens, 1-spec],\n", " [1-sens, spec]],\n", " evidence=['D'],\n", " evidence_card=[2],\n", " state_names={'D':['yes', 'no'],\n", " 'T1':['positive', 'negative']})\n", " t2 = TabularCPD(variable='T2', variable_card=2, \n", " values=[[sens, 1-spec],\n", " [1-sens, spec]],\n", " evidence=['D'],\n", " evidence_card=[2],\n", " state_names={'D':['yes', 'no'],\n", " 'T2':['positive', 'negative']})\n", " # Make the actual Bayesian DAG with the previously defined CPD's \n", " model = bn.make_DAG(DAG=edges, CPD=[d, t1, t2], verbose=1) # quiet those messages!\n", " return model" ] }, { "cell_type": "code", "execution_count": 33, "id": "ee07c5b7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >Variable Elimination..\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "fdf72908948d4f5b993b90bde28f9ce8", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a759925e004e49a6874ebb8d9f733e7c", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "+----+-----+----------+\n", "| | D | p |\n", "+====+=====+==========+\n", "| 0 | 0 | 0.867923 |\n", "+----+-----+----------+\n", "| 1 | 1 | 0.132077 |\n", "+----+-----+----------+\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = disease_test(0.026, sens=51/52, spec=2000/2008)\n", "bn.inference.fit(model, variables=['D'], evidence={'T1':'positive'})" ] }, { "cell_type": "markdown", "id": "7225da2b", "metadata": {}, "source": [ "Reassuringly, the probability of having the disease is about 87%, after having tested positive with this test (and assuming prevalence of 2.6%)." ] }, { "cell_type": "code", "execution_count": 34, "id": "94ea9ce6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >Variable Elimination..\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "3690e9e4f79d4f8f92b703fde7c5f139", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "85873fbc1de74c95a48a19c328c86c87", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "+----+-----+-------------+\n", "| | D | p |\n", "+====+=====+=============+\n", "| 0 | 0 | 0.999382 |\n", "+----+-----+-------------+\n", "| 1 | 1 | 0.000617783 |\n", "+----+-----+-------------+\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bn.inference.fit(model, variables=['D'], evidence={'T1':'positive', 'T2':'positive'})" ] }, { "cell_type": "markdown", "id": "2bfe27a8", "metadata": {}, "source": [ "And after 2 tests, the probability rises to 99.9%. Not bad, right?\n", "\n", "Now, let's take a quick look with a rapid at-home test like the Abbot BinaxNOW test. There seems to be some conflicting data on this (i.e. the results are better if you're symptomatic), but we can try to use some available overall numbers. Skimming the results of Table 2 from [this Nov 2020 study](https://www.cdc.gov/mmwr/volumes/70/wr/mm7003e3.htm) gives us a Sensitivity of 52.5% and Specificity of 99.9%:\n", "\n", "|
|\n", "|:---:|\n", "|![binax_sens_spec](binax_sens_spec.png)|\n", "| From: [Prince-Guerra JL, Almendares O, Nolen LD, et al. Evaluation of Abbott BinaxNOW Rapid Antigen Test for SARS-CoV-2 Infection at Two Community-Based Testing Sites — Pima County, Arizona, November 3–17, 2020. MMWR Morb Mortal Wkly Rep 2021](https://www.cdc.gov/mmwr/volumes/70/wr/mm7003e3.htm)|" ] }, { "cell_type": "code", "execution_count": 35, "id": "43b56ba8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >Variable Elimination..\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c667a113980c4dffb478d6a9e27a7771", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ec3023826c464b6cbfe49f3919427c45", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "+----+-----+-----------+\n", "| | D | p |\n", "+====+=====+===========+\n", "| 0 | 0 | 0.933397 |\n", "+----+-----+-----------+\n", "| 1 | 1 | 0.0666028 |\n", "+----+-----+-----------+\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = disease_test(0.026, sens=0.525, spec=0.999)\n", "bn.inference.fit(model, variables=['D'], evidence={'T1':'positive'})" ] }, { "cell_type": "markdown", "id": "89c118ce", "metadata": {}, "source": [ "So that's *quite* good, considering the study found that sensitivity was much higher among symptomatic patients. We can also take a look if you take two tests:" ] }, { "cell_type": "code", "execution_count": 36, "id": "c61c06ae", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >Variable Elimination..\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "3ee1cb958bd9490fa0ae9b2674c6f063", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9467ad9668bf408f89698fd0592949b4", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "+----+-----+-------------+\n", "| | D | p |\n", "+====+=====+=============+\n", "| 0 | 0 | 0.999864 |\n", "+----+-----+-------------+\n", "| 1 | 1 | 0.000135896 |\n", "+----+-----+-------------+\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bn.inference.fit(model, variables=['D'], evidence={'T1':'positive', 'T2':'positive'})" ] }, { "cell_type": "markdown", "id": "8ff65fbd", "metadata": {}, "source": [ "And now there's near certainty with two tests. What if you had one test positive and another negative, out of curiosity?" ] }, { "cell_type": "code", "execution_count": 37, "id": "d21449b2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bnlearn] >Variable Elimination..\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "3acc19dbbb67438585bef61421861e53", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "0eba0daac30f4fef9cd28aba4105733b", "version_major": 2, "version_minor": 0 }, "text/plain": [ "0it [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "+----+-----+----------+\n", "| | D | p |\n", "+====+=====+==========+\n", "| 0 | 0 | 0.869511 |\n", "+----+-----+----------+\n", "| 1 | 1 | 0.130489 |\n", "+----+-----+----------+\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bn.inference.fit(model, variables=['D'], evidence={'T1':'positive', 'T2':'negative'})" ] }, { "cell_type": "markdown", "id": "64667c1d", "metadata": {}, "source": [ "## References\n", "\n", "* [Phillip Loick on Bayesian Statistics (Udemy)](https://www.udemy.com/course/bayesian-statistics/)\n", "* `bnlearn`: https://github.com/erdogant/bnlearn, much newer and very promising \n", "* `pymc` thread on bayes nets: https://discourse.pymc.io/t/bayes-nets-belief-networks-and-pymc/5150/8\n", "* [FDA Emergency Use Authorization statistics for Covid Tests](https://www.fda.gov/medical-devices/coronavirus-disease-2019-covid-19-emergency-use-authorizations-medical-devices/eua-authorized-serology-test-performance)\n", "* [CDC Covid Data Tracker for disease prevalence](https://covid.cdc.gov/covid-data-tracker/#datatracker-home)\n", "* [Prince-Guerra JL, Almendares O, Nolen LD, et al. Evaluation of Abbott BinaxNOW Rapid Antigen Test for SARS-CoV-2 Infection at Two Community-Based Testing Sites — Pima County, Arizona, November 3–17, 2020. MMWR Morb Mortal Wkly Rep 2021](https://www.cdc.gov/mmwr/volumes/70/wr/mm7003e3.htm)" ] }, { "cell_type": "code", "execution_count": null, "id": "5c78553a", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }