{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Quasi-binomial regression\n", "\n", "This notebook demonstrates using custom variance functions and non-binary data\n", "with the quasi-binomial GLM family to perform a regression analysis using\n", "a dependent variable that is a proportion.\n", "\n", "The notebook uses the barley leaf blotch data that has been discussed in\n", "several textbooks. See below for one reference:\n", "\n", "https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_sect016.htm" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:10.322649Z", "iopub.status.busy": "2022-11-02T17:08:10.322374Z", "iopub.status.idle": "2022-11-02T17:08:12.218975Z", "shell.execute_reply": "2022-11-02T17:08:12.218257Z" } }, "outputs": [], "source": [ "import statsmodels.api as sm\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from io import StringIO" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The raw data, expressed as percentages. We will divide by 100\n", "to obtain proportions." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:12.223182Z", "iopub.status.busy": "2022-11-02T17:08:12.222616Z", "iopub.status.idle": "2022-11-02T17:08:12.227384Z", "shell.execute_reply": "2022-11-02T17:08:12.226826Z" } }, "outputs": [], "source": [ "raw = StringIO(\n", " \"\"\"0.05,0.00,1.25,2.50,5.50,1.00,5.00,5.00,17.50\n", "0.00,0.05,1.25,0.50,1.00,5.00,0.10,10.00,25.00\n", "0.00,0.05,2.50,0.01,6.00,5.00,5.00,5.00,42.50\n", "0.10,0.30,16.60,3.00,1.10,5.00,5.00,5.00,50.00\n", "0.25,0.75,2.50,2.50,2.50,5.00,50.00,25.00,37.50\n", "0.05,0.30,2.50,0.01,8.00,5.00,10.00,75.00,95.00\n", "0.50,3.00,0.00,25.00,16.50,10.00,50.00,50.00,62.50\n", "1.30,7.50,20.00,55.00,29.50,5.00,25.00,75.00,95.00\n", "1.50,1.00,37.50,5.00,20.00,50.00,50.00,75.00,95.00\n", "1.50,12.70,26.25,40.00,43.50,75.00,75.00,75.00,95.00\"\"\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The regression model is a two-way additive model with\n", "site and variety effects. The data are a full unreplicated\n", "design with 10 rows (sites) and 9 columns (varieties)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:12.230990Z", "iopub.status.busy": "2022-11-02T17:08:12.230548Z", "iopub.status.idle": "2022-11-02T17:08:12.247411Z", "shell.execute_reply": "2022-11-02T17:08:12.246859Z" } }, "outputs": [], "source": [ "df = pd.read_csv(raw, header=None)\n", "df = df.melt()\n", "df[\"site\"] = 1 + np.floor(df.index / 10).astype(int)\n", "df[\"variety\"] = 1 + (df.index % 10)\n", "df = df.rename(columns={\"value\": \"blotch\"})\n", "df = df.drop(\"variable\", axis=1)\n", "df[\"blotch\"] /= 100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the quasi-binomial regression with the standard variance\n", "function." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:12.250881Z", "iopub.status.busy": "2022-11-02T17:08:12.250504Z", "iopub.status.idle": "2022-11-02T17:08:12.295029Z", "shell.execute_reply": "2022-11-02T17:08:12.294418Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: blotch No. Observations: 90\n", "Model: GLM Df Residuals: 72\n", "Model Family: Binomial Df Model: 17\n", "Link Function: Logit Scale: 0.088778\n", "Method: IRLS Log-Likelihood: -20.791\n", "Date: Wed, 02 Nov 2022 Deviance: 6.1260\n", "Time: 17:08:12 Pearson chi2: 6.39\n", "No. Iterations: 10 Pseudo R-squ. (CS): 0.3198\n", "Covariance Type: nonrobust \n", "==================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "----------------------------------------------------------------------------------\n", "C(variety)[1] -8.0546 1.422 -5.664 0.000 -10.842 -5.268\n", "C(variety)[2] -7.9046 1.412 -5.599 0.000 -10.672 -5.138\n", "C(variety)[3] -7.3652 1.384 -5.321 0.000 -10.078 -4.652\n", "C(variety)[4] -7.0065 1.372 -5.109 0.000 -9.695 -4.318\n", "C(variety)[5] -6.4399 1.357 -4.746 0.000 -9.100 -3.780\n", "C(variety)[6] -5.6835 1.344 -4.230 0.000 -8.317 -3.050\n", "C(variety)[7] -5.4841 1.341 -4.090 0.000 -8.112 -2.856\n", "C(variety)[8] -4.7126 1.331 -3.539 0.000 -7.322 -2.103\n", "C(variety)[9] -4.5546 1.330 -3.425 0.001 -7.161 -1.948\n", "C(variety)[10] -3.8016 1.320 -2.881 0.004 -6.388 -1.215\n", "C(site)[T.2] 1.6391 1.443 1.136 0.256 -1.190 4.468\n", "C(site)[T.3] 3.3265 1.349 2.466 0.014 0.682 5.971\n", "C(site)[T.4] 3.5822 1.344 2.664 0.008 0.947 6.217\n", "C(site)[T.5] 3.5831 1.344 2.665 0.008 0.948 6.218\n", "C(site)[T.6] 3.8933 1.340 2.905 0.004 1.266 6.520\n", "C(site)[T.7] 4.7300 1.335 3.544 0.000 2.114 7.346\n", "C(site)[T.8] 5.5227 1.335 4.138 0.000 2.907 8.139\n", "C(site)[T.9] 6.7946 1.341 5.068 0.000 4.167 9.422\n", "==================================================================================\n" ] } ], "source": [ "model1 = sm.GLM.from_formula(\n", " \"blotch ~ 0 + C(variety) + C(site)\", family=sm.families.Binomial(), data=df\n", ")\n", "result1 = model1.fit(scale=\"X2\")\n", "print(result1.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot below shows that the default variance function is\n", "not capturing the variance structure very well. Also note\n", "that the scale parameter estimate is quite small." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:12.298368Z", "iopub.status.busy": "2022-11-02T17:08:12.297965Z", "iopub.status.idle": "2022-11-02T17:08:12.527953Z", "shell.execute_reply": "2022-11-02T17:08:12.524895Z" }, "lines_to_next_cell": 1 }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Residual')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVBklEQVR4nO3deXxTVd4/8E9S2pRiFwq0KVBoxQVKpQU67dQVhS6gHVHHHwgoMCwDY12oj0pHpaJiXUAZkQGXARfkgWfQEZFa6YCIYrGyzrALFnCgC1BoaittaO7vj05CQ5beJDfJvbmf9+ulL3Jz7r0np2nzzTnfc45GEAQBRERERCqk9XcFiIiIiPyFgRARERGpFgMhIiIiUi0GQkRERKRaDISIiIhItRgIERERkWoxECIiIiLV6uTvCsidyWTCqVOnEB4eDo1G4+/qEBERkQiCIKChoQE9e/aEVuu434eBUAdOnTqF+Ph4f1eDiIiI3PDzzz+jd+/eDp9nINSB8PBwAG0NGRER4ZV7GI1GbNiwAdnZ2QgODvbKPQId21AabEdpsB09xzaUhprb0WAwID4+3vI57ggDoQ6Yh8MiIiK8GgiFhYUhIiJCdW9UqbANpcF2lAbb0XNsQ2mwHdFhWguTpYmIiEi1GAgRERGRaikqENqyZQvy8vLQs2dPaDQafPrppx2es3nzZgwZMgQ6nQ5XXXUV3nvvPa/Xk4iIiJRBUYFQY2MjUlJSsHjxYlHlKysrcfvtt+PWW2/F7t278eijj2Lq1Kn48ssvvVxTIiIiUgJFJUuPHDkSI0eOFF1+6dKlSExMxIIFCwAAAwYMwLfffovXX38dOTk53qomERERKYSiAiFXlZeXY8SIEVbHcnJy8Oijjzo8p7m5Gc3NzZbHBoMBQFvmvdFo9Eo9zdf11vXVgG0oDbajNNiOnmMbSkPN7Sj2NQd0IFRdXY3Y2FirY7GxsTAYDPj111/RuXNnm3OKi4sxd+5cm+MbNmxAWFiY1+oKAGVlZV69vhqwDaXBdpQG29FzbENpqLEdm5qaRJUL6EDIHYWFhSgoKLA8Ni/IlJ2d7dV1hMrKypCVlaXadR48xTaUBttRGmxHz7ENpaHmdjSP6HQkoAMhvV6Pmpoaq2M1NTWIiIiw2xsEADqdDjqdzuZ4cHCw199EvrhHoGMbSoPtKA22o+fYhtJQYzuKfb0BHQhlZmaipKTE6lhZWRkyMzP9VCMioktaTQIqKutQ23ABMeGhSE+MRpCWmzsT+ZKiAqFffvkFR44csTyurKzE7t27ER0djT59+qCwsBAnT57EBx98AACYMWMG3nzzTTzxxBP4wx/+gE2bNuH//u//sH79en+9BCIiAEDp3irMXbcfVfUXLMfiIkNRlJeE3OQ4P9aMSF0UtY7Q9u3bMXjwYAwePBgAUFBQgMGDB2POnDkAgKqqKpw4ccJSPjExEevXr0dZWRlSUlKwYMECvPvuu5w6T0R+Vbq3CjNX7LQKggCguv4CZq7YidK9VX6qGZH6KKpHaNiwYRAEweHz9laNHjZsGHbt2uXFWhERiddqEjB33X7Y+0smANAAmLtuP7KS9BwmI/IBRfUIEREpXUVlnU1PUHsCgKr6C6iorPNdpYhUjIEQEZEP1TY4DoLcKUdEnmEgRETkQzHhoZKWIyLPMBAiIvKh9MRoxEWGwlH2jwZts8fSE6N9WS0i1WIgRETkQ0FaDYrykgDAJhgyPy7KS2KiNJGPMBAiIvKx3OQ4LJkwBPpI6+EvfWQolkwYwnWEiHxIUdPniYgCRW5yHLKS9FxZmsjPGAgRkdu4RYRngrQaZPbr5u9qEKkaAyEicgu3iCCiQMAcISJyGbeIIKJAwUCIiFzS0RYRQNsWEa0mx9vhEBHJBQMhInIJt4ggokDCQIiIXMItIogokDAQIiKXcIsIIgokDISIyCXcIoKIAgkDISJyCbeIIKJAwkCIiFzGLSKIKFBwQUUicgu3iCCiQMBAiIjcxi0iiEjpODRGREREqsVAiIiIiFSLgRARERGpFgMhIiIiUi0GQkRERKRaDISIiIhItRgIERERkWoxECIiIiLVYiBEREREqsVAiIiIiFSLgRARERGpFgMhIiIiUi0GQkRERKRaDISIiIhItRQXCC1evBgJCQkIDQ1FRkYGKioqnJZfuHAhrr32WnTu3Bnx8fGYNWsWLly44KPaEhERkZwpKhBavXo1CgoKUFRUhJ07dyIlJQU5OTmora21W37lypWYPXs2ioqKcODAAfztb3/D6tWr8ec//9nHNSciIiI5UlQg9Nprr2HatGmYPHkykpKSsHTpUoSFhWHZsmV2y3/33Xe44YYbMG7cOCQkJCA7Oxv33Xdfh71IREREpA6d/F0BsVpaWrBjxw4UFhZajmm1WowYMQLl5eV2z7n++uuxYsUKVFRUID09HT/99BNKSkpw//33O7xPc3MzmpubLY8NBgMAwGg0wmg0SvRqrJmv663rqwHbUBpsR2mwHT3HNpSGmttR7GtWTCB05swZtLa2IjY21up4bGwsDh48aPeccePG4cyZM7jxxhshCAIuXryIGTNmOB0aKy4uxty5c22Ob9iwAWFhYZ69iA6UlZV59fpqwDaUBttRGmxHz7ENpaHGdmxqahJVTjGBkDs2b96MF198EX/961+RkZGBI0eO4JFHHsHzzz+PZ555xu45hYWFKCgosDw2GAyIj49HdnY2IiIivFJPo9GIsrIyZGVlITg42Cv3CHRsQ2mwHaXBdvQc21Aaam5H84hORxQTCHXv3h1BQUGoqamxOl5TUwO9Xm/3nGeeeQb3338/pk6dCgC47rrr0NjYiOnTp+Opp56CVmubIqXT6aDT6WyOBwcHe/1N5It7BDq2oTTYjtJgO3qObSgNNbaj2NermGTpkJAQDB06FBs3brQcM5lM2LhxIzIzM+2e09TUZBPsBAUFAQAEQfBeZYmIiEgRFNMjBAAFBQWYOHEi0tLSkJ6ejoULF6KxsRGTJ08GADzwwAPo1asXiouLAQB5eXl47bXXMHjwYMvQ2DPPPIO8vDxLQERERETqpahAaMyYMTh9+jTmzJmD6upqpKamorS01JJAfeLECaseoKeffhoajQZPP/00Tp48iR49eiAvLw/z5s3z10sgIiIiGVFUIAQA+fn5yM/Pt/vc5s2brR536tQJRUVFKCoq8kHNiIiISGkUkyNEREREJDUGQkRERKRaDISIiIhItRgIERERkWoxECIiIiLVYiBEREREqsVAiIiIiFRLcesIERGROK0mARWVdahtuICY8FCkJ0YjSKvxd7WIZIWBEBFRACrdW4W56/ajqv6C5VhcZCiK8pKQmxznx5oRyQuHxoiIAkzp3irMXLHTKggCgOr6C5i5YidK91b5qWZE8sNAiIgogLSaBMxdtx+CnefMx+au249Wk70SROrDQIiIKIBUVNbZ9AS1JwCoqr+Aiso631WKSMYYCBERBZDaBsdBkDvliAIdAyEiogASEx4qaTmiQMdAiIgogKQnRiMuMhSOJslr0DZ7LD0x2pfVIpItBkJERAEkSKtBUV4SANgEQ+bHRXlJXE+I6L8YCBERBZjc5DgsmTAE+kjr4S99ZCiWTBjCdYSI2uGCikREASg3OQ5ZSXquLE3UAQZCREQBKkirQWa/bv6uBpGscWiMiIiIVIuBEBEREakWAyEiIiJSLQZCREREpFoMhIiIiEi1GAgRERGRajEQIiIiItViIERERESqxUCIiIiIVIuBEBEREakWAyEiIiJSLQZCREREpFoMhIiIiEi1GAgRERGRajEQIiIiItVSXCC0ePFiJCQkIDQ0FBkZGaioqHBa/vz583jwwQcRFxcHnU6Ha665BiUlJT6qLREREclZJ39XwBWrV69GQUEBli5dioyMDCxcuBA5OTk4dOgQYmJibMq3tLQgKysLMTExWLNmDXr16oXjx48jKirK95UnIiIi2VFUIPTaa69h2rRpmDx5MgBg6dKlWL9+PZYtW4bZs2fblF+2bBnq6urw3XffITg4GACQkJDg9B7Nzc1obm62PDYYDAAAo9EIo9Eo0SuxZr6ut66vBmxDabAdpcF2FKfVJGD78XOobWhGTLgOaX27IkirAcA2lIqa21Hsa9YIgiB4uS6SaGlpQVhYGNasWYPRo0dbjk+cOBHnz5/H2rVrbc4ZNWoUoqOjERYWhrVr16JHjx4YN24cnnzySQQFBdm9z7PPPou5c+faHF+5ciXCwsIkez1ERGq256wGnxzT4nyLxnIsKkTA3QkmpHRTxMcSyVxTUxPGjRuH+vp6REREOCynmB6hM2fOoLW1FbGxsVbHY2NjcfDgQbvn/PTTT9i0aRPGjx+PkpISHDlyBH/6059gNBpRVFRk95zCwkIUFBRYHhsMBsTHxyM7O9tpQ3rCaDSirKwMWVlZlp4rcg3bUBpsR2mwHZ37cl8NlpfvweXhTn2LBssPB2HR2BTcdk0021ACan4vmkd0OqKYQMgdJpMJMTExePvttxEUFIShQ4fi5MmTePXVVx0GQjqdDjqdzuZ4cHCw199EvrhHoGMbSoPtKA22o61Wk4B5XxyyCYIAQACgATDvi0MYMeAmAGxDqaixHcW+XsUEQt27d0dQUBBqamqsjtfU1ECv19s9Jy4uDsHBwVbDYAMGDEB1dTVaWloQEhLi1ToTEZG1iso6VNVfcPi8AKCq/gK2Hz/nu0qRqilm+nxISAiGDh2KjRs3Wo6ZTCZs3LgRmZmZds+54YYbcOTIEZhMJsuxw4cPIy4ujkEQEZEf1DY4DoKsyzV3XIhIAooJhACgoKAA77zzDt5//30cOHAAM2fORGNjo2UW2QMPPIDCwkJL+ZkzZ6Kurg6PPPIIDh8+jPXr1+PFF1/Egw8+6K+XQESkajHhoSLL2aYoEHmDYobGAGDMmDE4ffo05syZg+rqaqSmpqK0tNSSQH3ixAlotZdiu/j4eHz55ZeYNWsWBg0ahF69euGRRx7Bk08+6a+XQESkaumJ0YiLDEV1/QW7eUIaAPrIUKT17YovD/i6dqRGigqEACA/Px/5+fl2n9u8ebPNsczMTGzbts3LtSIiIjGCtBoU5SVh5oqd0ABWwZB5In1RXpJlPSEib1PU0BgRESlfbnIclkwYAn2k9TCZPjIUSyYMQW5ynJ9qRmqkuB4hIiJSvtzkOGQl6VFRWYfahguICQ9FemI0e4LI5xgIERGRXwRpNcjs183f1SCV49AYERERqRYDISIiIlItBkJERESkWgyEiIiISLUYCBEREZFqMRAiIiIi1WIgRERERKrFQIiIiIhUi4EQERERqRYDISIiIlItBkJERESkWgyEiIiISLUYCBEREZFqMRAiIiIi1WIgRERERKrFQIiIiIhUq5O/K0BERORMq0lARWUdahsuICY8FOmJ0QjSavxdLRtKqSdZYyBERESy9eW+Gsz74hCq6i9YjsVFhqIoLwm5yXF+rJm10r1VmLtuv+zrSbY4NEZERLK056wGD63aYxVcAEB1/QXMXLETpXur/FQza6V7qzBzxU7Z15PsYyBERESy02oS8MkxLQQ7z5mPzV23H60meyV8p9UkYO66/bKvJznGQIiIiGRn+/FzON/iOL9GAFBVfwEVlXW+q5QdFZV1Nj1B7cmlnuQYAyEiIpKd2oZmkeUcByG+IPb+/q4nOcZAiIiIZCcmXCeyXKiXayLN/f1dT3KMgRAREclOWt+uiAoR4GhwTIO2WVnpidG+rJaN9MRoxEWGyr6e5BgDISIikp0grQZ3J5gAwCbIMD8uykvy+zo9QVoNivKSAMi7nuQYAyEiIpKllG4CFo1NgT7SelhJHxmKJROGyGZ9ntzkOCyZMET29ST7uKAikYJxJVsKdDkDYzFyUC/Zv89zk+OQlaSXfT3JFgMhIoXiSrakFkFaDTL7dfN3NTqklHqSNQ6NESkQV7IlIpIGAyEiheFKtkRE0lFcILR48WIkJCQgNDQUGRkZqKioEHXeqlWroNFoMHr0aO9WkMjLuJItEZF0FBUIrV69GgUFBSgqKsLOnTuRkpKCnJwc1NbWOj3v2LFj+J//+R/cdNNNPqopkfdwJVsiIukoKhB67bXXMG3aNEyePBlJSUlYunQpwsLCsGzZMofntLa2Yvz48Zg7dy6uvPJKH9aWyDu4ki35QqtJQPnRs1i7+yTKj57lUCsFLMXMGmtpacGOHTtQWFhoOabVajFixAiUl5c7PO+5555DTEwMpkyZgm+++abD+zQ3N6O5+dIeNwaDAQBgNBphNBo9eAWOma/rreurgZracHDvcOgjdKgxNNvNE9IA0EfqMLh3uMvtoaZ29Calt+OX+2rwQslBVBsu/S3UR+jw9Kj+yBkY65M6KL0N5ULN7Sj2NSsmEDpz5gxaW1sRG2v9SxgbG4uDBw/aPefbb7/F3/72N+zevVv0fYqLizF37lyb4xs2bEBYWJhLdXZVWVmZV6+vBmppw1F6DZYZzB267dcpESAAGBnbhC9Lv3D7+mppR29TYjvuOavBssO2761qwwXkr9qNP1xjQko33/UOKbEN5UiN7djU1CSqnGICIVc1NDTg/vvvxzvvvIPu3buLPq+wsBAFBQWWxwaDAfHx8cjOzkZERIQ3qgqj0YiysjJkZWUhODjYK/cIdGprw1EAhtj51h4XGYqnRrr/rV1t7egtSm3HVpOA4gVbANjb+V0DDYAvasLwxPibvb5QoFLbUG7U3I7mEZ2OKCYQ6t69O4KCglBTU2N1vKamBnq93qb80aNHcezYMeTl5VmOmUxt+9Z06tQJhw4dQr9+/WzO0+l00Olsdz0ODg72+pvIF/cIdGpqwztSe3ttxV01taM3Ka0dtx89axVYX65tRmIzdv2nwWcLByqtDeVKje0o9vUqJhAKCQnB0KFDsXHjRssUeJPJhI0bNyI/P9+mfP/+/fHvf//b6tjTTz+NhoYG/OUvf0F8fLwvqk3kVVzJlqTEGYmkRooJhACgoKAAEydORFpaGtLT07Fw4UI0NjZi8uTJAIAHHngAvXr1QnFxMUJDQ5GcnGx1flRUFADYHCci8hYl7QfHGYmkRooKhMaMGYPTp09jzpw5qK6uRmpqKkpLSy0J1CdOnIBWq6gVAYgogH25rwbzvjjks/3gPA260hOjERcZiur6C05mJLZdlyhQKCoQAoD8/Hy7Q2EAsHnzZqfnvvfee9JXiIjIjj1nNVhevscmoDDvB7dkwhBJgyEpNuEN0mpQlJeEmSt2QgNY1d0cThXlJcm2R4vIHew+ISKSWKtJwCfHtD7bD07KTXhzk+OwZMIQ6COth7/0kaGSB29EcqC4HiEiIrnbfvwczrc47jVpvx+cp8nuHW3Cq0Fb0JWVpBfdk5ObHIesJL1icpuIPMFAiIhIYrUNjqegW5fzfPaVK5vwuhJ0XT4j0bzlBgMjCjQMhIiIJBYTbrsWmf1yns++8sWUdynyj4jkijlCREQSS+vbFVEhAhz1l2jQFkhIMfvK21Pepcw/IpIjBkJERBIL0mpwd0LbSvaXB0NSz74yT3n3RtDVUf4RIG3SN5E/MBAiIvKClG4CFo1N8frsK/OUd0D6oMuV/CMipWKOEBGRl+QMjPXafnDtmae8X57Ho/cwj4dbbpAaMBAiIvIiX+0H540p79xyg9SAgRARUYCQOujilhukBswRIiIiu7yZf0QkF6J7hD777DPRF/3d737nVmWIiEhevJV/RCQXogOh0aNHiyqn0WjQ2trqbn2IiEhmuOUGBTLRgZDJZPJmPYiISMZ8lfRN5GvMESIiIiLVcnvWWGNjI77++mucOHECLS0tVs89/PDDHleMiIiIyNvcCoR27dqFUaNGoampCY2NjYiOjsaZM2cQFhaGmJgYBkJERESkCG4Njc2aNQt5eXk4d+4cOnfujG3btuH48eMYOnQo5s+fL3UdiYiIiLzCrUBo9+7deOyxx6DVahEUFITm5mbEx8fjlVdewZ///Gep60hERETkFW4FQsHBwdBq206NiYnBiRMnAACRkZH4+eefpasdEclKq0lA+dGzWLv7JMqPnuWu40SkeG7lCA0ePBg//PADrr76atxyyy2YM2cOzpw5gw8//BDJyclS15GIZKB0b5XNonpxXFSPiBTOrR6hF198EXFxbX/45s2bh65du2LmzJk4ffo03n77bUkrSET+V7q3CjNX7LQKggCguv4CZq7YidK9VX6qGRGRZ9zqEUpLS7P8OyYmBqWlpZJViIh8r9UkoKKyDlXnG/FTvQatJgHB7Z6bu26/3U03BbTtOTV33X5kJem50jARKQ53nydSOdshryCsWbAFz/5uIHKT49oCpMt6gtoTAFTVX0BFZR1XHiYixXErEEpMTIRG4/ib308//eR2hYjId8xDXpf39tQYmjFzxU4smTAEzRfFba9T2+A4WCIikiu3AqFHH33U6rHRaMSuXbtQWlqKxx9/XIp6EZGXiR3ymv/7FFHXiwkPlbJ6REQ+4VYg9Mgjj9g9vnjxYmzfvt2jChGRb4gd8oKmbXZYdf0Fu0GTBoA+sm03ciIipZF009WRI0fi448/lvKSROQlYoeyzvzSjKK8JABtQU975sdFeUlMlCYiRZI0EFqzZg2io/mtkEgJxA5lxYSHIjc5DksmDIE+0vocfWQolkwYwnWEiEix3F5QsX2ytCAIqK6uxunTp/HXv/5VssoRkfekJ0a7NOSVmxyHrCQ9KirrUNtwATHhbc+xJ4iIlMytQGj06NFWj7VaLXr06IFhw4ahf//+UtSLiLwsSKtBUV4SZq7YCQ1gFQw5GvIK0mo4RZ6IAopbgVBRUZHU9SAiPzAPeV2+dYY+UoeivIEc8iKigCc6EDIYDKIvGhER4VZliMj32g95VZ1vxE/7diN/zM0I1YX4u2pERF4nOlk6KioKXbt2FfWfNy1evBgJCQkIDQ1FRkYGKioqHJZ95513cNNNN1nqNWLECKflidTKPOSVNygOV0cKzPshCjCtJgHlR89i7e6TKD96Fq0me5mB6iS6R+irr76y/PvYsWOYPXs2Jk2ahMzMTABAeXk53n//fRQXF0tfy/9avXo1CgoKsHTpUmRkZGDhwoXIycnBoUOHEBMTY1N+8+bNuO+++3D99dcjNDQUL7/8MrKzs7Fv3z706tXLa/UkIiKSiz1nNShesAXVhmbLsbjIUBTlJXH4Gy4EQrfccovl38899xxee+013HfffZZjv/vd73Ddddfh7bffxsSJE6Wt5X+99tprmDZtGiZPngwAWLp0KdavX49ly5Zh9uzZNuU/+ugjq8fvvvsuPv74Y2zcuBEPPPCAV+pIREQkF1/uq8Gyw1oAzVbHq+svWLbRUXsw5FaydHl5OZYuXWpzPC0tDVOnTvW4Uva0tLRgx44dKCwstBzTarUYMWIEysvLRV2jqakJRqPR6VpHzc3NaG6+9IYx50YZjUYYjUY3a++c+breur4asA2lwXaUBtvRc2xDz7WaBDxfctDuc5e20dmHYVd3C8jhcLHvHbcCofj4eLzzzjt45ZVXrI6/++67iI+Pd+eSHTpz5gxaW1sRGxtrdTw2NhYHD9r/QV/uySefRM+ePTFixAiHZYqLizF37lyb4xs2bEBYWJhrlXZRWVmZV6+vBmxDabAdpcF29Bzb0H0/1mtQYwiC7Zrwbdq20WnGm6tLcXVk4OUMNTU1iSrnViD0+uuv45577sEXX3yBjIwMAEBFRQV+/PFH2W6x8dJLL2HVqlXYvHkzQkMdr6hbWFiIgoICy2ODwYD4+HhkZ2d7bTac0WhEWVkZsrKyEBwc7JV7BDq2oTTYjtJgO3qObei5df+qAvb/u8NyVw5MxahBgTc8Jna2u1uB0KhRo3D48GEsWbLE0huTl5eHGTNmeK1HqHv37ggKCkJNTY3V8ZqaGuj1eqfnzp8/Hy+99BL++c9/YtCgQU7L6nQ66HQ6m+PBwcFe/2X0xT0CHdtQGmxHabAdPcc2dF9cVBfR5QKxjcW+JrcCIaBteOzFF19093SXhYSEYOjQodi4caNlZWuTyYSNGzciPz/f4XmvvPIK5s2bhy+//BJpaWk+qi0REZF/pSdGQx+hQ7XhAuwNj12+jY5aiQ6E/vWvfyE5ORlarRb/+te/nJbtqNfFXQUFBZg4cSLS0tKQnp6OhQsXorGx0TKL7IEHHkCvXr0sU/hffvllzJkzBytXrkRCQgKqq6sBAFdccQWuuOIKr9SRiIhIDoK0Gjw9qj/yV+12uo0OAJQfPavaPQRFB0Kpqamorq5GTEwMUlNTodFoIAi2yVUajQatra2SVtJszJgxOH36NObMmYPq6mqkpqaitLTUkkB94sQJaLWX1ohcsmQJWlpa8Pvf/97qOkVFRXj22We9UkciIiK5yBkYiz9cY0JJdZjVOkL6/64jBAA3vrzJaosdta0xJDoQqqysRI8ePSz/9pf8/HyHQ2GbN2+2enzs2DHvV4iIvKrVJHDHeyIPpHQT8MT4m7HrPw1Wv0dl+6sxc8VOXN6lobY1hkQHQn379rX7byIibyndW2WzIazavq0SScG8jY5Zq0nA3HX7bYIgoP0aQ/uRlaQP+C8eovcaa+/999/H+vXrLY+feOIJREVF4frrr8fx48clqxwRqVfp3irMXLHTKggCLn1bLd1b5aeaESlfRWWdze9We21rDF1ARWWd7yrlJ24FQi+++CI6d+4MoG2V6TfffBOvvPIKunfvjlmzZklaQSJSn46+rQJt31a5cSSRe2obHAdB7pRTMremz//888+46qqrAACffvopfv/732P69Om44YYbMGzYMCnrR0ROBGr+jCvfVtt395PvBOp7Ty1iwh0vLOxOOSVzKxC64oorcPbsWfTp0wcbNmywrMQcGhqKX3/9VdIKEpF9gZw/w2+r8hbI7z2zQA/00hOjERcZiur6C3Z7XtW0xpBbgVBWVhamTp2KwYMH4/Dhwxg1ahQAYN++fUhISJCyfkRkhzl/JlBne/DbqnwF+nsPUEegF6TVoCgvCTNX7HS6xlAgBX+OuJUjtHjxYmRmZuL06dP4+OOP0a1bW9f0jh07cN9990laQSKypqT8mVaTgPKjZ7F290mUHz0ruk7mb6uO/gRr0PbBpIZvq55wt/2dXU8p7z13qSlJPzc5DksmDIE+0voLhT4yNCACWrHc6hGKiorCm2++aXPc3q7tRCQtpeTPePKtmt9WPeeNXg2lvPfcpcYp5bnJcchK0gf0MGBH3OoRAoBvvvkGEyZMwPXXX4+TJ08CAD788EN8++23klWOiGwpIX9Gim/V/LbqPm/1aijhvecJtU4pN68xdGdqL2T266aqIAhws0fo448/xv3334/x48dj586daG5uW7a7vr4eL774IkpKSiStJBFdIvf8GSm/VQfCt1VfJ916s1dDqveeXBORAz3QI/vcCoReeOEFLF26FA888ABWrVplOX7DDTfghRdekKxyRGRL7rM9pB4+uXxFXCX5cl8N5n1xyKdJt94cvpLivSfnRGS5f8kg73BraOzQoUO4+eabbY5HRkbi/PnzntaJiJww588AsEkmlkP+DL9Vt9lzVoOHVu3xedKtN9vf0/ee3BORmaSvTm4FQnq9HkeOHLE5/u233+LKK6/0uFJE5Jyc82f4rbpt6OeTY1q/zK7ydvu7+95TwowzuX/JIO9wa2hs2rRpeOSRR7Bs2TJoNBqcOnUK5eXleOyxxzBnzhyp60hEdsg1f0buQ3e+sP34OZxvcfxz8ObsKl+0vzvvPaXMODMHepcP3+llMnxH0nMrEJo9ezZMJhOGDx+OpqYm3HzzzdDpdHj88ccxdepUqetIRA7IMX9G7NR3ACg/elZWQZxUahuaRZaTfnjQV0sPuPreU9KQqVy/ZJB3uDU0ptFo8NRTT6Gurg579+7Ftm3bcPr0aURGRiIxMVHqOhKRwnQ0fAIAN768Cfe9sw2PrNqN+97Zhhtf3uT3HBGxOlqoMCZcJ+o63hoelOPQqdKGTNU+pVxNXOoRam5uxrPPPouysjJLD9Do0aOxfPly3HXXXQgKCuLu80QEwPG36rL91YreokHMrKe0vl0RFSKgvkXjt+FBufVqBMKQqVyn/ZNnXAqE5syZg7feegsjRozAd999h3vvvReTJ0/Gtm3bsGDBAtx7770ICgryVl2JSGEuHz5R+sq9YvfZCtJqcHeCCcsPB/l1ZWw5DZ0qfbVwOU/7J8+4NDT297//HR988AHWrFmDDRs2oLW1FRcvXsSePXswduxYBkFEAUjK/aqUvHKvq7OeUroJWDQ2RVbDU/4mxyE7MeQ+7Z8841KP0H/+8x8MHToUAJCcnAydTodZs2ZBo5FnBE9EHWvf3d8trBPaxzlSfwtWUsLs5VwJ4tL6RAAAcgbGYuSgXhxOaUduQ3YdUXovJnXMpUCotbUVISEhl07u1AlXXHGF5JUiIt+wF+hEhQQhOKEGnToFSZ7Lo4SEWUd5IK4FcRGWx3IanpILJbWJUqb9k/tcCoQEQcCkSZOg07XNiLhw4QJmzJiBLl26WJX75JNPpKshEXmFo3yX8y3AQ6v2IDIsWPJvwXJPmHXWA6aEII6kJzYALttfzUBIoVzKEZo4cSJiYmIQGRmJyMhITJgwAT179rQ8Nv9HRJ5xlJcjVb6Os+5+oG2m0/kmo8Pz3c3lkfPKvR3lgZxrbOH2CyokNrBdtvUYc4UUyqUeoeXLl3urHkT0XyX/OoWn1+5FXeOlQCQuMhS/S4nDZ3uqJMnX6ai7Xyx3cnnkuHKvmDyQ59fvxzO3D8CDK3d1OOvJ1OrtGpOvmHsxO/p9Ya6Qcrm1sjQReUdxyX68taXS5nhV/QW7x93N15EqGdmT/ao6Spj15ZotYvNAunbRyS6II+8y92LOWLHTaTnmCikXAyEimSj5V5XdYMcZd/N1PM1jkSKXx1nCrK/XbHElEfrO1F6KmvVEnstNjsOUGxLwt63HOiwrxxmP5JxbW2wQkbRaTQKeXrvXrXPdydcxd/c7y3fpGhZs+fflzwHey+Xxx5otriZCc/sF9RmRpBdVjsnyysNAiPxGyoX6lK6isg51jS0eXcOVb6LOkpbN2S/Fd1+HpT5e/M7VRQulIiYwZCK0uvE9Erg4NEZ+weXqrUnRne7qN1FHSctRIcALd6dYfg6+HAby15otSt/+gbyP75HAxUCIfE7sfk3+5OvNFT3tTnf3m+jlScvdwjrh9P5tyBkYaynjy8Xv/LnytBxns5G88D0SmBgI+cHlH7KDe4f7u0o+o4Tl6v3RWyV2iq4jz9w+wO32ah/oGI1GlBxw6zKS8PeihUrb/oF8j++RwMNAyMfsfcjqI3QYpddglB/r5YjUPSNyX67eX71V7bvd3cl+6dpFJ3md/EEOK08rafsHMXzdu6kGgfYeUTsGQj7k6EO2xtCMZQYthuyrwR2pvf1SN3u80TMi5003/d1b5ajbXYxAmbLLPAxxxAY3zMUj6pjiZo0tXrwYCQkJCA0NRUZGBioqKpyW//vf/47+/fsjNDQU1113HUpKSnxUU2tiZsPM++KgbGZOeWsKs7+HPpxxpbfKW3KT4/Dtk7fhf6f9Fn8Zm4pnbh8g6rxAmrJrDgh9OVtNSUr3VuHGlzfhvne24ZFVu3HfO9tw48ubbH4n/bEMAZESKapHaPXq1SgoKMDSpUuRkZGBhQsXIicnB4cOHUJMTIxN+e+++w733XcfiouLcccdd2DlypUYPXo0du7cieTkZJ/WveMtDTSoqm+Wxaqk3uwZkcPQhyP+7K1y9A2/1STg3W8rZdlernB1eIZ5GPaJHbr1d+8mkRhyGbZVVCD02muvYdq0aZg8eTIAYOnSpVi/fj2WLVuG2bNn25T/y1/+gtzcXDz++OMAgOeffx5lZWV48803sXTpUp/WXc5DQpfzZh6PnIc+/NVb1dHwhVzbSyx3h2eYh2HNleBG7rl4RHIatlVMINTS0oIdO3agsLDQckyr1WLEiBEoLy+3e055eTkKCgqsjuXk5ODTTz91eJ/m5mY0NzdbHhsMBgBts2mMRse7cXekW5i4pu4aGoRvD9egtqEZMeE6pPXt6vMPuarzjaLLGY0RLl9/+LXdsWhsCl4oOYhqw6W21kfq8NTI/hh+bXeX29pc3pOf0eDe4dBH6FBjaHbS+6LD4N7hHt2nvS/31eChVXscfsNfNDYFOQNjJW8vR6Rox/bEvr5AI3U7AsD3IoOb8iO1qG1odliuPXd/h33BaDTCJABbf6xF3a+tfvt7qHTeeC96yld/F8S+ZsUEQmfOnEFraytiY60bJzY2FgcPHrR7TnV1td3y1dXVDu9TXFyMuXPn2hzfsGEDwsLC3Kh5G5MARIUE4XwLYG8tX0BAWCfgoZXbUd9y6fmoEAF3J5iQ0s13uUM/1WsABHVcbt9ulPxnl9v3eTIJOGrQwGAEIoKBfhGNaD2+AyXH3b4kysrK3D8ZwCi9BssM5tS59j8nAQKA1PAmvLii9L/1FeDJ32STAMzdGfTfPwbWFxL++/+nP9kN47FWaDXeaS9HPG1HwPXXF4ikaEezHWfE/V5u+OZ7RARDVFlPf4fdYRIufx/b/z3ac1aDT44F4fy23ZZj/vh7GCikfC96wpd/F5qamkSVU0wg5CuFhYVWvUgGgwHx8fHIzs5GRIRn35yCE9qiYMB2iEMA0HTR9qde36LB8sNBHkfIrSYB24+fE9XT1GoSsGbBlg57RvLH3Cybb2dGoxFlZWXIyspCcHCw29cZBWDIvhqb3peosGBA0OCL/1z6hqGP0OHpUf3d/rl8X1mH89u2OymhwfkWoEfSb5HhYg6QKz/v9qRqR8C7r0/upGxHs26VdfjgR2ft2Sb7pgyk9e0qy9/hL/fVoPjynk07v0df7qvB8vI9//36cYlUfw/VxBvvRU/48u+CeUSnI4oJhLp3746goCDU1NRYHa+pqYFeb38zPL1e71J5ANDpdNDpbNdkCQ4O9vhNdEdqb3TqFGQzLhobEQJDUzOaLtqeYx77n/fFIYwc1MutP1qujsUGA3j2dwM7yEsZiFBdiMt1EcvdJDqpfk4jB/Wy3P/YmSYs/Odhu8sePLRqj9szmc7a+4E7KOfKa5Ji7F2KdvTW61MSKdrRLPOqGFETDTKvikGQVuP33+HLle6tsjsccvnvUatJwLwvDjnsMfD076FaSfle9IQv/y6IPV8x0+dDQkIwdOhQbNy40XLMZDJh48aNyMzMtHtOZmamVXmgrXvQUXlfuHx69P9O+y1evvs6u71BZp5M23Z3Cq0/pzCLnR7sTeZE3TsG9cSqH054ZRNQbyRny2nKtJyXSlAiZxvl2kucl9MyBK5spiuHZSzIe+T4d0ExPUIAUFBQgIkTJyItLQ3p6elYuHAhGhsbLbPIHnjgAfTq1QvFxcUAgEceeQS33HILFixYgNtvvx2rVq3C9u3b8fbbb/vzZdjMhvlkxwlR57k6o8zTKbS+mMJ8ec/PucYWPLhSPvuQeXP2jdRLCchtyrScl0pQKlf3upLLMgSu/B4paYYtuU6OfxcUFQiNGTMGp0+fxpw5c1BdXY3U1FSUlpZaEqJPnDgBrfZSJ9f111+PlStX4umnn8af//xnXH311fj00099voZQR2LCxW2P4GqELMWHuDenMNsbwtFqIJsPcsC7yx5IvZSA2J/3e1sr0T1c5/KHoqvDlXJeKkHJXA1u5LAMgSu/R3LsMSDpyPHvgqICIQDIz89Hfn6+3ec2b95sc+zee+/Fvffe6+VaeSatb1dEhQiob9FIGiHL+ZuVo4XhnI0wSbH2ifnDvNpwAXW/NCO6Swj0kZ0dfpB4+4+ylLtZi/05Pr/+0q6qYnOH3M074m7d3iGH4MYVrvweybHHgKQlt78LiguEAlGQVoO7E0xYfjhI0ghZrt+snA3hiOFu4Gbvw9zM0Ye6L/4oSzV84c7PUcyQo6cb0cpleIb8x5Xfo/Y9Bpf6gi+VA9iTGAjk9HdBMcnSgS6lm4BFY1MkTWw0//Fx9LbSoC0A8PU3q463G3HOnQ98R0nEZlUOkoldTVB1l/kb/p2pvZDZr5tb1+vo521PRwnfriS5OiPF6yPlcifRe9HYFERdNqmN+80FFrn8XWCPkIzkDIy1mrbtaYQsx7FYwP0eHXd7X8T2QAmwn4Mkt25cR5z9vJ1pP+SY1sd6rSxu1UBScfX3KGdgLIzHWtEj6bc423SRPYnkNQyEZEbqsX85foi706PjSeDmSg+Uow91OXXjOuPo5y1GW4AaYeeY2HPlTS4bPKqZq79HWg2QkRgti/VvKHAxEFIBuX2Id5QvALT9AWw/2uJJ4Obqh7Sj8kpJUM1NjsNt/WPxYfkxHK9rgiAI+HBbx0s02AtQ5Zpn5io5bfCodkr5PSL1YCCkEnL64yNmyO7N+wajaxedJIGbqx/Scv9Q74ijZQkcpfG0H3I0tVqv+hoIM3g8TfYOdOwpI7VjIER+4cshOzE9UGb+SB6XkqvLElw+5GhqtX5ernlmYsltkUm5YU8ZEQMh8iNfDdlZT8d1TAN5f6h3RExSuDtDjv7KM5Oip0Jssve2n87ihqu6e1hjZVFqTxl7sEhqDITIr3w1ZNdREnEgfAsWkxRuEoBnbh/g8srSvs4zk6qnQmx+2PQPtuPV3w/CqEE9Xa6rEvmzp8yTQIY9WOQNDIRINdp/mItdWVpJxH7odw/X4c7UXi5f31dBq5Q9FWLzvRpbWvGnlbvwx/+cR+GoJNF1VWrvhL+WRfAkkFFqDxbJHwMhUhU5JY1LLRBmeEndU+FKfhgAvLWlEim9u2LUoI4/UJ19qA+/Vt7DbP5YFsGTQIa5XuRNXFmaKEDIdSVxV7jSUyFG+xWNxXpm7d4OV8p2tFK5+UP9y301Lt3T13wdNHu6QrnU7wui9hgIEQUId7cDaTUJKD96Fmt3n8T3lXVON771Nm/0VJjzw6I6i1uU72xji9MPVDEf6vO+OOjXduyIr4NmTwOZQFrYk+SHQ2NEAcTVGV72hneiQoIQnFCDO1J7+6zeZt7qqchNjkO4Lhjj//a9qPLOPlDFfag346jBt0M0ruQr+XpZBE8DmUAY9iX5YiBEFGDEzvBylLNxvgV4aNUedOoUJFnyqdgPaW8u4Pjbft0Q3SUYdY3GDss6+0AV+6Fu6Pg2knEnCdmXyyJ4GsgEwsKeJF8MhIgCUEdJ4c7XHGoLUKRKPnXlQ9qbPRVBWg1euDMZf1q5y2m5joaExH6oR/hoeyxPkpB9tSyCp4GM0hf2JHljjhCRCvkq+bSjpOLSvVU255h7KvSR1gGHPjLU4ynSowb1xB9vTnT4vJhFNcXl1+jQL8L7SUKeJiEDl4LmO1N7IbNfN68EE+7mr7XnzfcFqRt7hIhUyBfJp55MeXa1p6LlosmyyWzf6DDcn5mAkE72v+cVjkpCSu+ueHrtXtQ1tliOi13PRkzvxFMj+6P1+A6n15GCv9YDcocUQ3Fy20CaAgMDISIV8kXyqacf0mLXfCou2Y93vqm0mqU1r+QApt2U6HBxxFGD4pCT7P4Hakcf6sOv7Y6S46Iu5RGlzaaSIpAJ5LXAyD8YCBGpkC+ST33xIV1csh9vbam0OW4SYDnuKBjy9APV2Ye60eibTGklzqZiIENywxwhIhVylrNhHuzxNPnU2x/SLRdNeOcb2yCovXe+qUTLRZNb1xfDF/k1zgTCIppE/sZAiEilHCWfRoUAi8am2M3ZaL/4YvnRs06TcL39If1h+bEOFy00CW3lPOHKa/Y1KZKQidSOQ2NEKnb58E63sE44vX8bcgbG2pR1da0ab095Pl7XJGk5e5Sw27kv1wMiCkQMhIhUrn3OhtFoRMkB2zLurlXjzQ/pvtFhkpa7nJJ2O+dsKiL3MRAiIqc83fnbWx/S92cmYF7JAafDY1pNWzlXKXG3cyYhE7mHOUJE5JQUiy96I6k4pJMW025yvDgiAEy7KdHhekLOqGG3cznnPhH5EnuEiMgpOa9VY54af/k6QloNnK4j1BE5v2YpKCH3ichXGAgRkVNyX6umcFQSHsvuL3plaTHk/po9oaTcJyJfYCBERE4pYefvkE5aTLnpSsmup4TX7A4l5j4ReRtzhIjIKTWuVePP1+zN3B015D4RuYo9QkTUITmuVdNqErw6Xdwfr9nbuTuBnvtE5A4GQkQkipzWqvFVsq8vX7MvcncCOfeJyF0MhIhINF+sVdNRT4+vk3199Zp9kbsTqLlPRJ5gIEREstFRT0+gJvu6krvjSVDm7W1PiJRIMcnSdXV1GD9+PCIiIhAVFYUpU6bgl19+cVr+oYcewrXXXovOnTujT58+ePjhh1FfX+/DWhORWOaenssDAnNPT+neqoBN9vVl7o6jzXb1kaGcOk+qpJgeofHjx6OqqgplZWUwGo2YPHkypk+fjpUrV9otf+rUKZw6dQrz589HUlISjh8/jhkzZuDUqVNYs2aNj2tPRM6I7el5IudaUddTWrKvr3N35JTvReRvigiEDhw4gNLSUvzwww9IS0sDACxatAijRo3C/Pnz0bNnT5tzkpOT8fHHH1se9+vXD/PmzcOECRNw8eJFdOpk/6U3NzejubnZ8thgMABo24zSaDRK+bIszNf11vXVgG0oDX+14/cie3pON/wq6nrdwjp59BpaTQK2Hz+H2oZmxITrkNa3q0tBgqvtOLh3OPQROtQYmp3k7ugwuHe4pD+btD4RACIAAKbWizC1SnZpj/F3Whpqbkexr1kRgVB5eTmioqIsQRAAjBgxAlqtFt9//z3uuusuUdepr69HRESEwyAIAIqLizF37lyb4xs2bEBYmHu7WItVVlbm1eurAdtQGr5uxx1nNACCOiz3848HEBWixfkWwHaFHwAQEBUCnN6/DSUH3KvLnrMafHJMi/Mtl64fFSLg7gQTUrq5tqaPK+04Sq/BMoM5W6H9axMgABgZ24QvS79w6f6BgL/T0lBjOzY1NYkqp4hAqLq6GjExMVbHOnXqhOjoaFRXV4u6xpkzZ/D8889j+vTpTssVFhaioKDA8thgMCA+Ph7Z2dmIiIhwvfIiGI1GlJWVISsrC8HBwV65R6BjG0rDX+3YrbIOH/y4vcNyOTdnIKPJiIdW7QFgL9lXgxfuTkHOwFi36vHlvhosL99j0ytT36LB8sNBWDRW3LXdacdRAIbsq8ELJQdRbbjUKx0XGYqnRvZ3+zUpFX+npaHmdjSP6HTEr4HQ7Nmz8fLLLzstc+CAm1/r2jEYDLj99tuRlJSEZ5991mlZnU4HnU5nczw4ONjrbyJf3CPQsQ2l4et2zLwqRtS07syrYhCk1aBTpyDJFzpsNQmY98Uhp3lK8744hJGDeokeJnO1He9I7Y2Rg3oxd6cd/k5LQ43tKPb1+jUQeuyxxzBp0iSnZa688kro9XrU1tZaHb948SLq6uqg1+udnt/Q0IDc3FyEh4fjH//4h+reCERK4Oq0bm8k+/pqCntHfLFuERFd4tdAqEePHujRo0eH5TIzM3H+/Hns2LEDQ4cOBQBs2rQJJpMJGRkZDs8zGAzIycmBTqfDZ599htBQrpZKJFeubmkhdcDA7SeI1EkROUIDBgxAbm4upk2bhqVLl8JoNCI/Px9jx461zBg7efIkhg8fjg8++ADp6ekwGAzIzs5GU1MTVqxYAYPBYBkv7NGjB4KCOk7MJCLf8ue0bm4/QaROigiEAOCjjz5Cfn4+hg8fDq1Wi3vuuQdvvPGG5Xmj0YhDhw5ZssR37tyJ77//HgBw1VVXWV2rsrISCQkJPqs7EYnnr6Ehbj/hGW9vgkvkLYoJhKKjox0unggACQkJEIRLf76GDRtm9ZiIyBlvbz8RyIGCrzbBJfIGxQRCRETe5mqekliBHCj4ehNcIqkxECIiakfqPKVADhQCdRNcUhcGQkREl5EqT6mjQAFQdqAglyUHiDyhmN3niYiUZvvxc04DBaAtUHhz048+qpG0uOQABQL2CBEReUltQ3PHhQC8/s8fca0+XHFDZFxywLcCOeHenxgIERF5SUy47XY9jihxiIxLDvhOICfc+xuHxoiIvCStb1fERYrrDTHn0iiJeckB4NISA2ZSLDlAbcwJ95cPs5oT7kv3VvmpZoGBgRARkZe0DxTEUGIujXnJAf1lAZ8+MlTRM+LkQmzCfauJ6+a5i0NjRERelJsch1kjrsHr/zzcYVml5tL4c2uUQMeZed7HQIiIyMvyb7sK/1txHNUG+8nTgZBL46+tUQIdZ+Z5H4fGiIi8LEirwbO/GwgNmEtDruHMPO9jIERE5APMpSF3mGfmOQqRNWibPabk3kR/49AYEZGPMJeGXOXtzYCJgRARkU8xl4Zc5a3NgKkNAyEiIiKZY2+i9zAQIiIiUgD2JnoHk6WJiIhItRgIERERkWpxaIyIyEXcBZwocDAQIiJyAXcBJwosHBojIhLJ0S7gVfUXMGPFTpT865RX7ttqElB+9CzW7j6J8qNnucEmkYTYI0REJIKzXcDN8v93F96EBlkDukt2X/ZAEXkXe4SIiEToaBdwADAJwJ9W7sSX+2okuaejHqjq+guYuWInSvdWSXIfIjVjIEREJIIru3vP++IgPB29ctYDZT42d91+DpMReYiBEBGRCK7s7l1V34yjBs9mkXXUAyWgLTeporLOo/sQqR0DISIiEcy7gItlMHp2P7E9UK70VBGRLQZCREQimHcBFysi2LP7ie2BcqWniohsMRAiIhIpNzkOfx03GM7WTtQAiIvUoV+EZ7k75h4oR7dqu0/bYo5E5D4GQkRELhg1qCfevG+I3efMQctTI/s7DZbEaN8DdfmlzI+L8pK4ojWRhxgIERG5aNSgOCydMMQmZ0gfGYolE4YgZ2CsJPfJTY7DkglDoHdwH64jROQ5LqhIROSG3OQ4ZCXp7e45ZjR6mCkt8j5E5DkGQkREbgrSapDZr5si78ONY32PbS5PDISIiFSG23b4HttcvhSTI1RXV4fx48cjIiICUVFRmDJlCn755RdR5wqCgJEjR0Kj0eDTTz/1bkWJiGSM23b4Httc3hQTCI0fPx779u1DWVkZPv/8c2zZsgXTp08Xde7ChQuh0bD7kYjUjdt2+B7bXP4UEQgdOHAApaWlePfdd5GRkYEbb7wRixYtwqpVq3Dq1Cmn5+7evRsLFizAsmXLfFRbIiJ54rYdvsc2lz9F5AiVl5cjKioKaWlplmMjRoyAVqvF999/j7vuusvueU1NTRg3bhwWL14MvV4v6l7Nzc1obm62PDYYDAAAo9Eo6UyQ9szX9db11YBtKA22ozTk2o5V5xtFlzMaI7xcG+fk2oau8nebB0o7ukPsa1ZEIFRdXY2YmBirY506dUJ0dDSqq6sdnjdr1ixcf/31uPPOO0Xfq7i4GHPnzrU5vmHDBoSFhYmvtBvKysq8en01YBtKg+0oDbm140/1GgBBHZfbtxsl/9nl/QqJILc2dJVc2lzp7eiOpqYmUeX8GgjNnj0bL7/8stMyBw4ccOvan332GTZt2oRdu1x7YxUWFqKgoMDy2GAwID4+HtnZ2YiI8M43JKPRiLKyMmRlZSE42MMNilSKbSgNtqM05NqOrSYBaxZsQY2h2W7OigaAPlKH/DE3+31at1zb0FX+bvNAaUd3mEd0OuLXQOixxx7DpEmTnJa58sorodfrUVtba3X84sWLqKurczjktWnTJhw9ehRRUVFWx++55x7cdNNN2Lx5s93zdDoddDqdzfHg4GCvv4l8cY9AxzaUBttRGnJrx2AAz/5uIGau2AkNYPXBfGnbjoEI1YX4vnIOuNqGclurRy5tLrf3oi+Ifb1+DYR69OiBHj16dFguMzMT58+fx44dOzB06FAAbYGOyWRCRkaG3XNmz56NqVOnWh277rrr8PrrryMvL8/zyhMRKZB5247L17TRB8CaNnJdqyeQ2zwQKCJHaMCAAcjNzcW0adOwdOlSGI1G5OfnY+zYsejZsycA4OTJkxg+fDg++OADpKenQ6/X2+0t6tOnDxITE339EoiIZCMQt+0wr9Vz+fCTea0ef+/NFohtHigUEQgBwEcffYT8/HwMHz4cWq0W99xzD9544w3L80ajEYcOHRKdHEVEpGa+2h7EFzpaq0eDtrV6spL0fg08AqnNA4liAqHo6GisXLnS4fMJCQkQBOcLUnX0PBERKY8ra/UwEKHLKWJBRSIiIkdqGxwHQe6UI3VhIERERIoWEx4qaTlSFwZCRESkaOmJ0YiLDIWj7B8N2maPpSdG+7JapBAMhIiISNGCtBoU5SUBgE0wdGmtniTO0CK7GAgREZHimdfq0UdaD3/pI0P9PnWe5E0xs8aIiIic4Vo95A4GQkREFDC4Vg+5ikNjREREpFoMhIiIiEi1GAgRERGRajFHiIhIhVpNApOKicBAiIhIdUr3VmHuuv1W+3PFRYaiKC+J08xJdTg0RkSkIqV7qzBzxU6bTUqr6y9g5oqdKN1b5aeaWWs1CfixXoN1/6pC+dGzaDVx02zyDvYIERGpRKtJwNx1+2EvpBDQtgrz3HX7kZWk9+swWeneKjz72T5UG4KA/f8GwB4r8h72CBERqURFZZ1NT1B7AoCq+guoqKzzXaUuY+6xqjY0Wx2XW48VBQ4GQkREKlHb4DgIcqec1DrqsQLaeqw4TEZSYiBERKQSMeGhHRdyoZzUlNBjRYGHgRARkUqkJ0YjLjLUZod2Mw3acnHSE6N9WS0LufdYUWBiIEREpBJBWg2K8pIAwCYYMj8uykvyW6K03HusKDAxECIiUpHc5DgsmTAE+kjrYEIfGYolE4b4dVaW3HusKDBx+jwRkcrkJschK0kvu5WlzT1WM1fshAawSpqWQ48VBSYGQkREKhSk1SCzXzd/V8OGuceqbR2hS1Po9VxHiLyEgRAREVnIYQ+y3OQ4DLu6G95cXYorB6YiLqqLLHqsKDAxECIiIgDy2oMsSKvB1ZECRg2KQ3BwsE/vTerCZGkiIlLMHmREUmMgRESkclzRmdSMgRARkcpxRWdSMwZCREQqxxWdSc0YCBERqRxXdCY1YyBERKRyXNGZ1IyBEBGRysl9DzIib2IgREREst6DjMibuKAiEREBkO8eZETepJgeobq6OowfPx4RERGIiorClClT8Msvv3R4Xnl5OW677TZ06dIFERERuPnmm/Hrr7/6oMZERMpj3oPsztReyOzXjUEQBTzFBELjx4/Hvn37UFZWhs8//xxbtmzB9OnTnZ5TXl6O3NxcZGdno6KiAj/88APy8/Oh1SrmZRMREZEXKWJo7MCBAygtLcUPP/yAtLQ0AMCiRYswatQozJ8/Hz179rR73qxZs/Dwww9j9uzZlmPXXnutT+pMRERE8qeIQKi8vBxRUVGWIAgARowYAa1Wi++//x533XWXzTm1tbX4/vvvMX78eFx//fU4evQo+vfvj3nz5uHGG290eK/m5mY0NzdbHhsMBgCA0WiE0WiU8FVdYr6ut66vBmxDabAdpcF29BzbUBpqbkexr1kRgVB1dTViYmKsjnXq1AnR0dGorq62e85PP/0EAHj22Wcxf/58pKam4oMPPsDw4cOxd+9eXH311XbPKy4uxty5c22Ob9iwAWFhYR6+EufKysq8en01YBtKg+0oDbaj59iG0lBjOzY1NYkq59dAaPbs2Xj55Zedljlw4IBb1zaZTACAP/7xj5g8eTIAYPDgwdi4cSOWLVuG4uJiu+cVFhaioKDA8thgMCA+Ph7Z2dmIiIhwqy4dMRqNKCsrQ1ZWFoKDg71yj0DHNpQG21EabEfPsQ2loeZ2NI/odMSvgdBjjz2GSZMmOS1z5ZVXQq/Xo7a21ur4xYsXUVdXB71eb/e8uLi2NS+SkpKsjg8YMAAnTpxweD+dTgedTmdzPDg42OtvIl/cI9CxDaXBdpQG29FzbENpqLEdxb5evwZCPXr0QI8ePTosl5mZifPnz2PHjh0YOnQoAGDTpk0wmUzIyMiwe05CQgJ69uyJQ4cOWR0/fPgwRo4c6XnliYiISPEUMY98wIAByM3NxbRp01BRUYGtW7ciPz8fY8eOtcwYO3nyJPr374+KigoAgEajweOPP4433ngDa9aswZEjR/DMM8/g4MGDmDJlij9fDhEREcmEIpKlAeCjjz5Cfn4+hg8fDq1Wi3vuuQdvvPGG5Xmj0YhDhw5ZJUc9+uijuHDhAmbNmoW6ujqkpKSgrKwM/fr188dLICIiIplRTCAUHR2NlStXOnw+ISEBgiDYHJ89e7bVOkKuMl9TbNKVO4xGI5qammAwGFQ3hisVtqE02I7SYDt6jm0oDTW3o/lz215s0J5iAiF/aWhoAADEx8f7uSZERETkqoaGBkRGRjp8XiN0FCqpnMlkwqlTpxAeHg6Nxjt77pin6P/8889em6If6NiG0mA7SoPt6Dm2oTTU3I6CIKChoQE9e/Z0urUWe4Q6oNVq0bt3b5/cKyIiQnVvVKmxDaXBdpQG29FzbENpqLUdnfUEmSli1hgRERGRNzAQIiIiItViICQDOp0ORUVFdle0JnHYhtJgO0qD7eg5tqE02I4dY7I0ERERqRZ7hIiIiEi1GAgRERGRajEQIiIiItViIERERESqxUBIZg4fPow777wT3bt3R0REBG688UZ89dVX/q6WIq1fvx4ZGRno3LkzunbtitGjR/u7SorU3NyM1NRUaDQa7N6929/VUZRjx45hypQpSExMROfOndGvXz8UFRWhpaXF31WTvcWLFyMhIQGhoaHIyMhARUWFv6ukGMXFxfjNb36D8PBwxMTEYPTo0Th06JC/qyVbDIRk5o477sDFixexadMm7NixAykpKbjjjjtQXV3t76opyscff4z7778fkydPxp49e7B161aMGzfO39VSpCeeeAI9e/b0dzUU6eDBgzCZTHjrrbewb98+vP7661i6dCn+/Oc/+7tqsrZ69WoUFBSgqKgIO3fuREpKCnJyclBbW+vvqinC119/jQcffBDbtm1DWVkZjEYjsrOz0djY6O+qyZNAsnH69GkBgLBlyxbLMYPBIAAQysrK/FgzZTEajUKvXr2Ed999199VUbySkhKhf//+wr59+wQAwq5du/xdJcV75ZVXhMTERH9XQ9bS09OFBx980PK4tbVV6Nmzp1BcXOzHWilXbW2tAED4+uuv/V0VWWKPkIx069YN1157LT744AM0Njbi4sWLeOuttxATE4OhQ4f6u3qKsXPnTpw8eRJarRaDBw9GXFwcRo4cib179/q7aopSU1ODadOm4cMPP0RYWJi/qxMw6uvrER0d7e9qyFZLSwt27NiBESNGWI5ptVqMGDEC5eXlfqyZctXX1wMA33cOMBCSEY1Gg3/+85/YtWsXwsPDERoaitdeew2lpaXo2rWrv6unGD/99BMA4Nlnn8XTTz+Nzz//HF27dsWwYcNQV1fn59opgyAImDRpEmbMmIG0tDR/VydgHDlyBIsWLcIf//hHf1dFts6cOYPW1lbExsZaHY+NjWWKgBtMJhMeffRR3HDDDUhOTvZ3dWSJgZAPzJ49GxqNxul/Bw8ehCAIePDBBxETE4NvvvkGFRUVGD16NPLy8lBVVeXvl+F3YtvRZDIBAJ566incc889GDp0KJYvXw6NRoO///3vfn4V/iW2DRctWoSGhgYUFhb6u8qyJLYd2zt58iRyc3Nx7733Ytq0aX6qOanNgw8+iL1792LVqlX+ropscYsNHzh9+jTOnj3rtMyVV16Jb775BtnZ2Th37hwiIiIsz1199dWYMmUKZs+e7e2qyprYdty6dStuu+02fPPNN7jxxhstz2VkZGDEiBGYN2+et6sqW2Lb8P/9v/+HdevWQaPRWI63trYiKCgI48ePx/vvv+/tqsqa2HYMCQkBAJw6dQrDhg3Db3/7W7z33nvQavkd1JGWlhaEhYVhzZo1VjM9J06ciPPnz2Pt2rX+q5zC5OfnY+3atdiyZQsSExP9XR3Z6uTvCqhBjx490KNHjw7LNTU1AYDNH0mtVmvp5VAzse04dOhQ6HQ6HDp0yBIIGY1GHDt2DH379vV2NWVNbBu+8cYbeOGFFyyPT506hZycHKxevRoZGRnerKIiiG1HoK0n6NZbb7X0TDIIci4kJARDhw7Fxo0bLYGQyWTCxo0bkZ+f79/KKYQgCHjooYfwj3/8A5s3b2YQ1AEGQjKSmZmJrl27YuLEiZgzZw46d+6Md955B5WVlbj99tv9XT3FiIiIwIwZM1BUVIT4+Hj07dsXr776KgDg3nvv9XPtlKFPnz5Wj6+44goAQL9+/dC7d29/VEmRTp48iWHDhqFv376YP38+Tp8+bXlOr9f7sWbyVlBQgIkTJyItLQ3p6elYuHAhGhsbMXnyZH9XTREefPBBrFy5EmvXrkV4eLgltyoyMhKdO3f2c+3kh4GQjHTv3h2lpaV46qmncNttt8FoNGLgwIFYu3YtUlJS/F09RXn11VfRqVMn3H///fj111+RkZGBTZs2MemcfKqsrAxHjhzBkSNHbAJIZiU4NmbMGJw+fRpz5sxBdXU1UlNTUVpaapNATfYtWbIEADBs2DCr48uXL8ekSZN8XyGZY44QERERqRYHq4mIiEi1GAgRERGRajEQIiIiItViIERERESqxUCIiIiIVIuBEBEREakWAyEiIiJSLQZCREREpFoMhIhIEhqNBp9++qm/qyELmzdvhkajwfnz5wEA7733HqKiovxaJyKyj4EQEYkyadIkq93AL1dVVYWRI0f6rkIKMmbMGBw+fFh0+YSEBCxcuNB7FSIiC+41RkSSkMMmooIgoLW1FZ06ef6nTcprde7c2S+bXba0tCAkJMTn9yVSEvYIEZEk2g+NHTt2DBqNBp988gluvfVWhIWFISUlBeXl5VbnfPvtt7jpppvQuXNnxMfH4+GHH0ZjY6Pl+Q8//BBpaWkIDw+HXq/HuHHjUFtba3nePAT1xRdfYOjQodDpdPj2229t6mauz6pVq3D99dcjNDQUycnJ+Prrrzu8lslkQnFxMRITE9G5c2ekpKRgzZo1VtcvKSnBNddcg86dO+PWW2/FsWPHrJ63NzS2bt06/OY3v0FoaCi6d++Ou+66C0DbRpnHjx/HrFmzoNFooNFoLOd8/PHHGDhwIHQ6HRISErBgwQKrayYkJOD555/HAw88gIiICEyfPt3BT4uILAQiIhEmTpwo3HnnnQ6fByD84x//EARBECorKwUAQv/+/YXPP/9cOHTokPD73/9e6Nu3r2A0GgVBEIQjR44IXbp0EV5//XXh8OHDwtatW4XBgwcLkyZNslzzb3/7m1BSUiIcPXpUKC8vFzIzM4WRI0danv/qq68EAMKgQYOEDRs2CEeOHBHOnj1rUzdzfXr37i2sWbNG2L9/vzB16lQhPDxcOHPmjNNrvfDCC0L//v2F0tJS4ejRo8Ly5csFnU4nbN68WRAEQThx4oSg0+mEgoIC4eDBg8KKFSuE2NhYAYBw7tw5QRAEYfny5UJkZKSlPp9//rkQFBQkzJkzR9i/f7+we/du4cUXXxQEQRDOnj0r9O7dW3juueeEqqoqoaqqShAEQdi+fbug1WqF5557Tjh06JCwfPlyoXPnzsLy5cst1+3bt68QEREhzJ8/Xzhy5Ihw5MgRcT9cIhVjIEREorgTCL377ruW5/ft2ycAEA4cOCAIgiBMmTJFmD59utU1vvnmG0Gr1Qq//vqr3Xv88MMPAgChoaFBEIRLwcunn37qtO7m+rz00kuWY0ajUejdu7fw8ssvO7zWhQsXhLCwMOG7776zut6UKVOE++67TxAEQSgsLBSSkpKsnn/yySedBkKZmZnC+PHjHda3b9++wuuvv251bNy4cUJWVpbVsccff9zq3n379hVGjx7t8LpEZItDY0TkNYMGDbL8Oy4uDgAsQ1t79uzBe++9hyuuuMLyX05ODkwmEyorKwEAO3bsQF5eHvr06YPw8HDccsstAIATJ05Y3SctLU1UfTIzMy3/7tSpE9LS0nDgwAGH1zpy5AiampqQlZVlVc8PPvgAR48eBQAcOHAAGRkZDu9jz+7duzF8+HBRdTY7cOAAbrjhBqtjN9xwA3788Ue0trbarT8RdYzJ0kTkNcHBwZZ/m3NdTCYTAOCXX37BH//4Rzz88MM25/Xp0weNjY3IyclBTk4OPvroI/To0QMnTpxATk4OWlparMp36dJFsjq3v9Yvv/wCAFi/fj169eplVU6n07l9D28mTkvZFkRqwECIiPxiyJAh2L9/P6666iq7z//73//G2bNn8dJLLyE+Ph4AsH37do/uuW3bNtx8880AgIsXL2LHjh3Iz893WD4pKQk6nQ4nTpyw9EZdbsCAAfjss89s7uPMoEGDsHHjRkyePNnu8yEhIVa9POb7bN261erY1q1bcc011yAoKMjp/YjIMQZCRCRafX09du/ebXWsW7dulkDFFU8++SR++9vfIj8/H1OnTkWXLl2wf/9+lJWV4c0330SfPn0QEhKCRYsWYcaMGdi7dy+ef/55j+q/ePFiXH311RgwYABef/11nDt3Dn/4wx8clg8PD8f//M//YNasWTCZTLjxxhtRX1+PrVu3IiIiAhMnTsSMGTOwYMECPP7445g6dSp27NiB9957z2k9ioqKMHz4cPTr1w9jx47FxYsXUVJSgieffBJA2+yvLVu2YOzYsdDpdOjevTsee+wx/OY3v8Hzzz+PMWPGoLy8HG+++Sb++te/etQmRKrn7yQlIlKGiRMnCgBs/psyZYogCPaTpXft2mU5/9y5cwIA4auvvrIcq6ioELKysoQrrrhC6NKlizBo0CBh3rx5ludXrlwpJCQkCDqdTsjMzBQ+++wzq+uaE5zNScmOmOuzcuVKIT09XQgJCRGSkpKETZs2Wco4upbJZBIWLlwoXHvttUJwcLDQo0cPIScnR/j6668tZdatWydcddVVgk6nE2666SZh2bJlTpOlBUEQPv74YyE1NVUICQkRunfvLtx9992W58rLy4VBgwYJOp1OaP9nes2aNUJSUpIQHBws9OnTR3j11VetrmkvyZqInNMIgiD4IwAjIvKVY8eOITExEbt27UJqaqq/q0NEMsJZY0RERKRaDISIiIhItTg0RkRERKrFHiEiIiJSLQZCREREpFoMhIiIiEi1GAgRERGRajEQIiIiItViIERERESqxUCIiIiIVIuBEBEREanW/wcqGQZbG1GyNgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.clf()\n", "plt.grid(True)\n", "plt.plot(result1.predict(linear=True), result1.resid_pearson, \"o\")\n", "plt.xlabel(\"Linear predictor\")\n", "plt.ylabel(\"Residual\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An alternative variance function is mu^2 * (1 - mu)^2." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:12.532112Z", "iopub.status.busy": "2022-11-02T17:08:12.531694Z", "iopub.status.idle": "2022-11-02T17:08:12.537871Z", "shell.execute_reply": "2022-11-02T17:08:12.537184Z" }, "lines_to_next_cell": 1 }, "outputs": [], "source": [ "class vf(sm.families.varfuncs.VarianceFunction):\n", " def __call__(self, mu):\n", " return mu ** 2 * (1 - mu) ** 2\n", "\n", " def deriv(self, mu):\n", " return 2 * mu - 6 * mu ** 2 + 4 * mu ** 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the quasi-binomial regression with the alternative variance\n", "function." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:12.542501Z", "iopub.status.busy": "2022-11-02T17:08:12.541183Z", "iopub.status.idle": "2022-11-02T17:08:12.602534Z", "shell.execute_reply": "2022-11-02T17:08:12.601782Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: blotch No. Observations: 90\n", "Model: GLM Df Residuals: 72\n", "Model Family: Binomial Df Model: 17\n", "Link Function: Logit Scale: 0.98855\n", "Method: IRLS Log-Likelihood: -21.335\n", "Date: Wed, 02 Nov 2022 Deviance: 7.2134\n", "Time: 17:08:12 Pearson chi2: 71.2\n", "No. Iterations: 25 Pseudo R-squ. (CS): 0.3115\n", "Covariance Type: nonrobust \n", "==================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "----------------------------------------------------------------------------------\n", "C(variety)[1] -7.9224 0.445 -17.817 0.000 -8.794 -7.051\n", "C(variety)[2] -8.3897 0.445 -18.868 0.000 -9.261 -7.518\n", "C(variety)[3] -7.8436 0.445 -17.640 0.000 -8.715 -6.972\n", "C(variety)[4] -6.9683 0.445 -15.672 0.000 -7.840 -6.097\n", "C(variety)[5] -6.5697 0.445 -14.775 0.000 -7.441 -5.698\n", "C(variety)[6] -6.5938 0.445 -14.829 0.000 -7.465 -5.722\n", "C(variety)[7] -5.5823 0.445 -12.555 0.000 -6.454 -4.711\n", "C(variety)[8] -4.6598 0.445 -10.480 0.000 -5.531 -3.788\n", "C(variety)[9] -4.7869 0.445 -10.766 0.000 -5.658 -3.915\n", "C(variety)[10] -4.0351 0.445 -9.075 0.000 -4.907 -3.164\n", "C(site)[T.2] 1.3831 0.445 3.111 0.002 0.512 2.255\n", "C(site)[T.3] 3.8601 0.445 8.681 0.000 2.989 4.732\n", "C(site)[T.4] 3.5570 0.445 8.000 0.000 2.686 4.428\n", "C(site)[T.5] 4.1079 0.445 9.239 0.000 3.236 4.979\n", "C(site)[T.6] 4.3054 0.445 9.683 0.000 3.434 5.177\n", "C(site)[T.7] 4.9181 0.445 11.061 0.000 4.047 5.790\n", "C(site)[T.8] 5.6949 0.445 12.808 0.000 4.823 6.566\n", "C(site)[T.9] 7.0676 0.445 15.895 0.000 6.196 7.939\n", "==================================================================================\n" ] } ], "source": [ "bin = sm.families.Binomial()\n", "bin.variance = vf()\n", "model2 = sm.GLM.from_formula(\"blotch ~ 0 + C(variety) + C(site)\", family=bin, data=df)\n", "result2 = model2.fit(scale=\"X2\")\n", "print(result2.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the alternative variance function, the mean/variance relationship\n", "seems to capture the data well, and the estimated scale parameter is\n", "close to 1." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2022-11-02T17:08:12.606322Z", "iopub.status.busy": "2022-11-02T17:08:12.605937Z", "iopub.status.idle": "2022-11-02T17:08:12.804158Z", "shell.execute_reply": "2022-11-02T17:08:12.803463Z" } }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Residual')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGwCAYAAABRgJRuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBBElEQVR4nO3deXzU1b3/8fckJpMEkiBLSJQAKaIYI6AgGLUWl7BVlNp6VUoLKlq5oBa8Llz9sdQFbWvxXrVosYJKufoQF0QxSl3ABaQsQQMohgaxkrAESTBIGDLf3x/pBIZkMpPJzHyXeT0fDx8yM9/MfObM5Pv95JzPOcdlGIYhAAAAm0swOwAAAIBIIKkBAACOQFIDAAAcgaQGAAA4AkkNAABwBJIaAADgCCQ1AADAEU4wO4BY8nq92rlzp9LT0+VyucwOBwAAhMAwDB04cEAnnXSSEhIC98fEVVKzc+dO5ebmmh0GAAAIwzfffKNu3boFfDyukpr09HRJDY2SkZFhcjRHeTwevfPOOxo6dKiSkpLMDsdWaLvw0Xbho+3CR9uFJ97braamRrm5uY3X8UDiKqnxDTllZGRYLqlJS0tTRkZGXH5Z24K2Cx9tFz7aLny0XXhotwbBSkcoFAYAAI5AUgMAAByBpAYAADgCSQ0AAHAEkhoAAOAIJDUAAMARSGoAAIAjkNQAAABHIKkBAACOEFcrCgOwnnqvoTXl+7T7wCFlpadoUF5HJSaw4SyA1iOpAWCa4tIKzVq6WRXVhxrvy8lM0YxR+RpekGNiZADsyDbDT3PnzlXfvn0b920qLCzUW2+9ZXZYAMJUXFqhiQvX+yU0klRZfUgTF65XcWmFSZEBsCvbJDXdunXTQw89pHXr1mnt2rW6+OKLdcUVV2jTpk1mhwagleq9hmYt3Syjmcd8981auln13uaOAIDm2SapGTVqlEaOHKnevXvr1FNP1QMPPKD27dtr9erVZocGoJXWlO9r0kNzLENSRfUhrSnfF7ugANieLWtq6uvr9dJLL6m2tlaFhYUBj6urq1NdXV3j7ZqaGkkNW7h7PJ6oxxkqXyxWiskuaLvwmdl2FftrQz7O48mIcjStx/cufLRdeOK93UJ93y7DMGzTv/v555+rsLBQhw4dUvv27bVo0SKNHDky4PEzZ87UrFmzmty/aNEipaWlRTNUAC34qtqlxzcnBj1ucn69emfa5hQFIEoOHjyoMWPGqLq6WhkZgf/QsVVSc/jwYe3YsUPV1dVavHixnn76aa1YsUL5+fnNHt9cT01ubq727t3bYqPEmsfj0fLly1VUVKSkpCSzw7EV2i58ZrZdvdfQkEdWaldNXbN1NS5J2ZluvT/1QktO7+Z7Fz7aLjzx3m41NTXq3Llz0KTGVsNPycnJOuWUUyRJAwYM0D/+8Q/9z//8j5566qlmj3e73XK73U3uT0pKsuSXwqpx2QFtFz4z2i5J0szLz9DEhevlkvwSG18KM2PUGUpxJ8c0rtbiexc+2i488dpuob5n2xQKN8fr9fr1xACwj+EFOZo79mxlZ6b43Z+dmaK5Y89mnRoArWabnppp06ZpxIgR6t69uw4cOKBFixbpgw8+0Ntvv212aADCNLwgR0X52awoDCAibJPU7N69W7/+9a9VUVGhzMxM9e3bV2+//baKiorMDg1AGyQmuFTYq5PZYQBwANskNX/961/NDgEAAFiYrWtqAAAAfEhqAACAI5DUAAAARyCpAQAAjkBSAwAAHIGkBgAAOAJJDQAAcASSGgAA4AgkNQAAwBFIagAAgCOQ1AAAAEcgqQEAAI5AUgMAAByBpAYAADgCSQ0AAHAEkhoAAOAIJDUAAMARSGoAAIAjkNQAAABHIKkBAACOQFIDAAAcgaQGAAA4AkkNAABwBJIaAADgCCeYHQAAOF2919Ca8n3afeCQstJTNCivoxITXGaHBTgOSQ0ARFFxaYVmLd2siupDjfflZKZoxqh8DS/IMTEywHkYfgKAKCkurdDEhev9EhpJqqw+pIkL16u4tMKkyABnIqkBgCio9xqatXSzjGYe8903a+lm1XubOwJAOEhqACAK1pTva9JDcyxDUkX1Ia0p3xe7oACHI6kBgCjYfSBwQhPOcQCCI6kBgCjISk+J6HEAgiOpAYAoGJTXUTmZKQo0cdulhllQg/I6xjIswNFIagAgChITXJoxKl+SmiQ2vtszRuWzXg0QQSQ1ABAlwwtyNHfs2crO9B9iys5M0dyxZ7NODRBhLL4HAFE0vCBHRfnZrCgMxABJDQBEWWKCS4W9OpkdBuB4DD8BAABHIKkBAACOQFIDAAAcgaQGAAA4AkkNAABwBJIaAADgCCQ1AADAEUhqAACAI7D4HgA4VL3XYCVjxBWSGgBwoOLSCs1aulkV1Yca78vJTNGMUfnsOQXHYvgJABymuLRCExeu90toJKmy+pAmLlyv4tIKkyIDooukBgAcpN5raNbSzTKaecx336ylm1Xvbe4IwN5IagDAQdaU72vSQ3MsQ1JF9SGtKd8Xu6CAGCGpAQAH2X0gcEITznGAnZDUAICDZKWnRPQ4wE5IagDAQQbldVROZooCTdx2qWEW1KC8jrEMC4gJkhoAcJDEBJdmjMqXpCaJje/2jFH5rFcDRyKpAQCHGV6Qo7ljz1Z2pv8QU3ZmiuaOPZt1auBYLL4HAA40vCBHRfnZrCiMuEJSAwAOlZjgUmGvTmaHAcQMw08AAMARSGoAAIAjMPwExIHmdmsGAKchqQEcLtBuzfeMOM3EqAAg8hh+Ahyspd2ab3lhozZWMRMGgHOQ1AAOFcpuza9sT2C3ZgCOQVIDOFQouzXvP+zS2q+/i11QABBFJDWAQ4W+W3NdlCMBgNiwTVIze/ZsnXPOOUpPT1dWVpZGjx6tL7/80uywAMsKfbdmd5QjAYDYsE1Ss2LFCk2aNEmrV6/W8uXL5fF4NHToUNXW1podGmBJoezW3CHZ0MAeJ8YyLACIGttM6S4uLva7vWDBAmVlZWndunW68MILm/2Zuro61dUd7VqvqamRJHk8Hnk8nugF20q+WKwUk13Qdi27Z8RpuuWFjXJJfgXDvkTnyp5eeeuPyONhFlRr8L0LH20Xnnhvt1Dft8swDFtOfSgrK1Pv3r31+eefq6CgoNljZs6cqVmzZjW5f9GiRUpLS4t2iIAlbKxy6ZXtCdp/+Gji0iHZ0JU9verXyZa//gDizMGDBzVmzBhVV1crIyMj4HG2TGq8Xq8uv/xy7d+/Xx999FHA45rrqcnNzdXevXtbbJRY83g8Wr58uYqKipSUlGR2OLZC24Wm3mto7dffafeBOmWluzWwx4ny1h+h7cLE9y58tF144r3dampq1Llz56BJjW2Gn441adIklZaWtpjQSJLb7Zbb3bQIMikpyZJfCqvGZQe0XcuSJF1wale/+3xDTrRd+Gi78NF24YnXdgv1PdsuqZk8ebLeeOMNrVy5Ut26dTM7HAAAYBG2SWoMw9Att9yiV199VR988IHy8vLMDgkAAFiIbZKaSZMmadGiRVqyZInS09NVWVkpScrMzFRqaqrJ0QEAALPZZp2auXPnqrq6WkOGDFFOTk7jfy+++KLZoQEAAAuwTU+NDSdpAQCAGLJNTw0AAEBLSGoAAIAjkNQAAABHIKkBAACOQFIDAAAcgaQGAAA4AkkNAABwBJIaAADgCCQ1AADAEUhqAACAI5DUAAAARyCpAQAAjkBSAwAAHIGkBgAAOAJJDQAAcASSGgAA4AgnmB0AAOuq9xpaU75Puw8cUlZ6igbldVRigsvssACgWSQ1AJpVXFqhWUs3q6L6UON9OZkpmjEqX8MLckyMDD4knYA/khoAfuq9hh5/r0xz/r61yWOV1Yc0ceF6zR17NomNyUg6gaaoqQHQqLi0Quc/9G6zCY0kGf/+/6ylm1XvNZo9BtFXXFqhiQvX+yU00tGks7i0wqTIAHOR1ACQdPRCWVlT1+JxhqSK6kNaU74vNoHBT73X0Kylm9VcSknSiXhHUgOgxQtlILsPHAp+ECJuTfm+Jj00xyLpRDwjqQEQ9ELZnKz0lChFg5aEmkySdCIeUSgMoFUXQJek7MyGmTaIvVCTSZJOxCN6agC0+gI4Y1Q+U4dNMiivo3IyUxSo9V1qmAVF0ol4RFIDIOiF0icnM4Xp3CZLTHBpxqh8SWryefluk3QiXpHUAGjxQukz5dLe+uiui0loLGB4QY7mjj1b2Zn+PWzZJJ2Ic9TUAJB09ELJgm72MLwgR0X52awoDByDpAZAIy6U9pKY4FJhr05mhwFYBkkNAD9cKAHYFTU1AADAEeipAQDYCruTIxCSGgCAbbA7OVrC8BMAwBbYnRzBkNQAACyP3ckRCpIaAIDlsTs5QkFSAwCwPHYnRyhIagAAlsfu5AgFSQ0AwPLYnRyhIKkBAFgeu5MjFCQ1ABCieq+hVduqtPSzCn1V7WKmTYyxOzmCYfE9AAhB00XfErX4kZWaefkZXExjiE1X0RKSGgAIwrfo2/H9Mrtq6jRx4Xp6CWKMTVcRCMNPANACFn0D7IOkBgBawKJvgH2Q1ABAC1j0DbAPkhoAaAGLvgH2QVIDAC1g0TfAPkhqAKAFLPoG2AdJDQAEEXjRNzfTuQELYZ0aAAjBsYu+Veyv1T83lWjy1RcqxZ1sdmgA/o2kBgBC5Fv0zePJ0LJ/bWDICbAYhp8AAIAjkNQAAABHYPgJQNTVew02IAQQdSQ1AKKq6e7WDeu6zBiVz6whABHF8BOAqPHtbn383kmV1Yc0ceF6FZdWmBQZACciqQEQFexuDSDWSGoARAW7WwOINWpqAAeyQmEuu1sDiLWQk5rXX3895Ce9/PLLwwoGQNtZpTCX3a0BxFrISc3o0aNDOs7lcqm+vj7ceACEqLnemOWbKzVx4fomdSy+wtxY7lPk2926svpQs3U1LknZ7G4NIIJCTmq8Xm804wDQCs31xmRnpOjQkfqAhbkuNRTmFuVnx2Qoyre79cSF6+WS/OKKxu7WVhhyA2AuWxUKr1y5UqNGjdJJJ50kl8ul1157zeyQgJgLOE265pD2H/QE/DkzCnMD726dEtFeo+LSCl3w8Hu6dt5q3fZCia6dt1oXPPweU8aBOBN2oXBtba1WrFihHTt26PDhw36P3XrrrW0OLNBr9uvXT9dff72uvPLKqLwGYGUtTZMOVawLc4/d3ToavSi+JM8KQ24AzBVWUrNhwwaNHDlSBw8eVG1trTp27Ki9e/cqLS1NWVlZUUtqRowYoREjRkTluQE7CDZNOhRmFOb6dreOtGBr4cR6yA2AucJKaqZMmaJRo0bpySefVGZmplavXq2kpCSNHTtWt912W6RjDFtdXZ3q6uoab9fU1EiSPB6PPJ7A3fSx5ovFSjHZRby1XcX+2rB/tqEw162zuqX7/Q7Yue0+DXEtnFVluzU4ggXJTmg7s9B24Yn3dgv1fbsMw2h1T3aHDh306aef6rTTTlOHDh20atUqnX766fr00081btw4ffHFF60OuLVcLpdeffXVFmdlzZw5U7NmzWpy/6JFi5SWlhbF6IDo+Krapcc3J4ZwpK+f4tjb0vWnetWvk3NW8F2316XnvgreHr/uXa8BnZ3zvoF4c/DgQY0ZM0bV1dXKyMgIeFxYPTVJSUlKSGioMc7KytKOHTt0+umnKzMzU9988014EUfBtGnTNHXq1MbbNTU1ys3N1dChQ1tslFjzeDxavny5ioqKlJSUZHY4thJvbVfvNbT4kZXaVVMXcJp0ZlqS3Iku7TpwtNYtJzNF94zoo2FndG28zwlt16l8n577am3Q44b+eHDEe2rs3nZmoe3CE+/t5htpCSaspOass87SP/7xD/Xu3Vs/+clPNH36dO3du1fPP/+8CgoKwnnKqHC73XK73U3uT0pKsuSXwqpx2UG8tF2SpJmXn9HiNOmHrjyzVYW5dm67wlOyQloLp/CUrKjU1Ni57cxG24UnXtst1Pcc1pTuBx98UDk5DbMJHnjgAZ144omaOHGi9uzZo7/85S/hPCWAEIUyTdpXmHtF/5NV2KuTY4tkfWvhSP6DbcfejuRaOACsLayemoEDBzb+OysrS8XFxRELqCXff/+9ysrKGm+Xl5erpKREHTt2VPfu3WMSA2AF0Z4mbSe+JK/JYoQmbA0BwFy22tBy7dq1uuiiixpv++plxo0bpwULFpgUFWCOaE2TtiOSPABSmElNXl6eXK7AJ4t//vOfYQfUkiFDhiiMyVoA4gBJHoCwkprf/va3frc9Ho82bNig4uJi3XHHHZGICwAAoFXCSmoCLbD3xBNPaO3a4NMrAbABIwBEWkRrakaMGKFp06Zp/vz5kXxawHGa22U7h8JWAGiTiO7SvXjxYnXsGLkFrgAnCrjL9r83YLTTztL1XkOrtlVpScm3WrWtSvVeat4AmCfsxfeOLRQ2DEOVlZXas2eP/vznP0csONgfQyz+nLQBI71NAKwmrKTm+P2WEhIS1KVLFw0ZMkR9+vSJRFxwAC56TQXbZdu3AeOa8n2Wnsnj6206Pjnz9Tb5FgEEgFgKK6mZMWNGpOOAw3DRa97uA4ETmnCOM4OTepsAOEvISU2om0lJstRmkYg9LnqBZaWnBD+oFceZwSm9TbAuhq0RrpCTmg4dOrS44N6x6uvrww4I9sdFL7BBeR1D2oBxUAR3lI40J/Q2wboYtkZbhJzUvP/++43/3r59u+6++26NHz9ehYWFkqRVq1bp2Wef1ezZsyMfJWyFi15gvg0YW9pl2+obMDqhtwlHWalXhGFrtFXISc1PfvKTxn//7ne/05/+9Cdde+21jfddfvnlOvPMM/WXv/xF48aNi2yUsBUuei2z+waMTuhtQgMr9YowbI1ICKtQeNWqVXryySeb3D9w4EBNmDChzUHB3rjoBWfnDRid0NsE6/WKMGyNSAhr8b3c3FzNmzevyf1PP/20cnNz2xwU7M130ZOOXuR8uOgd5duA8Yr+J6uwVydbtYevtyk707+3LTszhSECGwjWKyI19IrEcjFFhq0RCWH11MyZM0c///nP9dZbb2nw4MGSpDVr1uirr77Syy+/HNEAYU92H2JBcHbubYp3VuwVYdgakRBWUjNy5Eht3bpVc+fO1RdffCFJGjVqlG6++WZ6atCIi57z+XqbYC9W7BVh2BqREPaGlrm5uXrwwQcjGQsciIseYD1W7BWhVguREHJS89lnn6mgoEAJCQn67LPPWjy2b9++bQ4MAEJhpSnJdmHVXhGGrdFWISc1/fv3V2VlpbKystS/f3+5XC4ZRtNfB5fLxeJ7AGLCSlOS7cTKvSIMW6MtQk5qysvL1aVLl8Z/A4CZrDYl2W6s3CvCsDXCFXJS06NHj2b/DQCxxkJtkUGvCJwmrHVqnn32Wb355puNt++880516NBB5513nr7++uuIBQcAzWnNlGS0zM7rJQHHCyupefDBB5WamiqpYXXhxx9/XL///e/VuXNnTZkyJaIBAsDxrDglGYD5wprS/c033+iUU06RJL322mv6xS9+oZtuuknnn3++hgwZEsn4AKAJK05JBmC+sHpq2rdvr6qqKknSO++8o6KiIklSSkqKfvjhh8hFBwDN8E1JDjRQ4lLDLCgWagPiS1hJTVFRkSZMmKAJEyZo69atGjlypCRp06ZN6tmzZyTjA4Am2F8sfPVeQ6u2VWlJybdata0qpvs7AdEW1vDTE088oXvvvVfffPONXn75ZXXq1DD1bt26dbr22msjGiAANMfKU5KtKpx1fVjc0H7i+TMLK6np0KGDHn/88Sb3z5o1q80BAUComJIcunDW9WFxQ/uJ988srOEnSfrwww81duxYnXfeefr2228lSc8//7w++uijiAUHAMEwJTm4YOv6SA3r+hw7FOVLgo6fOu9LgopLK6IXMMLCZxZmUvPyyy9r2LBhSk1N1fr161VXVydJqq6uZpNLALCY1q7rE04SBHPxmTUIK6m5//779eSTT2revHlKSkpqvP/888/X+vXrIxYcACA8xxYEf1y2N6Sf8a3rw+KG9sNn1iCsmpovv/xSF154YZP7MzMztX///rbGBABog+bqKkLhW9eHxQ3tJ9TP4q1/D0E5tfYsrJ6a7OxslZWVNbn/o48+0o9+9KM2BwUACE+guoqWHL+uD4sb2k+on8Vzq77WtfNW64KH33NkjU1YSc2NN96o2267TZ9++qlcLpd27typv/3tb7r99ts1ceLESMcIAAhBS3UVgTS3rg+LG9pPsM/seE4tHg4rqbn77rs1ZswYXXLJJfr+++914YUXasKECZo4caImTJgQ6RiBuMHCaGiLYHUVzcnOTGkynZvFDe2npc+sOU4tHg6rpsblcumee+7RHXfcobKyMn3//ffKz8/XU089pby8PFVWVkY6zrhQ7zW0dlsV623EqViuL1HvNfRp+T6t2+tSp/J9Kjwli++aA4RaVzH5ol7q3TW9xfMMixvaT6DPLJBji4cLe3WKfoAx0Kqkpq6uTjNnztTy5cvldrt1xx13aPTo0Zo/f75+9rOfKTExkV26w7SxyqXZj6xUZU1d433xtGBSvAtnYbS2vNbRk16invtqLd81hwi1ruL8U7qEdBFjcUP7OfYze6u0Qs+t+jrozzip4LtVw0/Tp0/X3Llz1bNnT5WXl+uqq67STTfdpDlz5uiRRx5ReXm57rrrrmjF6lhvb9qlZ7Ym+CU0knPHPJ0iUkNFsVxfgsW5nC0atTAsbmg/vs9sRIh/pDip4LtVPTUvvfSSnnvuOV1++eUqLS1V3759deTIEW3cuFEuF1/0cNR7Dd2/7ItmHzPUcBKatXSzivKzOZlEWFv2R4nkUFFr1pdoSxdxsOSJ75r9+eoqJi5cL5fk91lTCxN/fEluZfWhZn/vXWoYTnRSwXeremr+9a9/acCAAZKkgoICud1uTZkyhYSmDdaU7/t3D03zbRgvCybFWnFphS54+D1dO2+1bnuhpFVTHCPd2xGrNUFYnCs++OoqsjP9//puriAYzhaPBd+t6qmpr69XcnLy0R8+4QS1b98+4kHFExa5ir221K9Eo7cjVmuC8F2LH9TCwCfeCr5bldQYhqHx48fL7XZLkg4dOqSbb75Z7dq18zvulVdeiVyEDsciV7HV1qQkGkNFseoi5rsWX3x1FUA8JbmtSmrGjRvnd3vs2LERDSYeDcrrqOwMtyprDqm5ISgnjnmaqa1JSTR6O2JVBxGP4+sAGsRLktuqpGb+/PnRiiNuJSa4dO/IPpr8QgmFfTHQ1qQkWr0dsegipogUgNOFtfgeImvYGV11/aleLatM85vW7dQxTzO1NSmJZm9HLLqI4218HUB8IamxiH6dDN35ywu14V8HHD/maaa2JiXR7u2IRRexL3laVbZb73z4qYb+eDArCgNwhLD2fkJ0sMhV9EViiqMTpswmJrg0OK+jBnQ2NJjk2ZLYBwxoPXpqEHciMQQTT7MJEHux3AcMcBKSGsSlSCQl8TKbALEVy33A4lVbVhOHtZHUIG6RlMBq2Moi+ugFczZqagCLo7YifrCVRXQF2uKkovqQbmZDV0egpwawMP6qjC9sZRE9LfWC+dz9yuf0gtkcPTWARUV640xYH1tZRE+wXjBJ2n/Qo8ffK4tRRIgGkhrAgoLVVkgNtRUMRTmLbx2lQP0ELjX01LGVReuF2rs1/5Nyfq9sjKQGsCBqK+JTJNZRQvNC7d3af9DD75WNkdQAFkRtRfxywuKOVjQor6M6pCaFdCy/V/ZFoTBgQdRWxDcWd4y8xASXrju/p+b8/augx/J7ZV8kNYAFRXPjTNgD6yhF3uSLe2v+J9u1/6Cn2cf5vbI/hp8AC6K2AlbhpHWSEhNceujKM5stxOb3yhnoqUGLWE7cPJHYo8qO+M5ZhxPXSWrp9+r//TRfmanJWlLyLd89myKpQUBOPKHZTbzVVkTjO0eSFB4n70HV3O/Vd7V1uu9Nznd2R1KDZjn5hGY38VJbEY3vHIl5eOJhD6pjf6+KSys0adEGzncOQE0NmmDhN8RaNL5zrMgcvnhaJ4nznbOQ1KCJeDqhwRoi/Z3jQtU28bROEuc7ZyGpQRPxdEKDNUT6O8eFqm3iaZ0kznfOYruk5oknnlDPnj2VkpKiwYMHa82aNWaH5DjxdEKDNYT6Xerczh3ScVyo2iae9qDifOcstkpqXnzxRU2dOlUzZszQ+vXr1a9fPw0bNky7d+82OzRHiacTGqwh2HfO5/aXNoZUC8OFqm3iaZ0kznfOYqvZT3/6059044036rrrrpMkPfnkk3rzzTf1zDPP6O67725yfF1dnerq6hpv19TUSJI8Ho88nuZXlDSDLxYrxXTPiNN0ywsb5ZL86hJcxzzurT8ib70JwR3Dim1nF1Zru0DfuWPtqmko8n3smn4adkbXgM91Vrd0ZWe4taumroUVmd06q1t6WO/fam0XDZec1lmPXdNP9y/7QpU1R8+j2Zlu3TOijy45rbNj2s4O5zsrtlsshfq+XYZh2KJS7vDhw0pLS9PixYs1evToxvvHjRun/fv3a8mSJU1+ZubMmZo1a1aT+xctWqS0tLRohusIG6tcemV7gvYfPvo3TIdkQ1f29KpfJ1t8bWACryFtq3GpxiNlJEm9MgyF+gf9xiqXXi5PULWnpR8w1CFZmnF2fYvPu7HKpWe2+jqjjz2w4bt7/al8j0PRls/TTjjfWdvBgwc1ZswYVVdXKyMjI+Bxtklqdu7cqZNPPlmffPKJCgsLG++/8847tWLFCn366adNfqa5nprc3Fzt3bu3xUaJNY/Ho+XLl6uoqEhJSaHtIhsr9V5Da7/+TrsP1Ckr3a2BPU60VJezldvO6qLRdm9v2tX0L/sMt+4d2afFnpVjfbKtSuMWrAt63MLrB2pwkCGB5uLJ+XdPQ6jxNIfvXfis3HZWPt9Zud1ioaamRp07dw6a1Nhq+Km13G633O6mhYVJSUmW/FJYMa4kSRecGv7JP1as2HZ2Eam2Ky6t0C0vbGwy3LOrpk63vLAx5AXM9h8KrY+/6uCRoHFf1r+bRvQ9OWorCvO9C58V284O5zsrtpsU/ZW7Q33PtklqOnfurMTERO3atcvv/l27dik7O9ukqABIkV2BNtJFvvGyIrPd1HsNrd1WxfYVDmCllbttM/spOTlZAwYM0Lvvvtt4n9fr1bvvvus3HAUg9iK5LgyzUZxvY5VLQx5ZqWvnrdZtL5To2nmrdc4Df9eyz1jl2W6stnK3bZIaSZo6darmzZunZ599Vlu2bNHEiRNVW1vbOBsKgDkiuS5MPE0njkdvb9qlZ7Ym+NU5SdK+2sP6z0Xr9cCbm02KDK1lxZW7bTP8JElXX3219uzZo+nTp6uyslL9+/dXcXGxuna19hgoYBXHjnt3SjtBkTrXRHrIaHhBjuaOPbtJl3a2AzejjKddxOu9hu5f9kWLx8z7sFySoXt+ekZsgkLYWtNDG6shYFslNZI0efJkTZ482ewwANtpbty7Q3Kiknru0mX9u7XpuX1DRpXVh1pYF6Z1Q0bDC3JUlJ/t6Au+lWoRYmFN+b5/99C0/BnO+3C7zso9USP7nhSbwBAWK67cbavhJwDhCTTuvf+wdMsLoa3S25JoDRn5inyv6H+yCnt1clxCY6VahFhozcXt3iWlbDhqcVZcuZukBnC4lsa9fSlHJMa9fUNG2Zn+J7DszJSQp3PHCyvWIsRCay5u+2o9bDhqcVYs6rfd8BMQKfFSyxDLce94GDKKBCvWIsTCoLyO6piWpH0HQ1vyng1Hrc3XQztx4fqAW0zEuqifpAZxKZ5qGWI97s26MMHF8jOxUvKemODSzFGn69YXNypYXY3EhqN2YLWifpIaxB1fLcPxHfu+WganDZVYcdw73sXqM7Fi8j6iIFsXrdyg9ysSWzyOtYjsw0o9tNTUIK7EYy2DFce9410sPhMrFyKP7mno+vO6B3zcJdYishurFPWT1CCuRHLlW7toaWaSL5Wz+wWk3mto1bYqLSn5Vqu2VVk+KY32AoN2SN6njeijP485Sx3b+e/pk0NhOdqA4SfEFSuuqxALwwty9MSYs3TvklLtqz1apNkhWbr/yn62voBYcYglFNGsRbBLIfLIvidpWEGOJYYt4AwkNYgrTq0vCVYMWlxaofve3OKX0HRMS9IV3Q5p2Bn2XZHb7vVR0apFsFPyTmE5IomkBnElGivfmi1YT0WgC/93Bz2avzVBAza1fUVhM0RyZ3AzReOi7tTkHQiGmhrEFadtlhisGHTZZxVBayseeOsLy9egNCce66NC1dZCZLvVKAE+9NQg7lhtXYVwhdJT8f+WlKqq9nALz+JSRXWd6bUV4bDTEEustWVRNLvWKAESSQ3ilJXWVQhXKD0VLSc0R9nxws8QS8vCSd5jVaPkNaRPy/ep6uARW/7uwbpIahC37F6gGMlExI4XfifWR0Vaa5L3WNUovb1pl2atT9T+1Wsb76MnCJFCTQ1gU6EmIh3bJbWwIL2hnEy3LS/8TquPipZQF0WLRY1ScWmFbnlho/Yf14FohQUB4QwkNUCMRaoIM9Ri0PuvKGi8ffzjknTPiD62vfCzM3jkRLtGyb8nyP/7ZpUFAWF/DD8BMRTJIsxQi0GHF+RoboKrmdoKt0Z0PWjrdWokZ9RHWUG0a5TssiAg7I2kBoiRaBRhhloM2tyF/6xu6Xq7+K22vi1LsHt9lBVEu0aJ2WqIBZIaIAaiWYQZak/F8Rd+j8dz/FMhjrVlGngomK2GWKCmBoiBaBdhWmWHXNhbNGuU2C0esUBPDRADdL3Hr2D7cllNtGqUju0JOto/2YDZaogUkhoEZLeTcSy1tm3oeo9Pdl2dN1o1SsMLcvTYNf107yslftO67baaN6yLpAbNsuvJOBbCaRsWios/dt9BPFqGndFVnu316pJ/LisKI+KoqUETwTZJjOcFssJtGxaKiy/BCsOl+F6TJcElDc7rSA0YIo6kBn44GQfW1rZhobjQOGGHaHYQB8zB8BP8sEBWYJFoGxaKa5kZw57RqB2jMBwwB0kN/HAyDixSbcNCcc0zowaluSSqY7tk3X9FgUb2Df+1rFoYTvE/nI6kBn6sejK2AtomemK1Q/SxAiVR+2oP6z8Xrddv/pWnaSPzw3puKxaGU/yPeEBNDfywQFZgtE30xLoGpaUkyuepleV6o+TbsJ7faoXhFP8jXpDUwI/VTsZWQttET6yHPYMlUT63vFiiZZ+Fd8G3SmE4xf+IJww/oYlQN0mMR7RNdMR6aC/U5MgwpP9ctF5PJoSXhFihMJzif8QTkho0ywonY6uibSIv1jUorU2O/vvVz3Vxn65KPqH1ndtmF4ZT/I94QlKDgMw+GVsZbRNZ0d4h+niD8jqqY7tk7as9HPxgSftqPTp39rt68GcFtuuN2773YEjHUeAOJ6CmBoAlxLIGJTHBpfuvKGjVz+yrPWy7otri0go9+vetQY/r2C5JlTWHbLvYIeBDTw0Ay2jt0F5b1l0Z2TdHv/lXnp5aWd6qGGct3awhvX/cqp8xQygzvHz21Xo05cUSSUzzhr2R1CBqWOjLmaL9uYY6tBeJhfOmjczXmSdl6pYXS2SEcPX3FdU+t/prdbF4h0aoM7yOF+8bbsLeSGoQFSz05UxW+VwjuXDeZf1PVkJCgv5z0fqQX//Bt7aqQ3Kiknru0mX9u7Ui8tgJt/A3WosdArFATQ0ijoW+nMkqn2uoC+ct+2xnyM85sm+Onhx7tjq2Swr5Z/Yflm55YaNlv89tKfxlw03YFUkNIirYQl+GGqbHvrrBvjswxyMrLeAW6rDKvUtKWxXP8IIcrZ52qTq2Sw7xJxp6MKy6cF2wFbBDwTTvyHPCLvRWxvATIiqUCw5FifZjpQXcQr3Q7qv1tDqe5BMS9ODPCjRxYcNQVLDLTSzed7g1TC1Nkw8V07wjyyrDt05GTw0iqrV/2TEkZQ9WWsCtNRfacOIJNLU80q8TiuLSCl3w8Hu6dt5q3fZCia6dt1oXPPxeyL8vgd5LTmaKOqQlsY9ZDFll+Nbp6KlBRLX2LzuKEu2hNdsYRHt2VGsWzgu3p8E3tXzBx+W6780tUXudlgQqhm7t7KRA0+SXb66M2WKH8c6MXejjFT01iKhwxvEpSrS+UHco/662rk09C6EIdeG8tvY0JCa4NP78PFN2Zo90DZNvmvwV/U9WYa9OSkxwWWbDzXgQ613o4xlJDSKqpZ2sg6Eo0bpC2aH88n45mrRoQ0y610f2zdFvLswL+LhL/j0N4RZntvx9bniOaPRoxOoiOLwgRx/ddbH+78Zz9T/X9Nf/3XiuPrrrYhKaCLPS8K3TkdQg4sKpSZAoSrS6lv6yf2LM2Xp9Y0VMZ0dNG5mvP485q8k07JzjehqiVZfSIVl67Jp+UUkAYnkRbK4XB5EV613o4xk1NYiKY8fxK6t/0H1vbtF3tYdjsgMzoidQfYZZs6NG9j1JwwpyAtbwRKsupVPaCdqzebWGndE1Yu/lWFwEnSXWu9DHM5IaRM2xy92nJidSlOgQzW1jYGb3eqBtFSJdnHns63g8Hi0LXj8cNjtfBNkepalY70Ifzxh+QkxQlOhsVuxZsHNxZig1TFa8CLZ1qM/JOAfGBj01iJnW7sAM+7Biz4LdizN9F8HjF2vLtuhibZEa6nMyzoHRR1KDmAp1B2bYixW7163Ye9RadrkItmaoL95xDowuhp8ARES0u9dbOy071LV1rFiXciw7zE6y81AfnIWeGgARE62ehXD2zLFi75FTtW6oLyO6wSCu0VMDIKIi3bPQlj1zKM6MDScM9cEZ6KkBYFmRmJZtl7oUO2tNobi3/kisw0McoacGgGVFqlbDDnUpdmbXKehwHpIaAJZl92nZ8YShPlgBw08ALItaDXthqA9mI6kBYFlWXNQPLWMdFpiJ4ScAlkWtBoDWIKkBYGnUagAIFcNPACyPWg0AoSCpAWAL8VqrUe81SOaAEJHUtFG0TjhOOpE56b1YldltbPbrO1Vz20N0SE3SdefnafLFp9DGwHFsk9Q88MADevPNN1VSUqLk5GTt37/f7JDC2o/GzOc1g5Pei1WZ3cZmv75T+baHOH7W1/4fPJrz962a/0m5HrryTNoYOIZtCoUPHz6sq666ShMnTjQ7FElt24/GjOc1g5Pei1WZ3cZmv75TtbQ9hM/+gx7dTBsDfmyT1MyaNUtTpkzRmWeeaXYoQfejkRr2o6n3tnRKit3zmsFJ78WqzG7jaLx+vdfQqm1VWlLyrVZtq4rb70ew7SGOxe8RcJRthp/CUVdXp7q6usbbNTU1kiSPxyOPxxP2834a4n40q8p2a3AIi4L5Ylm9bU9En9dMkW6jQHxt15bP067a2sZtbbtIf8Zvb9ql+5d9ocqao7+z2Rlu3Tuyj4ad0TWsGKMl2t+7iv21oR9rk3OCT6Tart5raO3X32n3gTplpbs1sMeJjq4xiudznRT6+3Z0UjN79mzNmjWryf3vvPOO0tLSwn7edXtdkhKDHvfOh5+qakvof0G9t2pdVJ7XDNFqo0CWL1/e5uewm0i1cbhtF8nPeGOVS89s9XUcH70wVdYc0uQXSnT9qV7162S973y0vnf/rA6tbX3scE44XlvabmOVS69sT9D+w0e/Kx2SDV3Z05rfk0iKx3OdJB08eDCk40xNau6++249/PDDLR6zZcsW9enTJ6znnzZtmqZOndp4u6amRrm5uRo6dKgyMjLCek5J6lS+T899tTbocUN/PDjknprly5fr4sIBeu6rkog9r5ki3UaB+NquqKhISUlJYT+PHbW1jdvadpH6jOu9hmY/slJSXTOPuuSS9NauNN35ywst85d4tL939V5Dix9Z6ddr1RI7nBN82tp2b2/apfmrNjYZ9qw+7NL8rYl67Jp+luvZi4R4PtdJR0dagjE1qbn99ts1fvz4Fo/50Y9+FPbzu91uud3uJvcnJSW16UtReEpWSPvRFJ6S1aqT8Lm9ukTlec0QrTYKpK2fqR1Fqo3DbbtIvf7abVUtXrwbhrHqtOFfByy3Tk20vndJkmZefkazs5+OZadzwvHCabt6r6EH3voyYB2XS9IDb32pEX1Ptl17hCoez3WSQn7PphYKd+nSRX369Gnxv+TkZDNDbFa09qNx0j43TnovVmV2G0fq9XcfCK0gNtTjnMK3PUSHtOZP5vH4exSsgNpXx7WmfF/sgoKl2Gb2044dO1RSUqIdO3aovr5eJSUlKikp0ffff29KPNHaj8ZJ+9w46b1YldltHInXz0pPCXpMa45zkuEFOVp3b5GmXNpbHVL9k5t4/D0iAUYwtikUnj59up599tnG22eddZYk6f3339eQIUNMiSla+9E4aZ8bJ70XqzK7jdv6+oPyOoY0jDXIJjUjkZaY4NJtl56qyRf3jvvfIxJgBGObpGbBggVasGCB2WE0Ea39aJy0z42T3otVmd3GbXl93zDWxIXr5ZL8Ept4HGIJxOzP2ApIgBGMbYafADiX2cNosAez68hgfbbpqQHgbGYPo0UDG31Gni8BPn6/sWz2G4NIagBYiJOGWNjoM3qcmAAjMkhqACDCAu2w7dvokyG1tnNSAozIoaYGACLI7I1GgXhGUgMAEcQCcYB5SGoAIIJYIA4wD0kNAEQQC8QB5qFQGAAiiAXiwFR+85DUAEAEsUJy29k5KWAqv7lIagAgwlggLnx2TgqYym8+khoAiAIWiGs9OycFwabyu9Qwlb8oP5vvQBSR1ABAlLBAXOjsnhS0Zio/34noYfYTAMB0dl/fh6n81kBSAwAwnd2TAqbyWwNJDQDAdHZPCnxT+QMNjLnUUPDMVP7oIqkBAJjO7kmBbyq/pCbvgan8sUNSAwAwnROSAt9U/uxM/96k7MwUS8/cchJmPwEALMEJ6/swld9cJDUAAMtwQlLAVH7zkNQAACyFpADhoqYGAAA4AkkNAABwBJIaAADgCCQ1AADAESgUBgC0Sr3XsPXsJDgXSQ0AIGTFpRVN1pHJsdE6MnA2hp8AACEpLq3QxIXrm+ymXVl9SBMXrldxaYVJkQENSGoAAEHVew3NWrpZRjOP+e6btXSz6r3NHQHEBkkNACCoNeX7mvTQHMuQVFF9SGvK98UuKOA4JDUAgKB2Hwic0IRzHBANJDUAgKCy0lOCH9SK44BoIKkBAAQ1KK+jcjJTFGjitksNs6AG5XWMZViAH5IaAEBQiQkuzRiVL0lNEhvf7Rmj8lmvBqYiqQEAhGR4QY7mjj1b2Zn+Q0zZmSmaO/Zs1qmB6Vh8DwAQsuEFOSrKz2ZFYVgSSQ0AoFUSE1wq7NXJ7DCAJhh+AgAAjkBSAwAAHIGkBgAAOAJJDQAAcASSGgAA4AgkNQAAwBFIagAAgCOQ1AAAAEcgqQEAAI4QVysKG4YhSaqpqTE5En8ej0cHDx5UTU2NkpKSzA7HVmi78NF24aPtwkfbhSfe28133fZdxwOJq6TmwIEDkqTc3FyTIwEAAK114MABZWZmBnzcZQRLexzE6/Vq586dSk9Pl8tlnc3XampqlJubq2+++UYZGRlmh2MrtF34aLvw0Xbho+3CE+/tZhiGDhw4oJNOOkkJCYErZ+KqpyYhIUHdunUzO4yAMjIy4vLLGgm0Xfhou/DRduGj7cITz+3WUg+ND4XCAADAEUhqAACAI5DUWIDb7daMGTPkdrvNDsV2aLvw0Xbho+3CR9uFh3YLTVwVCgMAAOeipwYAADgCSQ0AAHAEkhoAAOAIJDUAAMARSGosZuvWrbriiivUuXNnZWRk6IILLtD7779vdli28eabb2rw4MFKTU3ViSeeqNGjR5sdkq3U1dWpf//+crlcKikpMTscy9u+fbtuuOEG5eXlKTU1Vb169dKMGTN0+PBhs0OzpCeeeEI9e/ZUSkqKBg8erDVr1pgdkuXNnj1b55xzjtLT05WVlaXRo0fryy+/NDssyyKpsZjLLrtMR44c0Xvvvad169apX79+uuyyy1RZWWl2aJb38ssv61e/+pWuu+46bdy4UR9//LHGjBljdli2cuedd+qkk04yOwzb+OKLL+T1evXUU09p06ZNmjNnjp588kn993//t9mhWc6LL76oqVOnasaMGVq/fr369eunYcOGaffu3WaHZmkrVqzQpEmTtHr1ai1fvlwej0dDhw5VbW2t2aFZkwHL2LNnjyHJWLlyZeN9NTU1hiRj+fLlJkZmfR6Pxzj55JONp59+2uxQbGvZsmVGnz59jE2bNhmSjA0bNpgdki39/ve/N/Ly8swOw3IGDRpkTJo0qfF2fX29cdJJJxmzZ882MSr72b17tyHJWLFihdmhWBI9NRbSqVMnnXbaaXruuedUW1urI0eO6KmnnlJWVpYGDBhgdniWtn79en377bdKSEjQWWedpZycHI0YMUKlpaVmh2YLu3bt0o033qjnn39eaWlpZodja9XV1erYsaPZYVjK4cOHtW7dOl166aWN9yUkJOjSSy/VqlWrTIzMfqqrqyWJ71gAJDUW4nK59Pe//10bNmxQenq6UlJS9Kc//UnFxcU68cQTzQ7P0v75z39KkmbOnKl7771Xb7zxhk488UQNGTJE+/btMzk6azMMQ+PHj9fNN9+sgQMHmh2OrZWVlemxxx7Tb37zG7NDsZS9e/eqvr5eXbt29bu/a9euDK23gtfr1W9/+1udf/75KigoMDscSyKpiYG7775bLperxf+++OILGYahSZMmKSsrSx9++KHWrFmj0aNHa9SoUaqoqDD7bZgi1Lbzer2SpHvuuUc///nPNWDAAM2fP18ul0svvfSSye/CHKG23WOPPaYDBw5o2rRpZodsGaG23bG+/fZbDR8+XFdddZVuvPFGkyKHk02aNEmlpaV64YUXzA7FstgmIQb27NmjqqqqFo/50Y9+pA8//FBDhw7Vd99957e1fO/evXXDDTfo7rvvjnaolhNq23388ce6+OKL9eGHH+qCCy5ofGzw4MG69NJL9cADD0Q7VMsJte3+4z/+Q0uXLpXL5Wq8v76+XomJifrlL3+pZ599NtqhWk6obZecnCxJ2rlzp4YMGaJzzz1XCxYsUEICfy8e6/Dhw0pLS9PixYv9ZiSOGzdO+/fv15IlS8wLziYmT56sJUuWaOXKlcrLyzM7HMs6wewA4kGXLl3UpUuXoMcdPHhQkpqcEBMSEhp7IuJNqG03YMAAud1uffnll41Jjcfj0fbt29WjR49oh2lJobbd//7v/+r+++9vvL1z504NGzZML774ogYPHhzNEC0r1LaTGnpoLrroosbeQRKappKTkzVgwAC9++67jUmN1+vVu+++q8mTJ5sbnMUZhqFbbrlFr776qj744AMSmiBIaiyksLBQJ554osaNG6fp06crNTVV8+bNU3l5uX7605+aHZ6lZWRk6Oabb9aMGTOUm5urHj166A9/+IMk6aqrrjI5Omvr3r273+327dtLknr16qVu3bqZEZJtfPvttxoyZIh69OihP/7xj9qzZ0/jY9nZ2SZGZj1Tp07VuHHjNHDgQA0aNEiPPvqoamtrdd1115kdmqVNmjRJixYt0pIlS5Sent5Yg5SZmanU1FSTo7MekhoL6dy5s4qLi3XPPffo4osvlsfj0RlnnKElS5aoX79+ZodneX/4wx90wgkn6Fe/+pV++OEHDR48WO+99x5F1oia5cuXq6ysTGVlZU0SQEb2/V199dXas2ePpk+frsrKSvXv31/FxcVNiofhb+7cuZKkIUOG+N0/f/58jR8/PvYBWRw1NQAAwBEY/AUAAI5AUgMAAByBpAYAADgCSQ0AAHAEkhoAAOAIJDUAAMARSGoAAIAjkNQAAABHIKkB4Mflcum1114zOwxL+OCDD+RyubR//35J0oIFC9ShQwdTYwIQGEkNEGfGjx/vt1Py8SoqKjRixIjYBWQjV199tbZu3Rry8T179tSjjz4avYAA+GHvJwB+rLARo2EYqq+v1wkntP0UFcnnSk1NNWUTwcOHDys5OTnmrwvYDT01APwcO/y0fft2uVwuvfLKK7rooouUlpamfv36adWqVX4/89FHH+nHP/6xUlNTlZubq1tvvVW1tbWNjz///PMaOHCg0tPTlZ2drTFjxmj37t2Nj/uGed566y0NGDBAbrdbH330UZPYfPG88MILOu+885SSkqKCggKtWLEi6HN5vV7Nnj1beXl5Sk1NVb9+/bR48WK/51+2bJlOPfVUpaam6qKLLtL27dv9Hm9u+Gnp0qU655xzlJKSos6dO+tnP/uZpIYNCL/++mtNmTJFLpdLLper8WdefvllnXHGGXK73erZs6ceeeQRv+fs2bOn7rvvPv36179WRkaGbrrppgCfFgA/BoC4Mm7cOOOKK64I+Lgk49VXXzUMwzDKy8sNSUafPn2MN954w/jyyy+NX/ziF0aPHj0Mj8djGIZhlJWVGe3atTPmzJljbN261fj444+Ns846yxg/fnzjc/71r381li1bZmzbts1YtWqVUVhYaIwYMaLx8ffff9+QZPTt29d45513jLKyMqOqqqpJbL54unXrZixevNjYvHmzMWHCBCM9Pd3Yu3dvi891//33G3369DGKi4uNbdu2GfPnzzfcbrfxwQcfGIZhGDt27DDcbrcxdepU44svvjAWLlxodO3a1ZBkfPfdd4ZhGMb8+fONzMzMxnjeeOMNIzEx0Zg+fbqxefNmo6SkxHjwwQcNwzCMqqoqo1u3bsbvfvc7o6KiwqioqDAMwzDWrl1rJCQkGL/73e+ML7/80pg/f76RmppqzJ8/v/F5e/ToYWRkZBh//OMfjbKyMqOsrCy0DxeIcyQ1QJwJJ6l5+umnGx/ftGmTIcnYsmWLYRiGccMNNxg33XST33N8+OGHRkJCgvHDDz80+xr/+Mc/DEnGgQMHDMM4moi89tprLcbui+ehhx5qvM/j8RjdunUzHn744YDPdejQISMtLc345JNP/J7vhhtuMK699lrDMAxj2rRpRn5+vt/jd911V4tJTWFhofHLX/4yYLw9evQw5syZ43ffmDFjjKKiIr/77rjjDr/X7tGjhzF69OiAzwugeQw/AQiqb9++jf/OycmRpMbho40bN2rBggVq375943/Dhg2T1+tVeXm5JGndunUaNWqUunfvrvT0dP3kJz+RJO3YscPvdQYOHBhSPIWFhY3/PuGEEzRw4EBt2bIl4HOVlZXp4MGDKioq8ovzueee07Zt2yRJW7Zs0eDBgwO+TnNKSkp0ySWXhBSzz5YtW3T++ef73Xf++efrq6++Un19fbPxAwgNhcIAgkpKSmr8t682xOv1SpK+//57/eY3v9Gtt97a5Oe6d++u2tpaDRs2TMOGDdPf/vY3denSRTt27NCwYcN0+PBhv+PbtWsXsZiPfa7vv/9ekvTmm2/q5JNP9jvO7XaH/RrRLBqOZFsA8YKkBkCbnH322dq8ebNOOeWUZh///PPPVVVVpYceeki5ubmSpLVr17bpNVevXq0LL7xQknTkyBGtW7dOkydPDnh8fn6+3G63duzY0dhLdLzTTz9dr7/+epPXaUnfvn317rvv6rrrrmv28eTkZL/eF9/rfPzxx373ffzxxzr11FOVmJjY4usBaBlJDRCHqqurVVJS4ndfp06dGpOO1rjrrrt07rnnavLkyZowYYLatWunzZs3a/ny5Xr88cfVvXt3JScn67HHHtPNN9+s0tJS3XfffW2K/4knnlDv3r11+umna86cOfruu+90/fXXBzw+PT1d//Vf/6UpU6bI6/XqggsuUHV1tT7++GNlZGRo3Lhxuvnmm/XII4/ojjvu0IQJE7Ru3TotWLCgxThmzJihSy65RL169dI111yjI0eOaNmyZbrrrrskNcxiWrlypa655hq53W517txZt99+u8455xzdd999uvrqq7Vq1So9/vjj+vOf/9ymNgEgZj8B8WbcuHGGpCb/3XDDDYZhNF8ovGHDhsaf/+677wxJxvvvv99435o1a4yioiKjffv2Rrt27Yy+ffsaDzzwQOPjixYtMnr27Gm43W6jsLDQeP311/2e11fc6yvIDcQXz6JFi4xBgwYZycnJRn5+vvHee+81HhPoubxer/Hoo48ap512mpGUlGR06dLFGDZsmLFixYrGY5YuXWqccsophtvtNn784x8bzzzzTIuFwoZhGC+//LLRv39/Izk52ejcubNx5ZVXNj62atUqo2/fvobb7TaOPd0uXrzYyM/PN5KSkozu3bsbf/jDH/yes7kCYwDBuQzDMMxIpgCgtbZv3668vDxt2LBB/fv3NzscABbD7CcAAOAIJDUAAMARGH4CAACOQE8NAABwBJIaAADgCCQ1AADAEUhqAACAI5DUAAAARyCpAQAAjkBSAwAAHIGkBgAAOML/B7EYZsk/P4vfAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.clf()\n", "plt.grid(True)\n", "plt.plot(result2.predict(linear=True), result2.resid_pearson, \"o\")\n", "plt.xlabel(\"Linear predictor\")\n", "plt.ylabel(\"Residual\")" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "main_language": "python", "notebook_metadata_filter": "-all" }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.8" } }, "nbformat": 4, "nbformat_minor": 4 }