example/19-hypothesis-first-t.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example of mean hypothesis testing using the *Expected Encompassing Intrinsic Bayes Factor* (EEIBF) method described in [1].\n",
    "\n",
    "The data set comes from a clinical trial by paramacologists Cushny and Peebles ([2]) at the Michigan Asylum for Insane at Kalamazoo (now Kalamazoo Psychiatric Hospital).\n",
    "\n",
    "In the study ([3]), Cushny and Peebles dosed patients with small amount of different but related drugs and measured sleep activity to determine if any of the drugs could be an effective at inducing sleep. The data set was also studied by Student ([4]) and Fisher ([5]).\n",
    "\n",
    "In this notebook, we'll use the same model that Student and Fisher used, but apply EEIBF hypothesis testing.\n",
    "\n",
    "\n",
    "### References\n",
    "1: Berger, J. and J. Mortera (1999). Default bayes factors for nonnested hypothesis testing. Journal of the American Statistical Association 94 (446), 542–554.\n",
    "postscript for paper: https://www2.stat.duke.edu/~berger/papers/mortera.ps\n",
    "\n",
    "2: Cushny, A. R. and Peebles, A. R. ‘The action of optical isomers. 11: Hyoscines’, Journal of Physiology, 32,\n",
    "501-510 (1905).\n",
    "\n",
    "3: Senn S, Richardson W. The first t-test. Stat Med. 1994 Apr 30;13(8):785-803. doi: 10.1002/sim.4780130802. PMID: 8047737.\n",
    "pdf: http://www.medicine.mcgill.ca/epidemiology/hanley/bios601/Mean-Quantile/First.t.test.pdf\n",
    "\n",
    "4: Student. ‘The probable error of a mean’, Biometrika VI (1908); reprinted in Pearson E. S. and Wishart J.\n",
    "‘Student’s Collected Papers, The University Press, Cambridge, 1958 (citation uses reprint page numbers).\n",
    "\n",
    "5: Fisher R. A. Statistical Methods for Research Workers, Oliver and Boyd, Edinburgh, 1925."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy\n",
    "from IPython.display import display, Math\n",
    "from bbai.stat import NormalMeanHypothesis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>Avg Sleep w/ L-hyoscyamine HBr</th>\n",
       "      <th>Avg Sleep w/ L-hyoscine HBr</th>\n",
       "      <th>Delta</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Patient</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.3</td>\n",
       "      <td>2.5</td>\n",
       "      <td>-1.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.4</td>\n",
       "      <td>3.8</td>\n",
       "      <td>-2.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.5</td>\n",
       "      <td>5.8</td>\n",
       "      <td>-1.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.3</td>\n",
       "      <td>5.6</td>\n",
       "      <td>-1.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6.1</td>\n",
       "      <td>6.1</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6.6</td>\n",
       "      <td>7.6</td>\n",
       "      <td>-1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>6.2</td>\n",
       "      <td>8.0</td>\n",
       "      <td>-1.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>3.6</td>\n",
       "      <td>4.4</td>\n",
       "      <td>-0.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.1</td>\n",
       "      <td>5.7</td>\n",
       "      <td>-4.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>4.9</td>\n",
       "      <td>6.3</td>\n",
       "      <td>-1.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>6.3</td>\n",
       "      <td>6.8</td>\n",
       "      <td>-0.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Avg Sleep w/ L-hyoscyamine HBr  Avg Sleep w/ L-hyoscine HBr  Delta\n",
       "Patient                                                                    \n",
       "1                                   1.3                          2.5   -1.2\n",
       "2                                   1.4                          3.8   -2.4\n",
       "3                                   4.5                          5.8   -1.3\n",
       "4                                   4.3                          5.6   -1.3\n",
       "5                                   6.1                          6.1    0.0\n",
       "6                                   6.6                          7.6   -1.0\n",
       "7                                   6.2                          8.0   -1.8\n",
       "8                                   3.6                          4.4   -0.8\n",
       "9                                   1.1                          5.7   -4.6\n",
       "10                                  4.9                          6.3   -1.4\n",
       "11                                  6.3                          6.8   -0.5"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# See page 791 of [3] for a copy of the original data.\n",
    "sleep_a = np.array([1.3, 1.4, 4.5, 4.3, 6.1, 6.6, 6.2, 3.6, 1.1, 4.9, 6.3 ])\n",
    "sleep_b = np.array([2.5, 3.8, 5.8, 5.6, 6.1, 7.6, 8.0, 4.4, 5.7, 6.3, 6.8 ])\n",
    "delta = sleep_a - sleep_b\n",
    "n = len(delta)\n",
    "data = pd.DataFrame({\n",
    "    'Patient': range(1, n+1),\n",
    "    'Avg Sleep w/ L-hyoscyamine HBr' : sleep_a,\n",
    "    'Avg Sleep w/ L-hyoscine HBr' : sleep_b,\n",
    "    'Delta': delta\n",
    "}).set_index('Patient')\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>Avg Sleep w/ L-hyoscyamine HBr</th>\n",
       "      <th>Avg Sleep w/ L-hyoscine HBr</th>\n",
       "      <th>Delta</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>11.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>11.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>4.209091</td>\n",
       "      <td>5.690909</td>\n",
       "      <td>-1.481818</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>2.111613</td>\n",
       "      <td>1.618304</td>\n",
       "      <td>1.211460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>1.100000</td>\n",
       "      <td>2.500000</td>\n",
       "      <td>-4.600000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>2.500000</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>-1.600000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>4.500000</td>\n",
       "      <td>5.800000</td>\n",
       "      <td>-1.300000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>6.150000</td>\n",
       "      <td>6.550000</td>\n",
       "      <td>-0.900000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>6.600000</td>\n",
       "      <td>8.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Avg Sleep w/ L-hyoscyamine HBr  Avg Sleep w/ L-hyoscine HBr      Delta\n",
       "count                       11.000000                    11.000000  11.000000\n",
       "mean                         4.209091                     5.690909  -1.481818\n",
       "std                          2.111613                     1.618304   1.211460\n",
       "min                          1.100000                     2.500000  -4.600000\n",
       "25%                          2.500000                     5.000000  -1.600000\n",
       "50%                          4.500000                     5.800000  -1.300000\n",
       "75%                          6.150000                     6.550000  -0.900000\n",
       "max                          6.600000                     8.000000   0.000000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Hypothesis\n",
    "\n",
    "Test the three hypotheses regarding the mean of delta\n",
    "$$\n",
    "\\begin{align*}\n",
    "H_0&: \\mu = 0, \\quad\\quad\\textrm{L-hyoscyamine HBr and L-hyscine HBr have the same effect}\\\\\n",
    "H_-&: \\mu < 0, \\quad\\quad\\textrm{L-hyoscyamine HBr is less effective at inducing sleep than L-hyscine HBr}\\\\\n",
    "H_+&: \\mu > 0, \\quad\\quad\\textrm{L-hyoscyamine HBr is more effective at inducing sleep than L-hyscine HBr}\n",
    "\\end{align*}\n",
    "$$\n",
    "We'll show how the posterior probabilities for the hypotheses evolve as more data is observed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>delta</th>\n",
       "      <th>H_0</th>\n",
       "      <th>H_-</th>\n",
       "      <th>H_+</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>n</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-1.2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-2.4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-1.3</td>\n",
       "      <td>0.336883</td>\n",
       "      <td>0.472277</td>\n",
       "      <td>0.190840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.3</td>\n",
       "      <td>0.193446</td>\n",
       "      <td>0.734020</td>\n",
       "      <td>0.072534</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.213953</td>\n",
       "      <td>0.704815</td>\n",
       "      <td>0.081232</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.132711</td>\n",
       "      <td>0.827604</td>\n",
       "      <td>0.039685</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>-1.8</td>\n",
       "      <td>0.061951</td>\n",
       "      <td>0.923300</td>\n",
       "      <td>0.014749</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>-0.8</td>\n",
       "      <td>0.034256</td>\n",
       "      <td>0.958680</td>\n",
       "      <td>0.007064</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>-4.6</td>\n",
       "      <td>0.073049</td>\n",
       "      <td>0.911635</td>\n",
       "      <td>0.015316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>-1.4</td>\n",
       "      <td>0.041497</td>\n",
       "      <td>0.950811</td>\n",
       "      <td>0.007692</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.034546</td>\n",
       "      <td>0.959591</td>\n",
       "      <td>0.005863</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    delta       H_0       H_-       H_+\n",
       "n                                      \n",
       "1    -1.2       NaN       NaN       NaN\n",
       "2    -2.4       NaN       NaN       NaN\n",
       "3    -1.3  0.336883  0.472277  0.190840\n",
       "4    -1.3  0.193446  0.734020  0.072534\n",
       "5     0.0  0.213953  0.704815  0.081232\n",
       "6    -1.0  0.132711  0.827604  0.039685\n",
       "7    -1.8  0.061951  0.923300  0.014749\n",
       "8    -0.8  0.034256  0.958680  0.007064\n",
       "9    -4.6  0.073049  0.911635  0.015316\n",
       "10   -1.4  0.041497  0.950811  0.007692\n",
       "11   -0.5  0.034546  0.959591  0.005863"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = len(sleep_a)\n",
    "probs_H_null = []\n",
    "probs_H_minus = []\n",
    "probs_H_plus = []\n",
    "for i in range(n):\n",
    "    if i < 2:\n",
    "        probs_H_null.append(None)\n",
    "        probs_H_minus.append(None)\n",
    "        probs_H_plus.append(None)\n",
    "        continue\n",
    "    probs = NormalMeanHypothesis().test(delta[:i+1])\n",
    "    probs_H_null.append(probs.equal)\n",
    "    probs_H_minus.append(probs.left)\n",
    "    probs_H_plus.append(probs.right)\n",
    "pd.DataFrame({\n",
    "    'n': range(1, n+1),\n",
    "    'delta': delta,\n",
    "    'H_0': probs_H_null,\n",
    "    'H_-': probs_H_minus,\n",
    "    'H_+': probs_H_plus\n",
    "}).set_index('n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}