{
 "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": "2021-02-02T06:54:49.261749Z",
     "iopub.status.busy": "2021-02-02T06:54:49.260792Z",
     "iopub.status.idle": "2021-02-02T06:54:50.675737Z",
     "shell.execute_reply": "2021-02-02T06:54:50.676432Z"
    }
   },
   "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": "2021-02-02T06:54:50.684663Z",
     "iopub.status.busy": "2021-02-02T06:54:50.683808Z",
     "iopub.status.idle": "2021-02-02T06:54:50.687970Z",
     "shell.execute_reply": "2021-02-02T06:54:50.689001Z"
    }
   },
   "outputs": [],
   "source": [
    "raw = StringIO(\"\"\"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\"\"\")"
   ]
  },
  {
   "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": "2021-02-02T06:54:50.693554Z",
     "iopub.status.busy": "2021-02-02T06:54:50.692228Z",
     "iopub.status.idle": "2021-02-02T06:54:50.714471Z",
     "shell.execute_reply": "2021-02-02T06:54:50.715479Z"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(raw, header=None)\n",
    "df = df.melt()\n",
    "df[\"site\"] = 1 + np.floor(df.index / 10).astype(np.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": "2021-02-02T06:54:50.720221Z",
     "iopub.status.busy": "2021-02-02T06:54:50.718789Z",
     "iopub.status.idle": "2021-02-02T06:54:50.776744Z",
     "shell.execute_reply": "2021-02-02T06:54:50.777725Z"
    }
   },
   "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:                Tue, 02 Feb 2021   Deviance:                       6.1260\n",
      "Time:                        06:54:50   Pearson chi2:                     6.39\n",
      "No. Iterations:                    10                                         \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(\"blotch ~ 0 + C(variety) + C(site)\",\n",
    "          family=sm.families.Binomial(), data=df)\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": "2021-02-02T06:54:50.782390Z",
     "iopub.status.busy": "2021-02-02T06:54:50.780954Z",
     "iopub.status.idle": "2021-02-02T06:54:51.016812Z",
     "shell.execute_reply": "2021-02-02T06:54:51.017897Z"
    },
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Residual')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnA0lEQVR4nO3de5hddX3v8feXYWIGBhgu7UiG2KCF1CBiygAWTuvESwOoJEWpokWx2py04sFbNAGpWtpDbNSW81RLU8oB0ToKpEOUSFRg0xYMhRBCCGkgB2vIhFPlMuCQoZmEb//Yayd7dtZe+74ue39ez5Nn9rrsvX6/nZn1Xb+7uTsiIiLlHJR0AkREJN0UKEREJJIChYiIRFKgEBGRSAoUIiIS6eCkE9AKxxxzjM+aNavu97/44osceuihzUtQiihv2dXO+VPekrd+/fqn3f1Xwo61ZaCYNWsWDzzwQN3vz+VyDA0NNS9BKaK8ZVc75095S56Z/azcsUSrnszsOjP7uZk9Uua4mdn/MbNtZvawmf1m3GkUEel0SbdRXA+cHXH8HOCE4N8i4G9jSJOIiBRJNFC4+z8Dz0acsgD4huetA/rM7Nh4UiciIpD+NooB4Mmi7R3BvqdKTzSzReRLHfT395PL5eq+6Pj4eEPvTzPlLbvaOX/KW7qlPVBYyL7QyancfSWwEmBwcNAbaTzKSuNTPZS37Grn/Clv6Zb2QLEDmFm0fRywM6G0iKTSyIZRVqzdys6xCWb09bBk/mwWzh1IOlnSRpJuzK5kNfCBoPfTG4Hn3f2AaieRTjWyYZRlqzYxOjaBA6NjEyxbtYmRDaNJJ03aSNLdY78N/ASYbWY7zOzDZrbYzBYHp6wBngC2AX8P/ElCSRVJpRVrtzIxuXfKvonJvaxYuzWhFEk7SrTqyd0vrHDcgY/GlByRzNk5NlHTfpF6pL3qSUQizOjrqWm/SD0UKEQybMn82fR0d03Z19PdxZL5sxNKkbSjtPd6EpEIhd5N6vUkraRAIR2hnbuQLpw70DZ5kXRSoJC2V+hCWugdVOhCCugGK1IFtVFI21MXUpHGKFBI21MXUpHGKFBI21MXUpHGKFBI21MXUpHGqDFb2p66kIo0RoFCOoK6kIrUT1VPIiISSYFCREQiKVCIiEgkBQoREYmkQCEiIpEUKEREJFLSS6GebWZbzWybmS0NOX6EmX3PzDaa2WYz+1AS6RQR6WSJBQoz6wK+BpwDzAEuNLM5Jad9FHjU3U8BhoCvmNm0WBMqItLhkixRnA5sc/cn3H03MAwsKDnHgcPMzIBe4FlgT7zJFBHpbObuyVzY7N3A2e7+kWD7IuAMd7+k6JzDgNXAbwCHAe9x99vKfN4iYBFAf3//qcPDw3WnbXx8nN7e3rrfn2bKW3a1c/6Ut+TNmzdvvbsPhh1LcgoPC9lXGrXmAw8BbwZeA/zIzP7F3V844I3uK4GVAIODgz40NFR3wnK5HI28P82Ut+xq5/wpb+mWZNXTDmBm0fZxwM6Scz4ErPK8bcBPyZcuREQkJkkGivuBE8zs+KCB+r3kq5mKbQfeAmBm/cBs4IlYUyki0uESq3py9z1mdgmwFugCrnP3zWa2ODh+DXAlcL2ZbSJfVfVZd386qTSLiHSiRKcZd/c1wJqSfdcUvd4J/G7c6RIRkf00MltERCIpUIiISCQFChERiaSlUEU63MiGUa0nLpEUKEQ62MiGUZat2sTE5F4ARscmWLZqE4CCheyjqieRDrZi7dZ9QaJgYnIvK9ZuTShFkkYKFCIdbOfYRE37pTMpUIh0sBl9PTXtl86kQCHSwZbMn01Pd9eUfT3dXSyZPzuhFEkaqTFbpIMVGqzV60miKFCIdLiFcwcUGCSSqp5ERCSSAoWIiERSoBARkUgKFCIiEkmBQkREIilQiIhIpEQDhZmdbWZbzWybmS0tc86QmT1kZpvN7O640ygi0ukSG0dhZl3A14C3ATuA+81stbs/WnROH/B14Gx3325mv5pIYkVEOliSJYrTgW3u/oS77waGgQUl57wPWOXu2wHc/ecxp1FEpOMlOTJ7AHiyaHsHcEbJOScC3WaWAw4Drnb3b8STPBFpJi2QlF3m7slc2OwCYL67fyTYvgg43d0/VnTO3wCDwFuAHuAnwNvd/bGQz1sELALo7+8/dXh4uO60jY+P09vbW/f700x5y64s5+/enZNc/8hudr+8f9+0g+Di103jzBndmc5bJVnJ27x589a7+2DYsSRLFDuAmUXbxwE7Q8552t1fBF40s38GTgEOCBTuvhJYCTA4OOhDQ0N1JyyXy9HI+9NMecuuLOfv8uV3TgkSALtfhtu2d3HZ+4YynbdK2iFvSbZR3A+cYGbHm9k04L3A6pJzbgV+28wONrNDyFdNbYk5nSLSIC2QlG2JlSjcfY+ZXQKsBbqA69x9s5ktDo5f4+5bzOx24GHgZeBad38kqTSLSH1m9PUwGhIUtEBSNiQ6zbi7rwHWlOy7pmR7BbAiznSJSHMtmT+bZas2TVmfWwskZYfWoxCRltMCSdmmQCEisdACSdmluZ5ERCSSAoWIiERSoBARkUgKFCIiEkmBQkREIilQiIhIJAUKERGJpEAhIiKRNOBORBJ3785JLl9+Z2yjtrU2Rm0UKEQkUSMbRqesVTE6NsGyVZsAWnLzHtkwOmXeqVZfrx2o6klEErVi7dYD1qqYmNzLirVbW3a94skJW329dqBAISKJinutCq2NUTsFChFJVLk1KVq1VkXc12sHChQikqgl82czreRO1Mq1KpbMn01Pd1ds12sHasyW1FGPlM6ycO4Aj255lNu2d8Xyf661MWqXaKAws7OBq8kvhXqtuy8vc95pwDrgPe5+c4xJlJipR0pnOnNGN5e9byi262ltjNokVvVkZl3A14BzgDnAhWY2p8x5XyK/tra0OfVIEUmfJEsUpwPb3P0JADMbBhYAj5ac9zHgFuC0eJMnSVCPlPRSlWDnSrIxewB4smh7R7BvHzMbAH4PuCbGdEmC1CMlnQpVgqNjEzj7qwRHNowmnTSJQZIlCgvZ5yXbfw181t33moWdXvRhZouARQD9/f3kcrm6EzY+Pt7Q+9Ms7Xl7+6v2cv0LTBmANe2g/P5K6U573hqVZP6uzO1iYnLqn+fE5F6uvHUjfc8/3vDnt/P/XTvkLclAsQOYWbR9HLCz5JxBYDgIEscA55rZHncfKf0wd18JrAQYHBz0oaGhuhOWy+Vo5P1plva8DQFz6qziSHveGpVk/p69/bbw/S95U9LUzv937ZC3JAPF/cAJZnY8MAq8F3hf8QnufnzhtZldD3w/LEhIe1GPlHDNnDiv1vaGGX09jIa0E6lKsDMk1kbh7nuAS8j3ZtoCfNfdN5vZYjNbnFS6RNKoMHFeM9oI6mlv0CC1zpboOAp3XwOsKdkX2nDt7hfHkSaRNIqaOK/WUkVUF+Ryn6VBap1NU3iIZEAzuw3X+1kL5w5wz9I381fveQMAn/jOQ5y1/E71fOoAChQiGdDMbsONfJa6yXamyKonM/tk1HF3/2pzkyMiYZbMn81nbnpoSvVTvW0ES+bPnjJNSi2fVU+1lWRfpTaKw2JJhYhEaubEeY20N2jkfGeKDBTu/sW4EiIi0Zo5cV69XZDVTbYzVdXrycymAx8GTgKmF/a7+x+2KF0ikkKNVFtJdlXbmH0j8EpgPnA3+VHUv2xVokQknRbOHeCq809moK8HAwb6erjq/JPVPtHmqh1H8evufoGZLXD3G8zsH9G03yIdSSPnO0+1JYrJ4OeYmb0OOAKY1ZIUiYhIqlRbolhpZkcCVwCrgV7gT1uWKhERSY2qAoW7Xxu8vBt4deuSI9JcWmxHpHHV9noKLT24+581NzkizaP1t0Wao9o2iheL/u0lv871rBalSaQmIxtGOWv5nVx8+4tT5h7S+tsizVFt1dNXirfN7Mvk2ypEEhVVatAoYpHmqHdSwENQW4WkQFSpQetvizRHVYHCzDaZ2cPBv83AVuDq1iZNpLKoUoMW2xFpjmq7x76j6PUe4D+DFepEEhU195AW2xFpjkrTjB8VvCydruNwM8Pdn23k4mZ2NvmSSRdwrbsvLzn+fuCzweY48MfuvrGRa0p7qTT3kEYRizSuUoliPeCAAa8Cngte9wHbgePrvbCZdQFfA94G7ADuN7PV7v5o0Wk/Bd7k7s+Z2TnASuCMeq8p7ae41DA6NsGASg2SMvfunOTy5XdmulRbaZrx4wHM7BpgdbDGNcFN+60NXvt0YJu7PxF85jCwANgXKNz93qLz15GfjFBkikKpIZfLMTQ0lHRyRPYZ2TDK9Y/s3rfgVFbH8lTb6+m0QpAAcPcfAG9q8NoDwJNF2zuCfeV8GPhBg9cUEYnNirVbp6xKCNkcy1NtY/bTZvY54Jvkq6L+AHimwWtbyD4PPdFsHvlA8T/KfpjZImARQH9/P7lcru6EjY+PN/T+NFPesqud89eueQvraFHYn6X8VhsoLgQ+D/xTsP3Pwb5G7ABmFm0fB+wsPcnMXg9cC5zj7mWDk7uvJN+GweDgoDdSBdHOVRjKW3a1c/7aNW8D6+4MDRYDfT2Zym9VVU/u/qy7X+ruc4N/lzba4wm4HzjBzI43s2nAeykZ7W1mrwJWARe5+2MNXk9EJFZL5s9mWsldtvsgY9fuPRy/9LYpU86kWaXusX/t7h83s+8RUi3k7ufVe2F332Nml5BfAKkLuM7dN5vZ4uD4NeSnMj8a+LqZAexx98F6rynSCM1EK7VaOHeAR7c8ym3bu9g5NsERPd28uHsPz+3KL/GTlcbtSlVPNwY/v9yKiwcN5GtK9l1T9PojwEdacW2RWmgmWqnXmTO6uex9QwCctfxOxiYmpxwvNG6n+feoUvfY9cHPuwv7ggWMZrr7wy1Om0hqRM0pleY/8FZTKas2WZ2ostr1KHLAecH5DwG/MLO73f2TrUuadKK03niy+gfeSmkuZaX19yhqypk0q3YcxRHu/gJwPvB/3f1UGh9wJzJF4cYzOjaBs//G08rGvsJaFpUaFtt5Jtpqv4NSaV3vI4nfo2pldaLKagPFwWZ2LPD7wPdbmB7pYHHfeGq5oWT1D7ySRm6qaS1lpTWAQb6kddX5JzPQ14OR7yZ71fknp6K0E6XacRR/Rr530j3ufr+ZvRp4vHXJkk4U942nlnaHNMxE24o5gxppe6mlGiXOqqC0BrCCLE5UWe0KdzcBNxVtPwG8q1WJks4Ud/1trTeUJP/AWzVnUCM31Uoz9xbE3ZaR1XaANKt24aITzewOM3sk2H59MKWHSNPEXb2TpXaHVs0Z1Mh3UG01StxVQe1aTZikaque/h5YAvwdgLs/bGb/CPx5qxImnSfu6p2wJ+LiUbNJVC+Vq6JpVXVKtaWCcqopZcVdFZSGasJ2U22gOMTd/y0YHV2gFe6k6eKs3im9oSQ9ajaqiqZV1Slx3FSTqAoq/j0qBN9PfOchBY061TJ77GsIpvEws3cDT7UsVSINqrbxtPiGkvSo2agqmiXzZ/OZmx6aUv3UrOqUVgfnRkstjUjzWI8sqTZQfJT8zKy/YWaj5Feee3/LUiVSpeKAcNR044oj8t0667k5xFlFEhbIoq5fOmdQlp6Mk6wKKhd8v7B6cya+u7SottfTE8BbzexQ8g3gE8B7gJ+1MG0ikUqfFp95yVm2ahOvOPigurp8xlVFUu4p94ie7gNKNMXXL54zKGuS6jFWLviOTUwysmFUwaJKlWaPPZx8aWIAuBX4cbD9aWAj8K1WJ1Daz8iGUb74vc372gL6erp5xynHcte//6KmJ85yT4ul+woqlQziqiIpl+7p3QfR092VSBVNuyoX/IGOn6erFpW6x94IzAY2AX8E/BC4AFjo7gtanDZpQyMbRlly88Z9QQLyT3ffXLe95tHBtVYJVSoZRHX3rHeai1rSPbZrMpOjdtMsKsimZQBeFlSqenq1u58MYGbXAk8Dr3L3X7Y8ZdKWVqzdyuTe0BVvp2ikqujIQ7p5afLlup7Mw6pImt0gGlXFlcVRu2m2cO7AlNJrsTSOl0mrSiWKfd+uu+8FfqogkQ7NfMKNUy1PcdVUFYUNrPr8O09q6pN5sweMaUBYvD7/zpP0fTeoUoniFDN7IXhtQE+wbYC7++EtTV2MwnrPpPXJLu4uf82cpyeqzrjUET3dkcdLe9McNd24YsH+gNCs76LZvaE0ICxe+r4bV2nhoq6o440ys7OBq8kvhXqtuy8vOW7B8XOBXcDF7v5gs9NRrvcMtLavdb034DgX0Wl2UFoyfzZLbt5YVfXT1PGd4YqranK5HEMZmTsoTVVMaV27oZnS9H1nUbXTjDedmXUBXwPOAeYAF5rZnJLTzgFOCP4tAv62FWlJYlrirEzv3OzvZuHcAVa8+xSOPCS6tAD5xt00aIeqonJVlWleu0HSo9oBd61wOrAtGKOBmQ0DC4BHi85ZAHzD3R1YZ2Z9Znasuzd1VHgS0xLHNb1zo5r53ZQ+uX7+nSexcO4AZy2/M9HZPis9UWe96iKqVKglXttDq0uFSQaKAeDJou0dwBlVnDNAk6cPKXfjPcisZZPDxTG9czM0KyhF3ayyMMVDlqsuooJBXA9JUTeykQ2jXJnbxbO3JzMRY9bF0WZp+Yf1+JnZBcB8d/9IsH0RcLq7f6zonNuAq9z9X4PtO4DPuPv6kM9bRL56iv7+/lOHh4erTsu9OyenzPUfZtpBcPHrpnHmjMpVJvfunOSWxyZ55iXn6OnGu07sPuB9n8rt4pmXDvzuj55ufGXokKZcI8z4+Di9vb0Vzyu+TrnvppbrVspvVH6qzWuteasmXWlST/4ALr79xbLHjp5uLc9/2O9Q4e8JKHusmt+rLKj3/61azfodnjdv3np3Hww7lmSJYgcws2j7OGBnHecA4O4ryc9HxeDgoA8NDVWdkCFgTtETjwGl98XdL8Nt27sqTqEwsmGUG+/YxMRk/j/umZecG7fsZc5r50yJ7lccMRr6FH3FgpOrapAdAi4LuXal4mcul6Pe72Y0+G4Kv5Ll8hbm2dtvC9//kjM0NBSan0Keqvk+68lbNelKk3ryBzCwLrxqbyD4HWnk97Aaly+/84AHjcLfU+F12LGsTldSqt7/t2rF8TucWGM2cD9wgpkdb2bTgPcCq0vOWQ18wPLeCDzf7PaJgoVzB7hn6Zv56fK3HxAkCqopjlfb+NuMtXOLGyjf8MUfsuTmjS1plCx8NwN9PZQ+t1TbsF3vAjnlvs9PfXdjU8aQZGnxonpFNcbHsYZzVPVW2pctzYI4focTK1G4+x4zu4T8WtxdwHXuvtnMFgfHrwHWkO8au41899gPxZG2csXxar74Wn7xG6n3Lq2XDJtMrpZGybD5l75w3klT3ptEu0q5z94bVJkW18f2lZxTTQkryfaRuFRqjG91+0uldi4tW9qYOH6Hk6x6wt3XkA8GxfuuKXrt5CchjNW7Tuzmxi176/ri4+qRFPakHaaam3hh/qXisQ1jE5MsuWkjsP9G00je6u05VM0AvUJA/Is37i8g19JIXU+6sibJxvhKN7J2D9StFsfvcKKBIq3OnNHNnNfOqeuLj+sJtdqieTU38XLzL02+7FNKJHEsm1kq7Jph8t/Hofu2a+n22aqbaCcMZKtGNTeyK2/dyLMveUd/T41o9YOAAkUZ9X7xcT2hVvOkXe1NPCroFB9L4ul74dwBHvjZs3z7vif3VTeFKQ2ISdd9a2W1qaL+nhbOHaDv+cdT13lA9lOgaIE4ivlhT9rdBxm90w9mbNdkTTfxqKBTegOOuwpjZMMot6wfjQwS+wLi84/v25fEOs3Fsj6QTaUhKaZAkVHNfLovN/9S90GWeF1xubaYLjNe9qlVFbnc/kDR7CrAWm+c5Uouo2MTvGbZGi48YyZ/vvDkutLSaq0oDdXy/SlIpY8CRYY16+m+8BmVej0lodwN92V3frr87WXf18xAWs+NM6qUttedb67bDlA2WITdLPtqTnl9ml0aquX7U5VdOilQCJDeKSoa7WnVjDzVc+OsphH+2/c9GRooyt0sL3ptF0P1Z6NqzW7fqeX7y3qVXbtKcsCdSEWVZm4tHnT4qdyulsx6Ws+Ns3ggWznl2l3K3Sxveaw5s+lWWvSq2QO4avn+ku6EIOEUKCTVokYOl06RXVhHpNZg0aobZ2FEe1eZhTXK7S93UwwbBFqraqYVb/a06rV8f50wUj6LFCgk9YqnV7ln6ZuntD80ulZGM26cnxvZxGuWrWHW0tt4zbI1fG5k05RzLzxjJmHK7S93Uzx6ehUrOVVQzXfW7Gk9agk87bD2RztSG4VkVjOqKaqpE49qGP/cyKZ9DdMQ3lBd+FkYC9JlFtnrqVyPrXed2PiCk9V+Z81ss6qlY0GnjJTPGgUKyaxmjJVo9Mb57fuePGBfYX9xIPjzhSdX3R223M2yr2icSL2SGl9SS+BJa8eKTqZAIZlVaaxENf3xG71xlmuQjhogWKpcOg+cRr3xQNEJkyBK86mNQjKrtC796OlWtqG73LTrjdaJ19pQXSruNavjmFZc2o9KFJJpxU/euVxu32I71fbHb7RO/MIzZk5poyjeX40kxg1UqtrRyGgppUAhbSmudUFqbahuJJ1x0MhoCaNAIW0pzkbbWhqqSyU9eWEpjYyWMGqjkLYUR3/8SgP1qtHsdDaaprSVcCQdEgkUZnaUmf3IzB4Pfh4Zcs5MM7vLzLaY2WYzuzSJtEo2tbrRtlmN0M1MZzPSpJHREiapqqelwB3uvtzMlgbbny05Zw/wKXd/0MwOA9ab2Y/c/dG4EyvZ1Iz++OUadptZRZPk5IWl1H1WwiQVKBbAvokwbwBylAQKd38KeCp4/Usz2wIMAAoUEouoht00VtE0I00aGS1hkgoU/UEgwN2fMrNfjTrZzGYBc4H7YkibCBD9hN6MRuhmd0NtVsO4RkZLKfMaRpDW9MFmPwZeGXLocuAGd+8rOvc5dz+gnSI41gvcDfyFu6+KuN4iYBFAf3//qcPDw3WnfXx8nN7e3rrfn2bKW/Uuvv3FsscWvX4a1z+ym90v79837SC4+HXTOHNGd8XPvnfnZM3vr5S/ej4zLfR7mbx58+atd/fBsGMtCxRRzGwrMBSUJo4Fcu5+QCWomXUD3wfWuvtXq/38wcFBf+CBB+pOXy6Xa9uF3pW36p21/M7QJ/SBvh7uWfrmhkoElT47TDX5y+pgOf1eJs/MygaKpKqeVgMfBJYHP28tPcHMDPgHYEstQUKkWSo17DZSRdOqNg5VG0krJBUolgPfNbMPA9uBCwDMbAZwrbufC5wFXARsMrOHgvdd5u5rEkivdKBWNuw20p6Q1lJDWtMljUskULj7M8BbQvbvBM4NXv8r0PhKLSINaNUTer3dUNM6xUZa0yXNoZHZIgmod6BduZ5YX/ze5hamtrJmrDYo6aW5nkQSUk9ppVwbxnO7JhnZMJrY03sax5W0WidVtalEIZIhUW0YST69d9rUH3GvI5I0BQqRDIlqw0jy6T2OSRjTpNOq2hQoRDJk4dwB+nrCB88l+fTeaSvndVpVm9ooRDLmC+edlMqJ+zppDEfa1hFpNZUoRDKm057e06jTqtpUohDJoE56ek+jTptlV4FCRKQOnRSsVfUkIiKRVKIQSZFOGsQl2aFAIZISYfMlffw7D/HF723m8+88ib46Pk9BR5pBVU8iKRE2iAvy03MsW7WJe3dOVv1ZnTZyWFpLgUIkJaIGa01M7uWWx6oPFJ02clhaS4FCJCUqDdZ65qXqV6PstJHD0loKFCIpETaIq9jR06tfnqXTJumT1lKgEEmJwojrsLmcerq7eNeJ4XM8hem0kcPSWur1JJIihUFcYT2W+p5/vKbPgfpGDndyb6lOznuURAKFmR0FfAeYBfwH8Pvu/lyZc7uAB4BRd39HXGkUSVLYqN9crvpAUe4zKunkJU07Oe+VJFX1tBS4w91PAO4Itsu5FNgSS6pEOlwn95bq5LxXklSgWADcELy+AVgYdpKZHQe8Hbg2nmSJdLZO7i3VyXmvxNyr73LXtIuajbl7X9H2c+5+ZMh5NwNXAYcBn46qejKzRcAigP7+/lOHh4frTt/4+Di9vb11vz/NlLfsiiN/n8rtCu2Ge/R04ytDh7Tsumn4v2tV3tOQt2rMmzdvvbsPhh1rWRuFmf0YeGXIocurfP87gJ+7+3ozG6p0vruvBFYCDA4O+tBQxbeUlcvlaOT9aaa8ZVcc+bviiNHQRZGuWHAyQy2sp4/KW1wNzK3Kezv8XrYsULj7W8sdM7P/NLNj3f0pMzsW+HnIaWcB55nZucB04HAz+6a7/0GLkizS8dK2zkKcDcxpy3uaJNU9djXwQWB58PPW0hPcfRmwDCAoUXxaQUKk9dK0zkJUA3Mr0pimvKdJUo3Zy4G3mdnjwNuCbcxshpmtSShNIpIyamBOh0RKFO7+DPCWkP07gXND9ueAXMsTJiKpMqOvh9GQoKCpSOKlKTxEJLU0FUk6aAoPEUktNTCngwKFiKSaGpiTp6onERGJpBKFiEyhGVSllAKFiOwT9wyqhaA0OjbBwLo7FZRSSlVPIrJPnDOoFoJSoftrISiNbBht+rWkMQoUIrJPnAPcNK13dihQiMg+ca61rVHX2aFAISL7xDnALc6gJI1RoBCRfRbOHeCq809moK8HAwb6erjq/JNb0sCsUdfZoV5PIjJF6QC3kQ2jnLX8zqZ3ly0edT06NsGAuuKmlgKFiJTV6u6yhaDUDov7tDNVPYlIWeqZJKBAISIR1DNJQIFCRCKoZ5KAAoWIRFDPJIGEAoWZHWVmPzKzx4OfR5Y5r8/MbjazfzezLWb2W3GnVaSTxdldVtIrqV5PS4E73H25mS0Ntj8bct7VwO3u/m4zmwYcEmciRUTrQUhyVU8LgBuC1zcAC0tPMLPDgd8B/gHA3Xe7+1hM6RMRkYC5e/wXNRtz976i7efc/ciSc94ArAQeBU4B1gOXuvuLZT5zEbAIoL+//9Th4eG60zc+Pk5vb2/d708z5S272jl/ylvy5s2bt97dB0MPuntL/gE/Bh4J+bcAGCs597mQ9w8Ce4Azgu2rgSurufapp57qjbjrrrsaen+aKW/Z1c75U96SBzzgZe6pLWujcPe3ljtmZv9pZse6+1Nmdizw85DTdgA73P2+YPtm8m0ZIiISo6TaKFYDHwxefxC4tfQEd///wJNmVuiH9xby1VAiIhKjpNoojga+C7wK2A5c4O7PmtkM4Fp3Pzc47w3AtcA04AngQ+7+XBWf/wvgZw0k8Rjg6Qben2bKW3a1c/6Ut+T9mrv/StiBRAJF2pnZA16uUSfjlLfsauf8KW/pppHZIiISSYFCREQiKVCEW5l0AlpIecuuds6f8pZiaqMQEZFIKlGIiEgkBQoREYmkQFGGmb3BzNaZ2UNm9oCZnZ50mprJzD5mZlvNbLOZ/WXS6Wk2M/u0mbmZHZN0WprFzFYEU+4/bGb/ZGZ9SaepUWZ2dvB7uC2YSbotmNlMM7srWB5hs5ldmnSaGqFAUd5fAl909zcAfxpstwUzm0d+zq3Xu/tJwJcTTlJTmdlM4G3kB3O2kx8Br3P31wOPAcsSTk9DzKwL+BpwDjAHuNDM5iSbqqbZA3zK3V8LvBH4aJbzpkBRngOHB6+PAHYmmJZm+2Ngubv/F4C7h821lWV/BXyG/P9h23D3H7r7nmBzHXBckulpgtOBbe7+hLvvBobJP8Bknrs/5e4PBq9/CWwBMruohwJFeR8HVpjZk+SfuDP99FbiROC3zew+M7vbzE5LOkHNYmbnAaPuvjHptLTYHwI/SDoRDRoAniza3kGGb6blmNksYC5wX4VTUyupFe5Swcx+DLwy5NDl5Cch/IS732Jmv09+AaWyM+KmTYW8HQwcSb5IfBrwXTN7tWekr3SFvF0G/G68KWqeqLy5+63BOZeTr9r4VpxpawEL2ZeJ38FqmVkvcAvwcXd/Ien01EvjKMows+eBPnd3MzPgeXc/vNL7ssDMbidf9ZQLtv8f8EZ3/0WiCWuQmZ0M3AHsCnYdR77K8PRgNuLMM7MPAouBt7j7rkrnp5mZ/RbwBXefH2wvA3D3qxJNWJOYWTfwfWCtu3816fQ0QlVP5e0E3hS8fjPweIJpabYR8nnCzE4kPztvFma3jOTum9z9V919lrvPIl+V8ZttFCTOJr+2/HlZDxKB+4ETzOx4M5sGvJf8EgSZFzxc/gOwJetBAjq86qmCPwKuNrODgZcIllltE9cB15nZI8Bu4INZqXbqcH8DvAL4Uf4+xDp3X5xskurn7nvM7BJgLdAFXOfumxNOVrOcBVwEbDKzh4J9l7n7muSSVD9VPYmISCRVPYmISCQFChERiaRAISIikRQoREQkkgKFiIhEUqCQtmNm4yH7FpvZB5JIT7MU8mVmM8zs5grnftzMDoknZdLu1D1W2o6Zjbt7b8zXNPJ/Ty/X+L6Diyb6q3Ru1fkys/8ABt296oGUZtbl7nurPV86h0oU0hHM7Atm9ungdc7MvmRm/2Zmj5nZbwf7u4I1H+4P1nz4n8H+XjO7w8weNLNNZrYg2D8rWG/g68CDwMySa/5H0XX+zcx+Pdh/vZl91czuAr5kZq8xs9vNbL2Z/YuZ/UZw3vFm9pMgPVcWfe6sYLBkIc1fDtL1cLDOyP8CZgB3BdfAzC4MznnEzL5U9FnjZvZnZnYf8Fut+fYl6zQyWzrVwe5+upmdC3ye/ISPHyY/p9dpZvYK4B4z+yH5GU5/z91fsPxCSOvMrDDVxGzgQ+7+J2Wu80JwnQ8Afw28I9h/IvBWd99rZncAi939cTM7A/g6+SlWrgb+1t2/YWYfLfP5i4DjgbnBSOej3P1ZM/skMM/dnzazGcCXgFOB54AfmtlCdx8BDgUecfc/rf0rlE6hQCGdalXwcz0wK3j9u8DrzezdwfYRwAnk54z632b2O8DL5KfC7g/O+Zm7r4u4zreLfv5V0f6bgiDRC5wJ3BRMywH5aTogPw3Eu4LXN5K/2Zd6K3BNofrK3Z8NOec0IFeY9NHMvgX8Dvk5v/aSn91UpCwFCulU/xX83Mv+vwMDPubua4tPNLOLgV8BTnX3yaD+f3pw+MUK1/EyrwvvOwgYC1ZSrPT+MFblOeW8pHYJqURtFCL7rQX+OJgeGjM70cwOJV+y+HkQJOYBv1bDZ76n6OdPSg8GaxT81MwuCK5pZnZKcPge8jOqAry/zOf/EFgcTF6JmR0V7P8lcFjw+j7gTWZ2jOWXH70QuLuGPEiHU6CQdnSIme0o+vfJKt93LfAo8GDQWPx35Esb3wIGzewB8jfsf68hLa8IGoovBT5R5pz3Ax82s43AZvYvB3op+bWW7ycfrMqleTvwcPD+9wX7VwI/MLO73P0p8is03gVsBB4sLIIkUg11jxVpkXq6qIqkkUoUIiISSSUKERGJpBKFiIhEUqAQEZFIChQiIhJJgUJERCIpUIiISKT/BlTTC9AX5FOpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": "2021-02-02T06:54:51.022657Z",
     "iopub.status.busy": "2021-02-02T06:54:51.021131Z",
     "iopub.status.idle": "2021-02-02T06:54:51.027933Z",
     "shell.execute_reply": "2021-02-02T06:54:51.028907Z"
    },
    "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": "2021-02-02T06:54:51.033503Z",
     "iopub.status.busy": "2021-02-02T06:54:51.032053Z",
     "iopub.status.idle": "2021-02-02T06:54:51.160002Z",
     "shell.execute_reply": "2021-02-02T06:54:51.160950Z"
    }
   },
   "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:                Tue, 02 Feb 2021   Deviance:                       7.2134\n",
      "Time:                        06:54:51   Pearson chi2:                     71.2\n",
      "No. Iterations:                    25                                         \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": "2021-02-02T06:54:51.165268Z",
     "iopub.status.busy": "2021-02-02T06:54:51.163984Z",
     "iopub.status.idle": "2021-02-02T06:54:51.362613Z",
     "shell.execute_reply": "2021-02-02T06:54:51.363004Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Residual')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAffUlEQVR4nO3df5xddX3n8dcnYUImjDj8atqMYKBKFAgyy0AtbOsMZRssFqZYRbRWtrZZumoBMW2AVX5ol9AI2MeubaUV6bI8GvkZVJSgwNgWCpKQQIAYtFLQiV2FdIQJEzMkn/3j3hsmk3vnnnvvOed7fryfj0cemTlz55zv99y538/5/jZ3R0REymdW6ASIiEgYCgAiIiWlACAiUlIKACIiJaUAICJSUvuETkArDj74YF+4cGHi19m2bRv77bdf4tdJk/KUD0XMExQzX3nK07p1615w90OmH89VAFi4cCFr165N/DojIyMMDg4mfp00KU/5UMQ8QTHzlac8mdlz9Y6rCUhEpKQUAERESkoBQESkpBQARERKSgFARKSkcjUKSCQNq9ePsnLNZraMTbCgt5tlSxYx3N8XOlkisQsWAMxsLvCPwL7VdNzm7peFSo8IVAr/i+/YyMTkTgBGxya4+I6NAAoCUjghm4B+Dpzi7m8DjgNOM7O3B0yPCCvXbN5d+NdMTO5k5ZrNgVIkkpxgNQCvbEQwXv22q/pPmxNIUFvGJlo6LpJnFnJDGDObDawD3gR83t3/rM5rlgJLAebPn3/8qlWrEk/X+Pg4PT09iV8nTcpTNBeNvMKL2/f+TBw017hmcF6s16qniO8TFDNfecrT0NDQOncfmH48aADYnQizXuBO4GPu/mSj1w0MDLiWgmiP8hTN9D4AgO6u2Vx11uJU+gCK+D5BMfOVpzyZWd0AkIlhoO4+BowAp4VNiZTdcH8fV521mL7ebgzo6+1OrfAXSVvIUUCHAJPuPmZm3cCpwNWh0iNSM9zfpwJfSiHkPIBfAv6+2g8wC7jF3b8WMD0iIqUSchTQE0B/qOuLiJRdJvoAREQkfQoAIiIlpQAgIlJSCgAiIiWlACAiUlIKACIiJaUAICJSUgoAIiIlpQAgIlJSCgAiIiWlACAiUlIKACIiJRVyNVARadHq9aOsXLOZLWMTLOjtZtmSRVq6WtqmACCSE9N3Kxsdm+DiOzYCKAhIW9QEJJITK9ds3mOrSoCJyZ2sXLM5UIok7xQARHJiy9hES8dFmlEAEMmJBb3dLR0XaUYBQCQnli1ZRHfX7D2OdXfNZtmSRYFSJHmnTmCRnKh19GoUkMRFAUAkR4b7+1TgS2zUBCQiUlIKACIiJaUAICJSUgoAIiIlpU5gkRLT2kLlpgAgUlJaW0jUBCRSUlpbSBQAREpKawuJAoBISWltIVEAECkprS0k6gQWKSmtLSQKACIlprWFyk0BQHJj+pj10w/byWDoRInkmPoAJBdqY9ZHxyZwKmPWb3xyB6vXj4ZOmkhuBQsAZnaomT1gZpvM7CkzOz9UWiT76o1Z37ELjVkX6UDIJqBXgYvc/TEzex2wzsy+6e5PB0yTZJTGrIvEL1gNwN1/7O6PVb9+GdgEqDdK6tKYdZH4ZaIPwMwWAv3AI4GTIhlVb8z6nFlozLpIB8zdwybArAf4NvDn7n5HnZ8vBZYCzJ8///hVq1Ylnqbx8XF6enoSv06aipCnh7ZMcvszk7y43TlorvGuw3YxdES+8zRdEd6neoqYrzzlaWhoaJ27D0w/HjQAmFkX8DVgjbtf2+z1AwMDvnbt2sTTNTIywuDgYOLXSZPylA9FzBMUM195ypOZ1Q0AIUcBGfBFYFOUwl9EROIVsg/gZOCDwClmtqH677cCpkdEpFSCDQN1938GLNT1RUTKLhOjgEREJH0KACIiJaUAICJSUgoAIiIlpQAgIlJSCgAiIiWlACAiUlIKACIiJaUAICJSUtoTWApl+r7By5Ys0qbn6L5IfQoAUgir149y+VeeYmxicvex0bEJLr5jI0CpC7vafsq1LTV1X6RGTUCSe7UCbmrhXzMxubP0+wbX209Z90VAAUAKoF4BN1XZ9w3WfsrSiAKA5F6zgqzs+wZrP2VpRAFAcm+mgqy7a3bp9w2ut5+y7ouAAoAUQL0CDuCAeV1cddbi0nd0Dvf3cdVZi+nr7caAvt5u3RcBNApICqBWkGmYY2PD/X26H7IXBQApBBVwxaD5CulSABCRTNB8hfSpD0BEMkHzFdKnACAimaD5CulTABCRTNB8hfQpAIhkwOr1o5y84n7OvWcbJ6+4n9XrR0MnKXWar5A+dQKLBKbOzwoN502fAoBIYDN1fpat8NNw3nSpCUgkMHV+SigKACKBqfNTQlEAEAlMnZ8SivoARAKb2vk5OjZBnzo/JSUKACIZUOv8HBkZYXBwMHRypCTUBCQiUlIz1gDM7OMz/dzdr403OSLp0uqTUmbNmoBel0oqRALQBCwpuxkDgLtfkVZCROpJ8gldE7Ck7CJ1ApvZXODDwNHA3Npxd/+DhNIlkvgTuiZgSdlFHQV0E/BdYAlwJfABYFOnFzezG4B3AT9x92M6PZ/k2/Sn/W0/fzXRJ/QFvd2M1insNQFLyiLqKKA3ufsngW3u/vfA6cDiGK5/I3BaDOeRnKs97Y+OTeBUnvbHJibrvjauJ/S4JmDVVvI8fPndpV3JU/Ipag2g9kkcM7NjgH8HFnZ6cXf/RzPr+DySf/Xa4xuJ6wk9jtUn1ZEseWbu3vxFZn8I3A4cC3wJ6AE+5e5/03ECKgHga42agMxsKbAUYP78+cevWrWq00s2NT4+Tk9PT+LXSVPW83TuPdsivW7OLDj3mDmctKArE3m6aOQVXty+92fooLnGNYPzWj5fFvKUhCLmK095GhoaWufuA9OPRwoASWoWAKYaGBjwtWvXJp6mIs7GzHqeTl5xf932+APmdTFvzj51n9CzkKfDl99NvU+QAc+uOL3l82UhT0koYr7ylCczqxsAoo4C+lS94+5+ZacJE4FKe/zUphSotMdf9ttHZ7opRR3JkmdRO4G3Tfm3E3gnMfQBiNQM9/dx1VmL6evtxoC+3m6uOmtxpgt/0Eqekm+RagDufs3U783ss8BXOr24mf0DMAgcbGY/Ai5z9y92el7JpzzuBqVtDCXP2l0NdB5wRKcXd/dzOj2HZFsZ1trJY+ASgeh9ABthd1/XbOAQKhPCRBrSEEmRbItaA3jXlK9fBf6fu7+aQHqkjrw+RWdlrZ283j+RpDVbDvrA6pcvT/vR/maGu29NJllSk+en6CystZPn+yeStGajgNYBa6v//xR4Bvhe9et1ySZNYOan6KzLwmbneb5/RaclNMKbMQC4++HufgSwBvhtdz/Y3Q+i0iR0RxoJLLssPEW3KwtDJPN8/0JKunCut/bTxXdsVBBIWdR5ACe4+9dr37j7N4B3JJMkmSoLT9HtysLY/jzfv1DSKJxVM8uGqJ3AL5jZ/wD+L5XRQL8HvJhYqmS3RjNk8zLRKPQQybzfvxDS6LxXzSwbotYAzqEy9PNOYDXwC9VjkrAsPEXnme5f69IonFUzy4aoM4G3AucnnBZpIPRTdN4lef+KOMQ0jfWNVDPLhmbDQD/n7heY2Vdh70UP3f2MxFImknFFHWKaRuGsJTSyoVkN4Kbq/59NOiEieZOViW5xS6twVs02vBkDgLuvq/7/7doxMzsAONTdn0g4bSKZVuSOTBXO5RCpE9jMRsxs/+rM4MeBL5nZtckmTSTb1JEpeRd1GOjr3f2l6taQX3L3y8xMNQAptaJ1ZDbr0C5ih3fc8naPogaAfczsl4D3ApcmmB6R3ChSR2azDu2idnjHKY/3KGoAuJLKchAPuvujZnYElTWBREqtKG3lzTq0i9rhHac83qOo8wBuBW6d8v0PgHcnlSgRSd7U5op6G9vDax3aRe7wjkujezE6NsHhy+/OZA0xaifwkWZ2n5k9Wf3+2OrSECKSQ9PX+2mk1qGtDu/mZroXWV3wLupSEH8LXAxMAlSHgL4vqURJcWkJ4Gyo11wx3dQO7Sys7Jp19e7RdFlb8C5qH8A8d/+OmU09VtgdwR7aMsmlK+7Pfcde1sTZSbZ6/SifHnmFrfdks2qddTM13RjsdU+L1OGdlOn3qFmzWha0shroL1NdDsLMfhf4cWKpCmj1+lFufHIHO3ZVvs9DT34I7Qx3i6uT7LVAUvmI6T1qXaP1fvp6u3lw+Sl1f6coHd5JmnqPTl5xf+JrKnUqahPQR4AvAG8xs1HgAuC8pBIV0so1m3cX/jVZq7YlJWrzTLvrxcfVkai15DunJp3k5eEeRwoA7v4Ddz+VypLQbwEGgf+cYLqCKetoh1YK9XYL4Lg6Esv6HsVJy2QnLw/3uNlqoPtTefrvA+4CvlX9/hNUloS4OekEpi2NpXCzqJXmmXYL4Lhmzpb1PYqbmnSSl/V73KwGcBOwCNgI/BFwL/AeYNjdz0w4bUEsW7KIOdPuStaqbUlopVBv90k+rieiPFStRfKgWSfwEe6+GMDM/g54ATjM3V9OPGWBDPf38fSmp7n7+dmlGu3QylN1J0/ycTwR1X7/03c9ztbtXpr3KIq8rUUjYTULAJO1L9x9p5k9W+TCv+akBV1c8v7B0MlIVSuFehaGBA7399H7s+8xODiY2jWzLo9r0XRCwa5zzQLA28zsperXBnRXvzfA3X3/RFMnqWm1UM9622YZ5XEtmnbVC3YXfnkDa5/bymeGFwdOXX402xBm5mltUihpF+p6gotXmUZH1Qt2Dtz88PMMvPFA/R1FFHUegEis2p1LII2Vab2eRkHNQfNBWqAAIEFoMlf8yjQ6aqagVsQaT1IUACSIMjVXpCUPE4/ismzJIqzBz4pY40lK1LWARGKlyVzJSLIfJ0t9NsP9fax9bis3P/z8HouuFbXGkxTVAAIp+7LIeWmuaOd9KuJ7m8U+m88ML+a6s4/bXeM5YF4X++4ziwu/vKEw9z1pCgABZPHDlLY8NFe08z4V9b3Nap/NcH8fDy4/hevOPo7tk7sYm5gs1H1PWtAAYGanmdlmM/u+mS0PmZY0ZfXDlLbah/fZFafz4PJTMlX4Q3vvU1Hf26z32RT1victWAAws9nA54F3AkcB55jZUaHSk6asf5ikYqY9Xhs9WRb1vc36ENOi3vekhawBnAh8v7rU9A5gFVDIBeamy/qHSSpmej8aNS8U9b3Nep9NUe970sx9pi2hE7xwZVex09z9D6vffxD4FXf/6LTXLQWWAsyfP//4VatWJZ628fFxenp6Ejv/Q1sm99h1DGDOLDj3mDmctKArkWsmnacQ2s3TQ1smuf2ZSV7c7hw013j3kV1173u992mqg+Ya1wzOa/o7rby3WX6fot63evSZCmtoaGiduw9MPx5yGGi9Ybx7RSN3vx64HmBgYMDTWPxrZGQk0UXGBoGjUh5Sl3SeQmgnT6vXj3LTfa9tJ/niduemTTs56q1H7XX/B6m8Txd8eUPdc23d7ntdv/Y77b63WX6fBoFL2vxdfabal+Tw25AB4EfAoVO+fwOwJVBaUqfF1MJodcG04f4+Vq7Z3NKcBb23e3toyyQXXnkv//FKZYHh3u4uLj/j6FjvUxHve9IrvIbsA3gUeLOZHW5mc4D3AV8JmB7Jgalj7C8aeaXlYX7tdBZmuf07D3MOVq8f5YaNO3YX/gBjE5N8/JYNmUxvliQ9uilYAHD3V4GPAmuATcAt7v5UqPRI9k0fY//idm95rHc7nYVZnbOQlzkHK9ds5tU6XY27HK74qj7yM0l6dFPQpSDc/evA10OmQfIjjvXu293NLIvNC3lZ/3+mwmpqrUD2lvSSKZoJLLtlvTkhjqehrD7NtyMvY981FLN9STc/ajE4AfKxnWBcT0NZfJpvRxz3I40F3pYtWdRwJFVvdzJDNIsi6e1XFQAEyEdzQieb0Scp1CqZnd6PtIL+cH8fdz20kQd+uOffV9cs4/Izjo7tOkWV5AOLAoAA2WtOqFeoAsztmrW7wNpvH/jzwM03IWtOnT4dphn0P3T0XM486c2ZWU5aKhQABMjW+vz1CtVltz0ODpO7XhtOMtlghm6aQtecOnk6TDvoF6XprUjUCSxAtsa61ytUJ3f6HoU/wI5d4fd/zVrNqRVRh8RmfXCAtE8BQIBsjY5ppfAMXdDmeRGyKEE/rrkGD22ZVBDJIDUByW5ZqaI3ao5q9NqQstoxHUWUPoQ4mrhWrx/dY6G2LI4wKysFAElMu6Nj6hWqXbNtrz6AObMIXtAmPUwvac2CfhxNXCvXbN5rRdWsjTArKwUASUQno2MaFarTj51+2M5MFCBZqTklIY7BAXnuJyk6BQBJRKdNB40K1anHRkZGOk5nnoSYbxBHE1eWRpjJnhQAAgg1cSguUdKvp754hZpvEEcT17Ili/jTWzfs0QyUl36SolMASFkellyYSdT066kvXiHnG3TaxDXc38fTm57m7udn5/ahp6gUAFIWeuJQp6KmP8+jY2rirql1cr6816hOWtDFJe8fDJ0MmUYBIGV5/yBHTX/eR8fEXVNbvX6UZbc9zuRO332+Zbc9Hvl8Sdeo8t4sKe1RAEhZ3ptGWkl/nkfHxF1Tu+KrT+0u/GsmdzqX3PFE20Nj46pR5b1ZUtqnmcApy9KSC+3Ie/qjirum1mjjk1cmd0WaFZvkTO2ktx2U7FINIGV5bxrJe/qjSrOmdtEt0ZqCkqpRNZp1nZdmSWmfAkAAeW4agfynP4q4m1x6u7sYm6hfC9jpvrvJpbets7dv9fpRDKizZS+zzDh8+d2FDfKiAJB56pxrXyf3LkpNZ/X6Ua746lO7m3d6u7u4/Iyj617j8jOObrgrFlSaXC748gYOmmt88vWjqb3HK9dsrlv4QyUwgfoEikwBIMPUOde+OO7dTDWd6aN6AMYmJll2a/3mnOH+PtY+t5WbH36+YYEL8OJ2T/U9jtrMk6ehyhKdOoEzrFHn3EW3PK5ldZtIumNz5ZrNe43qgcpidY2u8ZnhxVx39nHMNpvx3Gl2wMaxpk+RlG3vAwWADGv0gdvp3tHa7GWQ9HyLmc4z08+G+/u45r1v22skVSvniCJqQVZvVFej8JSXocrtimvvgzxRAMiwKB84Dderb6aNWuJ4ypvpvWn2vk0d0tnuOWbSSkFWb3jpB95+WCmG+k5XxuGwCgAZVu/prJ4yVM1b1Wi+wtBbDonlKW/ZkkWVPQqm6ZplLFuyqGmQGe7v48Hlp/C5s4+LvbBttSCrpeXZFafz4PJT+Mzw4szsDpemvM/Sb4c6gTNs+kiUWWa7R2ZMVfSqeTsajeKJa4Zv7bX1RgEBkTugp6fzwLnGJ8/srLCNoyArw1Df6fI+S78dCgAZN/WDOH1kC5Sjat6ueoXYhQ2GYrbzlNeokDx5xf0tBZmp5xkZGWGww4I3VEGW9yHLRVjAsFVqAsqRLG3cnldpbOIeuikhxHIdRehALePnSzWAnClj1TxOnT7lRXnKDd2UEGK5jrwvc15Tts+XAoCUSieFY9TJZVloSki7IAtd65H2KABI6bRbOEZ9yi3LgnlTha71SHsUAEQiauUpt2xNCVmo9Ujr1AksElEaHch5VcYO1CJQDUAkIj3lzqxstZ4iUAAQiaiMbftSbEECgJm9B7gceCtworuvTepanUxOycrElqykI01J5DmOc+blKbeVvQqkvELVAJ4EzgK+kORFOlkTPitr8WclHWlKIs9luo+t7lUg5RWkE9jdN7l74kvsdbK6X1ZWBsxKOtKURJ7bOWde14ZvZ6+CLMnrfc8j8zqLi6V2cbMR4BMzNQGZ2VJgKcD8+fOPX7VqVeTzn3vPtoY/u/G0/Rr+bHx8nI/+c+NNO2b63bi1m4fpxsfH6enpiSNJiYua51by1Op9fGjLJDc+uYMdu147NmcWnHvMHE5a0BXpmu2I432aKa+Q7t9vTdR8hbrv7cjTZ2poaGiduw9MP55YE5CZfQv4xTo/utTd74p6Hne/HrgeYGBgwAcHByOnoe/h++tOTunr7Wam84yMjNDXu6ut341bu3mYbmRkJNV0dyJqnlvJU6v38dIV9+9RCAHs2AV3Pz+bS94f7ZrtiON9apRXSP/vtyZqvkLd93bk6TPVSGJNQO5+qrsfU+df5MK/U50sihViQa0spyNNSeS51XPmeWmDZnsVZFme73seFXoYaCfD9rIy5C8r6UhTEnlu9Zx5Xtpgpr0Ksv53k+f7nkehhoH+DvC/gEOAu81sg7svSeJanQzby8qQv6ykI01J5LmVc+Z90lde/2byft/zJkgAcPc7gTtDXFskiqzVvMoyFyRr973oCt0EJNKJrDxFl2kOA2TnvpeBFoMTybgyzgWRdCgAiGScRsZIUtQEJJJxRR0ZE7pfI/T1s0A1AJGMK+JckNCbyIe+flYoAIhkXBE3WwndrxH6+lmhJiCRHCjayJjQ/Rqhr58VqgGISOpCb68Z+vpZoQAgIqkL3a8R+vpZoSYgEUld6Bm/oa+fFQoAIhJE6H6N0NfPAjUBiYiUlAKAiEhJqQlIRFqmWbTFoAAgIi0p2+qkRaYmIBFpiWbRFocCgIi0RLNoi0MBQERaolm0xaEAICIt0Sza4lAnsIi0RLNoi0MBQERaplm0xaAmIBGRklIAEBEpKQUAEZGSUgAQESkpBQARkZIydw+dhsjM7KfAcylc6mDghRSukyblKR+KmCcoZr7ylKc3uvsh0w/mKgCkxczWuvtA6HTESXnKhyLmCYqZryLkSU1AIiIlpQAgIlJSCgD1XR86AQlQnvKhiHmCYuYr93lSH4CISEmpBiAiUlIKACIiJaUAUIeZHWdmD5vZBjNba2Ynhk5TXMzsY2a22cyeMrO/CJ2euJjZJ8zMzezg0GnplJmtNLPvmtkTZnanmfWGTlO7zOy06t/b981seej0xMHMDjWzB8xsU/VzdH7oNLVLAaC+vwCucPfjgE9Vv889MxsCzgSOdfejgc8GTlIszOxQ4L8Az4dOS0y+CRzj7scCzwAXB05PW8xsNvB54J3AUcA5ZnZU2FTF4lXgInd/K/B24CN5zZcCQH0O7F/9+vXAloBpidMfAyvc/ecA7v6TwOmJy3XAn1J533LP3e9191er3z4MvCFkejpwIvB9d/+Bu+8AVlF5AMk1d/+xuz9W/fplYBOQy80RFADquwBYaWY/pPKUnMsnsDqOBH7NzB4xs2+b2QmhE9QpMzsDGHX3x0OnJSF/AHwjdCLa1Af8cMr3PyKnBWUjZrYQ6AceCZyUtpR2RzAz+xbwi3V+dCnwG8CF7n67mb0X+CJwaprpa1eTfO0DHECl2noCcIuZHeEZHwvcJE+XAL+Zboo6N1Oe3P2u6msupdLccHOaaYuR1TmW6b+1VphZD3A7cIG7vxQ6Pe3QPIA6zOxnQK+7u5kZ8DN337/Z72Wdmd1DpQlopPr9vwJvd/efBk1Ym8xsMXAf8Er10BuoNNed6O7/HixhMTCzDwHnAb/h7q80e30WmdmvApe7+5Lq9xcDuPtVQRMWAzPrAr4GrHH3a0Onp11qAqpvC/CO6tenAN8LmJY4raaSH8zsSGAO+VnNcC/uvtHdf8HdF7r7QipNDP+pAIX/acCfAWfktfCvehR4s5kdbmZzgPcBXwmcpo5VHwq/CGzKc+EPJW4CauKPgL80s32A7cDSwOmJyw3ADWb2JLAD+FDWm39K6n8D+wLfrJQ1POzu54VNUuvc/VUz+yiwBpgN3ODuTwVOVhxOBj4IbDSzDdVjl7j718MlqT1qAhIRKSk1AYmIlJQCgIhISSkAiIiUlAKAiEhJKQCIiJSUAoDkipmN1zl2npn9foj0xKWWLzNbYGa3NXntBWY2L52USZFpGKjkipmNu3tPytc0Kp+VXS3+3j5TFnVr9trI+TKzfwMG3D3yJD4zm+3uO6O+XspBNQDJPTO73Mw+Uf16xMyuNrPvmNkzZvZr1eOzq+vsP1pdZ/+/VY/3mNl9ZvaYmW00szOrxxdW13v/K+Ax4NBp1/y3Kdf5jpm9qXr8RjO71sweAK42s182s3vMbJ2Z/ZOZvaX6usPN7F+q6fn0lPMurE7Uq6X5s9V0PVHdy+FPgAXAA9VrYGbnVF/zpJldPeVc42Z2pZk9AvxqMndf8kwzgaWI9nH3E83st4DLqCzk92EqazqdYGb7Ag+a2b1UVqv8HXd/ySqbyTxsZrXlChYB/9Xd/3uD67xUvc7vA58D3lU9fiRwqrvvNLP7gPPc/Xtm9ivAX1FZjuMvgb929/9jZh9pcP6lwOFAf3VW7YHuvtXMPg4MufsLZrYAuBo4HvgP4F4zG3b31cB+wJPu/qnWb6GUgQKAFNEd1f/XAQurX/8mcKyZ/W71+9cDb6ayftD/NLNfB3ZRWa54fvU1z7n7wzNc5x+m/H/dlOO3Vgv/HuAk4Nbqkg5QWeIBKssJvLv69U1UCvHpTgX+ptaM5O5b67zmBGCktqCfmd0M/DqVdZ92UlmtUqQuBQApop9X/9/Ja3/jBnzM3ddMfaGZnQscAhzv7pPV9vW51R9va3Idb/B17fdmAWPVneWa/X49FvE1jWxXu7/MRH0AUhZrgD+uLuOLmR1pZvtRqQn8pFr4DwFvbOGcZ0/5/1+m/7C6RvyzZvae6jXNzN5W/fGDVFbHBPhAg/PfC5xXXZQQMzuwevxl4HXVrx8B3mFmB1tlC8ZzgG+3kAcpMQUAyZt5ZvajKf8+HvH3/g54Gnis2sn6BSq1g5uBATNbS6Ug/m4Ladm32sF6PnBhg9d8APiwmT0OPMVrWyKeT2Uv2UepBKFGaX4eeKL6+++vHr8e+IaZPeDuP6ayY90DwOPAY7UNZUSa0TBQkTa0MxRTJGtUAxARKSnVAERESko1ABGRklIAEBEpKQUAEZGSUgAQESkpBQARkZL6/wd+wHhyK/nXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}