{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weighted Generalized Linear Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:55.716556Z",
     "iopub.status.busy": "2022-11-02T17:04:55.716267Z",
     "iopub.status.idle": "2022-11-02T17:04:56.969964Z",
     "shell.execute_reply": "2022-11-02T17:04:56.968805Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weighted GLM: Poisson response data\n",
    "\n",
    "### Load data\n",
    "\n",
    "In this example, we'll use the affair dataset using a handful of exogenous variables to predict the extra-marital affair rate. \n",
    "\n",
    "Weights will be generated to show that `freq_weights` are equivalent to repeating records of data. On the other hand, `var_weights` is equivalent to aggregating data.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:56.975084Z",
     "iopub.status.busy": "2022-11-02T17:04:56.974655Z",
     "iopub.status.idle": "2022-11-02T17:04:56.982557Z",
     "shell.execute_reply": "2022-11-02T17:04:56.981648Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "::\n",
      "\n",
      "    Number of observations: 6366\n",
      "    Number of variables: 9\n",
      "    Variable name definitions:\n",
      "\n",
      "        rate_marriage   : How rate marriage, 1 = very poor, 2 = poor, 3 = fair,\n",
      "                        4 = good, 5 = very good\n",
      "        age             : Age\n",
      "        yrs_married     : No. years married. Interval approximations. See\n",
      "                        original paper for detailed explanation.\n",
      "        children        : No. children\n",
      "        religious       : How relgious, 1 = not, 2 = mildly, 3 = fairly,\n",
      "                        4 = strongly\n",
      "        educ            : Level of education, 9 = grade school, 12 = high\n",
      "                        school, 14 = some college, 16 = college graduate,\n",
      "                        17 = some graduate school, 20 = advanced degree\n",
      "        occupation      : 1 = student, 2 = farming, agriculture; semi-skilled,\n",
      "                        or unskilled worker; 3 = white-colloar; 4 = teacher\n",
      "                        counselor social worker, nurse; artist, writers;\n",
      "                        technician, skilled worker, 5 = managerial,\n",
      "                        administrative, business, 6 = professional with\n",
      "                        advanced degree\n",
      "        occupation_husb : Husband's occupation. Same as occupation.\n",
      "        affairs         : measure of time spent in extramarital affairs\n",
      "\n",
      "    See the original paper for more details.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sm.datasets.fair.NOTE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the data into a pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:56.988995Z",
     "iopub.status.busy": "2022-11-02T17:04:56.987733Z",
     "iopub.status.idle": "2022-11-02T17:04:57.001414Z",
     "shell.execute_reply": "2022-11-02T17:04:57.000865Z"
    }
   },
   "outputs": [],
   "source": [
    "data = sm.datasets.fair.load_pandas().data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The dependent (endogenous) variable is ``affairs``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.006097Z",
     "iopub.status.busy": "2022-11-02T17:04:57.004699Z",
     "iopub.status.idle": "2022-11-02T17:04:57.043867Z",
     "shell.execute_reply": "2022-11-02T17:04:57.043317Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rate_marriage</th>\n",
       "      <th>age</th>\n",
       "      <th>yrs_married</th>\n",
       "      <th>children</th>\n",
       "      <th>religious</th>\n",
       "      <th>educ</th>\n",
       "      <th>occupation</th>\n",
       "      <th>occupation_husb</th>\n",
       "      <th>affairs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "      <td>6366.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>4.109645</td>\n",
       "      <td>29.082862</td>\n",
       "      <td>9.009425</td>\n",
       "      <td>1.396874</td>\n",
       "      <td>2.426170</td>\n",
       "      <td>14.209865</td>\n",
       "      <td>3.424128</td>\n",
       "      <td>3.850141</td>\n",
       "      <td>0.705374</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.961430</td>\n",
       "      <td>6.847882</td>\n",
       "      <td>7.280120</td>\n",
       "      <td>1.433471</td>\n",
       "      <td>0.878369</td>\n",
       "      <td>2.178003</td>\n",
       "      <td>0.942399</td>\n",
       "      <td>1.346435</td>\n",
       "      <td>2.203374</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>17.500000</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>9.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>4.000000</td>\n",
       "      <td>22.000000</td>\n",
       "      <td>2.500000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>12.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>4.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>14.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>5.000000</td>\n",
       "      <td>32.000000</td>\n",
       "      <td>16.500000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>16.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>0.484848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>5.000000</td>\n",
       "      <td>42.000000</td>\n",
       "      <td>23.000000</td>\n",
       "      <td>5.500000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>20.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>57.599991</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       rate_marriage          age  yrs_married     children    religious  \\\n",
       "count    6366.000000  6366.000000  6366.000000  6366.000000  6366.000000   \n",
       "mean        4.109645    29.082862     9.009425     1.396874     2.426170   \n",
       "std         0.961430     6.847882     7.280120     1.433471     0.878369   \n",
       "min         1.000000    17.500000     0.500000     0.000000     1.000000   \n",
       "25%         4.000000    22.000000     2.500000     0.000000     2.000000   \n",
       "50%         4.000000    27.000000     6.000000     1.000000     2.000000   \n",
       "75%         5.000000    32.000000    16.500000     2.000000     3.000000   \n",
       "max         5.000000    42.000000    23.000000     5.500000     4.000000   \n",
       "\n",
       "              educ   occupation  occupation_husb      affairs  \n",
       "count  6366.000000  6366.000000      6366.000000  6366.000000  \n",
       "mean     14.209865     3.424128         3.850141     0.705374  \n",
       "std       2.178003     0.942399         1.346435     2.203374  \n",
       "min       9.000000     1.000000         1.000000     0.000000  \n",
       "25%      12.000000     3.000000         3.000000     0.000000  \n",
       "50%      14.000000     3.000000         4.000000     0.000000  \n",
       "75%      16.000000     4.000000         5.000000     0.484848  \n",
       "max      20.000000     6.000000         6.000000    57.599991  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.048157Z",
     "iopub.status.busy": "2022-11-02T17:04:57.047068Z",
     "iopub.status.idle": "2022-11-02T17:04:57.060475Z",
     "shell.execute_reply": "2022-11-02T17:04:57.059910Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rate_marriage</th>\n",
       "      <th>age</th>\n",
       "      <th>yrs_married</th>\n",
       "      <th>children</th>\n",
       "      <th>religious</th>\n",
       "      <th>educ</th>\n",
       "      <th>occupation</th>\n",
       "      <th>occupation_husb</th>\n",
       "      <th>affairs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.111111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>3.230769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.400000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   rate_marriage   age  yrs_married  children  religious  educ  occupation  \\\n",
       "0            3.0  32.0          9.0       3.0        3.0  17.0         2.0   \n",
       "1            3.0  27.0         13.0       3.0        1.0  14.0         3.0   \n",
       "2            4.0  22.0          2.5       0.0        1.0  16.0         3.0   \n",
       "\n",
       "   occupation_husb   affairs  \n",
       "0              5.0  0.111111  \n",
       "1              4.0  3.230769  \n",
       "2              5.0  1.400000  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following we will work mostly with Poisson. While using decimal affairs works, we convert them to integers to have a count distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.065032Z",
     "iopub.status.busy": "2022-11-02T17:04:57.063918Z",
     "iopub.status.idle": "2022-11-02T17:04:57.078112Z",
     "shell.execute_reply": "2022-11-02T17:04:57.077512Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rate_marriage</th>\n",
       "      <th>age</th>\n",
       "      <th>yrs_married</th>\n",
       "      <th>children</th>\n",
       "      <th>religious</th>\n",
       "      <th>educ</th>\n",
       "      <th>occupation</th>\n",
       "      <th>occupation_husb</th>\n",
       "      <th>affairs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   rate_marriage   age  yrs_married  children  religious  educ  occupation  \\\n",
       "0            3.0  32.0          9.0       3.0        3.0  17.0         2.0   \n",
       "1            3.0  27.0         13.0       3.0        1.0  14.0         3.0   \n",
       "2            4.0  22.0          2.5       0.0        1.0  16.0         3.0   \n",
       "\n",
       "   occupation_husb  affairs  \n",
       "0              5.0      1.0  \n",
       "1              4.0      4.0  \n",
       "2              5.0      2.0  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[\"affairs\"] = np.ceil(data[\"affairs\"])\n",
    "data[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.082371Z",
     "iopub.status.busy": "2022-11-02T17:04:57.081273Z",
     "iopub.status.idle": "2022-11-02T17:04:57.089184Z",
     "shell.execute_reply": "2022-11-02T17:04:57.088609Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6775054979579014"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(data[\"affairs\"] == 0).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.093435Z",
     "iopub.status.busy": "2022-11-02T17:04:57.092354Z",
     "iopub.status.idle": "2022-11-02T17:04:57.101342Z",
     "shell.execute_reply": "2022-11-02T17:04:57.100763Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4313,  934,  488,  180,  130,  172,    7,   21,   67,    2,    0,\n",
       "          0,   17,    0,    0,    0,    3,   12,    8,    0,    0,    0,\n",
       "          0,    0,    2,    2,    2,    3,    0,    0,    0,    0,    0,\n",
       "          0,    0,    0,    0,    0,    0,    1,    1,    0,    0,    0,\n",
       "          0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,\n",
       "          0,    0,    0,    1])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.bincount(data[\"affairs\"].astype(int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Condensing and Aggregating observations\n",
    "\n",
    "We have 6366 observations in our original dataset. When we consider only some selected variables, then we have fewer unique observations. In the following we combine observations in two ways, first we combine observations that have values for all variables identical, and secondly we combine observations that have the same explanatory variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset with unique observations\n",
    "\n",
    "We use pandas's groupby to combine identical observations and create a new variable `freq` that count how many observation have the values in the corresponding row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.106245Z",
     "iopub.status.busy": "2022-11-02T17:04:57.104698Z",
     "iopub.status.idle": "2022-11-02T17:04:57.138148Z",
     "shell.execute_reply": "2022-11-02T17:04:57.137544Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(476, 5)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>affairs</th>\n",
       "      <th>rate_marriage</th>\n",
       "      <th>age</th>\n",
       "      <th>yrs_married</th>\n",
       "      <th>freq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>17.5</td>\n",
       "      <td>0.5</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   affairs  rate_marriage   age  yrs_married  freq\n",
       "0      0.0            1.0  17.5          0.5     1\n",
       "1      0.0            1.0  22.0          2.5     3\n",
       "2      0.0            1.0  27.0          2.5     1\n",
       "3      0.0            1.0  27.0          6.0     5\n",
       "4      0.0            1.0  27.0          9.0     1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data2 = data.copy()\n",
    "data2[\"const\"] = 1\n",
    "dc = (\n",
    "    data2[\"affairs rate_marriage age yrs_married const\".split()]\n",
    "    .groupby(\"affairs rate_marriage age yrs_married\".split())\n",
    "    .count()\n",
    ")\n",
    "dc.reset_index(inplace=True)\n",
    "dc.rename(columns={\"const\": \"freq\"}, inplace=True)\n",
    "print(dc.shape)\n",
    "dc.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset with unique explanatory variables (exog)\n",
    "\n",
    "For the next dataset we combine observations that have the same values of the explanatory variables. However, because the response variable can differ among combined observations, we compute the mean and the sum of the response variable for all combined observations.\n",
    "\n",
    "We use again pandas ``groupby`` to combine observations and to create the new variables. We also flatten the ``MultiIndex`` into a simple index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.143475Z",
     "iopub.status.busy": "2022-11-02T17:04:57.141940Z",
     "iopub.status.idle": "2022-11-02T17:04:57.179529Z",
     "shell.execute_reply": "2022-11-02T17:04:57.178930Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(130, 6)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rate_marriage</th>\n",
       "      <th>age</th>\n",
       "      <th>yrs_married</th>\n",
       "      <th>affairs_mean</th>\n",
       "      <th>affairs_sum</th>\n",
       "      <th>affairs_count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>17.5</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>3.900000</td>\n",
       "      <td>39.0</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>3.400000</td>\n",
       "      <td>17.0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.900000</td>\n",
       "      <td>9.0</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.333333</td>\n",
       "      <td>4.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   rate_marriage   age  yrs_married  affairs_mean  affairs_sum  affairs_count\n",
       "0            1.0  17.5          0.5      0.000000          0.0              1\n",
       "1            1.0  22.0          2.5      3.900000         39.0             10\n",
       "2            1.0  27.0          2.5      3.400000         17.0              5\n",
       "3            1.0  27.0          6.0      0.900000          9.0             10\n",
       "4            1.0  27.0          9.0      1.333333          4.0              3"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gr = data[\"affairs rate_marriage age yrs_married\".split()].groupby(\n",
    "    \"rate_marriage age yrs_married\".split()\n",
    ")\n",
    "df_a = gr.agg([\"mean\", \"sum\", \"count\"])\n",
    "\n",
    "\n",
    "def merge_tuple(tpl):\n",
    "    if isinstance(tpl, tuple) and len(tpl) > 1:\n",
    "        return \"_\".join(map(str, tpl))\n",
    "    else:\n",
    "        return tpl\n",
    "\n",
    "\n",
    "df_a.columns = df_a.columns.map(merge_tuple)\n",
    "df_a.reset_index(inplace=True)\n",
    "print(df_a.shape)\n",
    "df_a.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After combining observations with have a dataframe `dc` with 467 unique observations, and a dataframe `df_a` with 130 observations with unique values of the explanatory variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.184042Z",
     "iopub.status.busy": "2022-11-02T17:04:57.182950Z",
     "iopub.status.idle": "2022-11-02T17:04:57.193154Z",
     "shell.execute_reply": "2022-11-02T17:04:57.192609Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of rows: \n",
      "original, with unique observations, with unique exog\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(6366, 476, 130)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"number of rows: \\noriginal, with unique observations, with unique exog\")\n",
    "data.shape[0], dc.shape[0], df_a.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "In the following, we compare the GLM-Poisson results of the original data with models of the combined observations where the multiplicity or aggregation is given by weights or exposure.\n",
    "\n",
    "\n",
    "### original data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.197525Z",
     "iopub.status.busy": "2022-11-02T17:04:57.196445Z",
     "iopub.status.idle": "2022-11-02T17:04:57.272124Z",
     "shell.execute_reply": "2022-11-02T17:04:57.271356Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                affairs   No. Observations:                 6366\n",
      "Model:                            GLM   Df Residuals:                     6362\n",
      "Model Family:                 Poisson   Df Model:                            3\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -10351.\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       15375.\n",
      "Time:                        17:04:57   Pearson chi2:                 3.23e+04\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):             0.2420\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.7155      0.107     25.294      0.000       2.505       2.926\n",
      "rate_marriage    -0.4952      0.012    -41.702      0.000      -0.518      -0.472\n",
      "age              -0.0299      0.004     -6.691      0.000      -0.039      -0.021\n",
      "yrs_married      -0.0108      0.004     -2.507      0.012      -0.019      -0.002\n",
      "=================================================================================\n"
     ]
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs ~ rate_marriage + age + yrs_married\",\n",
    "    data=data,\n",
    "    family=sm.families.Poisson(),\n",
    ")\n",
    "res_o = glm.fit()\n",
    "print(res_o.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.277676Z",
     "iopub.status.busy": "2022-11-02T17:04:57.276089Z",
     "iopub.status.idle": "2022-11-02T17:04:57.285369Z",
     "shell.execute_reply": "2022-11-02T17:04:57.284815Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.078702313363233"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_o.pearson_chi2 / res_o.df_resid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### condensed data (unique observations with frequencies)\n",
    "\n",
    "Combining identical observations and using frequency weights to take into account the multiplicity of observations produces exactly the same results. Some results attribute will differ when we want to have information about the observation and not about the aggregate of all identical observations. For example, residuals do not take ``freq_weights`` into account."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.290773Z",
     "iopub.status.busy": "2022-11-02T17:04:57.289195Z",
     "iopub.status.idle": "2022-11-02T17:04:57.324679Z",
     "shell.execute_reply": "2022-11-02T17:04:57.324128Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                affairs   No. Observations:                  476\n",
      "Model:                            GLM   Df Residuals:                     6362\n",
      "Model Family:                 Poisson   Df Model:                            3\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -10351.\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       15375.\n",
      "Time:                        17:04:57   Pearson chi2:                 3.23e+04\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):             0.9754\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.7155      0.107     25.294      0.000       2.505       2.926\n",
      "rate_marriage    -0.4952      0.012    -41.702      0.000      -0.518      -0.472\n",
      "age              -0.0299      0.004     -6.691      0.000      -0.039      -0.021\n",
      "yrs_married      -0.0108      0.004     -2.507      0.012      -0.019      -0.002\n",
      "=================================================================================\n"
     ]
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs ~ rate_marriage + age + yrs_married\",\n",
    "    data=dc,\n",
    "    family=sm.families.Poisson(),\n",
    "    freq_weights=np.asarray(dc[\"freq\"]),\n",
    ")\n",
    "res_f = glm.fit()\n",
    "print(res_f.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.329809Z",
     "iopub.status.busy": "2022-11-02T17:04:57.328287Z",
     "iopub.status.idle": "2022-11-02T17:04:57.337082Z",
     "shell.execute_reply": "2022-11-02T17:04:57.336539Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0787023133632"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_f.pearson_chi2 / res_f.df_resid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### condensed using ``var_weights`` instead of ``freq_weights``\n",
    "\n",
    "Next, we compare ``var_weights`` to ``freq_weights``. It is a common practice to incorporate ``var_weights`` when the endogenous variable reflects averages and not identical observations.\n",
    "I do not see a theoretical reason why it produces the same results (in general).\n",
    "\n",
    "This produces the same results but ``df_resid``  differs the ``freq_weights`` example because ``var_weights`` do not change the number of effective observations.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.342189Z",
     "iopub.status.busy": "2022-11-02T17:04:57.340681Z",
     "iopub.status.idle": "2022-11-02T17:04:57.365395Z",
     "shell.execute_reply": "2022-11-02T17:04:57.364849Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                affairs   No. Observations:                  476\n",
      "Model:                            GLM   Df Residuals:                      472\n",
      "Model Family:                 Poisson   Df Model:                            3\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -10351.\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       15375.\n",
      "Time:                        17:04:57   Pearson chi2:                 3.23e+04\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):             0.9754\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.7155      0.107     25.294      0.000       2.505       2.926\n",
      "rate_marriage    -0.4952      0.012    -41.702      0.000      -0.518      -0.472\n",
      "age              -0.0299      0.004     -6.691      0.000      -0.039      -0.021\n",
      "yrs_married      -0.0108      0.004     -2.507      0.012      -0.019      -0.002\n",
      "=================================================================================\n"
     ]
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs ~ rate_marriage + age + yrs_married\",\n",
    "    data=dc,\n",
    "    family=sm.families.Poisson(),\n",
    "    var_weights=np.asarray(dc[\"freq\"]),\n",
    ")\n",
    "res_fv = glm.fit()\n",
    "print(res_fv.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dispersion computed from the results is incorrect because of wrong ``df_resid``. \n",
    "It is correct if we use the original ``df_resid``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.369727Z",
     "iopub.status.busy": "2022-11-02T17:04:57.368657Z",
     "iopub.status.idle": "2022-11-02T17:04:57.375148Z",
     "shell.execute_reply": "2022-11-02T17:04:57.374611Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(68.45488160512008, 5.0787023133632)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_fv.pearson_chi2 / res_fv.df_resid, res_f.pearson_chi2 / res_f.df_resid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### aggregated or averaged data (unique values of explanatory variables)\n",
    "\n",
    "For these cases we combine observations that have the same values of the explanatory variables. The corresponding response variable is either a sum or an average.\n",
    "\n",
    "#### using ``exposure``\n",
    "\n",
    "If our dependent variable is the sum of the responses of all combined observations, then under the Poisson assumption the distribution remains the same but we have varying `exposure` given by the number of individuals that are represented by one aggregated observation.\n",
    "\n",
    "The parameter estimates and covariance of parameters are the same with the original data, but log-likelihood, deviance and Pearson chi-squared differ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.380077Z",
     "iopub.status.busy": "2022-11-02T17:04:57.379003Z",
     "iopub.status.idle": "2022-11-02T17:04:57.400559Z",
     "shell.execute_reply": "2022-11-02T17:04:57.399969Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:            affairs_sum   No. Observations:                  130\n",
      "Model:                            GLM   Df Residuals:                      126\n",
      "Model Family:                 Poisson   Df Model:                            3\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -740.75\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       967.46\n",
      "Time:                        17:04:57   Pearson chi2:                     926.\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):              1.000\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.7155      0.107     25.294      0.000       2.505       2.926\n",
      "rate_marriage    -0.4952      0.012    -41.702      0.000      -0.518      -0.472\n",
      "age              -0.0299      0.004     -6.691      0.000      -0.039      -0.021\n",
      "yrs_married      -0.0108      0.004     -2.507      0.012      -0.019      -0.002\n",
      "=================================================================================\n"
     ]
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs_sum ~ rate_marriage + age + yrs_married\",\n",
    "    data=df_a,\n",
    "    family=sm.families.Poisson(),\n",
    "    exposure=np.asarray(df_a[\"affairs_count\"]),\n",
    ")\n",
    "res_e = glm.fit()\n",
    "print(res_e.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.404872Z",
     "iopub.status.busy": "2022-11-02T17:04:57.403791Z",
     "iopub.status.idle": "2022-11-02T17:04:57.411459Z",
     "shell.execute_reply": "2022-11-02T17:04:57.410881Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.35078910917956"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_e.pearson_chi2 / res_e.df_resid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### using var_weights\n",
    "\n",
    "We can also use the mean of all combined values of the dependent variable. In this case the variance will be related to the inverse of the total exposure reflected by one combined observation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.415814Z",
     "iopub.status.busy": "2022-11-02T17:04:57.414740Z",
     "iopub.status.idle": "2022-11-02T17:04:57.446055Z",
     "shell.execute_reply": "2022-11-02T17:04:57.445452Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:           affairs_mean   No. Observations:                  130\n",
      "Model:                            GLM   Df Residuals:                      126\n",
      "Model Family:                 Poisson   Df Model:                            3\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -5954.2\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       967.46\n",
      "Time:                        17:04:57   Pearson chi2:                     926.\n",
      "No. Iterations:                     5   Pseudo R-squ. (CS):              1.000\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.7155      0.107     25.294      0.000       2.505       2.926\n",
      "rate_marriage    -0.4952      0.012    -41.702      0.000      -0.518      -0.472\n",
      "age              -0.0299      0.004     -6.691      0.000      -0.039      -0.021\n",
      "yrs_married      -0.0108      0.004     -2.507      0.012      -0.019      -0.002\n",
      "=================================================================================\n"
     ]
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs_mean ~ rate_marriage + age + yrs_married\",\n",
    "    data=df_a,\n",
    "    family=sm.families.Poisson(),\n",
    "    var_weights=np.asarray(df_a[\"affairs_count\"]),\n",
    ")\n",
    "res_a = glm.fit()\n",
    "print(res_a.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison\n",
    "\n",
    "We saw in the summary prints above that ``params`` and ``cov_params`` with associated Wald inference agree across versions. We summarize this in the following comparing individual results attributes across versions.\n",
    "\n",
    "Parameter estimates `params`, standard errors of the parameters `bse` and `pvalues` of the parameters for the tests that the parameters are zeros all agree. However, the likelihood and goodness-of-fit statistics, `llf`, `deviance` and `pearson_chi2` only partially agree. Specifically, the aggregated version do not agree with the results using the original data.\n",
    "\n",
    "**Warning**: The behavior of `llf`, `deviance` and `pearson_chi2` might still change in future versions.\n",
    "\n",
    "Both the sum and average of the response variable for unique values of the explanatory variables have a proper likelihood interpretation. However, this interpretation is not reflected in these three statistics. Computationally this might be due to missing adjustments when aggregated data is used. However, theoretically we can think in these cases, especially for `var_weights` of the misspecified case when likelihood analysis is inappropriate and the results should be interpreted as quasi-likelihood estimates. There is an ambiguity in the definition of ``var_weights`` because they can be used for averages with correctly specified likelihood as well as for variance adjustments in the quasi-likelihood case. We are currently not trying to match the likelihood specification. However, in the next section we show that likelihood ratio type tests still produce the same result for all aggregation versions when we assume that the underlying model is correctly specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.450857Z",
     "iopub.status.busy": "2022-11-02T17:04:57.449298Z",
     "iopub.status.idle": "2022-11-02T17:04:57.455506Z",
     "shell.execute_reply": "2022-11-02T17:04:57.454929Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "results_all = [res_o, res_f, res_e, res_a]\n",
    "names = \"res_o res_f res_e res_a\".split()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.460377Z",
     "iopub.status.busy": "2022-11-02T17:04:57.459287Z",
     "iopub.status.idle": "2022-11-02T17:04:57.479737Z",
     "shell.execute_reply": "2022-11-02T17:04:57.479152Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>res_o</th>\n",
       "      <th>res_f</th>\n",
       "      <th>res_e</th>\n",
       "      <th>res_a</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>2.715533</td>\n",
       "      <td>2.715533</td>\n",
       "      <td>2.715533</td>\n",
       "      <td>2.715533</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rate_marriage</th>\n",
       "      <td>-0.495180</td>\n",
       "      <td>-0.495180</td>\n",
       "      <td>-0.495180</td>\n",
       "      <td>-0.495180</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>-0.029914</td>\n",
       "      <td>-0.029914</td>\n",
       "      <td>-0.029914</td>\n",
       "      <td>-0.029914</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>yrs_married</th>\n",
       "      <td>-0.010763</td>\n",
       "      <td>-0.010763</td>\n",
       "      <td>-0.010763</td>\n",
       "      <td>-0.010763</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  res_o     res_f     res_e     res_a\n",
       "Intercept      2.715533  2.715533  2.715533  2.715533\n",
       "rate_marriage -0.495180 -0.495180 -0.495180 -0.495180\n",
       "age           -0.029914 -0.029914 -0.029914 -0.029914\n",
       "yrs_married   -0.010763 -0.010763 -0.010763 -0.010763"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.concat([r.params for r in results_all], axis=1, keys=names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.484370Z",
     "iopub.status.busy": "2022-11-02T17:04:57.482850Z",
     "iopub.status.idle": "2022-11-02T17:04:57.502196Z",
     "shell.execute_reply": "2022-11-02T17:04:57.501583Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>res_o</th>\n",
       "      <th>res_f</th>\n",
       "      <th>res_e</th>\n",
       "      <th>res_a</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>0.107360</td>\n",
       "      <td>0.107360</td>\n",
       "      <td>0.107360</td>\n",
       "      <td>0.107360</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rate_marriage</th>\n",
       "      <td>0.011874</td>\n",
       "      <td>0.011874</td>\n",
       "      <td>0.011874</td>\n",
       "      <td>0.011874</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>0.004471</td>\n",
       "      <td>0.004471</td>\n",
       "      <td>0.004471</td>\n",
       "      <td>0.004471</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>yrs_married</th>\n",
       "      <td>0.004294</td>\n",
       "      <td>0.004294</td>\n",
       "      <td>0.004294</td>\n",
       "      <td>0.004294</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  res_o     res_f     res_e     res_a\n",
       "Intercept      0.107360  0.107360  0.107360  0.107360\n",
       "rate_marriage  0.011874  0.011874  0.011874  0.011874\n",
       "age            0.004471  0.004471  0.004471  0.004471\n",
       "yrs_married    0.004294  0.004294  0.004294  0.004294"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.concat([r.bse for r in results_all], axis=1, keys=names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.507389Z",
     "iopub.status.busy": "2022-11-02T17:04:57.506297Z",
     "iopub.status.idle": "2022-11-02T17:04:57.526339Z",
     "shell.execute_reply": "2022-11-02T17:04:57.525742Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>res_o</th>\n",
       "      <th>res_f</th>\n",
       "      <th>res_e</th>\n",
       "      <th>res_a</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>3.756282e-141</td>\n",
       "      <td>3.756280e-141</td>\n",
       "      <td>3.756282e-141</td>\n",
       "      <td>3.756282e-141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rate_marriage</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>2.221918e-11</td>\n",
       "      <td>2.221918e-11</td>\n",
       "      <td>2.221918e-11</td>\n",
       "      <td>2.221918e-11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>yrs_married</th>\n",
       "      <td>1.219200e-02</td>\n",
       "      <td>1.219200e-02</td>\n",
       "      <td>1.219200e-02</td>\n",
       "      <td>1.219200e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       res_o          res_f          res_e          res_a\n",
       "Intercept      3.756282e-141  3.756280e-141  3.756282e-141  3.756282e-141\n",
       "rate_marriage   0.000000e+00   0.000000e+00   0.000000e+00   0.000000e+00\n",
       "age             2.221918e-11   2.221918e-11   2.221918e-11   2.221918e-11\n",
       "yrs_married     1.219200e-02   1.219200e-02   1.219200e-02   1.219200e-02"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.concat([r.pvalues for r in results_all], axis=1, keys=names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.531147Z",
     "iopub.status.busy": "2022-11-02T17:04:57.530057Z",
     "iopub.status.idle": "2022-11-02T17:04:57.541068Z",
     "shell.execute_reply": "2022-11-02T17:04:57.540531Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>res_o</th>\n",
       "      <th>res_f</th>\n",
       "      <th>res_e</th>\n",
       "      <th>res_a</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>llf</th>\n",
       "      <td>-10350.913296</td>\n",
       "      <td>-10350.913296</td>\n",
       "      <td>-740.748534</td>\n",
       "      <td>-5954.219866</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>deviance</th>\n",
       "      <td>15374.679054</td>\n",
       "      <td>15374.679054</td>\n",
       "      <td>967.455734</td>\n",
       "      <td>967.455734</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pearson chi2</th>\n",
       "      <td>32310.704118</td>\n",
       "      <td>32310.704118</td>\n",
       "      <td>926.199428</td>\n",
       "      <td>926.199428</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     res_o         res_f       res_e        res_a\n",
       "llf          -10350.913296 -10350.913296 -740.748534 -5954.219866\n",
       "deviance      15374.679054  15374.679054  967.455734   967.455734\n",
       "pearson chi2  32310.704118  32310.704118  926.199428   926.199428"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(\n",
    "    np.column_stack([[r.llf, r.deviance, r.pearson_chi2] for r in results_all]),\n",
    "    columns=names,\n",
    "    index=[\"llf\", \"deviance\", \"pearson chi2\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood Ratio type tests\n",
    "\n",
    "We saw above that likelihood and related statistics do not agree between the aggregated and original, individual data. We illustrate in the following that likelihood ratio test and difference in deviance agree across versions, however Pearson chi-squared does not.\n",
    "\n",
    "As before: This is not sufficiently clear yet and could change.\n",
    "\n",
    "As a test case we drop the `age` variable and compute the likelihood ratio type statistics as difference between reduced or constrained and full or unconstrained model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### original observations and frequency weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.545479Z",
     "iopub.status.busy": "2022-11-02T17:04:57.544399Z",
     "iopub.status.idle": "2022-11-02T17:04:57.582389Z",
     "shell.execute_reply": "2022-11-02T17:04:57.581759Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52.913431618813775, 45.726693322507344, -22.863346661253672)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs ~ rate_marriage + yrs_married\", data=data, family=sm.families.Poisson()\n",
    ")\n",
    "res_o2 = glm.fit()\n",
    "# print(res_f2.summary())\n",
    "res_o2.pearson_chi2 - res_o.pearson_chi2, res_o2.deviance - res_o.deviance, res_o2.llf - res_o.llf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.587543Z",
     "iopub.status.busy": "2022-11-02T17:04:57.586007Z",
     "iopub.status.idle": "2022-11-02T17:04:57.618365Z",
     "shell.execute_reply": "2022-11-02T17:04:57.617547Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52.91343161867189, 45.726693322505525, -22.863346661251853)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs ~ rate_marriage + yrs_married\",\n",
    "    data=dc,\n",
    "    family=sm.families.Poisson(),\n",
    "    freq_weights=np.asarray(dc[\"freq\"]),\n",
    ")\n",
    "res_f2 = glm.fit()\n",
    "# print(res_f2.summary())\n",
    "res_f2.pearson_chi2 - res_f.pearson_chi2, res_f2.deviance - res_f.deviance, res_f2.llf - res_f.llf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### aggregated data: ``exposure`` and ``var_weights``\n",
    "\n",
    "Note: LR test agrees with original observations, ``pearson_chi2`` differs and has the wrong sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.624923Z",
     "iopub.status.busy": "2022-11-02T17:04:57.622892Z",
     "iopub.status.idle": "2022-11-02T17:04:57.648925Z",
     "shell.execute_reply": "2022-11-02T17:04:57.648259Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-31.61852752510731, 45.72669332250655, -22.863346661253217)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs_sum ~ rate_marriage + yrs_married\",\n",
    "    data=df_a,\n",
    "    family=sm.families.Poisson(),\n",
    "    exposure=np.asarray(df_a[\"affairs_count\"]),\n",
    ")\n",
    "res_e2 = glm.fit()\n",
    "res_e2.pearson_chi2 - res_e.pearson_chi2, res_e2.deviance - res_e.deviance, res_e2.llf - res_e.llf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.652024Z",
     "iopub.status.busy": "2022-11-02T17:04:57.651722Z",
     "iopub.status.idle": "2022-11-02T17:04:57.674290Z",
     "shell.execute_reply": "2022-11-02T17:04:57.673618Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-31.618527525103218, 45.72669332250655, -22.863346661253672)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glm = smf.glm(\n",
    "    \"affairs_mean ~ rate_marriage + yrs_married\",\n",
    "    data=df_a,\n",
    "    family=sm.families.Poisson(),\n",
    "    var_weights=np.asarray(df_a[\"affairs_count\"]),\n",
    ")\n",
    "res_a2 = glm.fit()\n",
    "res_a2.pearson_chi2 - res_a.pearson_chi2, res_a2.deviance - res_a.deviance, res_a2.llf - res_a.llf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Investigating Pearson chi-square statistic\n",
    "\n",
    "First, we do some sanity checks that there are no basic bugs in the computation of `pearson_chi2` and `resid_pearson`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.677276Z",
     "iopub.status.busy": "2022-11-02T17:04:57.676987Z",
     "iopub.status.idle": "2022-11-02T17:04:57.684466Z",
     "shell.execute_reply": "2022-11-02T17:04:57.683821Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(894.5809002315173, 926.1994277566246, 894.5809002315173, 926.1994277566246)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_e2.pearson_chi2, res_e.pearson_chi2, (res_e2.resid_pearson ** 2).sum(), (\n",
    "    res_e.resid_pearson ** 2\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.687463Z",
     "iopub.status.busy": "2022-11-02T17:04:57.687167Z",
     "iopub.status.idle": "2022-11-02T17:04:57.694418Z",
     "shell.execute_reply": "2022-11-02T17:04:57.693713Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.617851152192228e-14,\n",
       " array([ 5.42753476, 46.42940306, 19.98971769, 38.50138978, 11.18341883]),\n",
       " array([ 5.42753476, 46.42940306, 19.98971769, 38.50138978, 11.18341883]))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_e._results.resid_response.mean(), res_e.model.family.variance(res_e.mu)[\n",
    "    :5\n",
    "], res_e.mu[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.697392Z",
     "iopub.status.busy": "2022-11-02T17:04:57.697016Z",
     "iopub.status.idle": "2022-11-02T17:04:57.704017Z",
     "shell.execute_reply": "2022-11-02T17:04:57.703349Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "926.1994277566246"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(res_e._results.resid_response ** 2 / res_e.model.family.variance(res_e.mu)).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.707008Z",
     "iopub.status.busy": "2022-11-02T17:04:57.706598Z",
     "iopub.status.idle": "2022-11-02T17:04:57.713826Z",
     "shell.execute_reply": "2022-11-02T17:04:57.713150Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9.9257354472334e-14,\n",
       " array([ 4.77165474, 44.4026604 , 22.2013302 , 39.14749309, 10.54229538]),\n",
       " array([ 4.77165474, 44.4026604 , 22.2013302 , 39.14749309, 10.54229538]))"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_e2._results.resid_response.mean(), res_e2.model.family.variance(res_e2.mu)[\n",
    "    :5\n",
    "], res_e2.mu[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.716726Z",
     "iopub.status.busy": "2022-11-02T17:04:57.716418Z",
     "iopub.status.idle": "2022-11-02T17:04:57.722657Z",
     "shell.execute_reply": "2022-11-02T17:04:57.721340Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "894.5809002315173"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(res_e2._results.resid_response ** 2 / res_e2.model.family.variance(res_e2.mu)).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.725624Z",
     "iopub.status.busy": "2022-11-02T17:04:57.725259Z",
     "iopub.status.idle": "2022-11-02T17:04:57.729891Z",
     "shell.execute_reply": "2022-11-02T17:04:57.729216Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(51204.8573783234, 47104.647795959645)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(res_e2._results.resid_response ** 2).sum(), (res_e._results.resid_response ** 2).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One possible reason for the incorrect sign is that we are subtracting quadratic terms that are divided by different denominators. In some related cases, the recommendation in the literature is to use a common denominator. We can compare pearson chi-squared statistic using the same variance assumption in the full and reduced model. \n",
    "\n",
    "In this case we obtain the same pearson chi2 scaled difference between reduced and full model across all versions. (Issue [#3616](https://github.com/statsmodels/statsmodels/issues/3616) is intended to track this further.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.733538Z",
     "iopub.status.busy": "2022-11-02T17:04:57.733172Z",
     "iopub.status.idle": "2022-11-02T17:04:57.738022Z",
     "shell.execute_reply": "2022-11-02T17:04:57.737331Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.4331417512193"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    (res_e2._results.resid_response ** 2 - res_e._results.resid_response ** 2)\n",
    "    / res_e2.model.family.variance(res_e2.mu)\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.741321Z",
     "iopub.status.busy": "2022-11-02T17:04:57.740957Z",
     "iopub.status.idle": "2022-11-02T17:04:57.745882Z",
     "shell.execute_reply": "2022-11-02T17:04:57.745197Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.43314175121901"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    (res_a2._results.resid_response ** 2 - res_a._results.resid_response ** 2)\n",
    "    / res_a2.model.family.variance(res_a2.mu)\n",
    "    * res_a2.model.var_weights\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.749527Z",
     "iopub.status.busy": "2022-11-02T17:04:57.748839Z",
     "iopub.status.idle": "2022-11-02T17:04:57.753819Z",
     "shell.execute_reply": "2022-11-02T17:04:57.753285Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.43314175121954"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    (res_f2._results.resid_response ** 2 - res_f._results.resid_response ** 2)\n",
    "    / res_f2.model.family.variance(res_f2.mu)\n",
    "    * res_f2.model.freq_weights\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.757347Z",
     "iopub.status.busy": "2022-11-02T17:04:57.756988Z",
     "iopub.status.idle": "2022-11-02T17:04:57.761943Z",
     "shell.execute_reply": "2022-11-02T17:04:57.761404Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.43314175121988"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    (res_o2._results.resid_response ** 2 - res_o._results.resid_response ** 2)\n",
    "    / res_o2.model.family.variance(res_o2.mu)\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remainder\n",
    "\n",
    "The remainder of the notebook just contains some additional checks and can be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.765566Z",
     "iopub.status.busy": "2022-11-02T17:04:57.765196Z",
     "iopub.status.idle": "2022-11-02T17:04:57.770236Z",
     "shell.execute_reply": "2022-11-02T17:04:57.769694Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1., 10.,  5., 10.,  3.]), array([ 1, 10,  5, 10,  3]))"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(res_e2.model.exposure)[:5], np.asarray(df_a[\"affairs_count\"])[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.773002Z",
     "iopub.status.busy": "2022-11-02T17:04:57.772645Z",
     "iopub.status.idle": "2022-11-02T17:04:57.777536Z",
     "shell.execute_reply": "2022-11-02T17:04:57.777011Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.664817945861138"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_e2.resid_pearson.sum() - res_e.resid_pearson.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.780846Z",
     "iopub.status.busy": "2022-11-02T17:04:57.780488Z",
     "iopub.status.idle": "2022-11-02T17:04:57.784860Z",
     "shell.execute_reply": "2022-11-02T17:04:57.784336Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.77165474, 44.4026604 , 22.2013302 , 39.14749309, 10.54229538])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_e2.mu[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.790037Z",
     "iopub.status.busy": "2022-11-02T17:04:57.789659Z",
     "iopub.status.idle": "2022-11-02T17:04:57.794724Z",
     "shell.execute_reply": "2022-11-02T17:04:57.794195Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(894.5809002315159, 926.1994277566191, -42.347207135187446, -32.68238918932903)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_a2.pearson_chi2, res_a.pearson_chi2, res_a2.resid_pearson.sum(), res_a.resid_pearson.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.797429Z",
     "iopub.status.busy": "2022-11-02T17:04:57.797075Z",
     "iopub.status.idle": "2022-11-02T17:04:57.802076Z",
     "shell.execute_reply": "2022-11-02T17:04:57.801535Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "894.5809002315159"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    (res_a2._results.resid_response ** 2)\n",
    "    / res_a2.model.family.variance(res_a2.mu)\n",
    "    * res_a2.model.var_weights\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.805322Z",
     "iopub.status.busy": "2022-11-02T17:04:57.804957Z",
     "iopub.status.idle": "2022-11-02T17:04:57.810574Z",
     "shell.execute_reply": "2022-11-02T17:04:57.810050Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "926.1994277566191"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    (res_a._results.resid_response ** 2)\n",
    "    / res_a.model.family.variance(res_a.mu)\n",
    "    * res_a.model.var_weights\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.813190Z",
     "iopub.status.busy": "2022-11-02T17:04:57.812827Z",
     "iopub.status.idle": "2022-11-02T17:04:57.819051Z",
     "shell.execute_reply": "2022-11-02T17:04:57.818522Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "850.1477584802968"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    (res_a._results.resid_response ** 2)\n",
    "    / res_a.model.family.variance(res_a2.mu)\n",
    "    * res_a.model.var_weights\n",
    ").sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.821708Z",
     "iopub.status.busy": "2022-11-02T17:04:57.821349Z",
     "iopub.status.idle": "2022-11-02T17:04:57.827748Z",
     "shell.execute_reply": "2022-11-02T17:04:57.827208Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0., 39., 17.,  9.,  4.]), array([ 0., 39., 17.,  9.,  4.]))"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_e.model.endog[:5], res_e2.model.endog[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.830477Z",
     "iopub.status.busy": "2022-11-02T17:04:57.830110Z",
     "iopub.status.idle": "2022-11-02T17:04:57.836405Z",
     "shell.execute_reply": "2022-11-02T17:04:57.835876Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.        , 3.9       , 3.4       , 0.9       , 1.33333333]),\n",
       " array([0.        , 3.9       , 3.4       , 0.9       , 1.33333333]))"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_a.model.endog[:5], res_a2.model.endog[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.839098Z",
     "iopub.status.busy": "2022-11-02T17:04:57.838740Z",
     "iopub.status.idle": "2022-11-02T17:04:57.844810Z",
     "shell.execute_reply": "2022-11-02T17:04:57.844287Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0., 39., 17.,  9.,  4.])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_a2.model.endog[:5] * np.exp(res_e2.model.exposure)[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.847621Z",
     "iopub.status.busy": "2022-11-02T17:04:57.847201Z",
     "iopub.status.idle": "2022-11-02T17:04:57.852991Z",
     "shell.execute_reply": "2022-11-02T17:04:57.852468Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0., 39., 17.,  9.,  4.])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_a2.model.endog[:5] * res_a2.model.var_weights[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.855795Z",
     "iopub.status.busy": "2022-11-02T17:04:57.855374Z",
     "iopub.status.idle": "2022-11-02T17:04:57.861720Z",
     "shell.execute_reply": "2022-11-02T17:04:57.861195Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.8390448369994542e-07, 6.931421143170174e-08)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "stats.chi2.sf(27.19530754604785, 1), stats.chi2.sf(29.083798806764687, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.864523Z",
     "iopub.status.busy": "2022-11-02T17:04:57.864102Z",
     "iopub.status.idle": "2022-11-02T17:04:57.871772Z",
     "shell.execute_reply": "2022-11-02T17:04:57.871246Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept        3.756282e-141\n",
       "rate_marriage     0.000000e+00\n",
       "age               2.221918e-11\n",
       "yrs_married       1.219200e-02\n",
       "dtype: float64"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_o.pvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.874610Z",
     "iopub.status.busy": "2022-11-02T17:04:57.874177Z",
     "iopub.status.idle": "2022-11-02T17:04:57.899462Z",
     "shell.execute_reply": "2022-11-02T17:04:57.898789Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:            affairs_sum   No. Observations:                  130\n",
      "Model:                            GLM   Df Residuals:                      127\n",
      "Model Family:                 Poisson   Df Model:                            2\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -763.61\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       1013.2\n",
      "Time:                        17:04:57   Pearson chi2:                     895.\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):              1.000\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.0754      0.050     41.512      0.000       1.977       2.173\n",
      "rate_marriage    -0.4947      0.012    -41.743      0.000      -0.518      -0.471\n",
      "yrs_married      -0.0360      0.002    -17.542      0.000      -0.040      -0.032\n",
      "=================================================================================\n",
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:            affairs_sum   No. Observations:                  130\n",
      "Model:                            GLM   Df Residuals:                      126\n",
      "Model Family:                 Poisson   Df Model:                            3\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -740.75\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       967.46\n",
      "Time:                        17:04:57   Pearson chi2:                     926.\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):              1.000\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.7155      0.107     25.294      0.000       2.505       2.926\n",
      "rate_marriage    -0.4952      0.012    -41.702      0.000      -0.518      -0.472\n",
      "age              -0.0299      0.004     -6.691      0.000      -0.039      -0.021\n",
      "yrs_married      -0.0108      0.004     -2.507      0.012      -0.019      -0.002\n",
      "=================================================================================\n"
     ]
    }
   ],
   "source": [
    "print(res_e2.summary())\n",
    "print(res_e.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:04:57.903757Z",
     "iopub.status.busy": "2022-11-02T17:04:57.903210Z",
     "iopub.status.idle": "2022-11-02T17:04:57.925149Z",
     "shell.execute_reply": "2022-11-02T17:04:57.924377Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                affairs   No. Observations:                  476\n",
      "Model:                            GLM   Df Residuals:                     6363\n",
      "Model Family:                 Poisson   Df Model:                            2\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -10374.\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       15420.\n",
      "Time:                        17:04:57   Pearson chi2:                 3.24e+04\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):             0.9729\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.0754      0.050     41.512      0.000       1.977       2.173\n",
      "rate_marriage    -0.4947      0.012    -41.743      0.000      -0.518      -0.471\n",
      "yrs_married      -0.0360      0.002    -17.542      0.000      -0.040      -0.032\n",
      "=================================================================================\n",
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                affairs   No. Observations:                  476\n",
      "Model:                            GLM   Df Residuals:                     6362\n",
      "Model Family:                 Poisson   Df Model:                            3\n",
      "Link Function:                    Log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -10351.\n",
      "Date:                Wed, 02 Nov 2022   Deviance:                       15375.\n",
      "Time:                        17:04:57   Pearson chi2:                 3.23e+04\n",
      "No. Iterations:                     6   Pseudo R-squ. (CS):             0.9754\n",
      "Covariance Type:            nonrobust                                         \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "Intercept         2.7155      0.107     25.294      0.000       2.505       2.926\n",
      "rate_marriage    -0.4952      0.012    -41.702      0.000      -0.518      -0.472\n",
      "age              -0.0299      0.004     -6.691      0.000      -0.039      -0.021\n",
      "yrs_married      -0.0108      0.004     -2.507      0.012      -0.019      -0.002\n",
      "=================================================================================\n"
     ]
    }
   ],
   "source": [
    "print(res_f2.summary())\n",
    "print(res_f.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}