{ "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": "2025-12-05T18:07:53.106017Z", "iopub.status.busy": "2025-12-05T18:07:53.105463Z", "iopub.status.idle": "2025-12-05T18:07:54.918198Z", "shell.execute_reply": "2025-12-05T18:07:54.917432Z" } }, "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": "2025-12-05T18:07:54.922532Z", "iopub.status.busy": "2025-12-05T18:07:54.922126Z", "iopub.status.idle": "2025-12-05T18:07:54.927589Z", "shell.execute_reply": "2025-12-05T18:07:54.926804Z" } }, "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": "2025-12-05T18:07:54.931265Z", "iopub.status.busy": "2025-12-05T18:07:54.930975Z", "iopub.status.idle": "2025-12-05T18:07:54.946603Z", "shell.execute_reply": "2025-12-05T18:07:54.944837Z" } }, "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": "2025-12-05T18:07:54.950467Z", "iopub.status.busy": "2025-12-05T18:07:54.949327Z", "iopub.status.idle": "2025-12-05T18:07:54.984874Z", "shell.execute_reply": "2025-12-05T18:07:54.983383Z" } }, "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: Fri, 05 Dec 2025 Deviance: 6.1260\n", "Time: 18:07:54 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": "2025-12-05T18:07:54.987576Z", "iopub.status.busy": "2025-12-05T18:07:54.987261Z", "iopub.status.idle": "2025-12-05T18:07:55.186696Z", "shell.execute_reply": "2025-12-05T18:07:55.185572Z" }, "lines_to_next_cell": 1 }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/hostedtoolcache/Python/3.10.19/x64/lib/python3.10/site-packages/statsmodels/genmod/generalized_linear_model.py:985: FutureWarning: linear keyword is deprecated, use which=\"linear\"\n", " warnings.warn(msg, FutureWarning)\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'Residual')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVQZJREFUeJzt3Xl8U1XeP/BPUtqUYhcKtClQaMUFSqUFOu3UFYUuoB1Rxx8IKDAsA2NdqI9KR6WiYl1AGZEBlwEX5IFn0BGRWumAiGKxss6wCxZwoAtQaGorbWju749OQkOW3iQ3yb25n/frpS9yc+69J6dp880533OORhAEAUREREQqpPV3BYiIiIj8hYEQERERqRYDISIiIlItBkJERESkWgyEiIiISLUYCBEREZFqMRAiIiIi1erk7wrInclkwqlTpxAeHg6NRuPv6hAREZEIgiCgoaEBPXv2hFbruN+HgVAHTp06hfj4eH9Xg4iIiNzw888/o3fv3g6fZyDUgfDwcABtDRkREeGVexiNRmzYsAHZ2dkIDg72yj0CHdtQGmxHabAdPcc2lIaa29FgMCA+Pt7yOe4IA6EOmIfDIiIivBoIhYWFISIiQnVvVKmwDaXBdpQG29FzbENpsB3RYVoLk6WJiIhItRgIERERkWopKhDasmUL8vLy0LNnT2g0Gnz66acdnrN582YMGTIEOp0OV111Fd577z2v15OIiIiUQVGBUGNjI1JSUrB48WJR5SsrK3H77bfj1ltvxe7du/Hoo49i6tSp+PLLL71cUyIiIlICRSVLjxw5EiNHjhRdfunSpUhMTMSCBQsAAAMGDMC3336L119/HTk5Od6qJhERESmEogIhV5WXl2PEiBFWx3JycvDoo486PKe5uRnNzc2WxwaDAUBb5r3RaPRKPc3X9db11YBtKA22ozTYjp5jG0pDze0o9jUHdCBUXV2N2NhYq2OxsbEwGAz49ddf0blzZ5tziouLMXfuXJvjGzZsQFhYmNfqCgBlZWVevb4asA2lwXaUBtvRc2xDaaixHZuamkSVC+hAyB2FhYUoKCiwPDYvyJSdne3VdYTKysqQlZWl2nUePMU2lAbbURpsR8+xDaWh5nY0j+h0JKADIb1ej5qaGqtjNTU1iIiIsNsbBAA6nQ46nc7meHBwsNffRL64R6BjG0qD7SgNtqPn2IbSUGM7in29AR0IZWZmoqSkxOpYWVkZMjMz/VQjIqJLWk0CKirrUNtwATHhoUhPjEaQlps7E/mSogKhX375BUeOHLE8rqysxO7duxEdHY0+ffqgsLAQJ0+exAcffAAAmDFjBt5880088cQT+MMf/oBNmzbh//7v/7B+/Xp/vQQiIgBA6d4qzF23H1X1FyzH4iJDUZSXhNzkOD/WjEhdFLWO0Pbt2zF48GAMHjwYAFBQUIDBgwdjzpw5AICqqiqcOHHCUj4xMRHr169HWVkZUlJSsGDBArz77rucOk9EflW6twozV+y0CoIAoLr+Amau2InSvVV+qhmR+iiqR2jYsGEQBMHh8/ZWjR42bBh27drlxVoREYnXahIwd91+2PtLJgDQAJi7bj+ykvQcJiPyAUX1CBERKV1FZZ1NT1B7AoCq+guoqKzzXaWIVIyBEBGRD9U2OA6C3ClHRJ5hIERE5EMx4aGSliMizzAQIiLyofTEaMRFhsJR9o8GbbPH0hOjfVktItViIERE5ENBWg2K8pIAwCYYMj8uyktiojSRjzAQIiLysdzkOCyZMAT6SOvhL31kKJZMGMJ1hIh8SFHT54mIAkVuchyykvRcWZrIzxgIEZHbuEWEZ4K0GmT26+bvahCpGgMhInILt4ggokDAHCEichm3iCCiQMFAiIhc0tEWEUDbFhGtJsfb4RARyQUDISJyCbeIIKJAwkCIiFzCLSKIKJAwECIil3CLCCIKJAyEiMgl3CKCiAIJAyEicgm3iCCiQMJAiIhcxi0iiChQcEFFInILt4ggokDAQIiI3MYtIohI6Tg0RkRERKrFQIiIiIhUi4EQERERqRYDISIiIlItBkJERESkWgyEiIiISLUYCBEREZFqMRAiIiIi1WIgRERERKrFQIiIiIhUi4EQERERqRYDISIiIlItBkJERESkWgyEiIiISLUUFwgtXrwYCQkJCA0NRUZGBioqKpyWX7hwIa699lp07twZ8fHxmDVrFi5cuOCj2hIREZGcKSoQWr16NQoKClBUVISdO3ciJSUFOTk5qK2ttVt+5cqVmD17NoqKinDgwAH87W9/w+rVq/HnP//ZxzUnIiIiOVJUIPTaa69h2rRpmDx5MpKSkrB06VKEhYVh2bJldst/9913uOGGGzBu3DgkJCQgOzsb9913X4e9SERERKQOnfxdAbFaWlqwY8cOFBYWWo5ptVqMGDEC5eXlds+5/vrrsWLFClRUVCA9PR0//fQTSkpKcP/99zu8T3NzM5qbmy2PDQYDAMBoNMJoNEr0aqyZr+ut66sB21AabEdpsB09xzaUhprbUexrVkwgdObMGbS2tiI2NtbqeGxsLA4ePGj3nHHjxuHMmTO48cYbIQgCLl68iBkzZjgdGisuLsbcuXNtjm/YsAFhYWGevYgOlJWVefX6asA2lAbbURpsR8+xDaWhxnZsamoSVU4xgZA7Nm/ejBdffBF//etfkZGRgSNHjuCRRx7B888/j2eeecbuOYWFhSgoKLA8NhgMiI+PR3Z2NiIiIrxST6PRiLKyMmRlZSE4ONgr9wh0bENpsB2lwXb0HNtQGmpuR/OITkcUEwh1794dQUFBqKmpsTpeU1MDvV5v95xnnnkG999/P6ZOnQoAuO6669DY2Ijp06fjqaeeglZrmyKl0+mg0+lsjgcHB3v9TeSLewQ6tqE02I7SYDt6jm0oDTW2o9jXq5hk6ZCQEAwdOhQbN260HDOZTNi4cSMyMzPtntPU1GQT7AQFBQEABEHwXmWJiIhIERTTIwQABQUFmDhxItLS0pCeno6FCxeisbERkydPBgA88MAD6NWrF4qLiwEAeXl5eO211zB48GDL0NgzzzyDvLw8S0BERERE6qWoQGjMmDE4ffo05syZg+rqaqSmpqK0tNSSQH3ixAmrHqCnn34aGo0GTz/9NE6ePIkePXogLy8P8+bN89dLICIiIhlRVCAEAPn5+cjPz7f73ObNm60ed+rUCUVFRSgqKvJBzYiIiEhpFJMjRERERCQ1BkJERESkWgyEiIiISLUYCBEREZFqMRAiIiIi1WIgRERERKrFQIiIiIhUS3HrCBERkTitJgEVlXWobbiAmPBQpCdGI0ir8Xe1iGSFgRARUQAq3VuFuev2o6r+guVYXGQoivKSkJsc58eaEckLh8aIiAJM6d4qzFyx0yoIAoDq+guYuWInSvdW+almRPLDQIiIKIC0mgTMXbcfgp3nzMfmrtuPVpO9EkTqw0CIiCiAVFTW2fQEtScAqKq/gIrKOt9VikjGGAgREQWQ2gbHQZA75YgCHQMhIqIAEhMeKmk5okDHQIiIKICkJ0YjLjIUjibJa9A2eyw9MdqX1SKSLQZCREQBJEirQVFeEgDYBEPmx0V5SVxPiOi/GAgREQWY3OQ4LJkwBPpI6+EvfWQolkwYwnWEiNrhgopERAEoNzkOWUl6rixN1AEGQkREASpIq0Fmv27+rgaRrHFojIiIiFSLgRARERGpFgMhIiIiUi0GQkRERKRaDISIiIhItRgIERERkWoxECIiIiLVYiBEREREqsVAiIiIiFSLgRARERGpFgMhIiIiUi0GQkRERKRaDISIiIhItRgIERERkWoxECIiIiLVUlwgtHjxYiQkJCA0NBQZGRmoqKhwWv78+fN48MEHERcXB51Oh2uuuQYlJSU+qi0RERHJWSd/V8AVq1evRkFBAZYuXYqMjAwsXLgQOTk5OHToEGJiYmzKt7S0ICsrCzExMVizZg169eqF48ePIyoqyveVJyIiItlRVCD02muvYdq0aZg8eTIAYOnSpVi/fj2WLVuG2bNn25RftmwZ6urq8N133yE4OBgAkJCQ4PQezc3NaG5utjw2GAwAAKPRCKPRKNErsWa+rreurwZsQ2mwHaXBdhSn1SRg+/FzqG1oRky4Dml9uyJIqwHANpSKmttR7GvWCIIgeLkukmhpaUFYWBjWrFmD0aNHW45PnDgR58+fx9q1a23OGTVqFKKjoxEWFoa1a9eiR48eGDduHJ588kkEBQXZvc+zzz6LuXPn2hxfuXIlwsLCJHs9RERqtuesBp8c0+J8i8ZyLCpEwN0JJqR0U8THEslcU1MTxo0bh/r6ekRERDgsp5geoTNnzqC1tRWxsbFWx2NjY3Hw4EG75/z000/YtGkTxo8fj5KSEhw5cgR/+tOfYDQaUVRUZPecwsJCFBQUWB4bDAbEx8cjOzvbaUN6wmg0oqysDFlZWZaeK3IN21AabEdpsB2d+3JfDZaX78Hl4U59iwbLDwdh0dgU3HZNNNtQAmp+L5pHdDqimEDIHSaTCTExMXj77bcRFBSEoUOH4uTJk3j11VcdBkI6nQ46nc7meHBwsNffRL64R6BjG0qD7SgNtqOtVpOAeV8csgmCAEAAoAEw74tDGDHgJgBsQ6mosR3Fvl7FBELdu3dHUFAQampqrI7X1NRAr9fbPScuLg7BwcFWw2ADBgxAdXU1WlpaEBIS4tU6ExGRtYrKOlTVX3D4vACgqv4Cth8/57tKkaopZvp8SEgIhg4dio0bN1qOmUwmbNy4EZmZmXbPueGGG3DkyBGYTCbLscOHDyMuLo5BEBGRH9Q2OA6CrMs1d1yISAKKCYQAoKCgAO+88w7ef/99HDhwADNnzkRjY6NlFtkDDzyAwsJCS/mZM2eirq4OjzzyCA4fPoz169fjxRdfxIMPPuivl0BEpGox4aEiy9mmKBB5g2KGxgBgzJgxOH36NObMmYPq6mqkpqaitLTUkkB94sQJaLWXYrv4+Hh8+eWXmDVrFgYNGoRevXrhkUcewZNPPumvl0BEpGrpidGIiwxFdf0Fu3lCGgD6yFCk9e2KLw/4unakRooKhAAgPz8f+fn5dp/bvHmzzbHMzExs27bNy7UiIiIxgrQaFOUlYeaKndAAVsGQeSJ9UV6SZT0hIm9T1NAYEREpX25yHJZMGAJ9pPUwmT4yFEsmDEFucpyfakZqpLgeISIiUr7c5DhkJelRUVmH2oYLiAkPRXpiNHuCyOcYCBERkV8EaTXI7NfN39UglePQGBEREakWAyEiIiJSLQZCREREpFoMhIiIiEi1GAgRERGRajEQIiIiItViIERERESqxUCIiIiIVIuBEBEREakWAyEiIiJSLQZCREREpFoMhIiIiEi1GAgRERGRajEQIiIiItViIERERESqxUCIiIiIVKuTvytARETkTKtJQEVlHWobLiAmPBTpidEI0mr8XS0bSqknWWMgREREsvXlvhrM++IQquovWI7FRYaiKC8JuclxfqyZtdK9VZi7br/s60m2ODRGRESytOesBg+t2mMVXABAdf0FzFyxE6V7q/xUM2ule6swc8VO2deT7GMgREREstNqEvDJMS0EO8+Zj81dtx+tJnslfKfVJGDuuv2yryc5xkCIiIhkZ/vxczjf4ji/RgBQVX8BFZV1vquUHRWVdTY9Qe3JpZ7kGAMhIiKSndqGZpHlHAchviD2/v6uJznGQIiIiGQnJlwnslyol2sizf39XU9yjIEQERHJTlrfrogKEeBocEyDtllZ6YnRvqyWjfTEaMRFhsq+nuQYAyEiIpKdIK0GdyeYAMAmyDA/LspL8vs6PUFaDYrykgDIu57kGAMhIiKSpZRuAhaNTYE+0npYSR8ZiiUThshmfZ7c5DgsmTBE9vUk+7igIpGCcSVbCnQ5A2MxclAv2b/Pc5PjkJWkl309yRYDISKF4kq2pBZBWg0y+3XzdzU6pJR6kjUOjREpEFeyJSKSBgMhIoXhSrZERNJRXCC0ePFiJCQkIDQ0FBkZGaioqBB13qpVq6DRaDB69GjvVpDIy7iSLRGRdBQVCK1evRoFBQUoKirCzp07kZKSgpycHNTW1jo979ixY/if//kf3HTTTT6qKZH3cCVbIiLpKCoQeu211zBt2jRMnjwZSUlJWLp0KcLCwrBs2TKH57S2tmL8+PGYO3currzySh/Wlsg7uJIt+UKrSUD50bNYu/skyo+e5VArBSzFzBpraWnBjh07UFhYaDmm1WoxYsQIlJeXOzzvueeeQ0xMDKZMmYJvvvmmw/s0NzejufnSHjcGgwEAYDQaYTQaPXgFjpmv663rq4Ga2nBw73DoI3SoMTTbzRPSANBH6jC4d7jL7aGmdvQmpbfjl/tq8ELJQVQbLv0t1Efo8PSo/sgZGOuTOii9DeVCze0o9jUrJhA6c+YMWltbERtr/UsYGxuLgwcP2j3n22+/xd/+9jfs3r1b9H2Ki4sxd+5cm+MbNmxAWFiYS3V2VVlZmVevrwZqacNReg2WGcwduu3XKREgABgZ24QvS79w+/pqaUdvU2I77jmrwbLDtu+tasMF5K/ajT9cY0JKN9/1DimxDeVIje3Y1NQkqpxiAiFXNTQ04P7778c777yD7t27iz6vsLAQBQUFlscGgwHx8fHIzs5GRESEN6oKo9GIsrIyZGVlITg42Cv3CHRqa8NRAIbY+dYeFxmKp0a6/61dbe3oLUptx1aTgOIFWwDY2/ldAw2AL2rC8MT4m72+UKBS21Bu1NyO5hGdjigmEOrevTuCgoJQU1NjdbympgZ6vd6m/NGjR3Hs2DHk5eVZjplMbfvWdOrUCYcOHUK/fv1sztPpdNDpbHc9Dg4O9vqbyBf3CHRqasM7Unt7bcVdNbWjNymtHbcfPWsVWF+ubUZiM3b9p8FnCwcqrQ3lSo3tKPb1KiYQCgkJwdChQ7Fx40bLFHiTyYSNGzciPz/fpnz//v3x73//2+rY008/jYaGBvzlL39BfHy8L6pN5FVcyZakxBmJpEaKCYQAoKCgABMnTkRaWhrS09OxcOFCNDY2YvLkyQCABx54AL169UJxcTFCQ0ORnJxsdX5UVBQA2BwnIvIWJe0HxxmJpEaKCoTGjBmD06dPY86cOaiurkZqaipKS0stCdQnTpyAVquoFQGIKIB9ua8G87445LP94DwNutIToxEXGYrq+gtOZiS2XZcoUCgqEAKA/Px8u0NhALB582an57733nvSV4iIyI49ZzVYXr7HJqAw7we3ZMIQSYMhKTbhDdJqUJSXhJkrdkIDWNXdHE4V5SXJtkeLyB3sPiEiklirScAnx7Q+2w9Oyk14c5PjsGTCEOgjrYe/9JGhkgdvRHKguB4hIiK52378HM63OO41ab8fnKfJ7h1twqtBW9CVlaQX3ZOTmxyHrCS9YnKbiDzBQIiISGK1DY6noFuX83z2lSub8LoSdF0+I9G85QYDIwo0DISIiCQWE267Fpn9cp7PvvLFlHcp8o+I5Io5QkREEkvr2xVRIQIc9Zdo0BZISDH7yttT3qXMPyKSIwZCREQSC9JqcHdC20r2lwdDUs++Mk9590bQ1VH+ESBt0jeRPzAQIiLygpRuAhaNTfH67CvzlHdA+qDLlfwjIqVijhARkZfkDIz12n5w7ZmnvF+ex6P3MI+HW26QGjAQIiLyIl/tB+eNKe/ccoPUgIEQEVGAkDro4pYbpAbMESIiIru8mX9EJBeie4Q+++wz0Rf93e9+51ZliIhIXryVf0QkF6IDodGjR4sqp9Fo0Nra6m59iIhIZrjlBgUy0YGQyWTyZj2IiEjGfJX0TeRrzBEiIiIi1XJ71lhjYyO+/vprnDhxAi0tLVbPPfzwwx5XjIiIiMjb3AqEdu3ahVGjRqGpqQmNjY2Ijo7GmTNnEBYWhpiYGAZCREREpAhuDY3NmjULeXl5OHfuHDp37oxt27bh+PHjGDp0KObPny91HYmIiIi8wq1AaPfu3Xjssceg1WoRFBSE5uZmxMfH45VXXsGf//xnqetIRERE5BVuBULBwcHQattOjYmJwYkTJwAAkZGR+Pnnn6WrHRHJSqtJQPnRs1i7+yTKj57lruNEpHhu5QgNHjwYP/zwA66++mrccsstmDNnDs6cOYMPP/wQycnJUteRiGSgdG+VzaJ6cVxUj4gUzq0eoRdffBFxcW1/+ObNm4euXbti5syZOH36NN5++21JK0hE/le6twozV+y0CoIAoLr+Amau2InSvVV+qhkRkWfc6hFKS0uz/DsmJgalpaWSVYiIfK/VJKCisg5V5xvxU70GrSYBwe2em7tuv91NNwW07Tk1d91+ZCXpudIwESkOd58nUjnbIa8grFmwBc/+biByk+PaAqTLeoLaEwBU1V9ARWUdVx4mIsVxKxBKTEyERuP4m99PP/3kdoWIyHfMQ16X9/bUGJoxc8VOLJkwBM0XxW2vU9vgOFgiIpIrtwKhRx991Oqx0WjErl27UFpaiscff1yKehGRl4kd8pr/+xRR14sJD5WyekREPuFWIPTII4/YPb548WJs377dowoRkW+IHfKCpm12WHX9BbtBkwaAPrJtN3IiIqWRdNPVkSNH4uOPP5bykkTkJWKHss780oyivCQAbUFPe+bHRXlJTJQmIkWSNBBas2YNoqP5rZBICcQOZcWEhyI3OQ5LJgyBPtL6HH1kKJZMGMJ1hIhIsdxeULF9srQgCKiursbp06fx17/+VbLKEZH3pCdGuzTklZsch6wkPSoq61DbcAEx4W3PsSeIiJTMrUBo9OjRVo+1Wi169OiBYcOGoX///lLUi4i8LEirQVFeEmau2AkNYBUMORryCtJqOEWeiAKKW4FQUVGR1PUgIj8wD3ldvnWGPlKHoryBHPIiooAnOhAyGAyiLxoREeFWZYjI99oPeVWdb8RP+3Yjf8zNCNWF+LtqREReJzpZOioqCl27dhX1nzctXrwYCQkJCA0NRUZGBioqKhyWfeedd3DTTTdZ6jVixAin5YnUyjzklTcoDldHCsz7IQowrSYB5UfPYu3ukyg/ehatJnuZgeokukfoq6++svz72LFjmD17NiZNmoTMzEwAQHl5Od5//30UFxdLX8v/Wr16NQoKCrB06VJkZGRg4cKFyMnJwaFDhxATE2NTfvPmzbjvvvtw/fXXIzQ0FC+//DKys7Oxb98+9OrVy2v1JCIikos9ZzUoXrAF1YZmy7G4yFAU5SVx+BsuBEK33HKL5d/PPfccXnvtNdx3332WY7/73e9w3XXX4e2338bEiROlreV/vfbaa5g2bRomT54MAFi6dCnWr1+PZcuWYfbs2TblP/roI6vH7777Lj7++GNs3LgRDzzwgFfqSEREJBdf7qvBssNaAM1Wx6vrL1i20VF7MORWsnR5eTmWLl1qczwtLQ1Tp071uFL2tLS0YMeOHSgsLLQc02q1GDFiBMrLy0Vdo6mpCUaj0elaR83NzWhuvvSGMedGGY1GGI1GN2vvnPm63rq+GrANpcF2lAbb0XNsQ8+1mgQ8X3LQ7nOXttHZh2FXdwvI4XCx7x23AqH4+Hi88847eOWVV6yOv/vuu4iPj3fnkh06c+YMWltbERsba3U8NjYWBw/a/0Ff7sknn0TPnj0xYsQIh2WKi4sxd+5cm+MbNmxAWFiYa5V2UVlZmVevrwZsQ2mwHaXBdvQc29B9P9ZrUGMIgu2a8G3attFpxpurS3F1ZODlDDU1NYkq51Yg9Prrr+Oee+7BF198gYyMDABARUUFfvzxR9lusfHSSy9h1apV2Lx5M0JDHa+oW1hYiIKCAstjg8GA+Ph4ZGdne202nNFoRFlZGbKyshAcHOyVewQ6tqE02I7SYDt6jm3ouXX/qgL2/7vDclcOTMWoQYE3PCZ2trtbgdCoUaNw+PBhLFmyxNIbk5eXhxkzZnitR6h79+4ICgpCTU2N1fGamhro9Xqn586fPx8vvfQS/vnPf2LQoEFOy+p0Ouh0OpvjwcHBXv9l9MU9Ah3bUBpsR2mwHT3HNnRfXFQX0eUCsY3Fvia3AiGgbXjsxRdfdPd0l4WEhGDo0KHYuHGjZWVrk8mEjRs3Ij8/3+F5r7zyCubNm4cvv/wSaWlpPqotERGRf6UnRkMfoUO14QLsDY9dvo2OWokOhP71r38hOTkZWq0W//rXv5yW7ajXxV0FBQWYOHEi0tLSkJ6ejoULF6KxsdEyi+yBBx5Ar169LFP4X375ZcyZMwcrV65EQkICqqurAQBXXHEFrrjiCq/UkYiISA6CtBo8Pao/8lftdrqNDgCUHz2r2j0ERQdCqampqK6uRkxMDFJTU6HRaCAItslVGo0Gra2tklbSbMyYMTh9+jTmzJmD6upqpKamorS01JJAfeLECWi1l9aIXLJkCVpaWvD73//e6jpFRUV49tlnvVJHIiIiucgZGIs/XGNCSXWY1TpC+v+uIwQAN768yWqLHbWtMSQ6EKqsrESPHj0s//aX/Px8h0Nhmzdvtnp87Ngx71eIiLyq1SRwx3siD6R0E/DE+Jux6z8NVr9HZfurMXPFTlzepaG2NYZEB0J9+/a1+28iIm8p3VtlsyGs2r6tEknBvI2OWatJwNx1+22CIKD9GkP7kZWkD/gvHqL3Gmvv/fffx/r16y2Pn3jiCURFReH666/H8ePHJascEalX6d4qzFyx0yoIAi59Wy3dW+WnmhEpX0Vlnc3vVnttawxdQEVlne8q5SduBUIvvvgiOnfuDKBtlek333wTr7zyCrp3745Zs2ZJWkEiUp+Ovq0Cbd9WuXEkkXtqGxwHQe6UUzK3ps///PPPuOqqqwAAn376KX7/+99j+vTpuOGGGzBs2DAp60dETgRq/owr31bbd/eT7wTqe08tYsIdLyzsTjklcysQuuKKK3D27Fn06dMHGzZssKzEHBoail9//VXSChKRfYGcP8Nvq/IWyO89s0AP9NIToxEXGYrq+gt2e17VtMaQW4FQVlYWpk6disGDB+Pw4cMYNWoUAGDfvn1ISEiQsn5EZIc5fyZQZ3vw26p8Bfp7D1BHoBek1aAoLwkzV+x0usZQIAV/jriVI7R48WJkZmbi9OnT+Pjjj9GtW1vX9I4dO3DfffdJWkEisqak/JlWk4Dyo2exdvdJlB89K7pO5m+rjv4Ea9D2waSGb6uecLf9nV1PKe89d6kpST83OQ5LJgyBPtL6C4U+MjQgAlqx3OoRioqKwptvvmlz3N6u7UQkLaXkz3jyrZrfVj3njV4Npbz33KXGKeW5yXHIStIH9DBgR9zqEQKAb775BhMmTMD111+PkydPAgA+/PBDfPvtt5JVjohsKSF/Ropv1fy26j5v9Woo4b3nCbVOKTevMXRnai9k9uumqiAIcLNH6OOPP8b999+P8ePHY+fOnWhublu2u76+Hi+++CJKSkokrSQRXSL3/Bkpv1UHwrdVXyfderNXQ6r3nlwTkQM90CP73AqEXnjhBSxduhQPPPAAVq1aZTl+ww034IUXXpCsckRkS+6zPaQePrl8RVwl+XJfDeZ9ccinSbfeHL6S4r0n50RkuX/JIO9wa2js0KFDuPnmm22OR0ZG4vz5857WiYicMOfPALBJJpZD/gy/VbfZc1aDh1bt8XnSrTfb39P3ntwTkZmkr05uBUJ6vR5HjhyxOf7tt9/iyiuv9LhSROScnPNn+K26bejnk2Nav8yu8nb7u/veU8KMM7l/ySDvcGtobNq0aXjkkUewbNkyaDQanDp1CuXl5XjssccwZ84cqetIRHbINX9G7kN3vrD9+Dmcb3H8c/Dm7CpftL877z2lzDgzB3qXD9/pZTJ8R9JzKxCaPXs2TCYThg8fjqamJtx8883Q6XR4/PHHMXXqVKnrSEQOyDF/RuzUdwAoP3pWVkGcVGobmkWWk3540FdLD7j63lPSkKlcv2SQd7g1NKbRaPDUU0+hrq4Oe/fuxbZt23D69GlERkYiMTFR6joSkcJ0NHwCADe+vAn3vbMNj6zajfve2YYbX97k9xwRsTpaqDAmXCfqOt4aHpTj0KnShkzVPqVcTVzqEWpubsazzz6LsrIySw/Q6NGjsXz5ctx1110ICgri7vNEBMDxt+qy/dWK3qJBzKyntL5dERUioL5F47fhQbn1agTCkKlcp/2TZ1wKhObMmYO33noLI0aMwHfffYd7770XkydPxrZt27BgwQLce++9CAoK8lZdiUhhLh8+UfrKvWL32QrSanB3ggnLDwf5dWVsOQ2dKn21cDlP+yfPuDQ09ve//x0ffPAB1qxZgw0bNqC1tRUXL17Enj17MHbsWAZBRAFIyv2qlLxyr6uznlK6CVg0NkVWw1P+JschOzHkPu2fPONSj9B//vMfDB06FACQnJwMnU6HWbNmQaORZwRPRB1r393fLawT2sc5Un8LVlLC7OVcCeLS+kQAAHIGxmLkoF4cTmlHbkN2HVF6LyZ1zKVAqLW1FSEhIZdO7tQJV1xxheSVIiLfsBfoRIUEITihBp06BUmey6OEhFlHeSCuBXERlsdyGp6SCyW1iVKm/ZP7XAqEBEHApEmToNO1zYi4cOECZsyYgS5duliV++STT6SrIRF5haN8l/MtwEOr9iAyLFjyb8FyT5h11gOmhCCOpCc2AC7bX81ASKFcyhGaOHEiYmJiEBkZicjISEyYMAE9e/a0PDb/R0SecZSXI1W+jrPufqBtptP5JqPD893N5ZHzyr0d5YGca2zh9gsqJDawXbb1GHOFFMqlHqHly5d7qx5E9F8l/zqFp9fuRV3jpUAkLjIUv0uJw2d7qiTJ1+mou18sd3J55Lhyr5g8kOfX78cztw/Agyt3dTjrydTq7RqTr5h7MTv6fWGukHK5tbI0EXlHccl+vLWl0uZ4Vf0Fu8fdzdeRKhnZk/2qOkqY9eWaLWLzQLp20ckuiCPvMvdizlix02k55gopFwMhIpko+VeV3WDHGXfzdTzNY5Eil8dZwqyv12xxJRH6ztReipr1RJ7LTY7DlBsS8LetxzosK8cZj+ScW1tsEJG0Wk0Cnl67161z3cnXMXf3O8t36RoWbPn35c8B3svl8ceaLa4mQnP7BfUZkaQXVY7J8srDQIj8RsqF+pSuorIOdY0tHl3DlW+izpKWzdkvxXdfh6U+XvzO1UULpSImMGQitLrxPRK4ODRGfsHl6q1J0Z3u6jdRR0nLUSHAC3enWH4OvhwG8teaLUrf/oG8j++RwMVAiHxO7H5N/uTrzRU97U5395vo5UnL3cI64fT+bcgZGGsp48vF7/y58rQcZ7ORvPA9EpgYCPnB5R+yg3uH+7tKPqOE5er90VsldoquI8/cPsDt9mof6BiNRpQccOsykvD3ooVK2/6BfI/vkcDDQMjH7H3I6iN0GKXXYJQf6+WI1D0jcl+u3l+9Ve273d3JfunaRSd5nfxBDitPK2n7BzF83bupBoH2HlE7BkI+5OhDtsbQjGUGLYbsq8Edqb39Ujd7vNEzIudNN/3dW+Wo212MQJmyyzwMccQGN8zFI+qY4maNLV68GAkJCQgNDUVGRgYqKiqclv/73/+O/v37IzQ0FNdddx1KSkp8VFNrYmbDzPvioGxmTnlrCrO/hz6ccaW3yltyk+Pw7ZO34X+n/RZ/GZuKZ24fIOq8QJqyaw4IfTlbTUlK91bhxpc34b53tuGRVbtx3zvbcOPLm2x+J/2xDAGREimqR2j16tUoKCjA0qVLkZGRgYULFyInJweHDh1CTEyMTfnvvvsO9913H4qLi3HHHXdg5cqVGD16NHbu3Ink5GSf1r3jLQ00qKpvlsWqpN7sGZHD0Icj/uytcvQNv9Uk4N1vK2XZXq5wdXiGeRj2iR269XfvJpEYchm2VVQg9Nprr2HatGmYPHkyAGDp0qVYv349li1bhtmzZ9uU/8tf/oLc3Fw8/vjjAIDnn38eZWVlePPNN7F06VKf1l3OQ0KX82Yej5yHPvzVW9XR8IVc20ssd4dnmIdhzZXgRu65eERyGrZVTCDU0tKCHTt2oLCw0HJMq9VixIgRKC8vt3tOeXk5CgoKrI7l5OTg008/dXif5uZmNDc3Wx4bDAYAbbNpjEbHu3F3pFuYuKbuGhqEbw/XoLahGTHhOqT17erzD7mq842iyxmNES5ff/i13bFobApeKDmIasOlttZH6vDUyP4Yfm13l9vaXN6Tn9Hg3uHQR+hQY2h20vuiw+De4R7dp70v99XgoVV7HH7DXzQ2BTkDYyVvL0ekaMf2xL6+QCN1OwLA9yKDm/IjtahtaHZYrj13f4d9wWg0wiQAW3+sRd2vrX77e6h03ngvespXfxfEvmbFBEJnzpxBa2srYmOtGyc2NhYHDx60e051dbXd8tXV1Q7vU1xcjLlz59oc37BhA8LCwtyoeRuTAESFBOF8C2BvLV9AQFgn4KGV21Hfcun5qBABdyeYkNLNd7lDP9VrAAR1XG7fbpT8Z5fb93kyCThq0MBgBCKCgX4RjWg9vgMlx92+JMrKytw/GcAovQbLDObUufY/JwECgNTwJry4ovS/9RXgyd9kkwDM3Rn03z8G1hcS/vv/pz/ZDeOxVmg13mkvRzxtR8D11xeIpGhHsx1nxP1ebvjme0QEQ1RZT3+H3WESLn8f2/892nNWg0+OBeH8tt2WY/74exgopHwvesKXfxeamppElVNMIOQrhYWFVr1IBoMB8fHxyM7ORkSEZ9+cghPaomDAdohDANB00fanXt+iwfLDQR5HyK0mAduPnxPV09RqErBmwZYOe0byx9wsm29nRqMRZWVlyMrKQnBwsNvXGQVgyL4am96XqLBgQNDgi/9c+oahj9Dh6VH93f65fF9Zh/PbtjspocH5FqBH0m+R4WIOkCs/7/akakfAu69P7qRsR7NulXX44Edn7dkm+6YMpPXtKsvf4S/31aD48p5NO79HX+6rwfLyPf/9+nGJVH8P1cQb70VP+PLvgnlEpyOKCYS6d++OoKAg1NTUWB2vqamBXm9/Mzy9Xu9SeQDQ6XTQ6WzXZAkODvb4TXRHam906hRkMy4aGxECQ1Mzmi7anmMe+5/3xSGMHNTLrT9aro7FBgN49ncDO8hLGYhQXYjLdRHL3SQ6qX5OIwf1stz/2JkmLPznYbvLHjy0ao/bM5nO2vuBOyjnymuSYuxdinb01utTEina0SzzqhhREw0yr4pBkFbj99/hy5XurbI7HHL571GrScC8Lw457DHw9O+hWkn5XvSEL/8uiD1fMdPnQ0JCMHToUGzcuNFyzGQyYePGjcjMzLR7TmZmplV5oK170FF5X7h8evT/TvstXr77Oru9QWaeTNt2dwqtP6cwi50e7E3mRN07BvXEqh9OeGUTUG8kZ8tpyrScl0pQImcb5dpLnJfTMgSubKYrh2UsyHvk+HdBMT1CAFBQUICJEyciLS0N6enpWLhwIRobGy2zyB544AH06tULxcXFAIBHHnkEt9xyCxYsWIDbb78dq1atwvbt2/H222/782XYzIb5ZMcJUee5OqPM0ym0vpjCfHnPz7nGFjy4Uj77kHlz9o3USwnIbcq0nJdKUCpX97qSyzIErvweKWmGLblOjn8XFBUIjRkzBqdPn8acOXNQXV2N1NRUlJaWWhKiT5w4Aa32UifX9ddfj5UrV+Lpp5/Gn//8Z1x99dX49NNPfb6GUEdiwsVtj+BqhCzFh7g3pzDbG8LRaiCbD3LAu8seSL2UgNif93tbK9E9XOfyh6Krw5VyXipByVwNbuSwDIErv0dy7DEg6cjx74KiAiEAyM/PR35+vt3nNm/ebHPs3nvvxb333uvlWnkmrW9XRIUIqG/RSBohy/mblaOF4ZyNMEmx9on5w7zacAF1vzQjuksI9JGdHX6QePuPspS7WYv9OT6//tKuqmJzh9zNO+Ju3d4hh+DGFa78Hsmxx4CkJbe/C4oLhAJRkFaDuxNMWH44SNIIWa7frJwN4YjhbuBm78PczNGHui/+KEs1fOHOz1HMkKOnG9HKZXiG/MeV36P2PQaX+oIvlQPYkxgI5PR3QTHJ0oEupZuARWNTJE1sNP/xcfS20qAtAPD1N6uOtxtxzp0PfEdJxGZVDpKJXU1QdZf5G/6dqb2Q2a+bW9fr6OdtT0cJ364kuTojxesj5XIn0XvR2BREXTapjfvNBRa5/F1gj5CM5AyMtZq27WmELMexWMD9Hh13e1/E9kAJsJ+DJLduXEec/bydaT/kmNbHeq0sbtVAUnH19yhnYCyMx1rRI+m3ONt0kT2J5DUMhGRG6rF/OX6Iu9Oj40ng5koPlKMPdTl14zrj6OctRluAGmHnmNhz5U0uGzyqmau/R1oNkJEYLYv1byhwMRBSAbl9iHeULwC0/QFsP9riSeDm6oe0o/JKSVDNTY7Dbf1j8WH5MRyva4IgCPhwW8dLNNgLUOWaZ+YqOW3wqHZK+T0i9WAgpBJy+uMjZsjuzfsGo2sXnSSBm6sf0nL/UO+Io2UJHKXxtB9yNLVar/oaCDN4PE32DnTsKSO1YyBEfuHLITsxPVBm/kgel5KryxJcPuRoarV+Xq55ZmLJbZFJuWFPGREDIfIjXw3ZWU/HdUwDeX+od0RMUrg7Q47+yjOToqdCbLL3tp/O4oaruntYY2VRak8Ze7BIagyEyK98NWTXURJxIHwLFpMUbhKAZ24f4PLK0r7OM5Oqp0Jsftj0D7bj1d8PwqhBPV2uqxL5s6fMk0CGPVjkDQyESDXaf5iLXVlaScR+6HcP1+HO1F4uX99XQauUPRVi870aW1rxp5W78Mf/nEfhqCTRdVVq74S/lkXwJJBRag8WyR8DIVIVOSWNSy0QZnhJ3VPhSn4YALy1pRIpvbti1KCOP1CdfagPv1bew2z+WBbBk0CGuV7kTVxZmihAyHUlcVe40lMhRvsVjcV6Zu3eDlfKdrRSuflD/ct9NS7d09d8HTR7ukK51O8LovYYCBEFCHe3A2k1CSg/ehZrd5/E95V1Tje+9TZv9FSY88OiOotblO9sY4vTD1QxH+rzvjjo13bsiK+DZk8DmUBa2JPkh0NjRAHE1Rle9oZ3okKCEJxQgztSe/us3mbe6qnITY5DuC4Y4//2vajyzj5QxX2oN+OowbdDNK7kK/l6WQRPA5lAGPYl+WIgRBRgxM7wcpSzcb4FeGjVHnTqFCRZ8qnYD2lvLuD4237dEN0lGHWNxg7LOvtAFfuhbuj4NpJxJwnZl8sieBrIBMLCniRfDISIAlBHSeHO1xxqC1CkSj515UPamz0VQVoNXrgzGX9auctpuY6GhMR+qEf4aHssT5KQfbUsgqeBjNIX9iR5Y44QkQr5Kvm0o6Ti0r1VNueYeyr0kdYBhz4y1OMp0qMG9cQfb050+LyYRTXF5dfo0C/C+0lCniYhA5eC5jtTeyGzXzevBBPu5q+15833Bakbe4SIVMgXyaeeTHl2taei5aLJssls3+gw3J+ZgJBO9r/nFY5KQkrvrnh67V7UNbZYjotdz0ZM78RTI/uj9fgOp9eRgr/WA3KHFENxcttAmgIDAyEiFfJF8qmnH9Ji13wqLtmPd76ptJqlNa/kAKbdlOhwccRRg+KQk+z+B2pHH+rDr+2OkuOiLuURpc2mkiKQCeS1wMg/GAgRqZAvkk998SFdXLIfb22ptDluEmA57igY8vQD1dmHutHom0xpJc6mYiBDcsMcISIVcpazYR7s8TT51Nsf0i0XTXjnG9sgqL13vqlEy0WTW9cXwxf5Nc4EwiKaRP7GQIhIpRwln0aFAIvGptjN2Wi/+GL50bNOk3C9/SH9YfmxDhctNAlt5Tzhymv2NSmSkInUjkNjRCp2+fBOt7BOOL1/G3IGxtqUdXWtGm9PeT5e1yRpOXuUsNu5L9cDIgpEDISIVK59zobRaETJAdsy7q5V480P6b7RYZKWu5ySdjvnbCoi9zEQIiKnPN3521sf0vdnJmBeyQGnw2NaTVs5Vylxt3MmIRO5hzlCROSUFIsveiOpOKSTFtNucrw4IgBMuynR4XpCzqhht3M55z4R+RJ7hIjIKTmvVWOeGn/5OkJaDZyuI9QROb9mKSgh94nIVxgIEZFTcl+rpnBUEh7L7i96ZWkx5P6aPaGk3CciX2AgREROKWHn75BOWky56UrJrqeE1+wOJeY+EXkbc4SIyCk1rlXjz9fszdwdNeQ+EbmKPUJE1CE5rlXTahK8Ol3cH6/Z27k7gZ77ROQOBkJEJIqc1qrxVbKvL1+zL3J3Ajn3ichdDISISDRfrFXTUU+Pr5N9ffWafZG7E6i5T0SeYCBERLLRUU9PoCb7upK740lQ5u1tT4iUSDHJ0nV1dRg/fjwiIiIQFRWFKVOm4JdffnFa/qGHHsK1116Lzp07o0+fPnj44YdRX1/vw1oTkVjmnp7LAwJzT0/p3qqATfb1Ze6Oo8129ZGhnDpPqqSYHqHx48ejqqoKZWVlMBqNmDx5MqZPn46VK1faLX/q1CmcOnUK8+fPR1JSEo4fP44ZM2bg1KlTWLNmjY9rT0TOiO3peSLnWlHXU1qyr69zd+SU70Xkb4oIhA4cOIDS0lL88MMPSEtLAwAsWrQIo0aNwvz589GzZ0+bc5KTk/Hxxx9bHvfr1w/z5s3DhAkTcPHiRXTqZP+lNzc3o7m52fLYYDAAaNuM0mg0SvmyLMzX9db11YBtKA1/teP3Int6Tjf8Kup63cI6efQaWk0Cth8/h9qGZsSE65DWt6tLQYKr7Ti4dzj0ETrUGJqd5O7oMLh3uKQ/m7Q+EQAiAACm1oswtUp2aY/xd1oaam5Hsa9ZEYFQeXk5oqKiLEEQAIwYMQJarRbff/897rrrLlHXqa+vR0REhMMgCACKi4sxd+5cm+MbNmxAWJh7u1iLVVZW5tXrqwHbUBq+bscdZzQAgjos9/OPBxAVosX5FsB2hR8AEBAVApzevw0lB9yry56zGnxyTIvzLZeuHxUi4O4EE1K6ubamjyvtOEqvwTKDOVuh/WsTIAAYGduEL0u/cOn+gYC/09JQYzs2NTWJKqeIQKi6uhoxMTFWxzp16oTo6GhUV1eLusaZM2fw/PPPY/r06U7LFRYWoqCgwPLYYDAgPj4e2dnZiIiIcL3yIhiNRpSVlSErKwvBwcFeuUegYxtKw1/t2K2yDh/8uL3Dcjk3ZyCjyYiHVu0BYC/ZV4MX7k5BzsBYt+rx5b4aLC/fY9MrU9+iwfLDQVg0Vty13WnHUQCG7KvBCyUHUW241CsdFxmKp0b2d/s1KRV/p6Wh5nY0j+h0xK+B0OzZs/Hyyy87LXPggJtf69oxGAy4/fbbkZSUhGeffdZpWZ1OB51OZ3M8ODjY628iX9wj0LENpeHrdsy8KkbUtO7Mq2IQpNWgU6cgyRc6bDUJmPfFIad5SvO+OISRg3qJHiZztR3vSO2NkYN6MXenHf5OS0ON7Sj29fo1EHrssccwadIkp2WuvPJK6PV61NbWWh2/ePEi6urqoNfrnZ7f0NCA3NxchIeH4x//+Ifq3ghESuDqtG5vJPv6agp7R3yxbhERXeLXQKhHjx7o0aNHh+UyMzNx/vx57NixA0OHDgUAbNq0CSaTCRkZGQ7PMxgMyMnJgU6nw2effYbQUK6WSiRXrm5pIXXAwO0niNRJETlCAwYMQG5uLqZNm4alS5fCaDQiPz8fY8eOtcwYO3nyJIYPH44PPvgA6enpMBgMyM7ORlNTE1asWAGDwWAZL+zRoweCgjpOzCQi3/LntG5uP0GkTooIhADgo48+Qn5+PoYPHw6tVot77rkHb7zxhuV5o9GIQ4cOWbLEd+7cie+//x4AcNVVV1ldq7KyEgkJCT6rOxGJ56+hIW4/4Rlvb4JL5C2KCYSio6MdLp4IAAkJCRCES3++hg0bZvWYiMgZb28/EciBgq82wSXyBsUEQkRE3uZqnpJYgRwo+HoTXCKpMRAiImpH6jylQA4UAnUTXFIXBkJERJeRKk+po0ABUHagIJclB4g8oZjd54mIlGb78XNOAwWgLVB4c9OPPqqRtLjkAAUC9ggREXlJbUNzx4UAvP7PH3GtPlxxQ2RccsC3Ajnh3p8YCBEReUlMuO12PY4ocYiMSw74TiAn3Psbh8aIiLwkrW9XxEWK6w0x59IoiXnJAeDSEgNmUiw5QG3MCfeXD7OaE+5L91b5qWaBgYEQEZGXtA8UxFBiLo15yQH9ZQGfPjJU0TPi5EJswn2rievmuYtDY0REXpSbHIdZI67B6/883GFZpebS+HNrlEDHmXnex0CIiMjL8m+7Cv9bcRzVBvvJ04GQS+OvrVECHWfmeR+HxoiIvCxIq8GzvxsIDZhLQ67hzDzvYyBEROQDzKUhd5hn5jkKkTVomz2m5N5Ef+PQGBGRjzCXhlzl7c2AiYEQEZFPMZeGXOWtzYCpDQMhIiIimWNvovcwECIiIlIA9iZ6B5OliYiISLUYCBEREZFqcWiMiMhF3AWcKHAwECIicgF3AScKLBwaIyISydEu4FX1FzBjxU6U/OuUV+7bahJQfvQs1u4+ifKjZ7nBJpGE2CNERCSCs13AzfL/dxfehAZZA7pLdl/2QBF5F3uEiIhE6GgXcAAwCcCfVu7El/tqJLmnox6o6voLmLliJ0r3VklyHyI1YyBERCSCK7t7z/viIDwdvXLWA2U+Nnfdfg6TEXmIgRARkQiu7O5dVd+MowbPZpF11AMloC03qaKyzqP7EKkdAyEiIhHMu4CLZTB6dj+xPVCu9FQRkS0GQkREIph3ARcrItiz+4ntgXKlp4qIbDEQIiISKTc5Dn8dNxjO1k7UAIiL1KFfhGe5O+YeKEe3artP22KOROQ+BkJERC4YNagn3rxviN3nzEHLUyP7Ow2WxGjfA3X5pcyPi/KSuKI1kYcYCBERuWjUoDgsnTDEJmdIHxmKJROGIGdgrCT3yU2Ow5IJQ6B3cB+uI0TkOS6oSETkhtzkOGQl6e3uOWY0epgpLfI+ROQ5BkJERG4K0mqQ2a+bIu/DjWN9j20uTwyEiIhUhtt2+B7bXL4UkyNUV1eH8ePHIyIiAlFRUZgyZQp++eUXUecKgoCRI0dCo9Hg008/9W5FiYhkjNt2+B7bXN4UEwiNHz8e+/btQ1lZGT7//HNs2bIF06dPF3XuwoULodGw+5GI1I3bdvge21z+FBEIHThwAKWlpXj33XeRkZGBG2+8EYsWLcKqVatw6tQpp+fu3r0bCxYswLJly3xUWyIieeK2Hb7HNpc/ReQIlZeXIyoqCmlpaZZjI0aMgFarxffff4+77rrL7nlNTU0YN24cFi9eDL1eL+pezc3NaG5utjw2GAwAAKPRKOlMkPbM1/XW9dWAbSgNtqM05NqOVecbRZczGiO8XBvn5NqGrvJ3mwdKO7pD7GtWRCBUXV2NmJgYq2OdOnVCdHQ0qqurHZ43a9YsXH/99bjzzjtF36u4uBhz5861Ob5hwwaEhYWJr7QbysrKvHp9NWAbSoPtKA25teNP9RoAQR2X27cbJf/Z5f0KiSC3NnSVXNpc6e3ojqamJlHl/BoIzZ49Gy+//LLTMgcOHHDr2p999hk2bdqEXbtce2MVFhaioKDA8thgMCA+Ph7Z2dmIiPDONySj0YiysjJkZWUhONjDDYpUim0oDbajNOTajq0mAWsWbEGNodluzooGgD5Sh/wxN/t9Wrdc29BV/m7zQGlHd5hHdDri10Dosccew6RJk5yWufLKK6HX61FbW2t1/OLFi6irq3M45LVp0yYcPXoUUVFRVsfvuece3HTTTdi8ebPd83Q6HXQ6nc3x4OBgr7+JfHGPQMc2lAbbURpya8dgAM/+biBmrtgJDWD1wXxp246BCNWF+L5yDrjahnJbq0cubS6396IviH29fg2EevTogR49enRYLjMzE+fPn8eOHTswdOhQAG2BjslkQkZGht1zZs+ejalTp1odu+666/D6668jLy/P88oTESmQeduOy9e00QfAmjZyXasnkNs8ECgiR2jAgAHIzc3FtGnTsHTpUhiNRuTn52Ps2LHo2bMnAODkyZMYPnw4PvjgA6Snp0Ov19vtLerTpw8SExN9/RKIiGQjELftMK/Vc/nwk3mtHn/vzRaIbR4oFBEIAcBHH32E/Px8DB8+HFqtFvfccw/eeOMNy/NGoxGHDh0SnRxFRKRmvtoexBc6WqtHg7a1erKS9H4NPAKpzQOJYgKh6OhorFy50uHzCQkJEATnC1J19DwRESmPK2v1MBChyyliQUUiIiJHahscB0HulCN1YSBERESKFhMeKmk5UhcGQkREpGjpidGIiwyFo+wfDdpmj6UnRvuyWqQQDISIiEjRgrQaFOUlAYBNMHRprZ4kztAiuxgIERGR4pnX6tFHWg9/6SND/T51nuRNMbPGiIiInOFaPeQOBkJERBQwuFYPuYpDY0RERKRaDISIiIhItRgIERERkWoxR4iISIVaTQKTionAQIiISHVK91Zh7rr9VvtzxUWGoigvidPMSXU4NEZEpCKle6swc8VOm01Kq+svYOaKnSjdW+WnmllrNQn4sV6Ddf+qQvnRs2g1cdNs8g72CBERqUSrScDcdfthL6QQ0LYK89x1+5GVpPfrMFnp3io8+9k+VBuCgP3/BsAeK/Ie9ggREalERWWdTU9QewKAqvoLqKis812lLmPusao2NFsdl1uPFQUOBkJERCpR2+A4CHKnnNQ66rEC2nqsOExGUmIgRESkEjHhoR0XcqGc1JTQY0WBh4EQEZFKpCdGIy4y1GaHdjMN2nJx0hOjfVktC7n3WFFgYiBERKQSQVoNivKSAMAmGDI/LspL8luitNx7rCgwMRAiIlKR3OQ4LJkwBPpI62BCHxmKJROG+HVWltx7rCgwcfo8EZHK5CbHIStJL7uVpc09VjNX7IQGsEqalkOPFQUmBkJERCoUpNUgs183f1fDhrnHqm0doUtT6PVcR4i8hIEQERFZyGEPstzkOAy7uhveXF2KKwemIi6qiyx6rCgwMRAiIiIA8tqDLEirwdWRAkYNikNwcLBP703qwmRpIiJSzB5kRFJjIEREpHJc0ZnUjIEQEZHKcUVnUjMGQkREKscVnUnNGAgREakcV3QmNWMgRESkclzRmdSMgRARkcrJfQ8yIm9iIERERLLeg4zIm7igIhERAZDvHmRE3qSYHqG6ujqMHz8eERERiIqKwpQpU/DLL790eF55eTluu+02dOnSBREREbj55pvx66+/+qDGRETKY96D7M7UXsjs141BEAU8xQRC48ePx759+1BWVobPP/8cW7ZswfTp052eU15ejtzcXGRnZ6OiogI//PAD8vPzodUq5mUTERGRFyliaOzAgQMoLS3FDz/8gLS0NADAokWLMGrUKMyfPx89e/a0e96sWbPw8MMPY/bs2ZZj1157rU/qTERERPKniECovLwcUVFRliAIAEaMGAGtVovvv/8ed911l805tbW1+P777zF+/Hhcf/31OHr0KPr374958+bhxhtvdHiv5uZmNDc3Wx4bDAYAgNFohNFolPBVXWK+rreurwZsQ2mwHaXBdvQc21Aaam5Hsa9ZEYFQdXU1YmJirI516tQJ0dHRqK6utnvOTz/9BAB49tlnMX/+fKSmpuKDDz7A8OHDsXfvXlx99dV2zysuLsbcuXNtjm/YsAFhYWEevhLnysrKvHp9NWAbSoPtKA22o+fYhtJQYzs2NTWJKufXQGj27Nl4+eWXnZY5cOCAW9c2mUwAgD/+8Y+YPHkyAGDw4MHYuHEjli1bhuLiYrvnFRYWoqCgwPLYYDAgPj4e2dnZiIiIcKsuHTEajSgrK0NWVhaCg4O9co9AxzaUBttRGmxHz7ENpaHmdjSP6HTEr4HQY489hkmTJjktc+WVV0Kv16O2ttbq+MWLF1FXVwe9Xm/3vLi4tjUvkpKSrI4PGDAAJ06ccHg/nU4HnU5nczw4ONjrbyJf3CPQsQ2lwXaUBtvRc2xDaaixHcW+Xr8GQj169ECPHj06LJeZmYnz589jx44dGDp0KABg06ZNMJlMyMjIsHtOQkICevbsiUOHDlkdP3z4MEaOHOl55YmIiEjxFDGPfMCAAcjNzcW0adNQUVGBrVu3Ij8/H2PHjrXMGDt58iT69++PiooKAIBGo8Hjjz+ON954A2vWrMGRI0fwzDPP4ODBg5gyZYo/Xw4RERHJhCKSpQHgo48+Qn5+PoYPHw6tVot77rkHb7zxhuV5o9GIQ4cOWSVHPfroo7hw4QJmzZqFuro6pKSkoKysDP369fPHSyAiIiKZUUwgFB0djZUrVzp8PiEhAYIg2ByfPXu21TpCrjJfU2zSlTuMRiOamppgMBhUN4YrFbahNNiO0mA7eo5tKA01t6P5c9tebNCeYgIhf2loaAAAxMfH+7kmRERE5KqGhgZERkY6fF4jdBQqqZzJZMKpU6cQHh4OjcY7e+6Yp+j//PPPXpuiH+jYhtJgO0qD7eg5tqE01NyOgiCgoaEBPXv2dLq1FnuEOqDVatG7d2+f3CsiIkJ1b1SpsQ2lwXaUBtvRc2xDaai1HZ31BJkpYtYYERERkTcwECIiIiLVYiAkAzqdDkVFRXZXtCZx2IbSYDtKg+3oObahNNiOHWOyNBEREakWe4SIiIhItRgIERERkWoxECIiIiLVYiBEREREqsVASGYOHz6MO++8E927d0dERARuvPFGfPXVV/6uliKtX78eGRkZ6Ny5M7p27YrRo0f7u0qK1NzcjNTUVGg0Guzevdvf1VGUY8eOYcqUKUhMTETnzp3Rr18/FBUVoaWlxd9Vk73FixcjISEBoaGhyMjIQEVFhb+rpBjFxcX4zW9+g/DwcMTExGD06NE4dOiQv6slWwyEZOaOO+7AxYsXsWnTJuzYsQMpKSm44447UF1d7e+qKcrHH3+M+++/H5MnT8aePXuwdetWjBs3zt/VUqQnnngCPXv29Hc1FOngwYMwmUx46623sG/fPrz++utYunQp/vznP/u7arK2evVqFBQUoKioCDt37kRKSgpycnJQW1vr76opwtdff40HH3wQ27ZtQ1lZGYxGI7Kzs9HY2OjvqsmTQLJx+vRpAYCwZcsWyzGDwSAAEMrKyvxYM2UxGo1Cr169hHfffdffVVG8kpISoX///sK+ffsEAMKuXbv8XSXFe+WVV4TExER/V0PW0tPThQcffNDyuLW1VejZs6dQXFzsx1opV21trQBA+Prrr/1dFVlij5CMdOvWDddeey0++OADNDY24uLFi3jrrbcQExODoUOH+rt6irFz506cPHkSWq0WgwcPRlxcHEaOHIm9e/f6u2qKUlNTg2nTpuHDDz9EWFiYv6sTMOrr6xEdHe3vashWS0sLduzYgREjRliOabVajBgxAuXl5X6smXLV19cDAN93DjAQkhGNRoN//vOf2LVrF8LDwxEaGorXXnsNpaWl6Nq1q7+rpxg//fQTAODZZ5/F008/jc8//xxdu3bFsGHDUFdX5+faKYMgCJg0aRJmzJiBtLQ0f1cnYBw5cgSLFi3CH//4R39XRbbOnDmD1tZWxMbGWh2PjY1lioAbTCYTHn30Udxwww1ITk72d3VkiYGQD8yePRsajcbpfwcPHoQgCHjwwQcRExODb775BhUVFRg9ejTy8vJQVVXl75fhd2Lb0WQyAQCeeuop3HPPPRg6dCiWL18OjUaDv//9735+Ff4ltg0XLVqEhoYGFBYW+rvKsiS2Hds7efIkcnNzce+992LatGl+qjmpzYMPPoi9e/di1apV/q6KbHGLDR84ffo0zp4967TMlVdeiW+++QbZ2dk4d+4cIiIiLM9dffXVmDJlCmbPnu3tqsqa2HbcunUrbrvtNnzzzTe48cYbLc9lZGRgxIgRmDdvnrerKlti2/D//b//h3Xr1kGj0ViOt7a2IigoCOPHj8f777/v7arKmth2DAkJAQCcOnUKw4YNw29/+1u899570Gr5HdSRlpYWhIWFYc2aNVYzPSdOnIjz589j7dq1/qucwuTn52Pt2rXYsmULEhMT/V0d2erk7wqoQY8ePdCjR48OyzU1NQGAzR9JrVZr6eVQM7HtOHToUOh0Ohw6dMgSCBmNRhw7dgx9+/b1djVlTWwbvvHGG3jhhRcsj0+dOoWcnBysXr0aGRkZ3qyiIohtR6CtJ+jWW2+19EwyCHIuJCQEQ4cOxcaNGy2BkMlkwsaNG5Gfn+/fyimEIAh46KGH8I9//AObN29mENQBBkIykpmZia5du2LixImYM2cOOnfujHfeeQeVlZW4/fbb/V09xYiIiMCMGTNQVFSE+Ph49O3bF6+++ioA4N577/Vz7ZShT58+Vo+vuOIKAEC/fv3Qu3dvf1RJkU6ePIlhw4ahb9++mD9/Pk6fPm15Tq/X+7Fm8lZQUICJEyciLS0N6enpWLhwIRobGzF58mR/V00RHnzwQaxcuRJr165FeHi4JbcqMjISnTt39nPt5IeBkIx0794dpaWleOqpp3DbbbfBaDRi4MCBWLt2LVJSUvxdPUV59dVX0alTJ9x///349ddfkZGRgU2bNjHpnHyqrKwMR44cwZEjR2wCSGYlODZmzBicPn0ac+bMQXV1NVJTU1FaWmqTQE32LVmyBAAwbNgwq+PLly/HpEmTfF8hmWOOEBEREakWB6uJiIhItRgIERERkWoxECIiIiLVYiBEREREqsVAiIiIiFSLgRARERGpFgMhIiIiUi0GQkRERKRaDISISBIajQaffvqpv6shC5s3b4ZGo8H58+cBAO+99x6ioqL8Wiciso+BEBGJMmnSJKvdwC9XVVWFkSNH+q5CCjJmzBgcPnxYdPmEhAQsXLjQexUiIgvuNUZEkpDDJqKCIKC1tRWdOnn+p03Ka3Xu3Nkvm122tLQgJCTE5/clUhL2CBGRJNoPjR07dgwajQaffPIJbr31VoSFhSElJQXl5eVW53z77be46aab0LlzZ8THx+Phhx9GY2Oj5fkPP/wQaWlpCA8Ph16vx7hx41BbW2t53jwE9cUXX2Do0KHQ6XT49ttvbepmrs+qVatw/fXXIzQ0FMnJyfj66687vJbJZEJxcTESExPRuXNnpKSkYM2aNVbXLykpwTXXXIPOnTvj1ltvxbFjx6yetzc0tm7dOvzmN79BaGgounfvjrvuugtA20aZx48fx6xZs6DRaKDRaCznfPzxxxg4cCB0Oh0SEhKwYMECq2smJCTg+eefxwMPPICIiAhMnz7dwU+LiCwEIiIRJk6cKNx5550Onwcg/OMf/xAEQRAqKysFAEL//v2Fzz//XDh06JDw+9//Xujbt69gNBoFQRCEI0eOCF26dBFef/114fDhw8LWrVuFwYMHC5MmTbJc829/+5tQUlIiHD16VCgvLxcyMzOFkSNHWp7/6quvBADCoEGDhA0bNghHjhwRzp49a1M3c3169+4trFmzRti/f78wdepUITw8XDhz5ozTa73wwgtC//79hdLSUuHo0aPC8uXLBZ1OJ2zevFkQBEE4ceKEoNPphIKCAuHgwYPCihUrhNjYWAGAcO7cOUEQBGH58uVCZGSkpT6ff/65EBQUJMyZM0fYv3+/sHv3buHFF18UBEEQzp49K/Tu3Vt47rnnhKqqKqGqqkoQBEHYvn27oNVqheeee044dOiQsHz5cqFz587C8uXLLdft27evEBERIcyfP184cuSIcOTIEXE/XCIVYyBERKK4Ewi9++67luf37dsnABAOHDggCIIgTJkyRZg+fbrVNb755htBq9UKv/76q917/PDDDwIAoaGhQRCES8HLp59+6rTu5vq89NJLlmNGo1Ho3bu38PLLLzu81oULF4SwsDDhu+++s7relClThPvuu08QBEEoLCwUkpKSrJ5/8sknnQZCmZmZwvjx4x3Wt2/fvsLrr79udWzcuHFCVlaW1bHHH3/c6t59+/YVRo8e7fC6RGSLQ2NE5DWDBg2y/DsuLg4ALENbe/bswXvvvYcrrrjC8l9OTg5MJhMqKysBADt27EBeXh769OmD8PBw3HLLLQCAEydOWN0nLS1NVH0yMzMt/+7UqRPS0tJw4MABh9c6cuQImpqakJWVZVXPDz74AEePHgUAHDhwABkZGQ7vY8/u3bsxfPhwUXU2O3DgAG644QarYzfccAN+/PFHtLa22q0/EXWMydJE5DXBwcGWf5tzXUwmEwDgl19+wR//+Ec8/PDDNuf16dMHjY2NyMnJQU5ODj766CP06NEDJ06cQE5ODlpaWqzKd+nSRbI6t7/WL7/8AgBYv349evXqZVVOp9O5fQ9vJk5L2RZEasBAiIj8YsiQIdi/fz+uuuoqu8//+9//xtmzZ/HSSy8hPj4eALB9+3aP7rlt2zbcfPPNAICLFy9ix44dyM/Pd1g+KSkJOp0OJ06csPRGXW7AgAH47LPPbO7jzKBBg7Bx40ZMnjzZ7vMhISFWvTzm+2zdutXq2NatW3HNNdcgKCjI6f2IyDEGQkQkWn19PXbv3m11rFu3bpZAxRVPPvkkfvvb3yI/Px9Tp05Fly5dsH//fpSVleHNN99Enz59EBISgkWLFmHGjBnYu3cvnn/+eY/qv3jxYlx99dUYMGAAXn/9dZw7dw5/+MMfHJYPDw/H//zP/2DWrFkwmUy48cYbUV9fj61btyIiIgITJ07EjBkzsGDBAjz++OOYOnUqduzYgffee89pPYqKijB8+HD069cPY8eOxcWLF1FSUoInn3wSQNvsry1btmDs2LHQ6XTo3r07HnvsMfzmN7/B888/jzFjxqC8vBxvvvkm/vrXv3rUJkSq5+8kJSJShokTJwoAbP6bMmWKIAj2k6V37dplOf/cuXMCAOGrr76yHKuoqBCysrKEK664QujSpYswaNAgYd68eZbnV65cKSQkJAg6nU7IzMwUPvvsM6vrmhOczUnJjpjrs3LlSiE9PV0ICQkRkpKShE2bNlnKOLqWyWQSFi5cKFx77bVCcHCw0KNHDyEnJ0f4+uuvLWXWrVsnXHXVVYJOpxNuuukmYdmyZU6TpQVBED7++GMhNTVVCAkJEbp37y7cfffdlufKy8uFQYMGCTqdTmj/Z3rNmjVCUlKSEBwcLPTp00d49dVXra5pL8maiJzTCIIg+CMAIyLylWPHjiExMRG7du1Camqqv6tDRDLCWWNERESkWgyEiIiISLU4NEZERESqxR4hIiIiUi0GQkRERKRaDISIiIhItRgIERERkWoxECIiIiLVYiBEREREqsVAiIiIiFSLgRARERGp1v8HKhkGWxtRsjYAAAAASUVORK5CYII=", "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": "2025-12-05T18:07:55.189941Z", "iopub.status.busy": "2025-12-05T18:07:55.189631Z", "iopub.status.idle": "2025-12-05T18:07:55.195301Z", "shell.execute_reply": "2025-12-05T18:07:55.193864Z" }, "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": "2025-12-05T18:07:55.198225Z", "iopub.status.busy": "2025-12-05T18:07:55.197954Z", "iopub.status.idle": "2025-12-05T18:07:55.238962Z", "shell.execute_reply": "2025-12-05T18:07:55.237438Z" } }, "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: Fri, 05 Dec 2025 Deviance: 7.2134\n", "Time: 18:07:55 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": "2025-12-05T18:07:55.243059Z", "iopub.status.busy": "2025-12-05T18:07:55.241811Z", "iopub.status.idle": "2025-12-05T18:07:55.404888Z", "shell.execute_reply": "2025-12-05T18:07:55.403192Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/hostedtoolcache/Python/3.10.19/x64/lib/python3.10/site-packages/statsmodels/genmod/generalized_linear_model.py:985: FutureWarning: linear keyword is deprecated, use which=\"linear\"\n", " warnings.warn(msg, FutureWarning)\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'Residual')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGwCAYAAABRgJRuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQQRJREFUeJzt3Xl81NW9//H3JCaTBJIgS0iUACmiGCOgIBi1FpewVZTaelVKCypauaAWvC5c/bHUBW1r8V61aLGCSrn6EBdEMUpdwAWkLEEDKIYGsZKwBEkwSBgy398f6QSGZDKTycx8l3k9Hw8fMjPfzHzmzOT7/eSczznHZRiGIQAAAJtLMDsAAACASCCpAQAAjkBSAwAAHIGkBgAAOAJJDQAAcASSGgAA4AgkNQAAwBFOMDuAWPJ6vdq5c6fS09PlcrnMDgcAAITAMAwdOHBAJ510khISAvfHxFVSs3PnTuXm5podBgAACMM333yjbt26BXw8rpKa9PR0SQ2NkpGRYXI0R3k8Hr3zzjsaOnSokpKSzA7HVmi78NF24aPtwkfbhSfe262mpka5ubmN1/FA4iqp8Q05ZWRkWC6pSUtLU0ZGRlx+WduCtgsfbRc+2i58tF14aLcGwUpHKBQGAACOQFIDAAAcgaQGAAA4AkkNAABwBJIaAADgCCQ1AADAEUhqAACAI5DUAAAARyCpAQAAjhBXKwoDsJ56r6E15fu0+8AhZaWnaFBeRyUmsOEsgNYjqQFgmuLSCs1aulkV1Yca78vJTNGMUfkaXpBjYmQA7Mg2w09z585V3759G/dtKiws1FtvvWV2WADCVFxaoYkL1/slNJJUWX1IExeuV3FphUmRAbAr2yQ13bp100MPPaR169Zp7dq1uvjii3XFFVdo06ZNZocGoJXqvYZmLd0so5nHfPfNWrpZ9d7mjgCA5tkmqRk1apRGjhyp3r1769RTT9UDDzyg9u3ba/Xq1WaHBqCV1pTva9JDcyxDUkX1Ia0p3xe7oADYni1raurr6/XSSy+ptrZWhYWFAY+rq6tTXV1d4+2amhpJDVu4ezyeqMcZKl8sVorJLmi78JnZdhX7a0M+zuPJiHI0rcf3Lny0XXjivd1Cfd8uwzBs07/7+eefq7CwUIcOHVL79u21aNEijRw5MuDxM2fO1KxZs5rcv2jRIqWlpUUzVAAt+Krapcc3JwY9bnJ+vXpn2uYUBSBKDh48qDFjxqi6uloZGYH/0LFVUnP48GHt2LFD1dXVWrx4sZ5++mmtWLFC+fn5zR7fXE9Nbm6u9u7d22KjxJrH49Hy5ctVVFSkpKQks8OxFdoufGa2Xb3X0JBHVmpXTV2zdTUuSdmZbr0/9UJLTu/mexc+2i488d5uNTU16ty5c9CkxlbDT8nJyTrllFMkSQMGDNA//vEP/c///I+eeuqpZo93u91yu91N7k9KSrLkl8KqcdkBbRc+M9ouSdLMy8/QxIXr5ZL8EhtfCjNj1BlKcSfHNK7W4nsXPtouPPHabqG+Z9sUCjfH6/X69cQAsI/hBTmaO/ZsZWem+N2fnZmiuWPPZp0aAK1mm56aadOmacSIEerevbsOHDigRYsW6YMPPtDbb79tdmgAwjS8IEdF+dmsKAwgImyT1OzevVu//vWvVVFRoczMTPXt21dvv/22ioqKzA4NQBskJrhU2KuT2WEAcADbJDV//etfzQ4BAABYmK1ragAAAHxIagAAgCOQ1AAAAEcgqQEAAI5AUgMAAByBpAYAADgCSQ0AAHAEkhoAAOAIJDUAAMARSGoAAIAjkNQAAABHIKkBAACOQFIDAAAcgaQGAAA4AkkNAABwBJIaAADgCCQ1AADAEUhqAACAI5DUAAAARyCpAQAAjkBSAwAAHIGkBgAAOAJJDQAAcASSGgAA4AgnmB0AADhdvdfQmvJ92n3gkLLSUzQor6MSE1xmhwU4DkkNAERRcWmFZi3drIrqQ4335WSmaMaofA0vyDExMsB5GH4CgCgpLq3QxIXr/RIaSaqsPqSJC9eruLTCpMgAZyKpAYAoqPcamrV0s4xmHvPdN2vpZtV7mzsCQDhIagAgCtaU72vSQ3MsQ1JF9SGtKd8Xu6AAhyOpAYAo2H0gcEITznEAgiOpAYAoyEpPiehxAIIjqQGAKBiU11E5mSkKNHHbpYZZUIPyOsYyLMDRSGoAIAoSE1yaMSpfkpokNr7bM0bls14NEEEkNQAQJcMLcjR37NnKzvQfYsrOTNHcsWezTg0QYSy+BwBRNLwgR0X52awoDMQASQ0ARFligkuFvTqZHQbgeAw/AQAARyCpAQAAjkBSAwAAHIGkBgAAOAJJDQAAcASSGgAA4AgkNQAAwBFIagAAgCOw+B4AOFS912AlY8QVkhoAcKDi0grNWrpZFdWHGu/LyUzRjFH57DkFx2L4CQAcpri0QhMXrvdLaCSpsvqQJi5cr+LSCpMiA6KLpAYAHKTea2jW0s0ymnnMd9+spZtV723uCMDeSGoAwEHWlO9r0kNzLENSRfUhrSnfF7uggBghqQEAB9l9IHBCE85xgJ2Q1ACAg2Slp0T0OMBOSGoAwEEG5XVUTmaKAk3cdqlhFtSgvI6xDAuICZIaAHCQxASXZozKl6QmiY3v9oxR+axXA0ciqQEAhxlekKO5Y89Wdqb/EFN2Zormjj2bdWrgWCy+BwAONLwgR0X52awojLhCUgMADpWY4FJhr05mhwHEDMNPAADAEUhqAACAIzD8BMSB5nZrBgCnIakBHC7Qbs33jDjNxKgAIPIYfgIcrKXdmm95YaM2VjETBoBzkNQADhXKbs2vbE9gt2YAjkFSAzhUKLs17z/s0tqvv4tdUAAQRSQ1gEOFvltzXZQjAYDYsE1SM3v2bJ1zzjlKT09XVlaWRo8erS+//NLssADLCn23ZneUIwGA2LBNUrNixQpNmjRJq1ev1vLly+XxeDR06FDV1taaHRpgSaHs1twh2dDAHifGMiwAiBrbTOkuLi72u71gwQJlZWVp3bp1uvDCC5v9mbq6OtXVHe1ar6mpkSR5PB55PJ7oBdtKvlisFJNd0HYtu2fEabrlhY1ySX4Fw75E58qeXnnrj8jjYRZUa/C9Cx9tF554b7dQ37fLMAxbTn0oKytT79699fnnn6ugoKDZY2bOnKlZs2Y1uX/RokVKS0uLdoiAJWyscumV7Qnaf/ho4tIh2dCVPb3q18mWv/4A4szBgwc1ZswYVVdXKyMjI+BxtkxqvF6vLr/8cu3fv18fffRRwOOa66nJzc3V3r17W2yUWPN4PFq+fLmKioqUlJRkdji2QtuFpt5raO3X32n3gTplpbs1sMeJ8tYfoe3CxPcufLRdeOK93WpqatS5c+egSY1thp+ONWnSJJWWlraY0EiS2+2W2920CDIpKcmSXwqrxmUHtF3LkiRdcGpXv/t8Q060Xfhou/DRduGJ13YL9T3bLqmZPHmy3njjDa1cuVLdunUzOxwAAGARtklqDMPQLbfcoldffVUffPCB8vLyzA4JAABYiG2SmkmTJmnRokVasmSJ0tPTVVlZKUnKzMxUamqqydEBAACz2Wadmrlz56q6ulpDhgxRTk5O438vvvii2aEBAAALsE1PjQ0naQEAgBiyTU8NAABAS0hqAACAI5DUAAAARyCpAQAAjkBSAwAAHIGkBgAAOAJJDQAAcASSGgAA4AgkNQAAwBFIagAAgCOQ1AAAAEcgqQEAAI5AUgMAAByBpAYAADgCSQ0AAHAEkhoAAOAIJ5gdAADrqvcaWlO+T7sPHFJWeooG5XVUYoLL7LAAoFkkNQCaVVxaoVlLN6ui+lDjfTmZKZoxKl/DC3JMjAw+JJ2AP5IaAH7qvYYef69Mc/6+tcljldWHNHHhes0dezaJjclIOoGmqKkB0Ki4tELnP/RuswmNJBn//v+spZtV7zWaPQbRV1xaoYkL1/slNNLRpLO4tMKkyABzkdQAkHT0QllZU9ficYakiupDWlO+LzaBwU+919CspZvVXEpJ0ol4R1IDoMULZSC7DxwKfhAibk35viY9NMci6UQ8I6kBEPRC2Zys9JQoRYOWhJpMknQiHlEoDKBVF0CXpOzMhpk2iL1Qk0mSTsQjemoAtPoCOGNUPlOHTTIor6NyMlMUqPVdapgFRdKJeERSAyDohdInJzOF6dwmS0xwacaofElq8nn5bpN0Il6R1ABo8ULpM+XS3vrorotJaCxgeEGO5o49W9mZ/j1s2SSdiHPU1ACQdPRCyYJu9jC8IEdF+dmsKAwcg6QGQCMulPaSmOBSYa9OZocBWAZJDQA/XCgB2BU1NQAAwBHoqQEA2Aq7kyMQkhoAgG2wOzlawvATAMAW2J0cwZDUAAAsj93JEQqSGgCA5bE7OUJBUgMAsDx2J0coSGoAAJbH7uQIBUkNAMDy2J0coSCpAQBYHruTIxQkNQAQonqvoVXbqrT0swp9Ve1ipk2MsTs5gmHxPQAIQdNF3xK1+JGVmnn5GVxMY4hNV9ESkhoACMK36Nvx/TK7auo0ceF6eglijE1XEQjDTwDQAhZ9A+yDpAYAWsCib4B9kNQAQAtY9A2wD5IaAGgBi74B9kFSAwAtYNE3wD5IagCgBSz6BtgHSQ0ABBF40Tc307kBC2GdGgAIwbGLvlXsr9U/N5Vo8tUXKsWdbHZoAP6NpAYAQuRb9M3jydCyf21gyAmwGIafAACAI5DUAAAAR2D4CUDU1XsNNiAEEHUkNQCiqunu1g3ruswYlc+sIQARxfATgKjx7W59/N5JldWHNHHhehWXVpgUGQAnIqkBEBXsbg0g1khqAEQFu1sDiDVqagAHskJhLrtbA4i1kJOa119/PeQnvfzyy8MKBkDbWaUwl92tAcRayEnN6NGjQzrO5XKpvr4+3HgAhKi53pjlmys1ceH6JnUsvsLcWO5T5NvdurL6ULN1NS5J2exuDSCCQk5qvF5vNOMA0ArN9cZkZ6To0JH6gIW5LjUU5hblZ8dkKMq3u/XEhevlkvziisbu1lYYcgNgLlsVCq9cuVKjRo3SSSedJJfLpddee83skICYCzhNuuaQ9h/0BPw5MwpzA+9unRLRXqPi0gpd8PB7unbeat32QomunbdaFzz8HlPGgTgTdqFwbW2tVqxYoR07dujw4cN+j916661tDizQa/br10/XX3+9rrzyyqi8BmBlLU2TDlWsC3OP3d06Gr0oviTPCkNuAMwVVlKzYcMGjRw5UgcPHlRtba06duyovXv3Ki0tTVlZWVFLakaMGKERI0ZE5bkBOwg2TToUZhTm+na3jrRga+HEesgNgLnCSmqmTJmiUaNG6cknn1RmZqZWr16tpKQkjR07VrfddlukYwxbXV2d6urqGm/X1NRIkjwejzyewN30seaLxUox2UW8tV3F/tqwf7ahMNets7ql+/0O2LntPg1xLZxVZbs1OIIFyU5oO7PQduGJ93YL9X27DMNodU92hw4d9Omnn+q0005Thw4dtGrVKp1++un69NNPNW7cOH3xxRetDri1XC6XXn311RZnZc2cOVOzZs1qcv+iRYuUlpYWxeiA6Piq2qXHNyeGcKSvn+LY29L1p3rVr5NzVvBdt9el574K3h6/7l2vAZ2d876BeHPw4EGNGTNG1dXVysjICHhcWD01SUlJSkhoqDHOysrSjh07dPrppyszM1PffPNNeBFHwbRp0zR16tTG2zU1NcrNzdXQoUNbbJRY83g8Wr58uYqKipSUlGR2OLYSb21X7zW0+JGV2lVTF3CadGZaktyJLu06cLTWLSczRfeM6KNhZ3RtvM8JbdepfJ+e+2pt0OOG/nhwxHtq7N52ZqHtwhPv7eYbaQkmrKTmrLPO0j/+8Q/17t1bP/nJTzR9+nTt3btXzz//vAoKCsJ5yqhwu91yu91N7k9KSrLkl8KqcdlBvLRdkqSZl5/R4jTph648s1WFuXZuu8JTskJaC6fwlKyo1NTYue3MRtuFJ17bLdT3HNaU7gcffFA5OQ2zCR544AGdeOKJmjhxovbs2aO//OUv4TwlgBCFMk3aV5h7Rf+TVdirk2OLZH1r4Uj+g23H3o7kWjgArC2snpqBAwc2/jsrK0vFxcURC6gl33//vcrKyhpvl5eXq6SkRB07dlT37t1jEgNgBdGeJm0nviSvyWKEJmwNAcBcttrQcu3atbrooosab/vqZcaNG6cFCxaYFBVgjmhNk7YjkjwAUphJTV5enlyuwCeLf/7zn2EH1JIhQ4YojMlaAOIASR6AsJKa3/72t363PR6PNmzYoOLiYt1xxx2RiAsAAKBVwkpqAi2w98QTT2jt2uDTKwGwASMARFpEa2pGjBihadOmaf78+ZF8WsBxmttlO4fCVgBok4ju0r148WJ17Bi5Ba4AJwq4y/a/N2C0087S9V5Dq7ZVaUnJt1q1rUr1XmreAJgn7MX3ji0UNgxDlZWV2rNnj/785z9HLDjYH0Ms/py0ASO9TQCsJqyk5vj9lhISEtSlSxcNGTJEffr0iURccAAuek0F22XbtwHjmvJ9lp7J4+ttOj458/U2+RYBBIBYCiupmTFjRqTjgMNw0Wve7gOBE5pwjjODk3qbADhLyElNqJtJSbLUZpGIPS56gWWlpwQ/qBXHmcEpvU2wLoatEa6Qk5oOHTq0uODeserr68MOCPbHRS+wQXkdQ9qAcVAEd5SONCf0NsG6GLZGW4Sc1Lz//vuN/96+fbvuvvtujR8/XoWFhZKkVatW6dlnn9Xs2bMjHyVshYteYL4NGFvaZdvqGzA6obcJR1mpV4Rha7RVyEnNT37yk8Z//+53v9Of/vQnXXvttY33XX755TrzzDP1l7/8RePGjYtslLAVLnots/sGjE7obUIDK/WKMGyNSAirUHjVqlV68sknm9w/cOBATZgwoc1Bwd646AVn5w0YndDbBOv1ijBsjUgIa/G93NxczZs3r8n9Tz/9tHJzc9scFOzNd9GTjl7kfLjoHeXbgPGK/iersFcnW7WHr7cpO9O/ty07M4UhAhsI1isiNfSKxHIxRYatEQlh9dTMmTNHP//5z/XWW29p8ODBkqQ1a9boq6++0ssvvxzRAGFPdh9iQXB27m2Kd1bsFWHYGpEQVlIzcuRIbd26VXPnztUXX3whSRo1apRuvvlmemrQiIue8/l6m2AvVuwVYdgakRD2hpa5ubl68MEHIxkLHIiLHmA9VuwVoVYLkRByUvPZZ5+poKBACQkJ+uyzz1o8tm/fvm0ODABCYaUpyXZh1V4Rhq3RViEnNf3791dlZaWysrLUv39/uVwuGUbTXweXy8XiewBiwkpTku3Eyr0iDFujLUJOasrLy9WlS5fGfwOAmaw2JdlurNwrwrA1whVyUtOjR49m/w0AscZCbZFBrwicJqx1ap599lm9+eabjbfvvPNOdejQQeedd56+/vrriAUHAM1pzZRktMzO6yUBxwsrqXnwwQeVmpoqqWF14ccff1y///3v1blzZ02ZMiWiAQLA8aw4JRmA+cKa0v3NN9/olFNOkSS99tpr+sUvfqGbbrpJ559/voYMGRLJ+ACgCStOSQZgvrB6atq3b6+qqipJ0jvvvKOioiJJUkpKin744YfIRQcAzfBNSQ40UOJSwywoFmoD4ktYSU1RUZEmTJigCRMmaOvWrRo5cqQkadOmTerZs2ck4wOAJthfLHz1XkOrtlVpScm3WrWtKqb7OwHRFtbw0xNPPKF7771X33zzjV5++WV16tQw9W7dunW69tprIxogADTHylOSrSqcdX1Y3NB+4vkzCyup6dChgx5//PEm98+aNavNAQFAqJiSHLpw1vVhcUP7iffPLKzhJ0n68MMPNXbsWJ133nn69ttvJUnPP/+8Pvroo4gFBwDBMCU5uGDr+kgN6/ocOxTlS4KOnzrvS4KKSyuiFzDCwmcWZlLz8ssva9iwYUpNTdX69etVV1cnSaqurmaTSwCwmNau6xNOEgRz8Zk1CCupuf/++/Xkk09q3rx5SkpKarz//PPP1/r16yMWHAAgPMcWBH9ctjekn/Gt68PihvbDZ9YgrJqaL7/8UhdeeGGT+zMzM7V///62xgQAaIPm6ipC4VvXh8UN7SfUz+Ktfw9BObX2LKyemuzsbJWVlTW5/6OPPtKPfvSjNgcFAAhPoLqKlhy/rg+LG9pPqJ/Fc6u+1rXzVuuCh99zZI1NWEnNjTfeqNtuu02ffvqpXC6Xdu7cqb/97W+6/fbbNXHixEjHCAAIQUt1FYE0t64PixvaT7DP7HhOLR4OK6m5++67NWbMGF1yySX6/vvvdeGFF2rChAmaOHGiJkyYEOkYgbjBwmhoi2B1Fc3JzkxpMp2bxQ3tp6XPrDlOLR4Oq6bG5XLpnnvu0R133KGysjJ9//33ys/P11NPPaW8vDxVVlZGOs64UO81tHZbFettxKlYri9R7zX0afk+rdvrUqfyfSo8JYvvmgOEWlcx+aJe6t01vcXzDIsb2k+gzyyQY4uHC3t1in6AMdCqpKaurk4zZ87U8uXL5Xa7dccdd2j06NGaP3++fvaznykxMZFdusO0scql2Y+sVGVNXeN98bRgUrwLZ2G0trzW0ZNeop77ai3fNYcIta7i/FO6hHQRY3FD+zn2M3urtELPrfo66M84qeC7VcNP06dP19y5c9WzZ0+Vl5frqquu0k033aQ5c+bokUceUXl5ue66665oxepYb2/apWe2JvglNJJzxzydIlJDRbFcX4LFuZwtGrUwLG5oP77PbESIf6Q4qeC7VT01L730kp577jldfvnlKi0tVd++fXXkyBFt3LhRLhdf9HDUew3dv+yLZh8z1HASmrV0s4ryszmZRFhb9keJ5FBRa9aXaEsXcbDkie+a/fnqKiYuXC+X5PdZUwsTf3xJbmX1oWZ/711qGE50UsF3q3pq/vWvf2nAgAGSpIKCArndbk2ZMoWEpg3WlO/7dw9N820YLwsmxVpxaYUuePg9XTtvtW57oaRVUxwj3dsRqzVBWJwrPvjqKrIz/f/6bq4gGM4WjwXfreqpqa+vV3Jy8tEfPuEEtW/fPuJBxRMWuYq9ttSvRKO3I1ZrgvBdix/UwsAn3gq+W5XUGIah8ePHy+12S5IOHTqkm2++We3atfM77pVXXolchA7HIlex1dakJBpDRbHqIua7Fl98dRVAPCW5rUpqxo0b53d77NixEQ0mHg3K66jsDLcqaw6puSEoJ455mqmtSUk0ejtiVQcRj+PrABrES5LbqqRm/vz50YojbiUmuHTvyD6a/EIJhX0x0NakJFq9HbHoIqaIFIDThbX4HiJr2Blddf2pXi2rTPOb1u3UMU8ztTUpiWZvRyy6iONtfB1AfCGpsYh+nQzd+csLteFfBxw/5mmmtiYl0e7tiEUXsS95WlW2W+98+KmG/ngwKwoDcISw9n5CdLDIVfRFYoqjE6bMJia4NDivowZ0NjSY5NmS2AcMaD16ahB3IjEEE0+zCRB7sdwHDHASkhrEpUgkJfEymwCxFct9wOJVW1YTh7WR1CBukZTAatjKIvroBXM2amoAi6O2In6wlUV0BdripKL6kG5mQ1dHoKcGsDD+qowvbGURPS31gvnc/crn9ILZHD01gEVFeuNMWB9bWURPsF4wSdp/0KPH3yuLUUSIBpIawIKC1VZIDbUVDEU5i28dpUD9BC419NSxlUXrhdq7Nf+Tcn6vbIykBrAgaiviUyTWUULzQu3d2n/Qw++VjZHUABZEbUX8csLijlY0KK+jOqQmhXQsv1f2RaEwYEHUVsQ3FneMvMQEl647v6fm/P2roMfye2VfJDWABUVz40zYA+soRd7ki3tr/ifbtf+gp9nH+b2yP4afAAuitgJW4aR1khITXHroyjObLcTm98oZ6KlBi1hO3DyR2KPKjvjOWYcT10lq6ffq//00X5mpyVpS8i3fPZsiqUFATjyh2U281VZE4ztHkhQeJ+9B1dzv1Xe1dbrvTc53dkdSg2Y5+YRmN/FSWxGN7xyJeXjiYQ+qY3+viksrNGnRBs53DkBNDZpg4TfEWjS+c6zIHL54WieJ852zkNSgiXg6ocEaIv2d40LVNvG0ThLnO2chqUET8XRCgzVE+jvHhapt4mmdJM53zmK7pOaJJ55Qz549lZKSosGDB2vNmjVmh+Q48XRCgzWE+l3q3M4d0nFcqNomnvag4nznLLZKal588UVNnTpVM2bM0Pr169WvXz8NGzZMu3fvNjs0R4mnExqsIdh3zuf2lzaGVAvDhapt4mmdJM53zmKr2U9/+tOfdOONN+q6666TJD355JN688039cwzz+juu+9ucnxdXZ3q6uoab9fU1EiSPB6PPJ7mV5Q0gy8WK8V0z4jTdMsLG+WS/OoSXMc87q0/Im+9CcEdw4ptZxdWa7tA37lj7appKPJ97Jp+GnZG14DPdVa3dGVnuLWrpq6FFZndOqtbeljv32ptFw2XnNZZj13TT/cv+0KVNUfPo9mZbt0zoo8uOa2zY9rODuc7K7ZbLIX6vl2GYdiiUu7w4cNKS0vT4sWLNXr06Mb7x40bp/3792vJkiVNfmbmzJmaNWtWk/sXLVqktLS0aIbrCBurXHple4L2Hz76N0yHZENX9vSqXydbfG1gAq8hbatxqcYjZSRJvTIMhfoH/cYql14uT1C1p6UfMNQhWZpxdn2Lz7uxyqVntvo6o489sOG7e/2pfI9D0ZbP004431nbwYMHNWbMGFVXVysjIyPgcbZJanbu3KmTTz5Zn3zyiQoLCxvvv/POO7VixQp9+umnTX6muZ6a3Nxc7d27t8VGiTWPx6Ply5erqKhISUmh7SIbK/VeQ2u//k67D9QpK92tgT1OtFSXs5Xbzuqi0XZvb9rV9C/7DLfuHdmnxZ6VY32yrUrjFqwLetzC6wdqcJAhgebiyfl3T0Oo8TSH7134rNx2Vj7fWbndYqGmpkadO3cOmtTYaviptdxut9zupoWFSUlJlvxSWDGuJEkXnBr+yT9WrNh2dhGptisurdAtL2xsMtyzq6ZOt7ywMeQFzPYfCq2Pv+rgkaBxX9a/m0b0PTlqKwrzvQufFdvODuc7K7abFP2Vu0N9z7ZJajp37qzExETt2rXL7/5du3YpOzvbpKgASJFdgTbSRb7xsiKz3dR7Da3dVsX2FQ5gpZW7bTP7KTk5WQMGDNC7777beJ/X69W7777rNxwFIPYiuS4Ms1Gcb2OVS0MeWalr563WbS+U6Np5q3XOA3/Xss9Y5dlurLZyt22SGkmaOnWq5s2bp2effVZbtmzRxIkTVVtb2zgbCoA5IrkuTDxNJ45Hb2/apWe2JvjVOUnSvtrD+s9F6/XAm5tNigytZcWVu20z/CRJV199tfbs2aPp06ersrJS/fv3V3Fxsbp2tfYYKGAVx457d0o7QZE610R6yGh4QY7mjj27SZd2tgM3o4ynXcTrvYbuX/ZFi8fM+7BckqF7fnpGbIJC2FrTQxurIWBbJTWSNHnyZE2ePNnsMADbaW7cu0NyopJ67tJl/bu16bl9Q0aV1YdaWBemdUNGwwtyVJSf7egLvpVqEWJhTfm+f/fQtPwZzvtwu87KPVEj+54Um8AQFiuu3G2r4ScA4Qk07r3/sHTLC6Gt0tuSaA0Z+Yp8r+h/sgp7dXJcQmOlWoRYaM3F7d4lpWw4anFWXLmbpAZwuJbGvX0pRyTGvX1DRtmZ/iew7MyUkKdzxwsr1iLEQmsubvtqPWw4anFWLOq33fATECnxUssQy3HveBgyigQr1iLEwqC8juqYlqR9B0Nb8p4NR63N10M7ceH6gFtMxLqon6QGcSmeahliPe7NujDBxfIzsVLynpjg0sxRp+vWFzcqWF2NxIajdmC1on6SGsQdXy3D8R37vloGpw2VWHHcO97F6jOxYvI+oiBbF63coPcrEls8jrWI7MNKPbTU1CCuxGMtgxXHveNdLD4TKxcij+5p6Przugd83CXWIrIbqxT1k9QgrkRy5Vu7aGlmki+Vs/sFpN5raNW2Ki0p+VartlVZPimN9gKDdkjep43ooz+POUsd2/nv6ZNDYTnagOEnxBUrrqsQC8MLcvTEmLN075JS7as9WqTZIVm6/8p+tr6AWHGIJRTRrEWwSyHyyL4naVhBjiWGLeAMJDWIK06tLwlWDFpcWqH73tzil9B0TEvSFd0OadgZ9l2R2+71UdGqRbBT8k5hOSKJpAZxJRor35otWE9FoAv/dwc9mr81QQM2tX1FYTNEcmdwM0Xjou7U5B0IhpoaxBWnbZYYrBh02WcVQWsrHnjrC8vXoDQnHuujQtXWQmS71SgBPvTUIO5YbV2FcIXSU/H/lpSqqvZwC8/iUkV1nem1FeGw0xBLrLVlUTS71igBEkkN4pSV1lUIVyg9FS0nNEfZ8cLPEEvLwkneY1Wj5DWkT8v3qergEVv+7sG6SGoQt+xeoBjJRMSOF34n1kdFWmuS91jVKL29aZdmrU/U/tVrG++jJwiRQk0NYFOhJiId2yW1sCC9oZxMty0v/E6rj4qWUBdFi0WNUnFphW55YaP2H9eBaIUFAeEMJDVAjEWqCDPUYtD7ryhovH3845J0z4g+tr3wszN45ES7Rsm/J8j/+2aVBQFhfww/ATEUySLMUItBhxfkaG6Cq5naCrdGdD1o63VqJGfUR1lBtGuU7LIgIOyNpAaIkWgUYYZaDNrchf+sbul6u/ittr4tS7B7fZQVRLtGidlqiAWSGiAGolmEGWpPxfEXfo/Hc/xTIY61ZRp4KJithligpgaIgWgXYVplh1zYWzRrlNgtHrFATw0QA3S9x69g+3JZTbRqlI7tCTraP9mA2WqIFJIaBGS3k3EstbZt6HqPT3ZdnTdaNUrDC3L02DX9dO8rJX7Tuu22mjesi6QGzbLryTgWwmkbFoqLP3bfQTxahp3RVZ7t9eqSfy4rCiPiqKlBE8E2SYznBbLCbRsWiosvwQrDpfhekyXBJQ3O60gNGCKOpAZ+OBkH1ta2YaG40Dhhh2h2EAfMwfAT/LBAVmCRaBsWimuZGcOe0agdozAcMAdJDfxwMg4sUm3DQnHNM6MGpbkkqmO7ZN1/RYFG9g3/taxaGE7xP5yOpAZ+rHoytgLaJnpitUP0sQIlUftqD+s/F63Xb/6Vp2kj88N6bisWhlP8j3hATQ38sEBWYLRN9MS6BqWlJMrnqZXleqPk27Ce32qF4RT/I16Q1MCP1U7GVkLbRE+shz2DJVE+t7xYomWfhXfBt0phOMX/iCcMP6GJUDdJjEe0TXTEemgv1OTIMKT/XLReTyaEl4RYoTCc4n/EE5IaNMsKJ2Orom0iL9Y1KK1Njv771c91cZ+uSj6h9Z3bZheGU/yPeEJSg4DMPhlbGW0TWdHeIfp4g/I6qmO7ZO2rPRz8YEn7aj06d/a7evBnBbbrjdu+92BIx1HgDiegpgaAJcSyBiUxwaX7ryho1c/sqz1su6La4tIKPfr3rUGP69guSZU1h2y72CHgQ08NAMto7dBeW9ZdGdk3R7/5V56eWlneqhhnLd2sIb1/3KqfMUMoM7x89tV6NOXFEklM84a9kdQgaljoy5mi/bmGOrQXiYXzpo3M15knZeqWF0tkhHD19xXVPrf6a3WxeIdGqDO8jhfvG27C3khqEBUs9OVMVvlcI7lw3mX9T1ZCQoL+c9H6kF//wbe2qkNyopJ67tJl/bu1IvLYCbfwN1qLHQKxQE0NIo6FvpzJKp9rqAvnLftsZ8jPObJvjp4ce7Y6tksK+Wf2H5ZueWGjZb/PbSn8ZcNN2BVJDSIq2EJfhhqmx766wb47MMcjKy3gFuqwyr1LSlsVz/CCHK2edqk6tksO8ScaejCsunBdsBWwQ8E078hzwi70VsbwEyIqlAsORYn2Y6UF3EK90O6r9bQ6nuQTEvTgzwo0cWHDUFSwy00s3ne4NUwtTZMPFdO8I8sqw7dORk8NIqq1f9kxJGUPVlrArTUX2nDiCTS1PNKvE4ri0gpd8PB7unbeat32QomunbdaFzz8Xsi/L4HeS05mijqkJbGPWQxZZfjW6eipQUS19i87ihLtoTXbGER7dlRrFs4Lt6fBN7V8wcfluu/NLVF7nZYEKoZu7eykQNPkl2+ujNlih/HOjF3o4xU9NYiocMbxKUq0vlB3KP+utq5NPQuhCHXhvLb2NCQmuDT+/DxTdmaPdA2Tb5r8Ff1PVmGvTkpMcFlmw814EOtd6OMZSQ0iqqWdrIOhKNG6Qtmh/PJ+OZq0aENMutdH9s3Rby7MC/i4S/49DeEWZ7b8fW54jmj0aMTqIji8IEcf3XWx/u/Gc/U/1/TX/914rj6662ISmgiz0vCt05HUIOLCqUmQKEq0upb+sn9izNl6fWNFTGdHTRuZrz+POavJNOyc43oaolWX0iFZeuyaflFJAGJ5EWyuFweRFetd6OMZNTWIimPH8Surf9B9b27Rd7WHY7IDM6InUH2GWbOjRvY9ScMKcgLW8ESrLqVT2gnas3m1hp3RNWLv5VhcBJ0l1rvQxzOSGkTNscvdpyYnUpToEM1tY2Bm93qgbRUiXZx57Ot4PB4tC14/HDY7XwTZHqWpWO9CH88YfkJMUJTobFbsWbBzcWYoNUxWvAi2dajPyTgHxgY9NYiZ1u7ADPuwYs+C3YszfRfB4xdry7boYm2RGupzMs6B0UdSg5gKdQdm2IsVu9et2HvUWna5CLZmqC/ecQ6MLoafAEREtLvXWzstO9S1daxYl3IsO8xOsvNQH5yFnhoAEROtnoVw9syxYu+RU7VuqC8jusEgrtFTAyCiIt2z0JY9cyjOjA0nDPXBGeipAWBZkZiWbZe6FDtrTaG4t/5IrMNDHKGnBoBlRapWww51KXZm1ynocB6SGgCWZfdp2fGEoT5YAcNPACyLWg17YagPZiOpAWBZVlzUDy1jHRaYieEnAJZFrQaA1iCpAWBp1GoACBXDTwAsj1oNAKEgqQFgC/Faq1HvNUjmgBCR1LRRtE44TjqROem9WJXZbWz26ztVc9tDdEhN0nXn52nyxafQxsBxbJPUPPDAA3rzzTdVUlKi5ORk7d+/3+yQwtqPxsznNYOT3otVmd3GZr++U/m2hzh+1tf+Hzya8/etmv9JuR668kzaGDiGbQqFDx8+rKuuukoTJ040OxRJbduPxoznNYOT3otVmd3GZr++U7W0PYTP/oMe3UwbA35sk9TMmjVLU6ZM0Zlnnml2KEH3o5Ea9qOp97Z0Sord85rBSe/Fqsxu42i8fr3X0KptVVpS8q1WbauK2+9HsO0hjsXvEXCUbYafwlFXV6e6urrG2zU1NZIkj8cjj8cT9vN+GuJ+NKvKdmtwCIuC+WJZvW1PRJ/XTJFuo0B8bdeWz9Ou2trGbW27SH/Gb2/apfuXfaHKmqO/s9kZbt07so+GndE1rBijJdrfu4r9taEfa5Nzgk+k2q7ea2jt199p94E6ZaW7NbDHiY6uMYrnc50U+vt2dFIze/ZszZo1q8n977zzjtLS0sJ+3nV7XZISgx73zoefqmpL6H9BvbdqXVSe1wzRaqNAli9f3ubnsJtItXG4bRfJz3hjlUvPbPV1HB+9MFXWHNLkF0p0/ale9etkve98tL53/6wOrW197HBOOF5b2m5jlUuvbE/Q/sNHvysdkg1d2dOa35NIisdznSQdPHgwpONMTWruvvtuPfzwwy0es2XLFvXp0yes5582bZqmTp3aeLumpka5ubkaOnSoMjIywnpOSepUvk/PfbU26HFDfzw45J6a5cuX6+LCAXruq5KIPa+ZIt1GgfjarqioSElJSWE/jx21tY3b2naR+ozrvYZmP7JSUl0zj7rkkvTWrjTd+csLLfOXeLS/d/VeQ4sfWenXa9USO5wTfNradm9v2qX5qzY2GfasPuzS/K2Jeuyafpbr2YuEeD7XSUdHWoIxNam5/fbbNX78+BaP+dGPfhT287vdbrnd7ib3JyUltelLUXhKVkj70RSektWqk/C5vbpE5XnNEK02CqStn6kdRaqNw227SL3+2m1VLV68G4ax6rThXwcst05NtL53SZJmXn5Gs7OfjmWnc8Lxwmm7eq+hB976MmAdl0vSA299qRF9T7Zde4QqHs91kkJ+z6YWCnfp0kV9+vRp8b/k5GQzQ2xWtPajcdI+N056L1ZldhtH6vV3HwitIDbU45zCtz1Eh7TmT+bx+HsUrIDaV8e1pnxf7IKCpdhm9tOOHTtUUlKiHTt2qL6+XiUlJSopKdH3339vSjzR2o/GSfvcOOm9WJXZbRyJ189KTwl6TGuOc5LhBTlad2+RplzaWx1S/ZObePw9IgFGMLYpFJ4+fbqeffbZxttnnXWWJOn999/XkCFDTIkpWvvROGmfGye9F6syu43b+vqD8jqGNIw1yCY1I5GWmODSbZeeqskX94773yMSYARjm6RmwYIFWrBggdlhNBGt/WictM+Nk96LVZndxm15fd8w1sSF6+WS/BKbeBxiCcTsz9gKSIARjG2GnwA4l9nDaLAHs+vIYH226akB4GxmD6NFAxt9Rp4vAT5+v7Fs9huDSGoAWIiThljY6DN6nJgAIzJIagAgwgLtsO3b6JMhtbZzUgKMyKGmBgAiyOyNRoF4RlIDABHEAnGAeUhqACCCWCAOMA9JDQBEEAvEAeahUBgAIogF4sBUfvOQ1ABABLFCctvZOSlgKr+5SGoAIMJYIC58dk4KmMpvPpIaAIgCFohrPTsnBcGm8rvUMJW/KD+b70AUkdQAQJSwQFzo7J4UtGYqP9+J6GH2EwDAdHZf34ep/NZAUgMAMJ3dkwKm8lsDSQ0AwHR2Twp8U/kDDYy51FDwzFT+6CKpAQCYzu5JgW8qv6Qm74Gp/LFDUgMAMJ0TkgLfVP7sTP/epOzMFEvP3HISZj8BACzBCev7MJXfXCQ1AADLcEJSwFR+85DUAAAshaQA4aKmBgAAOAJJDQAAcASSGgAA4AgkNQAAwBEoFAYAtEq917D17CQ4F0kNACBkxaUVTdaRybHROjJwNoafAAAhKS6t0MSF65vspl1ZfUgTF65XcWmFSZEBDUhqAABB1XsNzVq6WUYzj/num7V0s+q9zR0BxAZJDQAgqDXl+5r00BzLkFRRfUhryvfFLijgOCQ1AICgdh8InNCEcxwQDSQ1AICgstJTgh/UiuOAaCCpAQAENSivo3IyUxRo4rZLDbOgBuV1jGVYgB+SGgBAUIkJLs0YlS9JTRIb3+0Zo/JZrwamIqkBAIRkeEGO5o49W9mZ/kNM2Zkpmjv2bNapgelYfA8AELLhBTkqys9mRWFYEkkNAKBVEhNcKuzVyewwgCYYfgIAAI5AUgMAAByBpAYAADgCSQ0AAHAEkhoAAOAIJDUAAMARSGoAAIAjkNQAAABHIKkBAACOEFcrChuGIUmqqakxORJ/Ho9HBw8eVE1NjZKSkswOx1Zou/DRduGj7cJH24Un3tvNd932XccDiauk5sCBA5Kk3NxckyMBAACtdeDAAWVmZgZ83GUES3scxOv1aufOnUpPT5fLZZ3N12pqapSbm6tvvvlGGRkZZodjK7Rd+Gi78NF24aPtwhPv7WYYhg4cOKCTTjpJCQmBK2fiqqcmISFB3bp1MzuMgDIyMuLyyxoJtF34aLvw0Xbho+3CE8/t1lIPjQ+FwgAAwBFIagAAgCOQ1FiA2+3WjBkz5Ha7zQ7Fdmi78NF24aPtwkfbhYd2C01cFQoDAADnoqcGAAA4AkkNAABwBJIaAADgCCQ1AADAEUhqLGbr1q264oor1LlzZ2VkZOiCCy7Q+++/b3ZYtvHmm29q8ODBSk1N1YknnqjRo0ebHZKt1NXVqX///nK5XCopKTE7HMvbvn27brjhBuXl5Sk1NVW9evXSjBkzdPjwYbNDs6QnnnhCPXv2VEpKigYPHqw1a9aYHZLlzZ49W+ecc47S09OVlZWl0aNH68svvzQ7LMsiqbGYyy67TEeOHNF7772ndevWqV+/frrssstUWVlpdmiW9/LLL+tXv/qVrrvuOm3cuFEff/yxxowZY3ZYtnLnnXfqpJNOMjsM2/jiiy/k9Xr11FNPadOmTZozZ46efPJJ/fd//7fZoVnOiy++qKlTp2rGjBlav369+vXrp2HDhmn37t1mh2ZpK1as0KRJk7R69WotX75cHo9HQ4cOVW1trdmhWZMBy9izZ48hyVi5cmXjfTU1NYYkY/ny5SZGZn0ej8c4+eSTjaefftrsUGxr2bJlRp8+fYxNmzYZkowNGzaYHZIt/f73vzfy8vLMDsNyBg0aZEyaNKnxdn19vXHSSScZs2fPNjEq+9m9e7chyVixYoXZoVgSPTUW0qlTJ5122ml67rnnVFtbqyNHjuipp55SVlaWBgwYYHZ4lrZ+/Xp9++23SkhI0FlnnaWcnByNGDFCpaWlZodmC7t27dKNN96o559/XmlpaWaHY2vV1dXq2LGj2WFYyuHDh7Vu3TpdeumljfclJCTo0ksv1apVq0yMzH6qq6slie9YACQ1FuJyufT3v/9dGzZsUHp6ulJSUvSnP/1JxcXFOvHEE80Oz9L++c9/SpJmzpype++9V2+88YZOPPFEDRkyRPv27TM5OmszDEPjx4/XzTffrIEDB5odjq2VlZXpscce029+8xuzQ7GUvXv3qr6+Xl27dvW7v2vXrgytt4LX69Vvf/tbnX/++SooKDA7HEsiqYmBu+++Wy6Xq8X/vvjiCxmGoUmTJikrK0sffvih1qxZo9GjR2vUqFGqqKgw+22YItS283q9kqR77rlHP//5zzVgwADNnz9fLpdLL730ksnvwhyhtt1jjz2mAwcOaNq0aWaHbBmhtt2xvv32Ww0fPlxXXXWVbrzxRpMih5NNmjRJpaWleuGFF8wOxbLYJiEG9uzZo6qqqhaP+dGPfqQPP/xQQ4cO1Xfffee3tXzv3r11ww036O677452qJYTatt9/PHHuvjii/Xhhx/qggsuaHxs8ODBuvTSS/XAAw9EO1TLCbXt/uM//kNLly6Vy+VqvL++vl6JiYn65S9/qWeffTbaoVpOqG2XnJwsSdq5c6eGDBmic889VwsWLFBCAn8vHuvw4cNKS0vT4sWL/WYkjhs3Tvv379eSJUvMC84mJk+erCVLlmjlypXKy8szOxzLOsHsAOJBly5d1KVLl6DHHTx4UJKanBATEhIaeyLiTahtN2DAALndbn355ZeNSY3H49H27dvVo0ePaIdpSaG23f/+7//q/vvvb7y9c+dODRs2TC+++KIGDx4czRAtK9S2kxp6aC666KLG3kESmqaSk5M1YMAAvfvuu41Jjdfr1bvvvqvJkyebG5zFGYahW265Ra+++qo++OADEpogSGospLCwUCeeeKLGjRun6dOnKzU1VfPmzVN5ebl++tOfmh2epWVkZOjmm2/WjBkzlJubqx49eugPf/iDJOmqq64yOTpr6969u9/t9u3bS5J69eqlbt26mRGSbXz77bcaMmSIevTooT/+8Y/as2dP42PZ2dkmRmY9U6dO1bhx4zRw4EANGjRIjz76qGpra3XdddeZHZqlTZo0SYsWLdKSJUuUnp7eWIOUmZmp1NRUk6OzHpIaC+ncubOKi4t1zz336OKLL5bH49EZZ5yhJUuWqF+/fmaHZ3l/+MMfdMIJJ+hXv/qVfvjhBw0ePFjvvfceRdaImuXLl6usrExlZWVNEkBG9v1dffXV2rNnj6ZPn67Kykr1799fxcXFTYqH4W/u3LmSpCFDhvjdP3/+fI0fPz72AVkcNTUAAMARGPwFAACOQFIDAAAcgaQGAAA4AkkNAABwBJIaAADgCCQ1AADAEUhqAACAI5DUAAAARyCpAeDH5XLptddeMzsMS/jggw/kcrm0f/9+SdKCBQvUoUMHU2MCEBhJDRBnxo8f77dT8vEqKio0YsSI2AVkI1dffbW2bt0a8vE9e/bUo48+Gr2AAPhh7ycAfqywEaNhGKqvr9cJJ7T9FBXJ50pNTTVlE8HDhw8rOTk55q8L2A09NQD8HDv8tH37drlcLr3yyiu66KKLlJaWpn79+mnVqlV+P/PRRx/pxz/+sVJTU5Wbm6tbb71VtbW1jY8///zzGjhwoNLT05Wdna0xY8Zo9+7djY/7hnneeustDRgwQG63Wx999FGT2HzxvPDCCzrvvPOUkpKigoICrVixIuhzeb1ezZ49W3l5eUpNTVW/fv20ePFiv+dftmyZTj31VKWmpuqiiy7S9u3b/R5vbvhp6dKlOuecc5SSkqLOnTvrZz/7maSGDQi//vprTZkyRS6XSy6Xq/FnXn75ZZ1xxhlyu93q2bOnHnnkEb/n7Nmzp+677z79+te/VkZGhm666aYAnxYAPwaAuDJu3DjjiiuuCPi4JOPVV181DMMwysvLDUlGnz59jDfeeMP48ssvjV/84hdGjx49DI/HYxiGYZSVlRnt2rUz5syZY2zdutX4+OOPjbPOOssYP35843P+9a9/NZYtW2Zs27bNWLVqlVFYWGiMGDGi8fH333/fkGT07dvXeOedd4yysjKjqqqqSWy+eLp162YsXrzY2Lx5szFhwgQjPT3d2Lt3b4vPdf/99xt9+vQxiouLjW3bthnz58833G638cEHHxiGYRg7duww3G63MXXqVOOLL74wFi5caHTt2tWQZHz33XeGYRjG/PnzjczMzMZ43njjDSMxMdGYPn26sXnzZqOkpMR48MEHDcMwjKqqKqNbt27G7373O6OiosKoqKgwDMMw1q5dayQkJBi/+93vjC+//NKYP3++kZqaasyfP7/xeXv06GFkZGQYf/zjH42ysjKjrKwstA8XiHMkNUCcCSepefrppxsf37RpkyHJ2LJli2EYhnHDDTcYN910k99zfPjhh0ZCQoLxww8/NPsa//jHPwxJxoEDBwzDOJqIvPbaay3G7ovnoYcearzP4/EY3bp1Mx5++OGAz3Xo0CEjLS3N+OSTT/ye74YbbjCuvfZawzAMY9q0aUZ+fr7f43fddVeLSU1hYaHxy1/+MmC8PXr0MObMmeN335gxY4yioiK/++644w6/1+7Ro4cxevTogM8LoHkMPwEIqm/fvo3/zsnJkaTG4aONGzdqwYIFat++feN/w4YNk9frVXl5uSRp3bp1GjVqlLp376709HT95Cc/kSTt2LHD73UGDhwYUjyFhYWN/z7hhBM0cOBAbdmyJeBzlZWV6eDBgyoqKvKL87nnntO2bdskSVu2bNHgwYMDvk5zSkpKdMkll4QUs8+WLVt0/vnn+913/vnn66uvvlJ9fX2z8QMIDYXCAIJKSkpq/LevNsTr9UqSvv/+e/3mN7/Rrbfe2uTnunfvrtraWg0bNkzDhg3T3/72N3Xp0kU7duzQsGHDdPjwYb/j27VrF7GYj32u77//XpL05ptv6uSTT/Y7zu12h/0a0SwajmRbAPGCpAZAm5x99tnavHmzTjnllGYf//zzz1VVVaWHHnpIubm5kqS1a9e26TVXr16tCy+8UJJ05MgRrVu3TpMnTw54fH5+vtxut3bs2NHYS3S8008/Xa+//nqT12lJ37599e677+q6665r9vHk5GS/3hff63z88cd+93388cc69dRTlZiY2OLrAWgZSQ0Qh6qrq1VSUuJ3X6dOnRqTjta46667dO6552ry5MmaMGGC2rVrp82bN2v58uV6/PHH1b17dyUnJ+uxxx7TzTffrNLSUt13331tiv+JJ55Q7969dfrpp2vOnDn67rvvdP311wc8Pj09Xf/1X/+lKVOmyOv16oILLlB1dbU+/vhjZWRkaNy4cbr55pv1yCOP6I477tCECRO0bt06LViwoMU4ZsyYoUsuuUS9evXSNddcoyNHjmjZsmW66667JDXMYlq5cqWuueYaud1ude7cWbfffrvOOecc3Xfffbr66qu1atUqPf744/rzn//cpjYBIGY/AfFm3LhxhqQm/91www2GYTRfKLxhw4bGn//uu+8MScb777/feN+aNWuMoqIio3379ka7du2Mvn37Gg888EDj44sWLTJ69uxpuN1uo7Cw0Hj99df9ntdX3OsryA3EF8+iRYuMQYMGGcnJyUZ+fr7x3nvvNR4T6Lm8Xq/x6KOPGqeddpqRlJRkdOnSxRg2bJixYsWKxmOWLl1qnHLKKYbb7TZ+/OMfG88880yLhcKGYRgvv/yy0b9/fyM5Odno3LmzceWVVzY+tmrVKqNv376G2+02jj3dLl682MjPzzeSkpKM7t27G3/4wx/8nrO5AmMAwbkMwzDMSKYAoLW2b9+uvLw8bdiwQf379zc7HAAWw+wnAADgCCQ1AADAERh+AgAAjkBPDQAAcASSGgAA4AgkNQAAwBFIagAAgCOQ1AAAAEcgqQEAAI5AUgMAAByBpAYAADjC/wexGGbJPz+L3wAAAABJRU5ErkJggg==", "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.19" } }, "nbformat": 4, "nbformat_minor": 4 }