examples/Natality_Experiment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Natality Experiment on the DAME-FLAME Python Package\n",
    "\n",
    "I use the DAME-FLAME Python package, which is publicly available for download [here](https://pypi.org/project/dame-flame/), on a real dataset, aiming to examine the effect of \"extreme smoking\" during pregancy on a newborn's birthweight. Overall, the results of the DAME-FLAME package lead us to believe that smoking during pregnancy causes significantly lower birthweight of infants. It also leads us to believe that there is little to no heterogeneity in *conditional treatment effects* of a mother's smoking on birthweight. In other words, when we examine groupings of characteristics of mothers and infants (including race, gender of infant, etc), then we don't see evidence that there exists some mothers and infants for whom the impact of maternal smoking on infant birthweight is different than the outcomes for other types of mothers and infants. \n",
    "\n",
    "I use the 2010 US Natality dataset from the National Vital Statistics System of the National Center for Health Statistics (2010), available here: http://www2.nber.org/data/vital-statistics-natality-data.html.\n",
    "The accompanying user guide to the dataset is also linked on this URL, as \"documentation\". \n",
    "\n",
    "This code and analysis was written by Neha Gupta, following similar work done by Tianyu Wang, in the paper [FLAME: A Fast Large-scale Almost Matching Exactly Approach to Causal Inference](https://arxiv.org/abs/1707.06315). Notable differences between the two are that the latter uses a slightly modified algorithm (Wang uses FlameDB, I use FlameBit), and this experiment drops all units that have missing data rather than matching on units with missing data while skipping matching on missing covariates (Both options are available in the `dame-flame` package! Use the `missing_data_replace` parameter). I also cut back on the dataset size for processing time. \n",
    "\n",
    "This experiment aims to build upon the analysis on maternal smoking and infant birthweight shown in [Kondracki (2020)](https://link.springer.com/content/pdf/10.1186/s12884-020-02981-1.pdf). \n",
    "\n",
    "For questions on this code, or on the DAME-FLAME Python Package, please contact Neha at neha.r.gupta \"at\" duke \"dot\" edu. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import dame_flame\n",
    "import random\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Cleaning \n",
    "We begin by reading in the data file from the URL and cleaning it. We take a randomly chosen small sample, the variable `p`, for this experiment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/home1/nrg18/.local/lib/python3.6/site-packages/IPython/core/interactiveshell.py:3072: DtypeWarning: Columns (28,55,60,61,62,63,64,65,66,67,68,69,70,71,77,78,79,80,83,84,85,86,87,88,89,90,91,92,93,99,104,105,106,107,108,109,111,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,158,159,160) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n"
     ]
    }
   ],
   "source": [
    "p = 0.1\n",
    "\n",
    "df = pd.read_csv('http://data.nber.org/natality/2010/natl2010.csv.zip',\n",
    "         header=0, \n",
    "         skiprows=lambda i: i>0 and random.random() > p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict to singular births only. This is the DPLURAL variable\n",
    "df = df[df['dplural'] == 1]\n",
    "\n",
    "# need to also restrict to standard term births. This is 37 weeks+, as defined by Kondracki (2020)\n",
    "df = df[df['gestrec3'] == 2]\n",
    "\n",
    "# These are missing from the weight, I presume from the User Guide. \n",
    "df = df[df['dbwt'] != 9999]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, I recode cigarette usage. I follow [Kondracki (2020)](https://link.springer.com/content/pdf/10.1186/s12884-020-02981-1.pdf) and define \"extreme smoking\" as 10+ cigarettes per day for the duration of preganancy. I will use this as the treatment indicator, and the control group will be individuals who smoked 0 cigarettes per day for the duration of the pregnancy. The dataset encodes cigarette usage for each trimester of pregnancy differently, so I gather from each variable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, drop anyone for whom we don't have smoking data:\n",
    "df = df[df['cig_1'] != 99]\n",
    "df = df[df['cig_2'] != 99]\n",
    "df = df[df['cig_3'] != 99]\n",
    "\n",
    "df['smokes_10'] = np.where(((df['cig_1'] >= 10) & (df['cig_2'] >=10) & (df['cig_3'] >= 10)), 1, 0)\n",
    "df['nonsmoker'] = np.where(((df['cig_1'] == 0) & (df['cig_2'] == 0) & (df['cig_3'] == 0)), 1, 0)\n",
    "\n",
    "# if you're not smoking 10 and you're not a nonsmoker then these two values will both be 0. We remove where they're equal. \n",
    "df = (df[df['nonsmoker'] != df['smokes_10']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I remove many observational covariates from the original dataset. The variables names here are the exact variable names in the original dataset, so they can be searched for in the [User Guide](http://data.nber.org/natality/2010/natl2010.pdf). More details on the meaning are available there. Briefly, I summarize the variables below. Note that some of these variables are recoded and could be considered \"binned\" continuous variables. The `DAME-FLAME` package is only for discrete observational datasets, but there are acceptable forms of data binning for some research questions. In this case, we assume binning these variables is acceptable because of prior research that relies on binned data such as [Kondracki (2020)](https://link.springer.com/content/pdf/10.1186/s12884-020-02981-1.pdf). For example, data that is binned include fathers age into fixed width bins in intervals of 5 years, and mother's education into bins of degree obtained rather than years of study. There is substantial research precedence that uses educational degree obtained as a variable on outcome rather than years of study. \n",
    "\n",
    "| Variable Name | Meaning |\n",
    "| ------------- | ------- |\n",
    "| fagerec11 | Fathers Age Recoded |\n",
    "| meduc | Mothers Education |\n",
    "| dob_wk | Weekday |\n",
    "| mager9 | Mother's Age Recode 9 |\n",
    "| fbrace | Father's Race |\n",
    "| mar | Mother's marital status |\n",
    "| rf_diab | Pre-pregnancy diabetes |\n",
    "| rf_cesar | Previous Cesarean Deliveries |\n",
    "| tbo_rec | Total Birth Order |\n",
    "| rf_phyp | Prepregnancy hypertension |\n",
    "| urf_chyper | chronic hypertension |\n",
    "| lbo_rec | live birth order |\n",
    "| rf_ppterm | previous preterm birth |\n",
    "| sex | Sex of Infant |\n",
    "| mbrace | Mother's bridged race |\n",
    "| bfacil3 | Birth Place (home, hospital) |\n",
    "| dbwt | Infant's birthweight |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_trunc = df.loc[:, df.columns.intersection(['fagerec11', 'meduc', 'dob_wk', 'mager9',\n",
    "                                    'fbrace', 'mar','rf_diab', 'rf_cesar', \n",
    "                                  'tbo_rec', 'rf_phyp', 'urf_chyper', 'lbo_rec',\n",
    "                                  'rf_ppterm', 'sex', 'mbrace', 'bfacil3', 'smokes_10', 'dbwt'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Missing values are encoded differently throughout the dataset (see user guide), but I aim to encode them all with the 999 character\n",
    "\n",
    "df_trunc.loc[df_trunc['bfacil3'] == 9, 'bfacil3'] = np.nan\n",
    "df_trunc.loc[df_trunc['mar'] == 9, 'mar'] = np.nan\n",
    "df_trunc.loc[df_trunc['meduc'] == 9, 'meduc'] = np.nan\n",
    "df_trunc.loc[df_trunc['fagerec11'] == 11, 'fagerec11'] = np.nan\n",
    "df_trunc.loc[df_trunc['fbrace'] == 99, 'fbrace'] = np.nan\n",
    "df_trunc.loc[df_trunc['lbo_rec'] == 9, 'lbo_rec'] = np.nan\n",
    "df_trunc.loc[df_trunc['tbo_rec'] == 9, 'tbo_rec'] = np.nan\n",
    "df_trunc.loc[df_trunc['urf_chyper'] == 9, 'urf_chyper'] = np.nan\n",
    "\n",
    "# Need to transform all the \"N\" and \"Y\" to 0/1\n",
    "d = {\"N\": 0, \"Y\": 1, \"U\": np.nan}\n",
    "df_trunc['rf_diab'] = df_trunc['rf_diab'].map(d)\n",
    "df_trunc['rf_phyp'] = df_trunc['rf_phyp'].map(d)\n",
    "df_trunc['rf_ppterm'] = df_trunc['rf_ppterm'].map(d)\n",
    "df_trunc['rf_cesar'] = df_trunc['rf_cesar'].map(d)\n",
    "\n",
    "# Transforming Gender to boolean\n",
    "d = {\"F\":1, \"M\": 0}\n",
    "df_trunc['sex'] = df_trunc['sex'].map(d)\n",
    "\n",
    "# Make sure all clear for input to DAME-FLAME package\n",
    "df_trunc = df_trunc.dropna()\n",
    "df_trunc = df_trunc.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>dob_wk</th>\n",
       "      <th>bfacil3</th>\n",
       "      <th>mager9</th>\n",
       "      <th>mbrace</th>\n",
       "      <th>mar</th>\n",
       "      <th>meduc</th>\n",
       "      <th>fagerec11</th>\n",
       "      <th>fbrace</th>\n",
       "      <th>lbo_rec</th>\n",
       "      <th>tbo_rec</th>\n",
       "      <th>rf_diab</th>\n",
       "      <th>rf_phyp</th>\n",
       "      <th>rf_ppterm</th>\n",
       "      <th>rf_cesar</th>\n",
       "      <th>urf_chyper</th>\n",
       "      <th>sex</th>\n",
       "      <th>dbwt</th>\n",
       "      <th>smokes_10</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>19864</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>3544</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>19865</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>7</td>\n",
       "      <td>21</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>3487</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>19866</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2642</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>19867</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>5</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>3220</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>19869</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>3118</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       dob_wk  bfacil3  mager9  mbrace  mar  meduc  fagerec11  fbrace  \\\n",
       "19864       6        1       3       1    1      3          5       1   \n",
       "19865       1        2       7      21    1      4          8       1   \n",
       "19866       6        1       5       5    1      6          5       5   \n",
       "19867       7        1       5       8    1      7          5       8   \n",
       "19869       1        1       6       1    1      1          7       1   \n",
       "\n",
       "       lbo_rec  tbo_rec  rf_diab  rf_phyp  rf_ppterm  rf_cesar  urf_chyper  \\\n",
       "19864        2        2        0        0          0         0           2   \n",
       "19865        3        3        0        0          0         0           2   \n",
       "19866        1        1        0        0          0         0           2   \n",
       "19867        1        1        0        0          0         0           2   \n",
       "19869        3        3        0        0          0         0           2   \n",
       "\n",
       "       sex  dbwt  smokes_10  \n",
       "19864    0  3544          0  \n",
       "19865    0  3487          0  \n",
       "19866    0  2642          0  \n",
       "19867    0  3220          0  \n",
       "19869    1  3118          0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_trunc.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    169577\n",
       "1      5251\n",
       "Name: smokes_10, dtype: int64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is the number of smokers and non-smokers in the dataset:\n",
    "df_trunc['smokes_10'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create treatment/control matches with the DAME-FLAME package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration number:  1\n",
      "\tNumber of matched groups formed in total:  2408\n",
      "\tUnmatched treated units:  2052 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  139666 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  0\n",
      "\tNumber of matches made in this iteration:  33110\n",
      "\tNumber of matches made so far:  33110\n",
      "\tIn this iteration, the covariates dropped are:  set()\n",
      "Iteration number:  2\n",
      "\tNumber of matched groups formed in total:  3173\n",
      "\tUnmatched treated units:  1222 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  136968 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  608255.8314184071\n",
      "\tNumber of matches made in this iteration:  3528\n",
      "\tNumber of matches made so far:  36638\n",
      "\tIn this iteration, the covariates dropped are:  dob_wk\n",
      "Iteration number:  3\n",
      "\tNumber of matched groups formed in total:  3680\n",
      "\tUnmatched treated units:  670 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  130662 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  536873.2696106575\n",
      "\tNumber of matches made in this iteration:  10386\n",
      "\tNumber of matches made so far:  43496\n",
      "\tIn this iteration, the covariates dropped are:  fagerec11\n",
      "Iteration number:  4\n",
      "\tNumber of matched groups formed in total:  3987\n",
      "\tUnmatched treated units:  340 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  129197 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  481086.21810902137\n",
      "\tNumber of matches made in this iteration:  12181\n",
      "\tNumber of matches made so far:  45291\n",
      "\tIn this iteration, the covariates dropped are:  meduc\n",
      "Iteration number:  5\n",
      "\tNumber of matched groups formed in total:  4149\n",
      "\tUnmatched treated units:  160 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  128132 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  444174.4633542034\n",
      "\tNumber of matches made in this iteration:  13426\n",
      "\tNumber of matches made so far:  46536\n",
      "\tIn this iteration, the covariates dropped are:  mager9\n",
      "Iteration number:  6\n",
      "\tNumber of matched groups formed in total:  4224\n",
      "\tUnmatched treated units:  72 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  127258 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  428300.65459263563\n",
      "\tNumber of matches made in this iteration:  14388\n",
      "\tNumber of matches made so far:  47498\n",
      "\tIn this iteration, the covariates dropped are:  tbo_rec\n",
      "Iteration number:  7\n",
      "\tNumber of matched groups formed in total:  4268\n",
      "\tUnmatched treated units:  23 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  125200 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  413828.7703970227\n",
      "\tNumber of matches made in this iteration:  16495\n",
      "\tNumber of matches made so far:  49605\n",
      "\tIn this iteration, the covariates dropped are:  fbrace\n",
      "Iteration number:  8\n",
      "\tNumber of matched groups formed in total:  4286\n",
      "\tUnmatched treated units:  5 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  124668 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  405036.25523897505\n",
      "\tNumber of matches made in this iteration:  17045\n",
      "\tNumber of matches made so far:  50155\n",
      "\tIn this iteration, the covariates dropped are:  lbo_rec\n",
      "Iteration number:  9\n",
      "\tNumber of matched groups formed in total:  4289\n",
      "\tUnmatched treated units:  2 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  124360 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  402658.87978245097\n",
      "\tNumber of matches made in this iteration:  17356\n",
      "\tNumber of matches made so far:  50466\n",
      "\tIn this iteration, the covariates dropped are:  rf_diab\n",
      "Iteration number:  10\n",
      "\tNumber of matched groups formed in total:  4290\n",
      "\tUnmatched treated units:  1 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  124355 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  401409.0781930627\n",
      "\tNumber of matches made in this iteration:  17362\n",
      "\tNumber of matches made so far:  50472\n",
      "\tIn this iteration, the covariates dropped are:  rf_cesar\n",
      "Iteration number:  11\n",
      "\tNumber of matched groups formed in total:  4290\n",
      "\tUnmatched treated units:  1 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  124355 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  401412.0607323472\n",
      "\tNumber of matches made in this iteration:  17362\n",
      "\tNumber of matches made so far:  50472\n",
      "\tIn this iteration, the covariates dropped are:  bfacil3\n",
      "Iteration number:  12\n",
      "\tNumber of matched groups formed in total:  4290\n",
      "\tUnmatched treated units:  1 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  124355 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  401422.4167158159\n",
      "\tNumber of matches made in this iteration:  17362\n",
      "\tNumber of matches made so far:  50472\n",
      "\tIn this iteration, the covariates dropped are:  urf_chyper\n",
      "Iteration number:  13\n",
      "\tNumber of matched groups formed in total:  4291\n",
      "\tUnmatched treated units:  0 out of a total of  5251 treated units\n",
      "\tUnmatched control units:  124348 out of a total of  169577 control units\n",
      "\tPredictive error of covariates chosen this iteration:  401144.5790765865\n",
      "\tNumber of matches made in this iteration:  17370\n",
      "\tNumber of matches made so far:  50480\n",
      "\tIn this iteration, the covariates dropped are:  rf_phyp\n",
      "50480 units matched. We finished with no more treated units to match\n"
     ]
    }
   ],
   "source": [
    "# extreme smoking is the treatment column, birthweight is the outcome\n",
    "model_flame = dame_flame.matching.FLAME(repeats=False, verbose=3, adaptive_weights='decisiontreeCV')\n",
    "model_flame.fit(holdout_data=df_trunc, treatment_column_name='smokes_10', outcome_column_name='dbwt')\n",
    "result_flame = model_flame.predict(df_trunc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot and Discuss Outcome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The covariate dropping order, shown above using the `verbose=3` parameter, is something we can try to use to interpret covariate importance. The FLAME algorithm learns the importance of each covariate on the outcome when choosing the best covariate to drop. Thus, the earlier dropped covariates are less relevant to the outcome. In this case, FLAME would lead us to believe that some factors, such as mother's education and the day of the week of the infant birth, are less predictive of an infant's birth weight than other factors, such as the sex of the infant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-219.36514154344763"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dame_flame.utils.post_processing.ATE(model_flame)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, as we can see, the ATE is highly negative, implying that maternal extreme smoking causes an infant birthweight to be lower than not-smoking by this many grams. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "group_size_treated = []\n",
    "group_size_overall = []\n",
    "cate_of_group = []\n",
    "for group in model_flame.units_per_group:\n",
    "    \n",
    "    # find len of just treated units\n",
    "    df_mmg = df.loc[group]\n",
    "    treated = df_mmg.loc[df_mmg[\"smokes_10\"] == 1]\n",
    "    \n",
    "    group_size_treated.append(len(treated))\n",
    "    group_size_overall.append(len(group))\n",
    "    \n",
    "    cate_of_group.append(dame_flame.utils.post_processing.CATE(model_flame, group[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEaCAYAAADDgSq4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABUnUlEQVR4nO2dd5wdV3X4v2dmXt2qslpb1V22EQSIgm0gRmAMDgkxSQgdbDA4JJAAKbQktACBX+ihhRZMNSUQGzAQYxDVBmyKLdvIlm3Z0kpeldX212bm/P64s9LTastb6c3W8/183ue9uVPumXkz98y959xzRFUxDMMwjDTw5loAwzAMY/FiSsYwDMNIDVMyhmEYRmqYkjEMwzBSw5SMYRiGkRqmZAzDMIzUMCUzDSKyRURURFbOVZ0zlUFEdorIP6Yr5eJARM4WkRtFpCwiO+danulI834UkW+KyKebfdz5jIh4IvJfInIwua5b5lqmxcaSVTIicq2I3DDJunOSG+5JwM+Ak4GDsyrg0cxUhj8APpyeOEcQkW4Reb+I3CMiFRHpEZFvi8hTJtj270UkEpG31ZVtTa71ZJ+d02x39QmewluBUeBs3HWbNUTkchEZns06jWN4CvBC4Km4Z+xncysOiMibJrnXn1a3flsDx8mKyH4RGRaRjgnWjz1T/zrBui8l6z5YV/bpSeS6aSo5gobOenHySeDrInKKqu4ct+4K4H7ge6oaAw/OtnD1qGp1JjKo6v4UxTmMiJwC/BQYAl4H/Bb34nIR8FFg/bhdrgDeAVwuIm9Q1Qj4cyCbrF8O3A78BUce9qhu//8GXj/umKUTPI0zgGsmuAeMJiAiARDp/J31fQawV1UnVS4ikk2ewdlkO7BlXNmhGR7jacB9wADwHOAjE2yzC/c8vnXsPxKRFcClybrxfA94/riyqa+Nqi7JD07B7gXePK48A/QCb0iWtwAKrEyWO4DPAvuAMnAv8Mq6/RV4+rhj7gT+sW7574FbgRGgB/gE0Fm3fnyd45d3JsvjP6dMUp8CVwJfSeq8F3jeOBnPA36VnNOvcW94CmyZ4hpel8jfOsG6znHLFyTXNQPsAP5kgn1WTlYnsBX44Az/Yw/4V9zDUgFuAy4dd13qP2+a5DifBr4JvAan7AdwytID3pTcCw8Crxm336T/c91/ekz9OKX7dtyLTiX5v/5u3H4XAT/H9cJuBh45ru5HAz9M1vfgGpj2uvXF5LyGk//l9ck5fnqaa/oi4IHkuN8A/gbQuvVvArYBlwP34F4SWnEvHF/HvZAMAV8D1o7fb1xdlwPDExz7xYkMJeB/SZ6LZJuHAjcAg8m5/RZ4/BT/a/3131l3r30EeBewH/hlUn5hcs3LyTV7L5Add49+BHg30Jfs+wogB3wI6E/kfv401/iYazGT9XXbfRf4W5xSuHmSZ+q/knN5fF35K5J7Zyt1z1xyvb45k2dQVZfucJmqhsBVOC1efx2eimvs/nuSXd+Ku5H/BNiIe+h6Zlh9DLwSeAjuDeNRwH/OYP8/wHXtxz7fBH6Hu1km4w3ANcDvAV8CPiUi6wFEpLXuGL8PvBr4j6kEEJHlwCXAh1T1mCEfVe0fV/Ri4GpVrQGfS5bT5hXAP+GUw0NxjdzXROThyfqTcW+M705+v2uKY10InIpr5F+Ku0bX4RqQx+Ie/HeIyO/X7TPV//yzZN0oR/7HsfqvAl6AU1Ln4HqA/ePk+XfgtcAjccOonxcRARCRhwL/B1yL+7//HHg48Km6/d8FXIzrNV4EPCI5x0kRkQtwivJDyfGuBd48waanJuf7l0n9Vdy91w08PvmsBv53TOYZcArwPNyb9hOBM8ed1xdwL4+PSmR8E04pTMQrgLcAu3HXv3649HmAAH8IvEBE1gDfxr2APQL3nzwb9z/U81ycEj0P9yLyPpwivAvYjPtvPyEiJzd8xseBiGzA3atX4xT62SLyexNsWgM+g2vHxngRbqSnOcxUKy2mD+4GVeBJdWXfAr5dt7yFo3sR1wKfmuKY0/ZkJtjnEtwbqzdJnUctj9v3NcAB4PTJ6kv2/fe65QDXuD0vWf4r3JtXoW6b5zBFTwb3ECvwZw1c51bcg7c5WT4d1/CcNG676XoyVdzbaf3nb6aot4ekRzruOJ+rW97GJD2Yum0+jesN+XVlNwO/PcH/+XLq3tTH3ZOXTHKMsXvhyXVlj0nK1ibLnwE+OW6/hyfbrEr+jwrw3HH/UT9T9GSALwLfGVf2MY7tydSA7rqyi3E9mlPqyk7DKeEn1u3XSE8mAtbXlT02Oa8zk+VB4LLp7sm6/f+RpAcz7h65dVzZ24C7x/67OvkqQLFuvxvr1guuN3NtXVkGdx8/fQqZxs6z/j6/fdz6KXsyyTbfrFv+DONGAhJ5P4h7kRkB2nGKcADX093KsT2ZkGOfwXdOJcuS7ckAqOrduG7hiwBEZDXwZKbW4h8BnikivxWRd4nI42Zar4g8QUSuF5HdIjI2dJAFTprhcZ6Ke5P8C1W9Z5rNbx37oa4Xtx/X4IAzem9T1Xr7xs+nq34Goj4L2K2qNyf13wP8ErhsBscA1wN7+LjP5ycUTqQd97b803GrfgKcO8N6Ae5QZ0MaoxenoBhXNnZNj/d/fgSu8f3BNPLcWvd7T/I9VvfvA89LDL7DiXPB2HU4PflkgRvHDqCuN3rbNHWeDfxiXNlE98luVa3vVZ8D7NE6u5eq3pvIPdP/okdVHxhXf5zUAfAeXE/h+yLyzyJy9gyPP8Yt45bPAW5SZ6Md4ye463hGXVn9c6a4odTb6spqONvKKqbmHo6+z49xpJmMZGTmhbhh/TE+CzxXRPLjt1fVO3HDis/G9dCuVtXRSQ7/I459Bqcc9VjSSibhk8DTkuGfy3Fv9NdMtrGqfhvYgBtuWAl8S0Tqh9aUYxvgzNiPpBv7LeBO3HDC73Okq5qlQURkE66BfZmq/rCBXWrjlpUT+//vTo5xznQb4obGNopIOPbB2WiumGGdA6q6Y9xnYIbHACf3TJno+k16TZv1Pzcoz9j5eHXfn+DohuD3cL2k3zSh7ukYmcG2Y7LHTPHcNHww1TfhFNf/4uxSt4rIi6baZxKO5xxghvfJFFTH3ef3z0CeJ+FsYJ+ve96+DXTihkcn4lO4EY1nc/Tw43hGJ3gGD0wljCkZ+CpuzPZ5uEbgM8nbxqSo6gFV/ayqXo5rKC8TkVyyej9ufBdwLr71y7juaBZ4lareqKp34d64GyaZI/EN4OOq2oyx098Bm0SkUFf2qKl2UNU+nGHx5YlNZ7yMncn3Q3Dj00/i6EbvPOAUEZnSDnC8qOog7k35MeNWPRa4I406x9HI/1wF/HFlv8E9l48/gbp/BTxkgsZgR9JbvQfX8J0/toOItACbpjnu7zjWzXvK+yThTmB14o04Vt9puOsx9l/sB7rH2WgePsGx1ojIunH1e0kdgBuhUNUPqOof414im2H/uxM4f5z99rG4/3C6UYTZ5gpcr/nh4z4fZ/IXuy8BZ+F6odONYsyIpezCDICqlkTkC7gxzGVMY/ASkbfgHuLbcdfvz4F7VbWSbPJ94GUi8jPcuOrbOdrweDfuoXiliHwN96C/coZi/w/O3vBuEakfetk/bkinUb6Ac2j4uIi8Hffwj7kKT/XW/zLcMMzNia/9rbi30cfjXJrX4x7wX6vq98bvLG6e0otxXfBGKI47X3BvfH2TbP8fwFtE5G7c8MfzcIbcRzZY34nQyP+8E8iLyMU4g/Koqt4lIl/GDfm8AnevrcXZMz5LY7wTuElEPorzHhrCDXU9VVX/SlWHReSTwDtFZD9OGb+BYxXeeD4A/ERE/gnXU7gQ+LMG5Pke7t74fHJO4BwgfoV7XsCN/y8HXp/MfdoCPH2CY5WAq0Tk74ECzlX+W6p6d/KS9C6cF+VOnKPBY5l+6LcRPoz7/z4sIu/H2ZTegbNZTDa0lCb5OgeWMUZxQ3F/Cvylqh41nJv85zeKyOnjh9dVdShxbpiu/chN8AxGOsW0CevJOD6BUzA/S8Ynp6KCMwL+FtfAtuE80sb4B5zL6VZcL+kTuHFZAFT1VpxXy9/j3uJejDM+zoQLcW/oPThPmrHPuql2mgxVHUrO4SG4xu4/cEoXJvfMGRtXfyRwPa5huxXXaPwpcKWIZHEN+1cnOcRXgKdPNFFsEl7I0ee7F+eIMRkfSM7l/+HsJ3+Gs1/9tsH6jptG/md1czM+ijOo78d5rIHzLPtCIv/vcAbXRq/RWN0X4jyxfoi7V/+do70P/xFn9/l68r2NaZS9qt4IvAT4O9x//TTc/z7pPZLspzhvsP1JXT/AuXw/LVk3Zhf4a5yr/a04Z4G3T3C4nTiPqW/g7rV7cfcFuAZyGe56bU/O7Ubcf3BCqGoP8Ec4m9lvcENKX+TYeVuzxem4Z7X+8wWcu3IFN8ownl/gHFgm7M2o6oBO4Ck6jidy7DP466l2kOQ/NoyjEJFLcQ/pqunGXI2li4i8F+ch9tBZqOtNOK+s6Yb1jHnEkh8uMxwichnurXAXbmz+fcA3TMEY9SRDZdfjXFefiJszNFdv88YCwJSMMUY3zh36ZNxQxrdwc3AMo57NuKG2DlzIktcB759TiYx5jQ2XGYZhGKlhhn/DMAwjNWy4DFi5cqWecsopcy2GYRjGguKWW245oKpdU23TkJJJXFH/BTcbdDXOr/5q4G2qOqX74kLglFNO4eabb55rMQzDMBYUIjJtJIJGezIfwUUc/jtc+PENOI+SNRwdvdMwDMMwDtOoknkaLspvf7J8h4j8HJcXxJSMYRiGMSGNGv4fxIV+rqeAm+1pGIZhGBPSaE/ms8B3ROQ/cQl+1uHiVn1GRJ4wtpGqfn+S/Q3DMIwlSKNK5q+S7/Eze1+afMAFUjytGUIZhmEYi4OGlIyqnpqmEEno7s/gZp0r8DFVfX+S4+VLuEB/O4FnqOqhJBz4+3GJfEaBy1X1V8mxLsN5wgG8VVWvSkPm3l5l2/aQvsGI5e0+mzYGdHfPNJOsYRjG4ma+zJMJgX9Q1V+JSBtwi4hcj0sidoOqvkNEXovLaf4aXDTUM5PPeTjvt/MSpfRGXOgLTY5zraoeaqawvb3KNTeMMhCPEGpE0O9zz+4WLr2oaIrGMAyjjoYM/yKyS0QemOjTDCFUde9YTyQJO38nzj36UmCsJ3IVzsuNpPwz6rgJ6BSRk3Gpk69X1b5EsVyPy6veVH58U41dw/3ku0bo2lAi3zXCruF+fnzTlLnODMMwlhyN9mSeN275ZFyujKubKw4k2fMegUs01K2qYx5sD+KG08ApoF11u+1OyiYrn6ieK3G5K1i/fv2MZLztrhpdD61SLLq4b8Wi0rW+ym231Xh6UzLrGoZhLA4atckck0NeRLYC36GJEViTNL7/A7xSVQfrM7GqqopI06J5qurHgI8BbN68eUbHVQUi4aikkZFgsUYNwzCO5kQCZFaApjkEiEgGp2A+r6pfS4p7k2Ewku+xDJM9HJ0Fcm1SNll5U3noWRn27SwyOuShCqNDHvt2FnnoWZlmV2UYhrGgaTR22VvGFRVxnl3fboYQibfYJ4E7VfU9dauuBS7D5dK+DLimrvzlSS7w84ABVd0rIt8F3i4iy5LtnoTLd9FULrwgQ9/1LQw86DMYKxlPWNeR58ILTMkYhmHU06hNZnzu+BHgPbhJms3gMbjc1LeJyG+SstfjlMuXReQKXMy0ZyTrrsMpuR04F+YXAqhqn4j8G/DLZLu3qGpfk2Q8THe38LSLi2zbnjUXZsMwjClo1CbzwjSFUNWfAJO10BdNsL3iIg5MdKxPAZ9qnnSGYRjG8dLwPBkR2QK8AOet1QN8VlV/kI5Y8xubJ2MYhtEYjc6TeTHwZZwb8ddwgTG/KCIvSVG2eYvNkzEMw2iMRnsyrwYuVtXfjhWIyJdw3mAfT0Ow+YzNkzEMw2iMRl2YVwB3jCvbDixvrjgLgyPzZOqweTKGYRjH0GhP5ifAe0TkNao6KiItwL8DP0tPtPnLQ8/K8NNbi0SZCqEqgQh+LcdjzjEXZsMwjHoaVTIvxYWQGRCRPlwP5mfAs9MSbD5zzhkB3/t5gCwbxcvF1Coe1cGAc86YL/FGTwyLMG0YRrOYdrhMRHxc1OMn42b4PxU4VVUfp6p7UpZvXtLbF7H5vJi8n2X4YIa8n2XzeTG9fdFci3bC9PYqW28eJWobYs3ZI0RtQ2y9eZTeXhsLNAxj5kz76q2qkYi8J5l/sjv5LGnu2xWxT8ucvqnGpoJSKoXs2x8jI63Awh4y27Y9ZNnqCm3tTqm47wrbtmfp7l7Y52YYxuzTqOH/GyLy1FQlWUD09ytBIaJYVEScd1lQiOjvX/hv+32DEa1tR59Ha5vSN7jwe2mGYcw+jRoR8sBXReRGXCj9w62Qqr4gDcHmMx3tcM/uPD27S0TE+HgUyXP6IvC1W97uMzwkh3syAMNDwvJ2fw6lMgxjodKoktmWfAygs82nvFvQjE+sHp4I5ZrQuSG9hni2jPGbNgZsvTkHVGhtU4aHhEN7cmzZvDicGgzDmF0ajV325rQFWUiIgF+oQb4CSU/GL3tISg5YY8b4ZasrrFntGv6tN+fYsrn5YWy6u4Utm13wz54ep9C2bDbvMsMwjo8plYyIZIF1qnpPsvxcoP51/SuqWkpRvnlJ/3BMtQJ7dmUYLQvFvLJ6lStPg9k2xnd3S3JcM/QbhnFiTNeTeQUu8dcrkuWPAb9KfncDXcC70xFt/rLrgZj7H6xRyQ9TbgmJJeD+Pa0sj9NRMn2DEWtWH2uM7+mJMEVgGMZ8Zjol82zgmXXLVVX9QwARORM3QXPJKZmdD8Tc11el2F3Fz8WUyzGjvVWWVdNRMmaMNwxjoTKdC/M6Vb27bvk7Yz+S8vWpSDXPeWBvSGFFmWxLiJeJyLa45Qf2hqnUt2ljwKE9OYYGXXy0oUFnjN+00YzxhmHMb6ZTMjkR6RxbUNXDYWSS8lw6Ys1vKjVFNWJ4yGPwoM/wkIdqRKWWzjyZMWO8P9RGz+9a8IfaUjH6G4ZhNJvpXoVvBv6ciTNN/gVH7DNLipa8z933tZFbdxDJ19ByhsquTs5sSW/4yozxhmEsRKZTMv8OfEVE2nHJyh4ETsYpnjcCz0hXvPlJPiuIVtBqhrgaIAhSrJDPWs/CMAyjnimVjKpen2S/fDdHG/h7gCtV9f/SFG6+Uo1iCsWYchwRqeKLUCjGVFMy/BuGYSxUGgmQ+RVcb2YjsBI4CGxXXbopuoaGlHIeyAYQAx6UqzBUXrKXxDAMY0Iadk9S1e24bJhLnuEhYaScpX3jANnWCtFwjsHtXQzXbLjMMAyjHvOBPQ4qcYSXq3Fo+wriSPB8JZOrUalYpGLDMIx6TMkcB9WyUmuLyZw8iAYREvrUDrZSteEywzCMo5h0noyIvLzu9xmzI87CYHREiLyQMBSiUsZ9eyGjIzZcZhiGUc9UkzHfVvd7Sc6HmYzhkZh4JEemtUxu9SEyrWXikRzDI+ZdZhiGUc9Uw2X3isi7gduBjIi8aKKNkrTMS4pIlUL7KOWeZWg1g2RrZNpHqSxdhzvDMIwJmUrJPBN4NS5IZgZ4/gTbKBNHA1jcxEptoEhhwwGkpYqOZCndvxJiUzKGYRj1TKpkVPUu4MUAInKDql40a1LNc7IZH69QJa5lkIEAVcErVMlmLCqyYRhGPY1mxrxIRALg0cAaYDdwo6qmE3Z4npPJQRxESLaGBDGEHgRZMksyXKhhGMbkNKRkktn+3wQKwC5gHVAWkaeq6p0pyjcvKVdjsqrUDrSjKogoEFK2sDKGYRhHMV2o/zE+gsuKuU5VL1DVtcBHgQ+nJtk8pjYK1b2daCTgxWgkVPd2Uhuda8kMwzDmF41Oxnw4cPG4eGXvA/652QItBKIQvKEs1dgDXyBSopGAaEkOHhqGYUxOo0pmD/A44Pt1ZX+YlC89Ash0lchvOIjkQ7QcUL5/BbXB4lxLZhiGMa9odLjs9cC1InK1iLxTRK4Grk3Km4KIfEpE9onItrqy5SJyvYjcnXwvS8pFRD4gIjtE5FYReWTdPpcl298tIpc1S756vAxk1/QRtJcJciFBe5nsmj48yydmGIZxFA0pGVW9FngksA1oS75/X1WvaaIsnwYuGVf2WuAGVT0TuCFZBvgj4MzkcyXOZoSILMclUzsPeBTwxjHF1ExEFK9Qwy9U8VvL+IUqXqGWOAAYhmEYY8wk1P9dwFvTEkRVfyQip4wrvhTYkvy+CtgKvCYp/0xiI7pJRDpF5ORk2+tVtQ9ARK7HKa4vNlPWOBS0nEFjgchDY7cchxa7zDAMo575HoW5W1X3Jr8fBLqT32twrtRj7E7KJis/BhG5EtcLYv369TMSSiOo9CxHsjFevkpczlLpWY5apH/DMIyjmO9K5jCqqtLE8ShV/RjOLZvNmzfP6LgSgGpE+b4VEPvgRahGyIK5moZhGLPDfG8We0XkZFXdmwyH7UvKe3ATQsdYm5T1cGR4bax8a7OF0hDi4Txeaw3xYlSFeDhPmvEPenuVbdtD+gYjlrf7bNoY0N1tw3OGYcxvGjL8i8gHJil/X1OlOZZrgTEPscuAa+rKX5B4mZ0PDCTDat8FniQiyxKD/5OSsqYSh6A1H2I5/NGaT5ySkuntVbbePErUNsSas0eI2obYevMovb3maGAYxvymURfmyycpnygy83EhIl8EbgQ2ishuEbkCeAdwsYjcDTwxWQa4DrgX2AF8HPgbgMTg/2/AL5PPW8acAJqJF4DXXsZvq+B3uG+vvYyXUr9w2/aQZasrtLUrItDWrixbXWHbdpv9aRjG/GbKZrEuh0wwQT6Z04ADzRJEVZ89yapjoj8nXmUvm+Q4nyLl9APiQVCokeksIX6ERj7UBGlUZc+QvsGINauP7rW0tik9PREuC4NhGMb8ZLp377GeSpajey0K9HJkKGtp4SnBihGyXUNIEKGhj8YKXjrDV8vbfYaHhLb2I8cfHhKWt1tqAcMw5jdTKhlVfTyAiLxVVf9ldkSa/4gIxFA70IrGHuLFECflKbBpY8DWm3NAhdY2ZXhIOLQnx5bN891vwzCMpU6jAzw/EpGz6gtEZKOIXJyCTPOeOISwrw2/pUrmpAH8liphX1tqhv/ubmHL5iL+UBs9v2vBH2pjy+aieZcZhjHvafRV+EPAhePKhpLys47dfHEjANkQrXnocB5CgWxImk1+d7fQ3Z3BbDCGYSwkGlUyq+pm3o+xFzipyfIsCCRQ4hqUHlhxeLhMvRoSpOdSvG1bzHd+WOHBAyEnrQy45HE5Nm1KydPAMAyjSTTaSt0rIk8YV7YFuK+54iwM4lDQSha/rUxm2Qh+WxmtZFOLXbZtW8wHPtvPXQMHGSkc4q6Bg3zgs/1s22aZOA3DmN802pN5E/A1EfkkcA9wOvDC5LP0iMHPxuRWDeK1VIhHcsTDOUipzf/KN8r0VgYpDw9S1ZCsBOQrIV/5Rp5Nm5qfw8aiCxiG0SwaDfV/DW72fAvwx8n3k5sc6n/BIBnw2kYJh/NUezsJh/N4baNISuaSX/ymykHvEMGKYZatLRGsGOagd4hf/Kba9LosuoBhGM1kJqH+fwH8IkVZFhCKl40I2suIpy7Uf5THTR9qPgf6ldxpNUqjMDQoBAGQq3Ggv/n11UcXAJLvCtu2ZxPHA8MwjMZpSMmISA54A/BsYIWqdojIk4CzVPWDaQo4H9FQqB5sQzzwAiVOljUlm0xL1mPnrW1I5wjqR0jko/0tnJJtvuHfogsYhtFMGu3JvBeXl+W5wLeTstuT8iWnZMRTCD3iSoBqhFZ9CD1XngKB5zO0p0AQRmgQI6FHeKBAsLb5M/6Xt/vs2insOxAzNAxtrbBqpcdKiy5gGMZx0Oir8J8Bz1HVG0nM26rawyQJwRY9IniFMuCh5SzgueWUZvwfOBTitVTwszFB4JwOvJYKBw41f/Zn93Kfm26C0cwIXWcOMpoZ4aabXLlhGMZMabQnUx2/rYh0AQebLtECQEMlKuWQbBlPII4gKuXQMJ2ezKGhkOI5Q9A+fHi4LNPWyqE7VzS9rt6+iHM3xdx+T4bbf+uxvMPnIZtievsiNjX8TmIYhuFoVMl8BbhKRF4FkCQQex9wdUpyzW9UkFjwPBA/wos8JBbQdHoytUjwsiUKnWUyGaVWE0plnzhqfn337YrYp2VO31RjU0EplUL27Y+RkVbMJmMYxkxp9NX09biJl7cBncDdwB7gzemINc/xwWupku3uJ7e+j2x3P15LFVIaUcqIy7wZDmcoD/uEwxni4TyZFIbn+vuVoBBRLLrcNcWiW+5PwZPNMIzFz6Q9GRF5eZ3n2HpVfRXwqmSY7ECS02VJ4mUg2zWEl40hErxs7JZTetFvbRP27+ugdrCFWBRPBYkClrU1X8l0tENvf57RQolCa0xp2KPWn2d5e9OrMgxjCTBVT+Ztdb9/NfZDVfcvZQUDgCpx7KE1//Anjj1I6bK05n3C0awLwlmsQDYkHM3Smm9+1+m0DQFnbsjAYCv77mqHwVbO3JDhtA2WVsAwjJkzVctxr4i8G+eqnJkgMyZwOBPlkiKuCeGD7fjBISQXoRWf8MF24lo6NplKWQm9EtFgDo3zbgKoV6JSbr5S27Qx4JobcpTjESQXUQ59BvYV2HKRKRnDMGbOVC3HM4FX4yZgZjg6M+YYSsqpjucjGkM4koW9nUigaCiEI1k0pdhlu/bWoEsIugbxsiFxNaC2r82Vp4DGoFXXO1NP0Gwq1RiGsQSYSskUVPXFACJyg6peNEsyzXskUPAjJIiQsRErP0ot1P/QCLSeM0RudX9dumcY2rmy6XVt2x7SeVKF6gElHIZCUelcaWFlDMM4PqZSMj8Gxsy962dBlgWDiBC0lRERtOYjfnR4ORX8GL+tjBdESDZCAb+tDH7zu0733h/SW63RcXKJVaud4f/u+wsMZkMuMhfmeY9F0DbmG1MpmX4R+RPgDuAkETkVjk3+qKr3piXcfEVjJS7nyW44iF+oEpWylO9fgcbp9GRU3RycuJKBStLQq7jyJjMwCOV8hf77hNHRgGJRyecrDAw2vSqjyYxF0F62usKa1crwkLD15pyl6jbmlKmUzCtwEy434LzQ7plgGyW12SHzF/EEyVaoPdhOTX2QCMlWEC+dBzmuQXXvMoJiL9JSRUeyVPcuI07BJCPAHdsEf/kQXi5mf59H1NfOY85ofl1j2Nt3c7AI2sZ8ZFIXZlX9uqqeoaoZYFRVvQk+S07BAGikRMMF8GMkVwM/JhouoFFKATJFqPXnGL1nFaO3r2b0nlXU+nMEKQzPDQ4p2Y4S2dYILxORbY3IdpQYHErn3Cx/TfPoG4xobTs2gnbfYDRHEhlG42Flmh8kayGjAhK7137FfUucWliZXF6ohVDra0UQFEXjkFy++fWNlGL8IJkDhKB4+L7HSCkd17lt20PIVdi+fSzqs7LKHA2Oi+XtPsNDcrgnAzA8JCy3CNrGHDKlkhGRX6vqI1S1miz/t6q+sG79PlVdlbaQ8w4fgpaQXNcwkgnRWoDW/NQGDsNajJdXcifvx8tFxBWfyt5Owr7mN/yqUB4RhgahUhNyGWjz03NjngtHg8U6PLdpY8DWm3NAhdY2Z5M5tCfHls02x8mYO6a7+8aPxF86brnQRFkWDJ4PwYoh/GLFDThmIoIVQ3j+SanUFyoU2kbRco5oVMBTgrZRStrZ/MoiYe8eoXhyTLYlJg499u4RTlmXTiM8MAgHaxW23S0MDgW0tyknn1Qhn1InZjEbx7u7hS2bi2zbnqWnxynQLZsXhwI1Fi7TKZnxA+Pj79YlOnCuEHtIJoJsBFUfYo+0LoeqEpUzSCEETyGCqJwhjeg+o2UljoXBgz6hCIF6eLEwmkJ0AYCRIeWW7RGrzu3n5I0hAwcCbvltJ50b06lvsRvHu7slOY+Ffy7G4mCm/eglqlSORmNBqwHhaBapxmjoodUAjVPyLosEP6vk1x9C8lW0nKV038pUQv0fGFAKq4ap+RWoCdmMkokiDgw0f+InwN79MaduKuMVPEb7M7S1CC2byuzdk44NyNJLG8bsMp2SyYnIW+qWC+OWl2bAEVWikSx+Pu8mREYe0Ug2tQCZng/SUiEazUIp4/LZtFTwUrABVUoxByse0grqRVQqHjrs0Ran0+iLH1PMBLR1xmQLMdWSx1AYpDLRFCy9tGHMNtMpmS8A6+qWrx63/MWmS7QA0FiISwG1oSRYZbKcVk8GT4H4iCeb4pa95iu1Usnj0HCADmSJUTwEKQd0t6aTFfP0DRnuGc3wm1/A4HBEe6vP6WsynL4hnV5F93Kfb/4vnHTOCF1nRvQd8LnppjaueFo6SmaxOhkYRqNMqWTqPcmMI4wFyJRAwIsh9ghHMqkFyNRQqO5roVbL4TwNYjRTQcPmN1a9+yJGqgGZ7lHIhES1gFp/QO9oOnMtTjs5w+c+maP11GHaVkWMDvrcuLWVi69IybOsL+KCxyj7+/Mc3BXT1uJxwWM0lfTSvb3KNTeMMhCPEGpE0O9zz+4WLr1o4TsZjGFK1JgOS9p+PCh4WSVoL5FZPkrQXsLLamoWq6iiSJinuGkP7X/4O4qb9iBhnqjS/Ar39YV4bWWClhrZYkTQUsNrK7OvL2x6XQC/3V5lzYaQTFhg9MEWMmGBNRtCfru9mkp9fYMRLa1HX7eW1nQmLP74phq7hvvJd43QtaFEvmuEXcP9/PimdKJng2v0b/hRja98s8wNP6qlOqnVJtIajbAolYyIXCIi20Vkh4i8ttnH9/NQOG0f+fV95E4+RH59H4XT9uHnm11TUl9ByK3rRUtZqg+sREtZcut68QvNf2OsxUrh1ANkuwfxO0fIdg9SOPUAtZTist15b4315w6zdkPI6nXK2g0h688d5s5702mIJfK46dcV4uIoK9eXiIuj3PTrChI1/1G47a4aXeurR6Wy7lpf5ba70jm32W706z31RJyn3rLVFTfB1jASFt0sLRHxgQ8BFwO7gV+KyLWqekfT6vAUyYVo6OFmYLplScFGAuB5CuKT7RpA8jW0nKG8t8OVN5lMIESjGTLtFYJCRFT13XKQzhCIxsIvf5phoFahGsZkA4+OTI5TiynVpzAyELDjnjKRRPjq448EpDGlWBWIDhvRHJGk5R8y6+7Z5qlnNEJDSkZETlLVBxstn2MeBewYiw4tIlfjJpFOrmS2b4ctWxqu4Pr9NbzrKni+HrbFx5EQj+RgS/MfrusP1vBvLCH+kQdaIyEaLDS9vmsfCKn1VJDgSNgcDT0yUQ62NP+d5A13RAxWK3jZCPGcvSuu+rRnc7Cl+cb4x/bGPLwSUqkqcQyeB7ms0JoL4H3N7c28/IAyUgnJFCJ8H6IIaiWfllwA1zVfiT70YEyueOywX2XUhxXN76k9cUBRP8Kvuy+jSDgn8uFdZpcxHI3eeXdNUt603kETWQPsqlvenZQdhYhcKSI3i8jNtdrMhi9c5sjgsKF//HKzUQUNfRBFgghEXeKyFN6Ic1mXIyeuHvlozSeXTafRKFUU8RVVIY5cWgPxlVIK9iaAcgXCOKbQFtK2PKTQFhLGMeVK8+vqbBd8PEYGfA7t9xkZ8PHx6GxP51pmAiEaN3cqiiS1XmhLQQgrR+qMIrfcksIwrrFwafTV9Ji7RkTagZSa1fRR1Y8BHwPYvHmzsnVrw/tetHI/rec8SLByFPFjNPIIDxQZvvMkoq1dTZf1id37aTmnBy+vhz2Y47IwcucawibX99InDbCjtBdaq6gfIZEPw1nOKJzMjf/X0dS6AJ72kAP4G3eTOXngSGrpvR1E29eyZ2vzJ4B+8A3DPBj0QusIkcT46sFwCyeF3bz7La1NretQr/Kdeu8y8enw0vMuq9aFzDk6dlkRUqgvC2iv8utx3mVZ8y5bOjQQCX66AJm7cG1aQUQeGLd6BfNznkwPR8/lWZuUNQ0vowSrBvHEc+Fl/Ngt7+huZjWH0UiJKzkyy480xOFARyqpBfb0hoTtETqSQaMs4itCxJ7elIy5HpCvEBSqSBDh+TG1fCU1lxTxoFyJGDiUpVKBXA46ihGSgglh2/aQjlUVKgc8hoY98q3QkWKE6bmIXWZhbIzpmK4n8zxcL+Y64Pl15Qr0qur2tAQ7AX4JnJlk8uwBngU8p5kVaCxEQ0WyawaQbA2tZij1dKQ2GVMQyITEpSzxaA5EIRO68ibT1xdDLiB30gCSD9FyQOXBDleeAu2tHv0lH+9Qzg2bRUKt5NOZ0uRPrQm9e30Kq8oUWhQNhd69LZxycvOv5VxEmLZGv3nYHKDmMN1kzB8CiMhKVR2dHZFODFUNReTlwHdxrl+fUtXbm1pH7IzhGgnE7ltDLz2bjOdcb/2WMpKN0KpPOJhDU2iHq6HQ0jVEbm0/Xi4krgTEkcfIrnRSCi3v9Bip5inv9w47UWSjLMs701Eyw6NQGwmo9gmRp/ixICMBwync3ZbKurnMZqO/mKN1zzaN2mRCEbkSeDhw1MC1qr6g2UKdKKp6Ha73lQriK0hMOFRAhxVJkpjVe381tT4ByYXE5SyUk7Jc2Mhw6MzrysQEy4chFuJyACoEy4eRTDoa1Bef8oPtBCsHIYgh9CgfaMfvTifMS99ATMuKGkHRJ8bDQwiDGn0DzT8/Ae7ZIXSdMcyykyOGDvn07Ghj89qmV3WYxfr2PduN/mKP1j2b90mjSuYzwMOAbwC9qUiykFDBb6niZSJUPUQi/BZNLTOmiBIsK5E/aQDJRGjNp/xgByLNV2oaCdFwDokCNPYQLyYq+a7XlgIjozGSiQirPloNEBTJRIyMphSQ01OKq0ZZuaZKkFXCqnCgJ0R2p3AtgTPPiahohsF9AcWCcOY5ETrU9KqAxR3GxuYANY/ZVtiNKpknA6eqan/TJViIxBANtOAX+p0hvuYRDbSn5munKtQO5gj3rUMj52igfgVNQalFZaj0rCRzTg9B5yhRf5HKjm6ictOrAmBwJCKzepDc2kOHh+cqu5cxuKf5nmwA69d63FXJsn2bUq7E5HMeXcUs69c2f3ius1OoaEBXl1IoKKWSMLg/oLMznQb/xzfVuLNngNAPCSMI/BpBFLL8pgxPvzSdgOmz9UY8243+Yk5lPdsKu1El8wCQa3rtCxWB2kCWqNSFiKCqxFWZwNG7OURViA61k1t3EK8QEpcCKrtWEKUQ3ksCkKBK+YHl8MBKIEaCKpJSbIhaFJNZNYSfD0HBz4dkVg1R25WOxl69IuAnP8/Ssn6Ujq6I2qjPvgeybDmv+Sd46jqf9qjA/v4MB/a7YJxnrglY4afTUN14S42BfIUVq0fJFWMqox4H7y9y4y21VJTMbL4Rz3ajv2ljwDU35CZwP1/4QVJmW2HPZLjsGhF5P+OGy1T1+02Xap4ThyCxh5+tjgVFRstZ4rS8fAWCrmGyK0eRbIi2BETlHN6u9ubX5Sl+R4kgU5dXIC+phLAByGY8RgdyxLUjhv9oNEMxk47hf2A4pi2boe/eDso1JZ8RlhczDAw3X6lt2hiw9eY8Z55y9LyVTZvTaaj2HYzp/L0y+RZ3LvmWmJaVZfb9Nh2FPZtvxO5a5oDxc4DSa/TdJGufOPZQT9BFkj1rthV2o//Qy5Pvt48rV+C05omzMBDAa6mR6RpEsjFa9YijjrQ6MvhFJbN80M3yr/kgkFk+iF88uel1xbGHhj7SUj4SVqaUJ47TafQLOY+B3jbizlHnCxhB1F+kkEunvp27Q0oa4i8boeBFeLFPqeyxc3fz3xBme95K90rh0GCGQuuRBHCjgxm6V6ZT32y+Ec/2tdy2PeTUh1Roa4exIYqhwcVh+J9thd3QUVX11FRqX6BIBjIrh8iuHIEggtAnDj3k3nRSFGskaOwTtA8fNvyH+1pTMcZrFBMdaoUVI3gtVeKRLNGhVjRK5204E4Af+/gtNchUoZaFAZ9MSkPfe3shXneIk08t4wdKFAr77wvZm5KL9mzOWzn/kVm+d1uB+7Z5jJaVYl5YUcxx/iPTeQVf3u6zexfs7w8ZGnHDgV2dAStSeiOezWu5mA3/s62wG1ZdIpIBzgdWq+qXRKQFQFVHUpFsHiOeEkdCZU8nceTh+TGxxqlFYY5rEO7rICzW8Ao14lKGcF8HcQoR4yUr5E7uIx4qEPW1IZmQ3Ml9VPZ0Nr8yQFXJdA9AFBDXMnie4nUPoPvSUdieF6O1DGGpit8aEZY8tJbB8xZshKTDnH16wPW/KLByZQW/EBGVAuKDBc4+PZ03VJdlVDjpnDIr1rksozt+ml6W0dlkMRv+YXYVdqNRmB8KXAtUcGFavgQ8DrgMeGZq0s1T4pqgpRyZUw7itVSIR3KEO1cQ19JLvxyOBpR3LT88Kz4cDVJJv+yhaBSQX3sIr1glHs1S3r0ML6WMbLUI4nJA7MXEKngCXuxTSycRJ11dPstXBxw61MbBHqFYVE5d7eHrwm88evsinvB42HegnaFhaFsGqx5GKlk/x+o7/3zYd6CF/XdDWyucf3569c0mc2EDWqw0esU+ArxBVT8rIoeSsh8CH09HrHmOgPhJzyXyXH4ZP07Nu0xEyCwfhCjn5pJkQjLLB6nsWt70unzxiEOP6oOdaCyu1xZ6+JJOo1EpCSXJ03LaAaRYQUdzjNy7kmIpnYt54eYs3/1NJ2du6qdtRcjQwYA92zp58uaFb9XtG4xYd7ay/tQj105V6fldOkM8s13fbDIXceAWK40qmYcAn0t+K7hhMhEppCLVPMfPKl6xQq2vFY1d0EWvWMHPpjRcFinhoXYKGw4ihRpaylC6fwVxCgEy875Q7ulEV4wiQYyGHrWDRfJ+Og9XLYrJLysh4jz0RCC/rERtKJ3hqz99SpZDo+3cszPDgTuUtqLwqLMK/OlTFvY8Epj9IR4bUjIaodHX053A79cXiMijgB3NFmghEIcCEYAmQSoVoqQ8BTxfkBan1Kq7V1Dra0VaKngpNPy5FiXIOeeCuJxBaz5BLiLXklLWT18JlpfI5iCf9cnmIFhewkspRE93t/CUxxU4c3ULqzvznLm6hac8rpBKwz/b6ZA3bQw4tCfH0KDLvjk0mLhMb0xniGe26zMWJo3eDf8KfEtEPgpkReR1wEuBl6Qm2TzGzYpfRfHcHoKOYcKBVir3rUltVrzGisQeudX9h9Mvl3Z3onEKNhnxyC0vk1/T7+bkVAPKPZ14e9IZLmtr8ahGBYLiMJFXxY99MqMFsi3p1Nfbq9x+f4mHPbrCo9uU4aEat98f09XV/AmEsz2zetZdpm1IyWiARl2Yvykil+CUyg+BDcCfq+otaQo3X/FzoFQY+c06F69MFLIl/Fxzk16NIZ7gt5XAVwg98BW/rYR4zX+Y8wUl7izht9TccFlGCTpL5A+l8/b9sLMz3NKzjHyLh1+sEZUyjPR18LCz0xmimM2Gfy7cYGd7iMeGlIzpaLhfq6q/Bv4mRVkWDHGseDkonLYHKYRoKaDUs4w4hZ4F4NIu+4rfWjoc6l8GkrwyTabYAtKiVPqLRJHg+0pLi1JoaXpVADz1kizx1iIP9GYY2aO0FISzTs/w1C3p2Ej6BiMKWeXmX6rzwGqFU0+B0mDzG/7FbrMwjEZo1IU5AJ4NPIJjQ/1fmYJc8xrPF7IrhwmWlcGLIR+SrQxT3pnOhD4QiIRosODclmO3nIY724qODHfvbSPXPUSmEEHVp7S3jbUd6bypXnhBhvt6WmhpKTE06gzxp68ucOEFKb0Zhx433RLRfero4URiN91S5PfWNH94ztxgDaPxnszngIcC38ZC/SO+8y4TPwYvAiFZTqknoxCNZvEysRsyi4RoNEsaU1daix5eyTXwKupC6MQZWovpzXvoXAZndMR1gQhTq8rl5snU3HUE8BXJ1FLJzWM2C8NoXMlcAqxTTSsTxsJCI3G2GC+ZG+PFoJJazhUNoXYoT1zJMBaRMxr10RQCckaEtKyoMHqwjVrVI5ONaVlRISKd6J/btoecck6FtnYfF7ws3RhR6sec/4gcO3cHhyMjn/8In9HedFymzWZhLHUaVTK3A8sBUzK46KyV3g4kE+NlI+KqT6W3I7X0yy7qc5Zs1368YoV4NEf5vq5Uoj73DwKdStdD+pBciFYCRnZ30N/f/LpgbvKERIHHIx925IVgaNDsJIaRFo0qmecDnxCR/+PYUP+fabpU8xwB4lJArbftcJiXuBSkFoU5aIHMmv3EoznioTz4SmbNfoIDzfdmExViiSiXQaseEkMskUsxnQJzkSfE7CSGMXs0+mRdDvwhsAwo1ZUrLtfM0kIgaC8TtNQO50BRoJKSllFVRDyyJw0cjvpcPVhEtflGmVxeqY0G1Pws6sdI5MFoQC6f3gTC2Wz0zU5iGLNLo0/yK4BHqOqdaQqzUFAUYg8pVPESl2KGcq48DUJBAiXoKLl0z9WA2kABUogwoJHgtYQE4hOHLjpxmA/RkcUzoc/sJIYxezSqZHpxKZgNXMBKCKnu7TgyGdOvJuVpVKhEQ1nCQ0UkV0UrWaKhbCrzZBSQ4RYKZ+wjaK8QDuYY3rGKlEbLAGv0DWMx06iSeS/wORF5J7CvfoWq3tt0qeY5UUUJ+9spnNF7OLFXaUc3USWtyZhCrEp51wo09JAgJvaqpOF3m8vDyuUVDty/nMqoR64Ys3JNhVxf06syDGMJ0KiS+VDyfem4cmXM73QJ4WcFr32E2r5ONBI3b6Z9BD/bnkp9UQV0uIXcWb14rWXi4Tyjd3UTVZpf17I24d79QnHdITqKIbXRgAO7O9jQuXhsFrMZGdkwljqNxi5b2BmImo2neIGSWTaIeIrGQu1QIZUkYuBipXktZcKBAgzmQQWvpYyfa35dcSxk8zHVkQylwQDfd8txvDga4bHIyMtWV1iz2jkabL05x5bNzQ+QaRhGg6H+ReQDk5S/r6nSLBRU8PIRua4hcif3k+sawstHpGW48AIl2z1ApqNE0FYh01Ei2z2AFzRfqQ2MxGQ0Q74Y0bKiQr4YkdEMAyMLPz0xHB0gU8QFyFy2usK27elMNjWMpU6jPZTLJyl/fpPkWFgISBC6nosfg6duOa0X4ViIR7No6EHkoaFHPJp1McyaTLWmZDsrtAQ58rV2WoIc2c4K1VpK9qZZpm8worXt2MmffYMp5Xs2jCXOlMNlIvKise3qfo9xGnAgFanmOSIKsVLZ18bhWTKxuvIUKGaE2v4OMmfsw2srEQ8VqO1fQTHTfCXT1uoxWgvo2DBEri2kMhRwYEcbba2LY8TUIiMbxuwynU1mrKeS5ehei+Lcmi9LQ6j5joZCNNCClx1BgggNfaKBFjSlzJjtnR5DhSrV/e2wr91NBi1UaetsfsN/1mkemVKBwb0+Q/cruaywdlmWUwuLQ8nYjH/DmF2mfLJU9fEAIvJWVf2X2RFp/qMKUSmDN5J1ib1Cj6iUIYUJ+ICzyYgnZNrLkHXh92t9LanYZM5/ZJZ4ewtRXoiI8fHwy0XO35hOfpfZxmb8G8bs0qh32WEFI27GodStWxwW4RnitVTxCzXUU7yMELVUU6urVhGyxRpBIXKTJQsRUqxR621+w3jhBRn6hlsYqPrUYiXjCR0r8unld5kDbPKnYcwejSYtW42bK3Mh0Dlu9ZIbzPYySrBsBL+1dtgkExDjZdLpyoivxLUMGpWQTITWfOJahiCF/DXd3cLTLnZv+jaPxDCME6XRgej/AkaBi4Af4pTNm4Dr0hFrfqOI8+zy4sPzZIjFladANitoWwnRAC1lET8m01YiezC9eGL2pm8YRjNoVMk8GlivqiMioqr6WxG5AvgZ8PH0xJufaATRUB6/pXa4ZxEN5dGUvGDbWzweHM2TW3/wcOyyygMraG9ZHMZ4wzAWL422UhEcTo3YLyJdwAiw5kQFEJG/FJHbRSQWkc3j1r1ORHaIyHYReXJd+SVJ2Q4ReW1d+aki8vOk/Esiko61OoZwOEvtwXaqe5ZRe7CdcDgLKVmnii0Q9uUZ+O1qDv38VAZ+u5qwL0+xJZ36DMMwmkWjSubnwFOS398FvgR8Dbi5CTJsA/4c+FF9oYicCzwLeAgu/fOHRcQXER9nH/oj4Fzg2cm2AO8E3quqZwCHgCuaIN8xaKxoKYtGPlrz3Hcpi8bp2GQOHYA4DMh0lMmuGCbTUSYOAw4tyVlKhmEsJBpVMs/H2WIAXgl8H6ccnnOiAqjqnaq6fYJVlwJXq2pFVe8DdgCPSj47VPVeVa0CVwOXJl5vTwC+mux/FfC0E5VvIjxf8AplN/O+kkFDD69QxvPTsZEMjMRIxzCZ5cPkukbILB9GOoYXTagXwzAWL426MPfX/S4Bb01LoDrWADfVLe/myPDcrnHl5wErgH5VDSfY/hhE5ErgSoD169fPSLBYFfF9MisHkSRpWXVfO3FKE2VqUURh7SBBoUoMBHnn4VbpOSmV+gzDMJpFoy7MOeANwLOBFaraISJPAs5S1Q82sP/3gIlaxH9W1WtmInCzUNWPAR8D2Lx584y0g+cJmRXD5Nf24+VC4kqARh7evStSkdUPhOpgQDgaQOyDFxGHrtwwDGM+M5OkZWuA5wLfTspuT8qnVTKq+sTjkK0HWFe3vDYpY5Lyg0CniARJb6Z++6biZRS/YxQv44JkepkwWU6nJ5PFZ3DvClrO7cHvGCUaKFK+Yw3tS2+KkmEYC4xGlcyfAWckLswxgKr2iMgJe5dNwbXAF0TkPcBq4EzgF7jpj2eKyKk4JfIs4DmqqiLyA+DpODvNZUAqvSRVQUOPcDjP2GxMDT00pVD/xVYIu4bxRlvR4TY8Tyl2DVMknZ6TYRhGs2jU8F9lnEJK3JgPnqgAIvJnIrIbuAD4loh8F0BVbwe+DNwBfAd4mapGSS/l5TgvtzuBLyfbArwG+HsR2YGz0XzyROWbiLgMlZ4VLl4ZMVEpQ6VnBXE5jdqgvVNZc0pIW2dIvr1GW2fImlNC2jsXR/h9wzAWL432ZL4CXCUirwIQkZOB9+F6DCeEqn4d+Pok694GvG2C8uuYINqAqt6L8z5Llbii1HqLoCsRX9FIqO0rElfSafQ3npqlJyjQ1VUhRvEQyodyrDl1cQStNAxj8dJoT+b1wH3AbbjYZXcDe4A3pyPW/MZvFVRDKj3LDn9UQ/zWdIbLrnhmgeBAF5lqgfZclkzVLV/xzEIq9RmGYTSLaXsyyeTHfwFeq6qvSobJDqimFdh+/hNVlNxqaD33brzlw8R9rQzfsZrqwXQuyZYtPn19y/nv/8lz/6Eaq5ZleOlfFNmyxQz/hmHMb6ZVMqoaicjf4AJioqr70xZqvuPnhOJZewk6ShD5eB0limftpbKnPZX6enuVWq7C3/5jSGsbDA+FHNpTobfXt+jIhmHMaxodLvsM8NI0BVlI+FlF8lWXHxRAQfJV/Gw6PZlt20OWra7Q1q6IQFu7smx1hW3bw+l3NgzDmEOmVDIi8pjk56OA94vIThH5sYj8aOyTvojzj6gqRCN5573sxSAQjeSJqun0KvoGI1rbjlZgrW1K32BKYZ8NwzCaxHTDZd8G2nHh/JdcSP/JiMpQvqebTL6GLBtFB3OU7+kmSsmFeXm7z/CQ0NZ+RNEMDwnL280mYxjG/GY6JSMAqnrVLMiyYAgKSlyrcehnZyICquAVBwkK6QyXbdoYsPXmHFChtU0ZHhIO7cmxZXOjHuiGYRhzw7RKJplZP+k4UDI3ZUkRVYRc0af9EXfjtZeIBwuMbF9NpZJepsotm11K5J4elxJ5y2ZLiWwYxvxnOiVTxIXYn6w1U1h6AbS8LGS6+iH2iAdaQN2yd39nanVaSmTDMBYi0ymZEVVtmxVJFhASKOJHaHhEv4ofIcGSnTpkGIYxIdMpGWs1J0BrAvj4raPgKcRCWGpLyg3DMIwxppsnY63mBGgEtX3taM1H8lW05rtl8yg2DMM4iul6MufOihQLDD+reO2DlHZ2oaGHBDFe+yB+dtlci2YYhjGvmFLJqOquqdYvVbQmxMNF/FyEBCGoEA0XbbjMMAxjHI2GlTHqiGKIyzmCFUNk1/URrBgiLueI4rmWzDAMY35hSuY48POQX78fvxAikY9fCN1yfq4lMwzDmF/YlPHjwPMVrxDiZ2uMTfmPC4LnmzOeYRhGPZMqGRHZRQMuzKq6vqkSLQCk5hEOFsh0lJFMhNZ8wsECUrOOoWEYRj1T9WSeV/f7D4DLgA8A9wMbgJfjUgAsOQIfqruXQezj5WrElQzVPe0ESy72gWEYxtRMqmRU9Ydjv0XkQ8CTVbWnruzbwHeAd6cq4TzEy4DmK4R9RZxZK8bLVxCL+GIYhnEUjdpkVgPD48qGgTXNFWdhIEBmWY3C+r14hRpxKUPpgZWEvXMtmWEYxvyiUSVzLXCtiLwV2A2sA16XlC85CkUhztao9bfCIUDAz9bIFG2ejGEYRj2NKpmXAm8CPorr1ewBvgK8OR2x5jcdbR5Dnk/QVkL8GI08wmorbW1m+DcMw6inISWjqmXgtclnybNhncfdlQLZKAKJIfKoRgU2rDMlYxiGUU/D82RE5GLgWcAqVX2qiGwG2lX1+6lJN0/ZdE5AtKfA3n6PcjUmn/U4bVWOTatt2pFhGEY9Db16i8jfAh8B7gYuTIpLwFtTkmtec85pAV6UY1UXrF4Dq7rAi3Kcc5opGcMwjHoaHd95JfBEVX0HMBah63fAxjSEmu94PuDXiMtZ4tE8cTkLfs2VG4ZhGIdp9NW7DRiLyDwWBSADVJsu0QLgjh0hrR0RGpcgUrK+0OoVuGNHCOTmWjzDMIx5Q6M9mR9xrNH/74AfNFechcH9u2MGaxX8XEw2r/g5t3z/bgvDbBiGUU+jPZm/Bb4hIi8B2kRkOzAE/Elqks1jKhXlwKiQXxYhAWgI5UMZihYg0zAM4ygadWHeKyJ/gIthtgE3dPYLVV2Sr+6KoBJTCUpoECH4qGRRy1ZtGIZxFI16l12jjl+o6ldU9SZVjUXka2kLOB/xJMZrqeLlQ7xc6L5bqniyJHWuYRjGpDQ6XPb4Scq3NEmOBYUIVGoRlT05ogh8H3JBhFhHxjAM4yimVDIi8pbkZ7bu9xin4cL+nxAi8h/AU3GeavcAL1TV/mTd64ArgAj4O1X9blJ+CfB+wAc+kbhWIyKnAlcDK4BbgOeratM94Ab6YSjOkl/fh58P0XLA0APLGbAJ/4ZhGEcxXbO4Lvl4db/XAWtxdpm/bIIM1wObVPVhwF24wJuIyLm4CAMPAS4BPiwivoj4wIeAPwLOBZ6dbAvwTuC9qnoGLnTlFU2Q7xj2HVQCH7Jxlkw1TzbOEviu3DAMwzjClD0ZVX0hgIj8TFU/noYAqvp/dYs3AU9Pfl8KXK2qFeA+EdkBPCpZt0NV701kuxq4VETuBJ4APCfZ5ipcUM+PNFvm2ItoWVElrmQISxB40LKiStQXNbsqwzCMBU2j3mUfBxCRNmAlHHGjGmvsm8SLgC8lv9fglM4YuzmSv2bXuPLzcENk/aoaTrD9MYjIlcCVAOvXzyyD9Ip2n36JyXeEoAKiVAcDVrTblH/DMIx6GvUuO0dEfg0MADuSz93Jp5H9vyci2yb4XFq3zT8DIfD5GZ/FcaCqH1PVzaq6uaura0b7bjojAweWk89Ay4oK+QxwYLkrNwzDMA7TqHfZR3Cz+x8P3AecAvw78LNGdlbVJ061XkQux03svEhVxwwbPTj7zxhrkzImKT8IdIpIkPRm6rdvKpv/wKOab+GeXWsZKSktBeFRZ2fY/FCz/BuGYdTTqJL5PeBiVa2JiKjqgIj8E7AN+NyJCJB4ir0aeJyqjtatuhb4goi8B5co7UzgF7ihujMTT7IenHPAc1RVReQHOJvO1cBlwDUnIttknLYhoP3kgH0HPIaGoa0VVq30WJm1KMyGYRj1NPrqXcYFxAQ4ICLrk31XNEGGD+ICcF4vIr8RkY8CqOrtwJeBO4DvAC9T1Sjppbwc+C5wJ/DlZFuA1wB/nzgJrAA+2QT5jmHTxgAqOTZu9NiyRdi40YNKzpUbhmEYh2m0Vfwx8Azg08BXgW8DFeCEE5Yl7saTrXsb8LYJyq8Drpug/F6OeKClRne38JANBb7zQ48HD4SctDLgksfl6O622ZiGYRj1NOpd9oy6xdcDtwOtwGfSEGq+09ur3H5/iYc9usKj25ThoRq33x/T1VU0RWMYhlHHjMd3kqCYn01BlgXDtu0hy1ZXaGt3Pgruu8K27Vm6u83DzDAMY4yGlIyIdODyxzwC14M5jKo+KQW55jV9gxFrVh89u7+1TenpiThiujIMwzAa7cl8BRcn7OtAKT1xFgbL232Gh+RwTwZgeEhYbpMxDcMwjqJRJXM+sDKNYJMLkU0bA7benAMqtLYpw0PCoT05tmw27zLDMIx6GnVh/glwdpqCLCS6u4Utm4v4Q230/K4Ff6iNLZvN6G8YhjGeRl+9LweuE5GfA731K1R1fAqAJUF3tyRGfrPBGIZhTEajSuZtuFAuO4H2uvIlG9u+t1fZtj2kbzBiebvPpo2B9WQMwzDG0aiSeRZwlqruTVOYhUJvr7L15lGWra6wZrWzyWy9OWdDZoZhGONo1CZzL1BLU5CFRP08GRE3T2bZ6grbtofT72wYhrGEaLQn81ngWhH5T461yZxwaJmFhs2TMQzDaIxGlczLku+3jytX4LTmibMwsHkyhmEYjdFo7LJT0xZkIWHzZAzDMBrDsmwdBzZPxjAMozEmffUWkTtV9Zzk9y4mcVdW1fUpyTavsXkyhmEY0zPV+M5L6n4/L21BDMMwjMXHpEpGVX9St7hKVb8yfhsReXoqUhmGYRiLgkZtMpOlMf5YswQxDMMwFh9TukOJyJh7sicipwL1lu3TgHJags13LKyMYRjG9Eznc7sDZ/AX4J5x6x4E3pSCTPOe3l7lmhtGGYhHCDUi6Pe5Z3cLl15kHmaGYRj1TDlcpqqeqvrAj5Pf9Z/Vqrokh8t+fFONXcP95LtG6NpQIt81wq7hfn58k0XeMQzDqKchm4yqPq5+WUROE5FTUpFoAXDbXTW61lcpFl3ssmJR6Vpf5ba7TMkYhmHU05CSEZEvisijk98vBG4HbheRK9IUbr6iCkTjhsUiceWGYRjGYRr1LrsIuDn5/ffAE4FHAa9NQ6j5zkPPyrBvZ5HRIQ9VGB3y2LezyEPPsomZhmEY9TQabCurqlURWQMsV9WfAohId3qizV8uvCBD3/UtDDzoMxgrGU9Y15HnwgtMyRiGYdTTqJL5jYi8DtgAfAsgUTiDaQk2n+nuFp52cZFt27PmwmwYhjEFjSqZK4B/wyUu+6ek7ALg82kItRCw2GWGYRjT02io/3uA54wr+yrw1TSEMgzDMBYHUxr+ReQD45avGLf8P2kIZRiGYSwOpvMuu3zc8n+MW764eaIYhmEYi43plMx4S7ZZtg3DMIyGmU7JjJ9eaNMNDcMwjIaZzvAfiMjjOdKDGb/spybZLHLLLbccEJH751qOBlkJHJhrIVLCzm3hspjPz85tcjZMt4HoFLFQRGQn0/ReVPXUGYtlHDcicrOqbp5rOdLAzm3hspjPz87txJiyJ6Oqp6RZuWEYhrG4aTR2mWEYhmHMGFMyC4/FnMPHzm3hspjPz87tBJjSJmMYhmEYJ4L1ZAzDMIzUMCVjGIZhpIYpmQWAiKwTkR+IyB0icruIvGKuZWo2IuKLyK9F5JtzLUuzEZFOEfmqiPxORO4UkQvmWqZmISKvSu7JbUkG3fxcy3QiiMinRGSfiGyrK1suIteLyN3J97K5lPF4meTc/iO5L28Vka+LSGez6zUlszAIgX9Q1XOB84GXici5cyxTs3kFcOdcC5ES7we+o6pnA7/HIjnPJKfU3wGbVXUTbnL2s+ZWqhPm08Al48peC9ygqmcCN7BwMwJ/mmPP7Xpgk6o+DLgLeF2zKzUlswBQ1b2q+qvk9xCukVozt1I1DxFZC/wx8Im5lqXZiEgHcCHwSQBVrapq/5wK1VwCoCAiAVAE9syxPCeEqv4I6BtXfClwVfL7KuBpsylTs5jo3FT1/1Q1TBZvAtY2u15TMgsMETkFeATw8zkWpZm8D3g1EM+xHGlwKrAf+O9kOPATItIy10I1A1XtAd4FPADsBQZU9f/mVqpU6FbVvcnvB4HFmnb+RcC3m31QUzILCBFpBf4HeKWqLorU1yLyJ8A+Vb1lrmVJiQB4JPARVX0EMMLCHW45isQ2cSlOka4GWkTkeXMrVbqom/Ox6OZ9iMg/44blm57t2JTMAkFEMjgF83lV/dpcy9NEHgP8aRIn72rgCSLyubkVqansBnar6ljP86s4pbMYeCJwn6ruV9Ua8DXg0XMsUxr0isjJAMn3vjmWp6mIyOXAnwDP1RQmTpqSWQCIiODG9O9U1ffMtTzNRFVfp6prkzh5zwK+r6qL5m1YVR8EdonIxqToIuCOORSpmTwAnC8ixeQevYhF4tQwjmuBy5LflwHXzKEsTUVELsENVf+pqo6mUYcpmYXBY4Dn497yf5N8njLXQhkN87fA50XkVuDhwNvnVpzmkPTOvgr8CrgN154s6BAsIvJF4EZgo4jsTlLOvwO4WETuxvXe3jGXMh4vk5zbB4E24PqkXflo0+u1sDKGYRhGWlhPxjAMw0gNUzKGYRhGapiSMQzDMFLDlIxhGIaRGqZkDMMwjNQwJWOkjoh8WkTeOkd1i4j8t4gcEpFfzIUMix0R+baIXDb9ltMe57kishjD0ixpTMksQURkZxLyu6Wu7MUisnUOxUqLxwIXA2tV9VH1K0Tk9SIynHzKIhLVLd/eLAFE5BQR0SSI5KyTtpJX1T9S1auSui4XkZ8c53E+r6pPaq50xlxjSmbp4uPC6y8oRMSf4S4bgJ2qOjJ+haq+XVVbVbUVeClw49iyqj6krk4REXtWFglzpeyXKvbgLF3+A/jHiZIUTfTmLSJbReTFye/LReSnIvJeEekXkXtF5NFJ+a6klzR++GRlkvBpSER+KCIb6o59drKuT0S2i8gz6tZ9WkQ+IiLXicgI8PgJ5F0tItcm++8QkZck5Vfg0gdckPRO3tzoxUnO920i8lNgFDhtGjn/OImyPJhcgzfVHe5HyXd/IscFM72GIpITkXeJyAMi0isiHxWRQrJuSzKD+x+S/faKyAuTdVcCzwVendT9jQnOtZH/+ydJ/YdE5D4R+aPx24rIOcBH6653f7L+KeIS7g2JSI+I/OMk1/yoXlAi00vFJQvrF5EPiYhMsm9BRK5K5LtTRF4tIrvr1u8UkdeIi7owIiKBiPypuIRr/ck5nDOu7jPqlg/3Buuu9+tF5EBy7OdOJJdhSmYpczOwFZjwgW+A84BbgRXAF3DBLf8AOAN4HvBBcVGjx3gu8G/ASuA3JNFexQ3ZXZ8cYxUuftmH5eikbM8B3oYLfzHRUMzVuECUq4GnA28XkSeo6ic5uofyxhme4/OBK5N6908j5wjwAqATlxvnr0Xkacm6C5PvzkSOG5PlmVzDdwBn4cLSnIHLJ/SGOllPAjqS8iuAD4nIMlX9GO5a/7+k7qfO8BqMcR6wHff//T/gk+MbfFW9k6Ovd2ey6pPAX6lqG7AJ+P4M6v0T3DV5GPAM4MmTbPdG4BTgNNzw6ETx756N+286k+2+CLwS6AKuA74hItkG5ToJdy3W4OKZfUyOxKcz6jAls7R5A/C3ItJ1HPvep6r/raoR8CVgHfAWVa0kOUWquMZwjG+p6o9UtQL8M+5tdx2uEdmZHCtU1V/jok3/Zd2+16jqT1U1VtVyvRDJMR4DvEZVy6r6G1zv5QXHcU7j+bSq3p4kdbpkKjlVdauq3pbIeCuuAXvcNMdv6BomjfmVwKtUtS9JXPd2js5CWUv2ranqdcAw0MxG735V/Xgi61XAyTSeV6UGnCsi7ap6aCwBX4O8Q1X7VfUB4Ac4JTsRzwDenhx/N/CBCbb5gKruUtUS8EzcPXl9EkH6XUCBmUWR/tfkv/oh8K1EBmMcpmSWMKq6Dfgmx5ffpLfudyk53viy+p7Mrrp6h3EZ+lbjbCbnJUMW/ckQy3Nxb4rH7DsBq4GxhneM+2lO5tD6eqeUU0TOE5EfiMh+ERnAvdGvnOb4jV7DLlzWyVvq6v5OUj7GwboMh+CG+Oqv/4ny4NiPumi9jR7/L4CnAPeLGyq94HjqZepzWs3R/9dE90x92WrcfQKAqsbJ+kbvm0Pj7Hz3J8c0xmFKxngj8BKOfrjGHp5iXVl9o388rBv7kQwBLcel6t0F/FBVO+s+rar613X7ThXFdQ+wXETa6srWAz0nKO/4eqeT8wu4kPDrVLUDZ5uQCY5zPBzAKZyH1NXdkTgszPQ8JqKZ//cxdanqL1X1Utww4/8CXz7OY0/FXo5OHbxugm3qZduDe3EADqfTWMeR+2aUqa/HMjk6w+l6Fnjq6bQwJbPEUdUduKGav6sr24972J4nIr6IvAg4/QSreoqIPDYZ8/434CZV3YXrSZ0lIs8XkUzy+YN6I+w08u8Cfgb8u4jkReRhOJtEsxOfTSdnG65HVRaRR+HsSGPsx6WWPu14Kk7esj8OvFdEVgGIyBoRmcw+MZ7eqepu8v/dC6wds22ISFbc/JeOZFhqkHTSbH8ZeJ2ILBORNcDLG9j+j0XkInEJAf8BqODuJXB2w+ck1+MSJh76fHNyfn+IG/b9SjNOZLFhSsYAeAswPu/8S4B/Ag4CD+HIw3e8fAHXa+oDfp/EMJsMcz0JZ1/YgxseeSeQm8Gxn40z+u4Bvg68UVW/d4LyHkUDcv4N8BYRGcLZur5ct+8oznHhp8lw1/nHIcJrgB3ATSIyCHyPxm0un8TZRPpF5H8n2aZZ//f3gduBB0XkQFL2fGBnIvdLccOMzeYtOOeP+3DX5qs4pTEhqroddw/+J66n+FTgqapaTTZ5RVLWn8j7v+MO8SBwCHcvfB54qar+rjmnsriwfDKGYSw6ROSvgWep6nTOF8dz7C3A51R17TSbGlhPxjCMRYCInCwijxERL3El/gdcr9aYY2zmq2EYi4Es8F/AqbghrquBD8+lQIbDhssMwzCM1LDhMsMwDCM1TMkYhmEYqWFKxjAMw0gNUzKGYRhGapiSMQzDMFLj/wPCYivZ/3+IlgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(group_size_treated, cate_of_group, alpha=0.25, edgecolors='b')\n",
    "plt.axhline(y=0.0, color='r', linestyle='-')\n",
    "plt.xlabel('Number of Treatment units in group', fontsize=12)\n",
    "plt.ylabel('Estimated Treatment Effect of Group', fontsize=12)\n",
    "plt.title(\"Visualizing CATE of matched groups from FLAME\", fontsize=14)\n",
    "plt.savefig('interpretability.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot leads us to believe there is little heterogeneity in the Conditional Average Treatment Effect (CATE) of maternal smoking on infant birth weight. Most extreme values are observed in small matched groups, for which CATE estimation is difficult. Put simply, we don't see sufficient evidence that there exists some \"types of\" mothers and infants (classified by characteristics such as education and race) for whom the impact of maternal smoking on infant birthweight is different than the outcomes for other mothers and infants. "
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}