week07c_generalized_linear_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The generalized linear model ###\n",
    "\n",
    "#### Reference: \n",
    "\n",
    "http://statsthinking21.org\n",
    "\n",
    "* The general linear model\n",
    "* Comparing means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Common statistical tests as linear models ##\n",
    "\n",
    "Many statistical tests can be though of as implementations of the generalized linear model. Thinking of tests as part of a class of linear models can be more intuitive than thinking about how each test works individually.\n",
    "\n",
    "### Resources\n",
    "\n",
    "This approach is taken in Chapter 28 of *Statistical Thinking for the 21st Century* on comparing means:\n",
    "* http://statsthinking21.org/comparing-means.html\n",
    "\n",
    "A blog post by Jonas Kristoffer Lindeløv explains this approach for a wide array of statistical tests. Implementation of the statistical functions and linear models, with interpretations, are provided in both R and Python.\n",
    "* Original post (using R): https://lindeloev.github.io/tests-as-linear/\n",
    "* Python port: https://eigenfoo.xyz/tests-as-linear/\n",
    "\n",
    "### Examples\n",
    "\n",
    "The following examples show different ways of comparing means, using data from the 2007 West Coast Ocean Acidification cruise. The examples use quality controlled data from 0-10 dbar (upper 10m of the water column)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing two sample means (another example)\n",
    "\n",
    "The following shows the same calculations for temperature. In this case, the null hypothesis can be rejected at the 95% confidence level, and the 95% confidence intervals for the model slope do *not* overlap with zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load 2007 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import stats\n",
    "import statsmodels.formula.api as smf\n",
    "import pingouin as pg\n",
    "\n",
    "import PyCO2SYS as pyco2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename07 = 'data/wcoa_cruise_2007/32WC20070511.exc.csv'\n",
    "df07 = pd.read_csv(filename07,header=29,na_values=-999,parse_dates=[[6,7]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the PyCO2SYS package to calculate seawater carbon chemistry parameters.\n",
    "\n",
    "https://pyco2sys.readthedocs.io/en/latest/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/autograd/tracer.py:48: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return f_raw(*args, **kwargs)\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:99: RuntimeWarning: invalid value encountered in sqrt\n",
      "  lnKF = 1590.2 / TempK - 12.641 + 1.525 * IonS ** 0.5\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:577: RuntimeWarning: overflow encountered in power\n",
      "  K1 = 10.0 ** -(pK1)\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:583: RuntimeWarning: overflow encountered in power\n",
      "  K2 = 10.0 ** -(pK2)\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:603: RuntimeWarning: overflow encountered in power\n",
      "  K1 = 10.0 ** -pK1\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:611: RuntimeWarning: overflow encountered in power\n",
      "  K2 = 10.0 ** -pK2\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:636: RuntimeWarning: overflow encountered in power\n",
      "  K1 = 10.0 ** -pK1\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:641: RuntimeWarning: overflow encountered in power\n",
      "  K2 = 10.0 ** -pK2\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:653: RuntimeWarning: overflow encountered in power\n",
      "  K1 = 10.0 ** -pK1\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:658: RuntimeWarning: overflow encountered in power\n",
      "  K2 = 10.0 ** -pK2\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/equilibria/p1atm.py:715: RuntimeWarning: overflow encountered in power\n",
      "  K2 = 10.0 ** -pK2\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/solubility.py:41: RuntimeWarning: overflow encountered in power\n",
      "  KAr = 10.0 ** logKAr  # this is in (mol/kg-SW)^2\n",
      "/Users/tomconnolly/programs/miniconda3/envs/ms263-22/lib/python3.8/site-packages/PyCO2SYS/solubility.py:25: RuntimeWarning: overflow encountered in power\n",
      "  KCa = 10.0 ** logKCa  # this is in (mol/kg-SW)^2 at zero pressure\n"
     ]
    }
   ],
   "source": [
    "c07 = pyco2.sys(df07['ALKALI'], df07['TCARBN'], 1, 2,\n",
    "               salinity=df07['CTDSAL'], temperature=df07['CTDTMP'], \n",
    "                pressure=df07['CTDPRS'])\n",
    "\n",
    "df07['OmegaA'] = c07['saturation_aragonite']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing two sample means (another example)\n",
    "\n",
    "The following shows the same calculations for temperature. In this case, the null hypothesis can be rejected at the 95% confidence level, and the 95% confidence intervals for the model slope do *not* overlap with zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing one sample mean to a single value\n",
    "\n",
    "In this example the goal is to test whether the mean aragonite saturation state is different from a value of 1, a critical threshold for the ability of organisms to form calcium carbonate shells.\n",
    "\n",
    "* $H_0$: $\\bar{\\Omega}_A =$ 1\n",
    "* $H_A$: $\\bar{\\Omega}_A \\neq$ 1 \n",
    "\n",
    "The first step is to create a subset of good data from the upper 10m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "iisurf07 = ((df07['CTDPRS'] <= 10) &\n",
    "      (df07['NITRAT_FLAG_W'] == 2) & (df07['PHSPHT_FLAG_W'] == 2)\n",
    "      & (df07['CTDOXY_FLAG_W'] == 2) & (df07['CTDSAL_FLAG_W'] == 2) \n",
    "        & (df07['ALKALI_FLAG_W'] == 2) & (df07['TCARBN_FLAG_W'] == 2))\n",
    "\n",
    "df07surf = df07[iisurf07]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A box plot is one way of showing the distribution of $\\Omega_A$ values.\n",
    "* Orange line: median, or 50th percentile\n",
    "* Upper/lower limits on box: interquartile range, or 75th/25th percentiles\n",
    "* Whiskers: each have length 1.5*interquartile range (pyplot default)\n",
    "* Circles: extreme values\n",
    "* Green triangle: mean\n",
    "* Notches on box: 95% confidence intervals for median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, '$\\\\Omega_A$ - upper 10m 2007')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEMCAYAAADQ553CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAej0lEQVR4nO3df1TUVf4/8OfwW9HpHHOIEpHNX7FgammlyxHJk5iIKKnBYRUdcsstYdWFiOFomYih2ZqzphmY7BIbJ1KgTT0G6IqanUM/lmmhNQttwAF0VQSF4cf9/uHX+cw4yAwwOM2b5+OcOc373vfc92vq+PR23z9GJoQQICIiyXCydwFERGRbDHYiIolhsBMRSQyDnYhIYhjsREQSw2AnIpIYBjsRkcQw2KnP8vLyMHfuXEycOBEzZszAW2+9hZaWFnuXRTRgyXiDEvXF7t27UVhYiPT0dAQGBqK6uhrr16+Hm5sbsrKyIJPJ7F2izQkhkJycjHHjxiEuLg4AcOzYMbz99tvQ6/UYP348Nm/ejCFDhvT6GAUFBcjMzIRMJsOgQYOgUqkwYcIEi8fqqu+LL77Avn37DGNfv34ddXV1OH78OIYPH96HfxP0qyWIeunKlSti4sSJoqqqyqT98uXL4rHHHhPFxcV2qqz//Pjjj2Lp0qVi4sSJ4oMPPhBC3Pq+Tz31lPj555+FEEJkZGSIDRs29PoY586dE7/73e9EXV2dEEKIY8eOieDgYIvHsqYOvV4vlixZInJzc3tdH/36cSmGeu3rr7/G8OHDMX78eJP2YcOGYeLEiThz5oxNjnPmzBnMmzfPbPvMmTNYvHgxEhISEB4ejsWLF+PcuXOGfe7WBwAlJSVYvHgxFixYgKioKHzzzTc4c+YM5s+fj6ioKISHh0Ov15vVkpOTg8WLF2POnDmGtrKyMkyYMAF+fn4AgOjoaBQVFUEIgTNnzuD555/Hn/70J0RERCAqKgolJSVYsWIFZs6cic2bN5sdw83NDZs2bYKXlxcAIDAwEJcuXYJer+/2WN313bZ3714MGzYMUVFRvfuPQQ6BwU69duPGDQwdOrTLvqFDh+LmzZuG7bNnz8Lf3x86nc6mNWg0GixduhRFRUWIjIxEYmKixb7q6mq88847eP/993Hw4EG8+eabWL16NW7evImzZ8/i7bffRlFREdzc3MyOt379eoSHh5u06XQ6eHt7G7a9vb3R1NSE5uZmAEBFRQX+8Ic/oKCgAEOGDMH777+PPXv24NNPP8VHH32Euro6k/F8fHwwc+ZMALeWfdLT0/H000/Dzc2t22NZquN///sf9u3bh5SUlN78qyYHwmCnXhs5ciTOnz+P1tZWs74ff/wRo0aNMmzv3bsXERERJrNmW3jkkUcwZcoUAMBzzz2HyspKXLlypdu+kydPor6+HsuXL0dERAT+/Oc/QyaT4fz583jwwQcxYsSIHtXQ2dnZ5bkEJ6dbf7x8fHzw29/+FgDg6+uLJ598Em5ubhg2bBg8PT1x7dq1Lse9ceMGEhIScOHCBWzatMnisSzVkZeXh1mzZmHkyJE9+n7keBjs1GuPPvoovL298cknn5i0nz59GjU1NQgLCwMAVFVVQaFQICgoqMtg37FjByIiIhAREYEdO3aY9ctkMpPlhLa2NsN7Z2dns/1vt92tr7OzE9OmTUNBQYHhlZeXh3HjxmHw4MFWfvv/8+CDD6K+vt6wXVdXh/vuu88w1p0zfxcXF4tj1tbWIioqCs7OzsjOzoZcLrd4LEt1fP7554iMjOzx9yPHw2CnXpPJZNiyZQt27tyJo0ePArgV4snJyXjjjTcMywL79+/HypUrMWbMmC6DPSEhwRCwCQkJZv3Dhg1DbW0tLl++DCEE/vnPfxr6qqqqUFVVBQD4+OOPMXnyZEMI3q1v2rRpOHnypKGW48ePY/78+b2+RDMoKAjfffcdqqurAQD/+Mc/MGvWrF6NBQBNTU1YunQpZs+ejXfeeQceHh5WHau7vmvXruHChQuYPHlyr+six2F56kB0F7t378aePXvQ2tqKV155BT/88APS0tJQV1eH119/HUVFRVi3bh2+/vprrFu3Dp2dnejo6OjxccaMGYOoqCg899xzUCgUmDlzJioqKgAAw4cPx1/+8hfU1NRg2LBhyMjIMHzubn1jxozBxo0bsXbtWggh4OLigvfee69XtQHA/fffj/T0dMTHx6OtrQ2+vr546623ejUWcOsEbW1tLY4ePWr4CxMAPvzww26P1V3f+fPnoVAo4Orq2uu6yHHwOnbqV2vWrMHrr7+O++67DwCwaNEis6Wb3jpz5gzefPNNfPbZZz3qI5I6LsVQv9FoNPDw8DCEOgC4u7sbTm4SUf/gjJ2ISGI4YycikhgGOxGRxDDYiYgkxq6XO7a0tECj0UChUHR5MwkREZnr6OhAQ0MDAgMDTe5zuM2uwa7RaBATE2PPEoiIHFZOTo7hsRnG7BrsCoUCwK3ijB9eREREd6fT6RATE2PI0DvZNdhvL794e3vDx8fHnqUQETmcuy1h8+QpEZHEMNiJiCSGwU5EJDEMdiIiiWGwE90hNzcXgYGBcHZ2RmBgIHJzc+1dElGP8HnsREZyc3OhUqmQmZmJoKAglJWVIS4uDsCtH4cmcgScsRMZSUtLQ2ZmJkJCQuDq6oqQkBBkZmYiLS3N3qURWY3BTmSksrISQUFBJm1BQUGorKy0U0VEPcdgJzLi7++PsrIyk7aysjL4+/vbqSKinmOwExlRqVSIi4tDaWkp2traUFpairi4OKhUKnuXRmQ1njwlMnL7BOnq1atRWVkJf39/pKWl8cQpORQGO9EdoqOjGeTk0LgUQ3QHXsdOjo4zdiIjvI6dpIAzdiIjvI6dpIDBTmSE17GTFDDYiYzwOnaSAgY7kRFex05SwJOnREZ4HTtJgVXBvmPHDhw5cgQymQyLFi3CihUrTPrVajXy8/Mhl8sBAEuWLEFMTIztqyW6B3gdOzk6i8H+1Vdf4csvv0RhYSHa29sxd+5cBAcH4+GHHzbso9FosH37dkyePLlfiyUiIsssrrE/8cQTyM7OhouLCy5fvoyOjg4MHjzYZB+NRoM9e/YgPDwcGzduRGtra78VTERE3bPq5KmrqyveffddhIWFYdq0aXjggQcMfc3NzfD390diYiIOHDiAxsZG7Nq1y2yMxsZGaLVak5dOp7PdNyEiIgCATAghrN355s2beOmllzB37lw8//zzXe7zn//8BykpKTh48KBJ+86dO6FWq7v8THFxMXx8fKyvmohoANNqtZg1a9Zds9PiGvu5c+eg1+vh7++PQYMGYfbs2fjhhx8M/bW1tTh16hQWLVoEABBCwMXFfNjY2FgsXLjQpE2n0/EkKxGRjVlcitFqtUhNTYVer4der0dxcTEef/xxQ7+Hhwe2bt2KX375BUII5OTk4JlnnjEbRy6Xw8fHx+Tl7e1t229DZAN8CBg5Oosz9uDgYPz73//GggUL4OzsjNmzZyMsLAwrV65EfHw8JkyYgI0bN2LVqlVoa2vDY489ZnY5JJGjyM3NRUJCAjw9PSGEQHNzMxISEgDwIWDkOHq0xm5rltaJiO61kSNH4vLly2hvb0dbWxtcXV3h4uKC+++/H7/88ou9yyMCYDk7+UgBIiNarRatra3YsmULmpubsWXLFrS2tkKr1dq7NCKrMdiJ7vDCCy9g7dq1GDx4MNauXYsXXnjB3iUR9QiDnegOhYWFJg8BKywstHdJRD3Ch4ARGXFxccH169ehVCpx4cIF+Pr64vr1611ewkv0a8UZO5GRl156CTdv3sTNmzfR2dlpeP/SSy/ZuzQiq3EaQmRk586dAIC9e/cCAK5evYo//vGPhnYiR8AZO9Edpk+fjjFjxsDJyQljxozB9OnT7V0SUY9wxk5kJDc3FyqVCpmZmQgKCkJZWRni4uIA8AYlchycsRMZSUtLQ2ZmJkJCQuDq6oqQkBBkZmYiLS3N3qURWY0zdhowAgMD8f3331vc7+mnn+6yXSaTWfxsQEAANBpNj2sjsiXO2GnA0Gg0EEJ0+woICEBJSQluP2lDCIGSkhIEBARY/KwQgqFOvwoMdiIjKpUKcXFxKC0tBQCUlpYiLi4OKpXKzpURWY9LMURGbp8gXb16teGfaWlpPHFKDoXBTnSH6OhoREdHQyaTcWmFHBKXYoiIJIbBTkQkMQx2IiKJYbATEUkMg52ISGIY7EREEsNgJyKSGAY7EZHEMNiJiCSGwU5EJDEMdiIiibEq2Hfs2IG5c+ciLCwM+/btM+uvrKxEZGQkQkNDoVKp0N7ebvNCiYjIOhaD/auvvsKXX36JwsJC5Ofn429/+xt++uknk30SExOxfv16HDlyBEII5OXl9VvBRETUPYvB/sQTTyA7OxsuLi64fPkyOjo6MHjwYEN/TU0NWlpaMGnSJABAZGQkDh8+bDZOY2MjtFqtyUun09numxAREQArH9vr6uqKd999F1lZWZgzZw4eeOABQ199fT0UCoVhW6FQoK6uzmyM/fv3Q61W26BkIiLqjtUnT+Pj43H69GlcvHjRZKmls7PT5LcghRBd/jZkbGwsiouLTV45OTl9LJ+IiO5kccZ+7tw56PV6+Pv7Y9CgQZg9ezZ++OEHQ7+3tzcaGhoM25cuXYKXl5fZOHK5HHK53EZlExHR3VicsWu1WqSmpkKv10Ov16O4uBiPP/64oX/EiBFwd3dHeXk5AKCgoAAzZszov4qJiKhbFoM9ODgYM2fOxIIFC/Dcc89h8uTJCAsLw8qVK1FRUQEA2LZtG9LT0zFnzhzcuHEDy5Yt6/fCiYioazIhhLDXwbVaLWbNmoXi4mL4+PjYqwyiLslkMtjxjwfRXVnKTt55SkQkMQx2IiKJYbATEUkMg52ISGIY7EREEsNgJyKSGAY7EZHEMNiJiCSGwU5EJDEMdiIiiWGwExFJDIOdiEhiGOxERBLDYCcikhgGOxGRxDDYiYgkhsFORCQxDHYiIolhsBMRSQyDnYhIYhjsREQSw2AnIpIYBjsRkcQw2ImIJIbBTkQkMS7W7KRWq3Ho0CEAQHBwMJKSksz68/PzIZfLAQBLlixBTEyMjUslIiJrWAz2U6dOoaysDAcOHIBMJsMLL7yAo0eP4plnnjHso9FosH37dkyePLlfiyUiIsssBrtCoUBycjLc3NwAAKNHj0Ztba3JPhqNBnv27EFNTQ2mTp2KV199Fe7u7ib7NDY2orGx0aRNp9P1tX4iIrqDxWAfO3as4X11dTUOHTqE3NxcQ1tzczP8/f2RmJiIUaNGITk5Gbt27cKaNWtMxtm/fz/UarUNS6eB6vr161AqldDr9f1+rIiIiH4df/Dgwfjwww/NJkJEfWHVGjsAnD17Fi+++CKSkpLg5+dnaPf09MTevXsN20qlEikpKWbBHhsbi4ULF5q06XQ6rsVTj9XX1+PEiRPYs2dPvx7n22+/hVKp7NdjLFu2DNeuXYOXl1e/HocGFquCvby8HPHx8UhJSUFYWJhJX21tLU6dOoVFixYBAIQQcHExH1YulxtOrhL1laenZ7/Ppvt7fACcqVO/sHi548WLF/Hyyy9j27ZtZqEOAB4eHti6dSt++eUXCCGQk5NjcmKViIjuLYsz9szMTLS2tmLLli2GtqioKJSUlCA+Ph4TJkzAxo0bsWrVKrS1teGxxx7DihUr+rVoIiK6O4vBnpqaitTUVLP26Ohow/vQ0FCEhobatjIiIuoV3nlKRCQxDHYiIolhsBMRSQyDnYhIYhjsREQSw2AnIpIYBjsRkcQw2ImIJIbBTkQkMQx2IiKJYbCTw3F2dsbly5dRWFiItrY2e5fTK3q9Hvn5+WhqaoKzs7O9yyGJYbCTw/Hz88POnTuxdetW+Pr6IikpCZWVlTY9xpWb17ChZDuu3rxm03ErKiqwZs0a+Pj44K9//Ss+/PBD3H///TY9BhGDnRzS0qVLceLECRw7dgxOTk6YNWsWpk+fjg8++MDsJxh7I//7z1HV8CM++f7zPo919epVvPfee5g6dSrmzp2LIUOG4PTp0ygpKcGSJUv6PD7RnRjs5NDGjx+PLVu24MKFC0hJScHnn38OX19fLF++HP/6178ghOjxmFduXkNp9WkICByrPt2rWXtnZyeKi4sRExMDPz8/HDt2DJs2bUJ1dTXefPNNjB49usdjElmLwU6S4OLignnz5uHTTz/Ff//7Xzz66KNYtWoVxo0bh82bN6OmpsbqsfK//xxCdAIAOkVnj2bt58+fxxtvvIGHH34Y69atw1NPPYVz587h448/RmhoKNfT6Z5gsJPkeHl5Ye3atdBoNNi9ezfUajV8fX1x4sQJi5+9PVtv7+wAALR3dlg9az98+DB+85vfYN++ffj73/+Ob7/9FqtXr+YaOt1zVv+YNZGjEEKgrKwMWVlZOHjwIEJCQqBUKuH3n3z89K+/dPvZA4oh6Bw6CHCSGdo62vTIyo7HgktN3X52HICKjOVIPt6A+fPnY9GiRVAqlXjyySchk8m6/SyRLTHYSTJqamqQnZ2NrKwsuLm5IS4uDm+99Ra8vLz+/x7zLI5RdyQNHVe1Jm0dTjLo/Mbh4RdVVtVR9OdbP/KenZ2NZcuWwcXFBUqlEkuXLsUDDzzQ069F1GMMdnJoer0eRUVFyMrKwunTp7FkyRLk5ORg6tSpvZolZ4RaF96WPPTQQ0hOTsarr76KkydPIisrC4888giCg4OhVCrx7LPPwtXV1SbHIroT19jJIRlfD65WqxEVFQWtVovdu3fjiSee+NUsfchkMgQFBSErKwsXLlzA/PnzkZGRgZEjRyIpKQlVVVX2LpEkiMFODken02H69OmG68FLS0uxdOlSDB482N6ldWvo0KFQKpUoKyvD8ePH4eTkhKlTp+LaNdveBEXEYCeH09zcDC8vL4e+Hvz29feDBg1Ca2urvcshiWGwExFJDIOdiEhirAp2tVqNsLAwhIWFISMjw6y/srISkZGRCA0NhUqlQnt7u80LJSIi61gM9lOnTqGsrAwHDhzAwYMH8f333+Po0aMm+yQmJmL9+vU4cuQIhBDIy8vrt4KJiKh7FoNdoVAgOTkZbm5ucHV1xejRo1FbW2vor6mpQUtLCyZNmgQAiIyMxOHDh/utYCIi6p7FG5TGjh1reF9dXY1Dhw4hNzfX0FZfXw+FQmHYVigUqKurMxunsbHR7HGqOp2uV0UTEdHdWX3n6dmzZ/Hiiy8iKSkJfn5+hvbOzk6Tm0GEEF3eHLJ//36o1eq+VUtERBZZFezl5eWIj49HSkoKwsLCTPq8vb3R0NBg2L506ZLRszn+T2xsLBYuXGjSptPpEBMT05u6iYjoLiwG+8WLF/Hyyy/jnXfewbRp08z6R4wYAXd3d5SXl+Pxxx9HQUEBZsyYYbafXC6HXC63TdU0oLW1teGnn37Cs88+a+9S+uzq1atwcuJVx2RbFoM9MzMTra2t2LJli6EtKioKJSUliI+Px4QJE7Bt2zakpqaiqakJAQEBWLZsWb8WTQPb+PHj8cEHH+Chhx7q1+OsWrUK7733Xr8eIzk5GcOHD+/XY9DAIxO9+e0wG9FqtZg1axaKi4vh4+NjrzKIuiSTyXr103pE/c1SdvL/AYmIJIbBTkQkMQx2IiKJYbATEUkMg52ISGIY7EREEsNgJyKSGAY7EZHEMNiJiCSGwU5EJDEMdiIiiWGwExFJDIOdiEhiGOxERBLDYCcikhgGOxGRxDDYiYgkhsFORCQxDHYiIolhsBMRSQyDnYhIYhjsREQSw2AnIpIYBjsRkcQw2ImIJMaqYG9qasK8efOg1WrN+tRqNUJCQhAREYGIiAjk5OTYvEgiIrKei6UdvvvuO6SmpqK6urrLfo1Gg+3bt2Py5Mm2ro2IiHrB4ow9Ly8PGzZsgJeXV5f9Go0Ge/bsQXh4ODZu3IjW1labF0lERNazGOxpaWmYMmVKl33Nzc3w9/dHYmIiDhw4gMbGRuzatavLfRsbG6HVak1eOp2ub9UTEZEZi0sx3fH09MTevXsN20qlEikpKVizZo3Zvvv374dare7L4YiIyAp9Cvba2lqcOnUKixYtAgAIIeDi0vWQsbGxWLhwoUmbTqdDTExMX0ogIqI79CnYPTw8sHXrVjz55JPw8fFBTk4OnnnmmS73lcvlkMvlfTkcERFZoVfXsa9cuRIVFRUYNmwYNm7ciFWrVmHOnDkQQmDFihW2rpGIiHrA6hl7SUmJ4b3xunpoaChCQ0NtWxUREfUa7zwlIpIYBjsRkcQw2ImIJIbBTkQkMQx2IiKJYbATEUkMg52ISGIY7EREEsNgJyKSGAY7EZHEMNiJiCSGwU5EJDEMdiIiiWGwExFJDIOdiEhiGOxERBLDYCcikhgGOw0YgYGBkMlkVr8A9Gh/mUyGwMBAO39LIgY7DSAajQZCCIuvjz76CAEBAXByckJAQAA++ugjqz4nhIBGo7H31yRisBMZy83NRUJCApqbmyGEQHNzMxISEpCbm2vv0oisxmAnMpKUlARnZ2dkZWWhtbUVWVlZcHZ2RlJSkr1LI7Iag53IiFarRXZ2NkJCQuDq6oqQkBBkZ2dDq9XauzQiqzHYiYgkhsFOZMTHxwexsbEoLS1FW1sbSktLERsbCx8fH3uXRmQ1BjuRkYyMDLS3t0OpVMLDwwNKpRLt7e3IyMiwd2lEVrMq2JuamjBv3rwu1xkrKysRGRmJ0NBQqFQqtLe327xIonslOjoaO3bsgKenJwDA09MTO3bsQHR0tJ0rI7KexWD/7rvvEB0djerq6i77ExMTsX79ehw5cgRCCOTl5dm6RqJ7Kjo6GhqNBh0dHdBoNAx1cjgWgz0vLw8bNmyAl5eXWV9NTQ1aWlowadIkAEBkZCQOHz7c5TiNjY3QarUmL51O17fqiYjIjIulHdLS0u7aV19fD4VCYdhWKBSoq6vrct/9+/dDrVb3okQiIuqJPp087ezsNDxTAwCEECbbxmJjY1FcXGzyysnJ6cvhifpFbm4uAgMD4ezsjMDAQN51Sg7H4oy9O97e3mhoaDBsX7p0qcslGwCQy+WQy+V9ORxRv7v9SAFPT0+TRwoA4Fo7OYw+zdhHjBgBd3d3lJeXAwAKCgowY8YMmxRGZA98pABJQa+CfeXKlaioqAAAbNu2Denp6ZgzZw5u3LiBZcuW2bRAonuJjxQgKbB6KaakpMTwfu/evYb3jzzyCD755BPbVkVERL3GO0+JjPCRAiQFDHYiI3ykAEkBg53ICB8pQFLQp8sdiaQoOjqaQU4OjTN2ojvwBiVydJyxExnJzc2FSqVCZmYmgoKCUFZWhri4OAC8QYkcB2fsREbS0tKQmZlpch17ZmZmt89MIvq1YbATGamsrERQUJBJW1BQECorK+1UEVHPMdiJjPj7+6OsrMykraysDP7+/naqiKjnGOxERlQqFeLi4kxuUIqLi4NKpbJ3aURW48lTIiO3T5CuXr0alZWV8Pf3R1paGk+ckkNhsBPdgdexk6PjUgwRkcQw2ImIJIbBTkQkMQx2IiKJsevJ046ODgCATqezZxlERA7ldmbeztA72TXYb/8QdkxMjD3LICJySA0NDRg1apRZu0wIIexQDwCgpaUFGo0GCoUCzs7O9iqDyIxOp0NMTAxycnLg7e1t73KITHR0dKChoQGBgYHw8PAw67frjN3DwwNTpkyxZwlE3fL29ubP4tGvUlcz9dt48pSISGIY7EREEsNgJyKSGAY7URfkcjleeeUVyOVye5dC1GN2vSqGiIhsjzN2IiKJYbATEUkMn8dOA45arcahQ4cAAMHBwUhKSsKpU6eQnp6O1tZWPPvss1izZg2AW7+BqlKp0NzcjClTpuCNN97AtWvXoFQqDeNdv34dV65cwTfffGOX70NkRhANICdPnhTPP/+8aG1tFXq9XixbtkwUFRWJ4OBgceHCBdHW1iaUSqU4duyYEEKIsLAw8c033wghhHjttddETk6OyXgdHR3i97//vSgsLLzXX4XorrgUQwOKQqFAcnIy3Nzc4OrqitGjR6O6uhqjRo3CyJEj4eLigvDwcBw+fBg1NTVoaWnBpEmTAACRkZE4fPiwyXj5+fkYNGgQwsPD7fBtiLrGYKcBZezYsYagrq6uxqFDhyCTyaBQKAz7eHl5oa6uDvX19SbtCoUCdXV1hu2Ojg7s3r0b69atu2f1E1mDwU4D0tmzZ6FUKpGUlISRI0dCJpMZ+oQQkMlk6Ozs7LL9thMnTsDPzw/jx4+/p7UTWcJgpwGnvLwcy5cvx7p167Bw4UJ4e3sbHiEN3HoUqpeXl1n7pUuX4OXlZdj+4osvMHfu3HtaO5E1GOw0oFy8eBEvv/wytm3bhrCwMADAxIkT8fPPP+P8+fPo6OjAZ599hhkzZmDEiBFwd3dHeXk5AKCgoAAzZswwjPXtt9/y6aT0q8Q7T2lA2bRpE/Lz8+Hr62toi4qKgp+fn+Fyx+DgYLz22muQyWSoqqpCamoqmpqaEBAQgPT0dLi5uQG49RfCV199BXd3d3t9HaIuMdiJiCSGSzFERBLDYCcikhgGOxGRxDDYiYgkhsFORCQxDHYiIolhsBMRSQyDnYhIYv4ffx42O0T26MMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.boxplot(np.array(df07['OmegaA'][iisurf07]),labels=['2007'],showmeans=True,notch=True);\n",
    "plt.title('$\\Omega_A$ - upper 10m 2007')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Method 1: one sample t-test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A one-sample t-test can be used to test whether the null hypothesis can be rejected at the 95% confidence level ($\\alpha$ = 0.05). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_1sampResult(statistic=27.742288359723034, pvalue=4.457126416644256e-58)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ttest_1samp(np.array(df07surf['OmegaA']),popmean=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Method 2: generalized linear model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, this test can be framed in terms of a general linear model\n",
    "\n",
    "$$ \\hat{y} = \\hat{a}_1 + \\hat{a}_2 x ,$$\n",
    "\n",
    "In this application, $y$ represents the $\\Omega_A$ data. There is only one group, so we can set $x = 0$ for all values, making the slope parameter $\\hat{a}_2$ irrelevant. The model then reduces to\n",
    "\n",
    "$$ \\hat{y} = \\hat{a}_1 ,$$\n",
    "\n",
    "or a model for the intercept parameter only. This equation can also be expressed as\n",
    "\n",
    "$$ \\hat{y} = \\hat{a}_1 \\times 1.$$\n",
    "\n",
    "This model for the data in terms of a constant intercept can be implemented with the `statsmodels` library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = smf.ols(formula=\"OmegaA ~ 1\", data=df07surf).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The summary of the results shows that the intercept is 2.33, which is also the mean of the data. The 95% confidence intervals do not overlap with 1, which means that the null hypothesis can be rejected at $\\alpha$ = 0.05. This approach to hypothesis testing will give the same results as the one sample t-test for $N$ > 14.\n",
    "\n",
    "Notice that the test statistic $t$ is different from the one sample t-test. This is because `statsmodels` automatically tests whether parameters are different from *zero*, while in this case we are interested in whether the mean/intercept is different from *one*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>OmegaA</td>      <th>  R-squared:         </th> <td>   0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>     nan</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 08 Mar 2022</td> <th>  Prob (F-statistic):</th>  <td>   nan</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>14:09:53</td>     <th>  Log-Likelihood:    </th> <td> -102.08</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   138</td>      <th>  AIC:               </th> <td>   206.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   137</td>      <th>  BIC:               </th> <td>   209.1</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     0</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    2.2017</td> <td>    0.043</td> <td>   50.827</td> <td> 0.000</td> <td>    2.116</td> <td>    2.287</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 7.228</td> <th>  Durbin-Watson:     </th> <td>   0.655</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.027</td> <th>  Jarque-Bera (JB):  </th> <td>  11.596</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.177</td> <th>  Prob(JB):          </th> <td> 0.00303</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.375</td> <th>  Cond. No.          </th> <td>    1.00</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 OmegaA   R-squared:                       0.000\n",
       "Model:                            OLS   Adj. R-squared:                  0.000\n",
       "Method:                 Least Squares   F-statistic:                       nan\n",
       "Date:                Tue, 08 Mar 2022   Prob (F-statistic):                nan\n",
       "Time:                        14:09:53   Log-Likelihood:                -102.08\n",
       "No. Observations:                 138   AIC:                             206.2\n",
       "Df Residuals:                     137   BIC:                             209.1\n",
       "Df Model:                           0                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      2.2017      0.043     50.827      0.000       2.116       2.287\n",
       "==============================================================================\n",
       "Omnibus:                        7.228   Durbin-Watson:                   0.655\n",
       "Prob(Omnibus):                  0.027   Jarque-Bera (JB):               11.596\n",
       "Skew:                           0.177   Prob(JB):                      0.00303\n",
       "Kurtosis:                       4.375   Cond. No.                         1.00\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing two sample means\n",
    "\n",
    "We can also apply the generalized linear model approach when comparing two means. In this case, we will examine whether there is a statistically significant difference between the mean of $\\Omega_A$ to the north and south of Cape Mendocino. At a latitude of 40.4$^o$N, Cape Mendocino represents a sharp transition point in many oceanographic processes and water masses.\n",
    "\n",
    "The first steps are to make two subsets based on latitude, and then visualize the results in a box plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new boolean variable in the df07surf dataframe\n",
    "df07surf = df07[iisurf07]\n",
    "df07surf = df07surf.assign(is_northern = df07surf['LATITUDE'] > 40.4);\n",
    "\n",
    "iinorth = np.array(df07surf.is_northern == True)\n",
    "iisouth = np.array(df07surf.is_northern == False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, '$\\\\Omega_A$ - upper 10m 2007')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEXCAYAAAC59m+aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAl8ElEQVR4nO3de1RU5d4H8O8IcpHEwqApKTxp1iikpJUX8kaJC1S8YIEEJqRZCWYdSBlS01BD09TRvDQaGkyRloqZVmCdgxd607M4To29XsIcCBD1SCo4MOz3D1/mMAIOlxk37Pl+1mIt9n727P0b23x75tnP3iMTBEEAERFJRgexCyAiIutisBMRSQyDnYhIYhjsREQSw2AnIpIYBjsRkcQw2ImIJIbBTq2WmZmJ4OBg9O3bF0OHDsX777+PyspKscsislsy3qBErbFhwwbs2bMHS5cuha+vLwoKCjB//nw4OTlhy5YtkMlkYpdodYIgYO7cuejVqxdiY2MBAD/88AM++OADGAwGPProo1iyZAnuuuuuFh9j9+7dUKvVkMlkcHV1hVKphJ+fn8VjNdT2/fffY+vWraZ9//XXXygpKcGPP/6Ie++9txX/EtRmCUQtdPnyZaFv377CyZMnzdZfvHhReOKJJ4Ts7GyRKrOd06dPC1FRUULfvn2Fjz/+WBCEm+934MCBwu+//y4IgiCkpqYKCxYsaPExzpw5IwwZMkQoKSkRBEEQfvjhB2HYsGEWj9WUOgwGg/D8888LGo2mxfVR28ehGGqx48eP495778Wjjz5qtt7DwwN9+/ZFXl6eVY6Tl5eHMWPG1FvOy8vD5MmTMXv2bIwdOxaTJ0/GmTNnTNs01gYAOTk5mDx5MsaPH4/w8HD861//Ql5eHsaNG4fw8HCMHTsWBoOhXi3p6emYPHkyRo8ebVqXm5sLPz8/dO/eHQAQERGBrKwsCIKAvLw8vPDCC3jjjTcQGhqK8PBw5OTkYNq0aRg+fDiWLFlS7xhOTk5477334OXlBQDw9fVFWVkZDAbDbY91u7ZamzdvhoeHB8LDw1v2H4PaBQY7tdj169fRuXPnBts6d+6MiooK0/KpU6egUChQXFxs1Rq0Wi2ioqKQlZWFiRMnIiEhwWJbQUEBVq1ahU2bNmHXrl1YvHgx4uLiUFFRgVOnTuGDDz5AVlYWnJyc6h1v/vz5GDt2rNm64uJiyOVy07JcLsfVq1dx7do1AMCJEycwY8YM7N69G3fddRc2bdqEjRs34ssvv0RGRgZKSkrM9uft7Y3hw4cDuDnss3TpUowcORJOTk63PZalOi5duoStW7ciKSmpJf/U1I4w2KnFHnzwQZw7dw43btyo13b69Gn4+PiYljdv3ozQ0FCzXrM1PPbYYxgwYAAAYNKkSdDpdLh8+fJt2w4dOoTS0lK89NJLCA0Nxd///nfIZDKcO3cO999/P7p169asGmpqahq8ltChw80/L29vb/Tu3RsA8NBDD+Hpp5+Gk5MTPDw84ObmhitXrjS43+vXr2P27Nn4448/8N5771k8lqU6MjMzERgYiAcffLBZ74/aHwY7tdjjjz8OuVyOHTt2mK0/cuQICgsLERISAgA4efIkPD09ERAQ0GCwr169GqGhoQgNDcXq1avrtctkMrPhhKqqKtPvDg4O9bavXddYW01NDQYNGoTdu3ebfjIzM9GrVy906tSpie/+v+6//36UlpaalktKStClSxfTvm7t+Ts6OlrcZ1FREcLDw+Hg4IBt27bB3d3d4rEs1bFv3z5MnDix2e+P2h8GO7WYTCbDsmXLsHbtWnz33XcAbob43Llz8e6775qGBdLS0jB9+nT07NmzwWCfPXu2KWBnz55dr93DwwNFRUW4ePEiBEHA119/bWo7efIkTp48CQD4/PPP4e/vbwrBxtoGDRqEQ4cOmWr58ccfMW7cuBZP0QwICEB+fj4KCgoAAJ999hkCAwNbtC8AuHr1KqKiojBq1CisWrUKLi4uTTrW7dquXLmCP/74A/7+/i2ui9oPy10HokZs2LABGzduxI0bNzBr1iz89ttvSElJQUlJCRYuXIisrCy89dZbOH78ON566y3U1NTAaDQ2+zg9e/ZEeHg4Jk2aBE9PTwwfPhwnTpwAANx777348MMPUVhYCA8PD6Smpppe11hbz549sWjRIrz55psQBAGOjo746KOPWlQbAHTt2hVLly5FfHw8qqqq8NBDD+H9999v0b6Amxdoi4qK8N1335n+hwkAn3zyyW2Pdbu2c+fOwdPTEx07dmxxXdR+cB472dScOXOwcOFCdOnSBQAQFhZWb+impfLy8rB48WLs3bu3WW1EUsehGLIZrVYLFxcXU6gDgLOzs+niJhHZBnvsREQSwx47EZHEMNiJiCRG1FkxlZWV0Gq18PT0bHDOMRER1Wc0GnHhwgX4+vqaTYetJWqwa7VaREZGilkCEVG7lZ6ebrq7ui5Rg93T0xPAzeLqPuOCiIgaV1xcjMjISFOG3krUYK8dfpHL5fD29hazFCKidqexIWxePCUikhgGOxGRxDDYiYgkhsFORCQxDHaJ0Gg08PX1hYODA3x9faHRaMQuiYhEwsf2SoBGo4FSqYRarUZAQAByc3MRGxsL4Ob3XhKRfWGPXQJSUlKgVqsxYsQIdOzYESNGjIBarUZKSorYpRGRCBjsEqDT6RAQEGC2LiAgADqdTqSKiEhMDHYJUCgUyM3NNVuXm5sLhUIhUkVEJCYGuwQolUrExsbi4MGDqKqqwsGDBxEbGwulUil2aUQkAl48lYDaC6RxcXHQ6XRQKBRISUnhhVMiO8Vgl4iIiAgGOREB4FCMZHAeOxHVYo9dAjiPnYjqYo9dAjiPnYjqYrBLAOexE1FdDHYJ4Dx2IqqLwS4BnMdORHXx4qkEcB47EdXFYJcIzmMnolpNCvbVq1fjwIEDkMlkCAsLw7Rp08zaVSoVdu7cCXd3dwDA888/j8jISOtXS0REFlkM9p9++glHjx7Fnj17UF1djeDgYAwbNgwPP/ywaRutVouVK1fC39/fpsUSEZFlFoP9qaeewrZt2+Do6IiSkhIYjUZ06tTJbButVouNGzeisLAQTz75JN5++204OzubbVNeXo7y8nKzdcXFxVZ4C0REVFeThmI6duyINWvWYMuWLRg9ejTuu+8+U9u1a9egUCiQkJAAHx8fzJ07F+vXr8ecOXPM9pGWlgaVSmXd6omIqB6ZIAhCUzeuqKjAzJkzERwcjBdeeKHBbX799VckJSVh165dZusb67FHRkYiOzsb3t7eza+eiMgO6fV6BAYGNpqdFuexnzlzxnQHo6urK0aNGoXffvvN1F5UVIQdO3aYlgVBgKNj/Q8C7u7u8Pb2NvuRy+UtelNUHx8CRkS1LAa7Xq9HcnIyDAYDDAYDsrOz0b9/f1O7i4sLli9fjvPnz0MQBKSnp+O5556zadFkrvYhYGvXrkVlZSXWrl0LpVLJcKc2gZ2OO89isA8bNgzDhw/H+PHjMWnSJPj7+yMkJATTp0/HiRMn4OHhgUWLFuHVV1/F6NGjIQhCvemQZFspKSmYMmUK4uLi4OLigri4OEyZMoUPASPRaTQazJ49G9euXQNw85rc7NmzGe62Jojo/PnzQq9evYTz58+LWUa7J5PJhO7duws5OTmCwWAQcnJyhO7duwsymUzs0sjOeXt7C/fff7/ZuXn//fcL3t7eYpfWrlnKTj4rRgKcnJwQFxdn9tjeuLg4ODk5iV0a2Tm9Xo+0tDSzczMtLQ16vV7s0iSNwS4BBoMBKpXK7CFgKpUKBoNB7NKISAQMdgno3bt3g2PsvXv3Frs0snPe3t6Ijo4263RER0dzerONMdglQKlUIiMjw2xWTEZGBh/bS6JLTU2F0WhETEwMnJ2dERMTA6PRiNTUVLFLkzQ+3VEC+Nheaqtqz8GUlBTIZDK4ublhyZIlPDdtjMEuEXxsLxHV4lAMEdkM57GLg8FORDaTmJgIR0dHbNmyBZWVldiyZQscHR2RmJgodmmSxmAnIpvhPHZxMNiJiCSGwU5ENsN57OJgsBORzXAeuzgY7ERkMxEREVi9ejXc3NxM89hXr17Nqbk2xnnsRGRTvMfizmOPnYhsil+0ceexx05ENlP77V5qtRoBAQHIzc1FbGwsALAXb0PssRORzaSkpECtVpvNY1er1fx2LxtjsBORzeh0OgQEBJitCwgIgE6nE6ki+8BgJyKbUSgUyM3NNVuXm5sLhUIhUkX2gcFORDajVCoRGxtrdoNSbGwsvyvAxnjxlIhsht8VIA4GOxHZFOex33kciiEikhgGOxGRxDDYiYgkpknBvnr1agQHByMkJARbt26t167T6TBx4kQEBQVBqVSiurra6oUSEVHTWAz2n376CUePHsWePXuwc+dObN++HWfPnjXbJiEhAfPnz8eBAwcgCAIyMzNtVjAREd2exWB/6qmnsG3bNjg6OuLixYswGo3o1KmTqb2wsBCVlZXo168fAGDixInYv39/vf2Ul5dDr9eb/RQXF1vvnRAREYAmTnfs2LEj1qxZgy1btmD06NG47777TG2lpaXw9PQ0LXt6eqKkpKTePtLS0qBSqaxQMhER3U6TL57Gx8fjyJEj+PPPP82GWmpqaiCTyUzLgiCYLdeaOnUqsrOzzX7S09NbWT4REd3KYo/9zJkzMBgMUCgUcHV1xahRo/Dbb7+Z2uVyOS5cuGBaLisrg5eXV739uLu7w93d3UplExFRYyz22PV6PZKTk2EwGGAwGJCdnY3+/fub2rt16wZnZ2ccO3YMALB7924MHTrUdhUTEdFtWQz2YcOGYfjw4Rg/fjwmTZoEf39/hISEYPr06Thx4gQAYMWKFVi6dClGjx6N69evIzo62uaFExFRw2SCIAhiHVyv1yMwMBDZ2dnw9vYWqwwionbFUnbyzlMiIolhsBMRSQyDnYhIYhjsREQSw2Bvh3x9fSGTyZr14+vrK3bZZCeae37y3LQ+foNSO6TVahttk8lkEHGiE1Gj5yfPzTuHPXYiIolhsBMRSQyDnYhIYhjsREQSw2AnIpIYBjsRkcQw2ImIJIbBTkQkMQx2IiKJYbATEUkMg52ISGIY7EREEsNgJyKSGAY7EZHEMNiJiCSGwU5EJDEMdiIiiWGwExFJDIOdiEhimvSdpyqVCt988w0AYNiwYUhMTKzXvnPnTri7uwMAnn/+eURGRlq5VCIiagqLwX748GHk5ubiq6++gkwmw8svv4zvvvsOzz33nGkbrVaLlStXwt/f36bFEhGRZRaD3dPTE3PnzoWTkxMAoEePHigqKjLbRqvVYuPGjSgsLMSTTz6Jt99+G87OzmbblJeXo7y83GxdcXFxa+uXrJ9++gnDhg1DTU1Ns19767/97TzzzDP4/vvvm30MImq7LAb7I488Yvq9oKAA33zzDTQajWndtWvXoFAokJCQAB8fH8ydOxfr16/HnDlzzPaTlpYGlUplxdKl7dKlS3jmmWeQlZXVrNc98cQTOH78eJO2/eWXXxAdHd2S8oioDWvSGDsAnDp1Cq+88goSExPRvXt303o3Nzds3rzZtBwTE4OkpKR6wT516lRMmDDBbF1xcTHH4m+jQ4cOzep9AzfDuqlqP4URkbQ0KdiPHTuG+Ph4JCUlISQkxKytqKgIhw8fRlhYGABAEAQ4Otbfrbu7u+niKhER2Y7F6Y5//vknXn/9daxYsaJeqAOAi4sLli9fjvPnz0MQBKSnp5tdWCUiojvLYo9drVbjxo0bWLZsmWldeHg4cnJyEB8fDz8/PyxatAivvvoqqqqq8MQTT2DatGk2LZqIiBpnMdiTk5ORnJxcb31ERITp96CgIAQFBVm3MiIiahHeeUpEJDEMdiIiiWnydEciolp79+7FH3/80ezXrV+/vsnbOjg4ICoqCp06dWr2cewdg52Imu3111/H0KFD0blz5ya/5p577oFWq23y9l9++SV8fX0xZMiQlpRo1xjsRNQi7733Hnx8fJq8fXN66wCQn5/f3JLo/3GMnYhIYhjsREQSw2Bvo+RyOY4fP44jR47YZP/V1dVYt24dunXrZpP9E5F4GOxtVL9+/fDJJ58gNDQUGRkZVt33lStXMHbsWJw9exaZmZlW3TcRiY/B3oYFBwcjOzsbSUlJWLhwIQRBaPU+CwoKMGTIEDz88MP4+uuv0aVLFytUSkRtCYO9jfPz88PRo0exf/9+REZGorKystFtL1dcwYKclfhPxZUG248cOYLBgwdjxowZUKlUDT6Fk4jaPwZ7OyCXy3Hw4EEYjUaMGDECJSUlDW6385d9OHnhNHb8sq9em0ajwbhx47B582bEx8dDJpPZumwiEgmDvZ1wdXWFRqOBQqFAYGBgvfbLFVdwsOAIBAj4oeCIWa/9m2++wUsvvYQDBw40+OhlIluz9GmSrIufxduR/fv3IysrCwcTxuBsyiSztq8870JNZ1eggwzGKgO2bIvH+LKrAIBHAWS/FYx58+YhMzOT4+p0x9X9NPnygAjLL6BWYbC3E2vXrsWSJUuwe/du+A4ebNZ2ueIKjn/9DozGKgCAsYMMx7t2QUz0GtztejPEH6quxiNvvIHBgwdj7969+Nvf/nbH3wNJh7OzM86ePdukO09v/TQZ1ifYdF42pqKiAkVFRc3+aki6iUMxbVx1dTVmzZqFDRs24PDhwxh8S6gDN3tDglBjtq5GqDEba3d0dIRKpcLMmTMxePBgHD582Oa1k3StWLECUVFRKCwstLht3fPz1vOyIYIgICYmBoMHD0b//v2tUq+9YbC3YbXzzU+fPo3Dhw832sv+34tnUV1jNFtXXWPE/148W2/buLg4qNVqjB8/3urz48l+jBs3Dq+99homTJiAioqKRrer7a3Xnp/VNcZ614ButWzZMpw5cwYff/wxL/K3EIdi2qiysjIMHz4cw4cPx4cffnjbqYmpQcpm7bt2fvzYsWNRWFiIhISE1pZLdmjevHk4ceIEpk+fju3btzcYwrf7NNnQWPuePXuwbt065OXlwdXV1Wa1Sx177G3Uzz//jK5du9psvrmfnx+++OILpKWlWX3fZB9kMhnUajVOnjyJ5cuXN7hNcz5N/vLLL4iNjcXOnTv5qItWYo+9DbN1j4U9ImqtTp06YdeuXXj66afRp0+fetNpm/pp8uLFixg3bhxWrVqFp59+2hal2hX22ImoVby9vbFjxw5MmzYNOp2u2a+vqqrC5MmTERYWhhdffNEGFdofBjsRtdqgQYOwfPlyjBs3DpcuXWrWa9988024urpiyZIlNqrO/nAohoisIjg4GG+++SZK097Gf643PdxDna6jNDIKDg4ONqzOvjDYiajVDAYDwsLCMHPmTDw2J6VZr734P/+DKSEhePzxx+Hr62ujCu1Lk4ZiVCoVQkJCEBISgtTU1HrtOp0OEydORFBQEJRKJaqrq61eKBG1TYIgIC4uDl26dMHixYub/fonn3wSq1atQmhoKMrKymxQof2xGOyHDx9Gbm4uvvrqK+zatQu//PILvvvuO7NtEhISMH/+fBw4cACCIPDLG4jsyEcffYRDhw7h008/RYcOLbtsFxkZicmTJ2Py5MmoqqqycoX2x+JQjKenJ+bOnQsnJycAQI8ePVBUVGRqLywsRGVlJfr16wcAmDhxItasWYMpU6aY7ae8vBzl5eVm64qLi1tbv2QVFhbiwIEDpn9XW6ioqECnTp1stn+SvpycHCxatAiHDh2Cu7t7q/aVkpKC0NBQvPHGG1i3bp2VKrRPFoP9kUceMf1eUFCAb775BhqNxrSutLQUnp6epmVPT88GnxeelpYGlUrV2nrtxksvvYRevXqhc+fOzXpdWFgYduzY0eTt6/63I2qOs2fPYsqUKcjIyECPHj1avT8HBwdkZGRg4MCB2LBhA2bOnGmFKu1Tky+enjp1Cq+88goSExPRvXt30/qamhqzW4kFQWjw1uKpU6diwoQJZuuKi4sRGRnZgrKlz8HBAc8880yzX3fmzBmb9vKJAOCvv/7CuHHj8M4772DkyJFW26+7uzv27NmDIUOGQKFQYNiwYVbbtz1pUrAfO3YM8fHxSEpKqndnmVwux4ULF0zLZWVl8PLyqrcPd3f3Vn9UI6K2Yfr06RgyZAhee+01q++7Z8+eSE9PxwsvvID8/Hzcd999Vj+G1FkM9j///BOvv/46Vq1ahUGDBtVr79atG5ydnXHs2DH0798fu3fvxtChQ21SLBG1DUeOHMGPP/5os6cvPvvss3j44Ydx+vRpBnsLWAx2tVqNGzduYNmyZaZ14eHhyMnJQXx8PPz8/LBixQokJyfj6tWr6NOnD6Kjo21aNBGJz9aP1OUje1vOYrAnJycjOTm53vqIiP8+cvOxxx5r1gU7IiKyHT4rhohIYvhIASJqtj/++AN+fn6Qy+U2PUbt/TPUPAx2Imq2goIC3Lhxo1mvGTNmDPbu3dvk7R0dHfml6y3EYCeiZvPx8Wn2a06dOoVevXrZoBq6FcfYiYgkhsFORCQxDHYiIolhsBMRSQyDnYhIYhjsREQSw2AnIpIYBjsRkcQw2ImIJIbBTkQkMQx2IiKJYbATEUkMg52ISGIY7EREEsNgJyKSGAY7EZHEMNiJiCSGwU5EJDEMdiIiiWGwExFJDIOdiEhimhTsV69exZgxY6DX6+u1qVQqjBgxAqGhoQgNDUV6errViyQioqZztLRBfn4+kpOTUVBQ0GC7VqvFypUr4e/vb+3aiIioBSwGe2ZmJhYsWIDExMQG27VaLTZu3IjCwkI8+eSTePvtt+Hs7Fxvu/LycpSXl5utKy4ubmHZRETUGIvBnpKS0mjbtWvXoFAokJCQAB8fH8ydOxfr16/HnDlz6m2blpYGlUrVumqJiMgii8F+O25ubti8ebNpOSYmBklJSQ0G+9SpUzFhwgSzdcXFxYiMjGxNCUREdItWBXtRUREOHz6MsLAwAIAgCHB0bHiX7u7ucHd3b83hiIioCVo13dHFxQXLly/H+fPnIQgC0tPT8dxzz1mrNiIiaoEWBfv06dNx4sQJeHh4YNGiRXj11VcxevRoCIKAadOmWbtGIiJqhiYPxeTk5Jh+rzuuHhQUhKCgIOtWRURELcY7T4mIJIbBTkQkMQx2IiKJYbATEUkMg52ISGIY7EREEsNgb4d8fX0hk8ka/AHQ4HpfX1+RqyaiO6VVjxQgcWi1WrFLIKI2jD12IiKJYbATEUkMg52ISGIY7EREEsNgJyKramzWFsAZW3cKg52IrEqr1UIQBNNPRkYG+vTpgw4dOqBPnz7IyMgwa+csL+tjsEuERqOBr68vHBwc4OvrC41GI3ZJRNBoNFAqlVi7di0qKyuxdu1aKJVKnp82xmCXAP7xUFuVkpICtVqNESNGoGPHjhgxYgTUajVSUlLELk3SGOwSwD8eaqt0Oh0CAgLM1gUEBECn04lUkX1gsEuATqeDXq83G4rR6/X84yHRKRQK5Obmmq3Lzc2FQqEQqSL7wGCXgAceeACJiYlmQzGJiYl44IEHxC6N7JxSqURsbCwOHjyIqqoqHDx4ELGxsVAqlWKXJml8VoxEVFZWIiYmBufOnYOPjw8qKytx1113iV0W2bmIiAgAQFxcHHQ6HRQKBVJSUkzryTYY7BJQWFiIrl27AoBpvrCjoyMKCwvFLIsIwM1wZ5DfWRyKkQAnJyfMmzcPv//+O4xGI37//XfMmzcPTk5OYpdGRCJgsEuAwWCASqUyG8dUqVQwGAxil0ZEIuBQjAT07t0b48ePNxvHnDJlCnbt2iV2aUQkAvbYJUCpVCIjI8NsVkxGRgZnHhDZqSb12K9evYrw8HBs2LAB3t7eZm06nQ5KpRLXrl3DgAED8O6778LRkR8E7iTOPCCiuiz22PPz8xEREYGCgoIG2xMSEjB//nwcOHAAgiAgMzPT2jVSE0RERECr1cJoNEKr1TLUieyYxWDPzMzEggUL4OXlVa+tsLAQlZWV6NevHwBg4sSJ2L9/f4P7KS8vh16vN/spLi5uXfVERFSPxTGT2z1vpLS0FJ6enqZlT09PlJSUNLhtWloaVCpVC0qkptBoNEhJSTENxSiVSvbaiexUqwbDa2pqTDfEAIAgCGbLdU2dOhUTJkwwW1dcXIzIyMjWlEC4GeqzZ8+Gm5sbBEHAtWvXMHv2bABguBPZoVbNipHL5bhw4YJpuaysrMEhGwBwd3eHt7e32Y9cLm/N4en/JSYmwsHBAVu2bMGNGzewZcsWODg4IDExUezSiEgErQr2bt26wdnZGceOHQMA7N69G0OHDrVKYdR0er0e27ZtM3ts77Zt26DX68UujYhE0KJgnz59Ok6cOAEAWLFiBZYuXYrRo0fj+vXriI6OtmqBRETUPDJBEASxDq7X6xEYGIjs7Ox68+Op6R588EEYjUakp6cjICAAubm5iIyMhIODA86fPy92eURkZZayk3eeSkBqaiqqq6sRExMDFxcXxMTEoLq6GqmpqWKXRkQiYLBLQEREBFavXg03NzcAgJubG1avXs0ZMUR2ivf+SwSfeU1EtdhjJyKSGAY7EdmURqMx+6J1jUYjdkmSx6EYIrIZjUYDpVIJtVptmrEVGxsLgHdF2xJ77ERkMykpKVCr1WY3z6nV6ts+g4paj8FORDaj0+kQEBBgti4gIAA6nU6kiuwDg52IbEahUCA3N9dsXW5uLhQKhUgV2QcGOxHZjFKpRGxsrNkXrcfGxvJrG22MF0+JyGb4tY3iYLATkU3x5rk7j0MxREQSw2AnIpIYBjsRkcQw2ImIJEbUi6dGoxHAzS+1JiKipqnNzNoMvZWowV77RdiRkZFilkFE1C5duHABPj4+9daL+tV4lZWV0Gq18PT0hIODg1hlSEZxcTEiIyORnp4OuVwudjlEJjw3rctoNOLChQvw9fWFi4tLvXZRe+wuLi4YMGCAmCVIklwu53fIUpvEc9N6Guqp1+LFUyIiiWGwExFJDIOdiEhiGOwS4u7ujlmzZsHd3V3sUojM8Ny8s0SdFUNERNbHHjsRkcQw2ImIJIbBLhFZWVkIDg7GqFGjkJ6eLnY5RGauXr2KMWPGQK/Xi12KXWCwS0BJSQlWrVqFjIwM7Nq1C59//jlOnz4tdllEAID8/HxERESgoKBA7FLsBoNdAg4fPoyBAwfi7rvvRqdOnRAUFIT9+/eLXRYRACAzMxMLFiyAl5eX2KXYDX41ngSUlpbC09PTtOzl5YV///vfIlZE9F8pKSlil2B32GOXgJqaGshkMtOyIAhmy0RkXxjsEiCXy02PQAZuPsqTH3uJ7BeDXQIGDx6MI0eO4NKlS6ioqMC3336LoUOHil0WEYmEY+wScN9992HOnDmIjo5GVVUVwsLC8Pjjj4tdFhGJhI8UICKSGA7FEBFJDIOdiEhiGOxERBLDYG9DjEYjZs2ahb/++gvvvvsuxowZg7FjxyIyMrJVNxzNmzcPhYWFAICRI0c2+LyOZcuW4ddff23xMcg+1J6jqampCAoKQmVlpaktLy8PUVFRzdpfZmYm9u7dCwCYO3cuvvzyy3rbfPvtt/j0009bV7idYbC3IRqNBgEBAdBoNKipqUFWVhaysrLw5ptv4rXXXkNVVVWL9puXlwdL18hnzJiBJUuWtGj/ZD9qz1FXV1cUFhZi5cqVrdrf8ePHYTAYbrvNqFGj8O233+LixYutOpY9YbC3EYIgYPv27QgJCUFZWRmqqqpMQd6/f38sWbIENTU1AIANGzYgODgYY8eOxbJly2A0GqHX6zFy5EjT/tauXYu1a9di06ZNKC0txYwZM3D58mUAwLp16zB+/HgEBQUhPz8fAODh4QEPDw8cPXr0Dr9zai/qnqMA8MILL2Dfvn34+eef621bVlaGV155BWPHjsWECRPwj3/8A8DN8zI2NhbBwcHYvn07cnJysGbNGvzzn/8EAPzwww8ICwvDiBEj8Pnnn5v2x6eWNg+DvY04efIkOnfujM6dOyM6Ohr5+fkYNGgQXn31VWzbtg3+/v5wdnbGjz/+iJycHOzcuRNfffUVzp07h88++6zR/c6YMQNeXl7YtGkT7rnnHgBAz549sWvXLkRFRUGtVpu2HTBgAHJycmz+Xql9qnuOAsDdd9+NhQsXQqlUmg3JAMDixYsxcOBAZGVlYc2aNUhKSkJZWRkAwGAwYN++fYiKisLIkSMRHx+PZ555xtT2xRdfYOPGjVi1apVpfzw3m4fB3kYUFBRALpcDALy9vbF3715s3boVffv2xa5duxAaGory8nIcPXoUISEhcHV1haOjIyZNmoQjR44061jPPvssgJsBX9uLB4AHHngA586ds96bIkmpe47WevbZZ+Hr61tvSObo0aMICwsDADz44IPo27ev6dPh7W6eCwwMhEwmwyOPPGJ2bnbr1o3nZjMw2NsImUwGR8ebNwKvXLkSpaWlePzxxzFz5kx8+eWX8PLywqFDh0zDMXVVV1dDJpOZjaNXV1c3eiwHBwfTMetydHTkw8OoUXXP0bqSk5Oxb98+HDt2zLTu1ms6giDAaDQCAFxcXBo9Bs9N62CwtxE+Pj6mmSslJSVYt26d6aLShQsXcOnSJfTq1QsDBw7E119/jcrKSlRXV2Pnzp0YOHAg3N3d8Z///AeXLl2CwWAwjVkCN/9Yav+obqewsBA+Pj62eYPU7tU9R+u65557sHDhQqxfv960buDAgdixYwcA4Pz58zh+/Dj69etX77VNPTf1ej3PzWbgs2LaiMceewyXL1/GX3/9hXfeeQfvv/8+Ro8eDVdXV3Ts2BF///vf0aNHD/To0QM6nQ6TJk1CdXU1AgIC8OKLL8LR0REvv/wywsLCIJfL4efnZ9r38OHDMWPGDHz88ce3rSEvLw8vvviird8qtVN1z9FbPfvsswgKCkJpaSkAQKlUYv78+abpi++9916DTxwdPHgwVq5caRq3b0xeXh4CAwOt8C7sA58V04Zs27YNHTp0ECVcL168iFmzZkGj0dzxY1P7IdY5GhERAZVKha5du97R47ZXHIppQyIiInDo0CFUVFTc8WNv2LABSUlJd/y41L6IcY7u378fQUFBDPVmYI+diEhi2GMnIpIYBjsRkcQw2ImIJIbBTkQkMQx2IiKJ+T/cmwBjfXQHuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.boxplot([df07surf['OmegaA'][iisouth],df07surf['OmegaA'][iinorth]],\n",
    "            labels=['0\\n(South)','1\\n(North)'],showmeans=True,notch=True)\n",
    "plt.title('$\\Omega_A$ - upper 10m 2007')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Method 1: t-test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a difference of 0.104 in the mean of $\\Omega_A$ between the two regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.03056105006487364"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(df07surf['OmegaA'][iinorth]) - np.mean(df07surf['OmegaA'][iisouth])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is this difference statistically significant? A Student's t-test can be used to test whether the null hypothesis of no difference can be rejected at the 95% confidence level. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=-0.3235609219137568, pvalue=0.7467674907175257)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ttest_ind(df07surf['OmegaA'][iinorth],df07surf['OmegaA'][iisouth],equal_var=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Welch's t-test relaxes the assumption of equal variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=-0.30410900939221397, pvalue=0.7619693293723315)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ttest_ind(df07surf['OmegaA'][iinorth],df07surf['OmegaA'][iisouth],equal_var=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Method 2: generalized linear model\n",
    "\n",
    "This test can also be framed in terms of a general linear model\n",
    "\n",
    "$$ \\hat{y} = \\hat{a}_1 + \\hat{a}_2 x .$$\n",
    "\n",
    "Again, $\\hat{y}$ is a model for the aragonite saturation state data. In this case, we can think of the southern data points as having $x = 0$ and the northern data points as having $x = 1$.\n",
    "\n",
    "In this model, the intercept parameter $\\hat{a}_1$ is the mean of the points with $x = 0$, the southern points.\n",
    "\n",
    "The slope parameter $\\hat{a}_2$ is equal to the difference between the means of the two groups.\n",
    "\n",
    "$$ slope = \\frac{\\Delta\\bar{y}}{\\Delta x} = \\frac{\\Delta\\bar{y}}{1} = \\Delta\\bar{y}$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = smf.ols(formula=\"OmegaA ~ is_northern\", data=df07surf).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are summarized below. The slope parameter $\\hat{a}_2$ in our model is the coefficient for the `is_northern` variable. This is a Boolean variable that is equal to 0 (`False`) for southern points and 1 (`True`) for northern points. Notice that this coefficient is equal to 0.104, which the same as the difference between the two means.\n",
    "\n",
    "Also notice that the 95% confidence intervals (shown as the `[0.025 0.975]` interval) overlap 0 for this parameter. This means that the difference is not statistically significant at the 95% confidence ($\\alpha$ = 0.05) level. This summary also shows a t-statistic and p-value for this parameter, which are equivalent to the Student's t-test result shown above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>OmegaA</td>      <th>  R-squared:         </th> <td>   0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>  -0.007</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>  0.1047</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 08 Mar 2022</td> <th>  Prob (F-statistic):</th>  <td> 0.747</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>14:09:53</td>     <th>  Log-Likelihood:    </th> <td> -102.03</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   138</td>      <th>  AIC:               </th> <td>   208.1</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   136</td>      <th>  BIC:               </th> <td>   213.9</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "           <td></td>              <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>           <td>    2.2110</td> <td>    0.052</td> <td>   42.433</td> <td> 0.000</td> <td>    2.108</td> <td>    2.314</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>is_northern[T.True]</th> <td>   -0.0306</td> <td>    0.094</td> <td>   -0.324</td> <td> 0.747</td> <td>   -0.217</td> <td>    0.156</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 7.460</td> <th>  Durbin-Watson:     </th> <td>   0.655</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.024</td> <th>  Jarque-Bera (JB):  </th> <td>  12.047</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.189</td> <th>  Prob(JB):          </th> <td> 0.00242</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.397</td> <th>  Cond. No.          </th> <td>    2.42</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 OmegaA   R-squared:                       0.001\n",
       "Model:                            OLS   Adj. R-squared:                 -0.007\n",
       "Method:                 Least Squares   F-statistic:                    0.1047\n",
       "Date:                Tue, 08 Mar 2022   Prob (F-statistic):              0.747\n",
       "Time:                        14:09:53   Log-Likelihood:                -102.03\n",
       "No. Observations:                 138   AIC:                             208.1\n",
       "Df Residuals:                     136   BIC:                             213.9\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "=======================================================================================\n",
       "                          coef    std err          t      P>|t|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------------\n",
       "Intercept               2.2110      0.052     42.433      0.000       2.108       2.314\n",
       "is_northern[T.True]    -0.0306      0.094     -0.324      0.747      -0.217       0.156\n",
       "==============================================================================\n",
       "Omnibus:                        7.460   Durbin-Watson:                   0.655\n",
       "Prob(Omnibus):                  0.024   Jarque-Bera (JB):               12.047\n",
       "Skew:                           0.189   Prob(JB):                      0.00242\n",
       "Kurtosis:                       4.397   Cond. No.                         2.42\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing two sample means (another example)\n",
    "\n",
    "The following shows the same calculations for temperature. In this case, the null hypothesis can be rejected at the 95% confidence level, and the 95% confidence intervals for the model slope do *not* overlap with zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'CTDTMP - upper 10m 2007')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEJCAYAAACJwawLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAff0lEQVR4nO3de1RU5f4G8GdA5CJOSoBiqKc0IsP0mBVeU6C8ACJiKiKhiZeTYkoRKBBKEogmpySPeTIjtdQEM0RzpdAFr5lHDDWTH2KNOojXcRQEZvbvjxZTIyAz4wzMZp7PWqzlvPvd7/6OjM+8vrP3HokgCAKIiEi0rFq6ACIiejAMciIikWOQExGJHIOciEjkGORERCLHICciErk2LV0ANS+VSoXPPvsMubm5UKlUqKmpwfDhw/H666+jqqoK4eHhAIA7d+6gvLwcjz76KABg4MCBGDZsGGbMmKFpU6vV6NixI2bPno2BAweipKQEb7zxBgDg5s2buHXrFtzd3QEAwcHBkEqlWLhwIebMmYN58+ZpahIEAX5+frC3t8fOnTuRk5ODlJQUuLu7QyKRQBAE2NvbIzY2Fv/85z+b86/L6LZt24a9e/dizZo1Wm2ffPIJamtrMWDAACQkJMDGxsbgY3zyySfIzs6GtbU1nJyckJycjG7duqGyshIJCQk4deoU1Go1YmJi4OfnBwAoKipCcnIy7ty5A1dXVyxfvhyurq5YunQpfvrpJ83Y5eXlcHFxQW5uruF/CWR8AlmUhIQEISoqSlAoFIIgCMLt27eFf/3rX8Kbb76p1e/QoUOCv79/k22nT58WBg4cKBw/flyrPTs7W5g5c2a9tmHDhgm+vr5a7UeOHBEGDhyoGbuhffft2ycMGjRIqKmp0fMZm4fr168LiYmJQt++fbWe25kzZ4ShQ4cKV69eFVQqlbBgwQJh7dq1Bh9n//79wqhRo4Rbt24JgiAIGzduFCZPniwIgiAsW7ZMSEhIEARBEC5cuCAMHjxYuHTpknD37l1h6NChwtGjRwVBEIRNmzYJkZGR9cb+448/hKFDhwqnTp0yuD4yDS6tWBCZTIbc3Fy8++67aN++PQDAwcEBS5Ys0czM9OXp6Ynw8HB8+umnOvX38PCAg4MDjh07pmnbvn07xowZc9/9BgwYgIqKCigUCp1rk8lkWjP4vz/OycnBzJkzERkZidGjR2PatGkoLy8HAISHhyM1NRXjx4+Hr68vPvjgA80Yx44dw+TJkxEcHIyQkBAUFBRoxqtrr/tfzd/t3r0brq6uiI2N1Wrft28ffHx84OTkBCsrK0ycOBFff/21po60tDRMmjQJL730Ej7++GOkpaVh3LhxGDVqFM6cOVPvOM7Ozli8eDEcHR0BAL1798bFixcBAHv37sXLL78MAOjSpQsGDRqE3bt345dffoGjoyOeeeYZAMD48eNx8OBBXL9+XWvsxMRETJs2DU8++aSuvwJqJgxyC3Ly5En07NlT84+8jouLC0aMGGHwuJ6envjtt9907j927Fjs2LEDAFBZWYmff/4ZQ4YMabS/IAjYsmULPDw84OTkZHCd9/rpp58QHx+PXbt24amnnkJKSopm27lz5/DFF19g+/bt2LVrFwoKCnDz5k0sXLgQ6enp2L59O1avXo3FixdrgrKkpAQbNmzAhg0b6h0rNDQUc+fORdu2bbXaL126BDc3N83jzp07a95QAODChQvYvHkzli9fjuXLl+O5555DTk4OhgwZgo0bN9Y7joeHB5577jkAQHV1NVasWIGRI0c2eKxOnTpBLpdDLpejc+fOmva2bdvCyclJq47vv/8eFy9ebPBNiloe18gtiJWVFdRqtdHHlUgksLOz07l/YGAggoKCEB8fj2+//RY+Pj6wtrbW6nP06FEEBQVBIpGguroajz32mNbM2BgGDRqkWe+fMGECgoKCNNsmTpwIGxsb2NjYYOTIkSgsLISVlRUqKiowZ84cTT+JRKKZGT/xxBP13iSbItxzhwxBEGBl9df86sUXXwQAdO3aFQA0b3jdunXDkSNHGh332rVrmDdvHhwdHbFgwQLN2BKJRKtf3Wvi3nZBELR+J1lZWZg1a1a93xOZBwa5BXn66adRWloKpVKpFTjl5eVITEzEBx98oFcg1/nll1/g4eGhc38XFxf06tULP/zwA7766ivExcXV+298//798dFHHzU51owZM3D58mUAwLx58+Dr66vZVvdBaZ2amhqtff8eSmq1WutxmzZ//dOoC1eVSoUePXrgyy+/1GwrLy+Hk5MTcnNz4eDg0GS993Jzc9PUDwCXL1+uNzv+O10+BP3111/x2muvwc/PD7GxsZrnVXcsZ2dnzbE8PT3r1VBTU4MbN26gU6dOAP58UygqKkJmZqbez4+aB5dWLEinTp0QGBiIRYsWQalUAgCUSiUWL16MDh06GBTiJ06cwBdffIGIiAi99hs7dizWr1+PW7du6fUmcK///ve/2LFjB3bs2KEV4gAglUpRU1ODkpISAEBeXp7W9kOHDmmWDzZv3ozhw4drtn399ddQq9W4efMmdu/eDR8fH/Tt2xfnz5/XnMVx+vRpjBgxQmsJQl8+Pj7Iz8/H1atXNUtIhn5eAQByuRwRERF47bXXsGjRIq03J19fX2zZskXT78cff8Tw4cPRp08f3LhxQ/O5RXZ2Nvr27QupVArgz88FevfubdAbFTUPzsgtTFJSElavXo1JkybB2toa1dXV8PPzQ1RUlE77//7775olCCsrKzg6OmLFihXw9PTUqw4/Pz8kJSVp/ttvCu3bt0dMTAxmzJgBJycnzVpxnU6dOiEmJgYVFRXo2bMnkpOTNduqqqowfvx43L59G5MnT8aAAQMAAB988AHS09Nx9+5dCIKA9PR0uLu733eZ4348PT0xZ84cREREoKamBn369MGMGTMMfs6rV69GZWWl1lp927Zt8eWXXyIqKgqLFy+Gv78/VCoVYmJi0K1bNwBAZmYmkpOTUVlZiQ4dOmDZsmWaMcvKyvDII48YXBOZnkS4d5GOyALk5ORgz549DS7fhIeHIywsrF7wE5krLq0QEYkcZ+RERCLHGTkRkcgxyImIRK7Zz1qpqqpCcXExXFxceHEBEZGOVCoVKioq4OXlVe9U4WYP8uLiYoSFhTX3YYmIWoVNmzahf//+Wm3NHuQuLi6aYv5+BRsRETVOLpcjLCxMk6F/1+xBXrec0rlzZ829qomISDcNLUnzw04iIpFjkBMRiRyDnIhI5BjkREQixyAnIhI5BjkRkcgxyImIRI5BLhJeXl6QSCR6/Xh5ebV02UTUDPgNQSJRXFzc6LZ7v5uSiCwLZ+RERCLHICciEjkGORGRyDHIiYhETqcgVyqVCAgIgEwmAwAUFhZizJgxCAgIwFtvvYXq6mqTFklERI1rMsiLiooQGhqKsrIyTVt8fDwyMjKwc+dOVFVVYceOHaaskYiI7qPJIN+6dSuSkpLg6uqqaVOpVFAqlVCpVLh79y5sbW1NWiQRETWuyfPIU1JS6rUtXrwY4eHhcHR0hLu7O0aOHNngvgqFAgqFQqtNLpcbWCoRETVE7wuCKioqsGLFCuzcuRPu7u5ITU1FamoqkpKS6vXNyspCZmamUQolIqKG6R3kR48ehYeHB7p16wYAmDBhAubPn99g34iICAQHB2u11X3vHBERGYfeQe7h4YFly5bhypUrcHZ2xr59+9C7d+8G+0qlUkil0gcukoiIGqd3kPfo0QOvv/46XnnlFVhbW6N79+5ITk42RW1ERKQDnYM8Pz9f8+fg4OB6SyZERNQyeGUnEZHIMciJiESOQU5EJHIMciIikWOQExGJHIOciEjkGORERCLHICciEjkGORGRyDHIiYhEjkFORCRyDHIiIpFjkBMRiRyDnIhI5BjkREQixyAnIhI5BjkRkcgxyImIRI5BTkQkcgxyIiKRY5ATEYkcg5yISOQY5EREIscgJyISOQY5EZHIMciJiESOQU5EJHI6BblSqURAQABkMhkA4H//+x8mTJgAf39/REdHo7q62qRFEhFR45oM8qKiIoSGhqKsrAzAn6EeFRWF5ORk5OXlAQC2bdtm0iKJiKhxTQb51q1bkZSUBFdXVwDA/v370bdvX3h6egIAEhIS8OKLL5q2SiIialSbpjqkpKRoPT5//jwcHBywYMEClJaWol+/foiLi2twX4VCAYVCodUml8sfoFwiIrpXk0F+L5VKhcLCQmzZsgVdunRBfHw81q5di6ioqHp9s7KykJmZaZRCLUFxcTFWrVpl0L6zZs3Sue8LL7yAyZMnG3QcIjI/ege5s7Mz+vTpg65duwIARo0ahY0bNzbYNyIiAsHBwVptcrkcYWFhBpTa+hUWFuL06dN6//106dIF/fr106nviRMn8PnnnzPIiVoRvYN88ODBWLVqFS5dugQ3NzcUFBTgqaeearCvVCqFVCp94CItSa9evfSaXQP6zcZ37tyJNWvW6FsWEZkxvYPczc0NycnJmD17Nu7evYsnn3wSsbGxpqiNiIh0oHOQ5+fna/48bNgwDBs2zBT1EBGRnnhlJxGRyDHIiYhEjkFORCRyDHIiIpFjkBMRiRyDnIhI5BjkREQixyAnIhI5BjkRkcgxyImIRI5BTkQkcgxyIiKRY5ATEYkcg5yISOQY5EREIscgJyISOQY5EZHIMciJiESOQU5EJHIMciIikWOQExGJHIOciEjkGORERCLHICciEjkGORGRyDHIiYhEjkFuRqysrKBUKk16DKVSCSsr/tqJWhOd/kUrlUoEBARAJpNptW/cuBHh4eEmKcwSBQYGYvfu3SgtLTXJ+Gq1Gunp6Zg8ebJJxieiltFkkBcVFSE0NBRlZWVa7SUlJVi7dq2p6rJIbm5uWLBgAeLi4kwy/oYNG2Bra4uJEyeaZHwiahlNBvnWrVuRlJQEV1dXTVt1dTXefvttzJs3z6TFWaI33ngDhw4dQmFhoVHHvX37NuLj45GRkQGJRGLUsYmoZbVpqkNKSkq9tvfeew8hISFwd3e/774KhQIKhUKrTS6X61miZbG3t0dqaiqio6Nx6NAho61nL1++HEOGDIG3t7dRxiMi89FkkN9r//79uHTpEhYuXIjDhw/ft29WVhYyMzMNLs5ShYaG4v3338cXX3yBsLCwBx7vwoULWLVqFY4dO2aE6ojI3Ogd5Dt37sTZs2cRFBSEO3fu4MqVK5g/fz7+/e9/1+sbERGB4OBgrTa5XG6UcGrNrKyskJGRgdDQUAQHB8PBweGBxouPj8esWbPQvXt3I1VIROZE7yBPTU3V/Pnw4cPIzMxsMMQBQCqVQiqVGlycJRs0aBC8vb2xcuVKJCQkNNrveuVN/PvgOiwYMB0d7B+qt/3nn3/Gnj178Ntvv5myXCJqQTyh2IylpaUhIyMDFy9ebLRP9sld+LWiBNtO7qq3TRAEREdHIzk5Ge3btzdlqUTUgnQO8vz8/Hofbj7//PPYsGGD0YuiPz322GOIjIxEYmJig9uvV95EQdlBCBDwXdlB3Ki8qbX9q6++wvXr1/Hqq682R7lE1EI4IzdzixYtQl5eHo4fP15vW/bJXRAENQBALai1ZuXV1dWIiYnBypUrYW1t3VzlElELYJCbuYceeggJCQn11snrZuO1ahUAoFat0pqVr1+/Ho8//jj8/PyavWYial56f9hJze/s2bN4u387lKaEaNq2uzhC3d4esPrr4h5VTTU++Wwexl5R4kUAj3tZoaamBjY2Ni1QNRE1Fwa5mTtz5gw2bdqEhNOn4eLiomkv35MC1Q3te9+orCSQ/8MDj82KhyAI+NfIkTi5Zg2ioqKau2wiakYMcjMXExOD2NhYrRAHgPQR8ffdTyKRYMWKFfD19cWUKVPQsWNHU5ZJRC2Ia+RmbN++fSguLjb4nja9e/dGcHAwli5dauTKiMicMMjNlEqlwhtvvIH09HTY2toaPE5ycjKysrJw9uxZI1ZHROaEQW6mPv30U7Rv3x4hISFNd76PTp064c0330RsbKyRKiMic8MgN0O3bt1CYmIiVq5caZRbzs6fPx/Hjh3Dd9999+DFEZHZYZCbofT0dPj6+uLZZ581ynh2dnZIS0tDdHQ01Gq1UcYkIvPBIDczv//+O1avXo13333XqONOnDgRtra2vKUCUSvEIDczixYtwpw5c9C1a1ejjiuRSJCRkYH4+Hjcvn3bqGMTUctikJuRn3/+GQUFBXjrrbdMMr63tzeGDBmC9957zyTjE1HLYJCbkZ9++gn+/v5wdHQ02TFCQ0Nx5MgRk41PRM2PQW5mjPUdnS01PhE1P/6rJiISOQY5EZHIMciJiESOQU5EJHIMciIikWOQExGJHIOciEjkGORERCLHICciEjkGORGRyDHIieiBeXl5QSKR6Pzj5eXV0iW3Km1augAiEr/i4uIG2yUSCQRBaOZqLI9OM3KlUomAgADIZDIAwJYtWxAQEIDAwEAsXLgQ1dXVJi2SiIga1+SMvKioCAkJCSgrKwMAnDt3DuvWrUNOTg7atWuHuLg4fP7555g6daqJS2391Go1SktLUVpaarJjyOVyk41NRC2jySDfunUrkpKSNF920LZtWyQlJWnume3h4YGLFy+atkoLUV5ejm+//RZ+fn567SeTyeDu7q5z//Hjx+tbGhGZsSaDPCUlRevxI488gkceeQQAcO3aNWzatAmpqakN7qtQKKBQKLTaOCNs3JIlS7BkyRK995NIJCadxROReTP4w87y8nJERkYiJCQEzz//fIN9srKykJmZaXBxRETUNIOC/P/+7/8QGRmJ8PBwvPrqq432i4iIQHBwsFabXC5HWFiYIYclIqIG6B3kSqUS06dPx/z58zF27Nj79pVKpZBKpYbWRkREOtD7gqBt27bhypUrWL9+PYKCghAUFIT333/fFLUREZEOdJ6R5+fnAwCmTp3KUw2JiMwIL9EnIhI5BjkRkcgxyImIRI43zSKiJgmCgHnz5uHatWt676vP6cZ2dnb48MMPYWdnp/dxLBmDnIh0kpmZiY0bN+q1zw8//IDRo0fr3H/u3LlITk7WXD1OumGQE5HO9L2YT9/+dfd0Iv1wjZyISOQY5ESkEzs7O5w6dcpk48tkMigUCtja2prsGK0Vg5yImiSRSPDRRx8hICAAly9fNvr4SqUSgYGBSEhIgLOzs9HHb+0Y5ESkk1deeQVhYWEICgpCZWWl0cZVqVQIDQ3FM888wzVyAzHIiUhnycnJ+Mc//oGpU6dCrVYbZczo6GhUVlbiP//5DyQSiVHGtDQMciLSmUQiwfr16yGTyZCQkPDA461atQp79+7Ftm3bYGNjY4QKLRNPPyQivdjZ2eGrr76Ct7c3Hn/8cUybNs2gcfLy8pCamor9+/ejQ4cOxi3SwjDIiUhvLi4uyMvLwwsvvIDu3bvDx8dHr/2PHz+OqVOnIjc3F48++qiJqrQcXFohIoN4enpi8+bNmDBhAmQymc77KRQKBAQE4MMPP4S3t7cJK7QcnJETkcHOnDmDhx9+GKrcdJReu6DzfpvC+qHAhOekWxoGOREZZM+ePVi8eDH279+P7j166LVvu/JyTPX2Rs+ePTFlyhQTVWg5GOREpLdffvkF4eHhyMnJQQ89QxwAOnXqhLy8PAwbNgzdunXD0KFDTVCl5eAaORHpRS6XIzAwEO+//z4GDx7caL/rlTeRlL8SNypvNri9V69e+PzzzzFhwgScPXvWVOVaBAY5Eenszp07CAwMxPTp0xEaGnrfvtknd+HXihJsO7mr0T5+fn5YunQp/P39cfXqVWOXazEY5ESkE7VajSlTpqBXr15NXgx0vfImCsoOQoCA78oONjorB4DIyEgEBwcjODgYd+/eNXbZFoFBTkQ6iYuLw7Vr17B27domL6XPPrkLgvDnJfxqQX3fWTkApKamwtXVFZGRkRAEwWg1Wwp+2ElETRIEAStWrMCVK1eavM1s3Wy8Vq0CANSqVfiu7CDGPzUaHewfanAfKysrbNiwAV26dIFcLoebm5vRn0Nrxhk5EelEEAQ4OTk12e/vs/E6uszK7e3t4eDgYLSbcVkSBjkRGdVvV0s1s/E6tWoVfrta2kIVtX5cWiEio0ofEd/SJVgczsiJiEROpyBXKpUICAjQ3BjnwIEDCAwMxEsvvYSMjAyTFkhERPfXZJAXFRUhNDQUZWVlAICqqiosWrQIq1evxq5du1BcXIzvv//e1HUSEVEjmgzyrVu3IikpCa6urgCAEydOoHv37ujatSvatGmDwMBAfPPNNw3uq1AoIJPJtH7kcrlxnwERkYVr8sPOlJQUrceXL1+Gi4uL5rGrqyvKy8sb3DcrKwuZmZkPWCIREd2P3metqNVqrau6BEFo9CqviIgIBAcHa7XJ5XKEhYXpe1giImqE3kHeuXNnVFRUaB5XVFRoll3uJZVKIZVKDa+OiIiapHeQ9+nTB+fOncP58+fh7u6OnTt3IiQkxBS1EZGZ+fjjj006/u3bt006fmuld5Db2toiLS0NUVFRuHv3Ll544QWMHDnSFLURkRlp164dDh06pNc+2dnZek30pkyZAmdnZ31Ls3gSoZlvNSaTyeDr64t9+/bB3d29OQ/dakkkEt4xjswSX5vGc7/s5JWdREQixyAnIhI5BjkRkcgxyImIRI5BTkQkcgxyIiKRY5ATEYkcg5yISOQY5EREIscgJyISOQY5EZHIMciJiESOQU5EJHIMciIikWOQExGJHIOciEjkGORERCLHICciEjkGORGRyDHIiYhEjkFORCRyDHIiIpFjkBMRiRyDnIhI5BjkREQixyAnIhI5BjkRkcg9UJDv2LED/v7+8Pf3x7Jly4xVExER6cHgIK+srERKSgo2bNiAHTt24OjRozhw4IAxayMiIh20MXRHlUoFtVqNyspKODg4oLa2Fra2tlp9FAoFFAqFVptcLjf0kERE1ACDg9zR0RGvv/46Ro0aBXt7ezz77LPo16+fVp+srCxkZmY+cJFERNQ4g5dWfv31V2RnZ6OgoAA//vgjrKyssG7dOq0+ERER2Ldvn9bPpk2bHrhoIiL6i8FBXlhYiAEDBuDhhx9G27ZtMW7cOBw5ckSrj1Qqhbu7u9ZP586dH7hoIiL6i8FB7unpiQMHDuDOnTsQBAH5+fno3bu3MWsjIiIdGLxGPnjwYJw6dQrjxo2DjY0NevfujZkzZxqzNiIi0oHBQQ4AM2fOZHgTEbUwXtlJRA/My8sLEomk3g+ABtu9vLxauOLW5YFm5EREAFBcXNzSJVg0zshForEZD2c9RMQZuUhwxkNEjeGMnIhI5BjkREQixyAnIhI5BjkRkcgxyImIRI5BTkQkcgxyIiKRa/bzyFUqFQB+UxARkT7qMrMuQ/+u2YO8oqICABAWFtbchyYiEr2Kigp0795dq00iCILQnEVUVVWhuLgYLi4usLa2bs5Dt0pyuRxhYWHYtGkTv7SDzApfm8alUqlQUVEBLy8v2NnZaW1r9hm5nZ0d+vfv39yHbfU6d+4Md3f3li6DqB6+No3n3pl4HX7YSUQkcgxyIiKRY5ATEYkcg1zkpFIp5s6dC6lU2tKlEGnha7P5NPtZK0REZFyckRMRiRyDnIhI5BjkZuabb77BuHHjMGbMGAQGBuLjjz82aJyCggKsX78eALBq1SqsWrXKmGWSBZLJZHjiiSewf/9+rXYfHx/IZDKdxvjjjz+waNEiAMDhw4cRHh5u9DotEb+z04yUl5dj2bJlyMnJQceOHXH79m2Eh4fj0Ucfha+vr15j8Ts+yRRsbGyQmJiIr7/+Go6Ojnrvf/HiRfzxxx8mqMyycUZuRq5fv46amhpUVVUBANq1a4e0tDT07NkTx48fx8svv4wxY8YgIiIC58+fBwCEh4fj8OHDAP6cMfn4+KCkpASbN2/G5s2bkZ2dDQA4ceIEJk2ahOHDh3N2TgZzdXXFwIEDsWzZsnrb1qxZg9GjRyMwMBBpaWlQqVSQyWQYOXIkQkNDMW3aNCxduhTFxcVYsmQJAODatWuYMWMGRowYgdmzZ6O6urq5n1KrwCA3I56envD19YWfnx/Gjx+P5cuXQ61Ww83NDdHR0ZqZ0KRJkxAdHd3oOD179sSkSZMwadIkhISEAACuXr2Kzz77DNnZ2Vi3bh2USmVzPS1qZeLi4lBYWKi1xPLDDz8gPz8f2dnZ2L59O86fP4/NmzcDAM6dO4fly5dj/fr1SEhIgJeXF5KSkgD8OUN/++23sXv3bly5cgUHDhxokeckdgxyM7NkyRLk5+cjNDQUFy9exIQJE7B27VpIpVI8/fTTAIBRo0bh999/x61bt3Qed8iQIWjbti2cnJzQsWNH3Lx501RPgVo5R0dHvPPOO0hMTNRMCA4dOgR/f3/Y29ujTZs2CAkJwcGDBwEADz/8cKP3WvH09ETXrl1hZWWFHj164Pr16832PFoTrpGbke+++w537tzB6NGjERISgpCQEGzduhW5ubn1+gqCoLkvcd2lALW1tY2O3abNX79qiUQCXj5AD2Lw4MFaSyxqtbpen7rX47136vs7vi6NgzNyM2JnZ4f33ntPcwaAIAg4ffo0+vbtixs3buDEiRMAgF27dqFLly7o0KEDOnbsiJKSEgDA3r17NWNZW1vfN9iJHlTdEsvly5fh7e2NvLw8VFVVoba2FtnZ2fD29q63D1+XpsEZuRnx9vbG3LlzMXv2bNTU1AD4c0kkKioKPj4+eOedd1BZWYmHHnoIGRkZAIDIyEjExcUhOztb68yWZ599FrGxsXB2dm6R50KtX90Sy/Tp0zFs2DAoFAqEhISgtrYWgwcPxpQpU+p9E1iPHj1w69YtxMTEYPz48S1UeevDS/SJiESOSytERCLHICciEjkGORGRyDHIiYhEjkFORCRyDHIiIpFjkBMRiRyDnIhI5P4fnGn7R5VLWSQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.boxplot([df07surf['CTDTMP'][iisouth],df07surf['CTDTMP'][iinorth]],\n",
    "            labels=['South','North'],showmeans=True,notch=True);\n",
    "plt.title('CTDTMP - upper 10m 2007')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=-7.375689840798797, pvalue=1.4522086213500351e-11)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ttest_ind(df07surf['CTDTMP'][iinorth],df07surf['CTDTMP'][iisouth],equal_var=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=-9.799130523807385, pvalue=1.7823812194777424e-17)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ttest_ind(df07surf['CTDTMP'][iinorth],df07surf['CTDTMP'][iisouth],equal_var=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>CTDTMP</td>      <th>  R-squared:         </th> <td>   0.286</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.280</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   54.40</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 08 Mar 2022</td> <th>  Prob (F-statistic):</th> <td>1.45e-11</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>14:09:53</td>     <th>  Log-Likelihood:    </th> <td> -310.63</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   138</td>      <th>  AIC:               </th> <td>   625.3</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   136</td>      <th>  BIC:               </th> <td>   631.1</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "           <td></td>              <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>           <td>   13.6931</td> <td>    0.236</td> <td>   57.959</td> <td> 0.000</td> <td>   13.226</td> <td>   14.160</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>is_northern[T.True]</th> <td>   -3.1587</td> <td>    0.428</td> <td>   -7.376</td> <td> 0.000</td> <td>   -4.006</td> <td>   -2.312</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 4.510</td> <th>  Durbin-Watson:     </th> <td>   0.191</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.105</td> <th>  Jarque-Bera (JB):  </th> <td>   4.064</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.408</td> <th>  Prob(JB):          </th> <td>   0.131</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.202</td> <th>  Cond. No.          </th> <td>    2.42</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 CTDTMP   R-squared:                       0.286\n",
       "Model:                            OLS   Adj. R-squared:                  0.280\n",
       "Method:                 Least Squares   F-statistic:                     54.40\n",
       "Date:                Tue, 08 Mar 2022   Prob (F-statistic):           1.45e-11\n",
       "Time:                        14:09:53   Log-Likelihood:                -310.63\n",
       "No. Observations:                 138   AIC:                             625.3\n",
       "Df Residuals:                     136   BIC:                             631.1\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "=======================================================================================\n",
       "                          coef    std err          t      P>|t|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------------\n",
       "Intercept              13.6931      0.236     57.959      0.000      13.226      14.160\n",
       "is_northern[T.True]    -3.1587      0.428     -7.376      0.000      -4.006      -2.312\n",
       "==============================================================================\n",
       "Omnibus:                        4.510   Durbin-Watson:                   0.191\n",
       "Prob(Omnibus):                  0.105   Jarque-Bera (JB):                4.064\n",
       "Skew:                          -0.408   Prob(JB):                        0.131\n",
       "Kurtosis:                       3.202   Cond. No.                         2.42\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = smf.ols(formula=\"CTDTMP ~ is_northern\", data=df07surf).fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing three sample means (another example)\n",
    "\n",
    "We can also compare more than two means. This would be an ANOVA analysis, which can also be expessed as a linear model.\n",
    "\n",
    "##### Create a categorical variable for the region\n",
    "\n",
    "We will divide the data into three categories: north if the Columbia River, between the Columbia River and Golden Gate, and south of the Golden Gate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22      None\n",
       "46      None\n",
       "71      None\n",
       "95      None\n",
       "151     None\n",
       "        ... \n",
       "2275    None\n",
       "2287    None\n",
       "2310    None\n",
       "2334    None\n",
       "2335    None\n",
       "Name: region, Length: 138, dtype: object"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a new variable called \"region\" with no values\n",
    "df07surf = df07surf.assign(region = [None]*len(df07surf))\n",
    "df07surf['region']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22      north\n",
       "46      north\n",
       "71      north\n",
       "95      north\n",
       "151     north\n",
       "        ...  \n",
       "2275    south\n",
       "2287    south\n",
       "2310    south\n",
       "2334    south\n",
       "2335    south\n",
       "Name: region, Length: 138, dtype: object"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# assign string values to region based on latitude\n",
    "northern = (df07surf['LATITUDE'] > 46.2)\n",
    "central = (df07surf['LATITUDE'] <= 46.2) & (df07surf['LATITUDE'] >= 37.8)\n",
    "southern = (df07surf['LATITUDE'] < 37.8) \n",
    "\n",
    "df07surf.loc[northern,'region'] = 'north'\n",
    "df07surf.loc[central,'region'] = 'central'\n",
    "df07surf.loc[southern,'region'] = 'south'\n",
    "df07surf['region']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Box plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'region')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEXCAYAAACkpJNEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArx0lEQVR4nO3dd1gU1/4G8HdpAsJeRUFU1BgVUbFcNXLtIGjsBrsgosaO2LsgSIKVGxMlei1XxRYblqiYm59KNJLYYizYscSgLmIBRIq4e35/5LLXlQ5bhHk/z8PzsLMz53zHlXnnzMzOyIQQAkREJDlGhi6AiIgMgwFARCRRDAAiIoliABARSRQDgIhIohgAREQSZWLoAqh0UCqV2LJlCw4dOgSlUomsrCy4ublh8uTJyMjIgI+PDwAgLS0NCQkJqF27NgCgTZs2cHV1xejRo9XTVCoVKlasiHHjxqFNmzaIi4vD9OnTAQDJycl49eoVHBwcAACenp6Qy+WYO3cu/Pz8MGnSJHVNQgh4eHjAwsIChw8fxr59+xAaGgoHBwfIZDIIIWBhYYHZs2fj73//uz7/ubRu7969OHbsGP71r39pTNu4cSPevn2L1q1bIyAgAKampsXuY+PGjYiMjISxsTFsbGwQEhKCmjVrIj09HQEBAbh+/TpUKhVmzpwJDw8PAMDly5cREhKCtLQ02NnZYfny5bCzs8OXX36J8+fPq9tOSEiAra0tDh06VPx/BNI+QVQIAQEBwt/fX6SkpAghhHj9+rUYP368mDFjhsZ8Z86cET169Chw2o0bN0SbNm3EpUuXNKZHRkaKMWPG5Jjm6uoq3N3dNaafO3dOtGnTRt12bsseP35ctG3bVmRlZRVxjT8ML1++FIGBgaJZs2Ya63br1i3RoUMH8fz5c6FUKsXUqVPFunXrit1PTEyM6Natm3j16pUQQoht27YJLy8vIYQQS5cuFQEBAUIIIR49eiTatWsnnjx5IjIzM0WHDh3EhQsXhBBCbN++XYwaNSpH23/++afo0KGDuH79erHrI93gISAqUHx8PA4dOoRFixbB2toaAGBpaYmFCxeq9wSLysnJCT4+Pti8eXOh5nd0dISlpSUuXryonrZ//3707t073+Vat26NxMREpKSkFLq2+Ph4jRHDu6/37duHMWPGYNSoUejevTtGjBiBhIQEAICPjw8WL16M/v37w93dHStXrlS3cfHiRXh5ecHT0xP9+vVDdHS0ur3s6dmjqHcdPXoUdnZ2mD17tsb048ePo1OnTrCxsYGRkREGDRqE77//Xl3HkiVLMHjwYHTp0gUbNmzAkiVL0LdvX3Tr1g23bt3K0U/lypURHBwMKysrAEDjxo3x+PFjAMCxY8cwYMAAAEC1atXQtm1bHD16FFevXoWVlRVatGgBAOjfvz9+/fVXvHz5UqPtwMBAjBgxAg0aNCjsR0B6wgCgAl27dg1169ZVbxyy2dra4tNPPy12u05OTrh9+3ah5//ss89w8OBBAEB6ejp+++03tG/fPs/5hRDYtWsXHB0dYWNjU+w633f+/HnMnz8fUVFRaNSoEUJDQ9Xv3b9/H9999x3279+PqKgoREdHIzk5GXPnzsWyZcuwf/9+rF69GsHBweoNbFxcHLZu3YqtW7fm6GvIkCGYOHEizMzMNKY/efIEVatWVb+2t7dXBxEAPHr0CDt37sTy5cuxfPlytGrVCvv27UP79u2xbdu2HP04OjqiVatWAIA3b94gLCwMXbt2zbWvKlWqQKFQQKFQwN7eXj3dzMwMNjY2GnWcPHkSjx8/zjXcyPB4DoAKZGRkBJVKpfV2ZTIZzM3NCz1/r1690KdPH8yfPx//93//h06dOsHY2FhjngsXLqBPnz6QyWR48+YNPv74Y409cW1o27at+nzGwIED0adPH/V7gwYNgqmpKUxNTdG1a1ecPn0aRkZGSExMhJ+fn3o+mUym3hOvX79+jnAtiHjvDi5CCBgZ/W9/rnPnzgCAGjVqAIA6KGvWrIlz587l2e6LFy8wadIkWFlZYerUqeq2ZTKZxnzZ/yfeny6E0PhMIiIiMHbs2ByfE30YGABUoCZNmuDevXtITU3V2FAlJCQgMDAQK1euLNKGPNvVq1fh6OhY6PltbW3RsGFDnDp1CgcOHMCcOXNyHG5o2bIl1q5dW2Bbo0ePxtOnTwEAkyZNgru7u/q97BPI2bKysjSWfXdjplKpNF6bmPzvTyp7o6xUKlGnTh3s2bNH/V5CQgJsbGxw6NAhWFpaFljv+6pWraquHwCePn2aY2/8XYU5OXzz5k1MmDABHh4emD17tnq9svuqXLmyui8nJ6ccNWRlZSEpKQlVqlQB8FeYXL58GeHh4UVeP9IPHgKiAlWpUgW9evXCvHnzkJqaCgBITU1FcHAwKlSoUKyN/5UrV/Ddd9/B19e3SMt99tln2LRpE169elWk8Hjf+vXrcfDgQRw8eFBj4w8AcrkcWVlZiIuLAwAcOXJE4/0zZ86oD3Ps3LkTbm5u6ve+//57qFQqJCcn4+jRo+jUqROaNWuGP/74Q31VzI0bN/Dpp59qHCopqk6dOuHEiRN4/vy5+lBXcc/HAIBCoYCvry8mTJiAefPmaYSau7s7du3apZ7v559/hpubG5o2bYqkpCT1eZnIyEg0a9YMcrkcwF/nPRo3blysgCP94AiACiUoKAirV6/G4MGDYWxsjDdv3sDDwwP+/v6FWv7hw4fqQyVGRkawsrJCWFgYnJycilSHh4cHgoKC1IcndMHa2hozZ87E6NGjYWNjoz4Wnq1KlSqYOXMmEhMTUbduXYSEhKjfy8jIQP/+/fH69Wt4eXmhdevWAICVK1di2bJlyMzMhBACy5Ytg4ODQ76HY/Lj5OQEPz8/+Pr6IisrC02bNsXo0aOLvc6rV69Genq6xrkIMzMz7NmzB/7+/ggODkaPHj2gVCoxc+ZM1KxZEwAQHh6OkJAQpKeno0KFCli6dKm6zQcPHqB69erFrol0TybeP5hIRHnat28f/vOf/+R6mMnHxwfe3t45AoPoQ8VDQEREEsURABGRRHEEQEQkUQwAIiKJKjVXAWVkZCA2Nha2trb8UgkRUSEplUokJibC2dk5xyXbpSYAYmNj4e3tbegyiIhKpe3bt6Nly5Ya00pNANja2gL4ayXe/cYjERHlTaFQwNvbW70NfVepCYDswz729vbqe8UTEVHh5HbonCeBiYgkigFARCRRDAAiIoliABARSRQDgIhIohgAREQSxQAgIpIoBgARlXnOzs6QyWTF/nF2djb0KuhEqfkiGBFRccXGxub7/vvPgZYKjgCIiCSKAUBEJFEMACIiiWIAEBFJlE4DIDU1FT179kR8fDwA4PTp0+jduzd69uyJWbNm4c2bN7rsnoiI8qGzALh8+TKGDBmCBw8eqKfNnz8fK1aswOHDh5GRkYGDBw/qqnsiIiqAzgJg9+7dCAoKgp2dnXqaUqlEamoqlEolMjMzUa5cOV11T0REBdDZ9wBCQ0NzTAsODoaPjw+srKzg4OCArl275rpsSkoKUlJSNKYpFAqd1ElEJFV6+yJYYmIiwsLCcPjwYTg4OGDx4sVYvHgxgoKCcswbERGB8PBwfZVGRCRJeguACxcuwNHRETVr1gQADBw4EFOmTMl1Xl9fX3h6empMy36uJRERaYfeAsDR0RFLly7Fs2fPULlyZRw/fhyNGzfOdV65XA65XK6v0oiIJElvAVCnTh1MnjwZw4YNg7GxMWrVqoWQkBB9dU9ERO/ReQCcOHFC/bunp2eOQztERGQY/CYwEZFEMQCIiCSKAUBEJFEMACIiiWIAEBFJFAOAiEiiGABERBLFACAikigGABGRRDEAiIgkigFARCRRDAAiIoliABARSRQDgIhIohgAREQSxQAgIpIoBgARkUQxAIiIJIoBQEQkUQwAIiKJYgAQEUkUA4CISKIYAEREEsUAICKSKAYAEZFEMQCIiCSKAUBEJFE6DYDU1FT07NkT8fHxAIDff/8dAwcORI8ePTBt2jS8efNGl90TEVE+THTV8OXLlxEQEIAHDx4A+CsM/P39sWHDBjg5OWHatGnYu3cvvLy8dFUCEZVxcXFxmDZtGlQqVYnb6tmzZ7GXNTU1xbp162Bra1viOvRJZwGwe/duBAUFYdasWQCAmJgYNGvWDE5OTgCAgIAAKJVKXXVPRBJw+/ZtPH78GMHBwSVq5+rVqxg3blyxl580aRLi4+MZANlCQ0M1Xv/xxx+wtLTE1KlTce/ePTRv3hxz5szJddmUlBSkpKRoTFMoFLoqlYhKMTs7uxLtvQMl2/sHgMDAwBItbyg6C4D3KZVKnD59Grt27UK1atUwf/58rFu3Dv7+/jnmjYiIQHh4uL5KIyKSJL0FQOXKldG0aVPUqFEDANCtWzds27Yt13l9fX3h6empMU2hUMDb21vndRIRSYXeAqBdu3ZYtWoVnjx5gqpVqyI6OhqNGjXKdV65XA65XK6v0oiIJElvAVC1alWEhIRg3LhxyMzMRIMGDTB79mx9dU9EZZC1tTXi4uKQnp4OCwsLg9Tw4sULPHr0CFZWVgbpvyR0HgAnTpxQ/+7q6gpXV1ddd0lEEtGuXTs0b94cU6ZMwdq1a/Xev0qlgq+vL4YOHYp69erpvf+S4jeBiajUkslkWLduHaKjo/M8p6hL//znP5GYmIglS5bovW9t0NshICIiXZDL5dizZw88PDzQokULNGjQQC/9nj59GmFhYTh//jzMzMz00qe2cQRARKVe06ZNsXjxYgwYMACvX7/WeX+JiYkYMmQINm7ciJo1a+q8P11hABBRmfD555+jefPm8PPz02k/KpUKPj4+GDp0KHr06KHTvnSNAUBEZYJMJsOaNWtw/vx5bNq0SWf9LF68GGlpafjiiy901oe+8BwAEZUZ5cuXx549e9CxY0e0aNECTZo00Wr70dHRCA8Px4ULF2BiUvo3nxwBEFGZ0rBhQ3z11VcYMGAAXr16pbV2FQoFhg4dii1btqB69epaa9eQGABEVOb4+PigQ4cOGDNmDIQQJW5PqVTCy8sLo0ePRufOnbVQ4YeBAUBEZdLKlStx/fp1rXxBbOHChTAyMiq1d/3MS+k/iEVElAsLCwusWrUKXbt2xeeffw5TU9NitfP06VOEhobi1q1bMDY21nKVhsUAIKIy6e3bt1iwYAFmz54NxaaZyEr8s9ht/RLQDwEBAfjuu+8gk8m0WKVhMQCIqEwKCgqCqakpAgICSrznXi0jA7dat8aaNWswYcIELVVoeAwAIipzjh49ioiICFy8eFErh23Mzc2xZ88etGnTBi4uLmjRooUWqjQ8ngQmojLlzz//xIgRI7Bjxw7Y2dlprd26deti9erVGDhwIJKSkrTWriExAIiozMjKysKgQYMwZcoUdOjQQevt9+/fHz169MCIESO0cnmpoTEAiKjMmDdvHipWrIhZs2bprI/ly5fj0aNH+Prrr3XWh74wAIioTDh48CB2796NLVu2wMio8Ju2l+nJCDrxFZLSkws1f7ly5bB7924sWbIEZ86cKW65HwQGABGVevfv38fo0aOxa9cuVKpUqUjLRl6Lws3EOOy9FlXoZT766COsX78egwYNwvPnz4ta7geDAUBEpVpmZiYGDRqEuXPn4h//+EeRln2ZnozoB79CQOCnB78WehQAAL1798bAgQMxbNgwqFSqopb9QWAAEFGpNnPmTFSvXh1Tpkwp8rKR16IgxF8bb5VQFWkUAACLFi1CUlISli1bVuS+PwQMACIqtc6fP48DBw5g06ZNRf6Gbvbe/1uVEgDwVqUs8ijA1NQUu3btQlhYGP74448i9f8hYAAQUamVmJgIZ2dnVKhQocjLvrv3n604owAHBwfUqFEDL168KHINhsYAICJJuv38nnrvP9tblRK3n98zUEX6x1tBEJEkLft0vqFLMDiOAIiIJIoBQEQkUQwAIiKJ0mkApKamomfPnoiPj9eYvm3bNvj4+OiyayIiKoDOAuDy5csYMmQIHjx4oDE9Li4O69at01W3RERUSDoLgN27dyMoKEjjftxv3rzBggULMGnSJF11S0REhZTnZaDNmzfPd0EhBIyMjPDbb7/l+n5oaGiOaf/85z/Rr18/ODg45Nt2SkoKUlJSNKYpFIp8lyEioqLJMwCqVq2a76EaIQTGjh1b6I5iYmLw5MkTzJ07F2fPns133oiICISHhxe6bSIiKro8A2Dq1KmoXr16vgtPnTq10B0dPnwYd+7cQZ8+fZCWloZnz55hypQpuT5UwdfXF56enhrTFAoFvL29C90fERHlL88A8PDwQGRkJOrVq4cmTZoAAJYtW4Z69eqpN84eHh6F7mjx4sXq38+ePYvw8PA8n6gjl8shl8sL3TYRERVdnieB9+7di7Vr18LU1FQ9rUWLFlizZg0OHDigj9qIiEiH8hwB7NixA5s3b0a1atXU09zd3eHo6IjJkyfjs88+K1QHJ06cyDHNxcUFLi4uRa+WiIi0Js8RgBBCY+OfrUaNGlAqlbksQVR2OTs7QyaTlejH2dnZ0KtBpCHPEYBSqYRKpcrxcGWVSoW3b9/qvDCiD0lsbGy+78tkMggh9FQNkXbkOQJo1aoVNm/enGP6pk2b0LhxY13WREREepDnCGDy5MkYOnQojh07hubNm0OlUuHSpUtITU3NNRiIiKh0yTMArK2tsWfPHkRFRSE2NhYymQze3t7o0qWLxpVBRERUOuX7RDAzMzN89tlnhb7ih4iISg8+D4CISKIYAEREEsUAICKSqHzPAQDA3LlzNV7LZDJYWFigXr16GDBgAIyNjXVWHBER6U6hRgDXr19H/fr10aBBA8TFxeHx48c4ffo0Fi1apOv6iIhIRwocAdy9exfbt2+HlZUVAGDAgAEYOXIkduzYgZ49e+q8QCIi0o0CRwDJycnqjT8AmJubIzU1FTKZjN8HICIqxQocATRr1gwzZsxA//79IYTAvn370KRJE5w8eRIWFhb6qJGIiHSgwBHAwoULUbVqVSxevBjLly9HjRo1sGDBArx+/RohISH6qJGIiHSgwBGAubk5/Pz80KNHDzg6OiIzMxPlypVD9+7d9VEfERHpSIEBcOnSJUycOBEmJibYuXMn+vTpgzVr1qB58+b6qI+IKF/Pnz9HYmKiQWvIysoyaP/FVWAALFu2DJs3b8aMGTNgb2+PZcuWITQ0FJGRkfqoj4goT69fv8a5c+fQsGHDErXz8uVLVKxYsdjLm5qalmh5QykwADIyMlC3bl31644dO2LFihU6LYqIqDAGDBiglQfxyGQyg48iDKHAk8AmJiZITk6GTCYDANy7d0/nRRERke4VOAIYP348hg4dimfPnmHatGmIiYnh1T9ERGVAgQHg5uaGjz/+GDExMVCpVPDz80OdOnX0URsREelQngHw+PFj9e+mpqZwdXXVeK9atWo6LYyIiHQrzwDo0aMHZDIZhBDIyMhA+fLlYWxsjJSUFFSqVAmnT5/WZ51ERKRleQbA77//DgBYsGABXFxc0KNHDwDA8ePHcezYMf1UR0REOlPgVUCxsbHqjT8AuLu74+bNmzotioiIdK/AAFCpVDh79qz69alTp9SXhBIRUelV4FVAAQEBmDJlCkxNTSGEgBAC3377baEaT01NxeDBg/Gvf/0LDg4O2LVrF7Zu3QqZTAZnZ2csXLgQZmZmJV4JIiIqujwD4ObNm3ByckLLli0RHR2N27dvAwDq168PExMTjXlyc/nyZQQEBODBgwcAgPv37+Pf//439u3bh/Lly2POnDnYsWMHhg8frt01IiqiWbNm4c6dOyVux9PTs0TLe3l5YcCAASWug6iw8gyAuXPnYv/+/QD+ugy0UaNG+c7zvt27dyMoKAizZs0CAJiZmSEoKEj9cBlHR0eNS02JDGXHjh0ICgpC5cqVi93Gb7/9hmHDhhV7+R9++AEnTpxgAJBe5TsCyO+On0IIvHnzJs/3Q0NDNV5Xr14d1atXBwC8ePEC27dvx+LFi3NdNiUlBSkpKRrTFApFnn0RlVS3bt3g4OBQ7OVLuvevUChw5cqVErVBVFR5BoCuLvVMSEjAqFGj0K9fP7i4uOQ6T0REBMLDw3XSPxER/SXPAMjeW9emu3fvYtSoUfDx8cHIkSPznM/X1zfHHpVCoYC3t7fWayIikqoCLwPVltTUVHz++eeYPHlyvht/AJDL5XBwcND4sbe311OlJDUNGjTA3r17Dda/SqXC3r170aBBA4PVQNKktwDYu3cvnj17hk2bNqFPnz7o06cPvvnmG311T5SnNWvW4Msvv8StW7cM1v/r168xYcIEg/RP0lXg9wBK6sSJEwCA4cOH85JP+iDVrVsXwcHBGD58OE6fPg1jY2O99R0XF4egoCDExMSoL68m0he9jQCIPmQTJkyAhYUFwsLC9NanUqnEiBEjMH/+fNSvX19v/RJlYwAQATAyMsLGjRsRFhaG2NhYvfT5zTffQCaTYfLkyXrpj+h9DACi//roo4+waNEi+Pr6IisrS6d93bx5E4sWLcKmTZtgZMQ/QzIM/s8jeseoUaNgZ2eX55cUteHt27fw9fVFSEgIn65HBsUA0BNnZ2fIZLIS/Tg7Oxt6Nco8mUyG9evXIzw8XP1MDG1bvnw5rK2tMW7cOJ20T1RYvOxATwpzXDn7CWxkWA4ODggLC4Ovry/Onz+PcuXKaa3tq1ev4quvvsJvv/3GQz9kcPwfSJQLHx8f1K5dGyEhIVprMysrC76+vliyZAlq1qyptXaJiosBQJQLmUyGtWvXYsOGDTh37pxW2gwNDUXVqlUL/CY8kb7wEBBRHuzt7bFy5Ur4+vri4sWLsLCwKHZbFy9exOrVq3Hp0iU+UY8+GBwBEOVj0KBBKF++PCIjI0vUzpIlS+Dn54dq1appqTKikuMIgCgfP/74I54+fYr2r87hXmjuDz8qjCVNgXsvzuDly5eoWLGiFiskKj4GAFEekpOTMWrUKGzYsAG1unQpcXsr/P0RMWkStm7dqoXqiEqOh4CI8jB16lR0794dXbSw8Qf+Ogx05swZHDhwQCvtEZUUA4AoF4cPH8ZPP/2E5cuXFzjvy/RkBJ34CknpyfnOV758eWzevBkTJkzAs2fPtFUqUbExAIje8+LFC4wdOxYbN26EtbV1gfNHXovCzcQ47L0WVeC8bdu2hZeXF+/9Tx8EBgDRe/z9/dG/f3+4uroWOO/L9GREP/gVAgI/Pfi1wFEAAHzxxRe4evUqdu3apYVqiYqPJ4G1YOPGjZgxY4ZW2rKxsSn2sra2tgZ7qlVZsW/fPpw/fx6XLl0q1PyR16IghAoAoBIq7L0WhVEth+S7jIWFBSIiItCrVy907NiRjzslg2EAaMGjR48wcuRIzJs3r0TttG3bFjExMcVaVqlUckNSQomJifDz80NkZCQsLS0LnD977/+tSgkAeKtS4qcHv6J/o+6oYPG3fJdt1aoVRo0ahbFjx+LAgQP8cpiOOTs749q1a/nOk99n0KhRI709J0KfGABaYmFhUaK9dwC4ceNGsZdVKpUl6lvqhBAYP348fHx80KZNm0It8+7ef7bCjgIAYMGCBfjkk0+wdetWDBs2rFh1U+GUxY23NjAAiADs3LkTN27cwLZt2wq9zO3n99R7/9neqpS4/fxeoZYvV64ctmzZgi5duqBTp05wcHAoUs1EJcUAIMl78uQJpkyZgiNHjsDc3LzQyy37dH6J+27WrBkmTpyIUaNG4ejRozwURHrFq4C05O3bt5LuvzSbOHEixowZg5YtWxqk/7lz5+LZs2eIiIgwSP8kXQwALejUqRM2bNiAEydOGKT/jIwMDBo0CD179jRI/6Xd2bNnMWbMGIP1b2pqipEjR+Ls2bMGq4GkiQGgBW3btsWePXswaNAgHDlyRK99p6WloXfv3jAxMcGePXv02ndZYuhDL4bun6SJAaAlrq6uOHz4MEaOHKm3DXFKSgq6du2KqlWrYufOnTAzM9NLv0RUNjAAtMjFxQU//vgjJk2ahM2bNxdp2cLeTybb8+fP4e7uDmdnZ2zatAkmJjyfT0RFwwDQsqZNmyI6OhqBgYH49ttvC71cUe4no1Ao4OrqCjc3N3z77bd8uDgRFYtOtxypqano2bMn4uPjAQC//PILevXqhS5dumDFihW67NqgnJyccOrUKXz11VdYunRpgfMX5X4yDx8+RIcOHTBw4EAsXbqUx46JqNh0FgCXL1/GkCFD8ODBAwB/Xakyb948rF69GlFRUYiNjcXJkyd11b3B1a5dG6dOncLmzZsRGBgIIUSe8+Z2P5ncxMXFoUOHDhg3bhwCAwO58SeiEtFZAOzevRtBQUGws7MDAFy5cgW1atVCjRo1YGJigl69euGHH37IddmUlBTEx8dr/CgUCl2VqjPVq1fHyZMncfjw4TzvE5TX/WTeHwXcv38frq6umDdvHqZNm6bz2omo7NPZmcPQ0FCN10+fPoWtra36tZ2dHRISEnJdNiIiAuHh4boqTa9MTExgZmaGvub3cS+0X47399taQWVtARj9b29emfUGG7dMwmfPUjXm3TywMZ5YWOi8ZiKSBr1dOqJSqTQOWQgh8jyE4evrC09PT41pCoUC3t7eOq1R2xISEtClSxd4eHig5YKwXNc34T+hUCbFa0xTGsmg+MgRH4/VvNVAxvXrGN6lC16/fo1x48bptHYiKvv0FgD29vZITExUv05MTFQfHnqfXC6HXC7XV2k6ER8fD3d3dwwZMgRBQUF5hl1R7ifTsGFDnDx5Eh4eHkhNTdXaMwiISJr0dv1g06ZNcf/+ffzxxx9QKpU4fPgwOnTooK/u9eru3bto3749Ro8ejeDgYK2erK1Tpw5OnTqF9evXIzg4ON+Ty0RE+dHbCKBcuXJYsmQJ/P39kZmZiY4dO6Jr16766l5vbty4gS5dumDevHkYP368TvqoUaMGTp06hS5duuDVq1cIC8v98BIRUX50HgDv3iCtdevW+P7773XdpcH8/vvv6N69O5YuXarzB3xUqVIF0dHR6NatG8aPH4/Vq1fzC2FEVCTcYmjJr7/+iq5duyI8PFxvT3eysbHBsWPHcOPGDfj6+vKW0ERUJAwALTh37hx69+6NzZs3o1+/nJd66pK1tTWOHj2KZ8+elbqrpIjIsHgHMS34z3/+g5EjR6Jbt24G6d/S0hL79u2DlZWVQfonotKJIwAtMfStmA3dPxGVPhwBaMHFixdx4MABje856JtKpTJY30RUOjEAtGD69OmoX78+PvrooxK1ExISggULFhR7+c6dO5eofyKSFgaAFrRr1w7t2rUrcTvjx4/nLR4M4OnTp7h58yZSU1MLnllHSuPNDqn0YwCQ5GVlZcHX1xfW1tbFbuP+/fuoXbt2ierw8/Mr0fJERcUAIMnTxu00ZDIZbt68qYVqiPSHVwEREUkUA4CISKIYAEREEsUAICKSKAYAEZFEMQCIiCSKAUBEJFEMACIiiWIAEBFJFANAT5ydnSGTyfL9AZDv+87OzgZeCyIqS3grCD2JjY01dAlERBo4AiAikigGABGRRDEAiIgkigFARCRRDAAiIoliABARSRQDgIhIogwSAAcPHkSPHj3Qo0cPLF261BAlEBFJnt4DID09HaGhodi6dSsOHjyICxcu4JdfftF3GUREkqf3bwIrlUqoVCqkp6fD0tISb9++Rbly5TTmSUlJQUpKisY0hUKhzzKJiMo8vQeAlZUVJk+ejG7dusHCwgKffPIJmjdvrjFPREQEwsPD9V0aEZGk6P0Q0M2bNxEZGYno6Gj8/PPPMDIywr///W+NeXx9fXH8+HGNn+3bt+u7VCKiMk3vAXD69Gm0bt0alSpVgpmZGfr27Ytz585pzCOXy+Hg4KDxY29vr+9SiYjKNL0HgJOTE3755RekpaVBCIETJ06gcePG+i6DiEjy9H4OoF27drh+/Tr69u0LU1NTNG7cGGPGjNF3GUREkmeQ5wGMGTOGG30iIgPjN4GJiCSKAUBEJFEMAKJCKOiZzkD+z3PmM53pQ8RnAhMVAp/pTGURRwBERBLFACAikigGABGRRDEAiIgkigFARCRRDAAiIoliABARSVSp+R6AUqkEwCeDEREVRfY2M3sb+q5SEwCJiYkAAG9vbwNXQkRU+iQmJqJWrVoa02RCCGGgeookIyMDsbGxsLW1hbGxsaHL0TqFQgFvb29s376dD78pZfjZlW5l/fNTKpVITEyEs7MzzM3NNd4rNSMAc3NztGzZ0tBl6Jy9vT0cHBwMXQYVAz+70q0sf37v7/ln40lgIiKJYgAQEUkUA4CISKIYAB8IuVyOiRMnQi6XG7oUKiJ+dqWblD+/UnMVEBERaRdHAEREEsUAICKSKAbAB+rVq1fw8/MDAMTHx6NTp04GrogKcuXKFSxfvrxIy6xatQqrVq3SUUVUXNHR0di0aROAsv0ZlZovgklNcnIybty4YegyqAji4uLw/PlzQ5dBWiCVZ0AzAHTo7NmzWLt2LczNzXH37l3Ur18fYWFhOHToEDZt2gSZTIZGjRohMDAQ5cuXxz/+8Q84OzsjMTERtra2ePr0Kfz8/DB37lxkZGRg6tSpuHPnDuRyOb799ltUrFjR0KtYZgghEBYWhmPHjsHY2BiDBg2Cq6srgoODkZSUBHNzcwQGBqJhw4aYM2cOrKyscO3aNSQkJMDPzw+dO3fGypUrkZaWhjVr1qBKlSrYv38/kpKS4Obmhp49e+KLL75AWloaXrx4gTFjxmDIkCGGXu0ySaFQYMaMGUhLS4ORkRECAgIAAKGhocjMzETFihUREhKCWrVqwcfHBxMnToSLiwvi4+MxbNgwrFu3Djt37gQAVKtWDcBfo7vBgwcjISEBffv2hb+/v8HWT6sE6cyZM2dEs2bNxJMnT4RSqRT9+vUTW7ZsER4eHuLFixdCCCGCg4PFkiVLhBBCODo6ijNnzgghhPjzzz+Fm5ub+vf69euLy5cvCyGE8Pf3F9u2bTPAGpVdUVFRYvDgwSIzM1OkpqaK3r17C3d3d3Ht2jUhhBB37twRXbp0EUIIMXv2bOHn5ydUKpW4efOmaNWqlRBCiMjISDF79mz17507dxZZWVlCCCG+/PJL8csvvwghhHj48KFo1qyZEEKIlStXipUrV+p1Xcu6VatWifXr1wshhDh58qRYt26dcHNzU//9REVFib59+wohhBg6dGiuf3Pvfi4rV64Unp6eIjMzUzx//lw0bdpUvHr1St+rpRM8B6Bj9erVg729PYyMjFCnTh0kJyfDzc1Nvfc+aNAgnDlzRj1/06ZNc23Hzs4OTZo0AQDUrVsXL1++1H3xEnL+/Hl069YNZmZmKF++PHbs2AGFQoG5c+eiT58+mD59OtLS0tT/7m3btoVMJoOjoyOSkpJybbNhw4YwMflrkD1nzhxkZmZi7dq1+Prrr5GWlqavVZOc1q1bY+PGjZg+fTqSkpLQsWNHyOVy9d9Pt27d8PDhQ7x69arQbbZv3x5mZmawsbFBxYoVkZycrKvy9YqHgHSsXLly6t9lMhnkcjlSUlLU04QQePv2rfr1+3fry5a9IcluR/DrG1plYmICmUymfv3nn39CCIGDBw+qpykUClSoUAHA/z7Xd5d537uf5ZQpUyCXy+Hm5obu3bvj8OHDWl4DytaiRQscOXIEP/30E6KiorBnz54c8wgh1PfHz/5bevfv8H1l9e+PIwADOHHihHqvcffu3XBxcckxj4mJSb7/IUm7PvnkE/z444/IyspCeno6pkyZAktLS3UAxMTEFPgsCmNj4zw/s5iYGEyaNAkeHh44deoUgNwf0EElt2zZMnz//ffw9PTEggULcPPmTSQlJeHKlSsAgKioKFSrVg0VKlRAxYoVERcXBwA4duyYuo38PsuyhCMAPbOyssLYsWPh4+ODrKwsNGrUCAsXLswxX6VKlVCtWjX4+Phg8eLFBqhUWjp37ozY2Fj07dsXKpUKw4YNg4uLC4KDg7FhwwaYmppixYoV+e7xN2nSBOHh4QgLC8PHH3+s8Z6/vz+8vLxQrlw5ODk5oXr16oiPj9f1akmSj48Ppk+fjn379sHY2BjLly/H3/72N3zxxRdIT0/H3/72N6xYsQIAMGrUKMyZMweRkZFwd3dXt/HJJ59g9uzZqFy5sqFWQy94KwgiIoniISAiIoliABARSRQDgIhIohgAREQSxQAgIpIoBgCRlowePVp9TTlRacDLQImIJIpfBCPCX3duDQ0NhaWlJV6/fo3Jkydj7dq1yMrKgrm5OWbPno2///3vSE9PR1BQEC5fvgxra2vUrVsXALBkyRJ06tQJ33zzDRo3boxdu3Zh69atMDIyQuXKlREYGIjatWur7yR669YtKBQK1K9fH0uXLkX58uUN/C9AUsQAIPqvO3fu4NixY8jKyoK/vz+2bNmCihUr4s6dOxgxYgR+/PFHrFmzBkqlEkePHkVaWhq8vLzQsGFDjXZ+/fVXbNiwAbt27YKNjQ327dsHPz8/HDlyBMBf95rfsmULZDIZBg4ciB9++AH9+vUzxCqTxPEcANF/Va1aFdWrV0dMTAyePn2K4cOHo0+fPpgxYwZkMhkePnyIkydPon///jAyMoKVlRU8PT1ztPPzzz+je/fusLGxAQD07dsXCQkJ6ls/ZN9Z0tTUFI6OjmXmzpJU+nAEQPRflpaWAACVSoXWrVvj66+/Vr/35MkT2NnZwcTERONOkEZGOfehVCpVjmnv3vX13buElqU7S1LpwxEA0Xtat26NmJgY3L17FwBw8uRJ9O7dGxkZGejYsSMiIyOhUqmQnp6Ow4cP57hBXPv27REVFYUXL14AACIjI1GhQgXUqlVL7+tClB+OAIjeU7duXYSEhGDatGkQQsDExARr1qxB+fLlMXbsWISEhKBXr16wtrZGpUqVcjzDoW3bthg+fDh8fX2hUqlgY2ODtWvX5jpaIDIkXgZKVARHjhyBlZUVOnbsCJVKBX9/f7Rt2xZeXl6GLo2oyBgAREVw+/ZtLFiwAOnp6cjKyoKLiwvmzZsHU1NTQ5dGVGQMACIiieJBSSIiiWIAEBFJFAOAiEiiGABERBLFACAikigGABGRRP0/yITYpp50OlQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.boxplot([df07surf['CTDTMP'][df07surf['region']=='north'],\n",
    "             df07surf['CTDTMP'][df07surf['region']=='central'],\n",
    "             df07surf['CTDTMP'][df07surf['region']=='south']],\n",
    "            labels=['north','central','south'],showmeans=True,notch=True);\n",
    "plt.title('CTDTMP - upper 10m 2007')\n",
    "plt.ylabel('[deg C]')\n",
    "plt.xlabel('region')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Analysis\n",
    "\n",
    "Since we have three groups, we can perform a classic Fisher ANOVA on this data set, to test the null hypothesis of no difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "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>Source</th>\n",
       "      <th>ddof1</th>\n",
       "      <th>ddof2</th>\n",
       "      <th>F</th>\n",
       "      <th>p-unc</th>\n",
       "      <th>np2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>region</td>\n",
       "      <td>2</td>\n",
       "      <td>135</td>\n",
       "      <td>86.689</td>\n",
       "      <td>6.087664e-25</td>\n",
       "      <td>0.562</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Source  ddof1  ddof2       F         p-unc    np2\n",
       "0  region      2    135  86.689  6.087664e-25  0.562"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pg.anova(data=df07surf,dv='CTDTMP',between='region')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, looking at the box plot, it appears that the groups do not have equal variances (heteroscedesticity). This suggests that a Welch ANOVA may be more appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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>Source</th>\n",
       "      <th>ddof1</th>\n",
       "      <th>ddof2</th>\n",
       "      <th>F</th>\n",
       "      <th>p-unc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>region</td>\n",
       "      <td>2</td>\n",
       "      <td>79.858</td>\n",
       "      <td>117.841</td>\n",
       "      <td>1.489427e-24</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Source  ddof1   ddof2        F         p-unc\n",
       "0  region      2  79.858  117.841  1.489427e-24"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pg.welch_anova(data=df07surf,dv='CTDTMP',between='region')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both have very low p-values. Now we move on to a post-hoc test. Which post-hoc test to use? The Tukey test is a good choice for the classic ANOVA because it corrects for pairwise comparisons. However, it assumes equal variances. A Games-Howell post-hoc test is preferable in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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>A</th>\n",
       "      <th>B</th>\n",
       "      <th>mean(A)</th>\n",
       "      <th>mean(B)</th>\n",
       "      <th>diff</th>\n",
       "      <th>se</th>\n",
       "      <th>tail</th>\n",
       "      <th>T</th>\n",
       "      <th>df</th>\n",
       "      <th>pval</th>\n",
       "      <th>hedges</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>central</td>\n",
       "      <td>north</td>\n",
       "      <td>10.792</td>\n",
       "      <td>10.124</td>\n",
       "      <td>0.667</td>\n",
       "      <td>0.227</td>\n",
       "      <td>two-sided</td>\n",
       "      <td>2.080</td>\n",
       "      <td>60.999</td>\n",
       "      <td>0.095682</td>\n",
       "      <td>0.535</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>central</td>\n",
       "      <td>south</td>\n",
       "      <td>10.792</td>\n",
       "      <td>14.617</td>\n",
       "      <td>-3.825</td>\n",
       "      <td>0.255</td>\n",
       "      <td>two-sided</td>\n",
       "      <td>-10.588</td>\n",
       "      <td>95.680</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-2.048</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>north</td>\n",
       "      <td>south</td>\n",
       "      <td>10.124</td>\n",
       "      <td>14.617</td>\n",
       "      <td>-4.492</td>\n",
       "      <td>0.210</td>\n",
       "      <td>two-sided</td>\n",
       "      <td>-15.145</td>\n",
       "      <td>88.746</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-3.587</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         A      B  mean(A)  mean(B)   diff     se       tail       T      df  \\\n",
       "0  central  north   10.792   10.124  0.667  0.227  two-sided   2.080  60.999   \n",
       "1  central  south   10.792   14.617 -3.825  0.255  two-sided -10.588  95.680   \n",
       "2    north  south   10.124   14.617 -4.492  0.210  two-sided -15.145  88.746   \n",
       "\n",
       "       pval  hedges  \n",
       "0  0.095682   0.535  \n",
       "1  0.001000  -2.048  \n",
       "2  0.001000  -3.587  "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pg.pairwise_gameshowell(data=df07surf,dv='CTDTMP',between='region')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-values indicate that there is no statistically significant difference between the northern and central regions.\n",
    "\n",
    "Ho does this compare to pairwise t-tests? The t-test indicates that there is a statistically significant difference (barely) at $\\alpha$=0.05 between the northern and central regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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>Contrast</th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>Paired</th>\n",
       "      <th>Parametric</th>\n",
       "      <th>T</th>\n",
       "      <th>dof</th>\n",
       "      <th>Tail</th>\n",
       "      <th>p-unc</th>\n",
       "      <th>BF10</th>\n",
       "      <th>hedges</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>region</td>\n",
       "      <td>north</td>\n",
       "      <td>central</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>-2.080</td>\n",
       "      <td>61.00</td>\n",
       "      <td>two-sided</td>\n",
       "      <td>4.174356e-02</td>\n",
       "      <td>1.575</td>\n",
       "      <td>-0.447</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>region</td>\n",
       "      <td>north</td>\n",
       "      <td>south</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>-15.145</td>\n",
       "      <td>88.75</td>\n",
       "      <td>two-sided</td>\n",
       "      <td>2.478834e-26</td>\n",
       "      <td>4.026e+23</td>\n",
       "      <td>-2.401</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>region</td>\n",
       "      <td>central</td>\n",
       "      <td>south</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>-10.588</td>\n",
       "      <td>95.68</td>\n",
       "      <td>two-sided</td>\n",
       "      <td>8.500978e-18</td>\n",
       "      <td>2.726e+15</td>\n",
       "      <td>-1.945</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Contrast        A        B  Paired  Parametric       T    dof       Tail  \\\n",
       "0   region    north  central   False        True  -2.080  61.00  two-sided   \n",
       "1   region    north    south   False        True -15.145  88.75  two-sided   \n",
       "2   region  central    south   False        True -10.588  95.68  two-sided   \n",
       "\n",
       "          p-unc       BF10  hedges  \n",
       "0  4.174356e-02      1.575  -0.447  \n",
       "1  2.478834e-26  4.026e+23  -2.401  \n",
       "2  8.500978e-18  2.726e+15  -1.945  "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pg.pairwise_ttests(data=df07surf,dv='CTDTMP',between='region')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use statsmodels to perform this test as a linear model. Note that the F-statistic and Prob (F-statistic), a.k.a. p-value, are the same as the classic ANOVA.\n",
    "\n",
    "The regression coefficients in this case can be thought of as pairwise comparisons with the central region. Note that the confidence intervals of the regression coefficient with the north region overlap with zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "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>W</th>\n",
       "      <th>pval</th>\n",
       "      <th>equal_var</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>levene</th>\n",
       "      <td>9.766</td>\n",
       "      <td>0.000109</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            W      pval  equal_var\n",
       "levene  9.766  0.000109      False"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pg.homoscedasticity(data=df07surf,dv='CTDTMP',group='region')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>CTDTMP</td>      <th>  R-squared:         </th> <td>   0.562</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.556</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   86.69</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 08 Mar 2022</td> <th>  Prob (F-statistic):</th> <td>6.09e-25</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>14:09:54</td>     <th>  Log-Likelihood:    </th> <td> -276.85</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   138</td>      <th>  AIC:               </th> <td>   559.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   135</td>      <th>  BIC:               </th> <td>   568.5</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>   10.7920</td> <td>    0.284</td> <td>   37.991</td> <td> 0.000</td> <td>   10.230</td> <td>   11.354</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>region[T.north]</th> <td>   -0.6675</td> <td>    0.474</td> <td>   -1.409</td> <td> 0.161</td> <td>   -1.605</td> <td>    0.270</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>region[T.south]</th> <td>    3.8250</td> <td>    0.354</td> <td>   10.801</td> <td> 0.000</td> <td>    3.125</td> <td>    4.525</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 8.796</td> <th>  Durbin-Watson:     </th> <td>   0.310</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.012</td> <th>  Jarque-Bera (JB):  </th> <td>   4.392</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.204</td> <th>  Prob(JB):          </th> <td>   0.111</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.227</td> <th>  Cond. No.          </th> <td>    4.28</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 CTDTMP   R-squared:                       0.562\n",
       "Model:                            OLS   Adj. R-squared:                  0.556\n",
       "Method:                 Least Squares   F-statistic:                     86.69\n",
       "Date:                Tue, 08 Mar 2022   Prob (F-statistic):           6.09e-25\n",
       "Time:                        14:09:54   Log-Likelihood:                -276.85\n",
       "No. Observations:                 138   AIC:                             559.7\n",
       "Df Residuals:                     135   BIC:                             568.5\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept          10.7920      0.284     37.991      0.000      10.230      11.354\n",
       "region[T.north]    -0.6675      0.474     -1.409      0.161      -1.605       0.270\n",
       "region[T.south]     3.8250      0.354     10.801      0.000       3.125       4.525\n",
       "==============================================================================\n",
       "Omnibus:                        8.796   Durbin-Watson:                   0.310\n",
       "Prob(Omnibus):                  0.012   Jarque-Bera (JB):                4.392\n",
       "Skew:                          -0.204   Prob(JB):                        0.111\n",
       "Kurtosis:                       2.227   Cond. No.                         4.28\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = smf.ols(formula=\"CTDTMP ~ region\", data=df07surf).fit()\n",
    "res.summary()"
   ]
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}