{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1f25ba4a",
   "metadata": {},
   "source": [
    "# Bonate (2011): Pharmacokinetic-Pharmacodynamic Modeling and Simulations\n",
    "\n",
    "Here, we replicate some of the examples described in Bonate (2011). I recommend following along with the 2nd edition of the book. The purpose of this notebook is to illustrate the versatility of `delicatessen` by illustrating its application for pharmacokinetic modeling. This can easily be done by using the built-in estimating equations, as will be shown.\n",
    "\n",
    "Bonate PL. (2011). Pharmacokinetic-Pharmacodynamic Modeling and Simulations. 2nd Edition. Springer, New York, NY.\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4d95e437",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NumPy version:        2.3.5\n",
      "SciPy version:        1.16.3\n",
      "Pandas version:       2.3.3\n",
      "Delicatessen version: 4.1\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import delicatessen as deli\n",
    "from delicatessen import MEstimator\n",
    "from delicatessen.estimating_equations import ee_emax, ee_glm\n",
    "from delicatessen.utilities import aggregate_efuncs\n",
    "\n",
    "np.random.seed(80950841)\n",
    "\n",
    "print(\"NumPy version:       \", np.__version__)\n",
    "print(\"SciPy version:       \", sp.__version__)\n",
    "print(\"Pandas version:      \", pd.__version__)\n",
    "print(\"Delicatessen version:\", deli.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba72185b",
   "metadata": {},
   "source": [
    "## Chapter 4: Variance Models, Weighting, and Transformations\n",
    "\n",
    "The first example comes from Chapter 4. Data comes from Table 9 (pg 153) from a study by Byers et al. (1989) on XomaZyme-791 dose on percent change in albumin concentration among 17 patients. Note that this number of patients may be below what is considered sufficient for inference with the sandwich. \n",
    "\n",
    "The E-max model is described by\n",
    "$$ R_i = E_0 + (E_{m} - E_0) \\frac{D_i}{ED_{50} + D_i} $$\n",
    "where $D$ is the dose, $R$ is the response, $E_0$ is the response at a dose of zero, $E_m$ is the maximum response, and $ED_{50}$ is the halfway maximal dose. Here, $E_0 = 0$ so the E-max model reduces to\n",
    "$$ R_i = E_{m} \\frac{D_i}{ED_{50} + D_i} $$\n",
    "This is the model we will consider.\n",
    "\n",
    "The corresponding data is presented below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bd74b455",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = [5, 5, 13, 11, 14.5, 6.8, 42.5, 37.5, 25, 38, 40, 26.5, 27.7, 27.4, 45, 61.4, 52.8]\n",
    "r = [-12, -13, -28, -24, -45, -18, -26, -40, -45, -26, -29, -26, -22, -28, -36, -27, -48]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6eeb4f8",
   "metadata": {},
   "source": [
    "In the book, a two parameter E-max model is used. The first parameter is maximum response and the second is the 50% effective dose. The minimum dose will be set as zero (i.e., not estimated) since there is no data on a dose of zero. \n",
    "\n",
    "Here, we will apply the same linear model described in the book (i.e., prior to the Box-Cox transformation). This can be done using the `ee_emax` built-in estimating equation. But before that, we will plot the dose-response data. This is an important step, as it will help us decided on starting values for the root-finding procedure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1dfa2759",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALTVJREFUeJzt3QmU12W9P/DPsKeyKQqaIK6Iu0ILesWNhK6VWYlb4k6apYmZGF2XbopHux5JOy5dt1Pekrqm2UmScikTl0wFEXAXUlETWbIEg9//PN9zZ/4zLAMMv5nf7/fM63XON+e7zPTwzG9m3r/n+T6fb12pVCoFAAA1r0OlGwAAQHkIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZyDLY/fCHP4yBAwdGt27d4hOf+EQ8/vjjlW4SAECryy7Y3XHHHTFu3Li46KKL4i9/+UvsueeeMXLkyHj77bcr3TQAgFZVVyqVSpGRNEL3sY99LK699tpif8WKFdG/f//4+te/HuPHj6908wAAWk2nyMiyZcviySefjAsuuKDhWIcOHWLEiBExbdq01X7O0qVLi61eCoILFiyIzTbbLOrq6tqk3QBA+1QqlWLJkiWx1VZbFZllQ2UV7P72t7/F8uXLo2/fvk2Op/3Zs2ev9nMmTpwYl1xySRu1EABgVfPmzYutt946NlRWwa4l0uheuiev3qJFi2LAgAFFB/fo0aOibQMA8rZ48eLilrHu3buX5etlFez69OkTHTt2jLfeeqvJ8bTfr1+/1X5O165di21lKdQJdgBAWyjX7V9ZrYrt0qVLDBkyJH7/+983uWcu7Q8bNqyibQMAaG1ZjdglaVr1hBNOiKFDh8bHP/7xuPrqq+P999+Pk046qdJNAwBoVdkFu6OOOireeeeduPDCC2P+/Pmx1157xZQpU1ZZUAEAkJvs6tiV4ybGnj17Foso3GMHANRS7sjqHjsAgPZMsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAcAkAnBDgAgE4IdAEAmBDsAgEwIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAcAkAnBDgAgE4IdAEAmBDsAgEwIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAcAkAnBDgAgE4IdAEAmaibYXXrppbHvvvvGRhttFL169VrtNXPnzo3DDjusuGaLLbaI8847L/71r3+1eVsBACqhU9SIZcuWxZFHHhnDhg2Lm266aZXzy5cvL0Jdv3794pFHHok333wzxowZE507d47LLrusIm0GAGhLdaVSqRQ15NZbb41vfOMbsXDhwibH77333vjMZz4Tb7zxRvTt27c4dv3118f5558f77zzTnTp0mWdvv7ixYujZ8+esWjRoujRo0er/BsAAFojd9TMVOzaTJs2LXbfffeGUJeMHDmy6LCZM2dWtG0AAG2hZqZi12b+/PlNQl1Sv5/OrcnSpUuLrV4KggAAtaiiI3bjx4+Purq6ZrfZs2e3ahsmTpxYDIHWb/3792/V/z8AgCxH7M4999w48cQTm71mu+22W6evlRZNPP74402OvfXWWw3n1uSCCy6IcePGNRmxE+4AgFpU0WC3+eabF1s5pNWyqSTK22+/XZQ6SaZOnVrciLjLLrus8fO6du1abAAAta5m7rFLNeoWLFhQ/DeVNnn66aeL4zvssENssskmceihhxYB7vjjj48rrriiuK/uO9/5Tpx55pmCGwDQLtRMuZM0ZXvbbbetcvyBBx6IAw88sPj4tddeizPOOCMefPDB2HjjjeOEE06Iyy+/PDp1Wvf8qtwJANBWyp07aibYtRXBDgBoK+rYAQCQd4FiAID2TrADAMiEYAcAkAnBDgAgEzVTx649W7FiRVG/b8mSJdG9e/cYMGBAdOggkwMATQl2VW7WrFkxZcqUYjl0vVTnZtSoUTF48OCKtg0AqC6Gfao81E2ePLlJqEvSfjqezgMA1BPsqnj6NY3UNSedT9cBACSCXZVK99StPFK3snQ+XQcAkAh2VSotlCjndQBA/gS7KpVWv5bzOgAgf4JdlUolTdLq1+ak8+k6AIBEsKtSqU5dKmnSnHRePTsAoJ5gV8VSnbp999036urqmhxP++m4OnYAQGOCXRVLdeoeeeSRKJVKTY6n/XRcHTsAoDHBrkqpYwcArC/BrkqpYwcArC/BrkqpYwcArC/BrkqpYwcArC/BrkqpYwcArC/BrkqpYwcArC/BroqlOnWjR49e5QkUaT8dV8cOAGisU5M9qk4Kb4MGDSpWyaYFFeneuzRN64kTAMDKBLsakELcwIEDK90MAKDKmYoFAMiEYAcAkAnBDgAgE4IdAEAmBDsAgEwIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAcAkAnBDgAgE4IdAEAmOlW6AbSeFStWxNy5c2PJkiXRvXv3GDBgQHToIMsDQK4Eu0zNmjUrpkyZEosXL2441qNHjxg1alQMHjy4om0DAFqH4ZtMQ93kyZObhLok7afj6TwAkB/BLsPp1zRS15x0Pl0HAORFsMtMuqdu5ZG6laXz6ToAIC+CXWbSQolyXgcA1A7BLjNp9Ws5rwMAaodgl5lU0iStfm1OOp+uAwDyIthlJtWpSyVNmpPOq2cHAPkR7DKU6tSNHj16lZG7tJ+Oq2MHAHlSoDhTKbwNGjTIkycAoB0R7DKWplsHDhxY6WYAAG3EVCwAQCYEOwCATAh2AACZEOwAADJRE8Hu1VdfjVNOOSW23Xbb+MhHPhLbb799XHTRRbFs2bIm102fPj3233//6NatW/Tv3z+uuOKKirUZAKCt1cSq2NmzZ8eKFSvihhtuiB122CGeffbZOO200+L999+P73//+w0Ptj/00ENjxIgRcf3118eMGTPi5JNPjl69esXYsWMr/U8AAGh1daVSqRQ16Morr4zrrrsuXn755WI/fTxhwoSYP39+dOnSpTg2fvz4uOuuu4pguK5SQOzZs2csWrRorY/mAgDYEOXOHTUxFbs6qQM23XTThv1p06bF8OHDG0JdMnLkyJgzZ0689957FWolAEDbqclg9+KLL8Y111wTX/nKVxqOpZG6vn37Nrmufj+dW5OlS5cWabnxBgBQiyoa7NJUaV1dXbPbytOor7/+evEQ+yOPPLK4z25DTZw4sRgCrd/SogsAgFpU0Xvs3nnnnXj33XebvWa77bZrmF5944034sADD4xPfvKTceuttxaPzKo3ZsyYYrQt3VNX74EHHoiDDz44FixYEL17917jiF3a6qWvkcKde+wAgFq7x66iq2I333zzYlsXaaTuoIMOiiFDhsQtt9zSJNQlw4YNKxZPfPjhh9G5c+fi2NSpU2PQoEFrDHVJ165diw0AoNbVxD12KdSlkboBAwYU5U3SSF+6b67xvXPHHntsMbKX6t3NnDkz7rjjjpg0aVKMGzeuom0HAGgrNVHHLo28pQUTadt6662bnKufSU7DmPfdd1+ceeaZxahenz594sILL1TDDgBoN2q2jl1rUccOAKjV3FETU7EAAKydYAcAkAnBDgAgE4IdAEB7D3YvvfRSfOc734ljjjkm3n777eLYvffeW5QaAQCgRoLdQw89FLvvvns89thjceedd8bf//734vgzzzwTF110UbnbCNSwFStWxKuvvhozZswo/pv2AaiiOnbpGa/f+973iuK/3bt3bzieHt917bXXlrN9QA2bNWtWTJkypVjOXy8t50/Pex48eHBF2waQoxaN2KV33kccccQqx7fYYov429/+Vo52ARmEusmTJzcJdUnaT8fTeQCqINj16tUr3nzzzVWOP/XUU/HRj360HO0Caliabk0jdc1J503LAlRBsDv66KPj/PPPL57VWldXV/xy/tOf/hTf/OY3Y8yYMWVuIlBr5s6du8pI3crS+XQdABUOdpdddlnsvPPO0b9//2LhxC677BLDhw+Pfffdt1gpC7RvS5YsKet1ALTi4okuXbrEj370o7jwwguL++1SuNt7771jxx13bMmXAzLTeFFVOa4DoBWDXb00Ype25cuXFwHvvffei969e2/IlwQyMGDAgGL1a3PTsel8ug6ACk/FfuMb34ibbrqp+DiFugMOOCD22WefIuQ9+OCDZWweUIs6dOhQlDRpTjqfrgOgfFr0W/UXv/hF7LnnnsXH99xzT7z88ssxe/bsOOecc2LChAllbB5Qq1KdutGjRxcjc42l/XRcHTuA8qsrlUql9f2kbt26xYsvvhhbb711jB07NjbaaKO4+uqr45VXXikC39pWw1Wz1PaePXvGokWLVvmDBKy/tGo+rX5NCyXSPXVp+tVIHUDr5I4W3WPXt2/feO6552LLLbcsalFdd911xfF//OMf0bFjxw1uFJCPFOIGDhxY6WYAtAstCnYnnXRSMZWSgl2qYzdixIjieHp2bCqDAgBAjQS7iy++OHbbbbeYN29eHHnkkdG1a9fieBqtS8+RBQCgRu6xy5l77ACAdnWPXfL73/++2N5+++1Vnvd48803b3DDAABog2B3ySWXxHe/+90YOnRow312sDpWRAJAlQe766+/Pm699dY4/vjjy98isjFr1qxi1XTj8jdpmDkVplXDDACqpEDxsmXLYt999y1/a8gq1E2ePHmVmoZpPx1P5wGAKgh2p556avzP//xPmZtCTtOvaaSuOen8yvdmAgAVmIr94IMP4sYbb4zf/e53sccee0Tnzp2bnL/qqqs2sFnUsvSUgbU9fSSdT9cpXAsAFQ5206dPj7322qv4+Nlnn21yzkIK0qOjynkdANCKwe6BBx5oyafRTqTngZbzOgCgFe+xa+yvf/1rsUG99JD3tRVZTOfTdQBAhYNduuk91bFLlZK32WabYuvVq1f853/+pxviKR76nkqaNCedT9cBABWeip0wYULcdNNNcfnll8d+++1XHHv44YeLZ8imhRWXXnppGZtILUp16kaPHq2OHQBU+7Nit9pqq6JI8ec+97kmx+++++746le/Gq+//nrUKs+KLS9PngCAKn9W7IIFC2LnnXde5Xg6ls5BvTTdqqQJALSNFt3ktOeee8a11167yvF0LJ0DAKDttWjE7oorrojDDjusKFA8bNiw4ti0adNi3rx58Zvf/KbcbQQAoLVG7A444IB4/vnn44gjjoiFCxcW2xe+8IWYM2dO7L///i35kgAAVGLxRM4sngAA2tXiieS9994rSp7MmjWr2N9ll13ipJNOik033XSDGwXUBque21cfV1NbgDKO2P3hD3+Iz372s0XCHDp0aHHsySefLKZk77nnnhg+fHjUKiN2sG7Sm7opU6YUPzP10rvNVHw61TEkrz6uprZAThaXecSuRcFu9913LxZNXHfdddGxY8fi2PLly4sado888kjMmDEjapVgB+v2R37y5MlrPJ+KU/tjn08fV1NbIDeLyxzsWjSG/uKLL8a5557bEOqS9PG4ceOKc0C+0nRcGrlpTjqfrqP2+7ia2gK0UrDbZ599Gu6taywdU8cO8pbusWo8Hbc66Xy6jtrv42pqC9BKiyfOOuusOPvss4vRuU9+8pPFsUcffTR++MMfFs+PnT59esO1e+yxR0v+L4AqlW6cL+d1VHcfV1NbgFYKdsccc0zx329961urPVdXVxfp1r3033TvHZCPtBqynNdR3X1cTW0BWinYvfLKKy35NCADqcRFusG3uem5dD5dR+33cTW1BWilYLfNNtu05NOADKS6ZanERXOrJNN59c3y6ONqagv5UBOx9bSo3Mltt90Wffr0KZ4XWz8le+ONNxZFin/605/WdPBT7gTWjbpm7auPq6kt1DavpSqsYzdo0KCiht3BBx8c06ZNi0MOOSSuvvrq+PWvfx2dOnWKO++8M2qVYAfrzrvu9tXH1dQWapOaiFX6SLF58+bFDjvsUHx81113xZe+9KUYO3Zs7LfffnHggQducKOA2pD+qA8cOLDSzchaNfVxNbWF2rOuNRHT4JE3DC3Xordam2yySbz77rvFx/fdd1986lOfKj7u1q1b/POf/9yA5gAAOVITsW20aMQuBblTTz019t5773j++efj3//934vjM2fO9G4OAFiFmohVPGKXChGnZ8W+88478b//+7+x2WabFceffPLJhhp3AAD11ESs4hG7Xr16xbXXXrvK8UsuuaQcbQIAMqMmYtto8XKmP/7xj/HlL3859t1333j99deLYz/+8Y/j4YcfLmf7aOc32r766qsxY8aM4r8eMl5bfP+A1dVEbI6aiBUasUvTr8cff3wcd9xx8Ze//CWWLl1aHE9LdS+77LL4zW9+U4am0Z6pc1TbfP+A1Uk1D0ePHq0mYitqUR27tGjinHPOiTFjxhRz5s8880xst9128dRTT8WnP/3pmD9/ftQqdewqT52j2ub7B6yNmohVVsduzpw5MXz48FWOp4YtXLhwgxtF+6XOUW3z/QPWhZqIVXaPXb9+/eLFF19c5Xi6vy6N3EFLqXNU23z/AGow2J122mlx9tlnx2OPPRZ1dXXxxhtvxO233x7nnntunHHGGeVvJe2GOke1zfcPoLJaNBU7fvz4YsolPSP2H//4RzEt27Vr1zjvvPOKwsXQUuoc1TbfP4AaHLFLo3QTJkyIBQsWxLPPPhuPPvpoUaw43WO37bbblr+VtLs6R81J59N1VB/fP4AaCnaprMkFF1wQQ4cOjf32268oa7LLLrsUjxJLD+2dNGlSsVoWWvyCVOeopvn+VR/1BKF9Wa9yJ+eff37ccMMNMWLEiHjkkUeKUbqTTjqpGLH79re/HUceeWR07NgxaplyJ9VBHbTa5vtXHXwfoP3ljvUKdmnF69VXXx2f+9zniinYPfbYI0488cS46aabiunZHAh21UOdo9rm+1dZ6glCbahoHbu//vWvMWTIkOLj3XbbrVgwkaZecwl1VBd1jmqb71/lqCcI7dd63WO3fPny6NKlS8N+p06dYpNNNmmNdgHQQuoJQvu1XiN2adY2Tb2mkbrkgw8+iNNPPz023njjJtfdeeed5W0lAOtMPUFov9Yr2J1wwglN9r/85S+Xuz0AbCD1BKH9Wq9gd8stt7ReSwAoaz3BdFP2mqgHCXlqUYFiAKqXeoLQftVMsEslVtK70G7dusWWW24Zxx9/fPGM2samT58e+++/f3FN//7944orrqhYewEqafDgwTF69OhVyiek/XQ8nQfy06JnxVbCQQcdVBRBTqHu9ddfj29+85vxpS99qSiUnKQph0MPPbQonnz99dfHjBkz4uSTT45evXrF2LFjK918gDaXwlt6KlBaJZsWVKR779Ib5DSiB+RpvQoUV5Nf/epX8fnPf754zFnnzp3juuuuK55fO3/+/IaSLOPHj4+77rorZs+evc5fV4FiAKCtlDt31OTbtgULFsTtt98e++67bxHqkmnTpsXw4cOb1NkbOXJkzJkzJ9577701fq0UDFOnNt4AAGpRTQW79KzaVDNvs802K6YW7r777oZzaaSub9++Ta6v30/n1mTixIlFUq7f0r15AAC1qKLBLk2VpseRNbc1nkY977zz4qmnnor77rsvOnbsGGPGjCmKJm+ICy64oBj+rN/mzZtXhn8ZAEA7Wzxx7rnnFk+yaM52223X8HGfPn2KbaeddipuCk6ja48++mgMGzYs+vXrF2+99VaTz63fT+fWJD1Fo/5JGgAAtayiwW7zzTcvtpY+5Lr+Hrkkhbu0eOLDDz9suO9u6tSpxYqw3r17l7HVAADVqSbusXvsscfi2muvjaeffjpee+21uP/+++OYY46J7bffvgh0ybHHHlssnDjllFNi5syZcccdd8SkSZNi3LhxlW4+AECbqIlgt9FGG8Wdd94ZhxxySDECl8LbHnvsEQ899FDDNGpa+JDuvXvllVdiyJAhxTTvhRdeqIYdANBu1Gwdu9aijh0A0FbUsQMAoHanYgEAWDvBDgAgE4IdAEAmBDsAgEwIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAcAkAnBDgAgE4IdAEAmBDsAgEwIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCY6VboBALQ/K1asiLlz58aSJUuie/fuMWDAgOjQwVgDbCjBDoA2NWvWrJgyZUosXry44ViPHj1i1KhRMXjwYN8N2ADeHgHQpqFu8uTJTUJdkvbT8XQeaDnBDoA2m35NI3XNSefTdUDLCHYAtIl0T93KI3UrS+fTdUDLCHYAtIm0UKKc1wGrEuwAaBNp9Ws5rwNWJdgB0CZSSZO0+rU56Xy6DmgZwQ6ANpHq1KWSJs1J59Wzg5YT7ABoM6lO3ejRo1cZuUv76bg6drBhFCgGoE2l8DZo0CBPnoBWINgB0ObSdOvAgQP1PJSZqVgAgEwIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkolOlGwDA/7dixYqYO3duLFmyJLp37x4DBgyIDh28BwcyDXZLly6NT3ziE/HMM8/EU089FXvttVfDuenTp8eZZ54ZTzzxRGy++ebx9a9/Pb71rW9VtL0A62rWrFkxZcqUWLx4ccOxHj16xKhRo2Lw4ME6ElirmnsbmILaVltttcrx9Ivw0EMPjW222SaefPLJuPLKK+Piiy+OG2+8sSLtBFjfUDd58uQmoS5J++l4Og+QVbC7995747777ovvf//7q5y7/fbbY9myZXHzzTfHrrvuGkcffXScddZZcdVVV1WkrQDrM/2aRuqak86n6wCyCHZvvfVWnHbaafHjH/84Ntpoo1XOT5s2LYYPHx5dunRpODZy5MiYM2dOvPfee81O7aZ3xI03gLaU7qlb2++edD5dB1Dzwa5UKsWJJ54Yp59+egwdOnS118yfPz/69u3b5Fj9fjq3JhMnToyePXs2bP379y9z6wGalxZKlPM6oP2qaLAbP3581NXVNbvNnj07rrnmmuIX2gUXXFD2NqSvuWjRooZt3rx5Zf//AGhOWv1azuuA9quiq2LPPffcYiSuOdttt13cf//9xVRr165dm5xLo3fHHXdc3HbbbdGvX79iurax+v10bk3S11z56wK0pVTSJK1+bW46Np1P1wFUbbBLJUnStjY/+MEP4nvf+17D/htvvFHcP3fHHXcUpU+SYcOGxYQJE+LDDz+Mzp07F8emTp0agwYNit69e7fivwJgw6Q6damkSVr9uibpvHp2QBb32KV3qbvttlvDttNOOxXHt99++9h6662Lj4899thi4cQpp5wSM2fOLELfpEmTYty4cRVuPcDapTp1o0ePLkbmGkv76bg6dkCWBYrXJC18SKVQUoHiIUOGRJ8+feLCCy+MsWPHVrppAOskhbc0y+DJE0BL1ZXSklMapHtcUkhMCylWfucMAFDNuaMmpmIBAFg7wQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAcAkAnBDgAgE4IdAEAmBDsAgEwIdgAAmRDsAAAyIdgBAGRCsAMAyIRgBwCQCcEOACATgh0AQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJnoVOkGQLmsWLEi5s6dG0uWLInu3bvHgAEDokMH7130IUD7IdiRhVmzZsWUKVNi8eLFDcd69OgRo0aNisGDB1e0bbVCHwLUPsMZZBFIJk+e3CTUJWk/HU/n0YcA7YFgR81Pv6aRuuak8+k69CFA7gQ7alq6p27lkbqVpfPpOvQhQO4EO2paWihRzuvaI30IkA/BjpqWVr+W87r2SB8C5EOwo6alkiZp9Wtz0vl0HfoQIHeCHTUt1alLJU2ak86rZ6cPAdoDwY6al+rUjR49epWRu7Sfjqtjpw8B2ou6UqlUqnQjqklaQdmzZ89YtGjRWqf4qC6ePKEPAdp77vDkCbKRplsHDhxY6WbUNH0IUNtMxQIAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAcAkAnBDgAgE50q3QAAaA9WrFgRc+fOjSVLlkT37t1jwIAB0aGD8RXKS7ADgFY2a9asmDJlSixevLjhWI8ePWLUqFExePBg/U/ZeKsAAK0c6iZPntwk1CVpPx1P56FcBDsAaMXp1zRS15x0Pl0H5SDYAUArSffUrTxSt7J0Pl0H5SDYAUArSQslynkdrI1gBwCtJK1+Led1sDaCHQC0klTSJK1+bU46n66DchDsAKCVpDp1qaRJc9J59ewoF8EOAFpRqlM3evToVUbu0n46ro4d5aRAMQC0shTeBg0a5MkTtDrBDgDaQJpuHThwoL6mVZmKBQDIhGAHAJAJwQ4AIBOCHQBAJmom2KUbTuvq6ppsl19+eZNrpk+fHvvvv39069Yt+vfvH1dccUXF2gsA0NZqalXsd7/73TjttNNW+wiW9BDlQw89NEaMGBHXX399zJgxI04++eTo1atXjB07tkItBgBoOzUV7FKQ69ev32rP3X777bFs2bK4+eabo0uXLrHrrrvG008/HVdddZVgBwC0CzUzFZukqdfNNtss9t5777jyyivjX//6V8O5adOmxfDhw4tQV2/kyJExZ86ceO+99yrUYgCAtlMzI3ZnnXVW7LPPPrHpppvGI488EhdccEG8+eabxYhcMn/+/Nh2222bfE7fvn0bzvXu3Xu1X3fp0qXFVm/RokUNU7sAAK2pPm+USqXyfMFSBZ1//vnpX9HsNmvWrNV+7k033VTq1KlT6YMPPij2P/WpT5XGjh3b5JqZM2cWX+O5555bYxsuuuiitbbBpg+8BrwGvAa8BrwGvAaiFfvgpZdeKku2qiuVLSKuv3feeSfefffdZq/Zbrvtmkyv1ps5c2bstttuMXv27OL5e2PGjClS71133dVwzQMPPBAHH3xwLFiwYJ1H7BYuXBjbbLNN8Ty/nj17btC/r71I/Z5WIc+bN2+Vh1yj37zWKs/PqD7zWqteaaZwwIABxW1jacFnTU/Fbr755sXWEmlhRHru3hZbbFHsDxs2LCZMmBAffvhhdO7cuTg2derUIvStKdQlXbt2LbaVpVAnpKyf1F/6bP3pN33WVrzW9JnXWvVKmaYsXydqQFoYcfXVV8czzzwTL7/8crEC9pxzzokvf/nLDaHt2GOPLUb2TjnllGI074477ohJkybFuHHjKt18AIA2UROLJ9KI2s9+9rO4+OKLi2nTtEgiBbvGoS2NsN13331x5plnxpAhQ6JPnz5x4YUXKnUCALQbNRHs0mrYRx99dK3X7bHHHvHHP/5xg0PkRRddtNrpWfRZOXmt6bO24rWmz7zW2s/PZ0UXTwAAUD41cY8dAABrJ9gBAGRCsAMAyIRg18gPf/jDGDhwYHTr1i0+8YlPxOOPP16570wV+sMf/hCf/exnY6uttoq6uromxaCTdLtmWom85ZZbxkc+8pEYMWJEvPDCC9GeTZw4MT72sY9F9+7di5qLn//854vnFzf2wQcfFKu503OQN9lkk/jiF78Yb731VrRX1113XbEQqr7mWqpRee+99zac11/r/mzt9HP6jW98Q9+tQaq0kPqo8bbzzjvrr3Xw+uuvFyXH0u+t9Pt+9913jz//+c8N5/09WFXKFyu/3tKWfv+X83ebYPd/Ut27VD4lrUz5y1/+EnvuuWeMHDky3n777fXu1Fy9//77Rb+kALw6V1xxRfzgBz+I66+/Ph577LHYeOONiz5ML9b26qGHHip+UNOq7lQwOxXQPvTQQ4u+rJdK99xzzz3x85//vLj+jTfeiC984QvRXm299dZFKHnyySeLPxTp6TGHH354UZ8y0V9r98QTT8QNN9xQBOTG9N2qdt111+K54/Xbww8/rL/WIj0hYb/99iseBpDedD333HPxX//1X00eBuDvwep/Lhu/1tLfhOTII48s789nWR5MloGPf/zjpTPPPLNhf/ny5aWtttqqNHHixIq2q1qll84vf/nLhv0VK1aU+vXrV7ryyisbji1cuLDUtWvX0k9/+tMKtbL6vP3220XfPfTQQw191Llz59LPf/7zhmvS85HTNdOmTatgS6tL7969S//93/+tv9bBkiVLSjvuuGNp6tSppQMOOKB09tlnF8e91lb/rPA999xztf2ov5p/zvu//du/rfG8vwfrJv1sbr/99kV/lfP1ZsQuIpYtW1aMDqSpw8aP9kj76akXrN0rr7wS8+fPb9KHqWh0mtLWh02fCZhsuummxX/T6y6N4jXutzQVlJ4bqN8ili9fXhQnTyOcaUpWf61dGiE+7LDDmrymvNbWLN0ukm4vSc8lP+6444rnhOuv5v3qV7+KoUOHFiNN6RaTvffeO370ox/5e7CeueMnP/lJnHzyycV0bDl/twl2EfG3v/2t+APSt2/fJp2T9lNYYe3q+0kfrtmKFSuK+53SFMZuu+3W0G/pUXgrP/i5vb/2ZsyYUdxjkgp2nn766fHLX/4ydtllF/21FikEp1tJ0r2dK/NaW1V643nrrbfGlClTins70xvU/fffP5YsWaK/mpEe7Zn6a8cdd4zf/va3ccYZZ8RZZ50Vt912W8NrLfH3YM3SPeoLFy6ME088sew/nzXx5AnIZSTl2WefbXIPD6s3aNCgePrpp4sRzl/84hdxwgknFPecsGbz5s2Ls88+u7hvJy0AY+0+/elPN3yc7kdMQW+bbbaJyZMnFwsCWPOb1DRid9lllxX7acQu/W5L91enn1XW7qabbipef2m0uNyM2EUUz5Xt2LHjKqtP0n6/fv3K3uk5qu8nfbh6X/va1+LXv/51PPDAA8XigMb9lobk0zu3xtr7ay+9c91hhx2K5z6n0ae0aGfSpEn6qxlpKict9kqPYOzUqVOxpTCcFjSlj9M7f6+15qXRkp122ilefPFFr7VmpMoHaQS9scGDBzdMY/t70LzXXnstfve738Wpp57aKn8LBLv/+yOS/oD8/ve/b/KOJO2n+3pYu2233bZ48TXuw8WLFxerY9tzH6Z1JinUpanE+++/v+inxtLrLq0sa9xvqRxK+gXZnvttZenncenSpfqrGYccckgxhZ1GOuu3NKqS7hur/9hrrXl///vf46WXXiqCi5/NNUu3k6xctun5558vRjsTfw+ad8sttxT3JqZ7YeuV9fW2XkstMvazn/2sWMF56623lp577rnS2LFjS7169SrNnz+/0k2rqtV2Tz31VLGll85VV11VfPzaa68V5y+//PKiz+6+++7S9OnTS4cffnhp2223Lf3zn/8stVdnnHFGqWfPnqUHH3yw9OabbzZs//jHPxquOf3000sDBgwo3X///aU///nPpWHDhhVbezV+/Phi1fArr7xSvI7Sfl1dXem+++4rzuuvddd4Vay+W9W5555b/Gym19qf/vSn0ogRI0p9+vQpVq/rrzV7/PHHS506dSpdeumlpRdeeKF0++23lzbaaKPST37yk4Zr/D1YvVRxI/2+TyuLV1au322CXSPXXHNN0aldunQpyp88+uij692hOXvggQeKQLfydsIJJxTn05Lt//iP/yj17du3CMmHHHJIac6cOaX2bHX9lbZbbrml4ZoUfL/61a8WJT3SL8cjjjiiCH/t1cknn1zaZpttip/DzTffvHgd1Ye6RH+1PNjpu6aOOuqo0pZbblm81j760Y8W+y+++KL+Wgf33HNPabfddit+1++8886lG2+8scl5fw9W77e//W3xN2B1fxvL9fNZl/5n/cb4AACoRu6xAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAMiHYAQBkQrADAMiEYAfwf0488cSoq6srtvRA7r59+8anPvWpuPnmm2PFihX6Cah6gh1AI6NGjYo333wzXn311bj33nvjoIMOirPPPjs+85nPxL/+9S99BVQ1wQ6gka5du0a/fv3iox/9aOyzzz7x7W9/O+6+++4i5N16663FNXPnzo3DDz88Ntlkk+jRo0eMHj063nrrrYav8cwzzxSBsHv37sX5IUOGxJ///OeG8w8//HDsv//+8ZGPfCT69+8fZ511Vrz//vu+D8AGE+wA1uLggw+OPffcM+68885iSjaFugULFsRDDz0UU6dOjZdffjmOOuqohuuPO+642HrrreOJJ56IJ598MsaPH19M7SYvvfRSMSr4xS9+MaZPnx533HFHEfS+9rWv+T4AG6yuVCqVNvzLAORxj93ChQvjrrvuWuXc0UcfXQSxSZMmxac//el45ZVXitG25Lnnnotdd901Hn/88fjYxz5WjNJdc801ccIJJ6zydU499dTo2LFj3HDDDQ3HUrA74IADilG7bt26tfK/EsiZETuAdZDeA6dFFbNmzSoCXX2oS3bZZZfo1atXcS4ZN25cEeBGjBgRl19+eTFK13iaNk3ppmnc+m3kyJHFSGAKiwAbQrADWAcptG277bbr1FcXX3xxzJw5Mw477LC4//77i+D3y1/+sjj397//Pb7yla/E008/3bClsPfCCy/E9ttv73sBbJBOG/bpAPlL4WzGjBlxzjnnFPfOzZs3r9gaT8WmKdwU4OrttNNOxZY+55hjjolbbrkljjjiiGJBRrp+hx12qOC/CMiVYAfQyNKlS2P+/PmxfPnyYqXrlClTYuLEiUW5kzFjxkSHDh1i9913LxZIXH311UUJlK9+9avFPXJDhw6Nf/7zn3HeeefFl770pWKE769//WuxiCItlkjOP//8+OQnP1kslkjTtRtvvHER9NIijGuvvdb3Atgggh1AIynIbbnlltGpU6fo3bt3sRr2Bz/4QbEQIoW6JJU/+frXvx7Dhw8vjqVVrmmxRJIWRrz77rtFCEzBsE+fPvGFL3whLrnkkuL8HnvsUaymnTBhQlHyJN27l6ZgG6+qBWgpq2IBADJh8QQAQCYEOwCATAh2AACZEOwAADIh2AEAZEKwAwDIhGAHAJAJwQ4AIBOCHQBAJgQ7AIBMCHYAAJkQ7AAAIg//D/I5tnKctcrVAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(d, r, 'o', color='gray')\n",
    "plt.xlabel(\"Dose\")\n",
    "plt.ylim([-50, 0])\n",
    "plt.ylabel(\"Response\")\n",
    "plt.xlim([0, 70])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edc2bc87",
   "metadata": {},
   "source": [
    "Given this dose-response data, it appears that the dose-response begins to asymptote at some point beyond a dose of 10. Therefore, we should select a starting value of $ED_{50}$ below 10. Further, a maximum response value between -30 and -50 seems reasonable, so we will select -40 as the max effect starting value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7bfeef8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    lower = 0\n",
    "    vals = [0., ] + list(theta)\n",
    "    return ee_emax(theta=vals, dose=d, response=r)[1:, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "26a67c7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "estr = MEstimator(psi, init=[-40, 5])\n",
    "estr.estimate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7402fa8e",
   "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>Est</th>\n",
       "      <th>SE</th>\n",
       "      <th>LCL</th>\n",
       "      <th>UCL</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Params</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>E-max</th>\n",
       "      <td>-38.6</td>\n",
       "      <td>4.1</td>\n",
       "      <td>-46.7</td>\n",
       "      <td>-30.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ED50</th>\n",
       "      <td>6.2</td>\n",
       "      <td>2.4</td>\n",
       "      <td>1.5</td>\n",
       "      <td>10.9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Est   SE   LCL   UCL\n",
       "Params                       \n",
       "E-max  -38.6  4.1 -46.7 -30.5\n",
       "ED50     6.2  2.4   1.5  10.9"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setting up a Table for the results\n",
    "results = pd.DataFrame()\n",
    "results['Params'] = [\"E-max\", \"ED50\"]\n",
    "results = results.set_index('Params')\n",
    "results['Est'] = estr.theta\n",
    "results['SE'] = np.diag(estr.variance) ** 0.5\n",
    "ci = estr.confidence_intervals()\n",
    "results['LCL'] = ci[:, 0]\n",
    "results['UCL'] = ci[:, 1]\n",
    "results.round(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4a8f385",
   "metadata": {},
   "source": [
    "These results are close to what is reported in the book (note that a different variance estimator is used). \n",
    "\n",
    "We can also use these parameters to generate a plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "12f99316",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQnNJREFUeJzt3QecFPX9//H3wdE7Uo52NBHp0lRQUBEpGntEAwrYEMWKMaIoYiwQNQqiP9T8VUwkCpagGEVQsYKASG9KR+QoUg6Qfvt/fGayd3u9bZ19PR+PcXdn5s5hdg/e9y2fb4LP5/MJAAAAMa9EpC8AAAAAwUGwAwAA8AiCHQAAgEcQ7AAAADyCYAcAAOARBDsAAACPINgBAAB4BMEOAADAIwh2AAAAHkGwAwAA8AhPBrsXX3xRjRo1UtmyZXXGGWdo/vz5kb4kAACAkPNcsJsyZYqGDx+uRx55RD/++KPatWun3r17a8eOHZG+NAAAgJBK8Pl8PnmItdB17txZL7zwgvM6LS1NDRo00B133KERI0ZE+vIAAABCJlEecvToUS1cuFAPPPBA+r4SJUqoZ8+emjt3bo5fc+TIEWfzsyC4e/dunXTSSUpISAjLdQMAgPjk8/m0f/9+1a1b18ksxeWpYLdr1y6dOHFCtWvXzrTfXq9evTrHrxkzZoweffTRMF0hAABAdlu2bFH9+vVVXJ4KdkVhrXs2Js9v3759Sk5Odm5w5cqVI3ptAADA21JTU50hY5UqVQrK9/NUsKtRo4ZKliyp7du3Z9pvr5OSknL8mjJlyjhbVpUqVSbYAQCAsAjW8C9PzYotXbq0OnbsqM8//zzTmDl73aVLl0J9r5SUEFwgAABACHmqxc5Yt+qgQYPUqVMnnX766Ro3bpwOHjyo66+/vlDfZ+VKqXnzkF0mAABA0Hku2F199dXauXOnRo0apZSUFJ122mmaMWNGtgkVBQl2l18esssEAAAIOs/VsQvGIMYqVaqof/99mjyZyRMAACD0ucMmbwZj0qanxtgF06pVkb4CAACAwiHY5RHs0tIKeTcBAAAiiGCXi8OHpQ0bwvtmAAAAFAfBLg/LlhXr3gIAAIQVwS4Py5eH740AAAAoLoJdHgh2AAAglhDs8kCwAwAAsYRgl4c1a6SjR8P3ZgAAABQHwS4XlSpJx4+74Q4AACAWEOxy0aKF+0h3LAAAiBUEu1y0auU+EuwAAECsINjlghY7AAAQawh2uWjZ0n2kxQ4AAMQKgl0+wW79eunAgTC+IwAAAEVEsMvFSSdJSUnu85Uri3p7AQAAwodgl4fWrd1HumMBAEAsINjlgWAHAABiCcEuDwQ7AAAQSwh2eWjTxn1ctixM7wYAAEAxEOwKMDM2JUXatas4txkAACD0CHZ5qFhRatzYfb5iRRjeDQAAgGIg2OWDcXYAACBWEOwKOM6OkicAACDaEewK2GLHBAoAABDtCHaF6Ir1+cLwjgAAABQRwS4fzZtLpUpJ+/ZJGzcW9TYDAACEHsEuH6VLS6ed5j6fPz8M7wgAAEAREewK4Iwz3EeCHQAAiGYEuwI4/XT3cd68EL8bAAAAxUCwK0Sw+/FH6dix4txuAACA0CHYFUCzZlKVKtKhQ6xAAQAAohfBriA3qQTdsQAAIPoR7ArZHcsECgAAEK0IdgVEsAMAANGOYFfIYLdihbR/fwjfEQAAgCIi2BVQUpKUnOwuK7ZwYVFvNwAAQOgQ7AqB7lgAABDNCHaFQKFiAAAQzQh2hcDSYgAAIJoR7AqhQwe3pt0vv0i//hq6NwUAAKAoCHaFULGi1KqV+5x6dgAAINoQ7AqJ7lgAABCtCHaFxMxYAAAQrQh2RQx2CxZIaWkheEcAAACKiGBXSDbGrnx5KTVVWrOmqLcdAAAg+Ah2hZSYKHXs6D6fNy8E7wgAAEAREeyKgHF2AAAgGhHsioBgBwAAohHBrhglT5YskQ4dCvI7AgAAUEQEuyJITpZq1ZKOH5cWLy7qrQcAAAgugl0RJCRkdMcygQIAAEQLgl0RdeniPn7zTRDfDQAAgGIg2BXReee5j19+SaFiAAAQHQh2RdSpk1ShgrR7t7R0aXDfFAAAgKIg2BVRqVJS9+7u89mzi/pdAAAAgodgF4Tu2C++CNK7AQAAUAwEu2Lo0cN9/Pprt/QJAABAJBHsiuG006SqVaXUVGnRouC9KQAAAEVBsCuGkiWlc85xn9MdCwAAIo1gF6RxdkygAAAAkUawC1Kws0LFR48G4R0BAAAoIoJdMbVuLdWoIf3+u7RgQXG/GwAAQNER7IqpRAnp3HPd54yzAwAAkUSwC2LZE8bZAQCASCLYBXGc3Zw50uHDwfiOAAAAhUewC4LmzaU6daQjR6S5c4PxHQEAAAqPYBcECQmUPQEAAJFHsAsS1o0FAACRRrAL8gSK+fOlgweD9V0BAAA8GOyeeOIJde3aVeXLl1dVW6A1B5s3b9ZFF13knFOrVi3dd999On78eFiur3FjKTlZOnZM+u67sPwvAQAAYjPYHT16VFdddZVuvfXWHI+fOHHCCXV23pw5c/TGG29o0qRJGjVqVNjG2flb7ahnBwAAIiFmgt2jjz6qe+65R23atMnx+MyZM7Vy5Uq9+eabOu2009S3b1899thjevHFF52wFw6sGwsAACIpZoJdfubOneuEvtq1a6fv6927t1JTU7VixYqwBrsffpB27w7L/xIAAMB7wS4lJSVTqDP+13YsN0eOHHHCX+BWVA0aSNagmJYmffJJkb8NAABA7AW7ESNGKCEhIc9t9erVIb2GMWPGqEqVKulbA0tnxXDxxe7jhx8G5/oAAAAKKlERdO+992rw4MF5ntOkSZMCfa+kpCTNt1ojAbZv355+LDcPPPCAhg8fnv7aWuyKE+4uuUR68km3xc6G9pUuXeRvBQAAEDvBrmbNms4WDF26dHFKouzYscMpdWJmzZqlypUrq2XLlrl+XZkyZZwtWDp3ti5gC5XS119LPXsG7VsDAAB4Y4yd1ahbvHix82ilTey5bQcOHHCO9+rVywlw1113nZYsWaJPP/1UDz30kIYNGxbU4JafEiWkP/zBfU53LAAACKcEn8/nUwywLlurTZfV7Nmzde655zrPN23a5NS5+/LLL1WhQgUNGjRIY8eOVWJiwRsmrSvWxtrt27fPae0rCgt0l14qNWokrV/v1rgDAAAIRe6IyWAXSzf499+lk06SDh+Wli51Z8oCAACEOtjFTFdsLClfPmNsHd2xAAAgXAh2IeIvezJ9eqj+DwAAAJkR7ELEP4Fi3jwrkByq/wsAAEAGgl2I1K3rlj4x//1vqP4vAAAAGQh2IcQqFAAAIJwIdiFkq1CYWbOkQ4eK/n3S0tK0ceNGLVu2zHm01wAAAFG18oTXtW0rJSdbcWXp888zxt0VxqpVqzRjxgxnOrSfTYfu06ePWrRoEdwLBgAAMY0WuxCywsTF6Y61UDd16tRMoc7Ya9tvxwEAAPwIdmHqjv3oI+tSLfjXWXertdTlxY7TLQsAAPwIdiF2zjlSxYrStm3SwoUF/zpbEzdrS11WdtzOAwAAMAS7ECtTRurTx33+/vsF/7r9+/cH9TwAAOB9BLswuOoq9/GttwreHVupUqWgngcAALyPYBcGNoHC8temTdLcuQX7muTk5HwXA7bjdh4AAIAh2IVBuXLSFVe4zydPLtjXlChRwilpkhc7bucBAAA4+YHbEB79+7uPU6dKx44V7GusTl3Xrl2VYHVTAthr208dOwAAEIhgFyY9eki1a0u//SZ9+mnBvsbq1M2ZM0c+ny/Tfntt+6ljBwAAAhHswiQxUbrmGvf5v/+d//nUsQMAAIVFsAujAQPcxw8+kA4cyPtc6tgBAIDCItiFUadOUrNm0u+/S9Om5X0udewAAEBhEezCyOZA+CdR5Dc7ljp2AACgsAh2EeqOnTVL2rEj9/OoYwcAAAqLYBdm1hXbubN04oRb+iQ31LEDAACFRbCLYKtdft2xVqeuX79+2VagsNe2nzp2AAAgUIIva5G0OJeamqoqVapo3759+S7pVVQpKVK9eu66sWvXSk2b5l/6xGbJ2oQKG3tn3bSsOAEAQOxLDXLuoMUuApKSpPPPL3hNOwtxjRo1Ups2bZxHQh0AAMgJwS4KumNpMwUAAMFAsIuQK66QKlSQ1qyRvvkmUlcBAAC8hGAXIZUqZdS0mzgxUlcBAAC8hGAXQbfe6j6+9560fXskrwQAAHgBwS6C2reXzjhDOnZMeu21SF4JAADwAoJdlLTavfyyW7QYAACgqAh2Edavn1StmrRpkzRjRqSvBgAAxDKCXYSVKyddf737nEkUAACgOAh2UWDoUPfx44+lDRsifTUAACBWEeyiQLNmUs+ebqHiV16J9NUAAIBYRbCLskkUr74qHTkS6asBAACxiGAXJS65RKpbV9q5U3r//UhfDQAAiEUEuyiRmCjdfLP7nEkUAACgKAh2UcSCXcmS7tqxy5dH+moAAECsIdhFkXr13C5ZM25cpK8GAADEGoJdlLn3Xvfxn/+Utm6N9NUAAIBYQrCLMmedJXXr5q4f+9xzkb4aAAAQSwh2UWjEiIz1Y3fvjvTVAACAWEGwi0J9+0pt2kgHDkj/93+RvhoAABArCHZRKCEho9Vu/Hjp998jfUUAACAWEOyiVL9+UuPG0q5d7moUAAAA+SHYRXHB4vvuc58/84w7mQIAACAvBLsoNniwVKuWtHmz9Pbbkb4aAAAQ7Qh2UaxcOenuu93nY8dKaWmRviIAABDNCHZR7rbbpMqVpZUrpY8+ivTVAACAaEawi3JVqki33uo+f/JJyeeL9BUBAIBoRbCLAdYdW7asNG8erXYAACB3BLsYkJQk3XWX+/zBB6UTJyJ9RQAAIBoR7GLE/fdLVatKy5dLkycX7GvS0tK0ceNGLVu2zHm01wAAwLsSI30BKJhq1dzVKGwbNUq6+mqpTJncz1+1apVmzJih1NTU9H2VK1dWnz591KJFC247AAAeRItdDLnjDqluXWnTJumll/IOdVOnTs0U6oy9tv12HAAAeA/BLoaULy898oj7/PHHLahlP8e6W62lLi92nG5ZAAC8h2AXY264QTrlFHcN2WefzX588+bN2VrqsrLjdh4AAPAWgl0MriH7xBPu87//XdqxI/Px/fv3F+j7FPQ8AAAQOwh2MejKK6VOnaQDBzJCnl+lSpUK9D0Keh4AAIgdBLsYlJDgrh1rJk6UNmzIOJacnOzMfs2LHbfzAACAtxDsYtT550s9e0rHjkl//nPG/hIlSjglTfJix+08AADgLfzrHsOee04qWVJ6/32b6Zqx3+rU9evXL1vLnb22/dSxAwDAmxJ8PpaVzzpjtEqVKtq3b1++XZrR4N573dmxJ5/srkoRWLTYSprY7FebKGFj6qz7lZY6AAC8mztosYtxVteuTh1p7VrpmWcyH7MQ16hRI7Vp08Z5JNQBAOBtBLsYZ+Heyp4YmyG7cWOkrwgAAEQKwc4DrrlGOvdc6dAh6Z57In01AAAgUgh2Hil/8sILbvHiadOkjz+O9BUBAIBIINh5RKtW0t13u8/vvFM6fDjSVwQAAMItJoLdxo0bdeONN6px48YqV66cmjZtqkceeURHjx7NdN7SpUvVrVs3lS1bVg0aNNBTTz2leDJqlFS3rrRunfS3v0X6agAAQLjFRLBbvXq1U7rj5Zdf1ooVK/Tcc8/ppZde0oMPPphpunCvXr3UsGFDLVy4UE8//bRGjx6tV155RfHCVgmz0ifmySellSsjfUUAACCcYraOnQW3iRMnav369c5rez5y5EilpKSodOnSzr4RI0Zo2rRpTjD0ah27rOzdvPhi6b//ddeTnTvXHXsHAACiD3Xs/seCV/Xq1dNvzNy5c9W9e/f0UGd69+6tNWvWaM+ePYqniRTWSFm1qvTDD3TJAgAQT2KiKzartWvXasKECbrlllvS91lLXe3atTOd539tx3Jz5MgRJy0HbrHOxtlNmOA+f/RRG3sY6SsCAACeD3bWVZqQkJDnlrUbdevWrc4i9ldddZVuvvnmYl/DmDFjnK5X/2aTLrxgwADpssukY8ekQYOkLPNMAACAB0V0jN3OnTv122+/5XlOkyZN0rtXf/31V5177rk688wzNWnSpExLZA0cONBpbbMxdX6zZ89Wjx49tHv3blWrVi3XFjvb/Ox7WLiL1TF2gbZvd8ug2C22pcdGj470FQEAgFCOsYvosPqaNWs6W0FYS915552njh076vXXX8+27mmXLl2cyRPHjh1TqVKlnH2zZs1S8+bNcw11pkyZMs7mRdYT/eKL7soUttzYJZdIHTpE+qoAAEBcj7GzUGctdcnJyXrmmWeclj4bNxc4dq5///5Oy57Vu7OSKFOmTNH48eM1fPhwxbN+/aQ//lE6flwaPNhaKCN9RQAAIFRiohCGtbzZhAnb6tevn+mYvyfZmjFnzpypYcOGOa16NWrU0KhRozRkyBDFM5sl+3//J331lbRsmY1rlJ57LtJXBQAAQiFm69iFSqzXscvN9OluV6yxYYiXXhrpKwIAANSxQ5FY0WJ/r7R1yW7cyI0EAMBrYmKMHYJjzBjp9NOlvXvdCRWUQAEAwFsIdnHEqsZMmeKuSjFvnvTAA5G+IgAAEEwEuzjTqJH0+uvu82eflT78MNJXBAAAIh7s1q1bp4ceekh/+tOftGPHDmffJ5984pQaQXSzFSnuvjtjvN2mTZG+IgAAELFg99VXX6lNmzaaN2+e3n//fR04cMDZv2TJEj1iSxwg6v3tb1LnztKePdJVV0mHDkX6iuBVaWlp2rhxo5YtW+Y82msAQBTVsbM1Xh9//HGn+G+lSpXS99vyXS+88EIwrw8hHm/XsaO0YIFky+7+619u3TsgWFatWqUZM2Y40/n9rIyQrffcokULbjQAREOLnf3mffnll2fbX6tWLe3atSsY14UwaNxYevddqWRJafJktxUPCGaomzp1aqZQZ+y17bfjAIAoCHZVq1bVtm3bsu1ftGiR6tWrF4zrQpj06CFNmOA+f/BB6YMPuPUoPututZa6vNhxumUBIAqC3TXXXKP777/fWas1ISHB+cv5u+++05///GcNHDgwyJeIULv1Vum222x5NmnAAGnpUu45imfz5s3ZWuqysuN2HgAgwsHuySef1KmnnqoGDRo4Eydatmyp7t27q2vXrs5MWcSeceOk88+XDh50lx7buTPSV4RYtn///qCeBwAI4eSJ0qVL6x//+IdGjRrljLezcNe+fXs1a9asKN8OUaBUKWnqVOmMM6S1a6UrrpA++0wqUybSV4ZYFDipKhjnAQDCUKDYWuwuvPBCXXnllTp48KD2WO0MxKzq1aXp06UqVaRvv5UGDbKxUpG+KsSi5ORkZ/ZrXuy4nQcAiHCwu/vuu/Xqq686z0+cOKFzzjlHHTp0cILel19+GcTLQ7ideqo7U9Za8Kwcyl13uWPvgMIoUaKEU9IkL3bczgMABE+R/lZ999131a5dO+f59OnTtX79eq1evVr33HOPRo4cGcTLQyT07JlR087KEj7+OO8DCs/q1PXr1y9by529tv3UsQOA4Evw+QrfHlO2bFmtXbtW9evX15AhQ1S+fHmNGzdOGzZscAJffrPhoplde5UqVbRv3758u5K8zkLdHXe4z196SbrllkhfEWKRzZq32a82UcLG1Fn3Ky11ABCa3FGkFrvatWtr5cqVTjes1aK64IILnP2///67Slq1W3jC7bdLDz+cURLFumiBwrIQ16hRI2cZQnsk1AFA6BQp2F1//fVOV0rr1q2dOnY9re9OctaOtTIo8I5HH3Vb6vw17r74ItJXBAAAglruZPTo0U6o27Jli6666iqV+V9NDGuts3Vk4R02zu7FFyVbKe6996SLL5Y++UTq3j3SVwYAAIIyxs7LGGOXs8OHpcsukz79VCpfnnAHAEA05o4itdiZzz//3Nl27NiRbb3H1157rdgXhuhStqw0bVpGuOvbl3AHAIAnxtg9+uij6tWrlxPsdu3a5RQmDtzg7XDXu7dNlHHD3ddf5/01Fvo3btzorFBijyz6DgBAlHXF1qlTR0899ZSuu+46eQ1dscHrll21apUzazqw/I01M1thWmqYAQCg6Ch3cvToUXXt2pX3I07l1HI3a1b2UDd16tRsNQ3tte234wAAILiKFOxuuukm/fvf/w7ypSCWw91FF0lTp7rHrLvVWuryYsfplgUAILiKNHni8OHDeuWVV/TZZ5+pbdu2KmULiwZ49tlng3V9iPJw98EHkvXIv/OOdM01blmUCy/cnO/qI3bcViOwgrUAACCCwW7p0qU67bTTnOfLly/PdMwKFiN+WAnDt96SatSQJk6Uhg2zVSrKqVYttwZeXmyJKQAAEOFgN3v27CBeAmKdrSJnRYwtzNlKFRMn1lbnzheqb99PVKJE7nNzbN1QAAAQ4TF2gX755RdnQ3yz1rnRo6UJE+y5TwsWdNY77/xRR4/m/LuDzfyxxeABAECEg50Nev/rX//qTM9t2LChs1WtWlWPPfYYA+Lj3O23W9dsghITfVq1qqUmTRqs/fsrZjvPSp6wGDwAAFEQ7EaOHKkXXnhBY8eO1aJFi5ztySef1IQJE/Twww8H+RIRa66+2lYmSVDVqsf166/19I9/3KRt25LSW+r69etHHTsAAKKlQHHdunX10ksv6ZJLLsm0/4MPPtBtt92mrVu3KlZRoDh41q2T/vAHn1avTlDZsmkaP36XbrqpBi11AABEU4Hi3bt369RTT8223/bZMcA0bSrNnZugCy6wEjklNHRoLT39dAkV/lcJAAAQsmDXrl07pys2K9tnxwC/qlWljz+WbrtNTqAbMcKtd0elEwAAoqQr9quvvtJFF13kzGrs0qWLs2/u3LnasmWLPv74Y3Xr1k2xiq7Y0LGSKHffLR0/LrVoIb33nvsIAEC8So2GrthzzjlHP/30ky6//HLt3bvX2a644gqtWbMmpkMdQsuKF3/1lY3RtLVkpdNPz1iGDAAARKjFzstosQu97dvd7tgvv3RfWyveU09JWVamAwDA81KD3GJX5GC3Z88evfrqq1plTS+SWrZsqeuvv17Vq1dXLCPYhYd1xz70kPS3v7mvzzxT+ve/pcaNw3QBCAqraWlr/trycLaSiA3PoD6hd+9xNF0L4BWp0RDsvv76a1188cXOhXTq1MnZt3DhQqdLdvr06erevbtiFcEuvKZNkwYNsvtuNe7c9Wb79w/zRaBI7Je6GTNmOD8zfvaXkhWfbsHgSc/d42i6FsBLUqMh2LVp08aZNDFx4kSVtIVCJZ04ccKpYTdnzhwtW7ZMsYpgF34bN0oDBkhz5rivr7vOZli7QQ/Ryf6Rn5rHAEmKUHvrHkfTtQBekxoNkyfWrl2re++9Nz3UGXs+fPhw5xhQGI0auZMqbK1Z69X517+k9u2lefO4j9HIuuOs5SYvdtzOQ+zf42i6FgAhCnYdOnRIH1sXyPZRxw5FkZgoPfKIdfNLDRtK69dLZ53ljsM7coR7Gk1sjFVgd1xO7Lidh9i/x9F0LQDyl6giuPPOO3XXXXc5rXNn2qh3Sd9//71efPFFZ/3YpUuXpp/btm3bovwvEKcszC1e7BY0fust6YknbKk6adIkqWPHSF8djA2cD+Z5iO57HE3XAiBEwe5Pf/qT8/iXv/wlx2MJCQmyoXv2aGPvgMKuVmEzZK+8Urr1Vmn5cumMM9xVKx5+WCpThvsZSTYbMpjnIbrvcTRdC4AQBbsNGzYU5cuAQrFgd8450u23S1OmZLTevfaa1LkzNzNSrMSFDfDNq3vOjtt5iP17HE3XAiBEY+waNmxY4A0ojho1pLfflt59V6pZM6P17o47pH37uLeRYHXLrMRFXuw49c28cY+j6VrgHTbZZuPGjU4VDXtk8k3wFKncyRtvvKEaNWo468X6u2RfeeUVp0jxW2+9FdOBjnIn0WvXLumee6Q333Rf16kjjR8v/fGPUkJCpK8u/lDXLL7ucTRdC2Ibn6UorGPXvHlzp4Zdjx49NHfuXJ1//vkaN26cPvroIyUmJur9999XrCLYRb/PP3fH3v38s/u6b1+37l2TJpG+svjDSgTxdY+j6VoQm6iJGKXBrnz58lq9erXzQ33//fdr27Zt+uc//6kVK1bo3HPP1c6dOxWrCHax4fBhaexYacwY6ehRqWxZ6b77pPvvlypUiPTVAQBy+sVg/Pjx+Y7XtKob8fQLQ2o0FCiuWLGifvvtN+f5zJkzdcEFFzjPy5Ytq0OHDhX7ooD8WJCzgsa2yMn557tB77HHpFNPdWfUFm0FZABAqFATMTyKFOwsyN10003O9tNPP+nCCy909luLXSNbRgAIk1NOkWbNkt57z13B4pdf3OXJzj5b+uEH3gYAiBbURIziYGeFiG2tWOtyfe+993TSSSc5+xcuXJhe4w4IF5s4ccUV0sqV0uOP21ABd93Z00+XBg6UNm3ivQCASKMmYngUaYydlzHGLvZt3eoWM/bPni1d2i2P8sAD0v9+BwEAhBlj7KJ4jJ355ptvdO2116pr167aav+SyhZv/5e+/fbbYl8UUJw6R/Xq2WdRmj9fOu88d3LF3/8uNW0q/e1vEsNAw4M6VQACURMxilvsrPv1uuuu04ABA5wwt3LlSjVp0kQvvPCCPv74Y2eLVbTYeavOkX26P/3UnS3rX8K4bl1p5EjpxhtZnixUqFMFgL8fYqjcSfv27XXPPfdo4MCBTp/5kiVLnGC3aNEi9e3bVykpKYpVBDtv1jmyJYtttuxDD9nMLHdfgwZuwLv+ere7FsFBnSoA+aEmYpR1xa5Zs0bdu3fPtt8ubO/evcW+KMT3D7u11OXFjhd2+ZmSJaXrrpN++sktZmytdlu2SEOHujNrX31VOnasmBePkL1/ALzXLWtVNNq0aeM8xlPdulAr0p1MSkrS2rVrs+238XXWcgdEa52jMmWkYcOkdevc5ciSktxZszfdJDVrZjO+GYNXHNSpAoAYDHY333yzUxl63rx5SkhI0K+//qrJkyfr3nvv1a221hMQ5XWOrMDxnXe6Ac8mVtSq5Qa822936+HZJIt88iUi+P4BAIIY7EaMGKH+/fs7a8QeOHDA6Za1YsUW6uwRiJU6R1bzbvhwaeNGt4u2YUNpxw63XEpysjsGL4aHjIYddaoAIAaDnbXSjRw5Urt379by5cv1/fffO8WKbYxd48aNg3+ViBu2/nB+g0ftuJ0XTOXKuV20P/8svfGGuzTZvn3Sk0+6Ye+GG6Tly4P6v/SkSL1/AIAiBLsjR47ogQceUKdOnXTWWWc5ZU1atmzpLCXWvHlzZ3Ffmy0LxGqdo1Kl3NUqVqyQ3n9f6trVrYP3+utSmzZS7962PjJr0Ubr+4fsqCcIxJdClTu5//779fLLL6tnz56aM2eO00p3/fXXOy12Dz74oK666iqVtOmHMYxyJ9Ehmuqgff+9Ow7Pgp5/Mqddgo3Hs5m2QeoV9pRoev/iGe8DEP0iWsfOZryOGzdOl1xyidMF27ZtWw0ePFivvvqq0z3rBQS76BFtdY42bHBn0lpplAMH3H0W6gYPlm67ze2+RfS+f/GGeoJAbIhosCtdurQ2bNigerZmkzMuqZzmz5/v1KHxCoId8v+MSP/8pzvZYs2ajP3nny/dcot06aUUPEZksSYnEDsiWqD4xIkTTrjzS0xMVMWKFYt9EUAssZ8764ZdtUqaNcsNctYQ9fnntiqGVL++u4RZDqUegbCgniAQvxILc7I17lnXaxmr8irp8OHDGjp0qCpUqJDpvPdtMBLgcTb6oGdPd7NyKdZFa9u2bdJTT7lbjx7ujNrLL3dLqwDhQD1BIH4VqsVu0KBBqlWrltNkaNu1116runXrpr/2b0C8saLGjz3mrkM7bZp04YVu8PviC+naa6U6daQhQ6S5c5lRi9CjniAQvwo1xi4eMMYOwWIhz8qkTJrktuj5NW9uvyRJ/fu7NfKAYGOMHRA7IjrGDkDBWQ3eRx5xly2bPdsNc9YdaxMuHnzQbeXr3l16+WVp927uLIKHeoJA/IqZYGclVqxcQtmyZVWnTh1dd911zhq1gZYuXapu3bo55zRo0EBP2SAnIMJsYsW557otd7Y82WuvuWPvrKv2m2+koUOlpCR3Esbbb2eUUgGKw+oF9uvXL1sLgL22/dQTBLwpZrpin3vuOXXp0sUJdVu3btWf//xnZ78VSvY3ZZ5yyilO8WRbHWPZsmW64YYbnLp7Q2xwUwHRFYtw2bpVeustafJkafHizMubXXSRO8PWHpl0geKgniAQ3SJaxy6afPjhh7rsssucZc5KlSqliRMnOuvXpqSkpJdkGTFihKZNm6bVq1cX+PsS7BAJtoSZhbwpUzKXSbFQ94c/SFdeKfXtyyoXAOA1qYyxs/FIuzV58mR17drVCXVm7ty56t69e6Y6e71799aaNWu0Z8+eXG+oBUO7qYEbEG6tWkmPPy799JP0449uHbzGjaXff5emTpWuvlqqWdPtrn3jDcbkAQBifIydf61aq5l30kknOQU4P/jgg/Rj1lJXu3btTOf7X9ux3IwZMyZTqRYbmwdEio27a99eGjvWnXQxf74b8k4+2X4JsZZqdwkz+2jbShe2xNn69bxfAIAoCHbWVWprzOa1BXaj3nfffVq0aJFmzpypkiVLauDAgU7R5OKw8XjWr+3ftmzZEoQ/GRCckNe5sxvyrCVv2TJp9GipbVvp+HG3Rt7dd0tNm0q2qt/IkdL339sKMdx9AIhXER1jt3PnTv322295ntOkSZNM3at+v/zyi9O6ZpMnbFKFhTzrRrUxdX6zZ89Wjx49nK7batWqFeiaGGOHWGCtedOnuy14X3+dOczVqCH16eMWSe7dW6pePZJXCgAIZ+4o1JJiwVazZk1nK+pML/8YOWPhziZPHDt2LH3c3axZs9S8efMChzogVlgrnbXW2WZDSD/5RLKRCZ9+Ku3aJb35prtZqZUuXdygZyGvY0d3HwDAm2JiVuy8efO0YMECnX322U5IW7dunR5++GFt375dK1ascNautaRrIa5Xr17OWLzly5c75U6sTArlThAvjh1zly37+GPpv/+Vli/PfPykk6QLLnBDXq9eUt26kbpSAEDcljuxmnR33XWXlixZooMHDzq17Pr06aOHHnpI9erVy1SgeNiwYU4IrFGjhu644w4n5BUGXbHw2rJm1ppnLXmff26f78zHW7SQevZ0J2JYEWWWegaA8IrLYBdOBDt4uTVv3jw35Nn2ww9S4E+/ddHaZI3zznO3s86SKlSI5BUDgPelEuxi6wYD0crWp/3yS+mzz9zt558zH7ehqqef7rbk2WZj9Qh6ABBcBLsQI9ghnrttrbvWwt7s2VLWyj+Jie7ki+7d3c1a9JiXBADFQ7ALMYId4HbRbtjgBjzbrKRK1qBndfZsxYyzz3ZDnm2NGrn7AQAFQ7ALMYIdkLNNm9yA59+saHJWdepIXbu63ba2degglS3LHQWA3BDsQoxgBxTM9u3SnDnSt99K333nrnFrEzSyjtOzcHfmmdIZZ7hj9po0oVUPAPwIdiFGsAOK5tAhacECt46ef9uxI/t5VkvPAp5tNgvXtlq1uOsA4lMqs2Jj6wYD8T5OzwKerWE7f760eLF09Gj2cxs0kDp1ytislc+WRgMAr0sl2MXWDQaQwVYAXLLEDXlWU89q6a1Zk7menl9yshvwbCauPbZvLyUl0Y0LwFtSCXaxdYMB5PczJy1a5IY8/7Z2bc7nWpftaae5Ic8ebWvWTCpZkrsMIDYR7GLsBgMovH373G5bm5Bh28KFbsteWlr2c23WbevWUtu27taundSmjTuWDwCiHcEuxm4wgOD4/Xdp+XK3dc9Cn21Ll7r7c2KlVyzgWeizR9tsbdzy5XlHAEQPgl2M3WAAoXPihLR+vRvwbOye/3HjxpzPt+LJjRu7hZX9W8uW0qmnEvgARAbBLsZuMIDw279fWrFCWrYsY7PWvl27cg98DRu6Ic9a9fybBb7q1cN99QDiSSqTJ2LrBgOIHlZXzwJf4LZypfTbb7l/Tc2absBr3jzzZi1/VoAZAIqDYBdiBDsg/uzcKa1a5YY8/6NN1si6Pm6gxER3FQ0LeTYz95RT3Efb6tWTSpQI558AQKxKpcUutm4wgNh14IC7Jq6FPAt8/uf2mNukDVOunNS0qXTyydm3+vUpzwIgA8EuxAh2APJjZVe2bs0IeT//7G723FbbOH4896+17ttGjdzg59+s5c82696tWJH7D8STVFrsYusGA4gvx465s3LXrXMLLdvmD362347nxYowW8CzoGcB0J77N1uNg3F9gLekEuxi6wYDQGB5ll9+cUNf4GatfFa2Zc+evO+Vjduz8Xs2g9dCn2323L9Z8LOCzQBiB8Euxm4wABTU3r0ZIc8eAzdr7Tt8OP/vYS1+gUHPtgYNMp7bLF8r7wIgOhDsYuwGA0Aw+HzS9u1uwNu0yX30b/batrwmdPiVKeNO4LCwZ1vW57bZcmyEPyA8CHYxdoMBIFzBb/duN+Bt3uw+WrkWe+7fUlLc8woa/myzrt+cNluyjfF+QPTljsQgXBMAIMKshc1a2mzr0CHnc44edWfzWuDLutl+G/9nRZyPHMkY/5fX/8+6dS3k1a2bsVngC3y0rmGr+QcgPPhxA4A4Ubp0xgzb3Fio+/VXN+T5w549Bm523Gb3Wgi0bdGivAOghTsLeoFbUlLG5n9doUJI/thAXCHYAQAydcPmF/6sjp8tw2YBzx/2tm1zX/sfbbOuXzvXxgbatnhx3jfaavjVrp059NnrnLby5XnTgJwQ7AAAhWJlV6wb1rZ27fIu72LLtVnYC9ws8AVuts8mfthKH7bl1QXsZ6171hJoIc8ec9vsGmvUoDsY8SPB5yvIUNr4weQJAAgv+1do//6Mlr3A0OffF7gVpOxLVtWrZ4TRwM1CX9bn9mjLwgHhwOQJAICn2Dg8mwxoW7NmBQuB/vF9FvT8j9Y66N/vf75rV8aMYdtsGbiCsK7ewKBnk1LsMfC5f7KKf6N7GNGArlgAQEyGwJNPzv986w62QOcPe4GBz577H/2bjR20iSHWNewvE1NQtupH1rBnm7UW+h/9m/91tWruuEYgWAh2AADPKlkyo6u1IPwtgoGhz8KePfc/+jd77d+OH3e7iP2TSQrDxgtawPOHPnvuf53Tc/9Wtar75wMCEewAAMihRbBp08KFwcCgZ5u1FAY+D9xsny0hZ7OGDx50NystU1h2nRbwAsOe/9G/ZX1tW5UqUqVKrDDiRQQ7AACCFAbzKhOTlYW61NTsoW/PnozHrM/9m80eNvb1thWmyzhwdrMFPNv8Yc//mNNz+/P5X/s3m2TC8nPRhWAHAEAEWLDyt6A1aVK4r7VxgBbwrNXPNv/znPb5n+/bl/Hcvt6CpT8oFpWtKuIPfAV59G/WWhj42sYnEhCDg2AHAECMsXV6/bX6Csu6jm08YGDY8z/6n+e3WSuhBUMbW+hvaSwOGyuYNfTl9FiQrUIFNzTHK4IdAABxxFrGrAvVNlvPtygsHFp3sD/kBT76nwduOe2zzcYm+mcvF7f1MOsqJhUrZoQ9//OcHrM+z7pZUIylLmeCHQAAKBQLOf7QVBzW6mcB0QJe1sDn35fT89y2tDT3+/pXMbEi18H682YNe1mf5/SY13PbQlH7kGAHAAi7tLQ0bd68Wfv371elSpWUnJysEvHcfxan7C33d8HWq1e872WtiIcOZYQ8f2AMfO4PfFmP26xk/zH/fttn9Qz939t/brAFu44hwQ4AEFarVq3SjBkznKWU/CpXrqw+ffqoRYsWvBsocqta+fLuZmsIB4N1EfvXMbagFxgCA8Ng1n2Bj4Fb4D7/gq5HjgT3DWet2CxYKxYAQhvqpk6dmuvxfv36Ee7geb7/TWCxgJeSkqo2bapo3759zi84xUW7NwAgbN2v1lKXFztu5wHxMIGlRg0pOTm435tgBwAICxtTF9j9mhM7bucBKBqCHQAgLGyiRDDPA5AdwQ4AEBY2+zWY5wHIjmAHAAgLK2mS3+BwO27nASgagh0AICysTp2VNMmLHaeeHVB0BDsAQNhYnToraZK15c5eU+oEKD4KFAMAwh7umjdvzsoTQAgQ7AAAYWfdrY0aNeLOA0FGVywAAIBHEOwAAAA8gmAHAADgEQQ7AAAAjyDYAQAAeATBDgAAwCMIdgAAAB5BsAMAAPAIgh0AAIBHEOwAAAA8gmAHAADgEQQ7AAAAjyDYAQAAeATBDgAAwCMIdgAAAB5BsAMAAPAIgh0AAIBHEOwAAAA8gmAHAADgEQQ7AAAAj0iM9AUAADKkpaVp8+bN2r9/vypVqqTk5GSVKMHv4AA8GuyOHDmiM844Q0uWLNGiRYt02mmnpR9bunSphg0bpgULFqhmzZq644479Je//CWi1wsABbVq1SrNmDFDqamp6fsqV66sPn36qEWLFtxIAPmKuV8DLajVrVs32377i7BXr15q2LChFi5cqKefflqjR4/WK6+8EpHrBIDChrqpU6dmCnXGXtt+Ow4Angp2n3zyiWbOnKlnnnkm27HJkyfr6NGjeu2119SqVStdc801uvPOO/Xss89G5FoBoDDdr9ZSlxc7bucBgCeC3fbt23XzzTfrX//6l8qXL5/t+Ny5c9W9e3eVLl06fV/v3r21Zs0a7dmzJ8+uXfuNOHADgHCyMXX5/d1jx+08AIj5YOfz+TR48GANHTpUnTp1yvGclJQU1a5dO9M+/2s7lpsxY8aoSpUq6VuDBg2CfPUAkDebKBHM8wDEr4gGuxEjRighISHPbfXq1ZowYYLzF9oDDzwQ9Guw77lv3770bcuWLUH/fwBAXmz2azDPAxC/Ijor9t5773Va4vLSpEkTffHFF05Xa5kyZTIds9a7AQMG6I033lBSUpLTXRvI/9qO5ca+Z9bvCwDhZCVNbPZrXt2xdtzOA4CoDXZWksS2/Dz//PN6/PHH01//+uuvzvi5KVOmOKVPTJcuXTRy5EgdO3ZMpUqVcvbNmjVLzZs3V7Vq1UL4pwCA4rE6dVbSxGa/5saOU88OgCfG2Nlvqa1bt07fTjnlFGd/06ZNVb9+fed5//79nYkTN954o1asWOGEvvHjx2v48OERvnoAyJ/VqevXr5/TMhfIXtt+6tgB8GSB4tzYxAcrhWIFijt27KgaNWpo1KhRGjJkSKQvDQAKxMKb9TKw8gSAokrw2ZRTpLMxLhYSbSJF1t+cAQAAojl3xERXLAAAAPJHsAMAAPAIgh0AAIBHEOwAAAA8gmAHAADgEQQ7AAAAjyDYAQAAeATBDgAAwCMIdgAAAB5BsAMAAPAIgh0AAIBHEOwAAAA8gmAHAADgEQQ7AAAAjyDYAQAAeATBDgAAwCMIdgAAAB5BsAMAAPAIgh0AAIBHEOwAAAA8gmAHAADgEQQ7AAAAjyDYAQAAeATBDgAAwCMIdgAAAB5BsAMAAPAIgh0AAIBHEOwAAAA8gmAHAADgEYmRvgAgWNLS0rR582bt379flSpVUnJyskqU4HcX7iEAxA+CHTxh1apVmjFjhlJTU9P3Va5cWX369FGLFi0iem2xgnsIALGP5gx4IpBMnTo1U6gz9tr223FwDwEgHhDsEPPdr9ZSlxc7bueBewgAXkewQ0yzMXVZW+qysuN2HriHAOB1BDvENJsoEczz4hH3EAC8g2CHmGazX4N5XjziHgKAdxDsENOspInNfs2LHbfzwD0EAK8j2CGmWZ06K2mSFztOPTvuIQDEA4IdYp7VqevXr1+2ljt7bfupY8c9BIB4keDz+XyRvohoYjMoq1Spon379uXbxYfowsoT3EMAiPfcwcoT8Azrbm3UqFGkLyOmcQ8BILbRFQsAAOARBDsAAACPINgBAAB4BMEOAADAIwh2AAAAHkGwAwAA8AiCHQAAgEcQ7AAAADyCYAcAAOARBDsAAACPINgBAAB4BMEOAADAIxIjfQEAAMSDtLQ0bd68Wfv371elSpWUnJysEiVoX0FwEewAAAixVatWacaMGUpNTU3fV7lyZfXp00ctWrTg/iNo+FUBAIAQh7qpU6dmCnXGXtt+Ow4EC8EOAIAQdr9aS11e7LidBwQDwQ4AgBCxMXVZW+qysuN2HhAMBDsAAELEJkoE8zwgPwQ7AABCxGa/BvM8ID8EOwAAQsRKmtjs17zYcTsPCAaCHQAAIWJ16qykSV7sOPXsECwEOwAAQsjq1PXr1y9by529tv3UsUMwUaAYAIAQs/DWvHlzVp5AyBHsAAAIA+tubdSoEfcaIUVXLAAAgEcQ7AAAADyCYAcAAOARBDsAAACPiJlgZwNOExISMm1jx47NdM7SpUvVrVs3lS1bVg0aNNBTTz0VsesFAAAIt5iaFfvXv/5VN998c45LsNgiyr169VLPnj310ksvadmyZbrhhhtUtWpVDRkyJEJXDAAAED4xFewsyCUlJeV4bPLkyTp69Khee+01lS5dWq1atdLixYv17LPPEuwAAEBciJmuWGNdryeddJLat2+vp59+WsePH08/NnfuXHXv3t0JdX69e/fWmjVrtGfPnghdMQAAQPjETIvdnXfeqQ4dOqh69eqaM2eOHnjgAW3bts1pkTMpKSlq3Lhxpq+pXbt2+rFq1arl+H2PHDnibH779u1L79oFAAAIJX/e8Pl8wfmGvgi6//777U+R57Zq1aocv/bVV1/1JSYm+g4fPuy8vuCCC3xDhgzJdM6KFSuc77Fy5cpcr+GRRx7J9xrYuAd8BvgM8BngM8BngM+AQngP1q1bF5RsleALWkQsvJ07d+q3337L85wmTZpk6l71W7FihVq3bq3Vq1c76+8NHDjQSb3Tpk1LP2f27Nnq0aOHdu/eXeAWu71796phw4bOen5VqlQp1p8vXth9t1nIW7ZsybbINbhvfNYij59R7hmftehlPYXJycnOsDGb8BnTXbE1a9Z0tqKwiRG27l6tWrWc1126dNHIkSN17NgxlSpVytk3a9YsJ/TlFupMmTJlnC0rC3WElMKx+8U9KzzuG/csXPiscc/4rEUvyzRB+T6KATYxYty4cVqyZInWr1/vzIC95557dO2116aHtv79+zstezfeeKPTmjdlyhSNHz9ew4cPj/TlAwAAhEVMTJ6wFrW3335bo0ePdrpNbZKEBbvA0GYtbDNnztSwYcPUsWNH1ahRQ6NGjaLUCQAAiBsxEexsNuz333+f73lt27bVN998U+wQ+cgjj+TYPQvuWTDxWeOehQufNe4Zn7X4+fmM6OQJAAAABE9MjLEDAABA/gh2AAAAHkGwAwAA8AiCXYAXX3xRjRo1UtmyZXXGGWdo/vz5kXtnotDXX3+tiy++WHXr1lVCQkKmYtDGhmvaTOQ6deqoXLly6tmzp37++WfFszFjxqhz586qVKmSU3Pxsssuc9YvDnT48GFnNretg1yxYkVdeeWV2r59u+LVxIkTnYlQ/pprVqPyk08+ST/O/Sr42tr2c3r33Xdz73JhlRbsHgVup556KverALZu3eqUHLO/t+zv+zZt2uiHH35IP86/B9lZvsj6ebPN/v4P5t9tBLv/sbp3Vj7FZqb8+OOPateunXr37q0dO3YU+qZ61cGDB537YgE4J0899ZSef/55vfTSS5o3b54qVKjg3EP7sMarr776yvlBtVndVjDbCmj36tXLuZd+Vrpn+vTpeuedd5zzf/31V11xxRWKV/Xr13dCycKFC51/KGz1mEsvvdSpT2m4X/lbsGCBXn75ZScgB+LeZdeqVStn3XH/9u2333K/8mErJJx11lnOYgD2S9fKlSv197//PdNiAPx7kPPPZeBnzf5NMFdddVVwfz6DsjCZB5x++um+YcOGpb8+ceKEr27dur4xY8ZE9LqilX10/vOf/6S/TktL8yUlJfmefvrp9H179+71lSlTxvfWW29F6Cqjz44dO5x799VXX6Xfo1KlSvneeeed9HNsfWQ7Z+7cuRG80uhSrVo13//7f/+P+1UA+/fv9zVr1sw3a9Ys3znnnOO76667nP181nJeK7xdu3Y53kfuV97rvJ999tm5Huffg4Kxn82mTZs69yuYnzda7CQdPXrUaR2wrsPApT3sta16gfxt2LBBKSkpme6hFY22Lm3uYeY1AU316tWdR/vcWSte4H2zriBbN5D7Jp04ccIpTm4tnNYly/3Kn7UQX3TRRZk+U3zWcmfDRWx4ia1LPmDAAGedcO5X3j788EN16tTJaWmyISbt27fXP/7xD/49KGTuePPNN3XDDTc43bHB/LuNYCdp165dzj8gtWvXznRz7LWFFeTPf5+4h7lLS0tzxjtZF0br1q3T75sthZd14ed4/+wtW7bMGWNiBTuHDh2q//znP2rZsiX3Kx8Wgm0oiY3tzIrPWnb2i+ekSZM0Y8YMZ2yn/YLarVs37d+/n/uVB1va0+5Xs2bN9Omnn+rWW2/VnXfeqTfeeCP9s2b49yB3NkZ97969Gjx4cNB/PmNi5QnAKy0py5cvzzSGBzlr3ry5Fi9e7LRwvvvuuxo0aJAz5gS527Jli+666y5n3I5NAEP++vbtm/7cxiNa0GvYsKGmTp3qTAhA7r+kWovdk08+6by2Fjv7u83GV9vPKvL36quvOp8/ay0ONlrsJGdd2ZIlS2abfWKvk5KSgn7Tvch/n7iHObv99tv10Ucfafbs2c7kgMD7Zk3y9ptboHj/7NlvrieffLKz7rO1PtmknfHjx3O/8mBdOTbZy5ZgTExMdDYLwzahyZ7bb/581vJmrSWnnHKK1q5dy2ctD1b5wFrQA7Vo0SK9G5t/D/K2adMmffbZZ7rppptC8m8Bwe5//4jYPyCff/55pt9I7LWN60H+Gjdu7Hz4Au9hamqqMzs2nu+hzTOxUGddiV988YVznwLZ585mlgXeNyuHYn9BxvN9y8p+Ho8cOcL9ysP555/vdGFbS6d/s1YVGzfmf85nLW8HDhzQunXrnODCz2bubDhJ1rJNP/30k9Paafj3IG+vv/66MzbRxsL6BfXzVqipFh729ttvOzM4J02a5Fu5cqVvyJAhvqpVq/pSUlIifWlRNdtu0aJFzmYfnWeffdZ5vmnTJuf42LFjnXv2wQcf+JYuXeq79NJLfY0bN/YdOnTIF69uvfVWX5UqVXxffvmlb9u2benb77//nn7O0KFDfcnJyb4vvvjC98MPP/i6dOnibPFqxIgRzqzhDRs2OJ8je52QkOCbOXOmc5z7VXCBs2K5d9nde++9zs+mfda+++47X8+ePX01atRwZq9zv3I3f/58X2Jiou+JJ57w/fzzz77Jkyf7ypcv73vzzTfTz+Hfg5xZxQ37+95mFmcVrL/bCHYBJkyY4NzU0qVLO+VPvv/++0LfUC+bPXu2E+iyboMGDXKO25Tthx9+2Fe7dm0nJJ9//vm+NWvW+OJZTvfLttdffz39HAu+t912m1PSw/5yvPzyy53wF69uuOEGX8OGDZ2fw5o1azqfI3+oM9yvogc77l1mV199ta9OnTrOZ61evXrO67Vr13K/CmD69Om+1q1bO3/Xn3rqqb5XXnkl03H+PcjZp59+6vwbkNO/jcH6+Uyw/xSujQ8AAADRiDF2AAAAHkGwAwAA8AiCHQAAgEcQ7AAAADyCYAcAAOARBDsAAACPINgBAAB4BMEOAADAIwh2AAAAHkGwA4D/GTx4sBISEpzNFuSuXbu2LrjgAr322mtKS0vjPgGIegQ7AAjQp08fbdu2TRs3btQnn3yi8847T3fddZf+8Ic/6Pjx49wrAFGNYAcAAcqUKaOkpCTVq1dPHTp00IMPPqgPPvjACXmTJk1yztm8ebMuvfRSVaxYUZUrV1a/fv20ffv29O+xZMkSJxBWqlTJOd6xY0f98MMP6ce//fZbdevWTeXKlVODBg1055136uDBg7wPAIqNYAcA+ejRo4fatWun999/3+mStVC3e/duffXVV5o1a5bWr1+vq6++Ov38AQMGqH79+lqwYIEWLlyoESNGOF27Zt26dU6r4JVXXqmlS5dqypQpTtC7/fbbeR8AFFuCz+fzFf/bAIA3xtjt3btX06ZNy3bsmmuucYLY+PHj1bdvX23YsMFpbTMrV65Uq1atNH/+fHXu3NlppZswYYIGDRqU7fvcdNNNKlmypF5++eX0fRbszjnnHKfVrmzZsiH+UwLwMlrsAKAA7Hdgm1SxatUqJ9D5Q51p2bKlqlat6hwzw4cPdwJcz549NXbsWKeVLrCb1rp0rRvXv/Xu3dtpCbSwCADFQbADgAKw0Na4ceMC3avRo0drxYoVuuiii/TFF184we8///mPc+zAgQO65ZZbtHjx4vTNwt7PP/+spk2b8l4AKJbE4n05AHifhbNly5bpnnvuccbObdmyxdkCu2KtC9cCnN8pp5zibPY1f/rTn/T666/r8ssvdyZk2Pknn3xyBP9EALyKYAcAAY4cOaKUlBSdOHHCmek6Y8YMjRkzxil3MnDgQJUoUUJt2rRxJkiMGzfOKYFy2223OWPkOnXqpEOHDum+++7TH//4R6eF75dffnEmUdhkCXP//ffrzDPPdCZLWHdthQoVnKBnkzBeeOEF3gsAxUKwA4AAFuTq1KmjxMREVatWzZkN+/zzzzsTISzUGSt/cscdd6h79+7OPpvlapMljE2M+O2335wQaMGwRo0auuKKK/Too486x9u2bevMph05cqRT8sTG7lkXbOCsWgAoKmbFAgAAeASTJwAAADyCYAcAAOARBDsAAACPINgBAAB4BMEOAADAIwh2AAAAHkGwAwAA8AiCHQAAgEcQ7AAAADyCYAcAAOARBDsAAACPINgBAADIG/4/bXtXPi+Q9fIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, 70, 100)\n",
    "y = estr.theta[0] * x / (estr.theta[1] + x)\n",
    "\n",
    "plt.plot(d, r, 'o', color='gray')\n",
    "plt.plot(x, y, '-', color='blue')\n",
    "plt.xlabel(\"Dose\")\n",
    "plt.ylim([-50, 0])\n",
    "plt.ylabel(\"Response\")\n",
    "plt.xlim([0, 70])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c70987a8-11e3-48f9-b41d-4df457c3a297",
   "metadata": {},
   "source": [
    "While not discussed in the chapter, one might be concern about outlying response values. Particularly, there are some large responses at low doses. These few response values might be overly influential on the overall E-max parameter estimates. To go beyond the book, we can consider the use of outlier-robust estimating functions. These apply an additional constraint that limits the maximum influence an observation can have to the estimating function. The following code uses the Huber function to limit the influence. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "19f573c0-ec45-4f44-85ed-3e13fdbc7a7a",
   "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>Est</th>\n",
       "      <th>SE</th>\n",
       "      <th>LCL</th>\n",
       "      <th>UCL</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Params</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>E-max</th>\n",
       "      <td>-35.5</td>\n",
       "      <td>4.2</td>\n",
       "      <td>-43.7</td>\n",
       "      <td>-27.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ED50</th>\n",
       "      <td>6.3</td>\n",
       "      <td>2.1</td>\n",
       "      <td>2.3</td>\n",
       "      <td>10.3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Est   SE   LCL   UCL\n",
       "Params                       \n",
       "E-max  -35.5  4.2 -43.7 -27.4\n",
       "ED50     6.3  2.1   2.3  10.3"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def psi_robust(theta):\n",
    "    lower = 0.\n",
    "    vals = [lower, ] + list(theta)\n",
    "    return ee_emax(theta=vals, dose=d, response=r, robust='huber', k=5.)[1:, :]\n",
    "\n",
    "\n",
    "estr = MEstimator(psi_robust, init=[-40, 5])\n",
    "estr.estimate(solver='lm')\n",
    "\n",
    "results = pd.DataFrame()\n",
    "results['Params'] = [\"E-max\", \"ED50\"]\n",
    "results = results.set_index('Params')\n",
    "results['Est'] = estr.theta\n",
    "results['SE'] = np.diag(estr.variance) ** 0.5\n",
    "ci = estr.confidence_intervals()\n",
    "results['LCL'] = ci[:, 0]\n",
    "results['UCL'] = ci[:, 1]\n",
    "results.round(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "081fb932-01fc-4cf4-be72-8a21b2fe6560",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAV0FJREFUeJzt3Qd4VGXaxvE7CRBa6CX03puAiFhQASl2UbGgIBYUWURlXcWKqyu2VVjxs60Fu7gqCquUtWChCEhv0kFK6EloAZL5rucdJo2QxiRT8v9d1+s5c85JPJxMkjtvjfB4PB4BAAAg5EUG+gYAAADgHwQ7AACAMEGwAwAACBMEOwAAgDBBsAMAAAgTBDsAAIAwQbADAAAIEwQ7AACAMEGwAwAACBMEOwAAgDARlsHulVdeUf369VWyZEl17txZv/32W6BvCQAAoMCFXbD79NNPdd999+nxxx/X77//rnbt2qlXr17asWNHoG8NAACgQEV4PB6PwojV0HXq1Enjxo1zr1NSUlSnTh0NGzZMDz74YKBvDwAAoMAUUxg5cuSI5s+fr5EjR6Yei4yMVI8ePTRr1qwsPyYpKckVHwuCe/bsUeXKlRUREVEo9w0AAIomj8ejxMRE1axZ02WWUxVWwW7Xrl1KTk5W9erVMxy31ytXrszyY0aPHq0nnniikO4QAADgRJs3b1bt2rV1qsIq2OWH1e5Znzyf+Ph41a1b1z3gcuXKBfTeAABAeEtISHBdxmJiYvzy+cIq2FWpUkVRUVGKi4vLcNxex8bGZvkx0dHRrmRWukwZgh0AACgU/ur+FVajYkuUKKGOHTvqu+++y9Bnzl536dIlT59r9h9rCuAOAQAACk5YBTtjzapvvvmmxo8frxUrVmjIkCE6cOCABg0alKfPM33xogK7RwAAgIIQVk2x5tprr9XOnTv12GOPafv27TrttNM0ZcqUEwZU5GTu5sUFdo8AAAAFIezmsfNHJ8by5cur+j3naftLPwb6dgAAQBHIHfHx8X7p2x92TbH+sjNqsZtbBgAAIFQQ7E4ipfhebdy3uXC/GgAAAKeAYJeNKYsWnMqzBQAAKFQEu2z8uJJgBwAAQgfBLhuL4gh2AAAgdBDssrHpKMEOAACEDoJdNg4W36zdB3cX3lcDAADgFBDsTmZvA7eZ++fCU3m+AAAAhYZgdxJRe9q67fQlNMcCAIDQQLA7iRoR3mA3eyPBDgAAhAaC3Uk0r9jObVclEOwAAEBoINidRJeG3hq73Vqlg0cPFubXBAAAIF8IdifRuVV1KTFWikjR4rjF+Xu6AAAAhYhgdxItW0ra3t7tz95AcywAAAh+BLuTqFxZKpPoDXY//kGwAwAAwY9gl42GpbzBbtF2gh0AAAh+BLtsdKzlDXabjyzR0eSjhfU1AQAAyBeCXTbOatFAOlxOyRFJWrlrJW8xAAAQ1Ah22WjXNlLafprbX0BzLAAACHIEuxxHxnqD3cz19LMDAADBjWCXjbJlpSrHvP3s5mxcWFhfEwAAgHwh2OWgZSVvsFsVv0Aejyd/TxkAAKAQEOxy0KVxS+loSR3yxOuP3X8UxtcEAAAgXwh2OWjXuri0pbPb/2XTL/l7ygAAAIWAYJeD1q0lbTrH7RPsAABAMCPY5aBZMylqizfY/bCOGjsAABC8CHY5KFFCalupi+SJ0MbENYrbH1c4XxkAAIA8Itjlwtkdy0txbdz+r5t/zeszBgAAKBQEu1w44wz62QEAgOBHsMtjsPt5I/3sAABAcCLY5UKTJlLMXm+wW7D9dx04cqCgvy4AAAB5RrDLzUOKlM5sWUfaV1fJnmTN2TIn708aAACggBHscol+dgAAINgR7HKJYAcAAIIdwS4fwW7W5lk6lnKsIL8uAAAAeUawy6XYWKlOyVbS4fLaf3S/lsQtyfvTBgAAKEAEuzzofEaktPkst8+6sQAAINgQ7PLbz24z89kBAIDgQrDLg86dM65A4fF4CurrAgAAkGcEuzzo0EGK2NZJSi6urYlbtWHfhrw/cQAAgAJCsMuDsmWl1s1KSVtPd6/pZwcAAIIJwe4Um2MBAACCBcEujxhAAQAAghXBLj/B7viUJ8t3Ltfug7sL4usCAACQZwS7PGrVSiqtKtLO5u71zM0z8/7UAQAACgDBLo+KFZM6dkzrZzdj44yC+LoAAADkGcEuv82x63q4/Wlrp+XnUwAAAPgdwe5Ugp0nUkt2LNGWhC3+/8oAAADkEcEuv1OeHKosbenkXk9dOzU/nwYAAMCvCHb5ULeuVK2apDW93espa6b496sCAACQDwS7fIiION4cezzYTV83XcdSjuXnUwEAAPgNwS6funSRa4otfqyi9h3ep9+2/Oa/rwoAAEA+EOzy6YILJHmiFLGup3tNcywAAAg0gl0+nX66VKaMdGQ5/ewAAEBwINjlU/HiUteu1s+ul3s9b+s87Tq4y59fGwAAgDwh2J1qc+z+Gip3sJ088mj62umn8ukAAABOCcHuFHTr5t0eXnq8OXYt054AAIDAIdidgtNOkypUsH523ubYqWumKsWT4q+vDQAAQJ4Q7E5BVJR03nmSNp+tEiqjuANxWrR90al8SgAAgHwj2Pmjn11yCZXf0929ZtoTAAAQKAQ7fwQ7Sfvm0c8OAAAEFsHuFLVuLVWpIh1d4e1nN3PzTMUfjvfH1wYAACBPCHanKDJSOv98SXsbqrKaujVjv1///al+WgAAgDwj2Plx2pNSf7IKBQAACByCnR/72W3/Na2fncfj8cenBgAAyDWCnR80aybVqCEdW3OeoiNLaVP8Ji3YvsAfnxoAACDXCHZ+EBFxvNbuaGk1SL7IHfts2Wf++NQAAAC5RrDzc3NsypJr3Paz5Z/RHAsAAAoVwc7PAyjWTblYJYuV1Nq9a7Vw+0J/fXoAAIDwCXb/+Mc/dNZZZ6l06dKqYAu0ZmHTpk26+OKL3TXVqlXT/fffr2PHjhXK/TVoINWtKx07WFYdyx1vjl1OcywAACg8IRPsjhw5omuuuUZDhgzJ8nxycrILdXbdzJkzNX78eL377rt67LHHCq2fna/WrtJ2mmMBAEDhC5lg98QTT+jee+9VmzZtsjw/bdo0LV++XB988IFOO+009enTR08++aReeeUVF/YKs5/d1h8ucc2xa/as0aK4RYXy/wYAAAiZYJeTWbNmudBXvXr11GO9evVSQkKCli1bVqjBbsGcsupRz9scO2HZhEL5fwMAAIRNsNu+fXuGUGd8r+3cySQlJbnwl77kV506klUopqRI9fbTHAsAAIpQsHvwwQcVERGRbVm5cmWB3sPo0aNVvnz51FLH0tkpuPRS73brjzTHAgCAIhTsRowYoRUrVmRbGjZsmKvPFRsbq7i4uAzHfK/t3MmMHDlS8fHxqWXz5s2n9G+67DLv9n/flFWvhn3cPpMVAwCAwlBMAVS1alVX/KFLly5uSpQdO3a4qU7M9OnTVa5cObVs2fKkHxcdHe2Kv3TqZE3AFiqllrpGX+lLN+3JU92ecjWQAAAAKup97GyOuoULF7qtTW1i+1b279/vzvfs2dMFuJtuukmLFi3S1KlT9cgjj2jo0KF+DW45iYyULrnEu79n1iWKjorW6j2rtThucaHdAwAAKJpCJtjZfHTt27fX448/7sKc7VuZN2+eOx8VFaXJkye7rdXe3XjjjRowYID+/ve/F/q9+ppjp06KUZ/Gx5tjmawYAAAUsAiPx+Mp6P9JKLFRsTaIwvrbWTNufhw8KFWuLB0+LI3+70caObe/mlZuqpVDV9IcCwAA/Jo7QrLGLpSULi316OHdP/j7pa459o/df2jJjiWBvjUAABDGCHYFxDftybTJMerTxNsc+9GSjwrqfwcAAECwKyi+ARRz5kgX17nR7b+36D0dSznG2w4AABQIauwKSM2a3qlPTPLyS1W1dFVt279N367+tqD+lwAAoIgj2BVCc+w3k0poQLsBbv+tBW8V5P8SAAAUYQS7Qpj2ZPp0qX/LW93+5D8ma1vitjx9npSUFG3YsEFLlixxW3sNAAAQVCtPhLu2baW6dW1yZWnLwhbqUruLZv05y/W1e+CcB3L1OWxZtSlTprjh0D42HLp3795q0aJFAd49AAAINdTYFSBbQczXHPv119Kt7b21dm8vfFu5mT7QQt2ECRMyhDpjr+24nQcAAPAh2BVSc+zkydI1La9V2RJl3Zx2v2z6JduPs+ZWq6nLjp2nWRYAAPgQ7ArYeedJZctK27ZJq5aU1bWtrnXH/73g39l+nK2Jm7mmLjM7b9cBAAAYgl0Bi46Wevf27n/xRVpz7GfLPlP84fiTflxiYmKuPn9urwMAAOGPYFcIrrnGu/34Y+mMmmeqZdWWOnTskD5Z+slJPyYmJiZXnzu31wEAgPBHsCsENoDC8tfGjdLs2RGptXbZzWlXt27dHBcDtvN2HQAAgCHYFYJSpaS+fb37H34o3dT2JhWPLK65W+dqcdziLD8mMjLSTWmSHTtv1wEAALj8wGMoHDfc4N1OmCBVKFFVlzXzDpd96/eT19rZPHVnnXWWImzelHTstR1nHjsAAJAewa6QdOsmVa8u7d4tTZ0q3dbhNnf8/cXv6+DRg1l+jM1TN3PmzBPmvLPXdpx57AAAQHoEu0JSrJh03XXe/Y8+ki5seKEaVmyovYf3avzC8Sdczzx2AAAgrwh2hah/f+/2q6+kQwejdE/ne9zrl2a/pBRPxvVfmccOAADkFcGuEJ1+utSkiXTwoDRxojSo/SBVKFlBq/es1uQ/Jme4lnnsAABAXhHsCpGNgfANorDRsba82B0d73Cv/znrnxmuZR47AACQVwS7ADXHTp8u7dghDTtjmIpFFtNPG3/SvK3zUq9jHjsAAJBXBLtCZk2xnTpJycneqU9qlaul61p7R1W8OOvFtC8M89gBAIA8ItgFsNbOmmPNiC4j3HbCsgnaFL8p9Tqbp65fv34nrEBhr+0489gBAID0IjyZJ0kr4hISElS+fHnFx8fnuKRXfm3fLtWqZVOaSGvWSI0aSd3f667v13/vQt4LPV84YeoTGyVrAyqs750107LiBAAAoS/Bz7mDGrsAiI2VundPm9PO3HfmfW775u9vKiEpIcP1FuLq16+vNm3auC2hDgAAZIVgFwTNsVZn2qdJHzWv0tyFuuyWGQMAADgZgl2A9O0rlSkjrVol/fyzFBkRmVprN3bOWB1LORaoWwMAACGKYBcgMTFpc9q9+qp3e2PbG1W1dFVtjN+oz5Z9FqhbAwAAIYpgF0BDhni3n38uxcVJpYqXcvPamSd/elLJKcmBvD0AABBiCHYB1L691LmzdPSo9Pbb3mN3d75bFUtW1IpdK/Tx0o8DeXsAACDEEOyCpNbu9de9kxaXL1le9591vzv2xIwn6GsHAAByjWAXYP36SRUrShs3SlOmeI8N6zzM9bVbs2eN3lv0XqBvEQAAhAiCXYCVKiUNGpRxEEXZEmX1wNkPuP2/z/i7jiQfCeAdAgCAUEGwCwJ33undfvONtH69d39IpyGKLRvrRsi+veB4BzwAAIBsEOyCQJMmUo8e3omK33jDe6x08dJ66JyH3P5TPz2lw8cOB/YmAQBA0CPYBdkgirfekpKSvPuDOw5WnXJ1tCVxi96YfzzxAQAAnATBLkhcdplUs6a0c6f0xRfeY9HFovVI10fc/tM/P62DRw8G9iYBAEBQI9gFiWLFpNtvzziIwgw6bZAaVGiguANx+r+5/xew+wMAAMGPYBdELNhFRXnXjl261HuseFRxPXbeY27/mV+e0d5DewN7kwAAIGgR7IJIrVreJlkzZkzacVtDtlXVVtp9aLdG/TgqYPcHAACCG8EuyIwY4d2+9560ZYt3v1hkMY3tPdbtvzL3FS3bsSyAdwgAAIIVwS7InH22dO653vVjX3op7Xj3ht11RfMrlOxJ1r1T75XH5kYBAABIh2AXhB58MG392D170o7/s+c/FR0VrenrpuvrVV8H7P4AAEBwItgFoT59pDZtpP37pf9LNxC2YcWGGtHF21Z737T7lHTs+IR3AAAABLvgFBGRVms3dqx0MN30dSPPHamaMTW1bu86vTQ7XVstAAAo8qixC1L9+kkNGki7dnlXo/ApW6Ksnu3xbOpSY1sTtwbuJgEAQFAh2AXxhMX33+/df+EF72AKnxva3KAza5+pA0cPaOR3IwN2jwAAILgQ7ILYzTdL1apJmzZJn3ySdjwyIlL/6v0vt//eovc0a/OswN0kAAAIGgS7IFaqlHTPPd79Z56RUlLSznWq1cktN2Zum3QbAykAAADBLtjddZdUrpy0fLk0eXLGc89f+Lyqlamm5TuX6x8//yNQtwgAAIIENXZBrnx5acgQ7/7TT0vp5yWuXLqyXrnoFbc/+pfRWrR9UYDuEgAABAOCXQiw5tiSJaU5c06stbu65dXq26KvjqUc0y1f3+K2AACgaCLYhYDYWGn4cO/+Qw9JyckZz1utXcWSFfX7tt/1wswXAnKPAAAg8Ah2IeKBB6QKFaSlS6UPP8x4LrZsrMb0HuP2R/04Sit3rXT7KSkp2rBhg5YsWeK29hoAAISvCA+ryWeQkJCg8uXLKz4+XuVs1EIQefZZ74oU9epJq1ZJ0dFp5+zLeNFHF2nKmik6q85Zer3z65o+bbr79/jYv6d3795q0aJFYP4BAACgQHMHNXYhZNgwqWZNaeNG6bXXMp6LiIjQ65e8rpgSMZq5eabu/8/9GUKdsdcTJkzQihUrCvfGAQBAoSDYhZDSpaXHH/fuP/WUBbWM5+uWr6tnejzj9v+n/2mHdmT5eaZMmUKzLAAAYYhgF2JuuUVq2tS7huyLL554vneV3mqkRjqmY/pMn+mo0q1Flq7mbpMtZwEAAMIKwS4E15D9x/G5iP/5T2lHpkq5A/sP6EpdqTIqo53aqamamuXnSUxMLIS7BQAAhYlgF4Kuuko6/XRp//60kOcTExOjsiqrvurrXs/TPC3X8hM+h10HAADCC8EuBEVEeNeONa++Kq1fn3aubt26blSNNceerbPdsa/1tfZpX+o1dt6uAwAA4YVgF6K6d5d69JCOHpX++te045GRkW5KE9NN3VRLtXRYh/W5PleyvDMb23m7DgAAhBd+u4ewl16SoqKkL76wka5px22eun79+qliuYq6WlcrWtHarM2aXWK2O848dgAAhCeCXQhr3TptqTGb4y4pKe2chbfhw4frnoH36MkznnTH/nfkf9oSvSVAdwsAAAoawS7E2bx2NWpIa9ZIL2RaJtaaW+vXr6/7+9yv29rfJo886vdZP63dszZQtwsAAAoQwS7E2eojNu2JsRGyGzZkfd3LF72szrU6a+/hvbr8k8uVmMR0JwAAhBuCXRi47jrp/POlQ4eke+/N+pqSxUrqi2u/UI2yNbRs5zINnDhQKZ6Uwr5VAABQgAh2YTL9ybhx3smLJ06Uvvkm6+tqxtTUl9d+qRJRJfTlyi/15Axv3zsAABAeCHZholUr6Z57vPt33y0dPpz1dZ1rd9ZrF7/m9kfNGKWJKycW4l0CAAAV9WC3YcMG3XrrrWrQoIFKlSqlRo0a6fHHH9eRI0cyXLd48WKde+65KlmypOrUqaPnnntORcljj0k1a0pr10rPPnvy6wa1H6Thnb3DaW/68iYt3bG08G4SAAAU7WC3cuVKpaSk6PXXX9eyZcv00ksv6bXXXtNDDz2UYWH7nj17ql69epo/f76ef/55jRo1Sm+88YaKClsl7MUXvftPPy0tP3ElsVQv9HxB3Rp00/4j+3Xpx5dqW+K2QrtPAABQMCI8Ho9HIciC26uvvqp169a517b/8MMPa/v27SpRooQ79uCDD2rixIkuGOaWBcTy5csrPj7eLb0Vauyreeml0n//611PdtYsb9+7rOw+uFtnvnWm1uxZo9NiT9OMm2eoXHTo/ZsBAAhVCX7OHSFRY5cVewCVKlVKfT1r1ix17do1NdSZXr16adWqVdq7d6+K0kAKq6SsUEGaNy/7JtnKpStr6o1TVb1MdS3cvlB9P+2rpGPpZjkGAAAhJSSD3Zo1a/Tyyy/rjjvuSD1mNXXVq1fPcJ3vtZ07maSkJJeW05dQZ/3sXn7Zu//EE9b38OTXNqzYUN/0/0ZlS5TVd+u/081f3cw0KAAAhKiABjtrKo2IiMi2ZG5G3bJli1vE/pprrtHtt99+yvcwevRoVwXqKzboIhz07y9dcYV09Kg0cKCUaZxJBh1qdNAX/b5Qschi+mTpJxoxdYRCtIUeAIAiLaB97Hbu3Kndu3dne03Dhg1Tm1e3bt2q888/X2eeeabeffddt2SWz4ABA1xtm/Wp8/nhhx/UrVs37dmzRxUrVjxpjZ0VH/scFu5CtY9denFx3mlQ7BHb0mOjRmV//YeLP9SNX97o9p+/8Hn99ay/Fs6NAgBQRCX4uY/dSbrVF46qVau6khtWU3fBBReoY8eOeueddzKEOtOlSxc3eOLo0aMqXry4OzZ9+nQ1a9bspKHOREdHuxKOrCX6lVe8K1PYcmOXXSZ16HDy6/u37a/t+7frr9P/qvun36+KJSvq1g63FuYtAwCAcO9jZ6HOaurq1q2rF154wdX0Wb+59H3nbrjhBlezZ/Pd2ZQon376qcaOHav77rtPRVm/ftLVV0vHjkk332w1lNlfP+KsERrRZYTbv33S7Xp34buFc6MAACC0a+xyy2rebMCEldq1a2c452tJtmrMadOmaejQoa5Wr0qVKnrsscc0ePBgFWU2Svb//k+aMUNassT6NUovvZT9x1gzrI2OHTd3nG756hZFRkRqQLsBhXXLAACgqM1jV1BCfR67k5k0ydsUa6wb4uWXZ3+9vS2GfjNUr857VRGK0PtXvu+aagEAgP8wjx3yxSYt9rVKW5Pshg3ZX28jksddNE6DOwyWRx4NmDhAHy/5mKcPAEAQC4k+dvCP0aOlM86Q9u3zDqjIbgoUY02wr17yqm5rf5ub285GzH669FO+HAAABCmCXRFis8Z8+ql3VYo5c6SRI3P+GAt3r1/6ugadNsiFu+s/v15vzn+zMG4XAADkEcGuiKlfX3rnHe/+iy9KX3+du3D35qVvpjbLDp48WM/88gyTGAMAEC7Bbu3atXrkkUd0/fXXa8eOHe7Yt99+66YaQXCzFSnuuSetv93GjTl/TFRklF675DWNPMdbzTfyu5H62/S/Ee4AAAj1YDdjxgy1adNGc+bM0RdffKH9+/e744sWLdLjtsQBgt6zz0qdOkl790rXXCMdOpTzx9iAiqe7P60XLnzBvX5h1gu69etbdSzlWMHfMEJWSkqKNmzYoCVLlritvQYABNF0J7bKg63VapP/xsTEuEBnS3/99ttv6tu3r/7880+FqnCd7iQr69dLHTt6w52tLfv++95573LjnQXv6LZJ3kEVVzS/Qh/1/Uilipcq6FtGiFmxYoWmTJnivq987PvK1ntu0aJFQO8NAIJBUEx3Yn95X3nllSccr1atmnbt2nXKN4XC0aCB9J//SFFR0ocfemvxcmtQ+0H6zzX/UYmoEpq4cqIuGH+B4vbHFeTtIgRD3YQJEzKEOmOv7bidBwD4V76CXYUKFbRt27YTji9YsEC1atXyx32hkHTrJr38snf/oYekr77K/cde2eJKTbtxmltTds6WOer8785aumNpgd0rQoc1t1pNXXbsPM2yABAEwe66667TAw884NZqtX5X9sP5119/1V//+lcNGMDSU6FmyBDprrtstQlvk+zixbn/2PPqn6fZt81W40qNtTF+o8566yxNWZP9L3SEv02bNp1QU5eZnbfrAAABDnZPP/20mjdvrjp16riBEy1btlTXrl111llnuZGyCD1jxkjdu0sHDniXHtu5M/cf27RyU82+dbbOq3eeEo8k6uKPLta438YV5O0iyCUmJvr1OgBAAQa7EiVK6M0339S6des0efJkffDBB1q5cqXef/99RVmHLYSc4sWlCROkxo2905/07SslJeX+4yuXrqxpN01Lnch42LfDdNd/79KR5ByWt0BYskFV/rwOAFAIExRbjd1FF12kq666SgcOHNBeG16JkFWpkjRpklS+vPTLL9LAgdZXKvcfbwMp3rrsLT3bwzsK49V5r+q8d8/TnwmhO0oa+VO3bt0cR3fZebsOABDgYHfPPfforbfecvvJyck677zz1KFDBxf0fvzxRz/eHgpb8+bekbJWg2fLjw0f7u17l1vW5/JvZ/9Nk6+frAolK2j2n7PV4fUO+mH9DwV52wgykZGRbkqT7Nh5uw4A4D/5+qn6n//8R+3atXP7kyZNck2y1hR777336uGHH/bj7SEQevRIm9Nu3Djpqafy/jkubnqx5g+er9NiT9POgzvV4/0eeu7X51ipogixeer69et3Qs2dvbbjzGMHAEEyQXHJkiW1Zs0a1a5dW4MHD1bp0qU1ZswYrV+/3gW+nEbDBbOiNEFxTizUDRvm3X/tNemOO/L+OQ4dPaQh/x2i8YvGu9dXNr9Sb1/+tqvNQ9Fgo+Zt9KsNlLA+ddb8Sk0dAATRBMXVq1fX8uXLXTOszUV14YUXuuMHDx5k8EQY+ctfpEcfTZsSxZpo88pWo3jn8nf02sWvuT54X678Uu1ea6efN/7s9/tFcLIQV79+fbcMoW0JdQBQcPIV7AYNGuSaUlq3bu36VPWwtjvJrR1r06AgfDzxhLemzjfH3fff5/1z2HvkjtPv0K+3/KpGFRtpU/wmnT/+fD3+w+OsMwsAQKCbYn397DZv3uzWjLUmWTN+/Hi3KsXll1+uUEVT7ImSk6Vrr5U+/1wqXVr69lupa9f8Pd/EpEQ3FYqvabZL7S76sO+HalCxwal94QAACEH+zh35DnbhimCXtcOHpSuukKZOPfVwZz5e8rHu/O+dSkhKULnochrXZ5xubHujq90DAKCoSAiWYPfdd9+5smPHjhPWe3z77bcVqgh2hRfuNuzboP5f9NfMzTPd68uaXeb64tWIqZH/TwoAQAgJisETTzzxhHr27OmC3a5du9zExOkLwlPJktLEiVKvXjZQRurTR/rpp+w/xkL/hg0btGTJErdN/0dA/Qr1NePmGXrqgqdUPLK4vl71tVr+X0u9v+h9pkUBAKCwauxq1Kih5557TjfddJPCDTV2/qu5W7FihRs1nX76G/trxCamzTyH2ZK4Jbr5q5v1+7bf3etLm16q1y95ndo7AEBYSwiGGrsjR47orLPOOuX/OcKn5m769BND3YQJE06Y09Be23E7n16b6m00+9bZqbV3k/6Y5Grv3pj/hlt7FgAAFFCwu+222/TRRx/l50MRpuHu4oulCRO856y51WrqsmPnM/fNLB5VXA93fVi/3/G7OtboqH2H9+mOyXfonLfP0eK4xQX5zwEAICwUy88HHT58WG+88Yb+97//qW3btipuC4um8+KLL/rr/hDk4e6rryRrkf/sM+m666Rdu6SLLtqU4+ojdt5WI7AJazNrXa21Zt82W6/89ooe+eERzfpzlltv9t4z79Xj5z+usiXKFuC/CgCAItbH7oILLjj5J4yI0Pf5mcU2SNDHLn/z3NnSY6++6n09ZEicqlV7za01m52+ffu61QiysyVhi4ZPGa7PV3zuXtcpV0djeo9xS5MxNQoAINQlBMt0J+GKYJc/9i6yVSqsmE6d5qpPn28VGXnyt9fAgQOzrLHLyjerv9HQb4a6KVLM+fXP15heY9Qutl0+7xgAgMALisET6f3555+uoGiz2rlRo6SXX7Z9j+bO7aTPPrtaR45k3dpvb15bDD63LmpykZbdtUyPnPuIShYrqR83/KgOb3TQnZPv1M4DO/34LwEAIHTlK9hZp/e///3vLmHWq1fPFVtK7MknnzyhQzyKlr/8Rfr44wgVK+bRihUt9e67Nysx8cQ+cTblSV4Xgy9dvLSe7PakVg5dqX6t+rnRsq/Pf11NXm6iF2e9qKRjSX78lwAAUESC3cMPP6xx48bpmWee0YIFC1x5+umn9fLLL+vRRx/1/10ipNi6st99F6EKFY5p69ZaevPN27RtW2xqTV2/fv1OmMcuL+pVqKdPr/5UP938k9rHtld8UrxGTBuhZuOaucmNk1OS/fivAQAgdOSrj13NmjX12muv6bLLLstw/KuvvtJdd92lLVu2KFTRx85/1q6VLrnEo5UrI1SyZIrGjt2l226rkueauuxYiHt34bt67MfHtDVxqzvWtnpbje4+Wn0a92GABQAgqAVFH7s9e/aoefPmJxy3Y3YOMI0aSbNmRejCC22KnEjdeWc1Pf98pBto4S9RkVG6tcOtWj1stZ7p/ozKR5d3c95d/NHFOn/8+fp10698MQAARUa+gl27du1cU2xmdszOAT4VKkjffCPddZd35OyDD3rnu0tM9O8zsv53D5zzgNYNX6f7z7pf0VHR+mnjTzrnnXPU8/2emrV5Fl8UAEDYy1dT7IwZM3TxxRe7UY1dunRxx2bNmqXNmzfrm2++0bnnnqtQRVNswXnlFemee6RjxyTrYvf5595tQdgcv1lP/vSk3ln4jo6lHHPHejbqqVHnjVKXOt73LAAAgRY089ht3bpVr7zyilauXOleW2d4619n/e9CGcGuYM2cKV1zjb1/pLJlpbfekvr1K7j/3/q96/X0z0/r3UXvZgh4D53zkLrW60ofPABAQAVNsAtXBLuCFxfnbY798Ufva6vFe+45KdPKdH61bu86/eOnf2j8ovFK9nhHzXap3UUPnvOgLml6iSIj/DegAwCAkAt2e/fu1VtvvaUVK1a41y1bttSgQYNUqVIlhTKCXeGw5thHHpGefdb7+swzpY8+kho0KNj/rwW853993jXRJiV7571rVbWVHjj7AV3X+joVjyrAdBmGbN5KW/M3MTFRMTExrnuGP0c9I7iecTDdCxAuEoIh2P3000+69NJL3Y2cfvrp7tj8+fO1b98+TZo0SV27dlWoItgVrokTbWkxe+42x513vdkbbij4/++2xG0aM3uMXp33qhKPeEdy1C5XW8POGKbbO9yuiqUqFvxNhDj7o27KlCnue8bHfijZ5NOnMk8hgvMZB9O9AOEkIRiCnS3cboMmXn31VUVFRbljycnJro/dzJkztWTJEoUqgl3h27BB6t/f2//O3HSTjbD2Br2Ctu/wPr0691WNmTNGOw7scMfKFC+jQacN0vAzh6txpcYFfxMhyH7JT5gw4aTnT3USagTXMw6mewHCTUIwzGO3Zs0ajRgxIjXUGdu/77773DkgL+rXt5HW3rVmrVXn/fel9u2lOXMK/jlWKFlBI88dqY33bNTbl72tNtXa6MDRAxo3d5yavtxUl39yuaatneaWL0Nac5zV3GTHzrO8YHg842C6FwAFFOw6dOiQ2rcuPTvGPHbIj2LFpMcft2Z+qV49ad066eyzvf3wkgphCdiSxUpqUPtBWnTnIk2/abouanKRPPLo61Vfq9cHvdR8XHO9NOsl7T20V0Wd9bFK3xyXFTtv1yH0n3Ew3QuAAgp2d999t4YPH64XXnhBv/zyiyu2f++997qyePHi1ALkhYW5hQul66+35n3pH/+QrBvn/PmF8xwjIiLUo2EP/feG/2rF0BW6+4y7VS66nFbvWa37pt2nWi/W0q1f3aq5W+aqqA4ot47z/rwOwf2Mg+leABRQH7ucRkHZL0f7tLa1vnehhD52wcMmMB4yRNq505r6vatWPPqoFB1duPex/8h+fbj4Q70y9xUt2ZHWf9TWpLWBFv3b9C9Sgy02bNig8ePH53jdwIEDVd/a2RHSzziY7gUIRwnB0Mdu/fr12ZZ169alboH8uuoqafly6dprM9bezZ1buM+0bImyuuP0O1wz7c+DftYNbW5wS5bZmrTDvh2mmi/W1E1f3qQfN/xYJPri2RQXOf3wsfN2HUL/GQfTvQDIGRMUZ0KNXfDX3kVESEOHSk89JZUvH5j72XNoj6vFe/P3NzPU4tWvUF83tb1JA9oNCOsRtYySLFrPOJjuBeGBORGDbLoTq5avUqWKWy/W/O1vf9Mbb7zhJin++OOPVc96v4cogl3w2rVLuvde6YMPvK9r1JDGjpWuvtob9gLBvn3mbZ3nAt6nyz5VQlJaJ/Oz65ytge0G6uqWV4dlUy3zmhWtZxxM94LQxnspCINds2bN3Bx23bp106xZs9S9e3eNGTNGkydPVrFixfTFF18oVBHsgt9333lr71av9r7u08c7713DhoG9r0NHD2niyolu2bLp66anNsuWiCrhRtne0PoGt3xZqeKlFC74q7toPeNguheEJmp/gzTYlS5dWitXrnTf1A888IC2bdum9957T8uWLdP555+vndZeFqIIdqHh8GHpmWek0aOlI0ekkiWl+++XHnhAKlMm0HcnbU3c6ppq31v8npbuWJqhv96Vza/U9a2vV/eG3V3oA4CiwP4wGDt2bLbT51iwsVk3itIfDAnBMHiibNmy2r17t9ufNm2aLrzwQrdfsmRJHTp06JRvCsiJBTmb0NgWOene3Rv0nnxSat7cu+ZsoGciqRlTU/effb+WDFmixXcu1shzRqpe+XpuhO37i9/XRR9dpOovVNegrwbpv3/8V0eSjwT2hgGggDEnYuHIV7CzIHfbbbe58scff+iiiy5yx63GjuHuKExNm0rTp3sHV9hMC3/+6V2e7JxzpHnzguNr0aZ6Gz3d/WmtH75ev97yq4Z2GqrqZaq75czeXfiuLvn4ElV7vpoGThyor1Z+5Zp0ASDcMCdiEAe7V155xa0Va02un3/+uSpXruyOz58/X9fbzLJAIbKBE337eqdGsZGypUt715094wxpwABp48bg+HLYvI5n1TlL4y4apy33bdGPA390IS+2bKzik+L13qL3dMWnV6jK81V01YSr9MHiD1jpAkDYsH6Z/rwOWWO6k0zoYxf6tmzxTmbsGz1booQ0bJg0cqR0/G+QoJKckqxfN/+qz5d/ri9XfqnNCZtTzxWLLKbz6p2nS5te6gZeNKrUKKD3CgD5RR+7IB48YX7++We9/vrrbhLizz77TLVq1dL777+vBg0a6BxrBwtRBLvgcaoj8GwiYxtM8cMP3tc2552Fu7vvlkoF6cBU+3ZcsH2BvlzxpSaumphh4IVpUaVFasjrUqeLC37BihGUADJjVGyQBjtrfr3pppvUv39/F+aWL1+uhg0baty4cfrmm29cCVUEu/Ca58je3VOnegOeb+nimjWlhx+Wbr218Jcny6vVu1dr8h+TNemPSfp50886lnIs9VyFkhV0YcML1adxH/Vq3MsN2AgWzFMFgJ8PIRTs2rdvr3vvvVcDBgxwNSmLFi1ywW7BggXq06ePtm/frlBFsAvPv+hsSTIbLfvIIzYyy3usTh1vwBs0yNtcG+xssMXUNVNdyPt2zbdu9Yv02lVvp96Ne6tno55ucuToYoFJrfxFDiAn1OgH4Tx2VktnI2DTBztrlrXVJw7b3BMhimAX3n0wkpKkf/9bevppaetW7zFbKOXRR70DLYoXV0iwfnlzt87Vt6u/dSHPVr/wKO1buVSxUupar6sLeVar17paazd4o6DRhwYAQnAeu9jYWK1Zs+aE47/88osLeECwznNkTa+2zuzatd7lyGJjvaNmb7tNatLERnxLoTAVY1RklM6sfaaeuOAJ/Xb7b4r7a5w+uPIDt0atjbI9dOyQpq6dqhHTRqjta20V+89YXfef6/Tm/De1ds9a15evIDBPFQAEVr6C3e233+5qTObMmeNqAbZu3aoPP/xQI0aM0BBb6wkI8nmObIJjG0RhAe+f/5SqVfMGvL/8xTsf3rPPWoBUyKhapqr6t+2v8VeM19b7trqJkV/s+aJrmi1dvLR2HNjh1rIdPHmwGr/cWPXH1tfNE2/WOwve0bq96/wW9JinCgACK19NsfYhTz/9tEaPHq2DBw+6Y9HR0br//vs1cuRIlQrWIYe5QFNsYG3YsEHjx4/P8bqBAwf6dTJsq6V7+23p+efT5r2zUbRWu2dTpVjNXqiyVS3m/DlH36//Xt+t/06z/5ytoylHM1xTu1xtN62KlXPrnatmlZvlq+k2UF8/AAhVCcHQx87nyJEjrkl2//79rm+dTX/y/PPPM3gCIdtH6+hR6eOPvWvQrlzpPWYDK2w1i/vuk1q3Vsg7cOSAmzfvxw0/asbGGZq7Ze4JQa9K6So6p+45OqfOOW7boUYHFY8qHvRfPwAINQENdklJSRo1apSmT5+eWkN3xRVX6J133tEjjzyiqKgoDR06VA/Y3BIhihq7wAuGUZUpKdJXX0kvvOBdxcKnZ09pxAhbVs+74kU4OHj0oGZtnuVCnpXftvymw8cyDoCywRin1zzdrZxhpUvtLq75N1i/fkjD6EMguAU02Flgs1q5Hj16aObMmW5JsUGDBmn27Nl66KGHdM0117hwF8oIdsEhmOZBmz3b2w/viy+8gc/YLVh/vJtusuVvFFaSjiXp922/65dNv+iXzb+4beapVUzjSo3dAI7OtTq70i62nUpElQi6r19RxtcBCH4BDXY24nXMmDG67LLLtHTpUrVt21Y333yz3nrrrUKZSqEwEOyCR7DVNKxf7x1J+9Zb0v793mMW6m6+WbrrLql5c4WlFE+K/tj9h2Zunulq9mb+OVPLdy4/4broqGi1r9FeZ9Q8Q51qdVLH2I6KPhCtA/sPBMXXr6ih5hQIDQENdiVKlND69evd8mHGBkn89ttvatOmjcIFwQ45v0ek996Txo2TVq1KO969u3THHdLll4fGhMenYu+hvZqzZY4blOG2W+ZkWatXPrq8OtbsqE41O6ljjY5uv0GFBmHzh2Cwoq8jEDoCGuysmdVWlaha1du3xv4KX7x4sVsfNlwQ7JBb9p3z3XfegDdpUlozrX172GoWt98uNW5cNJ6n/RhZs2eNC3g2GMMmT7Y1bzP31TMVS1Z0gzF8pX1sezWp3ESREdTm+Qujk4HQEdBgZ80otmSYDZwwkyZNUrdu3VSmTJkM131hnZFCFMEO+bFhg7eJ1sq2bWnHu3WTbrlFuvJKW7GlaD3bo8lHtWznMhf0bGWM+dvma8mOJW76lczKFC/j+uidVv0015xry6O1qtbKzcGHvFuyZEmufg737ds3rFpcgFAU0GBnAyVyw0bJhiqCHU7FsWPSf/8rvfGG9O233lo9Y9+r117rrck788zwGVGbVxbqlu5Yqvlb57ugt3D7Qi2OW+xWysjMavCaVGriAp8FvbbV26pNtTaqW74uTbk5oMYOCB1BNY9dOCLYwV9s1TP7G+fdd701ej7NmtkEvdINN3jXqS3qjqUc0+rdq13T7YJtC7QwbqEWbV+knQd3Znl9uehybu1bC3lWrGbPXtvce/Cijx0QOgh2IfaAAet799NP3oD32WfS8cVanHPP9U5+fM01UqVKPKv0tu/f7gLeojhvWRK3RCt3rTxhMmWfamWqqVXVVt5SrZVaVm3pSlENfIyKBUJDkQ12NsXKwoULtWPHDlWsWNHNpffss8+qZs2aqdfYQA6bIHnu3LlugMewYcP0t7/9LU//H4IdCpItcfuf/0gffCD98ENaU23x4lKfPtL110uXXCKVLcvX4WRNuTb1ioU8669nZdmOZVq/b/1JH1jV0lXVomoLtazSUs2rNHf7trVl1MJ9wAbz2AHBr8gGu5deekldunRRjRo1tGXLFv31r391x22iZN+Dadq0qQt8tl6tdR6+5ZZb3Lx7gwcPzvX/h2CHwrJli3f5sg8/lBYuTDtuSy1ffLGt0ODdFrVBF/ldJm3FrhWu/54FveW7lmvFzhXZBj4bmGFr4jar0sxtm1ZumrqNiQ6fWaeDbT5IABkV2WCX2ddff+2WM7NlzooXL65XX31VDz/8sJuOxebbMw8++KAmTpyolb5FP3OBYIdAWLbMG/I+/VRasybtuIU6q8G76ipvjV64rXJRGIFv1e5VbkJlC3ord690zbnWp+9kTbqmRtkaLuDZ4A23rdzE7Teq1Egli5Us1H8DgPCWQLCT9uzZoyFDhriau19++cU9mAEDBriHY0HO54cffnDTsdj11nybFQuGVtI/4Dp16tDHDgFhf2ZZ7Z0FPFtu1Va78LFZhnr1sikqpEsvpU/eqQ7YWLd3nQt5q3atcs27FgCt7Diw46QfF6EI14RrQa9xxcZuWTULe40qNnLbsiVoQweQN0U62NlatePGjdPBgwd15plnavLkyapcubI717NnTzdRsq1l67N8+XK1atXKbU+2PuWoUaP0xBNPnHCcwRMINPvOnDdP+vxzb0lfk1esmNS1q/U99Ya8hg0DeafhZd/hfS7oWa3e6j2rvfvHtwlJaWvfnmwAhy/kNazQUA0qNlDDig1dqRlTM+z79AEo4sHOmkptAEROnX+bH1+Ec9euXa72bePGjS6M2YOwcGfLE+U32FFjh1Bg36XWXGsBz+adXbw44/nWrdNCXqdOtkpMoO40fNmPyl0Hd7kVNlLLXu927Z612n1od7YfXyKqhOqVr+fCni2rZqV+hfrutW1tkAdLrQFFT0I4BbudO3dq9+7sfxg2bNgwtc9cen/++adrMrXBEzaoIr9NsZnRxw6hYO1a7zJmX3/tnUolOTntXJUqUu/e0kUXeZtumUalcMQfjtfavWtdE68FPdva4A3bbozf6Jp/s1OqWCnVq1DPhTwLgK5UqOcmZLb9GjE1VCyyWCH9awAUlrAKdqfCRnnVq1fPhbfzzz8/dfBEXFycG0xhHnroIbesDoMnEM727vWucvHVV9LUqdaNIO2cDX7s0sUb9CzkdezoPYbCZaHuz4Q/tX7vem3Yt8EFPituf+96bU3cKo+y/1EcFRHl+vdZ0KtTvo7qlju+tdfl6rh9W4eXWj8gtBTJYDdnzhw3N90555zjat7Wrl2rRx991IW4ZcuWubVr7YE0a9bMNclaX7ylS5e66U5smhSmO0FRcfSoNGuW9M033qXNli7NeN66pF54oTfk9ewppZsGEgGen29z/GZXs2dhz4rtb4rfpI37NmpzwuYca/x8U7hYyLMAaEGvdkxtt+8rtcrVUuVSlQl/QBApksHO5qQbPny4Fi1apAMHDri57Hr37q1HHnlEtWrVynKC4ipVqrgJii3k5QVNsQi3Zc2sNs9q8r77zt7fGc9b19MePaTu3aXzz5fKlw/UnSI7ySnJbiUOC3sWAC3wWdizrW/f+v/lRnRUtAt4tWJquW3NsjVTX9sAD18pVbwUXxSgEBTJYFeYCHYI59q8OXO8Ic+KjbhN/91vTbQ28OKCC7zl7LOlMmUCecfIi8PHDrvmXgt+FvRsuyVxizvm22Y3lUtmFUpWSA15Nq9fhm1MDcWWjXWvy5TgTQKcCoJdASPYoajYs0f68Ufpf//zltWrM563rqpnnOGtybNiffUIeqEt6ViStu3fpi0JW1zYS90mbnH9/KzYsUPHDuX6c8aUiPGGvONhL7ZMrHd7vFQvW91tbdRv8Shv/2cAaQh2BYxgh6LcbGvNtRb2bB3bzZsznre582zwhc2fZ8Vq9HI52BwhxBpx4pPiXcCzEGhhb1vi8e3x19YsbPsHjx7M0+e2/n2+sGdz/lUvU90Vt3/8mK9Yf0GgKEigKTa0HjAQiqyJ1la9sIBnxaZUyRz0IiKkVq2kc87xhjwr9et7j6NoBMD9R/a7gGfBz8JeajmQth+3P841ASd70s3JkwtlipdxYc9q+izopW7LVHX7mbf0CUSoSiDYhdYDBsLFxo3egOcrf/xx4jU1akhnneVttrXSoYNUkqVVi7wUT4p2H9ytuANxGcKevfZtfcesJCWnLfOYlyDoC3lVSldJLb7XlUtXznC8UqlKzAuIoECwC7EHDISruDhp5kzJlmv+9Vfp99+9AzQy99OzcHfmmVLnzt4+e7b8GbV6yK4mMPFIojfw7Y/TzoM73f7OA8e3vtcHd7pjNhr4aEqmN14eBohY87Av9Ln9468zby0I2r41ETNXIPyJYFfACHZA/hw6JM2d651Hz1d2ZDEI0+bSs4BnxUbhWqlWjaeO/AdBW8PXgp6FPF/Yc/sW/g7udLWFvmNW9h7em+/HbUvD+UKebX3FJodO3S9VMcMxe10+uryiIlnrDyci2BUwgh3g3356FvBmz5Z++01auFA6cuTEa+vUkU4/Pa1YLZ8tjQYUBJvsee+hvW59Xwt6Fvwy77tyfH/PoT1uP781gyZCESoXXS5D6LP9CtEV3Nb32rZWk+gr7pqSFVygRHhKoI9daD1gAGmSkqRFi7whz+bUs7n0Vq3KOJ+eT9263oBnI3Ft2769FBtLMy4CVzN44OiB1JBngc/CodX+2THfcXvtO2bnbWsfd6psLeH0gS99sdpAty1Z/oR937ZsibI0IQcpgl2IPWAAOX3PSQsWeEOer6xZk/W11mR72mnekGdbK02aSFG0cCHIl4zbd3hfathz4S/d1s75AqG9tulmfMdt3x8iIyJdjWH6sOfbpj/u27et2z9+zFeoOfQ/gl0BI9gBgRcf7222tQEZVubP99bspaSceK2Num3dWmrb1lvatZPatPH25QPCYTk560NoAc+CXvriC4G+AOi2h+NP2M/NOsO5ZUvS+UJeTHRMhtBnk1Vn2EbHnLCffktI9CLYFTCCHRCcDh6Uli711u5Z6LOyeLH3eFZs6hULeBb6bGvF1sYtzby3KGJNyDaRtAU8X9BLv7XQ6AuO6bdu/3DasbxORp3bkOgLetZUnDn8uWPpzmX32krJYiVDsrmZYBdiDxhAwUlOltat8wY867vn227YkPX19jO/QQPvxMq+0rKl1Lw5gQ/IjtX6JSYluqBn09H4gp9v33fOV+x4+nPu9fFr8jNPYW6bm8seD3m+YvMbpu6XKKOyxdPtZ7rGjtm+75xv3/o3FmRgJNgVMIIdEPoSE6Vly6QlS9KK1fbt2pX19fYzu149b8izWj1fscBXqVJh3z0Q3o4mH00Neum3tpLJSfePX2cDUWxrx3zHC6I2MfOIZpu/MH3wy3J7fN9dm2nf9/G2zXws6UCSKlSo4LcKpQiP1dMiFcEOCF82r54FvvRl+XJp9+6Tf0zVqt6A16xZxmI1fzYBM4DAr2xi4c4X9vYfD4UWArM6duDI8eNH96fu+6611779w8cOF8r9RyRFyDPaQ7ArKAQ7oOjZuVNascIb8nxbG6yReX3c9IoV866iYSHPRuY2berdWqlVS4qMLMx/AYCCGLhy8OjB1DCYm6273nfs+GvfMd/n8h2z0dKO5cdnbNAYNXYFgmAHwGf/fu+auBbyLPD59m17skEbplQpqVEjqXHjE0vt2kzPAkCu36IFvO27tqtZ7WYEu4JCsAOQE5t2ZcuWtJC3erW32L6ttnEsm9klrPm2fn1v8PMVq/mzYs27Zcvy/IGiJIGVJ0LrAQMoWo4e9Y7KXbvWO9GyFV/ws+N2Pjs2CbMFPAt6FgBt31dsNQ769QHhJYFgF1oPGADST8/y55/e0Je+WC2fTduyN4e16a3fnvXfsxG8Fvqs2L6vWPCzCZsBhA6CXYg9YADIrX370kKebdMXq+07nItBelbjlz7oWalTJ23fRvmG4ByuQNhKoMYutB4wAPiDTUwVF+cNeBs3ere+Yq+tZDegwyc62juAw8Kelcz7Vmw5NsIfUDgIdiH2gAGgsILfnj3egLdpk3dr07XYvq9s3+69Lrfhz4o1/WZVbMk2+vsBwZc7ivnhngAAAWY1bFbTZqVDh6yvOXLEO5rXAl/mYset/59N4pyUlNb/L7v/nzXrWsirWTOtWOBLv7WmYZvzD0Dh4NsNAIqIEiXSRtiejIW6rVu9Ic8X9mybvth5G91rIdDKggXZB0ALdxb00pfY2LTie12mTIH8s4EihWAHAMjQDJtT+LN5/GwZNgt4vrC3bZv3tW9rxZp+7VrrG2hl4cLsH7TN4Ve9esbQZ6+zKqVL80UDskKwAwDkiU27Ys2wVtq1y356F1uuzcJe+mKBL32xYzbww1b6sJJdE7CP1e5ZTaCFPNuerNg9VqlCczCKjgiPJzddaYsOBk8AQOGy30KJiWk1e+lDn+9Y+pKbaV8yq1QpLYymLxb6Mu/b1paFAwoDgycAAGHF+uHZYEArTZrkLgT6+vdZ0PNtrXbQd9y3v2tX2ohhK7YMXG5YU2/6oGeDUmybft83WMVXaB5GMKApFgAQkiGwceOcr7fmYAt0vrCXPvDZvm/rK9Z30AaGWNOwb5qY3LJVPzKHPStWW+jb+orvdcWK3n6NgL8Q7AAAYSsqKq2pNTd8NYLpQ5+FPdv3bX3FXvvKsWPeJmLfYJK8sP6CFvB8oc/2fa+z2veVChW8/z4gPYIdAABZ1Ag2apS3MJg+6FmxmsL0++mLHbMl5GzU8IED3mJTy+SV3acFvPRhz7f1lcyvrZQvL8XEsMJIOCLYAQDgpzCY3TQxmVmoS0g4MfTt3Zu2zbzvKzZ62NjHW8lLk3H60c0W8Kz4wp5vm9W+/ft8r33FBpmw/FxwIdgBABAAFqx8NWgNG+btY60foAU8q/Wz4tvP6phvPz4+bd8+3oKlLyjml60q4gt8udn6itUWpn9t/RMJiP5BsAMAIMTYOr2+ufryypqOrT9g+rDn2/r2cypWS2jB0PoW+moaT4X1Fcwc+rLa5qaUKeMNzUUVwQ4AgCLEasasCdWKreebHxYOrTnYF/LSb3376UtWx6xY30Tf6OVTrT3MvIpJ2bJpYc+3n9U2837mYkExlJqcCXYAACBPLOT4QtOpsFo/C4gW8DIHPt+xrPZPVlJSvJ/Xt4qJTXLtr39v5rCXeT+rbXb7Vgpi7kOCHQCg0KWkpGjTpk1KTExUTEyM6tatq8ii3H5WRNmX3NcEW6vWqX0uq0U8dCgt5PkCY/p9X+DLfN5GJfvO+Y7bMZvP0Pe5fdf6m7/nMSTYAQAK1YoVKzRlyhS3lJJPuXLl1Lt3b7Vo0YKvBvJdq1a6tLfYGsL+YE3EvnWMLeilD4Hpw2DmY+m36Uv6Y74FXZOS/PsFZ63YTFgrFgAKNtRNmDDhpOf79etHuEPY8xwfwGIBb/v2BLVpU17x8fHuD5xTRb03AKDQml+tpi47dt6uA4rCAJYqVaS6df37uQl2AIBCYX3q0je/ZsXO23UA8odgBwAoFDZQwp/XATgRwQ4AUChs9Ks/rwNwIoIdAKBQ2JQmOXUOt/N2HYD8IdgBAAqFzVNnU5pkx84znx2QfwQ7AEChsXnqbEqTzDV39pqpToBTxwTFAIBCD3fNmjVj5QmgABDsAACFzppb69evz5MH/IymWAAAgDBBsAMAAAgTBDsAAIAwQbADAAAIEwQ7AACAMEGwAwAACBMEOwAAgDBBsAMAAAgTBDsAAIAwQbADAAAIEwQ7AACAMEGwAwAACBMEOwAAgDBBsAMAAAgTBDsAAIAwQbADAAAIEwQ7AACAMEGwAwAACBMEOwAAgDBBsAMAAAgTxQJ9AwCANCkpKdq0aZMSExMVExOjunXrKjKSv8EBhGmwS0pKUufOnbVo0SItWLBAp512Wuq5xYsXa+jQoZo7d66qVq2qYcOG6W9/+1tA7xcAcmvFihWaMmWKEhISUo+VK1dOvXv3VosWLXiQAHIUcn8GWlCrWbPmCcftB2HPnj1Vr149zZ8/X88//7xGjRqlN954IyD3CQB5DXUTJkzIEOqMvbbjdh4AwirYffvtt5o2bZpeeOGFE859+OGHOnLkiN5++221atVK1113ne6++269+OKLAblXAMhL86vV1GXHztt1ABAWwS4uLk6333673n//fZUuXfqE87NmzVLXrl1VokSJ1GO9evXSqlWrtHfv3mybdu0v4vQFAAqT9anL6WePnbfrACDkg53H49HNN9+sO++8U6effnqW12zfvl3Vq1fPcMz32s6dzOjRo1W+fPnUUqdOHT/fPQBkzwZK+PM6AEVXQIPdgw8+qIiIiGzLypUr9fLLL7sfaCNHjvT7PdjnjI+PTy2bN2/2+/8DALJjo1/9eR2Aoiugo2JHjBjhauKy07BhQ33//feuqTU6OjrDOau969+/v8aPH6/Y2FjXXJue77WdOxn7nJk/LwAUJpvSxEa/Ztcca+ftOgAI2mBnU5JYycm//vUvPfXUU6mvt27d6vrPffrpp27qE9OlSxc9/PDDOnr0qIoXL+6OTZ8+Xc2aNVPFihUL8F8BAKfG5qmzKU1s9OvJ2HnmswMQFn3s7K/U1q1bp5amTZu6440aNVLt2rXd/g033OAGTtx6661atmyZC31jx47VfffdF+C7B4Cc2Tx1/fr1czVz6dlrO848dgDCcoLik7GBDzYVik1Q3LFjR1WpUkWPPfaYBg8eHOhbA4BcsfBmrQysPAEgvyI8NuQUqayPi4VEG0iR+S9nAACAYM4dIdEUCwAAgJwR7AAAAMIEwQ4AACBMEOwAAADCBMEOAAAgTBDsAAAAwgTBDgAAIEwQ7AAAAMIEwQ4AACBMEOwAAADCBMEOAAAgTBDsAAAAwgTBDgAAIEwQ7AAAAMIEwQ4AACBMEOwAAADCBMEOAAAgTBDsAAAAwgTBDgAAIEwQ7AAAAMIEwQ4AACBMEOwAAADCBMEOAAAgTBDsAAAAwgTBDgAAIEwQ7AAAAMIEwQ4AACBMEOwAAADCBMEOAAAgTBQL9A0A/pKSkqJNmzYpMTFRMTExqlu3riIj+duFZwgARQfBDmFhxYoVmjJlihISElKPlStXTr1791aLFi0Cem+hgmcIAKGP6gyERSCZMGFChlBn7LUdt/PgGQJAUUCwQ8g3v1pNXXbsvF0HniEAhDuCHUKa9anLXFOXmZ2368AzBIBwR7BDSLOBEv68rijiGQJA+CDYIaTZ6Fd/XlcU8QwBIHwQ7BDSbEoTG/2aHTtv14FnCADhjmCHkGbz1NmUJtmx88xnxzMEgKKAYIeQZ/PU9evX74SaO3ttx5nHjmcIAEVFhMfj8QT6JoKJjaAsX7684uPjc2ziQ3Bh5QmeIQAU9dzByhMIG9bcWr9+/UDfRkjjGQJAaKMpFgAAIEwQ7AAAAMIEwQ4AACBMEOwAAADCBMEOAAAgTBDsAAAAwgTBDgAAIEwQ7AAAAMIEwQ4AACBMEOwAAADCBMEOAAAgTBDsAAAAwkSxQN8AAABFQUpKijZt2qTExETFxMSobt26ioykfgX+RbADAKCArVixQlOmTFFCQkLqsXLlyql3795q0aIFzx9+w58KAAAUcKibMGFChlBn7LUdt/OAvxDsAAAowOZXq6nLjp236wB/INgBAFBArE9d5pq6zOy8XQf4A8EOAIACYgMl/HkdkBOCHQAABcRGv/rzOiAnBDsAAAqITWlio1+zY+ftOsAfCHYAABQQm6fOpjTJjp1nPjv4C8EOAIACZPPU9evX74SaO3ttx5nHDv7EBMUAABQwC2/NmjVj5QkUOIIdAACFwJpb69evz7NGgaIpFgAAIEwQ7AAAAMIEwQ4AACBMEOwAAADCRMgEO+twGhERkaE888wzGa5ZvHixzj33XJUsWVJ16tTRc889F7D7BQAAKGwhNSr273//u26//fYsl2CxRZR79uypHj166LXXXtOSJUt0yy23qEKFCho8eHCA7hgAAKDwhFSwsyAXGxub5bkPP/xQR44c0dtvv60SJUqoVatWWrhwoV588UWCHQAAKBJCpinWWNNr5cqV1b59ez3//PM6duxY6rlZs2apa9euLtT59OrVS6tWrdLevXsDdMcAAACFJ2Rq7O6++2516NBBlSpV0syZMzVy5Eht27bN1ciZ7du3q0GDBhk+pnr16qnnKlasmOXnTUpKcsUnPj4+tWkXAACgIPnyhsfj8c8n9ATQAw88YP+KbMuKFSuy/Ni33nrLU6xYMc/hw4fd6wsvvNAzePDgDNcsW7bMfY7ly5ef9B4ef/zxHO+BwjPgPcB7gPcA7wHeA7wHVIDPYO3atX7JVhEev0XEvNu5c6d2796d7TUNGzbM0Lzqs2zZMrVu3VorV6506+8NGDDApd6JEyemXvPDDz+oW7du2rNnT65r7Pbt26d69eq59fzKly9/Sv++osKeu41C3rx58wmLXIPnxnst8Pge5ZnxXgte1lJYt25d123MBnyGdFNs1apVXckPGxhh6+5Vq1bNve7SpYsefvhhHT16VMWLF3fHpk+f7kLfyUKdiY6OdiUzC3WElLyx58UzyzueG8+ssPBe45nxXgtelmn88nkUAmxgxJgxY7Ro0SKtW7fOjYC99957deONN6aGthtuuMHV7N16662uNu/TTz/V2LFjdd999wX69gEAAApFSAyesBq1Tz75RKNGjXLNpjZIwoJd+tBmNWzTpk3T0KFD1bFjR1WpUkWPPfYYU50AAIAiIySCnY2GnT17do7XtW3bVj///PMph8jHH388y+ZZ8Mz8ifcaz6yw8F7jmfFeKzrfnwEdPAEAAAD/CYk+dgAAAMgZwQ4AACBMEOwAAADCBMEunVdeeUX169dXyZIl1blzZ/3222+B+8oEoZ9++kmXXnqpatasqYiIiAyTQRvrrmkjkWvUqKFSpUqpR48eWr16tYqy0aNHq1OnToqJiXFzLl5xxRVu/eL0Dh8+7EZz2zrIZcuW1VVXXaW4uDgVVa+++qobCOWbc83mqPz2229Tz/O8cr+2tn2f3nPPPTy7k7CZFuwZpS/NmzfneeXCli1b3JRj9nPLft63adNG8+bNSz3P74MTWb7I/H6zYj///fmzjWB3nM17Z9On2MiU33//Xe3atVOvXr20Y8eOPD/UcHXgwAH3XCwAZ+W5557Tv/71L7322muaM2eOypQp456hvVmLqhkzZrhvVBvVbRNm2wTaPXv2dM/Sx6bumTRpkj777DN3/datW9W3b18VVbVr13ahZP78+e4Xha0ec/nll7v5KQ3PK2dz587V66+/7gJyejy7E7Vq1cqtO+4rv/zyC88rB7ZCwtlnn+0WA7A/upYvX65//vOfGRYD4PdB1t+X6d9r9jvBXHPNNf79/vTLwmRh4IwzzvAMHTo09XVycrKnZs2antGjRwf0voKVvXW+/PLL1NcpKSme2NhYz/PPP596bN++fZ7o6GjPxx9/HKC7DD47duxwz27GjBmpz6h48eKezz77LPUaWx/Zrpk1a1YA7zS4VKxY0fPvf/+b55ULiYmJniZNmnimT5/uOe+88zzDhw93x3mvZb1WeLt27bJ8jjyv7Nd5P+ecc056nt8HuWPfm40aNXLPy5/vN2rsJB05csTVDljTYfqlPey1rXqBnK1fv17bt2/P8Axt0mhr0uYZZlwT0FSqVMlt7X1ntXjpn5s1Bdm6gTw3KTk52U1ObjWc1iTL88qZ1RBffPHFGd5TvNdOzrqLWPcSW5e8f//+bp1wnlf2vv76a51++umupsm6mLRv315vvvkmvw/ymDs++OAD3XLLLa451p8/2wh2knbt2uV+gVSvXj3Dw7HXFlaQM99z4hmeXEpKiuvvZE0YrVu3Tn1uthRe5oWfi/p7b8mSJa6PiU3Yeeedd+rLL79Uy5YteV45sBBsXUmsb2dmvNdOZH94vvvuu5oyZYrr22l/oJ577rlKTEzkeWXDlva059WkSRNNnTpVQ4YM0d13363x48envtcMvw9Ozvqo79u3TzfffLPfvz9DYuUJIFxqUpYuXZqhDw+y1qxZMy1cuNDVcP7nP//RwIEDXZ8TnNzmzZs1fPhw12/HBoAhZ3369Endt/6IFvTq1aunCRMmuAEBOPkfqVZj9/TTT7vXVmNnP9usf7V9ryJnb731lnv/WW2xv1FjJ7l1ZaOiok4YfWKvY2Nj/f7Qw5HvOfEMs/aXv/xFkydP1g8//OAGB6R/blYlb3+5pVfU33v2l2vjxo3dus9W+2SDdsaOHcvzyoY15dhgL1uCsVixYq5YGLYBTbZvf/nzXsue1ZY0bdpUa9as4b2WDZv5wGrQ02vRokVqMza/D7K3ceNG/e9//9Ntt91WIL8LCHbHf4nYL5Dvvvsuw18k9tr69SBnDRo0cG++9M8wISHBjY4tys/QxplYqLOmxO+//949p/TsfWcjy9I/N5sOxX5AFuXnlpl9PyYlJfG8stG9e3fXhG01nb5itSrWb8y3z3ste/v379fatWtdcOF78+SsO0nmaZv++OMPV9tp+H2QvXfeecf1TbS+sD5+fb/laahFGPvkk0/cCM53333Xs3z5cs/gwYM9FSpU8Gzfvj3QtxZUo+0WLFjgir11XnzxRbe/ceNGd/6ZZ55xz+yrr77yLF682HP55Zd7GjRo4Dl06JCnqBoyZIinfPnynh9//NGzbdu21HLw4MHUa+68805P3bp1Pd9//71n3rx5ni5durhSVD344INu1PD69evd+8heR0REeKZNm+bO87xyL/2oWJ7diUaMGOG+N+299uuvv3p69OjhqVKlihu9zvM6ud9++81TrFgxzz/+8Q/P6tWrPR9++KGndOnSng8++CD1Gn4fZM1m3LCf9zayODN//Wwj2KXz8ssvu4daokQJN/3J7Nmz8/xAw9kPP/zgAl3mMnDgQHfehmw/+uijnurVq7uQ3L17d8+qVas8RVlWz8vKO++8k3qNBd+77rrLTelhPxyvvPJKF/6KqltuucVTr149931YtWpV9z7yhTrD88p/sOPZZXTttdd6atSo4d5rtWrVcq/XrFnD88qFSZMmeVq3bu1+1jdv3tzzxhtvZDjP74OsTZ061f0OyOp3o7++PyPsP3mr4wMAAEAwoo8dAABAmCDYAQAAhAmCHQAAQJgg2AEAAIQJgh0AAECYINgBAACECYIdAABAmCDYAQAAhAmCHQAAQJgg2AHAcTfffLMiIiJcsQW5q1evrgsvvFBvv/22UlJSeE4Agh7BDgDS6d27t7Zt26YNGzbo22+/1QUXXKDhw4frkksu0bFjx3hWAIIawQ4A0omOjlZsbKxq1aqlDh066KGHHtJXX33lQt67777rrtm0aZMuv/xylS1bVuXKlVO/fv0UFxeX+jkWLVrkAmFMTIw737FjR82bNy/1/C+//KJzzz1XpUqVUp06dXT33XfrwIEDfB0AnDKCHQDkoFu3bmrXrp2++OIL1yRroW7Pnj2aMWOGpk+frnXr1unaa69Nvb5///6qXbu25s6dq/nz5+vBBx90Tbtm7dq1rlbwqquu0uLFi/Xpp5+6oPeXv/yFrwOAUxbh8Xg8p/5pACA8+tjt27dPEydOPOHcdddd54LY2LFj1adPH61fv97Vtpnly5erVatW+u2339SpUydXS/fyyy9r4MCBJ3ye2267TVFRUXr99ddTj1mwO++881ytXcmSJQv4XwkgnFFjBwC5YH8D26CKFStWuEDnC3WmZcuWqlChgjtn7rvvPhfgevTooWeeecbV0qVvprUmXWvG9ZVevXq5mkALiwBwKgh2AJALFtoaNGiQq2c1atQoLVu2TBdffLG+//57F/y+/PJLd27//v264447tHDhwtRiYW/16tVq1KgRXwsAp6TYqX04AIQ/C2dLlizRvffe6/rObd682ZX0TbHWhGsBzqdp06au2Mdcf/31euedd3TllVe6ARl2fePGjQP4LwIQrgh2AJBOUlKStm/fruTkZDfSdcqUKRo9erSb7mTAgAGKjIxUmzZt3ACJMWPGuClQ7rrrLtdH7vTTT9ehQ4d0//336+qrr3Y1fH/++acbRGGDJcwDDzygM8880w2WsObaMmXKuKBngzDGjRvH1wLAKSHYAUA6FuRq1KihYsWKqWLFim407L/+9S83EMJCnbHpT4YNG6auXbu6YzbK1QZLGBsYsXv3bhcCLRhWqVJFffv21RNPPOHOt23b1o2mffjhh92UJ9Z3z5pg04+qBYD8YlQsAABAmGDwBAAAQJgg2AEAAIQJgh0AAECYINgBAACECYIdAABAmCDYAQAAhAmCHQAAQJgg2AEAAIQJgh0AAECYINgBAACECYIdAABAmCDYAQAAKDz8P5Z7Y0ntDQ3nAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y2 = estr.theta[0] * x / (estr.theta[1] + x)\n",
    "\n",
    "plt.plot(d, r, 'o', color='gray')\n",
    "plt.plot(x, y, '-', color='blue')\n",
    "plt.plot(x, y2, '-', color='green')\n",
    "plt.xlabel(\"Dose\")\n",
    "plt.ylim([-50, 0])\n",
    "plt.ylabel(\"Response\")\n",
    "plt.xlim([0, 70])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93e100ae-962c-40f6-85e0-5c39d142fdc0",
   "metadata": {},
   "source": [
    "Here, we observe a minor change in the estimated coefficients. In the plot, we can see that the dose-response curve estimated using the outlier-robust approach (green) did not drop as low as the non-robust method (blue). This is expected given the apparent outliers around doses of 15, 25, and 52. \n",
    "\n",
    "## Chapter 9: Nonlinear Mixed Effects Models: Case Studies\n",
    "\n",
    "This chapter presents some case studies of nonlinear mixed effects models. Here, we will review some of the case studies but modify them to use estimating equations instead.\n",
    "\n",
    "### Pharmacodynamic Modeling of Acetylcholinesterase Inhibition\n",
    "\n",
    "Data comes from Cutler et al. 1995 and is provided in Table 1 of this chapter (pg 360). Data comes from 12 participants with Zifrosilone dose on inhibition of acetylcholinesterase. Here, multiple measurements were obtained from the same participants over different timings. Here, we will explore the use of E-max models again but with some caveats. First, we load the data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "94400bae-c29c-4c72-9d65-e3a4195b432e",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = pd.read_csv(\"data/cutler1995_pd.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baf055a2-347e-44c8-8e40-a2f8ba6bd63f",
   "metadata": {},
   "source": [
    "In this data, participants have multiple measures over time. We can easily observe this by plotting our data (which has many more points that 12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4b4894cd-d0f5-464c-ac4b-c39e6cf5ffdd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQfFJREFUeJzt3Ql8lOW1+PETWcJmQRQSEWRRBBFEgcrmAgpEiveqgBE/KlFBWwVlqwvVirRWXHrFjcUFQe31IlDBqiAgAioEBBSNiFSWKBQCKoKyI87/c55/J00mycw7M+/Mu/2+n88Y5t3mZZwkh+d5zjkZoVAoJAAAAPC845y+AQAAANiDwA4AAMAnCOwAAAB8gsAOAADAJwjsAAAAfILADgAAwCcI7AAAAHyCwA4AAMAnKjt9A27zyy+/yPbt2+X444+XjIwMp28HAAD4UCgUkp9++kkaNGggxx1n3zgbgV0EDeoaNWpk2xsMAABQka1bt0rDhg3FLgR2EXSkLvxG/+pXv7LtjQYAAAj78ccfzUBSOO6wC4FdhPD0qwZ1BHYAACCV7F72RfIEAACATxDYAQAA+ASBHQAAgE8Q2AEAAPgEgR0AAIBPENgBAAD4BOVOAACOdvv55ptvTAV+red16qmn2lqFHwgaAjsAgCPWr18v77zzjinUGqb1Qy+99FI588wz+b8CJIB/FgEAHAnqZsyYUSqoU/pct+t+APEjsAMApH36VUfqotH9ehyA+BDYAQDSStfURY7URdL9ehyA+BDYAQDSShMl7DwOwH8Q2AEA0kqzX+08DsB/ENgBANJKS5po9ms0ul+PAxAfAjsAQFppnTotaRKN7qeeHRA/AjsAQNppnbrc3NwyI3f6XLdTxw7weYHiJk2ayNdff11m+2233SYTJkyQQ4cOyahRo2T69Oly+PBhycnJkYkTJ0pWVpYj9wsAiE6DtxYtWtB5AghiYLdq1So5duxY8fPPP/9cevbsKVdddZV5PmLECHn77bdl5syZUrt2bRk6dKj07dtXli1b5uBdAwCi0elW/Yc7AHtkhEKhkHjQ8OHD5a233pKvvvrK1DuqV6+evPrqq9K/f3+z/8svvzT/GszPz5dOnTpZvq5eSwPDvXv3xlzcCwAAkIhUxRueXGN35MgR+dvf/iY33XSTZGRkyJo1a+To0aPSo0eP4mNatmxpMqo0sAMAAAgCz0zFljRnzhzZs2eP3HDDDeZ5UVGRVK1aVerUqVPqOF1fp/ui0fV4+giLVQ0dAADArTw5YjdlyhTp3bu3NGjQIOlrjRs3zgyFhh+NGjWy5R4BAADSzXOBnWbGvvvuuzJ48ODibdnZ2WZ6VkfxStq5c6fZF83o0aPN/Hb4sXXr1pTdOwAAQCp5LrCbOnWq1K9fX/r06VO8rX379lKlShVZtGhR8bYNGzaYFPrOnTtHvV5mZqZZtFjyAQAA4EWeWmP3yy+/mMAuLy9PKlf+z63rFOqgQYNk5MiRUrduXROc3X777SaoiycjFgAAwMs8FdjpFKyOwmk2bKTx48ebekj9+vUrVaAYAAA/0sEO/Z34008/yfHHH28qQdCGDZ6tY5cq1LEDALjd+vXr5Z133ilVyUFnq7THLu3YvIE6dgAAwAR1M2bMKFOeS5/rdt2P4PJc8gQAAEGeftWRumh0vx6HYCKwAwDAI3RNXaxC+rpfj0MwEdgBAOARmihh53HwHwI7AAA8QrNf7TwO/kNgBwCAR2hJk1iF9HW/HodgIrADAMAjtE6dljSJRvdTzy64COwAAPAQrVOXm5tbZuROn+t26tgFm6c6TwAAgP8f3LVo0YLOEyiDwA4AAA/S6dYmTZo4fRtwGaZiAQAAfILADgAAwCcI7AAAAHyCNXYA4EPaK1TbSmkHAi1Wq3XNKIEB+B+BHQD4zPr1600j+JI9RbUUhtY3oxQG4G9MxQKAz4K6GTNmlGkUr891u+4H4F8EdgDgo+lXHamLRvfrcQD8icAOAHxC19RFjtRF0v16HAB/Yo0dkAQWqMNNNFHCzuMAeA+BHZAgFqjDbTT71c7jAHgPU7FAAligDjfSkiaRjeEj6X49DoA/EdgBcWKBOtxK69RpSZNodD/17AD/IrBDIAOzwsJCKSgoMF/jzRBkgTrcTOvU5ebmlhm50+e6nTp27vuZAtiJNXYIFDvWxbFAHW6nn+UWLVrQeSINWGsLt2HEDoFh17o4FqjDC3S6tUmTJtKmTRvzlelX+7HWFm5EYIdAsHNdHAvUAbDWFm5FYIdAsHNdHAvUAbDWFm5FYIdAsHtdHAvUgWBjrS3ciuQJBEIq1sWxQB0ILtbawq0I7BAI4XVx0aZjEyncGl6gDiBYUvUzBUgWU7EIBNbFAeBnCoKAwA6Bwbo4APxMgd9lhEKhkNM34SY6rF67dm3Zu3dvzJ6L8G6ZAs1o08XPuk5Gp0rirfFlxzX8gPcB4PsA7oo3WGOHwEl2XRyV5nkfADt/pgB2Ct4QA5AEKs3zPgCAmxHYARZRaZ73AQDcjsAOsIhK87wPAOB2ngrs/vWvf8l1110nJ554olSvXt00t169enXxfs0Duf/+++Xkk082+3v06CFfffWVo/cM/3C60ryOGBYWFkpBQYH5Gu5rW9F2v74PAAAfJE/88MMP0rVrV+nevbvMmzdP6tWrZ4K2E044ofiYRx99VJ566il56aWXpGnTpvLHP/5RcnJy5IsvvpBq1ao5ev/wPicrzVeUsNG6dWv5/PPPy2y/9NJLTXmXVKDiPgC4l2cCu0ceeUQaNWokU6dOLd6mwVvJ0bonnnhC7rvvPrn88svNtpdfflmysrJkzpw5MmDAAEfuG/7hVKX5cMJGJL2P5cuXl7tdj8/NzU1JcEfFfQBwL89Mxf7jH/+QDh06yFVXXSX169eXc889V55//vni/Vu2bJGioiIz/Rqm9WE6duwo+fn5Dt01/MSJ7hVWEjYqouelYlqWLh4A4F6eCew2b94skyZNkubNm8v8+fPl1ltvlTvuuMNMuyoN6pSO0JWkz8P7ynP48GEzwlHyAbile4WVhI2K6Hl6firQxQMA3MkzU7E68qAjdg899JB5riN2urZo8uTJkpeXl/B1x40bJ2PHjrXxTuF3GtS0aNEiLZ0nkk1ASGUCQzrfBwCANZ75CayZrq1atSrziyU8IpGdnW2+7ty5s9Qx+jy8rzyjR4827TzCj61bt6bk/uHPSvOama1fUxXMJJuIkYpEDifeBwCANZ75KawZsRs2bCi17Z///Kc0bty4OJFCA7hFixaVmopauXKldO7cucLrZmZmmmm0kg/ALcKJColIRSIHAMDdPBPYjRgxQlasWGGmYjdu3CivvvqqPPfcczJkyBCzPyMjQ4YPHy4PPvigSbTQml4DBw6UBg0ayBVXXOH07QMpS1RIVyIHAMD9MkJaJ8Qj3nrrLTN1qvXrdIRu5MiRcvPNNxfv17/KmDFjTMC3Z88eOf/882XixIlyxhlnWH4NHeXTbFqdlmX0Dm7hpjp2AIDkpSre8FRglw4EdnBzAlF5iQoVbQcABC/e8ExWLBB04UQFq9sBAMHDP+sBAAB8gsAOAADAJ5iKhe+w5gwAEFQEdghE9ihZogCAIGAqFr4K6mbMmFGmt6o+1+26HwAAPyOwg2+mX3WkLhrdr8cBAOBXBHbwBa3jFjlSF0n3h3sLAwDgR6yxg+uTHxo2bCjbtm2LWoBX91mxefNmCvgCAHyLwA6uT37QPsAlG6SUlwyhAZ8VH3zwgXz66ackUwAAfImpWLg++SGy6115yRA6ile9enVLr0MyBQDArwjs4JnkB7uTIUimAAD4DYEdPJP8EC0ZQr8ePHgw4fMBAPADAju4gtXkh4rOS/Z8AAD8gMAOrmA1+aGi85I9HwAAPyCwgyto8oNmu8ZDj9fz7DgfAAA/ILCDK2hdOi1hEg89PlzPLtnzAQDwA36rwTVatGgh3bp1K1O2ROvYRY605ebmlqpjp/S5bo8cubN6PgAAXkeBYri2MLEGeB07dpSuXbvG7DwRpsGaBojxdq4AAMAPCOzgmsLEkbR8yZIlS6R+/fpxja5p0NakSZNS2yKfAwDgRwxbwPWFiSkkDACANQR2cH1hYgoJAwBgDVOxSOvoXMm1b7rWzWqB4C+++MJ8ZX0cAAAVI7CDY8kRmp3arl07S+evWrXKPPQcLVNCRisAAGUxFYu0JUdETrnqc02OiCxvEo2eo9fSawIAgNII7OB4ckQiSKgAAKAsAjs4nhyhZU20MHE8LcFIqAAAoCzW2CGlCRKxgrqwunXryrBhw8y5miih6+lisZp4AQBAUBDYIaUJEjVq1LB0rgaBJQsLWwns9BwAAPAfTMUipQkSBw4ciHmuTsFqGZMw/XOsadnIcwAAAIEdXJAgoeVLSvZu1T/rtnjOAQAABHZIU4JEedOyOuqWm5tbbk063ab7Ikfuop0DAEDQscYOSSdJfPvtt5bOycnJMYFZyc4T0UbdNHhr3ry5rF69Wnbv3m0SLDp06CCVK/OxBQCgPPyGhC1JElZoUBdOjkj0dfLz8+k8AQBABVikBFuSJGKJN9khWrcKOk8AAFA+AjukJUkinmQHK69D5wkAAMoisIPtSRLJJjtYeR06TwAA4OHA7oEHHpCMjIxSj5YtWxbvP3TokAwZMkROPPFEqVWrlvTr10927tzp6D37jdVODxdccIH07dtX8vLyTDeJeDNYrb6OdqgoLCw0I3wAAMBjyRNnnXWWvPvuu8XPS2ZHjhgxQt5++22ZOXOm1K5dW4YOHWqCi2XLljl0t/5jtdNDs2bN4kqSSPR1tDuFPnRUUKd6KYECAAg6z4zYhQO57Ozs4sdJJ51ktu/du1emTJkijz/+uFx88cXSvn17mTp1qixfvlxWrFjh9G37Rro6Qlh5nZJIqAAAwIOB3VdffSUNGjQwI0LXXnutWYul1qxZI0ePHpUePXoUH6vTtBogaHkM2CNdHSGsvE55SKgAAASdZwK7jh07yrRp08wv70mTJsmWLVvMWi5dj1VUVCRVq1aVOnXqlDonKyvL7Ivm8OHDZsSn5ANl6To2Xc+mAfQ555wj1apVKzdJokWLFua4goKCpNa/VdR5IhoSKgAAQeeZNXa9e/cu/vPZZ59tAr3GjRubmmbVq1dP+Lrjxo2TsWPH2nSXwStIrO+9/r/QIHvDhg3y5JNPljoumfVveo4Gijoyq4kSup4uli+//DKp9X0AAHiZZ0bsIuno3BlnnCEbN2406+2OHDkie/bsKXWMZsXqvmhGjx5t1uiFH1u3bk3xnfurIPHBgwdlyZIlsmjRopQUFNZpWQ3UWrVqZen4lStXJvxaAAB4nWcDu3379smmTZvk5JNPNskSVapUMcFFmI4e6UhP586do14nMzPTjCqVfCD+gsSx1jImu/4tnoQK1toBAILKM4Hd73//e1m6dKlZt6XZrldeeaVUqlRJrrnmGlPeZNCgQTJy5EhZvHixSaa48cYbTVDXqVMnp289EAWJQ6FQSte/xZNQwVo7AEBQeWaN3bZt20wQ9/3330u9evXk/PPPN6VM9M9q/Pjx5pe/FibWhIicnByZOHGi07ftaVYLBafrerrmTgN1KyVs7L53AAC8wDOB3fTp06Pu1yzNCRMmmAfsYbVQcDqvp8kUVgI7u+8dAAAv8MxULNIv3kLB0Wj2bLKFi9NZJBkAAC8isIMt69q0jqCfiiQDAOBF/PZDUoWCdXu3bt1MuZlotCxKMskTVu4pXCSZnrEAgKDyzBo7OKdkoWDNOD1w4IDUqFGjeMpz3bp1lq5jZ0JDyXvS6+qaOr0XJ0bqtIyLG+4DAAACO8RVKDiZRAW7Exqi3ZOTXTmS6bYBAEAyGFZA0vbv3x/zGD8mNFTUlSPZbhsAACSKwA5JT0MuWLAg5nG9evXy1fSkla4cdMAAAKSbf37TwtXdKWrWrClB+3vTAQMAkG6ssUNSrCZEJJo44dbEhFT/vQEASASBHZKSysQJNycmOJUwAgBANM4PfcDTUtUJwu2JCXTAAAC4EYEdXNcJwguJCXTAAAC4EYEdXNcJwiuJCXTAAAC4DWvsYEvSggY5zZs3l9WrV8vu3bulbt260qFDB6lcuXLcr+OlxAQ3dcAAAIDADrYkLZR3Tn5+ftzn6Ou0a9fOU4kJbuiAAQCAYlgBSSct2H3OkiVLpHr16oHrZAEAQLII7JBU0kKqzokl3oQMAACCgN+MiDtpQUfUCgsLi9fHxZvoYOWcgwcPymmnnSa1atWyJSEDAIAgYI0d4k5G+OCDD8xDgyyrAVbJa1t9nU2bNklGRoacddZZJkGBxAQAAKJjxA4JJyPoqNvKlSvjvnY8rxMKhWTdunVSVFRkEhSYfgUAoGIEdoirm0J5dFQtnkSHRF5HM2x//vln/m8BABAFgR3i6iJR0ahaPIkOibyOvobWyAMAABUjsIOlbgqxdOrUKa7OE4m8jhY+ToQmeWiyR0FBQXHSBwAAfkTyBKJ2U9i8ebNJlIhFj+/Zs2fc3Sr0PC19smrVqpivod0s0lFsGQAAr2LEDlG7KXTr1i3mqFp4DV34nDZt2lhOdNBjevXqFXOdnu7XFmXxSKRwMgAAXkZgh+gfEAvr4ZItFqz9ZDt37hz1GN0fq+9sSYkUTgYAwOsI7BBTRevh7CwWrNO4Xbp0KTNyp891u+6PRyKFkwEA8DrW2CHudXdW19DFS4O37t27m+xXTZTQNXU6/RrPSF28RZCtHgcAgBcQ2MGy8Bq6VNIgTjNsk2W1CHK8RZkBAHAzpmLhS1aKIEcWTgYAwOsI7OBL6Uj6AADAbfitBt9KR9IHAABuwho7VEhLgaQyWSIdr5mOpA8AANyCwA6u6diQqtdMR9IHAABuwLAFXNGxgS4RAAAkj8AOjndsoEsEAAD2ILBzgAYyhYWFUlBQYL66qa2VEx0b6BIBt3Hz9ygA+HKN3cMPPyyjR4+WYcOGyRNPPGG2HTp0SEaNGiXTp0+Xw4cPS05OjkycOFGysrIkyGvX3N6xgS4RcBO3f48CgO9G7FatWiXPPvusnH322aW2jxgxQt58802ZOXOmLF26VLZv3y59+/YVt/DCOjInOjbQJQJu4YXvUQDwVWC3b98+ufbaa+X555+XE044oXj73r17ZcqUKfL444/LxRdfLO3bt5epU6fK8uXLZcWKFeI0r6wjc6JjA10i4AZe+R4FAF8FdkOGDJE+ffpIjx49Sm1fs2aNHD16tNT2li1bmqAhPz+/wuvplK3+a7zkIxW8so7MiY4NdImAG3jlexQAfBPY6dq5jz/+WMaNG1dmX1FRkVStWlXq1KlTaruur9N9FdFr1a5du/jRqFGjlNy7l9aROdGxgS4RcJqXvkcBwPPJE1u3bjWJEgsXLpRq1arZdl1NwBg5cmSpf5GnIrjz2joyJzo20CUCTvLa9ygAeDqw06nWXbt2Sbt27Yq3HTt2TN5//3155plnZP78+XLkyBHZs2dPqVG7nTt3SnZ2doXXzczMNI9UC68jizbVY/fatWQ50bGBLhFwihe/RwHAs1Oxl1xyiakptXbt2uJHhw4dTCJF+M9VqlSRRYsWFZ+zYcMGM+LUuXNncRrryAB343sUgB94ZsROpz9at25dalvNmjXlxBNPLN4+aNAgM61at25d8y/r22+/3QR1nTp1EjcIryOjRhbgTnyPAvA6zwR2VowfP978q7tfv36lChS7idfWkf3888+yevVq2b17twmYdWS0cmV7PjZaNiLW+2DlGCDI36MAUFJGKBQKldoScLq+RrNjtS5erHpufqeJKloqpuRHJCMjw4yC9uzZM+XV/ekAAADwqx9TFG8k/E/QTZs2yX333SfXXHONSWpQ8+bNk3Xr1tl2c3A2qNPizpFxvz7X7bo/ldX96QAAAECaAjtt19WmTRtZuXKlvP7666YbhPr0009lzJgxiVwSLpt+jVbUWel+PS4V1f31Hwj6iIYOAAAA2BTY3XPPPfLggw+aURstChymrbzc0L4LydE1dbFm6HW/HpeK6v66rilWEVg6AAAAYFNgp2VHrrzyyjLb69evL999910il4SLaKKEncelqmo/HQAAALAhsNMCwDt27Ciz/ZNPPpFTTjklkUvCRTT71c7jUlW1nw4AAADYENgNGDBA7r77btODVbMkdd3UsmXL5Pe//70MHDgwkUvCRbSkif5/jUb363GJVvePFbDFCtroAAAAgE2B3UMPPSQtW7Y0PVU1caJVq1Zy4YUXSpcuXUymLLxN69TF6tah+xOpZ2elun/v3r3NIxq9BnXFAACwsY7d1q1bzXo7De7OPfdcad68uXgddez+gzp2AAB4K96wpUDxsWPHTIDXuHFjOeGEE8TLCOxS33ki3E1C3+sDBw5IjRo1iqdW6TwBAAiCH1MU2CX0G3r48OGmjp32ZtWg7qKLLjJFa/UX9FtvvSXdunWz7QbhLA3i7Oy1G62bRHlTq7qtSZMmtr0+AAB+ltAau1mzZknbtm3Nn998803ZvHmzfPnllzJixAi599577b5H+ATdJAAAcGFgp7XqsrOzzZ/nzp0rubm5csYZZ8hNN91kpmSBRDpO0E0CAAAHArusrCz54osvzDSs/jION4TX9VKVKlVK8pbgZPBVWFhognP9qs/tYqXjBN0kAABwYI3djTfeaEbpTj75ZFPPrEePHma79o7VMijwnmhr384888y0dYmgmwQAAGkO7B544AFp3bq1KXdy1VVXSWZmptmuo3XaRxbeXPsWSYM83a5BfLLBndUuEXSTAAAgcQnXrejfv3+ZbXl5eUncCty89q1FixZJFQQOd5yINh1LNwkAABwK7BYtWmQeu3btKrMW68UXX0zytpAu8ax9S6bsSLjjRHkjg2F0kwAAIDkJDcGMHTtWevXqZQI7zZD94YcfSj3gHelc+6bTuTqtG1mIUZ/bMd0btIQUAABsGbGbPHmyTJs2Ta6//vpEToeLpHvtmwZvOq2rI4AaLOp1y+s44QepTkgBACBSQr9Njxw5Il26dEnkVLhMeO1bNHavfQt3k9DuJfrVr0GdTjtHTnOHE1J0PwAAdkvoN+rgwYPl1Vdftf1mkH7htW/RsPYtPhRjBgB4air20KFD8txzz8m7774rZ599tlSpUqXU/scff9yu+0MahNe+MW3orYQUAABsCew+++wzOeecc8yfP//881L7tGAxnBspSnTtWjrWviVzf156fYoxAwA8FdgtXrzY/juB4wv1w2vf/JhIkM7XpxgzAMApSQ9XbNu2zTzgHLcv1Hf6/tL9+k4kpAAAkHBgp1Naf/rTn6R27drSuHFj86hTp478+c9/pk5Xmrl9ob7T9+fE65OQAgDwVGB37733yjPPPCMPP/ywfPLJJ+bx0EMPydNPPy1//OMf7b9L2LJQP4j359TrB7EYMwDAo2vsXnrpJXnhhRfkv//7v4u3aXbsKaecIrfddpv85S9/sfMe4eGF+k7fn5OvH6RizAAADwd2u3fvlpYtW5bZrtt0H9LH7Qv1nb4/p18/lQkpAABESmjooG3btmYqNpJu031IH7cv1Hf6/px+fQAAXB/YPfroo/Liiy9Kq1atZNCgQeahf9b+sY899pj9dwnPLtR3+v6cfn0AANIpIxQKhRI5cfv27TJhwgT58ssvi9cT6fq6Bg0aiJfpQnrN9t27d2/MkR43cbpOnNvvz+nXBwAgHfFGwoGdX3k1sEtVZwU7r1nRtdLVEcLpzhcAAKQ63kgoeUL98MMPMmXKlOLirjoVe+ONN0rdunVtuzk4u1Df7lGu8u4vnSNpJDIAAPwuoeGK999/3/yCfuqpp0yApw/9c9OmTc0+eF86ujU43ZECAAC/SSiwGzJkiFx99dWyZcsWef31181j8+bNMmDAALMP3paObg1Od6QAAMCPEgrsNm7cKKNGjZJKlSoVb9M/jxw50uyD/TTAKSwslIKCAvM1MuD5+eefZcWKFTJ37lzzVZ+7uVuD0x0pAADwo4TW2LVr185Mk2lV/ZJ0W6rq2E2aNMk8NKhRZ511ltx///3Su3dv8/zQoUMm2Jw+fbocPnxYcnJyZOLEiZKVlSVeF2sd2sKFCyU/P19K5sEsWLBAOnfuLD179nRltwanO1IAAOBHCQV2d9xxhwwbNsyMznXq1Mls01EiLX+i/WM/++yzUq3G7NCwYUNz7ebNm5sARtuaXX755aZPrQZ5I0aMkLfffltmzpxpskyGDh0qffv2lWXLlomXhdehRQqvQ9PgesOGDWX263u0fPly8+d4g7t0dGtwuiMEAAB+lFC5k1glIjIyMkxgoV+PHTsmqaIZuFoQuX///lKvXj159dVXzZ+V1tfT0SwdyQoHn14rd6LTrU8++WTMKcto9P/BH/7wB6lcubKtr6vvjQb3yZQ+SfVrAADgVq4qd6JJE07SYFFH5vbv32+mG9esWSNHjx6VHj16lOpbq3XKYgV2Om2rj7Bkgii7WVmHFosG2KtXr44ruA13ayhvpNCubg3peA0AAIImocCucePG4gRNHNBATtfT1apVS2bPnm3q561du1aqVq0qderUKXW8rq8rKiqKes1x48bJ2LFjxY3sWl/2/fffx32Ojnbm5uamtMZcOl4DAIAgSSiw0/VtJ510kvTp08c8v+uuu+S5554zQdb//d//pSzw0/VkGsTpsOWsWbMkLy9Pli5dmtQ1R48ebbJ5wzTAaNSokbiBXevLNCBu1qxZ3IGSHq/veSq7NaTjNQAACIqEfns+9NBDUr16dfNnnep85pln5NFHHzXBniYxpIqOyp1++unSvn17M9KmGbi6Tis7O1uOHDkie/bsKXX8zp07zb5oMjMzzQhRyYdbaIBjx/3oVHOiBX/D3RratGljvqYi4ErHawAAEAQJ/QbdunWrCbDUnDlzTMLCLbfcYoKtDz74QNJFF+Br0KKBXpUqVWTRokXF+zRTVEeBdOrWq8Lr0KKJLDkTDQV/AQDwt4QCO13fFl63pfXSwuU0qlWrJgcPHpRU0ClTbVcWLtKrz5csWSLXXnutySoZNGiQmVJdvHixSabQvrUa1MWTNOBG4XVokSN3+ly3a7ePLl26mOzXWCj4CwCAvyW0xk4DucGDB8u5554r//znP+U3v/mN2b5u3Tpbm9CXtGvXLhk4cKDs2LHDBHJaH2/+/PnFQeX48ePNCFe/fv1KFSj2g1jr0PQ9qF+/vhk9TTQhQ0c/WecGAEAAAzstRHzfffeZKdm///3vcuKJJ5rtOlJ2zTXXSCpMmTIl6n4dLdT70ocfhdehVUSD3UQTMmJ1tgAAAD4uUOxnbipQnI6CvxV1tgjT6V6COwAAvBFvJJx+qEkS1113nVnf9a9//ctse+WVV+TDDz+07eZgb6JFZMFfDQZ1pC4aEi4AAPCOhAI7nX7VNWxa8uTjjz8u7tygUaeWQoE7Ey0iR96sdLYg4QIAAJ+vsXvwwQdl8uTJJplh+vTpxdu7du1q9sE58RT8tdquzK4OGIkisQMAgBQGdloj7sILLyyzXeeKI4sEw32JFuG1dZpVnM4OGIkgsQMAgBRPxWo3h40bN5bZruvrtHUV3C2cMHHgwIGYx+o0ro74OXmfkSOL+jzRThoAAPhZQoHdzTffbLIrV65caQrjbt++Xf73f/9XRo0aJbfeeqv9dwnbWEmYiJZwkS4kdgAAkKap2Hvuucf84r3kkkvMqI9Oy2rP1TvvvNMULoZ7151ZSZhQNWrUkMsuu8yxUifxJHakqig2AACBCOx0lO7ee+81gZxOye7bt09atWolzz77rDRt2lSKiorsv1PYsu7MaiKEZj07Wb/O6n06ndgBAICbxDXHpmVNtEdrhw4dTAbs3LlzTUCnrcQ0E1ML5I4YMSJ1d4uk151ZTYRwujiz1ft0MrEDAABPB3b333+/TJo0yUx9bdmyRa666iq55ZZbTJ/W//mf/zHb7r777tTdbcDZse5Mp2xjBW1OJkx47T4BAPBsYDdz5kx5+eWXZdasWbJgwQI5duyY/Pzzz/Lpp5/KgAEDpFKlSqm7U9hSUDiRDhVO8Mp9AgDgJnH9Vty2bZu0b9/e/Ll169YmYUKnXnXNHeynI2+FhYVSUFBgvtpVUDjeDhVO8cp9BvkzGW10GADg8uQJHaGrWrXqf06uXFlq1aqVivsKvPISJDRT1a51Z/F0qHCSV+4zCCgWDQA+C+xCoZDccMMNZqROHTp0SH73u99JzZo1Sx33+uuv23uXAU2QiGR3QWErHSrcwCv3GcTPZDhphxFUAPBgYJeXl1fq+XXXXWf3/QRevAWEI7HuDE4l7ejIKiOpAOChwG7q1KmpuxPEXUC45Aie1Tp2QLwoFg0APi9QjNSJp4CwBnPxrDtLtFsFgo1i0QDgHQR2LhNPAeF41p2x8B2p/kxSLBoAnMdwjcukojBvst0qEGwUiwYA7yCwcxm7C/Pa0a0CwUaxaADwDgI7F7KzMK8d3SoAikUDgDewxs6lyivM27BhQ9P9Q6v+W01+sLrwXYM77SRAYgXi+UySgAMA7kJg55HCvLoO7umnny41+malxInVBe3z58+nfAri+kwCANyHqVgPSCb5wcrC9/K6WpBYAQCA9xDYuVyyyQ9WFr4nem0AAOAuBHYuZ0fyQ0UL37V7RSwkVgAA4B2ssXOBaB0h7Kr6X97Cdw3aZs+enfS1AQCAOxDYOSxWR4jdu3dbuo6V4yIXvmsWrBV0FAAAwBuYinVxUsTChQtlyZIllq6lx8XbQYKOAgAA+AuBnYuTIvLz81Oa6EBHAQAA/IXAzsVJEaFQKK5rJpLooOvuunXrJtWrV0+6ywUAAHAWa+wckqqEhHiuW976Pg3wOnbsKBdccIHlfrQAAMAd+M3tkFQlJFi9bkXr+w4ePGjW623YsCEl9wcAAFKHwM4hVhIXMjIy4rqmXk+vm+qixwAAwJ0I7Jx64y10hOjcuXNc19TrWZk+taPoMQAAcB8COwdF6wjRv39/6dmzZ7n7I0fy4k10sKvoMQAAcBfPJE+MGzdOXn/9dfnyyy/NAv8uXbrII488YrI6ww4dOiSjRo2S6dOny+HDhyUnJ0cmTpwoWVlZ4lYajOmU59y5c+XAgQNmm35dsGCBGX0rr2NEw4YNZdu2beV2qrBzHR6FiQEA8BbPjNgtXbpUhgwZIitWrDCFe48ePSq9evWS/fv3Fx8zYsQIefPNN2XmzJnm+O3bt0vfvn3FzTSJYdasWcVBXWSRYt0f7hjRpk0b87Vy5cqlnsebvUphYgAA/CkjFG+xNJf49ttvpX79+iaAu/DCC2Xv3r1Sr149efXVV800ptLRPR3x0kK/nTp1snRdDahq165trhcruSFZOlL35JNPRl3vpvcwbNgw20uPhLNiK0INOwAAUidV8YZnRuwi6Ruh6tata76uWbPGjOL16NGj+JiWLVua0aloHRx0ylbf3JKPdHEyiaGi9X0UJgYAwLs8s8YucqRr+PDh0rVrV2ndurXZVlRUJFWrVpU6deqUOlbX1+m+aGv3xo4dK05IdRKDvk8l1+ZFrsUrb/1evOv1AACAe3gysNO1dp9//rl8+OGHSV9r9OjRMnLkyFIjZI0aNZJ0SGUSQ3ldJXQ0TkuilMyeDa/fAwAA3ue5oZmhQ4fKW2+9JYsXLzbZoWHZ2dly5MgR2bNnT6njd+7cafZVJDMz0wQ8JR/pkqokhoq6SpRMyAAAAP7jmcBOczw0qJs9e7a899570rRp01L727dvL1WqVJFFixYVb9O2WDrNGG+hXzcVKbZadDiMrhIAAARXZS9Nv2rG6xtvvGGmJsPr5jSjROva6ddBgwaZaVVNqNCRrttvv90EdVYzYp0QTmKwMm1qd0IGU7AAAPiLZwK7SZMmma/dunUrtX3q1Klyww03mD+PHz/ejG7169evVIFit0s2iSGcJKEB28aNGy2dQ1cJAAD8x7N17FIlnXXs7FBekoQVeXl5jNgBAOCzeMMzI3aIv8iwnQkZAADA/TyTPIH4kyTsSsgAAADewIidR1lJkoiUaEIGAADwBgI7j4o3+eGCCy4wiSeM1AEA4F/Mx3lUvN0omjVrRlAHAIDPEdh5lJWuFWEkSwAAEAwEdh5lpWtFGMkSAAAEA4Gdh4W7VlQ0cqfbdX84WUIzaQsLC6WgoMB81ecAAMA/SJ7wOA3aNECbO3euHDhwoHh7jRo1pFevXsVBXXmFjMmSBQDAXwjsPE4DtlmzZpXZrkGebg9nwZZXyFiDPN1eclQPAAB4F4Gdz4sUz5s3L+Z19Braq5ZSKAAAeBtr7HxepFjr3cWqeafX0GsBAABvY8QuQEWK03UtN4xkaqCqfyet96elYRiNBAAEAYFdgIoUp+taTiJJBAAQZEzF+rxIsQZssYI2vxQw1qBOk0Eip6fDSSK6HwAAPyOw83mR4t69e5uH3wsYW0kk0f3U7gMA+Jm3f5ujwiLFJYsTWzkmCIkkJIkAAPyONXY+oIGZliuJljBg5Rgvs5r84ackEQAAIhHY+YQGaE2aNEn6GK+ymvzhlyQRAADK44/hGgSelUQSvySJAABQEQI7BCaRxA9JIgAARMNvOfhGEJJEAACIhjV2DqE7Qmr4PUkEAIBoCOwcQHeE1PJzkggAANEwjJFmdEcAAACpQmCXRnRHAAAAqcRUbBrXylntjlBYWGiuwRoxAAAQDwK7NK6Vs9r1YNasWXLw4EFL1wQAAAhjKjaNa+Wsdj0oGdTFuiYAAEAYgV0a18pZ6Y4Q7zUBAADCCOwSYHWtnB4Xb3eEeK8JAAAQRmCXAKtr5co7rqLuCNWrV7f1tQEAQPCQPJEAq2vlKjquvO4IOsX6yiuv2PbaAAAgeAjsEhBeKxdtOlb363FWuyNoYJfsNQEAQLAxFZvIm2ZhrZzuj6c/aSquCQAAgoUoIUEVrZXT57o9kZpzqbgmAAAIjoxQKBRy+ibcRKdCa9euLXv37rVUmiSezhNWz0v0mgAAwJ/xhi/X2L3//vvy2GOPyZo1a2THjh0ye/ZsueKKK4r3a4w6ZswYef7552XPnj3StWtXmTRpkjRv3jxl9xS5Vs6OjhWJXBMAAMBTw0D79++Xtm3byoQJE8rd/+ijj8pTTz0lkydPlpUrV0rNmjUlJydHDh06JF7uWAEAAOC7EbvevXubR3l0tO6JJ56Q++67Ty6//HKz7eWXX5asrCyZM2eODBgwQLzSsUJLoTD1CgAAfD1iF82WLVukqKhIevToUbxN5647duwo+fn5FZ53+PBhM1pW8uG2jhUAAACBCuw0qFM6QleSPg/vK8+4ceNMABh+NGrUKGX3aLVrhAZ3hYWFUlBQYL7SHxYAAPhuKjYVRo8eLSNHjiwVVKUquLPaNWL+/Ply4MCBchMrAAAAfD9il52dbb7u3Lmz1HZ9Ht5XnszMTBM4lXykumNFLCWDOkViBQAACFRg17RpUxPALVq0qFRApNmxnTt3Fjew0l0iVmIF07IAAMAXgd2+fftk7dq15hFOmNA/a7JBRkaGDB8+XB588EH5xz/+YdanDRw4UBo0aFCq1p3TKuouUaNGjZjnklgBAAB8s8Zu9erV0r179+Ln4bVxeXl5Mm3aNLnrrrtMrbtbbrnFFCg+//zzzShXtWrVxE20nIlOAWtihNJixBq0asFluxIwAABA8NBSLE0tPmJ1nWjXrp0sWbIk5vndunWTiy66yPb7AgAA3o83PDUV63XRuk5oUFe9evWY19Dj6E4BAADKQ2Dnoq4TVpFEAQAAykNglyZWuk4cPHhQ6tev72gShQagFEcGAMCbPJU84WVWkx527dpl6/XsWP9HcWQAALyBEbs0sdp1wqnrRVv/p9tZ1wcAgPsR2KWJ1a4TVuh19HrpXP/Huj4AANyPwM4jXSdK0uvo9dK5/o/iyAAAuB+BnQu6Tlil5+n5eh07EyGsrtejODIAAO5G8kSaaVCmnSc++ugjmT9/vqVzLrjgAmnWrJmZfrU6UhdPIoTV9Xp2r+sDAAD2YsTOARqcnXfeeZYCJT1Gu01o27F4grp4EiGsrP+ze10fAACwH4GdQzRI6927d8zj9Jh41tMlkghhZf2f3ev6AACA/fhN7YI1d+W1EtNt8aynSzYRoqL1f/Gu6wMAAM5hjZ1L1txpcoM+VHgt3f79+822eNbWJZMIEb4XDfp0v04Dx/PaAADAWQR2LqCBkyZH6EPXv7355psJd39INhFC70XX8wEAAO9hKMZF7Oj+QCIEAADBRWDnEnZ1fyARAgCA4CKwcwk7uz+QCAEAQDCxxs4l7O7+QCIEAADBQ2DnEqno/kAiBAAAwcJUrEuQ9AAAAJJFYOcSJD0AAIBkEdi5CEkPAAAgGayxS5KWH4m3U0PJc2rWrGm2aZcJPV87P9D9AQAAJILALglaMFhry8XTJaK8c0qKp8sEAABASUzFprFLREXnWD0fAAAgGgK7NHWJsHJOtPMBAABiIbBLU5cIK+dEO1+DvMLCQikoKDBfCfoAAEAk1tilqUuE1XPKOz+RtXwAACB4GLFLU5eIeDpGlDwnkbV8AAAgmAjsEtCwYUPJyMiIeozu1+PiOackHZHTc+JdywcAAIKLwC4B27Ztk1AoFPUY3a/HxXNOSTrNqufEu5YPAAAEF4Gdy9bYVa9eXXJzc83auUReBwAABBfJEy5bY9e/f39p1qxZwq8DAACCixG7BGjbMF0DF43u1+PiPadJkyZJvQ4AAAguArtE3rTjjjNr4KLR/SV7xqbrHAAAEFxEBAnSNXC6Fi5yRE2fh9fIOXUOAAAIpoxQPKmaHjFhwgR57LHHpKioSNq2bStPP/20nHfeeZbO1SzT2rVry969e2NOg6rt27fLCy+8YDJetZzJ4MGDpUGDBlHP0fIkmsmqSQ+6Pk6nUmONuiVyDgAAcKd4443ABnavvfaaDBw4UCZPniwdO3aUJ554QmbOnCkbNmyQ+vXr2/pGjx07tsJ9Y8aMSej+AQCA//2YosDOd0M+jz/+uNx8881y4403SqtWrUyAV6NGDXnxxRdtfZ1oQZ2V/QAAAHbzVWB35MgRWbNmjfTo0aN4m05X6vP8/HzbXkeneO08DgAAwA6+Cuy+++47OXbsmGRlZZXars8rCrIOHz5shkNLPmJ5/vnnLd2P1eMAAADs4KvALhHjxo0zc9zhR6NGjWKeY7U3Kz1cAQBAOvkqsDvppJOkUqVKsnPnzlLb9Xl2dna554wePdosXAw/tm7dGvN1rGajkrUKAADSyVeBXdWqVaV9+/ayaNGiUqNm+rxz587lnpOZmWmyUUo+YtHkDCusHgcAAGAH3/WKHTlypOTl5UmHDh1M7Totd7J//36TJWuXikb/Ej0OAADADr4asVNXX321/PWvf5X7779fzjnnHFm7dq288847ZRIqkhWrTh117AAAQLr5rkBxugsGaratZr/qlK+uqdPpV0bqAABANHSe8PgbDQAAEEbnCQAAAARrjR0AAEBQ+S4r1ot0fd4333wjP/30kxx//PFy6qmnUgMPAADEjcDOYevXrzdZuyVbmenavksvvVTOPPNMR+8NAAB4C1OxDgd1M2bMKNOfVp/rdt0PAABgFYGdg9OvOlIXje6n3ywAALCKwM4huqYucqQuku7X4wAAAKxgjZ1DiQ96vBXh40iwAAAAsRDYOZT4oEGgFXocCRYAAMAKpmIdSnzQkb1YnS10//79+0mwAAAAlhDYOZT4oNO1OrIXTa9evWTBggVJvQ4AAAgOAjsHEx90ujY3N7fMyJ0+1+01a9YkwQIAAFjGGrsExJv4ECu4a968uaxevVp2794tdevWlQ4dOkjlypWloKDA1vsBAAD+RmCXgHgSH2IpLzEiPz/fTNPa+ToAAMD/mIpNgNXEBz0umQQMTZyw43UAAEAwENgl8qZZSHzQ/dHq2VlJwNDECU2gSOZ1AABAcBARJChW4kOsOnZWEzA0gSKZ1wEAAMHBGrskaFDVokWLlHeeaNOmTcKvAwAAgoPALkkaXDVp0iTu8+JNjEj0dQAAQHAw5OPxBAwAAIAwAjsPJ2AAAACURNTg4QQMAACAklhj5+EEDAAAgJII7FyAxAgAAGAHhoUAAAB8gsAOAADAJ5iKTZK2BmN9HAAAcAMCuySsX7/e9Hst2RpMM1q1TAkZrQAAIN2Yik0iqJsxY0aZfq/6XLfrfgAAgHQisEtw+lVH6qLR/XocAABAuhDYJUDX1EWO1EXS/XocAABAuhDYJUALCdt5HAAAgB0I7BKg3SHsPA4AAMAOBHYJ0JZfkf1dI+l+PQ4AACBdCOwSedOOO86UNIlG99PvFQAApBOBXYK0Tl1ubm6ZkTt9rtupYwcAANLNMwWK//KXv8jbb78ta9eulapVq8qePXvKHKNZqLfeeqssXrxYatWqJXl5eTJu3DipXDk1f00N3lq0aEHnCQAA4AqeCeyOHDkiV111lXTu3FmmTJlSZv+xY8ekT58+kp2dLcuXL5cdO3bIwIEDpUqVKvLQQw+l7L50urVJkyYpuz4AAIBVGaFQKCQeMm3aNBk+fHiZEbt58+bJZZddJtu3b5esrCyzbfLkyXL33XfLt99+a0b5rND6c7Vr15a9e/fGTJAAAABIRKriDd+sscvPz5c2bdoUB3UqJyfHvHHr1q1z9N4AAADSwTNTsbEUFRWVCupU+Lnuq8jhw4fNI0wjZxWrswQAAECiwnGG3ROnjgZ299xzjzzyyCNRj1m/fr20bNkyZfegyRVjx44ts71Ro0Ype00AAAD1/fffmylZXwR2o0aNkhtuuCHqMc2aNbN0LU2a+Oijj0pt27lzZ/G+iowePVpGjhxZ/PyXX36R3bt3y4knnigZGRmWo24NBLdu3RrodXm8D7wPfBb4nuBnAz8j+V1hjc4QaiODunXrip0cDezq1atnHnbQbFktibJr1y6pX7++2bZw4UITaLVq1arC8zIzM82jpDp16iR0D/paQQ7swngfeB/4LPA9wc8Gfkbyu8Iau5sZeGaNndao05E0/aqlTbSenTr99NNNzbpevXqZAO7666+XRx991Kyru++++2TIkCFlAjcAAAA/8kxgd//998tLL71U/Pzcc881X7UYcbdu3aRSpUry1ltvmQLFOnpXs2ZNU6D4T3/6k4N3DQAAkD6VvVS/Th/RNG7cWObOnSvppiOCY8aMCfzIIO8Dnwe+J/jZwM9IflfwO9PZ35meK1AMAAAAnxcoBgAACDoCOwAAAJ8gsAMAAPAJAjuLJkyYIE2aNJFq1apJx44dyxRDjjRz5kzTMUOP1x62TiR1OP0+aLKLFnku+dDzvOz999+X//qv/5IGDRqYv8+cOXNinrNkyRJp166dWSCr5XliJQH58X3Q9yDys6CPaO3+3E671vz617+W448/3tTOvOKKK2TDhg0xz/Pbz4ZE3gc//myYNGmSnH322cV1PLU6w7x58wL1WUjkffDjZ6E8Dz/8sPm7DR8+XFL9mSCws+C1114z3Sk0e+Xjjz+Wtm3bSk5OjimGXJ7ly5fLNddcI4MGDZJPPvnE/KDTx+effy5Beh+UfmPv2LGj+PH111+Ll+3fv9/8vTXAtWLLli3Sp08f6d69u6m9qN/UgwcPlvnz50uQ3ocw/YVf8vMQLibuRUuXLjV1MlesWGGKoR89etTU09T3piJ+/NmQyPvgx58NDRs2NL+816xZI6tXr5aLL75YLr/8clm3bl1gPguJvA9+/CxEWrVqlTz77LMm4I3Gts+EZsUiuvPOOy80ZMiQ4ufHjh0LNWjQIDRu3Lhyj8/NzQ316dOn1LaOHTuGfvvb3wbqfZg6dWqodu3aIb/Sb5/Zs2dHPeauu+4KnXXWWaW2XX311aGcnJxQkN6HxYsXm+N++OGHkF/t2rXL/B2XLl1a4TF+/dkQ7/vg958NYSeccELohRdeCOxnwcr74PfPwk8//RRq3rx5aOHChaGLLrooNGzYsAqPteszwYhdDEeOHDH/8ujRo0ep9h/6PD8/v9xzdHvJ45WObFV0vF/fB7Vv3z5TX1B76cb6V5sf+fGzkIxzzjlHTj75ZOnZs6csW7ZM/Nb3UUXr+xiEz4OV98HvPxu0O9L06dPNqKVORQb1s2DlffD7Z2HIkCFm1iby/3UqPxMEdjF899135sOZlZVVars+r2h9kG6P53i/vg8tWrSQF198Ud544w3529/+Jr/88ot06dJFtm3bJkFR0Wfhxx9/lIMHD0pQaDA3efJk+fvf/24e+gNcO8bolL4f6Gdbp9m7du0qrVu3rvA4P/5sSOR98OvPhoKCAtPiUtfT/u53v5PZs2dX2Kvcz5+FeN4Hv34WlAa1+jNO16FaYddnwjOdJ+A9+i+0kv9K02/WM88806w1+POf/+zovSG99Ie3Pkp+FjZt2iTjx4+XV155xRf/Ktd1MB9++KEEmdX3wa8/G/QzrmtpddRy1qxZpq2lrkGsKKjxq3jeB79+FrZu3SrDhg0z607TnQxCYBfDSSedZPrQ7ty5s9R2fZ6dnV3uObo9nuP9+j5EqlKliunxu3HjRgmKij4Luli4evXqEmTnnXeeLwKhoUOHmj7VmimsC8ej8ePPhkTeB7/+bKhatarJfFft27c3i+affPJJE6QE6bMQz/vg18/CmjVrTGKhVkQI01kv/f545pln5PDhw+Z3aio+E0zFWviA6gdz0aJFxdt0qFifV7RmQLeXPF5p1B5tjYEf34dI+qHWIXqdlgsKP34W7KL/ovfyZ0HzRjSY0Wmm9957T5o2bRrIz0Mi70NQfjboz0j9BR6Uz0Ii74NfPwuXXHKJ+Xvoz7nwo0OHDnLttdeaP0cGdbZ+JhJO9QiQ6dOnhzIzM0PTpk0LffHFF6FbbrklVKdOnVBRUZHZf/3114fuueee4uOXLVsWqly5cuivf/1raP369aExY8aEqlSpEiooKAgF6X0YO3ZsaP78+aFNmzaF1qxZExowYECoWrVqoXXr1oW8nOH0ySefmId++zz++OPmz19//bXZr39/fR/CNm/eHKpRo0bozjvvNJ+FCRMmhCpVqhR65513Ql4W7/swfvz40Jw5c0JfffWV+T7QzLDjjjsu9O6774a86tZbbzXZfEuWLAnt2LGj+HHgwIHiY4LwsyGR98GPPxv076eZwFu2bAl99tln5nlGRkZowYIFgfksJPI++PGzUJHIrNhUfSYI7Cx6+umnQ6eeemqoatWqpuzHihUrSv3PysvLK3X8jBkzQmeccYY5XstdvP3226GgvQ/Dhw8vPjYrKyv0m9/8JvTxxx+HvCxctiPyEf5761d9HyLPOeecc8z70KxZM5Pe73Xxvg+PPPJI6LTTTjM/sOvWrRvq1q1b6L333gt5WXl/f32U/P8bhJ8NibwPfvzZcNNNN4UaN25s/k716tULXXLJJcXBTFA+C4m8D378LFgN7FL1mcjQ/yQ75AgAAADnscYOAADAJwjsAAAAfILADgAAwCcI7AAAAHyCwA4AAMAnCOwAAAB8gsAOAADAJwjsAAAAfILADgAAwCcI7ABARG644QbJyMgwjypVqkhWVpb07NlTXnzxRdPEHAC8gMAOAP7t0ksvlR07dkhhYaHMmzdPunfvLsOGDZPLLrtMfv75Z94nAK5HYAcA/5aZmSnZ2dlyyimnSLt27eQPf/iDvPHGGybImzZtmjnmm2++kcsvv1xq1aolv/rVryQ3N1d27txZ/B5++umnJiA8/vjjzf727dvL6tWri/d/+OGHcsEFF0j16tWlUaNGcscdd8j+/fv5fwDAFgR2ABDFxRdfLG3btpXXX3/dTMlqULd7925ZunSpLFy4UDZv3ixXX3118fHXXnutNGzYUFatWiVr1qyRe+65x0ztqk2bNplRwX79+slnn30mr732mgn0hg4dyv8DALbICIVCIXsuBQDeXmO3Z88emTNnTpl9AwYMMIHYk08+Kb1795YtW7aY0Tb1xRdfyFlnnSUfffSR/PrXvzajdE8//bTk5eWVuc7gwYOlUqVK8uyzzxZv08DuoosuMqN21apVS/HfEoDfMWIHADHov381qWL9+vUmoAsHdapVq1ZSp04ds0+NHDnSBHA9evSQhx9+2IzSlZym1SldncYNP3JycsxIoAaLAJAsAjsAiEGDtqZNm1p6nx544AFZt26d9OnTR9577z0T+M2ePdvs27dvn/z2t7+VtWvXFj802Pvqq6/ktNNO4/8DgKRVTv4SAOBfGpwVFBTIiBEjzNq5rVu3mkfJqVidwtUALuyMM84wDz3nmmuukalTp8qVV15pEjL0+NNPP93BvxEAPyOwA4B/O3z4sBQVFcmxY8dMpus777wj48aNM+VOBg4cKMcdd5y0adPGJEg88cQTpgTKbbfdZtbIdejQQQ4ePCh33nmn9O/f34zwbdu2zSRRaLKEuvvuu6VTp04mWUKna2vWrGkCPU3CeOaZZ/j/ACBpBHYA8G8ayJ188slSuXJlOeGEE0w27FNPPWUSITSoU1r+5Pbbb5cLL7zQbNMsV02WUJoY8f3335sgUAPDk046Sfr27Stjx441+88++2yTTXvvvfeakie6dk+nYEtm1QJAMsiKBQAA8AmSJwAAAHyCwA4AAMAnCOwAAAB8gsAOAADAJwjsAAAAfILADgAAwCcI7AAAAHyCwA4AAMAnCOwAAAB8gsAOAADAJwjsAAAAfILADgAAQPzh/wHZfVUxrAA22wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(d['zifro'], d['acetyl'], 'o', color='gray')\n",
    "plt.xlabel(\"Dose\")\n",
    "plt.ylim([-10, 75])\n",
    "plt.ylabel(\"Response\")\n",
    "plt.xlim([-0.1, 4])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1c7ed44-3645-412b-91ca-4dbaeefd298f",
   "metadata": {},
   "source": [
    "While we will ignore the time-since-dose in this analysis, we should not ignore the fact that a single unit contributes multiple times to estimation. Our variance needs to account for this. Let's first ignore this aspect and treat all the observations as if they were independent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c19d28df-c6e5-4af6-ab85-a2a2ff78bcc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    return ee_emax(theta=theta, dose=d['zifro'], response=d['acetyl'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "35b4c18b-72df-4172-aeb9-1f6adb591a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "estr = MEstimator(psi, init=[0, 60, 0.5])\n",
    "estr.estimate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "103dd97a-c73d-4caa-b55e-d46126c0572b",
   "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>Est</th>\n",
       "      <th>SE</th>\n",
       "      <th>LCL</th>\n",
       "      <th>UCL</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Params</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Zero-Dose</th>\n",
       "      <td>-2.3</td>\n",
       "      <td>1.3</td>\n",
       "      <td>-4.9</td>\n",
       "      <td>0.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>E-max</th>\n",
       "      <td>84.5</td>\n",
       "      <td>6.6</td>\n",
       "      <td>71.4</td>\n",
       "      <td>97.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ED50</th>\n",
       "      <td>0.6</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Est   SE   LCL   UCL\n",
       "Params                          \n",
       "Zero-Dose  -2.3  1.3  -4.9   0.3\n",
       "E-max      84.5  6.6  71.4  97.5\n",
       "ED50        0.6  0.1   0.4   0.8"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = pd.DataFrame()\n",
    "results['Params'] = [\"Zero-Dose\", \"E-max\", \"ED50\"]\n",
    "results = results.set_index('Params')\n",
    "results['Est'] = estr.theta\n",
    "results['SE'] = np.diag(estr.variance) ** 0.5\n",
    "ci = estr.confidence_intervals()\n",
    "results['LCL'] = ci[:, 0]\n",
    "results['UCL'] = ci[:, 1]\n",
    "results.round(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8370fba6-b94b-4adf-919e-9c229cf8e773",
   "metadata": {},
   "source": [
    "To compare, let's modify the estimating functions. Here, we will use the `aggregate_efuncs` to collapse the contributions by unique subject IDs. The dimension of this new estimating functions will be 3-by-12, to reflect the fact that there are only 12 independent observations. The following code applies this process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e6be5054-e12e-477e-bf66-6de644b4720d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi_c(theta):\n",
    "    ef_emax_i = ee_emax(theta=theta, dose=d['zifro'], response=d['acetyl'])\n",
    "    return aggregate_efuncs(ef_emax_i, group=d['subject'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "0c16ebed-0c33-4105-bbd1-f05c00975724",
   "metadata": {},
   "outputs": [],
   "source": [
    "estr = MEstimator(psi_c, init=[0, 60, 0.5])\n",
    "estr.estimate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "3a65b0bc-b68b-4714-ab47-7971c9c49884",
   "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>Est</th>\n",
       "      <th>SE</th>\n",
       "      <th>LCL</th>\n",
       "      <th>UCL</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Params</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Zero-Dose</th>\n",
       "      <td>-2.3</td>\n",
       "      <td>1.6</td>\n",
       "      <td>-5.3</td>\n",
       "      <td>0.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>E-max</th>\n",
       "      <td>84.5</td>\n",
       "      <td>7.4</td>\n",
       "      <td>70.0</td>\n",
       "      <td>99.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ED50</th>\n",
       "      <td>0.6</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.7</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Est   SE   LCL   UCL\n",
       "Params                          \n",
       "Zero-Dose  -2.3  1.6  -5.3   0.8\n",
       "E-max      84.5  7.4  70.0  99.0\n",
       "ED50        0.6  0.1   0.4   0.7"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = pd.DataFrame()\n",
    "results['Params'] = [\"Zero-Dose\", \"E-max\", \"ED50\"]\n",
    "results = results.set_index('Params')\n",
    "results['Est'] = estr.theta\n",
    "results['SE'] = np.diag(estr.variance) ** 0.5\n",
    "ci = estr.confidence_intervals()\n",
    "results['LCL'] = ci[:, 0]\n",
    "results['UCL'] = ci[:, 1]\n",
    "results.round(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79486e3e-85cf-449b-a5ef-eb998da3be82",
   "metadata": {},
   "source": [
    "As seen here, the point estimates are similar but the variance changes.\n",
    "\n",
    "However, 12 observations is still relatively few so we can further apply a finite-sample correction using the `finite_correction` optional argument as shown in the following code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "c372eb62-4524-43ed-9fb9-04dfad5f13c6",
   "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>Est</th>\n",
       "      <th>SE</th>\n",
       "      <th>LCL</th>\n",
       "      <th>UCL</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Params</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Zero-Dose</th>\n",
       "      <td>-2.3</td>\n",
       "      <td>1.8</td>\n",
       "      <td>-6.3</td>\n",
       "      <td>1.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>E-max</th>\n",
       "      <td>84.5</td>\n",
       "      <td>8.5</td>\n",
       "      <td>65.1</td>\n",
       "      <td>103.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ED50</th>\n",
       "      <td>0.6</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Est   SE   LCL    UCL\n",
       "Params                           \n",
       "Zero-Dose  -2.3  1.8  -6.3    1.8\n",
       "E-max      84.5  8.5  65.1  103.8\n",
       "ED50        0.6  0.1   0.4    0.8"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estr = MEstimator(psi_c, init=[0, 60, 0.5], finite_correction='HC1')\n",
    "estr.estimate()\n",
    "\n",
    "results = pd.DataFrame()\n",
    "results['Params'] = [\"Zero-Dose\", \"E-max\", \"ED50\"]\n",
    "results = results.set_index('Params')\n",
    "results['Est'] = estr.theta\n",
    "results['SE'] = np.diag(estr.variance) ** 0.5\n",
    "ci = estr.confidence_intervals()\n",
    "results['LCL'] = ci[:, 0]\n",
    "results['UCL'] = ci[:, 1]\n",
    "results.round(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5606b501",
   "metadata": {},
   "source": [
    "Here, we see the variance increase a bit further (due to the additional correction). This concludes the case study showing how to account for non-independent observations and small sample sizes.\n",
    "\n",
    "## Chapter 11: Generalized Linear Models and Its Extensions\n",
    "\n",
    "\n",
    "### Adverse Events Case Study\n",
    "\n",
    "The first example from Chapter 11 is the *Case Study: Assessing the Relationship Between Drug Concentrations and Adverse Events Using Logistic Regression*. Data comes from Table 2 of the book. In the book, a variety of different models for different adverse events are considered. Here, we only consider nausea (and vomiting) by AUC. For the one observations with a missing AUC value, they are dropped from the data set (same as the book). For ease of examining the coefficients, we will also divide the AUC value by 1000. This means the coefficients for AUC are rescaled from those reported in the book.\n",
    "\n",
    "First, we load the data set and transform the coefficients and add an intercept column to the data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1be76e83",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Index: 42 entries, 0 to 42\n",
      "Data columns (total 12 columns):\n",
      " #   Column     Non-Null Count  Dtype  \n",
      "---  ------     --------------  -----  \n",
      " 0   id         42 non-null     int64  \n",
      " 1   c_max      42 non-null     float64\n",
      " 2   auc        42 non-null     float64\n",
      " 3   age        42 non-null     int64  \n",
      " 4   sex        42 non-null     int64  \n",
      " 5   ps         42 non-null     int64  \n",
      " 6   myalgia    42 non-null     int64  \n",
      " 7   phlebitis  42 non-null     int64  \n",
      " 8   asthenia   42 non-null     int64  \n",
      " 9   diarrhea   42 non-null     int64  \n",
      " 10  nausea     42 non-null     int64  \n",
      " 11  intercept  42 non-null     int64  \n",
      "dtypes: float64(2), int64(10)\n",
      "memory usage: 4.3 KB\n"
     ]
    }
   ],
   "source": [
    "d = pd.read_csv(\"data/bonate.csv\").dropna()\n",
    "d['intercept'] = 1   # Adding intercept to data\n",
    "d['auc'] /= 1000     # Rescaling AUC\n",
    "d['c_max'] /= 1000   # Rescaling C_max\n",
    "d.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "143778ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "table = pd.DataFrame(columns=[\"Model\", \"Intercept\", \"AUC\", \"Sex\", \"Age\", \"PS\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02260fdc",
   "metadata": {},
   "source": [
    "To begin, we will fit a null (intercept-only) logistic regression model. This is easily done by using the built-in `ee_glm` estimating equation. For the logistic model, we specify a binomial distribution with the logit link. Below is code to setup the estimating equation and then estimate the parameters using `MEstimator`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "207e9a4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    # Estimating equation for null model\n",
    "    return ee_glm(theta=theta,\n",
    "                  y=d['nausea'],\n",
    "                  X=d[['intercept', ]],\n",
    "                  distribution='binomial',\n",
    "                  link='logit')\n",
    "\n",
    "\n",
    "# Estimate the parameters of the logit model\n",
    "estr_null = MEstimator(psi, init=[0., ])\n",
    "estr_null.estimate()\n",
    "\n",
    "# Adding results to the output table\n",
    "table.loc[len(table)] = [\"Null\", estr_null.theta[0], ] + [np.nan, ]*4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e0cb694",
   "metadata": {},
   "source": [
    "Next we fit a logistic regression model that includes linear terms for all the independent variables in the data set. This is easily done by modifying the previous design matrix (i.e., `X`). Below is code to fit the full model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8d17c1bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    # Estimating equation for full model\n",
    "    return ee_glm(theta=theta,\n",
    "                  y=d['nausea'],\n",
    "                  X=d[['intercept', 'auc', 'sex', 'age', 'ps']],\n",
    "                  distribution='binomial',\n",
    "                  link='logit')\n",
    "\n",
    "\n",
    "# Estimate the parameters of the logit model\n",
    "estr_full = MEstimator(psi, init=[0., ]*5)\n",
    "estr_full.estimate()\n",
    "\n",
    "# Adding results to the output table\n",
    "table.loc[len(table)] = [\"Full\", ] + list(estr_full.theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f579eda",
   "metadata": {},
   "source": [
    "In the book, Bonate performs some variable selection. In general, we would not recommend use of backwards-selection procedures (like those done in the book). Such procedures complicate inference (P-values and confidence intervals after these procedures are no longer valid). For comparison purposes, we estimate the reduced model reported in the book. Again, this is easily done by modifying the `X` argument for `ee_glm`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6558123b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    # Estimating equation for reduced model\n",
    "    return ee_glm(theta=theta,\n",
    "                  y=d['nausea'],\n",
    "                  X=d[['intercept', 'auc', 'sex']],\n",
    "                  distribution='binomial',\n",
    "                  link='logit')\n",
    "\n",
    "\n",
    "# Estimate the parameters of the logit model\n",
    "estr_redu = MEstimator(psi, init=[0., ]*3)\n",
    "estr_redu.estimate()\n",
    "\n",
    "# Adding results to the output table\n",
    "table.loc[len(table)] = [\"Reduced\", ] + list(estr_redu.theta) + [np.nan, ]*2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c8ff25",
   "metadata": {},
   "source": [
    "Finally, two alternative models are considered: a probit regression model and a complimentary log-log model. Again, these models are easily implemented using `ee_glm`. For the probit model, we set the link equal to `probit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "24a66642",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    # Estimating equation for reduced probit model\n",
    "    return ee_glm(theta=theta,\n",
    "                  y=d['nausea'],\n",
    "                  X=d[['intercept', 'auc', 'sex']],\n",
    "                  distribution='binomial',\n",
    "                  link='probit')\n",
    "\n",
    "\n",
    "# Estimate the parameters of the probit model\n",
    "estr_prob = MEstimator(psi, init=[0., ]*3)\n",
    "estr_prob.estimate()\n",
    "\n",
    "# Adding results to the output table\n",
    "table.loc[len(table)] = [\"Probit\", ] + list(estr_prob.theta) + [np.nan, ]*2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "643eb533",
   "metadata": {},
   "source": [
    "Similarly, the complimentary log-log model only requires setting the link to `cloglog`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "6858e649",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(theta):\n",
    "    # Estimating equation for reduced C-log-log model\n",
    "    return ee_glm(theta=theta,\n",
    "                  y=d['nausea'],\n",
    "                  X=d[['intercept', 'auc', 'sex']],\n",
    "                  distribution='binomial',\n",
    "                  link='cloglog')\n",
    "\n",
    "\n",
    "# Estimate the parameters of the cloglog model\n",
    "estr_clog = MEstimator(psi, init=[0., ]*3)\n",
    "estr_clog.estimate()\n",
    "\n",
    "# Adding results to the output table\n",
    "table.loc[len(table)] = [\"CLogLog\", ] + list(estr_clog.theta) + [np.nan, ]*2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c28fe617",
   "metadata": {},
   "source": [
    "Now we can view the results across the different models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d90c3a58",
   "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>Intercept</th>\n",
       "      <th>AUC</th>\n",
       "      <th>Sex</th>\n",
       "      <th>Age</th>\n",
       "      <th>PS</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Model</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Null</th>\n",
       "      <td>-0.587787</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Full</th>\n",
       "      <td>-5.602899</td>\n",
       "      <td>0.289785</td>\n",
       "      <td>1.730299</td>\n",
       "      <td>0.049538</td>\n",
       "      <td>0.220545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Reduced</th>\n",
       "      <td>-2.663522</td>\n",
       "      <td>0.303502</td>\n",
       "      <td>1.772238</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Probit</th>\n",
       "      <td>-1.629729</td>\n",
       "      <td>0.184030</td>\n",
       "      <td>1.080253</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLogLog</th>\n",
       "      <td>-2.233720</td>\n",
       "      <td>0.193108</td>\n",
       "      <td>1.212290</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Intercept       AUC       Sex       Age        PS\n",
       "Model                                                     \n",
       "Null     -0.587787       NaN       NaN       NaN       NaN\n",
       "Full     -5.602899  0.289785  1.730299  0.049538  0.220545\n",
       "Reduced  -2.663522  0.303502  1.772238       NaN       NaN\n",
       "Probit   -1.629729  0.184030  1.080253       NaN       NaN\n",
       "CLogLog  -2.233720  0.193108  1.212290       NaN       NaN"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table.set_index(\"Model\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4791334a",
   "metadata": {},
   "source": [
    "These point estimates match those reported in Tables 3 and 4 on pages 469 and 470 in the book (note the null model differs slightly, since we dropped the one observation with the missing AUC value to fit this model, but the book does not). These results highlight how `delicatessen` allows one to easily fit a variety of different models. \n",
    "\n",
    "This is the end of the current replication."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}