statistics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model summaries and statistical tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import itertools\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from statsmodels.formula.api import logit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_fit_results(model, subject=None):\n",
    "    \"\"\"Read a json file with cached fit results from a known location.\"\"\"\n",
    "    fpath = f\"fits/{model}.json\" if subject is None else f\"fits/{model}_{subject}.json\"\n",
    "    with open(fpath, \"r\") as fid:\n",
    "        return json.load(fid)\n",
    "\n",
    "\n",
    "def get_param_and_ci(res, key):\n",
    "    \"\"\"Extract parameter value and compute CI from bootstrap distribution.\"\"\"\n",
    "    param = res[\"params\"][key]\n",
    "    boots = np.fromiter((boot[\"params\"][key] for boot in res[\"bootstraps\"]), np.float)\n",
    "    ci = np.percentile(boots, [2.5, 97.5])\n",
    "    return param, ci"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "trial_data = pd.read_csv(\"data/trial_data.csv\")\n",
    "pulse_data = pd.read_csv(\"data/pulse_data.csv\")\n",
    "subjects = [\"S1\", \"S2\", \"S3\", \"S4\", \"S5\"]\n",
    "models = [\"linear\", \"extrema\", \"counting\", \"leaky\"]\n",
    "trial_grouper = [\"subject\", \"timing\", \"session\", \"run\", \"trial\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option(\"display.precision\", 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Basic dataset statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14,869 Total trials\n"
     ]
    }
   ],
   "source": [
    "print(f\"{len(trial_data):,d} Total trials\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6,500 shorter trials\n",
      "8,369 longer trials\n"
     ]
    }
   ],
   "source": [
    "shorter_trials = trial_data.query(\"timing == 'shorter'\")\n",
    "print(f\"{len(shorter_trials):,d} shorter trials\")\n",
    "longer_trials = trial_data.query(\"timing == 'longer'\")\n",
    "print(f\"{len(longer_trials):,d} longer trials\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S1: 3,059 trials\n",
      "S2: 2,882 trials\n",
      "S3: 2,913 trials\n",
      "S4: 2,987 trials\n",
      "S5: 3,028 trials\n"
     ]
    }
   ],
   "source": [
    "for subj, subj_data in trial_data.groupby(\"subject\"):\n",
    "    print(f\"{subj}: {len(subj_data):,d} trials\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Trial duration:\n",
      "  mean+/-std: 10.1+/-5.59 s\n",
      "  range: 2.2–34.0 s\n",
      "\n"
     ]
    }
   ],
   "source": [
    "timing = trial_data.trial_dur.describe()\n",
    "print(f\"\"\"\n",
    "Trial duration:\n",
    "  mean+/-std: {timing['mean']:.1f}+/-{timing['std']:.2f} s\n",
    "  range: {timing['min']:.1f}–{timing['max']:.1f} s\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {}\n",
    "for model, subj in itertools.product(models, subjects):\n",
    "    fit = load_fit_results(model, subj)\n",
    "    for param in fit[\"params\"]:\n",
    "        val, ci = get_param_and_ci(fit, param)\n",
    "        params[(model, param, subj)] = [val] + list(ci)\n",
    "params = pd.DataFrame(params, index=[\"MLE\", \"CI (low)\", \"CI (high)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>MLE</th>\n",
       "      <th>CI (low)</th>\n",
       "      <th>CI (high)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">counting</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">sigma</th>\n",
       "      <th>S1</th>\n",
       "      <td>0.323</td>\n",
       "      <td>0.298</td>\n",
       "      <td>0.347</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>0.348</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.377</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>0.379</td>\n",
       "      <td>0.349</td>\n",
       "      <td>0.411</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>0.447</td>\n",
       "      <td>0.412</td>\n",
       "      <td>0.484</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>0.324</td>\n",
       "      <td>0.298</td>\n",
       "      <td>0.348</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"10\" valign=\"top\">extrema</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">sigma</th>\n",
       "      <th>S1</th>\n",
       "      <td>0.395</td>\n",
       "      <td>0.363</td>\n",
       "      <td>0.426</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>0.428</td>\n",
       "      <td>0.395</td>\n",
       "      <td>0.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>0.421</td>\n",
       "      <td>0.383</td>\n",
       "      <td>0.458</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>0.497</td>\n",
       "      <td>0.457</td>\n",
       "      <td>0.535</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>0.406</td>\n",
       "      <td>0.374</td>\n",
       "      <td>0.436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">theta</th>\n",
       "      <th>S1</th>\n",
       "      <td>0.487</td>\n",
       "      <td>0.451</td>\n",
       "      <td>0.521</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>0.503</td>\n",
       "      <td>0.466</td>\n",
       "      <td>0.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>0.503</td>\n",
       "      <td>0.466</td>\n",
       "      <td>0.541</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>0.513</td>\n",
       "      <td>0.47</td>\n",
       "      <td>0.556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>0.527</td>\n",
       "      <td>0.493</td>\n",
       "      <td>0.559</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"15\" valign=\"top\">leaky</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">alpha</th>\n",
       "      <th>S1</th>\n",
       "      <td>1.38e-15</td>\n",
       "      <td>6.99e-22</td>\n",
       "      <td>0.00902</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>8.2e-15</td>\n",
       "      <td>9.72e-20</td>\n",
       "      <td>0.0259</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>8.58e-16</td>\n",
       "      <td>2.48e-20</td>\n",
       "      <td>0.0114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>1.43e-16</td>\n",
       "      <td>3.24e-20</td>\n",
       "      <td>5.72e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>0.0391</td>\n",
       "      <td>0.00828</td>\n",
       "      <td>0.0702</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">sigma_eps</th>\n",
       "      <th>S1</th>\n",
       "      <td>4.78e-13</td>\n",
       "      <td>1.08e-29</td>\n",
       "      <td>4.43e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>2.36e-17</td>\n",
       "      <td>5.28e-32</td>\n",
       "      <td>1.39e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>3.71e-10</td>\n",
       "      <td>2.99e-30</td>\n",
       "      <td>0.0801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>8.1e-08</td>\n",
       "      <td>1.23e-19</td>\n",
       "      <td>0.16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>2.43e-08</td>\n",
       "      <td>6.13e-09</td>\n",
       "      <td>0.0492</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">sigma_eta</th>\n",
       "      <th>S1</th>\n",
       "      <td>0.424</td>\n",
       "      <td>0.394</td>\n",
       "      <td>0.453</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>0.457</td>\n",
       "      <td>0.423</td>\n",
       "      <td>0.491</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>0.476</td>\n",
       "      <td>0.436</td>\n",
       "      <td>0.508</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>0.568</td>\n",
       "      <td>0.465</td>\n",
       "      <td>0.605</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>0.408</td>\n",
       "      <td>0.374</td>\n",
       "      <td>0.434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">linear</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">sigma</th>\n",
       "      <th>S1</th>\n",
       "      <td>0.424</td>\n",
       "      <td>0.395</td>\n",
       "      <td>0.453</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>0.457</td>\n",
       "      <td>0.424</td>\n",
       "      <td>0.489</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>0.476</td>\n",
       "      <td>0.442</td>\n",
       "      <td>0.51</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>0.568</td>\n",
       "      <td>0.523</td>\n",
       "      <td>0.615</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>0.412</td>\n",
       "      <td>0.383</td>\n",
       "      <td>0.438</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            MLE  CI (low) CI (high)\n",
       "counting sigma     S1     0.323     0.298     0.347\n",
       "                   S2     0.348      0.32     0.377\n",
       "                   S3     0.379     0.349     0.411\n",
       "                   S4     0.447     0.412     0.484\n",
       "                   S5     0.324     0.298     0.348\n",
       "extrema  sigma     S1     0.395     0.363     0.426\n",
       "                   S2     0.428     0.395      0.46\n",
       "                   S3     0.421     0.383     0.458\n",
       "                   S4     0.497     0.457     0.535\n",
       "                   S5     0.406     0.374     0.436\n",
       "         theta     S1     0.487     0.451     0.521\n",
       "                   S2     0.503     0.466      0.54\n",
       "                   S3     0.503     0.466     0.541\n",
       "                   S4     0.513      0.47     0.556\n",
       "                   S5     0.527     0.493     0.559\n",
       "leaky    alpha     S1  1.38e-15  6.99e-22   0.00902\n",
       "                   S2   8.2e-15  9.72e-20    0.0259\n",
       "                   S3  8.58e-16  2.48e-20    0.0114\n",
       "                   S4  1.43e-16  3.24e-20  5.72e-10\n",
       "                   S5    0.0391   0.00828    0.0702\n",
       "         sigma_eps S1  4.78e-13  1.08e-29  4.43e-05\n",
       "                   S2  2.36e-17  5.28e-32  1.39e-05\n",
       "                   S3  3.71e-10  2.99e-30    0.0801\n",
       "                   S4   8.1e-08  1.23e-19      0.16\n",
       "                   S5  2.43e-08  6.13e-09    0.0492\n",
       "         sigma_eta S1     0.424     0.394     0.453\n",
       "                   S2     0.457     0.423     0.491\n",
       "                   S3     0.476     0.436     0.508\n",
       "                   S4     0.568     0.465     0.605\n",
       "                   S5     0.408     0.374     0.434\n",
       "linear   sigma     S1     0.424     0.395     0.453\n",
       "                   S2     0.457     0.424     0.489\n",
       "                   S3     0.476     0.442      0.51\n",
       "                   S4     0.568     0.523     0.615\n",
       "                   S5     0.412     0.383     0.438"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params.astype(object).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log-likelihoods and model comparisons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "loglikes = {}\n",
    "for model in models:\n",
    "    agg, agg_cv = 0, 0\n",
    "    for subj in subjects:\n",
    "        fit = load_fit_results(model, subj)\n",
    "        loglikes[(model, subj)] = [fit[\"loglike\"], fit[\"cv_loglike\"]]\n",
    "        agg += fit[\"loglike\"]\n",
    "        agg_cv += fit[\"cv_loglike\"]\n",
    "    loglikes[(model, \"agg\")] = [agg, agg_cv]\n",
    "loglikes = pd.DataFrame(loglikes, index=[\"loglike\", \"cv_loglike\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>loglike</th>\n",
       "      <th>cv_loglike</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">linear</th>\n",
       "      <th>S1</th>\n",
       "      <td>-925.781</td>\n",
       "      <td>-927.538</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>-920.391</td>\n",
       "      <td>-922.830</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>-954.744</td>\n",
       "      <td>-955.920</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>-1123.667</td>\n",
       "      <td>-1125.917</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>-880.606</td>\n",
       "      <td>-882.946</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>agg</th>\n",
       "      <td>-4805.189</td>\n",
       "      <td>-4815.151</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">extrema</th>\n",
       "      <th>S1</th>\n",
       "      <td>-1099.216</td>\n",
       "      <td>-1101.284</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>-1077.112</td>\n",
       "      <td>-1079.987</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>-1085.826</td>\n",
       "      <td>-1089.015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>-1218.818</td>\n",
       "      <td>-1222.534</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>-1108.981</td>\n",
       "      <td>-1111.011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>agg</th>\n",
       "      <td>-5589.953</td>\n",
       "      <td>-5603.831</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">counting</th>\n",
       "      <th>S1</th>\n",
       "      <td>-940.396</td>\n",
       "      <td>-942.150</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>-919.868</td>\n",
       "      <td>-922.141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>-977.324</td>\n",
       "      <td>-979.239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>-1129.088</td>\n",
       "      <td>-1132.006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>-908.208</td>\n",
       "      <td>-911.157</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>agg</th>\n",
       "      <td>-4874.884</td>\n",
       "      <td>-4886.694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">leaky</th>\n",
       "      <th>S1</th>\n",
       "      <td>-925.781</td>\n",
       "      <td>-927.538</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S2</th>\n",
       "      <td>-920.391</td>\n",
       "      <td>-922.954</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S3</th>\n",
       "      <td>-954.744</td>\n",
       "      <td>-955.920</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S4</th>\n",
       "      <td>-1123.667</td>\n",
       "      <td>-1125.928</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S5</th>\n",
       "      <td>-875.545</td>\n",
       "      <td>-881.457</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>agg</th>\n",
       "      <td>-4800.128</td>\n",
       "      <td>-4813.797</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               loglike  cv_loglike\n",
       "linear   S1   -925.781    -927.538\n",
       "         S2   -920.391    -922.830\n",
       "         S3   -954.744    -955.920\n",
       "         S4  -1123.667   -1125.917\n",
       "         S5   -880.606    -882.946\n",
       "         agg -4805.189   -4815.151\n",
       "extrema  S1  -1099.216   -1101.284\n",
       "         S2  -1077.112   -1079.987\n",
       "         S3  -1085.826   -1089.015\n",
       "         S4  -1218.818   -1222.534\n",
       "         S5  -1108.981   -1111.011\n",
       "         agg -5589.953   -5603.831\n",
       "counting S1   -940.396    -942.150\n",
       "         S2   -919.868    -922.141\n",
       "         S3   -977.324    -979.239\n",
       "         S4  -1129.088   -1132.006\n",
       "         S5   -908.208    -911.157\n",
       "         agg -4874.884   -4886.694\n",
       "leaky    S1   -925.781    -927.538\n",
       "         S2   -920.391    -922.954\n",
       "         S3   -954.744    -955.920\n",
       "         S4  -1123.667   -1125.928\n",
       "         S5   -875.545    -881.457\n",
       "         agg -4800.128   -4813.797"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loglikes[models].T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Integration vs. Extrema Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>S1</th>\n",
       "      <th>S2</th>\n",
       "      <th>S3</th>\n",
       "      <th>S4</th>\n",
       "      <th>S5</th>\n",
       "      <th>agg</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>loglike</th>\n",
       "      <td>173.435</td>\n",
       "      <td>156.721</td>\n",
       "      <td>131.081</td>\n",
       "      <td>95.151</td>\n",
       "      <td>228.375</td>\n",
       "      <td>784.764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>cv_loglike</th>\n",
       "      <td>173.746</td>\n",
       "      <td>157.158</td>\n",
       "      <td>133.095</td>\n",
       "      <td>96.617</td>\n",
       "      <td>228.065</td>\n",
       "      <td>788.680</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 S1       S2       S3      S4       S5      agg\n",
       "loglike     173.435  156.721  131.081  95.151  228.375  784.764\n",
       "cv_loglike  173.746  157.158  133.095  96.617  228.065  788.680"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loglikes[\"linear\"] - loglikes[\"extrema\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Integration vs. Counting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>S1</th>\n",
       "      <th>S2</th>\n",
       "      <th>S3</th>\n",
       "      <th>S4</th>\n",
       "      <th>S5</th>\n",
       "      <th>agg</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>loglike</th>\n",
       "      <td>14.615</td>\n",
       "      <td>-0.523</td>\n",
       "      <td>22.579</td>\n",
       "      <td>5.421</td>\n",
       "      <td>27.602</td>\n",
       "      <td>69.695</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>cv_loglike</th>\n",
       "      <td>14.611</td>\n",
       "      <td>-0.688</td>\n",
       "      <td>23.320</td>\n",
       "      <td>6.089</td>\n",
       "      <td>28.211</td>\n",
       "      <td>71.543</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                S1     S2      S3     S4      S5     agg\n",
       "loglike     14.615 -0.523  22.579  5.421  27.602  69.695\n",
       "cv_loglike  14.611 -0.688  23.320  6.089  28.211  71.543"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loglikes[\"linear\"] - loglikes[\"counting\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Integration vs. Leaky Integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>S1</th>\n",
       "      <th>S2</th>\n",
       "      <th>S3</th>\n",
       "      <th>S4</th>\n",
       "      <th>S5</th>\n",
       "      <th>agg</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>loglike</th>\n",
       "      <td>-1.003e-07</td>\n",
       "      <td>-1.027e-07</td>\n",
       "      <td>-7.947e-08</td>\n",
       "      <td>-2.936e-07</td>\n",
       "      <td>-5.061</td>\n",
       "      <td>-5.061</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>cv_loglike</th>\n",
       "      <td>-2.156e-04</td>\n",
       "      <td>1.244e-01</td>\n",
       "      <td>1.519e-04</td>\n",
       "      <td>1.093e-02</td>\n",
       "      <td>-1.489</td>\n",
       "      <td>-1.354</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   S1         S2         S3         S4     S5    agg\n",
       "loglike    -1.003e-07 -1.027e-07 -7.947e-08 -2.936e-07 -5.061 -5.061\n",
       "cv_loglike -2.156e-04  1.244e-01  1.519e-04  1.093e-02 -1.489 -1.354"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loglikes[\"linear\"] - loglikes[\"leaky\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likeratio_test(test, null, dof, cv=False):\n",
    "\n",
    "    col = \"cv_loglike\" if cv else \"loglike\"\n",
    "    stat = 2 * (loglikes.loc[col, test] - loglikes.loc[col, null])\n",
    "    pval = stats.chi2(dof).sf(stat)\n",
    "    return pd.DataFrame(np.c_[stat, pval],\n",
    "                        columns=[\"stat\", \"pval\"],\n",
    "                        index=loglikes.loc[\"loglike\", test].index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>S1</th>\n",
       "      <th>S2</th>\n",
       "      <th>S3</th>\n",
       "      <th>S4</th>\n",
       "      <th>S5</th>\n",
       "      <th>agg</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>stat</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>10.122</td>\n",
       "      <td>10.122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pval</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.006</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       S1   S2   S3   S4      S5     agg\n",
       "stat  0.0  0.0  0.0  0.0  10.122  10.122\n",
       "pval  1.0  1.0  1.0  1.0   0.006   0.006"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "likeratio_test(\"leaky\", \"linear\", dof=2).round(3).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>S1</th>\n",
       "      <th>S2</th>\n",
       "      <th>S3</th>\n",
       "      <th>S4</th>\n",
       "      <th>S5</th>\n",
       "      <th>agg</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>stat</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.249</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>-0.022</td>\n",
       "      <td>2.978</td>\n",
       "      <td>2.708</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pval</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.000</td>\n",
       "      <td>0.226</td>\n",
       "      <td>0.258</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       S1     S2   S3     S4     S5    agg\n",
       "stat  0.0 -0.249 -0.0 -0.022  2.978  2.708\n",
       "pval  1.0  1.000  1.0  1.000  0.226  0.258"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "likeratio_test(\"leaky\", \"linear\", dof=2, cv=True).round(3).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Statistical test of accuracy improvements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_odd_vs_even(trial_data):\n",
    "    \"\"\"Use logistic regression to test \"\"\"\n",
    "    trials = trial_data.pulse_count < 5\n",
    "    tsub = trial_data.loc[trials, [\"pulse_count\", \"correct\"]]\n",
    "    y = tsub[\"correct\"]\n",
    "    df = pd.DataFrame(dict(\n",
    "        correct=tsub[\"correct\"],\n",
    "        parity=tsub[\"pulse_count\"] % 2 == 0,\n",
    "        magnitude=tsub[\"pulse_count\"] > 2,\n",
    "    )).astype(int)\n",
    "    m = logit(\"correct ~ magnitude + parity\", df).fit(disp=False)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>      <td>correct</td>     <th>  No. Observations:  </th>  <td> 13427</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>               <td>Logit</td>      <th>  Df Residuals:      </th>  <td> 13424</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>MLE</td>       <th>  Df Model:          </th>  <td>     2</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Wed, 14 Nov 2018</td> <th>  Pseudo R-squ.:     </th> <td>0.009477</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>10:34:42</td>     <th>  Log-Likelihood:    </th> <td> -7198.2</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>           <td>True</td>       <th>  LL-Null:           </th> <td> -7267.1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                      <td> </td>        <th>  LLR p-value:       </th> <td>1.232e-30</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    0.9589</td> <td>    0.029</td> <td>   32.867</td> <td> 0.000</td> <td>    0.902</td> <td>    1.016</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>magnitude</th> <td>    0.4483</td> <td>    0.046</td> <td>    9.842</td> <td> 0.000</td> <td>    0.359</td> <td>    0.538</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>parity</th>    <td>    0.2517</td> <td>    0.042</td> <td>    5.949</td> <td> 0.000</td> <td>    0.169</td> <td>    0.335</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                correct   No. Observations:                13427\n",
       "Model:                          Logit   Df Residuals:                    13424\n",
       "Method:                           MLE   Df Model:                            2\n",
       "Date:                Wed, 14 Nov 2018   Pseudo R-squ.:                0.009477\n",
       "Time:                        10:34:42   Log-Likelihood:                -7198.2\n",
       "converged:                       True   LL-Null:                       -7267.1\n",
       "                                        LLR p-value:                 1.232e-30\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      0.9589      0.029     32.867      0.000       0.902       1.016\n",
       "magnitude      0.4483      0.046      9.842      0.000       0.359       0.538\n",
       "parity         0.2517      0.042      5.949      0.000       0.169       0.335\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = test_odd_vs_even(trial_data)\n",
    "m.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P = 2.7e-09\n"
     ]
    }
   ],
   "source": [
    "print(f\"P = {m.pvalues.parity:.2g}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S1: P = 0.0023\n",
      "S2: P = 0.0069\n",
      "S3: P = 0.0088\n",
      "S4: P = 0.0031\n",
      "S5: P = 0.0011\n"
     ]
    }
   ],
   "source": [
    "for subj in subjects:\n",
    "    tdata = trial_data.query(\"subject == @subj\")\n",
    "    m = test_odd_vs_even(tdata)\n",
    "    print(f\"{subj}: P = {m.pvalues.parity / 2:.2g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "## Statistical test of differences between shorter and longer gaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_longer_vs_shorter(pulse_data):\n",
    "    pulse_data = pulse_data.merge(trial_data, on=trial_grouper)\n",
    "    fit_data = pulse_data.groupby(trial_grouper).mean().reset_index()\n",
    "    m = logit(\"response ~ pulse_llr * timing\", fit_data).fit(disp=False)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>     <td>response</td>     <th>  No. Observations:  </th>  <td> 14869</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>               <td>Logit</td>      <th>  Df Residuals:      </th>  <td> 14865</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>MLE</td>       <th>  Df Model:          </th>  <td>     3</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Wed, 14 Nov 2018</td> <th>  Pseudo R-squ.:     </th>  <td>0.5393</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>10:34:42</td>     <th>  Log-Likelihood:    </th> <td> -4747.5</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>           <td>True</td>       <th>  LL-Null:           </th> <td> -10304.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                      <td> </td>        <th>  LLR p-value:       </th>  <td> 0.000</td> \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "               <td></td>                  <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>                   <td>    0.1367</td> <td>    0.034</td> <td>    3.981</td> <td> 0.000</td> <td>    0.069</td> <td>    0.204</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>timing[T.shorter]</th>           <td>   -0.1552</td> <td>    0.053</td> <td>   -2.950</td> <td> 0.003</td> <td>   -0.258</td> <td>   -0.052</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>pulse_llr</th>                   <td>    5.4972</td> <td>    0.118</td> <td>   46.579</td> <td> 0.000</td> <td>    5.266</td> <td>    5.728</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>pulse_llr:timing[T.shorter]</th> <td>    0.2910</td> <td>    0.184</td> <td>    1.582</td> <td> 0.114</td> <td>   -0.070</td> <td>    0.652</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:               response   No. Observations:                14869\n",
       "Model:                          Logit   Df Residuals:                    14865\n",
       "Method:                           MLE   Df Model:                            3\n",
       "Date:                Wed, 14 Nov 2018   Pseudo R-squ.:                  0.5393\n",
       "Time:                        10:34:42   Log-Likelihood:                -4747.5\n",
       "converged:                       True   LL-Null:                       -10304.\n",
       "                                        LLR p-value:                     0.000\n",
       "===============================================================================================\n",
       "                                  coef    std err          z      P>|z|      [0.025      0.975]\n",
       "-----------------------------------------------------------------------------------------------\n",
       "Intercept                       0.1367      0.034      3.981      0.000       0.069       0.204\n",
       "timing[T.shorter]              -0.1552      0.053     -2.950      0.003      -0.258      -0.052\n",
       "pulse_llr                       5.4972      0.118     46.579      0.000       5.266       5.728\n",
       "pulse_llr:timing[T.shorter]     0.2910      0.184      1.582      0.114      -0.070       0.652\n",
       "===============================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = test_longer_vs_shorter(pulse_data)\n",
    "m.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S1: P = 0.16\n",
      "S2: P = 0.91\n",
      "S3: P = 0.42\n",
      "S4: P = 0.03\n",
      "S5: P = 0.52\n"
     ]
    }
   ],
   "source": [
    "for subj in subjects:\n",
    "    pdata = pulse_data.query(\"subject == @subj\")\n",
    "    m = test_longer_vs_shorter(pdata)\n",
    "    print(f\"{subj}: P = {m.pvalues['pulse_llr:timing[T.shorter]']:.2g}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pulses-psych (py36)",
   "language": "python",
   "name": "pulses-psych"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}