content/11_statistics/linear_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Models - a statistical Swiss Army Knife"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In this section\n",
    "\n",
    "In this section we will learn:\n",
    "\n",
    "- How variants of the familiar equation for a line (y = mx + b) underlie many different statistical models that you may have learned under different names.\n",
    "- Adding an error term to the equation (y = mx + b + e) allows us to incorporate measurement error into our prediction \n",
    "- Why error terms are often assumed to follow a normal (bell-curve) distribution\n",
    "- How expanding the equation for a line lets us handle multiple predictor variables \n",
    "- How categories (like 'red' vs 'blue' or 'country') can be represented or encoded as 'dummy variables' with 0's and 1's, thereby allowing them to be incorporate categorical data as predictors in the equation for a line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Content"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lining it all up - Linear Models, Linear Regression and the Equation for a Line\n",
    "\n",
    "While statistics is incredibly diverse, and can be extremely complex, many statistical methods can be thought of as variations on a simple equation that you have no doubt already been introduced to: the equation for a line `y = mx + b`. \n",
    "\n",
    "In this section, we will discuss many of the statistical tests commonly taught in introductory statistics courses, and show how despite having unique names and histories, many common statistical tests, —ranging from T-tests to ANOVA —  can be seen as rooted in the equation for a line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "In the equation for a line, we predict some variable `y` based on a predictor variable `x` times a slope `m`, offset by an added constant or intercept `b`. The steeper the slope `m`, the greater the effect of changes in our independent variable `x` on `y`. The bigger the intercept, the more the line is shifted upwards.\n",
    "\n",
    "#### A 'perfect' prediction\n",
    "Imagine you were predicting something real using such a line. For example, imagine that you were trying to estimate the amount of DNA that had been replicated by a DNA polymerase during cell division in a bacterium like *Escherichia coli*. \n",
    "\n",
    "You might reasonably predict — in the absence of contrary evidence — that this replication would occur at a fixed rate over time (rather than e.g. increasing exponentially over time or speeding up initially and then slowing down later). If that were right, and if you knew the starting amount of DNA in a cell, and the rate of DNA replication, you could use the equation for a line to predict the amount of DNA that you should expect in a cell after DNA replication had gone on for a certain amount of time. \n",
    "\n",
    "Your equation might look something like this: `dna_replicated = rate_of_dna_replication * time`+ `starting_dna`.\n",
    "\n",
    "If you plotted the predictions of your equation onto a line plot (letting the `x`-axis represent time and `y`-axis represent the amount of DNA, you might get something like this:\n",
    "\n",
    "<img src=\"./resources/line_plot.png\" width=\"400\"  description=\"A cartoon of a line plot with x and y axes. A black regression line runs diagonally up and to the left. \">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A simple linear model in python\n",
    "\n",
    "In statistics, we would refer to this as a simple *linear model* for the amount of dna replicated over time.\n",
    "\n",
    ">**Simple Linear Models** use the equation for a line to predict a response variable based on a constant response to a predictor variable. \n",
    "\n",
    "We can build a simple linear model in python using one function (`linear_simulation`, below) to predict y-values corresponding to some x values using the equation for a line.\n",
    "\n",
    "We can then plot this simple linear model using a `matplotlib` scatterplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAESCAYAAADTx4MfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGqxJREFUeJzt3Xu0ZGV55/HvzwYviA4gDQuBtsG0F5JZgJ4QFHXhXYxLookTHZRWyWpJzFIjOoA6anRIdOIlcWWEaQOKsUEd8YLRRBlEiBlBuwURBeUOLR1oxXs7OOAzf+x9pDic7lOnTl1PfT9r1arab+2qeoo68LD3+z7PTlUhSdJi3WfUAUiSJpMJRJLUExOIJKknJhBJUk9MIJKknphAJEk9MYFIknpiApEk9cQEIknqyU6jDmCQ9txzz1q9evWow5CkibJp06YfVNXKhfZb1glk9erVbNy4cdRhSNJESXJjN/t5CkuS1BMTiCSpJyYQSVJPTCCSpJ6YQCRJPRlZAkmyf5ILklyZ5NtJXt2O/02Sq5JcnuRTSXZrx1cn+WWSy9rbaaOKXZLG1vUb4NOr4az7NPfXbxjYR41yGe+dwAlV9Y0kDwI2JTkPOA84uaruTPJO4GTgxPY111bVISOKV5LG2/Ub4Gvr4K5tzfa2G5ttgAOO6fvHjewIpKq2VNU32sc/A64E9q2qL1bVne1uFwP7jSpGSZoo33zj3clj1l3bmvEBGIs5kCSrgUOBS+Y89XLgnzu2D0hyaZILkzxxO++1LsnGJBu3bt06kHglaSxtu2lx40s08gSSZFfgHOA1VfXTjvE30pzmmj2BtwVYVVWHAq8Fzkry4LnvV1Xrq2qmqmZWrlywEl+Slo9dVi1ufIlGmkCS7EyTPDZU1Sc7xtcCzwGOqaoCqKo7quqH7eNNwLXAI4YftSSNqYNPgRW73HNsxS7N+ACMchVWgNOBK6vqPR3jz6KZNH9uVW3rGF+ZZEX7+EBgDXDdcKOWpDF2wDFw2HrY5WFAmvvD1g9kAh1GuwrrCOAlwLeSXNaOvQF4H3A/4Lwmx3BxVR0PPAl4W5I7gbuA46vq9uGHLUlj7IBjBpYw5hpZAqmqrwCZ56nPb2f/c2hOd0mSxsDIJ9ElSZPJBCJJ6okJRJLGyRBbkSzVsr4ioSRNlCG3Ilkqj0AkaVwMuRXJUplAJGlcDLkVyVKZQCRpXAy5FclSmUAkaVwMuRXJUplAJGlcDLkVyVK5CkuSxskQW5EslUcgkqSemEAkST0xgUhSP01QJflSOQciSf0yYZXkS+URiCT1y4RVki+VCUSS+mXCKsmXygQiSf0yYZXkSzXKa6Lvn+SCJFcm+XaSV7fjeyQ5L8nV7f3u7XiSvC/JNUkuT/KYUcUuSfOasErypRrlEcidwAlV9WjgcOCVSQ4CTgLOr6o1wPntNsBRwJr2tg44dfghS9IOTFgl+VKN8proW4At7eOfJbkS2Bc4Gjiy3e1M4MvAie34h6uqgIuT7JZkn/Z9JGk8TFAl+VKNxRxIktXAocAlwN6zSaG936vdbV/g5o6XbW7H5r7XuiQbk2zcunXrIMOWpKk28gSSZFfgHOA1VfXTHe06z1jda6BqfVXNVNXMypUr+xWmJGmOkSaQJDvTJI8NVfXJdvjWJPu0z+8D3NaObwb273j5fsAtw4pV0pSYokrypRrlKqwApwNXVtV7Op46F1jbPl4LfKZj/Nh2NdbhwE+c/5DUV7OV5NtuBOruSnKTyLxGeQRyBPAS4ClJLmtvzwbeATw9ydXA09ttgM8D1wHXAB8A/mwEMUtazqasknypRrkK6yvMP68B8NR59i/glQMNStJ0m7JK8qUa+SS6JI2NKaskXyoTiCTNmrJK8qUygUjSrCmrJF8qrwciSZ2mqJJ8qTwCkST1xAQiSeqJCUSS1BMTiKTlxVYkQ+MkuqTlY7YVyWw1+WwrEnBifAA8ApG0fNiKZKhMIJKWD1uRDJUJRNLyYSuSoTKBSFo+bEUyVCYQScuHrUiGylVYkpYXW5EMjUcgkqSemEAkST0Z5TXRz0hyW5IrOsY+1nF52xuSXNaOr07yy47nThtV3JIGzEryiTHKOZAPAX8PfHh2oKr+ePZxkncDP+nY/9qqOmRo0UkaPivJJ8rIjkCq6iLg9vmeSxLgPwFnDzUoSaNlJflEGdc5kCcCt1bV1R1jByS5NMmFSZ64vRcmWZdkY5KNW7duHXykkvrHSvKJMq4J5EXc8+hjC7Cqqg4FXgucleTB872wqtZX1UxVzaxcuXIIoUrqGyvJJ8rYJZAkOwHPBz42O1ZVd1TVD9vHm4BrgUeMJkJJA2Ml+UQZuwQCPA24qqo2zw4kWZlkRfv4QGANcN2I4pM0KFaST5SRrcJKcjZwJLBnks3AW6rqdOCF3Hvy/EnA25LcCdwFHF9V807AS5pwVpJPjL4lkCQzwLFV9apu9q+qF21n/KXzjJ0DnLOkACVJfbWkU1hJHprkxCTfBr4GvLI/YUmSxt2ij0CSPIBmknst8GSaJLQReBPwmb5GJ0kaW10nkCRHAscCfwjsClxPkzyOqaqPDiQ6SZPn+g1N4d+2m5rltwef4pzGMrXDBJJkDU3SeDHwMOBm4FRgA/AzmpVQ27b7BpKmi61IpspCRyBXAT8GPgGcVVUXzj6R5GGDDEzSBNpRKxITyLKz0CR62vuiWT4rSdtnK5KpslACeRRwGnAUcGHbYv2UJAcNPjRJE8dWJFNlhwmkqr5XVW8EVgPPAC4EXgV8CziP5sjkgQOOUdKksBXJVOmqDqQa51fVWmBv4GXATTQJ5CNJvp7kjUn+4wBjlTTubEUyVVJVvb842Y+7V2k9iibXrOhTbEs2MzNTGzduHHUYkjRRkmyqqpmF9ltSJXpVba6qv6qqg4DH0SzxlSRNgb71wqqqS4BL+vV+kqTxtlAh4X9e7BtW1Vm9hyNp5KwkV5cWOgL5CM1EeRbYrzruTSDSpLKSXIuwUAI5qov32B14HfAY7k4kkiaRleRahB0mkKr6wvaeS7Ir8BftbTfgs8B/7Wt0kobLSnItwqJXYSV5QJLX03Tj/Uua64D8XlUdXVWX9ztASUNkJbkWoesEkuS+SV4FXAu8E7gCeFJVPauqvr7YD05yRpLbklzRMfbWJN9Pcll7e3bHcycnuSbJd5M8c7GfJ6kLVpJrERZMIElWJHkFcA3wt8ANwNOr6slV9ZUlfPaHgGfNM/7eqjqkvX2+jeEgmmul/3b7mvcnGZuCRWnZsJJci7DQMt61wJtpemFdBvxpVX2uHx9cVRclWd3l7kcDH62qO4Drk1wDHAZ8tR+xSOpwwDEmDHVloVVYH6RZWfUN4NPAwUkO3sH+VVV/vcSY/jzJsTSXyT2hqn4E7Atc3LHP5nbsXpKsA9YBrFrleVtJGpRuKtEDPLa9LaSApSSQU4G3t+/zduDdwMuZvw5l3iXDVbUeWA9NL6wlxCJJ2oGFEsijhxJFq6punX2c5APAP7Wbm4H9O3bdD7hliKFJk8NKcg3JQnUg3x1WIABJ9qmqLe3m82hWegGcC5yV5D3AQ4E1NMuHJXWyklxDtOhmikl2Ah4E/Kyq7uz1g5OcDRwJ7JlkM/AW4Mgkh9CcnroBeAVAVX07yceB7wB3Aq+sKi+xK81lJbmGqKsEkuRA4PU0S2hXdYzfBPwz8O6qunYxH1xVL5pn+PQd7H8K4GJ0aUesJNcQdVMHchTwTZqjgZ1pLmX7yfb+vsDxwGVJ5qvpkDRMVpJriHaYQJLsBWwAfgg8o6r2ayvPX9De70tzVHI7sCHJysGHLGm7rCTXEC10BPInwC7AM6vqf8+3Q1V9EXgmsCtwXH/Dk7QoVpJriBaaA3ka8NmFVmNV1VVJzgWeAbyjX8FJ6oGV5BqShY5AHk337UK+Chy0tHAkSZNioQSyG/CDLt/rB8B/WFo4kqRJsVACuR9N3UU37qJZlSVJmgLd1IGsSfL4LvZ7xFKDkYStSDQxukkgb25vCwleE11aGluRaIIslED+dChRSGrYikQTZKFmiv9zWIFIwlYkmihdXxNd0hDYikQTxAQijRNbkWiCmECkcWIrEk2QRV8PRNKA2YpEE8IjEElST0wgkqSejCyBJDkjyW1JrugY+5skVyW5PMmnkuzWjq9O8sskl7W300YVt7Sg6zfAp1fDWfdp7q/fMOqIpIHY7hxIks/38H5VVb/f5b4fAv4e+HDH2HnAyVV1Z5J3AicDJ7bPXVtVh/QQkzQ8VpJriuxoEn0PBtiapKouSrJ6ztgXOzYvBv5oUJ8vDYSV5Joi200gVXX4MAOZx8uBj3VsH5DkUuCnwJuq6l/ne1GSdcA6gFWrLL7SkFlJrikylpPoSd5I00Z+9uTxFmBVVR0KvBY4K8mD53ttVa2vqpmqmlm50ku0a8isJNcUGbsEkmQt8BzgmKoqgKq6o6p+2D7eBFyL7eM1jqwk1xRZVAJJ8rgkn0xyc5KfJ9k25/aLpQST5Fk0k+bPraptHeMrk6xoHx8IrAGuW8pnSQNhJbmmSNeV6EmeDHwBuBH4IvAy4BPAzsCzgSva8W7f72zgSGDPJJuBt9CsurofcF4SgIur6njgScDbktxJc+XD46vq9m4/SxoqK8k1JdKeJVp4x+QCYE/gd4FdgduAp1XVl5L8DvAVYG1VfWZQwS7WzMxMbdy4cdRhSNJESbKpqmYW2m8xp7AeC3y4qv4v8Ot2bAVAVV0BrAfetNhAJUmTaTEJJMCP2sez8xO7dTz/PeCgfgQlSRp/i0kgNwMHArRHITcDj+t4/mDuTjDS5LIVidSVxbRz/zLwXOAN7fYG4PVJ7kdzKus4mtNY0uSyFYnUtcUkkHcB/5bk/u0RyFuBvYGX0qyMOhs4qd8BSkNlKxKpa10nkKq6jo7ai6r6Fc1Rx3EDiEsaDVuRSF3reg4kyfuTbHdZV5LHJnl/f8KSRsRWJFLXFjOJfjw7bh/yW8ArlhaONGK2IpG61s9eWHsAd/Tx/aThsxWJ1LUdzoEkeTzwhI6h5yTZb55ddwdeDHynj7FJo2ErEqkrC02iP52mRxU0F5d6YXubz63A6/sUlyRpzC2UQN5P0zAxwOU0y3Q/N2efAn4O3FzdNtaSJE28HSaQqtoKbAVIchTwzar692EEJvXs+g1N3ca2m5rVUwef4ikpaQAWUwfyhdnHSXYHVtMcfdxYVbYw0XiwklwamsVeUOqxSS4CfgBsBDYBP0hyUZLHDiJAaVF2VEkuqa8Wc0GpxwAX0cyH/CPNiqsAjwZeAFyY5IlVdekgApW6YiW5NDSL6YX132i67R5RVTd2PpHkzcD/Ad5Ocz1zaTR2WdWctppvXFJfLeYU1uOBU+cmD4Cqugk4jXvWjCwoyRlJbktyRcfYHknOS3J1e797O54k70tyTZLL2yMi6Z6sJJeGZjEJZGfglzt4flu7z2J8CHjWnLGTgPOrag1wPnd3+D0KWNPe1gGnLvKzNA2sJJeGZjGnsC4HXprktKq6xyxlkl1o2rp/czEfXlUXJVk9Z/ho4Mj28Zk01yE5sR3/cFtrcnGS3ZLsU1VbFvOZmgJWkktDsZgE8tfAp4HLkpwKXEWzjPcgmkaLDwf+oA8x7T2bFKpqS5K92vF9aa6COGtzO3aPBJJkHc0RCqtWed5bkgZlMXUg5yZZC7y7vc1WnYdmWe9Lq+qz/Q/xNzJfWPcaqFpPe2XEmZkZK+MlaUAWcwRCVf1jko/RXAt9dTt8A/DV9gJT/XDr7KmpJPsAt7Xjm4H9O/bbD7ilT5+pcWIluTQRdjiJ3q6S+r3Osar6VVVdWFVntrcL+5g8AM4F1raP1wKf6Rg/tl2NdTjwE+c/lqHZSvJtNwJ1dyX59RtGHZmkORZahfVSmrmNgUhyNvBV4JFJNic5DngH8PQkV9N0A35Hu/vnaS6pew3wAeDPBhWXRshKcmliLOoUVr9V1Yu289RT59m3gFcONiKNnJXk0sTo5xUJpaXzmuTSxOjmCGSvJAd2+4ZVdd0S4tG0O/iUe3bTBSvJpTHVTQKZXbbbrRU9xiLdvdrKVVjS2OsmgXyapgpdGg4ryaWJ0E0COaeqzhp4JJKkieIkuiSpJyYQSVJPTCDqv+s3wKdXw1n3ae6tIpeWpR3OgVSVCUaLM9uKZHYZ7mwrEnBiXFpmTBDqL1uRSFPDBKL+shWJNDVMIOovW5FIU8MEov46+JSm9UgnW5FIy5IJRP11wDFw2HrY5WFAmvvD1juBLi1DI23nrmXKViTSVPAIRJLUExOIJKknY5dAkjwyyWUdt58meU2Styb5fsf4s0cd67JlJbmkLozdHEhVfRc4BCDJCuD7wKeAlwHvrap3jTC85c9KckldGrsjkDmeClxbVTeOOpCpYSW5pC6NewJ5IXB2x/afJ7k8yRlJdp/vBUnWJdmYZOPWrVuHE+VyYiW5pC6NbQJJcl/gucD/aodOBR5Oc3prC9u5zG5Vra+qmaqaWbly5VBiXVasJJfUpbFNIMBRwDeq6laAqrq1qu6qql8DHwAOG2l0y5WV5JK6NM4J5EV0nL5Ksk/Hc88Drhh6RNPASnJJXRq7VVgASXYBng68omP4vyc5BCjghjnPqZ+sJJfUhbFMIFW1DXjInLGXjCgcSdI8xvkUliRpjJlAJEk9MYEsR7YikTQEYzkHoiWwFYmkIfEIZLmxFYmkITGBLDe2IpE0JCaQ5cZWJJKGxASy3NiKRNKQmECWG1uRSBoSV2EtR7YikTQEHoFIknpiApEk9cQEMo6sJJc0AZwDGTdWkkuaEB6BjBsrySVNCBPIuLGSXNKEMIGMGyvJJU2IsU0gSW5I8q0klyXZ2I7tkeS8JFe397uPOs6+s5Jc0oQY2wTSenJVHVJVM+32ScD5VbUGOL/dXl6sJJc0ISZtFdbRwJHt4zOBLwMnjiqYgbGSXNIEGOcjkAK+mGRTknYdK3tX1RaA9n6vuS9Ksi7JxiQbt27dOsRwJWm6jPMRyBFVdUuSvYDzklzVzYuqaj2wHmBmZqYGGaAkTbOxPQKpqlva+9uATwGHAbcm2Qegvb9tdBFK0nQbywSS5IFJHjT7GHgGcAVwLrC23W0t8JnRRLgAW5FImgLjegprb+BTSaCJ8ayq+pckXwc+nuQ44CbgBSOMcX62IpE0JcYygVTVdcDB84z/EHjq8CNahB21IjGBSFpGxvIU1kSzFYmkKWEC6TdbkUiaEiaQfrMViaQpYQLpN1uRSJoSYzmJPvFsRSJpCngEIknqiQlEktQTE8h8rCSXpAU5BzKXleSS1BWPQObaUSW5JOk3TCBzWUkuSV0xgcxlJbkkdcUEMpeV5JLUFRPIXFaSS1JXXIU1HyvJJWlBHoFIknpiApEk9cQEIknqiQlEktQTE4gkqSepqlHHMDBJtgI3LuEt9gR+0KdwJpHf3+/v959OD6uqlQvttKwTyFIl2VhVM6OOY1T8/n5/v//0fv9ueApLktQTE4gkqScmkB1bP+oARszvP938/toh50AkST3xCESS1BMTyDySPCvJd5Nck+SkUcczaEn2T3JBkiuTfDvJq9vxPZKcl+Tq9n73Ucc6SElWJLk0yT+12wckuaT9/h9Lct9RxzgoSXZL8okkV7V/B4+bpt8/yV+0f/tXJDk7yf2n6ffvlQlkjiQrgP8BHAUcBLwoyUGjjWrg7gROqKpHA4cDr2y/80nA+VW1Bji/3V7OXg1c2bH9TuC97ff/EXDcSKIajr8D/qWqHgUcTPPPYSp+/yT7Aq8CZqrqd4AVwAuZrt+/JyaQezsMuKaqrquqXwEfBY4ecUwDVVVbquob7eOf0fzHY1+a731mu9uZwB+MJsLBS7If8PvAP7TbAZ4CfKLdZdl+/yQPBp4EnA5QVb+qqh8zRb8/zaUtHpBkJ2AXYAtT8vsvhQnk3vYFbu7Y3tyOTYUkq4FDgUuAvatqCzRJBthrdJEN3N8C/wX4dbv9EODHVXVnu72c/w4OBLYCH2xP4f1DkgcyJb9/VX0feBdwE03i+Amwien5/XtmArm3zDM2FUvVkuwKnAO8pqp+Oup4hiXJc4DbqmpT5/A8uy7Xv4OdgMcAp1bVocAvWKanq+bTzu0cDRwAPBR4IM0p7LmW6+/fMxPIvW0G9u/Y3g+4ZUSxDE2SnWmSx4aq+mQ7fGuSfdrn9wFuG1V8A3YE8NwkN9CcsnwKzRHJbu0pDVjefwebgc1VdUm7/QmahDItv//TgOuramtV/T/gk8DjmZ7fv2cmkHv7OrCmXYFxX5rJtHNHHNNAtef7TweurKr3dDx1LrC2fbwW+MywYxuGqjq5qvarqtU0v/eXquoY4ALgj9rdlvP3/3fg5iSPbIeeCnyHKfn9aU5dHZ5kl/bfhdnvPxW//1JYSDiPJM+m+T/QFcAZVXXKiEMaqCRPAP4V+BZ3zwG8gWYe5OPAKpp/yV5QVbePJMghSXIk8Lqqek6SA2mOSPYALgVeXFV3jDK+QUlyCM0CgvsC1wEvo/kfzKn4/ZP8JfDHNCsSLwX+hGbOYyp+/16ZQCRJPfEUliSpJyYQSVJPTCCSpJ6YQCRJPTGBSJJ6YgKR+izJl5N8edRxSINmApEWkKS6vH1o1LFKw7TTwrtIU+8lc7afDzwPeB1wa8f4te39M4YRlDRqFhJKi5TkrcBbgDVVdc2Iw5FGxlNYUp/NNwfSnuL6SJKnJfl6kl8m+V6SF7TPH5LkS0l+keSWJK/dzns/L8m/Jfl5ezs/yeOH8LWkezGBSMNzMLAB+BxwYjv20TaJfAHYCLyepjvuu5Pc41RYktfQdIq9HTiZ5ihoH+CCJEcM5RtIHTyFJS3SQqewZo8+qurIjrGiaVT5u7NXf0zy28AVNNeZOLqqPtuOP4QmiXy+qv6wHduPpsnh+6rqdR3vuytN59ibquoJ/f6u0o44iS4NzyWzyQOgqr6d5CfAL2aTRzv+wyTfBX6r47XPB3YGzkqy55z3PR84NskuVbVtgPFL92ACkYbnxnnGfkxzGdX5xh/esf2o9n7TPPvOeghgAtHQmECk4blrkeOdl9Wdna98HvDz7ey/tZegpF6ZQKTJMDvXcktVfW2kkUgtV2FJk+ETNFfLe2vHdbp/I8leww9J084jEGkCVNUNSU4A/g7YlOTjwG3AfsCRNCu8njy6CDWNTCDShKiq97Wrs06gqRe5P80E/CXAB0cZm6aTdSCSpJ44ByJJ6okJRJLUExOIJKknJhBJUk9MIJKknphAJEk9MYFIknpiApEk9cQEIknqiQlEktST/w9oJBvFMs6UeAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def linear_simulation(x_data,slope,intercept = 0):\n",
    "    \"\"\"Return y data given x data and parameters\n",
    "    x_data -- a list of x data points (ints or floats)\n",
    "    slope -- the slope of our linear model (m in y=mx+b+error)\n",
    "    intercept -- the intercept of our linear model (b in y=mx+b+error)\n",
    "    \"\"\"\n",
    "    y_predictions = []\n",
    "    for i,x in enumerate(x_data):\n",
    "        y = slope*x + intercept \n",
    "        y_predictions.append(y)\n",
    "        \n",
    "    return y_predictions\n",
    "\n",
    "\n",
    "def plot_scatterplot(x,y,xlabel=\"x\",ylabel=\"y\"):\n",
    "    \"\"\"Plot a simple scatterplot of x vs. y\"\"\"\n",
    "    plt.plot(x,y,\"o\",color = \"orange\")        \n",
    "    plt.xlabel(xlabel,size = \"xx-large\")\n",
    "    plt.ylabel(ylabel,size = \"xx-large\")\n",
    "\n",
    "#set some x data (0-100 in steps of 5)\n",
    "x_data = list(range(0,100,5))    \n",
    "\n",
    "#calculate y data for our x data\n",
    "y_data = linear_simulation(x_data = x_data,slope = 2.0,intercept =35.0)\n",
    "\n",
    "#plot the result\n",
    "plot_scatterplot(x_data,y_data,xlabel = \"Time\",ylabel = \"Total DNA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inferring linear models from data\n",
    "\n",
    "Of course, in real life we often don't know the right slope or intercept for a biological process like DNA replication – we have to infer them from the data.\n",
    "\n",
    "If we are instead inferring the slope and intercept of our line from the data, that is called **linear regression**. If you've ever added a regression line to a scatterplot in Excel, you have encountered linear regression.\n",
    "\n",
    "Several variants of linear regression are commonly used in statistics, including those that have more than one `x` predictor variable (a **multiple linear model**), and those that have more than one response `y` variable (a **multivariate** linear model). We'll talk more about multiple linear models down below.\n",
    "\n",
    "Central to the idea of linear regression is the notion of error. So let's talk about measurement errors, building on our simple example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Virtually all real measurements involve some amount of error\n",
    "\n",
    "Let's say that you had a theoretical estimate for the `rate of replication` and amount of `starting_dna` in your cells, but you weren't totally sure if it was right. You want to compare that estimate against measurents of the *actual* rate of DNA replication over time, to see if the estimate is reasonable.\n",
    "\n",
    "You might set up an experiment using either culture based experiments — in which many bacteria are grown together in a flask or on a plate — or single-cell experiments — in which individual cells are studied using sophisticated instruments — to collect data on the amount of DNA in cells at different points in time. Let's say this all worked and you got some data back.\n",
    "\n",
    "*If* DNA replication was indeed happening at a linear rate, and *if* the rate of replication you estimated in that experiment were correct, then after you got your measurements back, do you think the points would fall precisely on the line predicted by the equation?\n",
    "\n",
    "Of course not! If your experimental setup, instrumentation, and execution were outstanding, then they might fall *extremely close* to the line predicted by your equation. But they would not fall right on the line with perfect mathematical precision. Instead, if we plotted time on the x-axis and the amount of DNA replicated on the y-axis, they might look something like the diagram below:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./resources/linear_scatter_high_r2_linear_regression.png\" width=\"400\"  description=\"A cartoon of a scatter plot with x and y axes. A series of orange points is plotted. The points have some scatter to them, but roughly form a diagonal line from the lower left to the upper right of the plot. A black regression line runs diagonally through the points. \">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The process of measuring the rate of replication will inevitably be at least somewhat noisy, due to a variety of causes that might include:\n",
    "\n",
    "- imprecisions of the instruments used to measure DNA replication\n",
    "- imprecisions in when all the cells began replicating\n",
    "- imprecisions in when each culture was measured\n",
    "- etc\n",
    "\n",
    "These types of imprecisions in measurement together constitute the amount of \"error\" that separates the *measured* number from our prediction if our model is in fact 100% correct.  The word \"error\" can sound a little scary, but some amount of error is virtually inevitable in any measurement. Small errors may be quite harmless, while larger errors may make it look like the data do not fit our linear estimate (even if they do). \n",
    "\n",
    "Even if our model is 100% correct, the presence of measurement error will inevitably mean our actual results are at best slightly blurred versions of that biological truth. Our job is not to pretend that measurement error doesn't exist, but rather to quantify it as best we can. In real situations, we will often need to decide whether differences between what a given model predicts and what we observe are likely to be due to measurement error, or if the data contradict the model.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Try it out**: try rerunning the code several times without changing any parameters. Does the shape of the distribution change, or stay exactly the same?\n",
    "\n",
    "**Try it out**: try increasing the `scale` parameter, how does the shape of the normal distibution change? Do you think a higher `scale` parameter corresponds to more observation error or less?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Residuals can be used to quantify error\n",
    "\n",
    "The image below shows one way to start thinking about how we could estimate error. The grey lines show how different the y-value for each observation (orange point) is from the y-value we predicted using our simple linear model (black line). \n",
    "\n",
    "<img src=\"./resources/linear_scatter_high_r2_linear_regression_residuals.png\" width=\"400\"  description=\"A cartoon of a scatter plot with x and y axes, with orange points and a black regression line just as in the previous figure. Now, however, thin grey lines stretch vertically from each point to the regression line. These are the residuals — the difference between what the regression model predicts and the actual observed data. \">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These differences are known as the 'residuals', short for 'residual variation' in our `y` values, after accounting for the variation in `y` values predicted by our model.\n",
    "\n",
    "<img src=\"./resources/labelled_residuals.png\" width=\"400\"  description=\"A cartoon of a scatter plot with x and y axes, with orange points and a black regression line just as in the previous figure. Now, however, thin grey lines stretch vertically from each point to the regression line. These are the residuals — the difference between what the regression model predicts and the actual observed data.\">\n",
    "\n",
    "> **Residuals** measure the difference between model predictions and observations. They represent *residual* variation that isn't explained by the model, but might be explained by other factors not included in the model.\n",
    "\n",
    "If we wanted to simulate this whole process, and incorporate error into our simple linear model, we might write an equation that looks like this:\n",
    "\n",
    "y = mx + b + e\n",
    "\n",
    "The `e` parameter would not be a fixed value across the whole dataset, but rather represent the random error in each *individual* measurement. Those random errors  might be either positive or negative. Adding error to each y value will cause data points to fall either above or below the line representing the predictions of our model. \n",
    "\n",
    "If we wanted to polish this up to make this look a bit more 'mathy' — that is, to use conventions similar to those used in slightly more complex models — we could represent all the parameters of the model itself as &beta;s, and that funky observation-specific error term as &epsilon;. This would produce an equation that looks like the one below:\n",
    "\n",
    "y = &beta;<sub>0</sub>x<sub>0</sub> + &beta;<sub>1</sub> + &epsilon;\n",
    "\n",
    "This is often rearranged so the constant goes first:\n",
    "\n",
    "y = &beta;<sub>0</sub> + &beta;<sub>1</sub>x<sub>1</sub>  + &epsilon;\n",
    "\n",
    "In this equation &beta;<sub>0</sub> is the intercept, and &beta;<sub>1</sub> is the coefficient or slope saying how `y` changes with increasing values of x<sub>1</sub>.\n",
    "\n",
    "A final way you will see simple linear models represented in programming languages like `R` or in the `patsy` python module is in a shorthand that removes the &beta; coefficients (slopes) for each term, and leaves the error term implicit: \n",
    "\n",
    "`y ~ 1 + x`\n",
    "\n",
    "This notation can be a little confusing at first, but is very quick to write. Note that if we multiply each term by it's own beta, and add on the implicit error term, the shorthand notation `y ~ 1 + x`  becomes y ~ &beta;<sub>0</sub> + &beta;<sub>1</sub>x<sub>1</sub> + &epsilon; just like we saw above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating normally-distributed errors in python\n",
    "\n",
    "We've established that even if we our simple linear model were a perfect description of biology, we would still have some errors during measurement, due to imprecision, less than infinite sampling, etc.\n",
    "\n",
    "So let's see how our model looks with different levels of random error. First we'll simulate drawing random numbers from a normal distribution, then add this in to our plot. You may wonder why we want random numbers from a normal distribution, rather than 'regular' random numbers that are equally likely. This is due to a very important statistical result known as the Central Limit Theorem, which tells us that the sum of many errors is normally distributed (even if the individual errors are not), so a **normal distribution** (aka Bell Curve, Gaussian distribution, Gauss-Laplace distribution) is useful to represent the sum of many small errors. We'll discuss the Central Limit Theorem further below.\n",
    "\n",
    "Before we add random error to our linear model, let's simulate normally distributed random data on it's own, just so we can see what this type of data looks like.\n",
    "\n",
    "The `np.random.normal` function (i.e. the `normal` function from the `numpy` package's `random` module), generates random numbers by drawing from a normal distribution. It's `loc` parameter set's the **mean** value for the normal distribution. It's `scale` parameter sets the standard deviation for the data, which controls how 'flat' or 'spread out' the normal distribution will be. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Count')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAESCAYAAADTx4MfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGsRJREFUeJzt3XucZGV95/HPNwh4CyIyXIQhPSquoonKTrzuGiIaAS9glASDBpQEjUi8xCiaVVxXd72L7ipmIgquiLJeAlGiIorGRI2DV646IpcRHMagKCrI5bd/nNNaFtUzXaeru6q6P+/X67yq6jlPnfM7Xd316/M8z3lOqgpJkob1W+MOQJI0nUwgkqROTCCSpE5MIJKkTkwgkqROTCCSpE5MIJKkTkwgkqROTCCSpE5uN+4AFtPOO+9cMzMz4w5DkqbKeeed98OqWrW1ess6gczMzLB+/fpxhyFJUyXJ5fOpZxOWJKkTE4gkqRMTiCSpExOIJKkTE4gkqRMTiCSpExOIJKmTsSWQJO9Ock2S83vK3pDk4iTfTPLRJDv2rHtpkg1JLkny2PFELUmaNc4zkJOBA/rKzgbuX1W/B3wbeClAkn2Aw4D7te95R5Jtli5USVK/sSWQqvo8cG1f2aeq6ub25ZeAPdvnBwMfqKobq+p7wAbgwUsWrJatmdW7kaTTMrN6t3GHL43VJE9l8kzgg+3zPWgSyqyNbdltJDkaOBpgr732Wsz4tAxcvnETdWq39+bwTaMNRpoyE9mJnuTvgJuB2T/tDKhWg95bVeuqam1VrV21aqtzgUmSOpq4M5AkRwCPB/avqtkksRFY3VNtT+CqpY5NkvRrE3UGkuQA4CXAE6vq5z2rzgQOS7J9kjXA3sC/jyNGSVJjbGcgSU4D9gN2TrIROJ5m1NX2wNlJAL5UVc+uqguSnA5cSNO0dUxV3TKeyCVJMMYEUlVPHVB80hbqvwZ4zeJFJEkaxkQ1YUmSpocJRJLUiQlEktSJCUSS1IkJRJLUiQlE6mj7bek8j5ZzaWk5mLgr0aVpceNNdJ5HC5xLS9PPMxBJUicmEElSJyYQSVInJhBJUicmEElSJyYQSVInJhBpTBZyHYnXkGgSeB2INCYLuY7Ea0g0CTwDkSR1YgKRJHViApEkdWICkSR1YgKRJHViApEkdWICkSR1YgKRppAXIWoSeCGhNIW8CFGTwDMQSVInY0sgSd6d5Jok5/eU7ZTk7CTfaR/v2pYnyduSbEjyzST7jituSVJjnGcgJwMH9JUdB5xTVXsD57SvAQ4E9m6Xo4ETlyhGSdIcxpZAqurzwLV9xQcDp7TPTwEO6Sl/bzW+BOyYZPeliVSSNMik9YHsWlVXA7SPu7TlewBX9tTb2JZJksZk0hLIXDKgrAZWTI5Osj7J+s2bNy9yWJK0ck1aAtk02zTVPl7Tlm8EVvfU2xO4atAGqmpdVa2tqrWrVq1a1GAlaSWbtARyJnBE+/wI4Iye8j9vR2M9FLhutqlLkjQeY7uQMMlpwH7Azkk2AscDrwVOT3IUcAVwaFv9LOAgYAPwc+AZSx6wJOk3jC2BVNVT51i1/4C6BRyzuBFJkoYxaU1YkqQpYQKRJHViApEkdWICkSR1YgKRJHViApEkdWICkSR1YgKRJHViApEkdWICkSR1YgLR1JtZvRtJOi2SuhvbXFjSqFy+cRN1arf35vDRxiKtJJ6BSJI6MYFIkjoxgUiSOjGBSJI6MYFIkjoxgUiSOjGBSJI6MYFIkjoxgUiSOjGBSJI6MYFIkjoxgUiSOjGBSJI6mcgEkuQFSS5Icn6S05LcPsmaJF9O8p0kH0yy3bjj1Og4Jbs0fSZuOvckewB/DexTVb9IcjpwGHAQ8Jaq+kCSdwJHASeOMVSNkFOyS9NnIs9AaBLbHZLcDrgjcDXwKOBD7fpTgEPGFJskiQlMIFX1feCNwBU0ieM64Dzgx1V1c1ttI7DHeCKUJMEEJpAkdwUOBtYAdwfuBBw4oGrN8f6jk6xPsn7z5s2LF6gkrXATl0CARwPfq6rNVXUT8BHg4cCObZMWwJ7AVYPeXFXrqmptVa1dtWrV0kQsSSvQJCaQK4CHJrljmiE2+wMXAp8FntLWOQI4Y0zxSZKYwARSVV+m6Sz/KvAtmhjXAS8BXphkA3A34KSxBSlJmrxhvABVdTxwfF/xpcCDxxCOJGmAoc5AktyS5M+2sP5Pk9yy8LAkSZNu2CasrV32+1vMMTpKkrS8dOkD2VKCeAjw446xSJKmyFb7QJI8D3heT9EJSV4zoOqOwF2A00YUmyRpgs2nE/0nwPfb5zPAj4D+K/QKuBj4Cs1V5JKkZW6rCaSq3gO8ByDJ94AXV9WZix2YJGmyDTWMt6rWLFYgkqTp0vk6kCR3AnZiwMisqrpiIUFJkibfUAkkybbAfwOeBWxpoqltFhKUJGnyDXsG8jbgaOBjNHNT/WjkEUmSpsKwCeRPgPdW1TMWIxhJi2/7bel8K+Df2XNXLrvyByOOSNNq2ASyLfDFxQhE0tK48SYWcPvgTaMNRlNt2CvRPwOsXYxAJEnTZdgEciywX5Jjk2y/GAFJkqbDsAnkX4G7AicA1yf5fpIr+pbLRx+mpEkw23/SZZlZvdu4w9eIDdsHcinOtiutWPafqNewV6Lvt0hxSJKmzMTd0laSNB1MIJKkToadyuRW5tEHUlVOZSJJy9ywneiv4rYJZBvgnsDBwIXAx0cQlyRpwg3bif7KudYlmQG+BJy/oIgkSVNhZH0gVXUZcCJw/Ki2KUmaXKPuRL8G2HvE25QkTaCRJZB2apPDgatGtU1J0uQadhTWu+dYdVfgYcAuwF8tNKgkOwLvAu5P02n/TOAS4IPADHAZ8CdV5f1IJGlMhh2F9ShuOwqraG4sdS5wYlV9bgRxvRX4RFU9Jcl2wB2BlwHnVNVrkxwHHAe8ZAT7kiR1MOworJlFiuNXkuwAPBI4st3nL4FfJjkY2K+tdgpNwjKBSNKYTOKV6PcANgPvSfK1JO9Kcidg16q6GqB93GWcQUrSSjdsExYASe4HPAFYQ9OEdRnwT1V1wYhi2hc4tqq+nOStNM1V843taJr7trPXXnuNIBxJ0iBDJ5AkbwOOAfpvqvyaJP+nqp63wJg2Ahur6svt6w/RJJBNSXavqquT7E4zZPg2qmodsA5g7dq1Tj0vSYtkqCasJM8Hngv8I/BwYMd2eTjwUeC5SRaUQKrqB8CVSf5TW7Q/zRQpZwJHtGVHAGcsZD+SpIUZ9gzkL2lGRz25r/xLwFOSnEXTfPTWBcZ1LHBqOwLrUuAZNMnu9CRHAVcAhy5wH5KkBRg2gdwTePsW1n8MeFP3cBpV9XVg7YBV+y9025Kk0Rh2FNaPaTrO53IP4Cfdw5EkTYthE8gnaPo5Dulf0V6n8RzgrFEEJkmabMM2Yb0M+EPgw0k2ABfTDOO9L3Av4Mq2jiRpmRvqDKSqrgIeRNPPcSvwR8Bj2+dvAvadvdhPkrS8DX0dSFVdC7y4XSRJK9S8zkCS3D3JHvOos/towpIkTbqtJpAkD6S57uJpW6n6NOCKdpoTSdIyN58zkOfQzHX1hq3UexPwPZppTiRJy9x8EsgfAv+vqm7dUqWquoVm3qpHjyIwSdJkm08C2RP4zjy3twFY3T0cSdK0mE8C+SVw+3lu7/bATd3DkSRNi/kkkO/SzLY7Hw+nOQvRCjSzejeSdFokTZ/5XAdyBvCyJG+oqm/MVakdrXUo8D9HFZymy+UbN1GndntvDh9tLJIW33zOQE4ANgGfTvLnSbbtXZlkuyRHAme39U4YeZSSpImz1QRSVdcBB9LMsvse4Lr2XuWfS/I1mhl6TwJ+ChzU1pckLXPzmsqkqi5I8nvAs4BDgH2AHWiSynqaOxSuq6rrFytQSdJkmfdcWFX1M+DN7SJJWuGGvR+IJEmACUSS1JEJRJLUiQlEktSJCUSS1IkJRJLUiQlEktSJCUS/spDJEJ0QUVp55n0hoZa/hUyGCE6IKK00E3sGkmSbds6tj7Wv1yT5cpLvJPlgku3GHaMkrWQTm0CA5wEX9bx+HfCWqtob+BFw1FiikiQBE5pAkuwJPA54V/s6wKNo7rkOcArNpI6SpDGZyARCc0+RFwO3tq/vBvy4qm5uX28E9hj0xiRHJ1mfZP3mzZsXP1JJ87L9tixokMbM6t3GfQjqM3Gd6EkeD1xTVecl2W+2eEDVGvT+qloHrANYu3btwDqSlt6NN7HAQRqbRheMRmLiEgjwCOCJSQ4Cbk9z35ETgB2T3K49C9kTuGqMMUrSijdxTVhV9dKq2rOqZoDDgM9U1eHAZ4GntNWOoLlXuyRpTCYugWzBS4AXJtlA0ydy0pjjkaQVbRKbsH6lqs4Fzm2fXwo8eJzxSJJ+bZrOQCRJE8QEIknqxAQiSerEBCJJ6sQEIknqxAQiSerEBCJJ6sQEIknqxAQiSerEBCJJ6sQEIknqxAQiSerEBCJJ6sQEIknqxAQiSerEBCJJ6sQEIknqxAQiSerEBLLMzKzejSSdFkkaxkTfE13Du3zjJurUbu/N4aONRdLy5hmIJKkTE4gkqRMTiKSpsP22dO7fm1m927jDX5bsA5E0FW68iQX0720abTACPAORJHU0cQkkyeokn01yUZILkjyvLd8pydlJvtM+3nXcsUrSSjZxCQS4Gfibqrov8FDgmCT7AMcB51TV3sA57WtJ0phMXAKpqqur6qvt858CFwF7AAcDp7TVTgEOGU+EkiSYwATSK8kM8CDgy8CuVXU1NEkG2GV8kUmSJjaBJLkz8GHg+VX1kyHed3SS9UnWb968efEClKQVbiITSJJtaZLHqVX1kbZ4U5Ld2/W7A9cMem9VrauqtVW1dtWqVUsTsCStQBOXQNLM6ncScFFVvbln1ZnAEe3zI4Azljo2SdKvTeKFhI8Ang58K8nX27KXAa8FTk9yFHAFcOiY4pMkMYEJpKq+AMw1t/j+SxmLJGluE9eEJUmaDiYQScueEzEujolrwpKkUXMixsXhGcgE8ra0kqaBZyATyNvSSpoGnoFIkjoxgUiSOjGBSJI6MYFIkjoxgUiSOjGBSJI6MYFIkjoxgUiSOjGBSJI6MYFIkjoxgUjSFjiT79ycC0uStsCZfOfmGcgicUZdScudZyCLxBl1JS13noFIkjoxgUiSOjGBSJI6MYFI0iJZyBDgaRgGbCe6JC2ShQwBhskfBuwZiCRNqEm/iHHqzkCSHAC8FdgGeFdVvXYx9jOzejcu3zjZ2V/S8jbpFzFOVQJJsg3wduAxwEbgK0nOrKoLR72vhVzHAV7LIWn5m7YmrAcDG6rq0qr6JfAB4OAxxyRJK9K0JZA9gCt7Xm9syyRJSyxVNe4Y5i3JocBjq+ov2tdPBx5cVcf21DkaOLp9eX/g/CUPdHHtDPxw3EGMmMc0HTym6TCKY/qdqlq1tUpT1QdCc8axuuf1nsBVvRWqah2wDiDJ+qpau3ThLT6PaTp4TNPBY1qYaWvC+gqwd5I1SbYDDgPOHHNMkrQiTdUZSFXdnOS5wCdphvG+u6ouGHNYkrQiTVUCAaiqs4Cz5ll93WLGMiYe03TwmKaDx7QAU9WJLkmaHNPWByJJmhDLLoEkOTTJBUluTbK2p3wmyS+SfL1d3jnOOIcx1zG1616aZEOSS5I8dlwxLlSSVyb5fs/nc9C4Y+oqyQHt57EhyXHjjmcUklyW5FvtZ7N+3PF0keTdSa5Jcn5P2U5Jzk7ynfbxruOMcVhzHNOS/S0tuwRCc93HHwOfH7Duu1X1wHZ59hLHtRADjynJPjQj0e4HHAC8o53uZVq9pefzmW8/10TpmW7nQGAf4Knt57Qc/GH72UzrsNeTaf5Oeh0HnFNVewPntK+nycnc9phgif6Wll0CqaqLquqScccxSls4poOBD1TVjVX1PWADzXQvGh+n25lQVfV54Nq+4oOBU9rnpwCHLGlQCzTHMS2ZZZdAtmJNkq8l+VyS/zruYEZguU3t8twk32xPy6eqKaHHcvtMZhXwqSTntbM9LBe7VtXVAO3jLmOOZ1SW5G9pKhNIkk8nOX/AsqX/9K4G9qqqBwEvBN6fZIeliXjrOh5TBpRN7LC6rRzjicA9gQfSfFZvGmuw3U3VZzKER1TVvjRNc8ckeeS4A9KcluxvaequAwGoqkd3eM+NwI3t8/OSfBe4NzARHYJdjol5TO0ySeZ7jEn+AfjYIoezWKbqM5mvqrqqfbwmyUdpmuoG9TNOm01Jdq+qq5PsDlwz7oAWqqp+dSOQxf5bmsozkC6SrJrtYE5yD2Bv4NLxRrVgZwKHJdk+yRqaY/r3McfUSfvHO+tJTO8kmMtuup0kd0ry27PPgT9iej+ffmcCR7TPjwDOGGMsI7GUf0tTeQayJUmeBPxvYBXw8SRfr6rHAo8EXpXkZuAW4NlVNbbOp2HMdUxVdUGS04ELgZuBY6rqlnHGugCvT/JAmuaey4BnjTecbpbpdDu7Ah9NAs13xvur6hPjDWl4SU4D9gN2TrIROB54LXB6kqOAK4BDxxfh8OY4pv2W6m/JK9ElSZ2smCYsSdJomUAkSZ2YQCRJnZhAJEmdmEAkSZ2YQDRWSc5thx+uWO1M0ZXkyHHHIg3DBKKRSrJLkjcmubidPv/HSf4lyTOTrOjftyQvnMQk0ZPA5lq+MO4YNZmW3YWEGp8k+wKfAH6bZmbT84A70cx4ehLw5CR/3E4rsxK9kGbG5JP7yi8H7gDctNQB9TkD+NCA8qmf3kOLwwSikWgnpvxHmquvH1JV3+xZfUKSvwNeTXPl7wvGEOK8JLlDVf1iKfdZzdW8NyzlPudwflW9b9g3tVO23DJoFoQ0l69vX1ULOr52GqJt2inyNSFWdJOCRupomkkEX9KXPACoqtcA/0Yzk+vu/euT3DvJJ5Ncn2Rzkre38y711rlHkve3d1u7MclVST6e5AF99dYkeW+STW29S5L8bX8TWts8874kj0uyPskNwEuSfDTJte0XY3+cr2jfd+/29e8mOSnNHe1+3jbZfSLJ7/fvi2Za9z/oaRq6rF03sA8kyR5JTu45jguTvKD9Uu6td26SjUn2SvKPSX7axv/OJNv3H8NCJDmyjfWgJP+r7b+6AVjdcxyvbutdQDOB6WE97396kq+2zZvXJvlwkvv07WO/djt/meRFaSY+vRF4+CiPRQvnGYhG5RCaP/JTt1Dn3TRfAge2z2fdAfg08DngxcDDgOcAa4CDAJJsC3wKuDPNdNVX0szR9EjgvsA32nr3Ar4IXE8zf9hmmrmCXg/MAMf0xfSf232cCPx9u93z2+M5kNtOrvdUYH1Vfbt9/VjgAcAHaZqidgWOAj6XZN+qurit93TgbcAm4DVt2fVz/aCS3I0m4e5Gc4fDS4HHA2+mmar7uX1vuQPNHfU+B/wtzc/wWe3xv3yu/fRvI8nOA8p/XlU/7yt7PfBz4I3Atu2x3LlddzBwN5qf6Wbg4vaYXgS8gWbCz5cCOwHHAl9M8vtVtaFvH8+j+Y5a1+7r6nkeh5ZKVbm4LHihuSvaN7ZSZ1+aCd7e2FN2blv2ur66b2rLD2xfP6B9fehW9nEWzZftDn3lbwZuBfbuKat2eVhf3dsD1wEfnCP+5/eU3WlADDvT9Bu8s698I3DugPoz7XaP7Cl7fVv25J6yAB9py393wM/w+X3bPRO4Zh6f3ez+51pe2VP3yLbsG8B2c2znBuB3+tbdDfgFze0Ttu/7md4CfKinbL92O5uAu4z7d9tl7sUmLI3KDsBPtlJndv1dBqw7oe/17E1wntA+Xtc+HpDkzgyQZEea+0N/GNguyc6zC83suAEe1fe2r1XVF3sLqmmv/zDwhL59/RnNl90Heur+rGf/d2zPHKD5L/s3mrGG9ESaW+N+uGdfRfMfPPz65zLrVpozqF6fBValnYp9Hk4FHjNgee+AuifV3P0RZ1XV5X1lj6FJzCdUzyCKqvoqzdnnQUn6W0ROrarr0MQygWhUfkKTRLZkdn1/ovlptbcVnVXNDYx+StOMRVVdRpNUngn8MMlnkhyXpPfmTfemSRIvomk66V1mpx/vv2Xp9+aI9f00zUJPgl91Bv8p8Jmq+sFspSQ7JHlbkquBnwE/bPf3OGAhtxKdoW366XNh+7imr/yaum3n/4/ax53muc/LqurTA5ZB982Z6+c217qZ9vGiAesupPlZ7zrEPjQB7APRqFwIrN3KKKYHto/998eY654Cv9FZXFUvSvIumv++Hw28Enh5OzT4k/z6H6J30pxBDNL/ZThXrJ+haXP/M+D/An9Ac3fB/v6E04D9gbcAX6U5U7qVpo3/nnNse762dK+F/nVbug/MoNvsLtSWRqptad2gY8oc65Z0NJyGZwLRqJwJPILmC/ekOeo8g+Zah3/uK98h7W1FZwuS3J2mU/Y3/gutplP6YuAN7dnH12huovNJ4Lu07fZV9emFHExV3ZrkA8CxSVa1x3UDTR/EbIw70nTA//eqemXv+5P8j0GbHSKEy2gGB/S7b8/6aXJZ+7gPTaLtdR+aZOH1JlPGJiyNyt8D3wdel+R+/SuTHAf8F+Ad/c1Vref3vf6b9vFj7ft36G8jr6oraZqLdmpfb6YZiXREOxqrP4YdhhzWeirNP1mHA08G/qmqepvfbqFJCv3Dg/cDHjJge9cz/2atfwLuleZulLPbnW2em10/Tc6mScB/3Ts8uh2C/Rjgn6vq5nEFp248A9FIVNV17ZfdWcBXkvReif5Ems7rs4DjBrz9R8BT01wf8iWaIahPAz5VVWe1dR4FnJjkQ8AlNLfwfTzNf6+v6NnWXwH/Cnytbe66iKbT/v40SeD+zPO/96o6L8klwKtorq4/tW/9T5OcA7w4yR2AbwO/S3OmdUH7nl5fAZ6W5Pi27vVVNVcieC3wJ8BpSWaH8T6OZmjx26vqW/M5hiHdP8nTBpTfUFWDrlCft6r6jySvoBld9vk0t2KdHcb7Ewb/XmjSjXsYmMvyWmiuW3gzzRfkDTR9Al+guTZimwH1z6UZ3npvmmao64H/oLmG4M499dYA/9Bu92ftdr/Sbjd929yjff8VwC9phoP+C831EbfvqVfA+7ZyPC9v611L37DVdv0uNH0k17Rx/RtNn8jJNJ3S/XF9nOYLs2bXM2AYb0/9U2jOsm6kSYYvHHC85wIbB8R2ZLvdma0c4+z+51p+OGCbj97Cdl69hX39OU2z4w00/zh8BLhPX5392u38xbh/n122vHhPdElSJ/aBSJI6MYFIkjoxgUiSOjGBSJI6MYFIkjoxgUiSOjGBSJI6MYFIkjoxgUiSOjGBSJI6+f81cKl4msdrIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#generate our random data\n",
    "data = np.random.normal(loc=0.0,scale=5.0,size = 1000)\n",
    "\n",
    "#plot a histogram of the data\n",
    "\n",
    "#set how many 'bins' or rectangles we want to divide the data up into for our histogram\n",
    "n_bins = 20\n",
    "\n",
    "#plot the histogram\n",
    "count,bins,ignored = plt.hist(data,n_bins,density=False,color=\"orange\",edgecolor='black')\n",
    "\n",
    "#set the x and y axis labels\n",
    "plt.xlabel(\"Observation Error\",size = \"xx-large\")\n",
    "plt.ylabel(\"Count\",size = \"xx-large\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding error to our linear model in python\n",
    "\n",
    "Now that we've seen how to generate normally-distributed random numbers, we can use them to add error into our linear model.\n",
    "\n",
    "To incorporate error, we just add a random number drawn from a normal distribution centered on 0 (this means the average value is 0 as set using the `loc` parameter in `np.random.normal`) to each y value. If you like, you can imagine taking the above histogram, rotating it 90 degrees clockwise, putting it on top of your models prediction line, then drawing one random number from it. Where that number falls will be the *observed* y value. Because the data are normally distributed, numbers close to the correct number are the most commonly observed values, but error will cause you to see y-values that are both bigger and smaller than the model's prediction value.\n",
    "\n",
    "The `scale` parameter, reflecting the standard deviation of the normal distribution, controls how 'spread out' these errors are — a larger standard deviation value results in a flatter, broader distribution corresponding to greater random error. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAESCAYAAADTx4MfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGQpJREFUeJzt3X+w3XV95/Hni4DaK90FBBwkhAQbV3E7QclmENwOVlfFukXsMosbFR2ctDs4aletKK26s5ttO+uPLTMLNq0/sF60DLKILbNKU8Vtt4KJgoKIUkJCJIXgb71dXON7//h+7+Ym3OTe+73n173n+Zg5c873c77n3PfJSfK6n+/n8/18U1VIkrRQRwy7AEnS0mSASJI6MUAkSZ0YIJKkTgwQSVInBogkqRMDRJLUiQEiSerEAJEkdXLksAvop+OPP75Wr1497DIkaUnZvn37I1V1wlz7LesAWb16Ndu2bRt2GZK0pCTZOZ/9PIQlSerEAJEkdWKASJI6MUAkSZ0YIJKkTgwQSRolOybhhtVwzRHN/Y7JYVd0SMt6Gq8kLSk7JuG2TbBvqtme2tlsA6zZOLy6DsEeiCSNijsu3x8e0/ZNNe0jyACRpFExtWth7UM2tABJckqSzyW5O8ldSd7Ytr87ybeT3N7eXjLjNW9Pcm+Se5K8aFi1S1JfTKxaWPuQDXMM5GfAm6vqy0l+Edie5Ob2ufdX1Xtm7pzkdOAi4JnAU4C/SvK0qto30KolqV/WbT5wDARgxUTTPoKG1gOpqj1V9eX28Y+Au4GTD/OS84FPVNWjVbUDuBfY0P9KJWlA1myEDVtg4lQgzf2GLSM5gA4jMgsryWrgWcCtwDnA65O8GthG00v5Hk24fHHGy3Zz+MCRpKVnzcaRDYyDDX0QPcnRwCeBN1XVD4GrgKcCZwB7gPdO7zrLy2uW99uUZFuSbXv37u1T1ZKkoQZIkqNowmOyqq4HqKqHqmpfVf0c+BP2H6baDZwy4+UrgQcPfs+q2lJV66tq/QknzLmcvSSpo2HOwgrwQeDuqnrfjPaTZux2AXBn+/hG4KIkj0+yBlgL3DaoeiVJBxrmGMg5wKuAryW5vW17B/CKJGfQHJ66H/hNgKq6K8m1wNdpZnBd6gwsSRqeoQVIVf0Ns49r3HSY12wGRnM+mySNmaEPokuSliYDRJLUiQEiSerEAJEkdWKASJI6MUAkSZ0YIJKkTgwQSVInBogkqRMDRJLUiQEiSerEAJEkdWKASJI6MUAkSZ0YIJKkTgwQSVInBogkqRMDRJLUiQEiScvJjkm4YTVcc0Rzv2Oybz9qaNdElyT12I5JuG0T7Jtqtqd2NtsAazb2/MfZA5Gk5eKOy/eHx7R9U017HxggkrRcTO1aWPsiGSCStFxMrFpY+yIZIJK0XKzbDCsmDmxbMdG094EBIkm9NMBZUI+xZiNs2AITpwJp7jds6csAOjgLS5J6Z8CzoGa1ZuPAfpY9EEnqlQHPgho2A0SSemXAs6CGzQCRpF4Z8CyoYTNAJKlXBjwLatgMEEnqlQHPgho2Z2FJUi8NcBbUsNkDkSR1YoBIkjoxQCRJnRggkqRODBBJUicGiCSpk6EFSJJTknwuyd1J7kryxrb9uCQ3J/lWe39s254kVyS5N8lXkzx7WLVLkobbA/kZ8OaqegZwFnBpktOBy4CtVbUW2NpuA5wHrG1vm4CrBl+yJGna0AKkqvZU1Zfbxz8C7gZOBs4Hrm53uxp4Wfv4fOCj1fgicEySkwZctiSpNRJjIElWA88CbgWeXFV7oAkZ4MR2t5OBB2a8bHfbJkkagqEHSJKjgU8Cb6qqHx5u11naapb325RkW5Jte/fu7VWZkqSDDDVAkhxFEx6TVXV92/zQ9KGp9v7htn03cMqMl68EHjz4PatqS1Wtr6r1J5xwQv+Kl6QxN8xZWAE+CNxdVe+b8dSNwMXt44uBT81of3U7G+ss4AfTh7okSYM3zNV4zwFeBXwtye1t2zuAPwCuTXIJsAu4sH3uJuAlwL3AFPDawZYrSZppaAFSVX/D7OMaAM+fZf8CLu1rUZKkeevZIawk65Nc0av3kySNtkUFSJKnJHlbkruA27CHIEljY8GHsJL8AvBymgHu59GE0Dbgd9k/4C1JWubmHSBJzgVeDfwGcDSwgyY8NlbVJ/pSnaSF2zEJd1wOU7tgYhWs2zw2l1jVYB02QJKspQmNVwKn0pwJfhUwCfwIuI9mRpSkUbBjEm7bBPvaf5ZTO5ttMETUc3P1QL4BfB+4Drimqm6ZfiLJqf0sTFIHd1y+Pzym7Ztq2g0Q9dhcg+jT02wL2NfnWiQt1tSuhbVLizBXgDwd+ADNUuq3JLk/yeZ22XVJo2Zi1cLapUU4bIBU1Ter6nJgNfBC4BbgDcDXgJtpeiZP7HONkuZr3WZYMXFg24qJpl3qsXmdB9Jeg2NrVV0MPJlmGZFdNAHysSRfSnJ5kl/uY62S5rJmI2zYAhOnAmnuN2xx/GMhdkzCDavhmiOa+x2Tw65oZKVZIaTji5OV7J+l9XSarFnRo9oWbf369bVt27ZhlyFpqTh4Fhs0PbgxC+Ek26tq/Vz7LepM9KraXVX/papOB56Dl5mVtJQdbhabHqNniylW1a00VxSUpKXJWWwLMteJhP9uoW9YVdd0L0eShmhiVXPy5Wzteoy5eiAfoxkoP9Sy69Nqxr0BImlpWrd59jEQZ7HNaq4AOW8e73Es8Bbg2cxyjXJJWjKmB8pdS2xeDhsgVfWZQz2X5Gjgt9vbMcCngd/raXWSxs+wF4Ncs9HAmKeuy7m/Hvgd4EnAZ4Hfq6ov9bg2SePGxSCXlHlP403yuCRvAP4e+EPgTuBXqurFhoeknujFNFpPBByYOXsgSVYArwMuB1YCXwReVVVb+1ybpHGz2Gm09mAG6rA9kCQXA98ErgT2Av+6qs42PCT1xWIXg/REwIGaqwfyYZqZVV8GbgDWJVl3mP2rqn6/V8VJGjOLnUbriYADNZ9B9ABntre5FGCASOpmsdNoPRFwoOYKkGcMpApJmraYabSeCDhQc50Hcs+gCpE0IoZ9HsZieCLgQHU5D+RI4BeBH1XVz3pfkjRkS/k/0MVaDrOYPBFwYOZ1HkiS05JclWQH8CjwCPBokh1Jrkzy1L5WKQ3K9H+gUzuB2v8f6LicS+AsJi3AnAGS5DzgDuA3gaNoLmV7fXv/OOC3gNuTvLiPdUqDMe7/gTqLSQsw13LuJwKTwHeAC6rqr2bZ54XAnwCTSZ5eVXv7Uqk0COP+H6izmLQAc/VAXgdMAC+aLTwAquqzwIuAo4FLelueNGCLPZFtqVu3uZm1NJOzmHQIcwXIC4BPzzUbq6q+AdwIvLBXhUlDMe7/ga7Z2Fz/e+JUIM39Qq8H7lpUY2M+54H813m+19/RrNArLV1OA13cLKblMItL8zZXgBxDM+NqPh4B/uniypFGgNNAuzvcJAT/TJeduQ5hPR6Y77ke+2hmZUkaV+M+CWHMzOdEwrVJzp7Hfk9bbDGSljhncY2V+QTIO9vbXILXRJfGm2tRjZW5AuTfD6QKScuDkxDGylyLKf7xoAqRtEw4CWFszPua6L2W5ENJHk5y54y2dyf5dpLb29tLZjz39iT3JrknyYuGU7UkadrQAgT4CDDb+lnvr6oz2ttNAElOBy4Cntm+5sr2Wu2SpCEZWoBU1ReA785z9/OBT1TVo1W1A7gX2NC34iRJcxpmD+RQXp/kq+0hrmPbtpOBB2bss7tte4wkm5JsS7Jt717XdZSkfhm1ALkKeCpwBrAHeG/bnln2nXXKcFVtqar1VbX+hBNO6E+VkqTRCpCqeqiq9lXVz2mWiJ8+TLUbOGXGriuBBwddnyRpv5EKkCQnzdi8AJieoXUjcFGSxydZA6wFbht0fZKk/Q55HkiSmzq8X1XVr81nxyQfB84Fjk+yG3gXcG6SM2gOT91PcxVEququJNcCX6dZm+vSqtrXoT5JUo8c7kTC4+jj0iRV9YpZmj94mP03A66HIEkj4pABUlVnDbIQSdLSMlJjIJKkpcMAkSR1sqAASfKcJNcneSDJj5NMHXT7Sb8KlcaG1xTXEjHvAEnyPOAW4JeBzwITwF8AnwFWAHcDV/ShRml8TF9TfGonUPuvKW6IaAQtpAfyTuAemgB5W9v2gaq6ADiT5gzyL/a2PGnMHO6a4tKIWUiAnAl8tKr+D/Dztm0FQFXdCWwBfre35UljxmuKawlZSIAE+F77ePpXpGNmPP9N4PReFCWNrUNdO9xrimsELSRAHgBOA2h7IQ8Az5nx/Dr2B4ykLtZtbq4hPpPXFNeImuua6DN9Hvh14B3t9iTw1iSPpzmUdQnNYSxJXXlNcS0hCwmQ9wB/m+QJbQ/k3cCTgdcA+4CPA5f1ukBp7HhNcS0R8w6QqroPuG/G9k9peh2X9KEuSdKIW8h5IFcmWX+Y589McmVvypIkjbqFDKL/FvC0wzz/S7TLr0uSlr9eroV1HPBoD99PkjTCDjsGkuRs4Lkzml6aZOUsux4LvJLmgk+SpDEw1yD6v6K5UiA0F5e6qL3N5iHgrT2qS5I04uYKkCuB62jOQv8qzTTdvzxonwJ+DDxQVX27gqEkabQcNkCqai+wFyDJecAdVfUPgyhMkjTaFnIeyGemHyc5FlhN0/vYWVUuYSJJY2ahF5Q6M8kXgEeAbcB24JEkX0hyZj8KlCSNpnn3QJI8G/gCzXjIn9HMuArwDOBC4JYk/7KqvtKPQiVJo2Uha2H9Z5rVds+pqp0zn0jyTuB/A/8JeGnvypMkjaqFHMI6G7jq4PAAqKpdwAc48JwRSdIytpAAOQr4x8M8P9XuI0kaAwsJkK8Cr0kycfATbdtrgDt6VJckacQtZAzk94EbgNuTXAV8g2Ya7+k0Cy0+FXhZzyuUJI2khZwHcmOSi4H3trfps85DM633NVX16d6XKEkaRQvpgVBVf5bkz2muhb66bb4f+Lv2AlOSpDEx12q8HwL+uKpunW5rg+KW9iZJGlNzDaK/hmZsQ5KkA/TyglKSpDFigEiSOpnPIPqJSU6b7xtW1X2LqEeStETMJ0Cmp+3O14qOtUiSlpD5BMgNNGehS5qPHZNwx+UwtQsmVsG6zbBm47CrknpuPgHyyaq6pu+VSMvBjkm4bRPsm2q2p3Y222CIaNkZ2iB6kg8leTjJnTPajktyc5JvtffHtu1JckWSe5N8tb02iTR67rh8f3hM2zfVtEvLzDBnYX0EePFBbZcBW6tqLbC13QY4D1jb3jYBVw2oRmlhpnYtrF1awoYWIFX1BeC7BzWfD1zdPr6a/Yszng98tBpfBI5JctJgKpUWYGLVwtqlJeywAVJVRwx4/OPJVbWn/dl7gBPb9pOBB2bst7ttk0bLus2w4qArHqyYaNqlZWapnEiYWdpqljaSbEqyLcm2vXv39rks6SBrNsKGLTBxKpDmfsMWB9C1LC1oNd4BeCjJSVW1pz1E9XDbvhs4ZcZ+K4EHZ3uDqtoCbAFYv379rCEj9dWajQaGxsKo9UBuBC5uH18MfGpG+6vb2VhnAT+YPtQlSRqOofVAknwcOBc4Pslu4F3AHwDXJrkE2AVc2O5+E/AS4F6aa6+/duAFS5IOMLQAqapXHOKp58+ybwGX9rciSdJCjNohLEnSEmGASJI6MUC0/OyYhBtWwzVHNPc7JoddkbQsjdo0XmlxXMxQGhh7IDrQUv/t3cUMpYGxB6L9lsNv7y5mKA2MPRDtNyq/vS+mF+RihtLAGCDabxR+e5/uBU3tBGp/L2i+IeJihtLAGCDabxR+e19sL8jFDKWBcQxE+63bfOAYCHT77X0x1wTvRS/IxQylgbAHov168dv7Yg9BjUIvSNK82APRgRb72/vhDkHN53171QuS1Hf2QNRbiz0E5RiGtGTYA1luFjP+0AsTq9rDV7O0z5djGNKSYA9kOVns+EMvOI1WGhsGyHIyCicCeghKGhsewlpORuFEQPAQlDQm7IEsJ06BlTRABshy4viDpAEyQJYTxx8kDZBjIMuN4w+SBsQeiCSpEwNEktSJASJJ6sQAkSR1YoBIkjoxQCRJnRggkqRODBBJUicGiCSpEwNEktSJASJJ6sQAkSR1YoBIkjoxQCRJnRggkqRODBBJUicGiCSpk5EMkCT3J/laktuTbGvbjktyc5JvtffHDrvOvtgxCTeshmuOaO53TA67Ikma1UgGSOt5VXVGVa1vty8DtlbVWmBru7287JiE2zbB1E6gmvvbNhkikkbSKAfIwc4Hrm4fXw28bIi19Mcdl8O+qQPb9k017ZI0YkY1QAr4bJLtSTa1bU+uqj0A7f2JQ6uuX6Z2LaxdkoboyGEXcAjnVNWDSU4Ebk7yjfm+sA2cTQCrVq3qV339MbGqPXw1S7skjZiR7IFU1YPt/cPA/wA2AA8lOQmgvX/4EK/dUlXrq2r9CSecsPAfPsxB7HWbYcXEgW0rJpp2SRoxIxcgSZ6Y5BenHwMvBO4EbgQubne7GPhUz3/4sAex12yEDVtg4lQgzf2GLU27JI2YVNWwazhAktNoeh3QHGK7pqo2J3kScC2wCtgFXFhV3z3ce61fv762bds2/x9+w+pDHEI6FV52//zeY8dkM+g9tas59LRuswEgaUlJsn3GDNhDGrkxkKq6D1g3S/t3gOf39YcvdhB7ugczPZNqugcDhoikZWfkDmEN1aEGq+c7iO00XEljxACZabGD2E7DlTRGDJCZFjuIvdgejCQtISM3BjJ0azZ2H69Yt/nAMRBwGq6kZcseSC85DVfSGLEH0muL6cFI0hJiD0SS1IkBIknqxACRJHVigEiSOjFAJEmdjNxiir2UZC8wy+qI83I88EgPy1lqxv3zg38Gfv7x/fynVtWc18NY1gGyGEm2zWc1yuVq3D8/+Gfg5x/vzz8fHsKSJHVigEiSOjFADm3LsAsYsnH//OCfgZ9fh+UYiCSpE3sgkqRODJBZJHlxknuS3JvksmHX029JTknyuSR3J7kryRvb9uOS3JzkW+39scOutZ+SrEjylSR/0W6vSXJr+/n/PMnjhl1jvyQ5Jsl1Sb7R/j14zjh9/0l+u/27f2eSjyd5wjh9/10ZIAdJsgL478B5wOnAK5KcPtyq+u5nwJur6hnAWcCl7We+DNhaVWuBre32cvZG4O4Z238IvL/9/N8DLhlKVYPxR8D/rKqnA+to/hzG4vtPcjLwBmB9Vf1zYAVwEeP1/XdigDzWBuDeqrqvqn4KfAI4f8g19VVV7amqL7ePf0Tzn8fJNJ/76na3q4GXDafC/kuyEvg14E/b7QC/ClzX7rJsP3+SfwL8CvBBgKr6aVV9nzH6/mkubfELSY4EJoA9jMn3vxgGyGOdDDwwY3t32zYWkqwGngXcCjy5qvZAEzLAicOrrO/+G/A7wM/b7ScB36+qn7Xby/nvwWnAXuDD7SG8P03yRMbk+6+qbwPvAXbRBMcPgO2Mz/ffmQHyWJmlbSymqiU5Gvgk8Kaq+uGw6xmUJC8FHq6q7TObZ9l1uf49OBJ4NnBVVT0L+AnL9HDVbNqxnfOBNcBTgCfSHMI+2HL9/jszQB5rN3DKjO2VwINDqmVgkhxFEx6TVXV92/xQkpPa508CHh5WfX12DvDrSe6nOWT5qzQ9kmPaQxqwvP8e7AZ2V9Wt7fZ1NIEyLt//C4AdVbW3qv4vcD1wNuPz/XdmgDzWl4C17QyMx9EMpt045Jr6qj3e/0Hg7qp634ynbgQubh9fDHxq0LUNQlW9vapWVtVqmu/7r6tqI/A54N+0uy3nz/8PwANJ/lnb9Hzg64zJ909z6OqsJBPtv4Xpzz8W3/9ieCLhLJK8hOY30BXAh6pq85BL6qskzwX+F/A19o8BvINmHORaYBXNP7ILq+q7QylyQJKcC7ylql6a5DSaHslxwFeAV1bVo8Osr1+SnEEzgeBxwH3Aa2l+wRyL7z/JfwT+Lc2MxK8Ar6MZ8xiL778rA0SS1ImHsCRJnRggkqRODBBJUicGiCSpEwNEktSJASL1WJLPJ/n8sOuQ+s0AkeaQpOZ5+8iwa5UG6ci5d5HG3qsO2n45cAHwFuChGe1/396/cBBFScPmiYTSAiV5N/AuYG1V3TvkcqSh8RCW1GOzjYG0h7g+luQFSb6U5B+TfDPJhe3zZyT56yQ/SfJgkv9wiPe+IMnfJvlxe9ua5OwBfCzpMQwQaXDWAZPAXwJva9s+0YbIZ4BtwFtpVsd9b5IDDoUleRPNSrHfBd5O0ws6CfhcknMG8gmkGTyEJS3QXIewpnsfVXXujLaiWajyX0xf/THJM4E7aa4zcX5VfbptfxJNiNxUVb/Rtq2kWeTwiqp6y4z3PZpm5dhdVfXcXn9W6XAcRJcG59bp8ACoqruS/AD4yXR4tO3fSXIP8EszXvty4CjgmiTHH/S+W4FXJ5moqqk+1i8dwACRBmfnLG3fp7mM6mztT52x/fT2fvss+057EmCAaGAMEGlw9i2wfeZldafHKy8AfnyI/fd2KUrqygCRlobpsZYHq+q2oVYitZyFJS0N19FcLe/dM67T/f8lOXHwJWnc2QORloCquj/Jm4E/ArYnuRZ4GFgJnEszw+t5w6tQ48gAkZaIqrqinZ31ZprzRZ5AMwB/K/DhYdam8eR5IJKkThwDkSR1YoBIkjoxQCRJnRggkqRODBBJUicGiCSpEwNEktSJASJJ6sQAkSR1YoBIkjr5f3A8IVGPdV06AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def linear_simulation(x_data,slope,error_stdev = 20, intercept = 0):\n",
    "    \"\"\"Return y data given x data and parameters\n",
    "    x_data -- a list of x data points (ints or floats)\n",
    "    slope -- the slope of our linear model (m in y=mx+b+error)\n",
    "    error_stdev -- the standard deviation of normal error\n",
    "    intercept -- the intercept of our linear model (b in y=mx+b+error)\n",
    "    \"\"\"\n",
    "    y_predictions = []\n",
    "    for i,x in enumerate(x_data):\n",
    "        error = np.random.normal(loc=0.0,scale = error_stdev,size = None)\n",
    "        y = slope*x + error + intercept \n",
    "        y_predictions.append(y)\n",
    "        \n",
    "    return y_predictions\n",
    "\n",
    "\n",
    "x_data = list(range(0,100,5))    \n",
    "y_data = linear_simulation(x_data = x_data,slope = 2.0,intercept =35.0,error_stdev=20.0)\n",
    "plot_scatterplot(x_data,y_data,xlabel=\"Time\",ylabel=\"Total DNA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Linear regression\n",
    "\n",
    "In the last section we imagined what would happen if we had a perfect biological model of a process. We saw that due to measurement error, observations will always deviate from model predictions, and noted that we could measure the extent of this difference between prediction and observation using residuals.\n",
    "\n",
    "Now let's imagine that we are in a situation that is somewhat the opposite of our previous one. Imagine that we do not know for sure the rate of DNA replication in a cell, but we have collected some data on the amount of DNA in various cells at different points in time during the replication process. If we were to plot this data, we might see something like this:\n",
    "\n",
    "<img src=\"./resources/linear_scatter_high_r2_linear_regression_no_line.png\" width=\"400\"  description=\"A cartoon of a scatter plot with x and y axes, with orange points and a black regression line just as in the previous figure. Now, however, thin grey lines stretch vertically from each point to the regression line. These are the residuals — the difference between what the regression model predicts and the actual observed data.\">\n",
    "\n",
    "Now it looks like a line might be a reasonable approximation for these points. But if the real trend is linear, how do we best estimate the slope and intercept of that line? This is the question addressed by **linear regression**. \n",
    "\n",
    "There are multiple approaches to how to pick the line that best explains the data in linear regression. We noted above that residuals represent differences between the model and the real data, so it probably makes sense that we somehow want to pick a line that keeps those residuals small. \n",
    "\n",
    "The simplest approach to picking which line is best, **ordinary least squares**, picks a line that minimizes the squared residuals (hence 'least squares'). This means small residuals have little effect, but large residuals are heavily penalized. But what is the justification for weighting by the square of the residual, rather than, for instance, just the residuals themselves?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normally distributed errors and the Central Limit Theorem\n",
    "\n",
    "**Ordinary least squares** regression assumes that measurement errors are **normally distributed** - that is, that they follow a particular kind of bell curve. If so, then the chances of getting small deviations from our expected value is large, whereas the chances of getting very large deviations from our expected value falls very quickly. For many people, this may be at least a little bit intuitive: you may notice by eye that small deviations from the regression line don't \"look bad\", whereas larger deviations are very noticable. \n",
    "\n",
    "But why on Earth should we think that errors are normally distributed? If you aren't already in the know, it may seem quite odd that everyone is always assuming normally distributed errors. We'll get back to linear regression for a moment, but this point is important enough that it deserves a detour. \n",
    "\n",
    "The reason linear regression typically assumes that errors are normally distributed comes from a very famous result in statistics known as the **Central Limit Theorem**. The [Central Limit Theorem](https://en.wikipedia.org/wiki/Central_limit_theorem) says that as the number of observations increases, the sum of independent and identically distributed <abbr title = \"random numbers drawn from a distribution\"> random variates </abbr> form a normal distribution - even if the distribution those measurements come from isn't itself normal.\n",
    "\n",
    "Here's a simple example: the chances of rolling each number on a fair die are equal — rolling 1,2,3,4,5,6 each have a 1/6th chance. Statisticians would call this one random variate from a uniform distribution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Possible results of rolling one die:<br>\n",
    "&nbsp;1:&#9861;<br>\n",
    "&nbsp;2:&#9861;<br>\n",
    "&nbsp;3:&#9861;<br>\n",
    "&nbsp;4:&#9861;<br>\n",
    "&nbsp;5:&#9861;<br>\n",
    "&nbsp;6:&#9861;<br>\n",
    "\n",
    "Notice how there is just one way to roll each number, so the results above look 'flat'. This is a *uniform distribution*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " But what happens when we roll more dice? The distribution formed by the *sum* of two dice rolls — that is two <abbr title=\"A uniform probability distribution is a flat distribution, meaning if you plotted all possible outcomes on the x-axis  against their probability on the y-axis, you would just get a straight line.\"> uniform distributions</abbr> — is not itself uniform. Instead, there are more ways to get a 7 — which can be formed by 1 + 6, 2+5, 3+4, etc — than there are a 12 — which can only be formed by rolling two 6's. If we make a histogram showing how often each total number comes up, rolling 1 die will look flat, but 2 dice will look triangular (peaked at 7), while we will see that rolling more dice will approximate a bell curve-shaped normal distribution.\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us illustrate this with a table of all 36 possible rolls of two dice and their sums:\n",
    " \n",
    "**Possible sums of rolling two dice:**<br>\n",
    "&nbsp;2:&#9856;+&#9856;<br>\n",
    "&nbsp;3:&#9856;+&#9857; or &#9857;+&#9856;<br>\n",
    "&nbsp;4:&#9856;+&#9858; or &#9857;+&#9857; or &#9858;+&#9856;<br>\n",
    "&nbsp;5:&#9856;+&#9859; or &#9857;+&#9858; or &#9858;+&#9857; or &#9859;+&#9856;<br>\n",
    "&nbsp;6:&#9856;+&#9860; or &#9857;+&#9859; or &#9858;+&#9858; or &#9859;+&#9857; or &#9860;+&#9856;<br>\n",
    "&nbsp;7:&#9856;+&#9861; or &#9857;+&#9860; or &#9858;+&#9859; or &#9859;+&#9858; or &#9860;+&#9857; or &#9861;+&#9856;<br>\n",
    "&nbsp;8:&#9857;+&#9861; or &#9858;+&#9860; or &#9859;+&#9859; or &#9860;+&#9858; or &#9861;+&#9857;<br>\n",
    "&nbsp;9:&#9858;+&#9861; or &#9859;+&#9860; or &#9860;+&#9859; or &#9861;+&#9858;<br>\n",
    "10:&#9859;+&#9861; or &#9860;+&#9860; or &#9861;+&#9859;<br>\n",
    "11:&#9860;+&#9861; or &#9861;+&#9860;<br>\n",
    "12:&#9861;+&#9861;<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    " \n",
    "If we do the same thing with 3 or more dice, the result will look progressively more and more like a Bell curve shaped normal distribution:\n",
    "\n",
    "**Possible sums of rolling three dice:**\n",
    "&nbsp;3:&#9856;&#9856;&#9856;<br>\n",
    "&nbsp;4:&#9856;&#9856;&#9857;,&#9856;&#9857;&#9856;,&#9857;&#9856;&#9856;<br>\n",
    "&nbsp;5:&#9856;&#9856;&#9858;,&#9856;&#9857;&#9857;,&#9856;&#9858;&#9856;,&#9857;&#9856;&#9857;,&#9857;&#9857;&#9856;,&#9858;&#9856;&#9856;<br>\n",
    "&nbsp;6:&#9856;&#9856;&#9859;,&#9856;&#9857;&#9858;,&#9856;&#9858;&#9857;,&#9856;&#9859;&#9856;,&#9857;&#9856;&#9858;,&#9857;&#9857;&#9857;,&#9857;&#9858;&#9856;,&#9858;&#9856;&#9857;,&#9858;&#9857;&#9856;,&#9859;&#9856;&#9856;<br>\n",
    "&nbsp;7:&#9856;&#9856;&#9860;,&#9856;&#9857;&#9859;,&#9856;&#9858;&#9858;,&#9856;&#9859;&#9857;,&#9856;&#9860;&#9856;,&#9857;&#9856;&#9859;,&#9857;&#9857;&#9858;,&#9857;&#9858;&#9857;,&#9857;&#9859;&#9856;,&#9858;&#9856;&#9858;,&#9858;&#9857;&#9857;,&#9858;&#9858;&#9856;,&#9859;&#9856;&#9857;,&#9859;&#9857;&#9856;,&#9860;&#9856;&#9856;<br>\n",
    "&nbsp;8:&#9856;&#9856;&#9861;,&#9856;&#9857;&#9860;,&#9856;&#9858;&#9859;,&#9856;&#9859;&#9858;,&#9856;&#9860;&#9857;,&#9856;&#9861;&#9856;,&#9857;&#9856;&#9860;,&#9857;&#9857;&#9859;,&#9857;&#9858;&#9858;,&#9857;&#9859;&#9857;,&#9857;&#9860;&#9856;,&#9858;&#9856;&#9859;,&#9858;&#9857;&#9858;,&#9858;&#9858;&#9857;,&#9858;&#9859;&#9856;,&#9859;&#9856;&#9858;,&#9859;&#9857;&#9857;,&#9859;&#9858;&#9856;,&#9860;&#9856;&#9857;,&#9860;&#9857;&#9856;,&#9861;&#9856;&#9856;<br>\n",
    "&nbsp;9:&#9856;&#9857;&#9861;,&#9856;&#9858;&#9860;,&#9856;&#9859;&#9859;,&#9856;&#9860;&#9858;,&#9856;&#9861;&#9857;,&#9857;&#9856;&#9861;,&#9857;&#9857;&#9860;,&#9857;&#9858;&#9859;,&#9857;&#9859;&#9858;,&#9857;&#9860;&#9857;,&#9857;&#9861;&#9856;,&#9858;&#9856;&#9860;,&#9858;&#9857;&#9859;,&#9858;&#9858;&#9858;,&#9858;&#9859;&#9857;,&#9858;&#9860;&#9856;,&#9859;&#9856;&#9859;,&#9859;&#9857;&#9858;,&#9859;&#9858;&#9857;,&#9859;&#9859;&#9856;,&#9860;&#9856;&#9858;,&#9860;&#9857;&#9857;,&#9860;&#9858;&#9856;,&#9861;&#9856;&#9857;,&#9861;&#9857;&#9856;<br>\n",
    "10:&#9856;&#9858;&#9861;,&#9856;&#9859;&#9860;,&#9856;&#9860;&#9859;,&#9856;&#9861;&#9858;,&#9857;&#9857;&#9861;,&#9857;&#9858;&#9860;,&#9857;&#9859;&#9859;,&#9857;&#9860;&#9858;,&#9857;&#9861;&#9857;,&#9858;&#9856;&#9861;,&#9858;&#9857;&#9860;,&#9858;&#9858;&#9859;,&#9858;&#9859;&#9858;,&#9858;&#9860;&#9857;,&#9858;&#9861;&#9856;,&#9859;&#9856;&#9860;,&#9859;&#9857;&#9859;,&#9859;&#9858;&#9858;,&#9859;&#9859;&#9857;,&#9859;&#9860;&#9856;,&#9860;&#9856;&#9859;,&#9860;&#9857;&#9858;,&#9860;&#9858;&#9857;,&#9860;&#9859;&#9856;,&#9861;&#9856;&#9858;,&#9861;&#9857;&#9857;,&#9861;&#9858;&#9856;<br>\n",
    "11:&#9856;&#9859;&#9861;,&#9856;&#9860;&#9860;,&#9856;&#9861;&#9859;,&#9857;&#9858;&#9861;,&#9857;&#9859;&#9860;,&#9857;&#9860;&#9859;,&#9857;&#9861;&#9858;,&#9858;&#9857;&#9861;,&#9858;&#9858;&#9860;,&#9858;&#9859;&#9859;,&#9858;&#9860;&#9858;,&#9858;&#9861;&#9857;,&#9859;&#9856;&#9861;,&#9859;&#9857;&#9860;,&#9859;&#9858;&#9859;,&#9859;&#9859;&#9858;,&#9859;&#9860;&#9857;,&#9859;&#9861;&#9856;,&#9860;&#9856;&#9860;,&#9860;&#9857;&#9859;,&#9860;&#9858;&#9858;,&#9860;&#9859;&#9857;,&#9860;&#9860;&#9856;,&#9861;&#9856;&#9859;,&#9861;&#9857;&#9858;,&#9861;&#9858;&#9857;,&#9861;&#9859;&#9856;<br>\n",
    "12:&#9856;&#9860;&#9861;,&#9856;&#9861;&#9860;,&#9857;&#9859;&#9861;,&#9857;&#9860;&#9860;,&#9857;&#9861;&#9859;,&#9858;&#9858;&#9861;,&#9858;&#9859;&#9860;,&#9858;&#9860;&#9859;,&#9858;&#9861;&#9858;,&#9859;&#9857;&#9861;,&#9859;&#9858;&#9860;,&#9859;&#9859;&#9859;,&#9859;&#9860;&#9858;,&#9859;&#9861;&#9857;,&#9860;&#9856;&#9861;,&#9860;&#9857;&#9860;,&#9860;&#9858;&#9859;,&#9860;&#9859;&#9858;,&#9860;&#9860;&#9857;,&#9860;&#9861;&#9856;,&#9861;&#9856;&#9860;,&#9861;&#9857;&#9859;,&#9861;&#9858;&#9858;,&#9861;&#9859;&#9857;,&#9861;&#9860;&#9856;<br>\n",
    "13:&#9856;&#9861;&#9861;,&#9857;&#9860;&#9861;,&#9857;&#9861;&#9860;,&#9858;&#9859;&#9861;,&#9858;&#9860;&#9860;,&#9858;&#9861;&#9859;,&#9859;&#9858;&#9861;,&#9859;&#9859;&#9860;,&#9859;&#9860;&#9859;,&#9859;&#9861;&#9858;,&#9860;&#9857;&#9861;,&#9860;&#9858;&#9860;,&#9860;&#9859;&#9859;,&#9860;&#9860;&#9858;,&#9860;&#9861;&#9857;,&#9861;&#9856;&#9861;,&#9861;&#9857;&#9860;,&#9861;&#9858;&#9859;,&#9861;&#9859;&#9858;,&#9861;&#9860;&#9857;,&#9861;&#9861;&#9856;<br>\n",
    "14:&#9857;&#9861;&#9861;,&#9858;&#9860;&#9861;,&#9858;&#9861;&#9860;,&#9859;&#9859;&#9861;,&#9859;&#9860;&#9860;,&#9859;&#9861;&#9859;,&#9860;&#9858;&#9861;,&#9860;&#9859;&#9860;,&#9860;&#9860;&#9859;,&#9860;&#9861;&#9858;,&#9861;&#9857;&#9861;,&#9861;&#9858;&#9860;,&#9861;&#9859;&#9859;,&#9861;&#9860;&#9858;,&#9861;&#9861;&#9857;<br>\n",
    "15:&#9858;&#9861;&#9861;,&#9859;&#9860;&#9861;,&#9859;&#9861;&#9860;,&#9860;&#9859;&#9861;,&#9860;&#9860;&#9860;,&#9860;&#9861;&#9859;,&#9861;&#9858;&#9861;,&#9861;&#9859;&#9860;,&#9861;&#9860;&#9859;,&#9861;&#9861;&#9858;<br>\n",
    "16:&#9859;&#9861;&#9861;,&#9860;&#9860;&#9861;,&#9860;&#9861;&#9860;,&#9861;&#9859;&#9861;,&#9861;&#9860;&#9860;,&#9861;&#9861;&#9859;<br>\n",
    "17:&#9860;&#9861;&#9861;,&#9861;&#9860;&#9861;,&#9861;&#9861;&#9860;<br>\n",
    "18:&#9861;&#9861;&#9861;<br>\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Central Limit Theorem tells us this isn't just a weird special property of rolling dice and adding up their values, but instead is a *general* property of adding up independent random values drawn from any kind of distribution. \n",
    "\n",
    "This means that if we have many different sources of error that influence our result, as long there are lots of them and our total measurement error can be approximated by the sum of those individual errors, then we can use a normal distribution to model the types of measurement errors we expect and their properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ordinary Least Squares helps find the model parameters that best explain the data if errors are normal\n",
    "\n",
    "In the previous section we saw that the Central Limit Theorem means that the sum of a large number of random variates from many statistical distributions will itself be normally distributed, even if the underlying statistical distribution(s) are not. \n",
    "\n",
    "This means that a normal distribution is often a pretty good approximation for **measurement error**. That's very good news when true, because if we have a particular normal distribution, we can calculate the exact chances of drawing a particular number from it. \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Linear regression in Python\n",
    "\n",
    "We can infer a regression model that best explains the data using **ordinary least squares** regression, which minimizes the sum of squared difference between each data point and the predicted line (x - mean)^2. This also maximizes the likelihood of the model, since when drawing random numbers from a normal distribution, the probability of getting a given level of error (distance from the mean) deceases with the square of the size of that error.\n",
    "\n",
    "The `statsmodels` python package provides a nice interface to do linear regression that can handle more complex cases as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.905\n",
      "Model:                            OLS   Adj. R-squared:                  0.899\n",
      "Method:                 Least Squares   F-statistic:                     170.6\n",
      "Date:                Thu, 15 Jul 2021   Prob (F-statistic):           1.28e-10\n",
      "Time:                        14:24:55   Log-Likelihood:                -87.304\n",
      "No. Observations:                  20   AIC:                             178.6\n",
      "Df Residuals:                      18   BIC:                             180.6\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept     38.4136      8.646      4.443      0.000      20.248      56.579\n",
      "x              2.0323      0.156     13.061      0.000       1.705       2.359\n",
      "==============================================================================\n",
      "Omnibus:                        4.859   Durbin-Watson:                   2.344\n",
      "Prob(Omnibus):                  0.088   Jarque-Bera (JB):                3.285\n",
      "Skew:                           0.988   Prob(JB):                        0.193\n",
      "Kurtosis:                       3.186   Cond. No.                         107.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "import statsmodels as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame({\"x\":x_data,\"y\":y_data})\n",
    "df.set_index(\"x\")\n",
    "\n",
    "model = smf.ols(formula='y ~ x', data=df)\n",
    "results = model.fit()\n",
    "summary = results.summary()\n",
    "print(summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting our results. \n",
    "\n",
    "**Estimated slope and intercept**. The number under coef. (coefficient) tells you the intercept (for the row labelled Intercept) and the slope on our one x variable (in the row labelled x). The final two rows also give you a 95% confidence interval on that number.\n",
    "\n",
    "Note that the slope and intercept from linear regression will almost always be at least a little bit wrong compared to the true values. When I ran our above code, I got an inferred slope of ~1.7 instead of 2, and an inferred intercept of ~52 rather than 35. However, both the true slope and true intercept are within the 95% confidence interval of the estimate.\n",
    "\n",
    "However, assuming the assumptions of linear regression are met (e.g. error is normal, the trend really is linear, data points are independent and identically distributed observations), the more data available, the closer linear regression will get to the right answer.\n",
    "\n",
    "**Effect Size**. The 'effect size' of a linear regression is reported in the R<sup>2</sup> value. It says what proportion of variance in the data is explained by your model. A value close to 1.0 means that your model's predictions predict virtually all the differences in `y` values that you observe. A value close to 0.0 means that your model predicts almost nothing about `y`. I got an R<sup>2</sup> value of ~0.9 indicating that about 90% of the variance in y in this dataset could be explained by the model. \n",
    "\n",
    "**Significance**. There are multiple ways to test how good of a job a model does at explaining data. We will explore these later. For now, note that there is a P-value column in the results. The p-value  represents the chances of getting as good of a fit as you actually got or better under the null hypothesis that the slope is 0 (e.g. a flat line where x has no effect on y). The intercept represents the same for the null hypothesis that the intecept is 0. p-values that are below 0.05 are traditionally considered significant. In my results both the intercept and slope p-values are 0.0, indicating the data are extremely unlikely under either null hypothesis.\n",
    "\n",
    "We will return in much greater detail to how model parameters (like slope and intercept), effect size (like R<sup>2</sup>) and significance iteract in another chapter.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple Linear Regression\n",
    "\n",
    "Multiple linear regression refers to linear regression with more than one predictor variable. These predictors might be continuous or categorical. Let's start with an example that includes categorical data, since these are easiest to draw. We'll apply the same ideas to continuous data afterwards. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To develop a graphical intuition for why multiple linear regression is often needed, let's imagine that when we did our DNA replication rate experiments in a second cell type, and colored our results from each type of cell (say skin vs. bone cells) a different color.  Our data look like this:\n",
    "\n",
    "<img src=\"./resources/linear_scatter_multivariate.png\" width=\"400\"  description=\"Another scatter plot. Now however there are two sets of points - one set orange, and one cyan. These sets of points roughly form separate diagonal lines, with the orange line offset so it is above the cyan one\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could, of course, ignore the type of cell we collected the data from and try to run a linear regression. However, the results would not be very satisfactory:\n",
    "\n",
    "<img src=\"./resources/linear_scatter_multivariate_with_single_regression.png\" width=\"400\"  description=\"Another scatter plot with two sets of points — one orange and one cyan — trending upwards and to the right. A black regression line is plotted. It does head up and to the right, but it falls between the orange and cyan points. It is clear that it will always underestimate the value of the orange points, and overestimate the cyan ones.\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While our line does fall roughly in the middle of all the data, it is clear that it will always underestimate the value of the orange points, and overestimate the cyan ones. To my eye, at least, it feels like we need *two* lines, one offset from the other by a fixed amount in order to explain these data:\n",
    "\n",
    "<img src=\"./resources/linear_scatter_multivariate_with_two_regression.png\" width=\"400\"  description=\"Another scatter plot with two sets of points — one orange and one cyan — trending upwards and to the right. Now each color of points has it's own regression line, and they seem to match the data pretty well.\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How could we incorporate this into our equation for a line? One way would be to let our line represent the value of the blue points, but add a positive number to the prediction of the line whenever we got to an orange point. This would have the effect of predicting the blue points with our line, and our orange points with our line plus some fixed value. But how can we represent that mathematically?\n",
    "\n",
    "The trick is to recognize that our categorical value (orange vs. blue) is binary in nature — each point is either orange or blue — and we can therefore represent one value with a `1` and the other with a `0`. This approach is sometimes called generating a 'dummy variable' to hold our categorical data.\n",
    "\n",
    "Let us rewrite our equation for a line using full words instead of variables. We will now include a *second* x value that represents whether a point is orange (1) or blue(0). That second x value, just like the first x value will be associated with it's own slope, which we'll for now call `effect_of_color`. \n",
    "\n",
    "`y = mx + is_orange * effect_of_being_orange + b + error`\n",
    "\n",
    "The is_orange parameter above is just 1 if the point is orange, or 0 otherwise. Since 0 times any number is just 0, this means that the `effect_of_color` value will only be added to the y value if we are dealing with an orange point! \n",
    "\n",
    "Now of course, conventionally, statistical equations aren't written out in words but rather with more concise symbols. Here's a critical observation: the parameter `is_orange` is effectively just another x value, and the parameter `effect of being orange` is effectively just another slope that gets multiplied by some x values to contribute to the result. \n",
    "\n",
    "That is, we now have a **multiple linear regression** with *two* x values that each can have different slopes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how the equation might look if we modified the most familiar form of the equation for a line:\n",
    "\n",
    "y = m<sub>1</sub>x<sub>1</sub> + m<sub>2</sub>x<sub>2</sub> + b + error\n",
    "\n",
    "Now in our DNA replication example we would have two x values: the amount of time that has passed (x<sub>1</sub>), and the type of cell (x<sub>2</sub>). The second x value would be just like our `is_orange` variable up above.\n",
    "\n",
    "If we extended our 'mathy' looking version of the same equation, it might look something like this:\n",
    "\n",
    "y = &beta;<sub>0</sub>x<sub>0</sub>  + &beta;<sub>1</sub>x<sub>1</sub> + &beta;<sub>2</sub> + &epsilon;\n",
    "    \n",
    "Remember that despite this now looking more mathematically formal, we're representing a really really simple idea: we have a linear model that predicts y from our first x value (x<sub>0</sub>). We also have two colors, one is offset from the other, so we just bump our estimate of y up a bit (specifically an amount described by &beta;<sub>1</sub>) if the point is orange (x<sub>1</sub> is not equal to 0). \n",
    "\n",
    "What if we had more than two categories (e.g. 3 colors, 4 cell types, etc)? If these categories aren't aranged in a particular order, we can't use a single binarry variable to represent this. Instead we would need to generate several 'dummy variables' that can represent different colors or cell types as a 0 or 1. So for instance, if we had measured 4 blood cell types (let's say basophils, thrombocytes, macrophage, and plasma cells) we might have variables is_basophil, is_thrombocyte, and is_macrophage, each of which is a 1 if a cell falls into that category. Note that we need 1 fewer dummy variable than the number of cell types, because in our analysis anything that isn't a basophil, thrombocyte or macrophage must be a plasma cell, since those are the only cell types we measured. Using this general trick of converting categorical data into many binary variables, we can represent any categorical data we choose in our linear models in a simple way.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating a multiple linear model in python\n",
    "\n",
    "We can simulate a multiple linear model in python by modifying our previous simulation.\n",
    "\n",
    "If we have two x values, one reflecting some quantitative measurement, and another a category. It will take a little more work to plot, but the idea is the same: we take in an extra x variable (let's call it `x2`), define a second slope (let's uncreatively call it `slope2`) and then just add an extra x2 * slope2 term to our equation that predicts `y`. Otherwise this is just the familiar equation for a line, with normal errors added:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_simulation(x_data,x2_data,slope,slope2,intercept = 0,error_stdev = 0):\n",
    "    \"\"\"Return y data given x data, x2 data and parameters\n",
    "    x_data -- a list of x data points (ints or floats)\n",
    "    x2_data -- a list of x2 data points (ints or floats), that must be of the same length as x_data\n",
    "    slope -- the slope of y in response to x our linear model \n",
    "    slope2 -- the slope of y in response to x2 in our linear model\n",
    "    intercept -- the intercept of our linear model (b in y=mx+b+error)\n",
    "    error_stdev -- the standard deviation of normal errors\n",
    "    \"\"\"\n",
    "    y_predictions = []\n",
    "    for i,x in enumerate(x_data):\n",
    "        x2 = x2_data[i]\n",
    "        error = np.random.normal(loc=0.0,scale = error_stdev,size = None)\n",
    "        y = slope*x + slope2*x2 + intercept + error\n",
    "        y_predictions.append(y)\n",
    "        \n",
    "    return y_predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's write a function to plot a scatterplot of both our actual data and lines representing the 'biological truth' according to our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_colored_scatterplot(x,y,x2,true_y=None,colors =[\"orange\",\"skyblue\"],markers=[\"^\",\"s\"],\\\n",
    "  labels=[\"Skin Cell\",\"Heart Cell\"],xlabel=\"x\",ylabel=\"y\"):\n",
    "    \"\"\"Plot a colored scatterplot of x vs. y\"\"\"\n",
    "    categories = np.array(x2)\n",
    "    colormap = np.array(colors)\n",
    "      \n",
    "    for i,color in enumerate(colors):\n",
    "            \n",
    "        x_of_this_color = [curr_x for j,curr_x in enumerate(x) if x2[j] == i]\n",
    "        y_of_this_color = [y[j] for j,x in enumerate(x) if x2[j]==i]\n",
    "        \n",
    "        curr_marker = markers[i]\n",
    "        curr_label = labels[i]\n",
    "        plt.scatter(x_of_this_color,y_of_this_color,c=color,marker=curr_marker,label=curr_label)\n",
    "        \n",
    "        if true_y:\n",
    "            x_of_this_color = [curr_x for j,curr_x in enumerate(x) if x2[j] == i]\n",
    "            true_y_of_this_color = [true_y[j] for j,x in enumerate(x) if x2[j]==i]\n",
    "            plt.plot(x_of_this_color,true_y_of_this_color,\"-\",color=color)\n",
    "            \n",
    "    plt.xlabel(xlabel, size = \"xx-large\")\n",
    "    plt.ylabel(ylabel, size = \"xx-large\")\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot the results of simulations from our multiple linear model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAESCAYAAADe2fNYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmclVX9wPHPdzYY9lVFBmQRRER2NUMTNUVNJSuttMRyDeyXUiZmpYmalmK2WZiIpGYuabilaKKpKILsoCzDNjCyCwwMs35/f5znztx757n7nbmzfN+v17zmzvM897nncZAv53zP+R5RVYwxxph0ycp0A4wxxjQvFliMMcaklQUWY4wxaWWBxRhjTFpZYDHGGJNWFliMMcaklQUWY4wxaWWBxRhjTFplLLCISGsRmS8iS0RkhYj8yjs+U0TWi8hi72u4d1xE5PcislZElorIyKB7TRCRNd7XhEw9kzHGGMjJ4GeXAWeoaomI5ALvisir3rmbVPXZsOvPBQZ4XycBDwEniUgX4DZgNKDAQhGZrap7In1wt27dtE+fPul9GmOMaeYWLly4U1W7x7ouY4FFXS2ZEu/HXO8rWn2Z8cAs730fiEgnEekBjAXmqOpuABGZA5wD/CPSjfr06cOCBQtSfwhjjGlBRGRjPNdlNMciItkishjYjgsOH3qn7vKGux4QkVbesZ7A5qC3F3nHIh03xhiTARkNLKpaparDgQLgRBEZAtwCDAJOALoAN3uXi98tohwPISLXiMgCEVmwY8eOtLTfGGNMXY1iVpiqfg7MBc5R1WJ1yoBHgRO9y4qAXkFvKwC2Rjke/hnTVXW0qo7u3j3mEKExxpgkZSzHIiLdgQpV/VxE8oEvA/eKSA9VLRYRAb4KLPfeMhu4XkSewiXv93rXvQbcLSKdvevOxvV6ElJRUUFRURGHDh1K9dGMp3Xr1hQUFJCbm5vpphhjGlAmZ4X1AB4TkWxcz+lpVX1JRP7rBR0BFgPXede/ApwHrAUOAt8DUNXdIjIV+Mi77o5AIj8RRUVFtG/fnj59+uBimkmFqrJr1y6Kioro27dvpptjjGlAmZwVthQY4XP8jAjXKzApwrkZwIxU2nPo0CELKmkkInTt2hXLZxnT8jSKHEtjYUElvey/pzEtkwUWY4wxaWWBpRG56667OO644xg6dCjDhw/nww/dsp4+ffqwc+fOOtd/8YtfTPgz7rvvPgYNGsSQIUMYNmwYs2bNinr9FVdcwbPPuiIIY8eOtYWlxjQy05bs5J5Fdb+mLan7d0ZDyWTy3gSZN28eL730Eh9//DGtWrVi586dlJeXR33P+++/n9Bn/OUvf2HOnDnMnz+fDh06sHfvXl544YVUmm2MybDy6sjH71m0k7wsmDysW4O2yXosqSgthtn9ofSzlG9VXFxMt27daNXKFRro1q0bRx55ZOjHlZZyzjnn8PDDDwPQrl07AObOncvYsWP5xje+waBBg7jssstwcx1C3X333fz5z3+mQ4cOAHTs2JEJE1zNzoULF3LaaacxatQoxo0bR3FxccrPZIzJvECAachejAWWVCybCiUbYPnUlG919tlns3nzZgYOHMjEiRN5++23Q86XlJRwwQUXcOmll3L11VfXef+iRYv43e9+x8qVKyksLOS9994LOb9//372799P//7967y3oqKCH/7whzz77LMsXLiQ73//+9x6660pP5MxpnEpr6ZBgosFlmSVFsP6R4FqKHw05V5Lu3btWLhwIdOnT6d79+5885vfZObMmTXnx48fz/e+9z0uv/xy3/efeOKJFBQUkJWVxfDhw9mwYUPIeVWNOEvr008/Zfny5Zx11lkMHz6cO++8k6KiopSexxjTOEUaOksny7Eka9lUUO83pFWu13LCn1K6ZXZ2NmPHjmXs2LEcf/zxPPbYY1xxxRUAjBkzhldffZVLL73UN0AEhtAC96msrAw536FDB9q2bUthYSH9+vULOaeqHHfcccybNy+l9htj0m/akp2+wSATuZN4WY8lGYHeSrWXXK8uT7nX8umnn7JmzZqanxcvXsxRRx1V8/Mdd9xB165dmThxYtKfccsttzBp0iT27dsHwL59+5g+fTrHHHMMO3bsqAksFRUVrFixIunPMcakT7TkPLgA09g0wiY1AcG9lYBAryVJJSUlTJgwgcGDBzN06FBWrlzJ7bffHnLN7373Ow4dOsRPf/rTpD7jBz/4AaeffjonnHACQ4YM4bTTTqNNmzbk5eXx7LPPcvPNNzNs2DCGDx+e8IwzY0xmTB7WjSkjupHbiNYji9/soeZu9OjRGr4eY9WqVRx77LHx3eD5AijdUvd4fk+4yHITwRL672qMqeOeRZGT7VNGdKNalaW7ynin+AAHK5UsIFYaZcqI5IbQRGShqo6OdZ3lWJJhwcMYkyap5FA27C/nzaID7DhURUHbHC7u15bHVu+N+p6GGDqzwGKMMRkUK4cSzVNr99ExL4uv9mnPMZ3yYtbnS7ankigLLMYY00hFGwYDOK1HG044LJ+crEaUYMECizHGNBn52UJplTKsaytW7i7j7eKDvF18MNPNqsMCizHGNBHd83M4s2dbDm+Tw5JdZZluTkQWWIwxpon49tEdkt7nqCHXu9g6lkYkUFQyYObMmVx//fVpufeGDRt48sknI55fvXo15513HkcffTTHHnssl1xyCdu2bYt6vyFDhgCuCOb555+flnYa09Ik8hd+vEFlyohudb4acpW+9VhagMrKyprAcumll9Y5f+jQIb7yla8wbdo0LrjgAgDeeustduzYweGHH97QzTWmWYk1nTjwF35VtbJw5yH+u+VAWj4zk+VeLLAkIRO1e3bs2MF1113Hpk2bALcKf8yYMcyfP58bbriB0tJS8vPzefTRRznmmGOYOXMmL7/8MocOHeLAgQMcPHiQVatWMXz4cCZMmMCNN95Yc+8nn3ySk08+uSaoAJx++ukAVFVVMWXKFObOnUtZWRmTJk3i2muvrZdnNKY5ijWdWFVZs7ect7YeYE9ZNQJEWrYe7/4qDVFoMpqMBRYRaQ28A7Ty2vGsqt4mIn2Bp4AuwMfAd1W1XERaAbOAUcAu4JuqusG71y3AlUAV8H+q+lp9tj2VeefRlJaWMnz48Jqfd+/ezYUXXgjAj370I2688UZOOeUUNm3axLhx41i1ahWDBg3inXfeIScnhzfeeIOf/exnPPfcc4DbPGzp0qV06dKFuXPnct999/HSSy/V+dzly5czatQo3zY98sgjdOzYkY8++oiysjLGjBnD2WefbfvZG5MG2w5W8uaWA2wqqaBr62wu7teB/h3zgMhTjWtqhGkZ5dLK95pMy2SPpQw4Q1VLRCQXeFdEXgUmAw+o6lMi8hdcwHjI+75HVY8WkW8B9wLfFJHBwLeA44AjgTdEZKCqVmXioVKRn5/P4sWLa36eOXNmzVbAb7zxBitXrqw5t2/fPvbv38/evXuZMGECa9asQUSoqKioueass86iS5cuKbXp9ddfZ+nSpTXbE+/du5c1a9YwcODAlO5rjIFHP/0cgLMK2jK8W2uy4/wH2z2LdkIjDSqQwcCirkhZifdjrvelwBlAIBHwGHA7LrCM914DPAv8Udw/m8cDT6lqGbBeRNYCJwLNqgZ8dXU18+bNIz8/P+T4D3/4Q04//XSef/55NmzYwNixY2vOtW3bNq57H3fccXU2FgtQVf7whz8wbty4kOPh+70Y0xKla1h8VPf82Bc1IRmdFSYi2SKyGNgOzAHWAZ+ramAzkSKgp/e6J7AZwDu/F+gafNznPc3G2WefzR//+MeanwM9m71799Kzp3vc4I3BwrVv3579+/f7nrv00kt5//33efnll2uO/ec//2HZsmWMGzeOhx56qKYntHr1ag4cSD25aExzUF/D4mmRhi3Tk5XRwKKqVao6HCjA9TL8yuAG8lh+fUSNcjyEiFwjIgtEZMGOHTuSbXLG/P73v2fBggUMHTqUwYMH85e//AWAn/70p9xyyy2MGTOGqqrIo39Dhw4lJyeHYcOG8cADD4Scy8/P56WXXuIPf/gDAwYMYPDgwcycOZPDDjuMq666isGDBzNy5EiGDBnCtddeW2cTMWNMZDmJpiNLi2F23S3EE5FXtT8tW6Ynq9GUzReR24CDwM3AEapaKSInA7er6jgRec17PU9EcoDPgO7AFABV/bV3n5rrIn1WqmXzm+KObpliZfNNcxatltek4zrzdvFBlu8uo02OcLAy8t+1IcUh50+EtX9l2qBNcSfnp6zqXvdgdj5cWAj5R8R1j3g0+rL5ItIdqFDVz0UkH/gyLiH/FvAN3MywCcC/vbfM9n6e553/r6qqiMwGnhSRabjk/QBgfn223YKHMSaW6av2UK3whcPyOfmIfB5Yujv2mwK701LN5DUDQgJD1IKUR/8A1v0NtHbyTrq2TE9GJmeF9QAeE5Fs3JDc06r6koisBJ4SkTuBRcAj3vWPAH/3kvO7cTPBUNUVIvI0sBKoBCY1xRlhxpimZdqS6JWH+3fIY+yRbenUKhtwIxqRRjooLYY5p0D3U2t3pw0LDHlZSnl13XG1vKoS2PTP0KACbsv0on+3rMCiqkuBET7HC3H5lvDjh4CLI9zrLuCuNLTJ1mekUWMZZjWmPkRL0OcIfLVvh5BjUUc65k+EkvVwYKMLKOACQ+GjMOQXkH8Ek8tug7V/pe7+kN6eka2PgGF3Q9/LISs7mUdKG1t572ndujW7du2ia9euFlzSQFXZtWsXrVu3znRTjElZpLxqJFniP3Tlm4etGf7S2qASUF3pei1Dfl4zRFZXNeS0hwtWQ277+BtZjyyweAoKCigqKqIpzhhrrFq3bk1BQUGmm2FMyhKdPpzQNORlU2uHv8JphRvOUvW/pt3RcOYb0PaoxBpYzyyweHJzc+nbt2+mm2GMaUkCvZXq8tpjWV4vv/qQey05UPR86DUBVaWNLqiABRZjTAuX0eUDfr2V4ABSXQ4HN0I7b11Lfk8Yfi/0+TZI4931xAKLMaZFizVstacssUmmkWZ/+doy26cnUl33dck6OPYmOP52yGmTUHsywQKLMcZE8N8tB1iwozSua4MXOUZbcxJ8Lq/f4tBe0fyJddejACBQeaBJBBWwwGKMaabSMcQ1f3spQ7u0YtWeMioSmD0fb6+lzjWbnvEJKgAaMvW4sbPAYoxpltJVIPKTz8v48XAXiKKufg8SHrhivm/fGlj0Eyjf6ZLx7QbC9rmNZiV9ohpv9scYYxqB4EAUaX963+OBYpJ7lkS8d6uqz+GZTvDycbDtLRj2azj/E9i3MvJK+ibAeizGmCYt0cWLwUorqwPr1uOS0CyxZVOhZAO8dxn0nBtySrSSEXse49Sd90LVXpBcOOtd6DzUXXBRUfyf0whZYDHGNGnJBJUqVT7ecYj3PjuIAiO6tWbRzkMRrw8fyoqZpwkqJsm+FSE7RPUteZMzt/2SbuWr2djmFI46+K7rnXwyDU6emfjDNEI2FGaMaXEeWfU5b245wBFtcvj+oE6M69UuoffHDGZh61PyqvbTtexTLt70Tb65+VtkaQXPFTzGcz0fq33PhsczujlXOlmPxRjTZKQy7BXuG/060L9DbvprA/qspp+8ul/IJV2yDvL14muhKqiXpFWweEqz6LVYj8UY02SkGlRaZwtf7tmWK4/txNEd80KCSqTEfMIi1f7qOAS+tgMuVThyXGhQCWgmvRbrsRhjWoxrB3cmP8c/gvjlTOKdXlyjtDhCaXugfA+09j5j03P+729CU4qjscBijGkxHlzmdnGMa5FkaTGQm9gHLPgRIUElKw/Gbwxd1FhaDEQpE5OhzbnSyYbCjDHNwpQR3bjimE5xXRvXkNqyqeRV7fc9lZcVtgy/dJubVrz5mboXL59a5741Q2VZeTBgohseC3w18anGYD0WY0wz8fLG/SzbXZaem3kJ+MlVD0F2PvS+BDb+wyXks/Kg/1XAn1ye5JPfwYq7XS0vBAgKOmG7QNZJ7Iefbyasx2KMaTKiJdhX7CnjpMPyuWFol9Q/KLhXUV3pkurBwWDdDFj7MLw4EJbcAl1GQVYOIUElIJA3Cb+v3/lmwnosxpiMS8c04mqFRTtLOb1n29RuFN6r8CsKWV0G86+BvK6AwKEdbpdH34aV1+ZN/MrkB59vJjIWWESkFzALOAKX7Zquqg+KyO3A1UBgj+Cfqeor3ntuAa7EZb7+T1Vf846fAzwIZAN/U9V7GvJZjDGpSdfalETuE5jxVSeRH22r4BoKOR2gssS93rci8qUDJtYGjWaQP4lHJofCKoEfq+qxwBeASSIy2Dv3gKoO974CQWUw8C3gOOAc4M8iki0i2cCfgHOBwcC3g+5jjGmBElmTUicY+W6+BSAuvzL4Zrh4L/S5jJqhr6w8yI7QU2oihSPTKWM9FlUtBoq91/tFZBUhFXXqGA88paplwHoRWQuc6J1bq6qFACLylHftynprvDGmUUt4Tcrs/nDWey6BHuhVaDVs+AcsmQIHi6DX12HEvdCun38SPjsbLipuVkn4ZDWKHIuI9AFGAB8CY4DrReRyYAGuV7MHF3Q+CHpbEbWBaHPY8ZN8PuMa4BqA3r17p/cBjDFRpbMUS70oWR+6MLHoRXj3YpdL6TIKvvgkHHZq7fXRkvDNKFeSrIzPChORdsBzwA2qug94COgPDMf1aO4PXOrzdo1yPPSA6nRVHa2qo7t3756Wthtj4tOogwpQs0Pjro/gvW/DOxe6oHL4GTBufmhQgehJeJPZHouI5OKCyhOq+i8AVd0WdP5h4CXvxyKgV9DbC4Ct3utIx40xLUheloYOayWiqgxeP9ntjSLZrgeycx4c2l73Xi0kCZ+sTM4KE+ARYJWqTgs63sPLvwBcBCz3Xs8GnhSRacCRwABgPq7HMkBE+gJbcAn+SxvmKYwxDW3KiCilWOZPdJtr+QxJ+e1DL1rF8D0zvZ+qgWzo+RXY8qILLIGKwzv+l1ywaqEy2WMZA3wXWCYii71jP8PN6hqOG87aAFwLoKorRORpXFK+EpikqlUAInI98BpuuvEMVY0y988YU98yklMJ3lzLZzV7SEL//QmwYZbPTbKg6AUXUMANb61/3N3T8idxy+SssHfxz4+8EuU9dwF3+Rx/Jdr7jDENq76CStRpxMEJdb9EemkxzDkFRj4YIajgvxgyUDCyGZZeqS8ZT94bY0y8IlYkjlSDK3hvk8U/h5JCl5iPKMpfieGlV0qLXT6nGeyfkm5p67GIyGjgclX9v3Td0xjT+CU67JWF724lMcXdWwkIBIJRD8In98P6GYET/vfI7xmalC8thtn9ajfkCi8YuWxqxHxOS5dSYBGRI3F5ksuBY3G/MQssxrQgiQ57Hd+1FUt2xa5CHDFJHxjSCk6mR5r+u/GfsP1t2LuCmsrDWXlw1LdjJ+SjBashP4+az2npEh4KE5F8EblMRF4HNgJ3AyXAz4GhaW6fMaaZObd3+9RuENxTCLioKGg/k63QpgAO/zKU73Ll7LNyqempBBLygUWRkURbq+KXzzE14u6xiMhYXM/k60A7YD0uMF2mqk/VS+uMMc1KyvvKx5j5Rdlu+O84V4Ll0HYY8VvYtwbWzwy7URwJ+UhrVQJDZM18T5VURA0sIjIAF0y+AxyFK53yEPAEsB8oBA7WcxuNMU3UGT3bMqpba7KzQieA+q0piUukmV/VFbD6z7Dsdqj43Ls4C/p8B/4zOkJRSerOHvMbZovWhkj3aeFi9Vg+AT4HngWeVNW3AydE5Kj6bJgxpn5ESrbHtQ98gk48LD/m5/rxXUHvN/Nr3Qzo+gVYcSfsXw35Ba6UvVZSs/Yk3Qn5FrKnSipiBZbAPzOUmr6jMaYpSHS2VjI9CFUlR6DSZ6JV+LBXtPvXSdT7raD36ylUl8EHl0OHY+DkWfDh1V5QwX+IKh0JeSvnElOswDIImIAbCrtaRDbjhsGeAA7Uc9uMMSlIJlBEC0bhPZr7Fu+MGFBS6vlEyqP47pOikNsRTn8dXh5CnanE4UNU0XobqtEXWJq4RU2lqepqVb0V6AOcDbyNm068DJiD+y2muA+oMaaxiBaMAudKKqp5ddN+36AS6x5xiTTj6sJ1MPxeyO0AkgPH/Ai+vgsu/hxW3AOV+2NXHA6ePRb8dc5HsRdYmrjFNStMVRV4E3hTRH4AfAOX1O8HPC4ik4EXgNmquqy+GmuMyax5nx1k3rZSKiPt7x4kqXphkfIonUfAyl+7lfPZbeDM/9aWsq/p4QDZ+XBhYeKzsywhn1YJT/5T1YOqOktVv4zryfwC12uZCiyO9l5jTNP2dvFByquV6thxJfVZXwHVZTD/ahdQjjzfJd43PlVbUmXxLamvKbH9VdIqpZX3qlqEWyB5t4ichFuFb4wxcQtJ8kfMo3SCsa/ASwOpyb1UHnCLHA9sDK1GnMyaEkvIp1XaaoWp6oe4rYWNMY1ArLUiU0Z0i74PfD2KWK7lgjWurtfKe1yQOOYGOO5WyOvoZooFeibVlbDhcUBrg0qADWFlXKwFkglvmKWqTybfHGNMukweFj1w1FdQSWp1vaob3lp8MxzcDL2+BsN/A+37uyGvF4a6lfSB3oxveXuPrSnJuFg9Fu+fBL77pgTToO8WWIxpYaLu6hjLzg9g4Y2w6wOXpD/573D4abXnl02Fg5uImhJONmlv6kWswHJuHPfoDPwEGEnEetTGmJYm0lBcXlUJlFa6MiyLp8DGJ6H1EXDSDOh7OWRl114cPOMrWrF9G/5qVKIGFlV9LdI5EWkH3Oh9dQJexM0QM8bUs3jLsiRbkyuR90Ua+gpZJDl/IhQ+4oapJBd2DoO9y925426FwVMgt13dmwTPEsvKg/5XuWGu0i2h19nwV6OScPJeRPKB64GfAl2B14FfqOpHaW6bMSaCSH/phx+PlWeJJNb7Ehr6Cl+bohWwewEUfBVG/Q7aRig7GGlXSBvyavQSKZufB1wHTAGOwK3C/4W3d33CRKQXMMu7VzUwXVUfFJEuwD9xa2Q2AJeo6h4REeBB4DxcReUrVPVj714TcPvBANypqo8l0yZjGpukFhk2NsumullcwSQH8o8MDSrhlYVt0WKTFTOwiEg2cBVwK1AAfAB8V1XfTPGzK4Efq+rHItIeWCgic4ArgDdV9R4RmYILZDfj8j0DvK+TcOX7T/IC0W3AaFyOZ6GIzFbVPSm2z5h6F2tIKx1BJdHhsMDQVsQcid/QV3BQQGtfV5W6XoaGBRatrDt0FV5Z2KoIN1mxphtPAH6J6z0sBn6gqi+n44NVtRgo9l7vF5FVQE9gPDDWu+wxYC4usIwHZnnlZT4QkU4i0sO7do6q7vbaPAc4B/hHOtppTH2KNqSV7HTgSO/rkJvF6T3bMqhTHm4AIAJvRfvk8JL1c06BU/8F//saDAzbryQ4KKi6hYtzz3NbAksOHH8HHPtjyGkT+TPDC0/aosUmK1aP5VFcL+BjXC2wYSIyLMr1qqq/TrQRItIHGIFbYHm4F3RQ1WIROcy7rCduo7GAIu9YpOPGmCBXD+5MblaslQP470kSOPbeZXXPBQeFtY/gdthQ2LMIel0Mox6ANjH+l4y0gZdpkuLJsQgwyvuKRYGEAos3u+w54AZV3RflX1J+JyKtsakz7VlErgGuAejdu3ciTTSmWYgrqPj1HNDaY/tWuOsKH4Wjr3G9l+6nBgWFstp7SQ607h47qERK0ttWv01WrMBybH1+uIjk4oLKE6r6L+/wNhHp4fVWegDbveNFQK+gtxcAW73jY8OOzw3/LFWdDkwHGD16tK23MU1eIA+T1gS/X88heJ+SAK3yei/rXQ/Gb42JVsYXICxJ3+zEWsfyaX19sDfL6xFglapOCzo1G7e52D3e938HHb9eRJ7CJe/3esHnNVwRzM7edWcDt9RXu41pLALBJG1BJVLJ+sDrYNXltb2XaOuiq0pd9eGTH418jSXpm51k1rHkAO2B/arhUz0SMgZXDXmZiATK7f8MF1CeFpErgU3Axd65V3BTjdfipht/D0BVd4vIVCCwjuaOQCLfmMYu2QWMAXOKStLXGN+S9eGVhpOw6ZnogcWS9M1OXIFFRPoBN+FmW/UOOr4JeBW4X1XXJfLB3vqXSIO+Z/pcr8CkCPeaAcxI5PONaQz8tvBNZDbYxzsOxXVdXIUhfUvWJxD1zl0CnYe6ns+/+7p9VMANa5V+ZvmSFiTmHzcRORdYAlwL5OK2JP6X9z2waHKxiJxTj+00xvj4/qBOMa+Jew/6SNv2XqrwlRXQw/tfPLejS8yHe98rhr5sqqsDFlBdntzmW6bJihpYvKm+TwC7gLNVtUBVz1HVi73vPXG9mN3AEyLSvf6bbIwJ6J4ffdBhyohu8QWVSA7thI8mwStDXRXikdMgu23dBY/g1qzsWQKFMwjt6VS7XI3tH99ixOqxXAW0Acap6ht+F6jq68A4oB1wZXqbZ4yJJHiFfLTzSakqh1XT4MWjYe1f4ejr3CZcg26Er22Bo3/gikIGy8pzM8WqffZKsV5LixIrx/Jl4MVYs8NU9RMRmY2bkXVPuhpnjAnlV/wxpR5JOFU3G2vRTVCy1g1/jbwfOg6u3WP+rPciz+Tatwr/vEy1zfJqQeJZx/LbOO81D1fx2BgTRbLrTgI9kHhL5vsKL/QYbM9i+HgybHsLOhzr9pg/MmhLpuAV+TaTy0QRK7B0AuKdorIT6Jhac4xpfpINJDkCPxleN1DEWzLfl1+5ltLPYOkvYN0j0KoLjP4jHH0tZAX99eC3It9meZkIYo3CtsJVIY5HFW6WmDEmSLLrVCqTqQ8RGK7yS5SHB4eSDbDi1/DiACicCcfc4PIoAyeFBhXwX5FvTATxrGMZICJfjOO6gak2xpimKtrwVEO655Nc6PUheavKmDwy7GRwcKgudzO9KvdDwXgY/lvoMMD/plbLyyQonsDyS+8rFsH2vDctVErDU/WgXFp5vRatLXcfsotjFVSWwCnPQO9v1L7RLwdjtbxMgmIFlh80SCuMMfEpLcatU45Dzd4oG+B/l7gpxMEkxyXqgwOLXw7GanmZBMUqQvnXhmqIMSY6N9wWZ1ABtyhRFaiGktV1z2tFaHA8S8+vAAAfzUlEQVSIlKC3GWAmQQ08AmyMCcjLcutS4l3gmNCwmla7Wl01+6NkQZ/v1C3VEhw0LEFv0iTh6sbGmOTV+wJHoOfB+Zy57eeEpjyrYfNzUPpb/4S7JehNGllgMSZO0WZ+pVr+Ph06VGxm7PY7GLzvBUpyDgfJdj2PgGgJd0vQmzSyoTBjovHWhUxbsiPqzK/Jw7px3eDOHNOp4Zdy5VWVMCXrASYWnszgA6/BkF/SLicrNKhAbcLdT7QEvTEJsh6LMdF4s6TKq6PvF/+bRTsRgXi2lU8X0SqO3/sUX9p+N1Rthz6XwbBfQ9teMPRXid3MEvQmjSywGBNJ8CypGKqBIZ1bcdqRbfjT8j310pzg4bbeB97jzG0/5/Cy5WzNP4F2p86GbifVy+cakygLLMZE4pd3iOL8o9oDkfMtqa7CnzysG+xf6yoPF70AbXrDF//BkUd9E6QBu0rGxBAxsIjIK0ncT1X1Kym0x5jGIXyWVALSPcsLgPLPYfmdsPr3bt+ToXfCoMmQk5/+zzImRdF6LF2wEi2mpUqwtwLR96oP9FYSLndfXQnrHoalv4SyXdDvChh2F+T3SKhtxjSkiIFFVb9Qnx8sIjOA84HtqjrEO3Y7cDWww7vsZ6r6infuFtwOlVXA/6nqa97xc4AHgWzgb6pqG42Z1HmzpKrIYVHnK9wK9hSGm6JNRY54butrsOjHbsvfw74EIx+ALuGVJY1pfDKZY5kJ/BGYFXb8AVW9L/iAiAwGvgUcBxwJvCEigWrKfwLOAoqAj0RktqqurM+Gm8Ynpc2v/FxUxLq95fx3ywF2lVVxVLtczixoy+OrP6//9Srb34G3zoGqUmjXD059DgousjyKaTIyFlhU9R0R6RPn5eOBp1S1DFgvImuBE71za1W1EEBEnvKutcDSwkRbYzJtyU7f4JLIBlwbSyqY8cnnNT/Xy4LIsl2w7Few+o+AQreT4cy3ILtVmj/ImPqV0DwVETlZRP4lIptFpEREDoZ9HUhDm64XkaUiMkNEOnvHegKbg64p8o5FOm5MjfooaZ/OoJKlFfDJg27DrdV/AvH+t9yzGMrrZ+qyMfUp7sAiIqcDbwPHA68DbYCXgNdw+Y1VwO9TbM9DQH9gOFAM3B/4eJ9rNcrxOkTkGhFZICILduzY4XeJMQ1Llf77X+fKwlPh4xugy2jo9Q1XigWsEKRpshIZCvsl8ClwAtAO+B7wF1X9r4gMAd4FPkilMaq6LfBaRB7GBS5wPZFeQZcWAFu915GOh997OjAdYPTo0TbbrZFIe24kic9KViLDYW5WmNas4O9+aCVnbP8lfQ+8ze68/nDaS9B5BLzY3wpBmiYvkcAyCpiqqodEpI13LBtAVZeLyHTg50DSxYVEpIeqFns/XgQs917PBp4UkWm45P0AYD6uxzJARPoCW3AJ/kuT/XzT8Bpy58V033PysG4wfyIUPhK63iUrD/pfVbd44/yJsPZhoDLkcJfKLdBlFCy7wwpBmmYhkRyLAIEB34Pe905B51cDg+O+mcg/gHnAMSJSJCJXAr8RkWUishQ4HbgRQFVXAE/jkvL/ASapapWqVgLX44bjVgFPe9eaFibWqvZpSyKvMUlJPMUbS4vh+d5uPUpYUPHe4IKHFYI0zUQiPZbNQD8Ar9eyGTgZeMY7P4zawBOTqn7b5/AjUa6/C7jL5/grQDJVAkwzMnlYt6hDXenurdQEsljFG1Vh3vegdHPkawLBwwpBmmYikcAyF7gQ+Jn38xPATSLSCjckdiVeDsOYTJg8rFvU1e+pSCrns/tj+GgS7IqSeszvaQHFNDuJBJb7gPdEpLWqHgJuBw4HrsCthv8HMCXdDTQmknQn4wP8dnlMSGkxLLkVCmdCVqvaDbci5V7C3zvnFDjrPUvYmyYr7hyLqhaq6hNeUEFVy1X1SlVtq6odVPVyVd1ff001zVG8+737STSoxFtd+J5FO+PPyXgbgVH6GVSWwvK73HqUDU/AgOsArd1wKzDLq/SzyPfz9n+xacamKYu7xyIifwZmqOqCCOdHAVeq6sR0Nc40f/VSCTjSZw2sqOkN3PNJ9D/6cQetZVOhZD3Muxz2fQoHN0Gvr8Pwe2HV/dRZVhVtllfw/i/B04ytF2OamERmhV0HDIxy/mjg2tSaY0z9yMsipDeQ6t4ogPsLv/ARQOGzOZDbAc6cC6c+C+37Jz7LK7iicvDiSOvFmCYmnXvedwHK0ng/Y1I2ZUQ3pozo5norQb2ByQMrY+dSog1ZHdgMc06rDRySDZ1Hwoffr33fRUVw0VZXSPKiYrhU3Zdfsj58/5fAsNmeJaG9mGhtMqaRiDoeICJfBE4JOnS+iBT4XNoZ+A5W/NHUs6QT9n69gViLDhdPgR3/Cx2CqiiBlffCqvug+lDttVrl8iqBNSmBewf3NqJ9nt/+L1oF712WeLuNybBYOZazgNu814pb2f6tCNduA25KU7tMCxRPeZekgkqk3sCQXxD1f4H1j1MTKEb/AdbPgiU/c/drdzQc2AhaEfQGL0lfc2/1z5n4iTRstm9F6M9W4sU0AbECy5+BZ3Gr7pfiphO/HHaNAiXAZlW1GlwmadEWN6a0PiVSb2D5VPLyfuUfzKr2UxMo1v3N9Vw+XwZdT4JTnoN3Lw4LKnXvjWr8vQ2/4TG/cjHWazFNgMQbC0RkHLBEVZv8IO/o0aN1wQLfyW0mg+pjceOUEd3g+QIo3VL3ZOsjIKeNG+r6z2j/awJy2sKJ0+Gob9WWtQ8oLYbZ/aAqaGgsq7X7Hjxclp0PFxbG39uI1G5bVGkyREQWquroWNfFPd04sBWwd/POQB9cb2WjqtqmEabR8S27Ejx1d9kdsPavrgcQfM2+1fDy4Nr1JwDVVdDxOLdGJXzar1+PKHxYCxLvbVjwME1UQjtIemtVHgDGhB1/D7hRVRemsW3GJCSuFfOBZPriKbDpn4TkP1p1d0NPC38UGlTAXffeZf6JeL/8CD7ja4GpxjaMZZq5RBZIjgTeweVb/o6bASbAscDFwNsicqqqLqqPhhoTTVzrUoIXIK5/HLKCNtSafw0c2ODyKFl5dd8bnEgPT6Bbz8KYEIn0WO7EVS8eo6obg0+IyC+B94GpwPnpa55p7NK1UVfhvnKECNt/RpFQXa+QIasqN7wFLmhseRHa9IZTnnEr5yVsc9LgRLol0I2JKpHA8kXgt+FBBUBVN4nIX4CfpK1lpklIZiZXcNDZeaiSt7YcYN2+CjrlZXFGz7YM6JiHiKSlyGTIPXLvgGPuAKB9eRGT1o2ovVCyocc50PsbdW8SbbqyTfs1po5EAksuUBrl/EHvGmOiKq+G0spq/ld8kEU7D5GXLZx+ZBtGdc8nJ6u2pxAIPknPFistprw69I9kllYwfM9jnLrzN6HXahVsDZ9J74kyXdl6LcbUlUhJl6XAFUHbEtfwjl0BLElTu0wz9+Cy3Xy88xAKVFcrJx3eJiSopMWy0Npa/Ure4PuFp3H2tlvY1up4dzC/Z/RSK2A7OxqToER6LL8GXgAWi8hDwCe4IfHBuAKV/YGvpr2FptmriJFYycvyH3LLy4ryxsDw1cA76Fb2CWds+yX9DrzF7ty+PFcwizXtzmHKyO7xNdCS88YkJJF1LLNFZAJwv/cV+L9agJ3AFar6YvqbaFq6OpMA5k90608GXAdEGIpaNhWqqzjrs58yYs8syrPa8eZhU1nY5ftUi8+sL2NM2iS0jkVV/y4i/8Ttdd/HO7wBmKeqPivCTHNSXzs2JiTSniXBDmyEdQ+DVjJizywWdb6Cd7vdRGlO14w02ZiWJlZ14xnAX1X1w8AxL4C87X0lzbv3+cB2VR3iHesC/BMXtDYAl6jqHhER4EHgPNwkgStU9WPvPROAn3u3vVNVH0ulXcZpFEHET3AivboSXjwGLvjUBRdVl/f4YAJoJQCzjnqFz9qMrHObtOzHYozxFbVWmIhUA99R1SfT/sEiX8IVr5wVFFh+A+xW1XtEZArQWVVvFpHzgB/iAstJwIOqepIXiBYAo3FDcwuBUbFKzFitsNjqo25XPKKuf/GryQXQdwIMugEW3gjb54Lk1AQWX1Zry5ikpL1WWLqp6jsi0ifs8HhgrPf6MWAucLN3fJZXPfkDEekkIj28a+eo6m4AEZkDnAP8o56bb8JccUwn3igqoehAlL/Q4xC1l+Q37RdcOfv1s6BVFxj9Jzj6GsjK2B9tY1q8xvZ/3+GqWgygqsUicph3vCewOei6Iu9YpON1iMg1wDUAvXv3TnOzzcxPP4/rusBK+bh7RMFFI31rcgEodBoGX34L8jrH2WJjTH2JJ7AcJiL94r2hqham0J5I/BY4aJTjdQ+qTgemgxsKS1/Tmo50lV9JVlJ5jeAdGC8qcnmUtdPho+tCr9u/GqpsZ2xjGoN4AktgenG8spNsC8A2Eenh9VZ6ANu940VAr6DrCoCt3vGxYcfnpvD5zVq08isNIeHgFT4DrOcFsOJO2PEehFcWs5XwxjQa8QSWF3Cr7hvCbGACcI/3/d9Bx68Xkadwyfu9XvB5Dbjb2x8G4GzglgZqa6OQrl5IppL1UQXnVKoOwdxzoVU3yO0EFWFDb1aS3phGI57A8lw9zQr7B6630U1EioDbcAHlaRG5EtiEK8cP8ApuRtha3HTj7wGo6m4RmQp85F13RyCR31LUZy8k4XwIkVfJx3ufmuGy0mIonBGUU/F6J0ecDWOeiLs9xpiGl8lZYd+OcOpMn2sVmBThPjOAGWlsmvEk04uJ1EuKdq86pe+1Gt7/LlT75Ew2/RNG3l+7KDI4uW+Vho1pFGyZmEmbtCw6LJoN/2wD2970Px/IpQQEJ/eNMY2CBZYWpL5Xm6c0s6xkA7x9Ibwz3vVUup0CbfvCuYshu3XotYWPQulndZP7pZ+l0nxjTJpEHQpTVQs8TVj48FNeVt1hp7Qm7Us/8x+OKi0m0lY9eVUlsHgafDItdLX8znnU7DMfaS8U1dpzNivMmEbDAkcTl0gvxC+pnpumLVDyqvZHHo5aVve4aBVDP3+ca9adCCt/DT3Pd4GiRhWgbp95v71QNv/Lf1dH67UYk3GNbeW9SZDf8FM8vZBqVZbsOkRutlBRqRzfpRXLdse/wLCm5xNcvys7v261YW+4Kq//TZRntweg94H3OHPbzzm8bDlbW4+kXet2kN0W8Il8WXnQ/6q6PZHAHvTBrNdiTKNggaWRq4/V8hv2lfPmlgPsOFRFQdscLunfjiPa5PDp52VxTVMO6SUFrzXx+4vdOz95tX/xhiMPfexeHNjo/2GR9pePtqujBRZjMsoCSyOXyDqVeErdP7NuL+v2VdAxL4uv9m3PMR3zcLsSpLAyPnw4KjgIRKzvFS5Kw/0CllUnNqbRshxLIzZtSWKJ9Xh6G5tLKhl7ZBuuPrYzgzq1qgkqSfGrNhw8Hbi6Eobc6lbLI9Dv+9A6ibUmtr+8MU2K9VgasViBItEZXVnAtYM70zY39N8TSQ+3RRuO6nkhLPox7F0Bh50GIx+ALiPcNYH8SPB7I+VSjDFNjgWWFqQa+MPy2oo3gcCRdFkYv+GovZ+4gDL3HGjXH079FxR8FYJ7RpYfMaZZs8DSgqW1qnHZLlh2O6x5CHLawojfwsAfQnarutdafsSYZs0Ci0lNVTms+TMs+xVU7oP+18DQX0Hrw2K/1xjTLFny3kQ3u3/oosPSYnfsYDEUvQivDIGPb4SuJ8C5S+D4X8LrJ9tCRWNaMAssjVh91/aKS3iBx2VToWQ9vHYivHMhSBac9hKc/hp0GmJFIY0xNhTW2MSzFiXdIu2hkldVQk2BxyG/gLIdsG46oFBaBEPvhME/hSyvDlh4UcjwRY3GmBbBAksjk0xQCS4smeimXBBhSnFwyZTqSnjna7BnoVunAiC5ULq1NqhA7FX4xpgWwQJLE5focFlcpWDCV9RrBeyaR8jIqVaE9kriWYVvjGkRLLA0IQM65nFGz7Z0bpWd1Pvr7NQYybKprpcSwmeFfnCvJNoqfOu1GNOiWGBpBOLNq3y9X4eY10TMl8Tbszm4FdY/Fro3CuD2nNfQQ8GLGm3RozHG0ygDi4hsAPbjNuWoVNXRItIF+CfQB9gAXKKqe8QVu3oQOA84CFyhqh9not3JSmeyfvLAiuT2gK8shU/uh5X3uKBy7E1w3K2Q1zG+99uiR2OMp1EGFs/pqhqciZ4CvKmq94jIFO/nm4FzgQHe10nAQ973ZiXuHkfwdN9YPYXSYvjPSVBVCtl5Lhnf6+sw/F5o3z/VJhtjWqjGHFjCjQfGeq8fA+biAst4YJaqKvCBiHQSkR6qWpyRVnriLeyoqnUvChJ3XgQSn+770SQo3exet+oGZ86Fw0+L//OMMcZHY1iC50eB10VkoYhc4x07PBAsvO+BmiE9gc1B7y3yjmVUPIUdi0oqeGz13vR9qN90Xz8HNrnpw0XPBzVsP3Q4Jn1tMca0WI01sIxR1ZG4Ya5JIvKlKNf6bShSpxsgIteIyAIRWbBjx450tTMpn5dV8cL6fTy+Zi8HKhJIsATKqfiVS4k03Tf42ooSWPILeOkYl2wPphW2Wt4YkxaNMrCo6lbv+3bgeeBEYJuI9ADwvm/3Li8CegW9vQDY6nPP6ao6WlVHd+/evT6bH9PDq/awdm85pxzhNtyKlD+pczxauZRo0321GgpnwksDYcWd0GMcED5luRrWPmI1vowxKWt0ORYRaQtkqep+7/XZwB3AbGACcI/3PbCl4GzgehF5Cpe035vp/Eosgzq14rQj29Ahz/3lHteWwLHyJ5Gm+258GnZ9CLsXQteT4IS/wvuX+UwnBrTM1p0YY1LWGHsshwPvisgSYD7wsqr+BxdQzhKRNcBZ3s8ArwCFwFrgYWBiwzc5MRf0aV8TVOIWK39yURFcqrVfF65zM7zKd8KhbTD6j3BoO2x+Dir3E3GP+c3/Svh5jDEmWKPrsahqITDM5/gu4Eyf4wpMaoCmJSRXoMJnwldSFYv98ifrZsCWl2HcB6E9l4p9sPwu+PR3IDlw/B1w7I/h45/AgQ2wYZO7LjsfLiyEZXfUbhOclQe9vpZEA40xplajCyxNVbyr58ur6xaKjFm/yy9/Ul0OBzfWDl1VV7kAsfQXrmfSdwIMuwva9AwaRtPaIpJaBYunwKZ/Wn0vY0xaNcahsCYpldXzgWAzbUmEysR++ZPAUFbho7DxGfjPSJh/LbQfCOM+gpNnuqACXu2vqrC3l8P6xyMn/I0xJknWY0mDwn3hf+knJ2JwCi+XEihpX10OVYfgvUugbR845RmXV5GgGdiB3opW+Ny4yj/gWH0vY0wKLLCkaN3ecp4p3NdwH1haDIUzgnow6nIpZ74F7frUvd5vGC1Yfk+r82WMSSsbCktR3w65nNu7XcN8WHUF/O8SqC4LPS5ZsOq3/u/xHUbDBZRL1YKKMSbtrMeSoiwRhnVtzaubSurnA0qLXbXi4+9wixv3fVL3mmjDVxY4jDENzAJLChpkf/oFN0JJIcz7jkvMf2k29Dy/bh5lzilu1bzN5jLGZJgFlhRECypTVh0OA67z7UVEq3wMuEDx+slupfzmp90xyYUz5kDb3rVvCASU7qfGXyrfGGPqmQWWehO5dH3UNSslG+ClQS6PcmAjLg1W7XooK+8NDRzLpkLJeu+6OEvlG2NMPbPkfX1KZE2IKmx+Hl4ZGpac97o24dWKIy16tDUoxpgMs8BSn/xK1/vZsxjePAP+9zWvjlcEwYEj0qLHeD7PGGPqkQWW+hYcDML3Uyn9DD68Cl4dCXuXQZcTo98rMPsr2qJH67UYYzLMAksKIu6jUhXU6wgEA6jdT2XpbbDibnhxABTOgtwOcOrzrrR9uKw8GDCxtmrxRUXRFz0Gf54xxmSAxNpzvTkaPXq0LliwIH03nD8R1v2ttgch2fDVotAkemkx/Luvlz8RQCG7jdt0q+jf0OFY2LfC//7hq+OfL4DSLbGvM8aYNBKRhao6OtZ1NissVX7DUoHKwSfPrD320fWhK+Db9HHVibfMBqrjDypgwcMY06hZYEmVXxIdYMPjMOhGePtC6HQ8bH056KTCwQ3ey6D39r0CTn60HhtrjDH1z3IsqdoyO3IS/c2z4OAm2PoKcf2nXv93m9FljGnyLLCkorQYsltB6wgLEst3eC+UiFsBh6iCxbekqXHGGJMZFlhSEZjl1etrbsbWWe9CFy+v1aq7K2cP1PxnHjARjv6Bm+kVyaZn6rPFxhhT7yzHkozSYnjtC24LYKph3SNwsMgNi+UfCaN+D4tuAq303uD1VtbNgLxO/mXsa1RbMUljTJPWbHosInKOiHwqImtFZEq9ftiyqS53Eii9Ul3m8ihDboMLVsPeVbjhrzDV5bW9m0vVv/diCxyNMU1cswgsIpIN/Ak4FxgMfFtEBtfLh5UWux4KEBI8JMdVM85pG3lzrUChyECC3u86W+BojGnimstQ2InAWlUtBBCRp4DxwMq0f9KH14FGCBqBsvWBdSbBe9MHBHokwdcZY0wz0ix6LEBPYHPQz0XesfTa8T5sne1/zq+nYT0SY0wL1Fx6LOJzLCTJISLXANcA9O7d2+fyOKx/nJpyLAFZedD/KtsW2BhjPM2lx1IE9Ar6uQDYGnyBqk5X1dGqOrp79+5Jfsrz1EnKWw/EGGNCNJfA8hEwQET6ikge8C0gwphVCgouCprFFbQ2xXomxhhTo1kEFlWtBK4HXgNWAU+raoSqjkkKFJusyZkErU2xMizGGFOjWQQWAFV9RVUHqmp/Vb0r7R8QaQ+U6nJbd2KMMUGaTWCpd/GuTTHGmBbOAku8Liqy1fLGGBMHCyyJsrUpxhgTVXNZx9JwbAaYMcZEZT0WY4wxaWWBxRhjTFpZYDHGGJNWFliMMcaklQUWY4wxaSWqPjsdNnMisgPYmMItugE709ScpqClPS+0vGduac8LLe+Z0/G8R6lqzCq+LTKwpEpEFqjq6Ey3o6G0tOeFlvfMLe15oeU9c0M+rw2FGWOMSSsLLMYYY9LKAktypme6AQ2spT0vtLxnbmnPCy3vmRvseS3HYowxJq2sx2KMMSatLLAkQETOEZFPRWStiEzJdHvqi4hsEJFlIrJYRBZ4x7qIyBwRWeN975zpdiZLRGaIyHYRWR50zPf5xPm99ztfKiIjM9fy5EV45ttFZIv3e14sIucFnbvFe+ZPRWRcZlqdPBHpJSJvicgqEVkhIj/yjjfb33OUZ27437Oq2lccX0A2sA7oB+QBS4DBmW5XPT3rBqBb2LHfAFO811OAezPdzhSe70vASGB5rOcDzgNeBQT4AvBhptufxme+HfiJz7WDvT/frYC+3p/77Ew/Q4LP2wMY6b1uD6z2nqvZ/p6jPHOD/56txxK/E4G1qlqoquXAU8D4DLepIY0HHvNePwZ8NYNtSYmqvgPsDjsc6fnGA7PU+QDoJCI9Gqal6RPhmSMZDzylqmWquh5Yi/vz32SoarGqfuy93g+sAnrSjH/PUZ45knr7PVtgiV9PYHPQz0VE/6U1ZQq8LiILReQa79jhqloM7g8wcFjGWlc/Ij1fc/+9X+8N/cwIGt5sVs8sIn2AEcCHtJDfc9gzQwP/ni2wxE98jjXXKXVjVHUkcC4wSUS+lOkGZVBz/r0/BPQHhgPFwP3e8WbzzCLSDngOuEFV90W71OdYc3nmBv89W2CJXxHQK+jnAmBrhtpSr1R1q/d9O/A8rnu8LTA04H3fnrkW1otIz9dsf++quk1Vq1S1GniY2mGQZvHMIpKL+wv2CVX9l3e4Wf+e/Z45E79nCyzx+wgYICJ9RSQP+BYwO8NtSjsRaSsi7QOvgbOB5bhnneBdNgH4d2ZaWG8iPd9s4HJv1tAXgL2BoZSmLiyHcBHu9wzumb8lIq1EpC8wAJjf0O1LhYgI8AiwSlWnBZ1qtr/nSM+ckd9zpmcyNKUv3MyR1bjZE7dmuj319Iz9cDNFlgArAs8JdAXeBNZ437tkuq0pPOM/cEMCFbh/tV0Z6flwwwV/8n7ny4DRmW5/Gp/5794zLfX+kukRdP2t3jN/Cpyb6fYn8byn4IZ1lgKLva/zmvPvOcozN/jv2VbeG2OMSSsbCjPGGJNWFliMMcaklQUWY4wxaWWBxRhjTFpZYDHGGJNWFliMaSAiMldE5ma6HcbUNwssxiRJRDTOr5mZbqsxDSkn0w0wpgn7btjPX8OtbP4JsC3o+Drv+9kN0ShjMs0WSBqTJiJyO3AbMEBV12a4OcZkjA2FGdNA/HIs3lDZ4yLyZRH5SERKRWS1iFzsnR8uIv8VkQMislVEJke490Ui8p6IlHhfb4rIFxvgsYypwwKLMZk3DHgCeBm42Tv2lBdcXgMWADfhanzdLyIhQ2oicgPwL9xGXrfgek09gLdEZEyDPIExQWwozJg0iTUUFuitqOrYoGMKVAMnqLf7n4gch6tAq8B4VX3RO94VF1xeUdWve8cKgELg96r6k6D7tgNWAptU9ZR0P6sx0Vjy3pjM+zAQVABUdYWI7AUOBIKKd3yXiHwKHB303q8BucCTItIt7L5v4krBt1HVg/XYfmNCWGAxJvM2+hz7HFfm3u94/6CfB3nfF0a5f1fAAotpMBZYjMm8qgSPB28pG8iTXgSURLh+RzKNMiZZFliMadoCuZytqtqkdnk0zZfNCjOmaXsWqARuF5E6/1AUkcMavkmmpbMeizFNmKpuEJEfAw8CC0XkaWA7UACMxc04Oz1zLTQtkQUWY5o4Vf29N1vsx7j1Lq1xif8PgUcz2TbTMtk6FmOMMWllORZjjDFpZYHFGGNMWllgMcYYk1YWWIwxxqSVBRZjjDFpZYHFGGNMWllgMcYYk1YWWIwxxqSVBRZjjDFpZYHFGGNMWv0/hf4t+fwGOeQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from random import choice\n",
    "\n",
    "#simulate observed data points\n",
    "N = 250 #number of data points\n",
    "\n",
    "\n",
    "#x values will be a number\n",
    "x = [i for i in range(N)]\n",
    "\n",
    "#x2 values in this case will be a category\n",
    "#represented by 0 or 1\n",
    "x2 = [choice([0,1]) for i in range(N)]\n",
    "\n",
    "x1_slope = 12  #increase in y with each x\n",
    "x2_slope = 300 #difference between cyan and orange\n",
    "intercept = 0 #intercept\n",
    "error_stdev = 100 #measure of 'scatter' or 'spread' of points around the model prediction\n",
    "\n",
    "\n",
    "observed_y = linear_simulation(x,x2,slope=x1_slope,slope2=x2_slope,\\\n",
    "  intercept=intercept,error_stdev=error_stdev)\n",
    "\n",
    "true_y = linear_simulation(x,x2,slope=x1_slope,slope2=x2_slope,error_stdev=0)\n",
    "\n",
    "plot_colored_scatterplot(x=x,y=observed_y,x2=x2,true_y=true_y,xlabel=\"Time\",ylabel=\"Total DNA\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stop and Think**: compare this linear simulation function to the one earlier in the chapter. Imagine running the single-variable linear model with `x = [0.0,1.0,3.0,4.0]`, `slope = 2.0`, `intercept=10.0`. Now imagine running this function with `x = [0.0,1.0,3.0,4.0]`, `x2 = [13.0,12.0,11.0,10.0]`, `slope = 2.0`, `slope2` = 0.0, `intercept=10.0`. Do you think the the two functions return the same or different y values? Why? \n",
    "\n",
    "**Try it out**. Test your prediction by generating y values using both the single and multiple linear models using the parameters described above. If you like, you can set up a new Jupyter Notebook in which to run this. Do the two sets of predictions differ? \n",
    "\n",
    "**Try it out** Now, repeat the exercise setting slope2 in the multiple linear model to `5.0`. Do you think the y values will differ? Do they?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting a Multiple Linear Model from Data in Python\n",
    "\n",
    "Look at the graph up above. Let's say we had measured this data and didn't know what model had produced it. We could use ordinary least squares model fitting to try to figure out how well the data fit to a multiple linear  model in which y (`Total DNA`) is determined by the amount of time that has passed (time or `x1`) and what type of cell we are studying (skin vs. heart cells, as represented by whether `x2` is 0 or 1).\n",
    "\n",
    "To fit this multiple linear model in `statsmodels.py`, we don't have to change much - we simply adjust our equation from `y ~ x` to `y ~ x + x2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.988\n",
      "Model:                            OLS   Adj. R-squared:                  0.988\n",
      "Method:                 Least Squares   F-statistic:                 1.006e+04\n",
      "Date:                Thu, 15 Jul 2021   Prob (F-statistic):          2.31e-237\n",
      "Time:                        16:19:29   Log-Likelihood:                -1499.0\n",
      "No. Observations:                 250   AIC:                             3004.\n",
      "Df Residuals:                     247   BIC:                             3015.\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept      8.4306     14.023      0.601      0.548     -19.189      36.051\n",
      "x             11.9593      0.086    139.502      0.000      11.790      12.128\n",
      "x2           294.2004     12.431     23.666      0.000     269.716     318.685\n",
      "==============================================================================\n",
      "Omnibus:                        3.103   Durbin-Watson:                   2.265\n",
      "Prob(Omnibus):                  0.212   Jarque-Bera (JB):                3.484\n",
      "Skew:                           0.016   Prob(JB):                        0.175\n",
      "Kurtosis:                       3.577   Cond. No.                         376.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "import statsmodels as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import pandas as pd\n",
    "\n",
    "#Set up our data table\n",
    "#Note that x, x2, and observed_y were generated in the previous cell. \n",
    "#(if you get a NameError, back up and generate them)\n",
    "df = pd.DataFrame({\"x\":x,\"y\":observed_y,\"x2\":x2})\n",
    "df.set_index(\"x\")\n",
    "\n",
    "#Here is where we actually fit our linear model to our data\n",
    "model = smf.ols(formula='y ~ x + x2', data=df)\n",
    "results = model.fit()\n",
    "summary = results.summary()\n",
    "print(summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interpreting our results.\n",
    "\n",
    "We interpet the results as before. Now, however, the `coef` column represents the value that the regression analysis inferred for the slope associated with x **and** x2. So for example, I got an x coefficient of 11.87 (yours may be slightly different due to differences in the random simulation). This implies that each time x goes up 1 unit, y goes up roughly 12 units. The coefficient labelled 'Intercept' is again just the inferred intercept (e.g. the b in the equation for a line y = mx + b). The effect size (R<sup>2</sup>) and p-values are interepreted as before.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The same linear model equation can handle both binary and continuous `x` variables\n",
    "\n",
    "<img src=\"./resources/linear_scatter_multivariate_with_two_regression.png\" width=\"400\"  description=\"Another scatter plot with two sets of points — one orange and one cyan — trending upwards and to the right. Now each color of points has it's own regression line, and they seem to match the data pretty well.\">\n",
    "\n",
    "We started our investigation of multiple linear models with orange and blue colors as our second indepedent or `x` variable. This was mostly because this example is easy to draw. However, as you can see from the above equations, there's nothing special about the math that requires our second x variable to be categorical. We can still just multiply `y` = `x1`*`slope1` + `x2`*`slope2` + `intercept` if x2 is a continuous variable rather than a categorical one.\n",
    "\n",
    "So using what we've already worked out, you could *also* build models that use two or more continuous x variables. For example, you might model the mating success of a male elk based on several continous factors like it's size, it's age, the number of competitors, etc. Each of those factors would be its own x variable (age could be x<sub>0</sub>,number of competitors could be x<sub>1</sub>, etc), and each factor would have it's own slope or coefficient that says how it affects mating success (e.g. the effect of age would be represented with &beta;<sub>0</sub>, the effect of number of competitors could be represented with &beta;<sub>1</sub>, etc). \n",
    "\n",
    "> To model extra independent predictor variables in linear regression, we just add extra slopes and x's to our equation for a line. If these predictors are categorical, we convert them into binary variables that can be represented by 1's or 0's."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Summary**. By extending the familiar equation for a line with an error term (&epsilon;) that typically represents errors drawn from a normal distribution, and adding extra variables and slopes, we have a pretty powerful way of representing the independent effects of several different parameters on a response variable that we care about. The major limitation of our approach so far is that the effects must be *independent*, because mathematically we are just adding together each factor's effects.  In the next section, we will elaborate on this to allow for situations in which x variables *interact* with one another to influence the y value. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiplication can model interactions between two predictor variables\n",
    "\n",
    "All of our above math stemmed from our observation in the cartoon that one line for some orange points looked offset from the line for other points. But, critically, in that example both sets of points seemed to have the same slope. What if instead we have data that looks more like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./resources/interaction_effect_no_lines.png\" width=\"400\"  description=\"Another scatter plot with two sets of points — one orange and one cyan. Now, however, the orange points trend downward while the cyan ones trend upward as the x value increases.\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the orange points have a positive slope, while the cyan points have a negative slope. That is, unlike in our last example, the color of the points doesn't just change the y value by a fixed amount — it actually reverses the slope! \n",
    "\n",
    "This is an example of an **interaction effect**. If you plot your data by category as in the plot up above, and your different categories appear to have different slopes, that suggests you should test for an interaction effect in your data. It may not be there — the other possibility is that any apparent differences are just due to chance — but it's worth testing. We can think about how to test the statistical significance of interaction effects in several different frameworks latter, but for now let's try to figure out how to incorporate them into our linear model. \n",
    "\n",
    "We can think of the situation depicted in the plot up above as an [interaction](https://en.wikipedia.org/wiki/Interaction_(statistics)) between the `color` and `x` variables. Here, the word interaction just means that the outcome of having different color and x variables can't be predicted by simply adding up the independent effects of color (orange is a little higher) and the `x` parameter. Since the math we developed in the previous section just added up terms, it can't handle interaction effects. However, a very simple extension of that math can handle this situation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can model interactions using multiplication. Specicially, we will add an extra term to our linear model that includes the interaction of x value and being orange, times `x` and times the `is_orange` value: `interaction_of_x_and_orange` * `x` * `is_orange`. The `interaction_of_x_and_orange` term will be inferred from the data using ordinary least squares, and the product of these values will be added into our linear model just like all the other terms.\n",
    "\n",
    "Before we add this to the equation, let's consider what we would get if we just modeled the cyan points. We'd get a negative slope, because y goes down as x increases for the cyan points. Let's imagine that the slope of the line describing the cyan points is -0.75 (just my guess based on the cartoon for now). Now we want the slope to be positive  if the points are orange, since those points are going up not down. If the `interaction_of_x_and_orange` is 2, then when a point is orange (and therefore `is_orange` is `1`), a prediction based on a positive slope is incorporated into the `y` value. Since the interaction term is bigger than the main slope (|2| > |-0.75|), our interaction term would overwhelm the negative trend in the case of orange points only, causing the orange slope to go up rather than down. This is what we wanted! Here's how all that would look if we wrote it out with descriptive variable names. The equation is long, but there's nothing super complex going on here:\n",
    "\n",
    "`y = x * effect_of_x + is_orange * effect_of_being_orange + interaction_of_x_and_orange * x * is_orange + intercept + error`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this is getting annoyingly long to write out, a more mathematical representation becomes a bit easier to read (consult the explanations above if this notation is still unfamiliar):\n",
    "\n",
    "y = &beta;<sub>0</sub> + &beta;<sub>1</sub>x<sub>1</sub>  + &beta;<sub>2</sub>x<sub>2</sub> + &beta;<sub>3</sub>x<sub>1</sub>x<sub>2</sub> + &epsilon;\n",
    "\n",
    "As with our 'dummy variable' approach to handling categorical data, there are several extensions that are possible to this multiplicative approach to handling interaction effects. If we are testing for more than one interaction effect, for instance, we can add more than one term. For instance, if studying the effects of age, diet, and microbiome diversity on the ability of Woodrats to put on weight, we might have a `diet*microbiome` interaction term to test whether microbiome diversity matters more for some diets vs. others, and also an `age * diet` term to check whether, for example, diet matters more for very young woodrats vs. mature ones. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating a multiple linear model with an interaction effect in Python\n",
    "\n",
    "We can model **interaction effects** by assigning a slope to the product of our two x values. Let's update our `linear_simulation` function to have a third slope that will be multiplied by `x*x2`. Thus, the results will be identical to those of our multiple linear model, except that we will add the product `slope3*x*x2` to them.\n",
    "\n",
    "First, we modify our simulation function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_simulation_with_interaction(x_data,x2_data,slope,slope2,\\\n",
    "  interaction_slope,intercept = 0,error_stdev = 0):\n",
    "    \"\"\"Return y data given x data and parameters\n",
    "    x_data -- a list of x data points (ints or floats)\n",
    "    slope -- the slope of our linear model (m in y=mx+b+error)\n",
    "    intercept -- the intercept of our linear model (b in y=mx+b+error)\n",
    "    error -- the standard deviation of normal errors\n",
    "    \"\"\"\n",
    "    y_predictions = []\n",
    "    for i,x in enumerate(x_data):\n",
    "        x2 = x2_data[i]\n",
    "        x_interaction = x * x2\n",
    "        \n",
    "        error = np.random.normal(loc=0.0,scale = error_stdev,size = None)\n",
    "        y = intercept + slope*x + slope2*x2 + interaction_slope * x_interaction +  error\n",
    "        y_predictions.append(y)\n",
    "        \n",
    "    return y_predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we plot the data for an example simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAESCAYAAADTx4MfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XuYVNWZ7/Hvy6UF71EwMTQEUBCRACp6TNARdRQ1BsYc4xhMgomGRIxjZHIUYuYkkZjHzPioueqQUdGTeAveOJoYL2eMM0ZFUJRbUAKIrR0FNVy0oaH7PX/sXU119a5rV9Xe1fX7PE8/dO29q2qxoevttd613mXujoiISLF6xd0AERGpTQogIiJSEgUQEREpiQKIiIiURAFERERKogAiIiIliS2AmNmtZvaOmS3POH6pma02sxVm9q9px+eY2Zrw3OTqt1hERNL1ifG95wM/B+5IHTCzk4CpwFh332FmB4XHRwPnAUcAHweeMLOR7t5W9VaLiAgQYw/E3Z8G3ss4fDFwrbvvCK95Jzw+Fbjb3Xe4+zpgDXBs1RorIiJdxNkDiTISOMHMrgG2A9929xeAQcBzadc1hce6MLMZwAyAvfba6+hRo0ZVtsUiIj3MkiVLNrn7wHzXJS2A9AE+AhwHHAPca2bDAYu4NrIGi7vPA+YBTJgwwRcvXlyhpoqI9Exm9noh1yVtFlYTcL8HFgHtwIDw+OC06xqBt2Jon4iIhJIWQB4ETgYws5FAA7AJWAicZ2Z7mNkwYASwKLZWiohIfENYZnYXMAkYYGZNwPeAW4Fbw6m9rcB0D8oFrzCze4GVwC7gEs3AEhGJl/Xkcu5ROZCdO3fS1NTE9u3bY2pVz9OvXz8aGxvp27dv3E0RkTIwsyXuPiHfdUlLoldcU1MT++yzD0OHDsUsKjcvxXB33n33XZqamhg2bFjczRGRKkpaDqTitm/fzoEHHqjgUSZmxoEHHqgenUiStDTDwkOg5a8VfZu6CyCAgkeZ6X6KJEQqcCydA9vWw/K5FX27ugwgIiI90rK5sG0drP810A5rb6toL0QBJAbXXHMNRxxxBGPHjmX8+PE8//zzAAwdOpRNmzZ1uf7Tn/500e9x3XXXMWrUKMaMGcO4ceO44447cl5/wQUXsGDBAgAmTZqEFmCK1JiWZlh3G+CQmqTqbRXthdRdEj1uzz77LA8//DAvvvgie+yxB5s2baK1tTXnc/70pz8V9R4333wzjz/+OIsWLWLfffdl8+bNPPjgg91ptogk3bK50J6xuqG9NeiFjPkX6P+xsr+leiCFKGNCqrm5mQEDBrDHHnsAMGDAAD7+8Y93fruWFk4//XR+9atfAbD33nsD8NRTTzFp0iTOOeccRo0axfnnn0/UNOwf/ehH/PKXv2TfffcFYL/99mP69OkALFmyhBNPPJGjjz6ayZMn09zc3O2/k4jELNX78J1dz1WwF6IAUohlc8uWkDrttNN44403GDlyJDNnzuSPf/xjp/Pbtm3js5/9LNOmTeNrX/tal+e/9NJL3HjjjaxcuZK1a9fyzDPPdDq/detWtm7dyiGHHNLluTt37uTSSy9lwYIFLFmyhK9+9atcddVV3f47iUjMls0Fb48+194KTQ9V5G0VQPLpGFcsT0Jq7733ZsmSJcybN4+BAwfyj//4j8yfP7/j/NSpU/nKV77Cl7/85cjnH3vssTQ2NtKrVy/Gjx/P+vXrO51396yzolavXs3y5cs59dRTGT9+PD/84Q9pamrq1t9HRBLgzYVBoMjUfxBMczi7Mj/nyoHkkx7ZU13BY37RrZfs3bs3kyZNYtKkSXzyk5/k9ttv54ILLgBg4sSJ/P73v2fatGmRgSA19JV6nV27dnU6v++++7LXXnuxdu1ahg8f3umcu3PEEUfw7LPPdqv9IlJlLc3w+PFw6jPRuYwKBYh81APJJdX7SEX2VEKqG72Q1atX89prr3U8Xrp0KZ/4xCc6Hl999dUceOCBzJw5s+T3mDNnDpdccglbtmwBYMuWLcybN4/DDjuMjRs3dgSQnTt3smLFipLfR0SqpIzD6OWkAJJL1LhiNxNS27ZtY/r06YwePZqxY8eycuVKvv/973e65sYbb2T79u1cccUVJb3HxRdfzEknncQxxxzDmDFjOPHEE9lzzz1paGhgwYIFXHnllYwbN47x48cXPcNLRKosfRj9L7fCg0MrvsK8UHVXTHHVqlUcfvjhhb3AA43Q8mbX4/0HxdZlTKqi7quI5Jcathp4Arx+VzgS0gtohxEzuz2UnouKKZaDgoSIxCW1qvyD13cvDCQcEfnLrRVb21EMDWGJiCRN1KrydO2ticiHKICIiMQh1wLlXOs6gGrUuSqEAoiISByyzazKnP0JQG/o1dD5ugrXuSpEbAHEzG41s3fC7Wszz33bzNzMBoSPzcx+amZrzOwVMzuq+i0WESlAIaWPci1Qjux9tHVdKFjBFeaFirMHMh84PfOgmQ0GTgU2pB0+AxgRfs0AbqpC+0REiperZ5EKLFELlFPyrSpP/4p5ok9sAcTdnwbeizh1A3AFkD6/eCpwhweeA/Y3s4Or0MyKSBVHTJk/fz7f/OY3y/La69ev584778x6/tVXX+XMM8/k0EMP5fDDD+fcc8/l7bffzvl6Y8aMAYJijmeddVZZ2inSI+XrWWxbD0tn516gfPoLsPdwOLs5UcEiSqJyIGY2BXjT3V/OODUIeCPtcVN4LOo1ZpjZYjNbvHHjxgq1NJl27dqVM4Bs376dz3zmM1x88cWsWbOGVatWcfHFF1Nv90mkYrL1LNIDy7pf516gnNBV51ESsw7EzPYErgJOizodcSxyBaS7zwPmQbCQsDttuv7lTbRGTIRo6AWzxg3ozktntXHjRr7xjW+wYUMwgnfjjTcyceJEFi1axLe+9S1aWlro378/t912G4cddhjz58/nkUceYfv27XzwwQd8+OGHrFq1ivHjxzN9+nQuv/zyjte+8847+dSnPsVnP/vZjmMnnXQSAG1tbcyePZunnnqKHTt2cMkll/D1r3+9In9HkR4pW+mjMf+Skddoi963o+khGPPdzj2YBKz1yCUxAQQ4BBgGvBwWEWwEXjSzYwl6HIPTrm0E3qp0g6KCR67jhWppaWH8+PEdj9977z2mTJkCwGWXXcbll1/O8ccfz4YNG5g8eTKrVq1i1KhRPP300/Tp04cnnniC73znO9x3331AsEnVK6+8wgEHHMBTTz3Fddddx8MPP9zlfZcvX87RRx8d2aZbbrmF/fbbjxdeeIEdO3YwceJETjvtNO13LpIpW2HDbKWPls6GDfd0zmv07g9T1nYNDotmlr14ayUlJoC4+zLgoNRjM1sPTHD3TWa2EPimmd0N/A9gs7vX7E5I/fv3Z+nSpR2P58+f37GF7BNPPMHKlSs7zm3ZsoWtW7eyefNmpk+fzmuvvYaZsXPn7o1jTj31VA444IButemxxx7jlVde6djWdvPmzbz22muMHDmyW68r0uOkDzGlf7hHJb/bW2HDguxDVunPz9WDSWgvJLYAYmZ3AZOAAWbWBHzP3W/JcvnvgDOBNcCHwFeq0sgYtLe38+yzz9K/f/9Oxy+99FJOOukkHnjgAdavX8+kSZM6zu21114FvfYRRxzRZQOrFHfnZz/7GZMnT+50PHO/EZG6lpkkT/9wz5bkjqqplxqySg8guYq3JrQXEucsrC+4+8Hu3tfdGzODh7sPdfdN4ffu7pe4+yHu/kl3Xxz9qrXvtNNO4+c//3nH41RPZfPmzQwaFMwbSN+AKtM+++zD1q1bI89NmzaNP/3pTzzyyCMdxx599FGWLVvG5MmTuemmmzp6Nq+++ioffPBBd/86Ij1Lrum32Zzd1HX6bdSsqmw9mJjXeuSSqFlYAj/96U9ZvHgxY8eOZfTo0dx8880AXHHFFcyZM4eJEyfS1hZRGyc0duxY+vTpw7hx47jhhhs6nevfvz8PP/wwP/vZzxgxYgSjR49m/vz5HHTQQVx00UWMHj2ao446ijFjxvD1r3+9y2ZVInWtAvsDdVJooEkQlXPPIY5ZWLVK5dylx1s0E9be0rmX0KsBDrkosUNMpVI59zJQkBCRDrmGmHpYACmUAoiISCESPJQUl7rMgfTkYbs46H5KTSqk6GEx19Whugsg/fr1491339WHXpm4O++++y79+vWLuykixSm0ZEi+6+o4wNTdEFZjYyNNTU2q/1RG/fr1o7GxMe5miBQu13qO1ErzE+6HP06B7e9EX5eSbWFhHai7ANK3b1+GDRsWdzNEJE5R6zlSH/6pgPDM+fDhBjoGanKtHq+R2lXlVndDWCJS53Kt50gPCFtWhE9o73pdSikLC3sQBRARqS+5Sobk24s8qkR7pRYW1gAFEBGpL1mLHv4W1twcvRtg+nWp0iK5AlGdqLsciIjUucz1HKmkufWFHVkm10StONfCQgUQEalzy+bCtnVk2aMu0N4Kb9wPzY/u3gdECws1hCUidawjaZ4WPHo1wIiZXYsaNp5dM1vNVosCiIjUl/SFf8vmRm8vm5kMz5yuW0eJ8lwUQESkvqTWeSydHQQF39n1msxkeJ1P180mtgBiZrea2Ttmtjzt2L+Z2Z/N7BUze8DM9k87N8fM1pjZajObHP2qIiI5pPck1v06+5Td9NlWmq6bVZw9kPnA6RnHHgfGuPtY4FVgDoCZjQbOA44In/NLM+tdvaaKSKIVUo+qpRn+72FpQaMtespu/0GdN3LSdN2s4tzS9mngvYxjj7l7ahu854BUgaWpwN3uvsPd1xHsjX5s1RorIslWSGHEl2bDrq2dg0bv/nB2c4/barZakjyN96vAPeH3gwgCSkpTeKwLM5sBzAAYMmRIJdsnIklQSD2qlmZY/5uuz42qb5UpwdN14941NZFJdDO7CtgFpP7FLeKyyEnb7j7P3Se4+4SBAwdWqokikhSFJLiXzQXauh6v8Z5EVPDIdbzcEtcDMbPpwFnAKb57044mYHDaZY3AW9Vum4gkTLYEd2Z59nW3dX5e7/4wZW1dVc6thET1QMzsdOBKYIq7f5h2aiFwnpntYWbDgBHAojjaKCIxyJYkLyTB3QOT4Ne/vIlrX9oUdzNincZ7F/AscJiZNZnZhcDPgX2Ax81sqZndDODuK4B7gZXAo8Al7h7RHxWRHikzSZ4KKE0P5E9w98AkeLWGqPKJbQjL3b8QcfiWHNdfA1xTuRaJSCJFJclTAWXEN/IXLkxoEjzuBHg5JGoIS0Ski8wkeWoFeY2XFalkAryhSp/sCiAiUn2FLPxLXZeZJE9fQV7juYxKmH3kgKr1YBRARKT6Cln4l7quS7mRtrotK5KU5HmKAoiIVFcxlW2jEuCZ6qgXkm94q1pDVykKICJSXbkW/mUObZ3d1LnMSP+IAhQ1PqOqXKo5dJWSuIWEItKD5Vv4lz60FTW7KqEzqkrR0Cu6R1HtXkR3KICISPXkWtQ35rv5a1r1ILUyVTcXBRARqZ5ci/rcuw5t5VvjIbFSABGR8mtphsePh1Of6dyLyDYE1dIMC4fnrmkliRv2UgARkfJL3zZ24391DSRR12cb2qrRXki2leaZill5nrRhrxpK14hITcjcNnbbuvzTbOu4XlVS6lqVQj0QESmvTr2JsOZpvuGoHjS7qp6oByIi5ZM5TTeljhb71RMFEBEpn8jSI9RdyZF6oSEsESmfXKVHCkiK94QS5/VEPRARKZ9U6ZESS47Evcd3ORU6tbaWVp5niq0HYma3Eux9/o67jwmPHQDcAwwF1gPnuvv7ZmbAT4AzgQ+BC9z9xTjaLSIFUFK8LnpMcca++cDpGcdmA0+6+wjgyfAxwBkE+6CPAGYAN1WpjSKSrtB9PKQuxBZA3P1p4L2Mw1OB28Pvbwf+Ie34HR54DtjfzA6uTktFpEOh+3hIXUhaEv2j7t4M4O7NZnZQeHwQ8EbadU3hsebMFzCzGQS9FIYMGVLZ1orUk6i9yVVmpIt6mghQK+kbizjmURe6+zx3n+DuEwYOHFjhZonUkVz7eJRJtoRyLSWae9JEgHyS1gN528wODnsfBwPvhMebgMFp1zUCb1W9dSL16v2lsOZmOn5vK7bYYbbiihl62m/oPV3S4vpCYHr4/XTgobTjX7bAccDm1FCXiFTBM1+kS6e/rQWWzins+cqd9EhxTuO9C5gEDDCzJuB7wLXAvWZ2IbAB+Hx4+e8IpvCuIZjG+5WqN1ikXrU0w5YV0ec2/BY+dVv+55c5d9KT8gy1/HeJLYC4+xeynDol4loHLqlsi0Skk9Sw08AToFdDMGzVqwE+8QXYcA+0bQfagym9hZZqL1OJ9lx5hmtf2tTxuBY+hGs5Z5K0HIiIJMWyuUEp9g9eDz74IQgi634NvXoHj/MFhHx7oBeo0L01MsXxIZy0TZ8qSQFEpN4UktDuGHby3cGjQxu0pwWUXAEhz0ZRhQ7fJPm38VoeguquHhgTRSSnQhLa2arqRsk1pTfPRlG1PHyT0hP+DqVSD0SknhSS0I7a06N3f5iyFh6dAC1vdr4+FRCihrFUE6tHUwARqSe5EtrpSfNsw041GBDKnVQvNR+TTS3nTBRAROpFvoR2amir5e3sw07dnD2VUs4P4dlHDij49crxnuUemqrlPIkCiEi9yJXQHvPd3UNbtMPZzRWtc1Xsh3C+39JTH8LpvQ2pPAUQkaQpsOxH0TIS2tePXEtr732CB38GRgb1ShvatjKrDGs1uiNz+KYWf0uvhSGo7lIAEUma9FlS5fwQz8hftGb5bb219z7wanzVdmcfWXvBIqWW216KOoiRIjUkc5ZUXBs3VajarvQs6oGIJEkFyn6U4trDwqm6L22qqQVx2XIl6ecr9R71MGSVSQFEJCnKVPaj3KI+LAud9ZQKPsVc3x3VCHS1EkyrQQFEJCnylP1IkkJnUaWuy3V9veUNehIFEJGkyFX2o4QAcv3LG2lt77qZZ6pXkG+4RyQfBRCRpCjzKu+o4BEcD/7MHIpJ0hqKei5QWEvqMO0jUgdacm/Yee1Lm7j+5eQEjEz1XKCwliSyB2JmlwMXEeyhuYxgB8KDgbuBA4AXgS+5e2vWFxHpYYr6rXzZXOh7dc7Xy3wtzS6SYiUugJjZIOCfgNHu3mJm9wLnEWxpe4O7321mNwMXAjfF2FSRqir4t/LUbK6RuQNIpmKGhgrNn6SCj4JTz5S4ABLqA/Q3s53AnkAzcDIwLTx/O/B9FEAkaSpVhqQYxezlUaJCak+lz65S3qJnSlwAcfc3zew6YAPQAjwGLAH+5u67wsuagEExNVEku7AMyfWr2mi1rh+ulUoCdypZ3vtKZrXfREPb1t21rkQqIHEBxMw+AkwFhgF/A34LnBFxqWd5/gxgBsCQIUMq1EqRCGllSFptj8hLWtsrP8OotddeXHv4xm6/Tpw05FUbEhdAgL8H1rn7RgAzux/4NLC/mfUJeyGNwFtRT3b3ecA8gAkTJkQGGZGKKHDoqOgZRqlhscHPl962CEku66Ehr9qQxACyATjOzPYkGMI6BVgM/CdwDsFMrOnAQ7G1UCRT1DawJbg2qvZUOCzW4Duy9mzySc9HpHpAre3d361PH/T1LXEBxN2fN7MFBFN1dwEvEfQoHgHuNrMfhsduia+VUq+yDj95G7PKlLju9Pppw2KzXhsBU9Zy7Z+792OrNRZSLkV1NM3sGjM7tFKNSXH377n7KHcf4+5fcvcd7r7W3Y9190Pd/fPuvqPS7RDJlPXD1/p1u/cRKao6r0hCFPurzGxgtpk9B9wG3OvuW8rfLJHaU0jiuqGXZy0x0kW26rwjf9CNVoqUT7EBZDhB/uGLBMNKPzGzB4D57v5EuRsnUssiq8wumsn1vf5X/um1LX/NXp23SJq5JJVS1H8td3/d3a9295HACcCvCVaI/8HMNpjZD81sRCUaKlLzwh7FrFeHM3vVwNzXLp+btTpvQ/sHBb1dQ68giCnRLZVScjbO3Z8BnjGzSwnWbVwEzAHmhENctwC/Ua5CJLR0NrRt73iYdaGfO9f2+QEM7TxUlZolNSvj8myrwbPla7TGQsqlHLOwPknQGzkaMOBVYD/gP4AfmNm57v5sGd5HJHbd2kNjw32dHs56dXjwTe+9wHd29Day5VLKMe0WNPVWyqekAGJmBxPkQaYDhwMfAPcCt4U9E8zseIIgcjMwriytFYlZyXtotDQDYf6id3+YsnZ3rawHGqHlzaLbomm3Erdip/GeZ2a/J1js92PgXeCrwMfc/aJU8ABw9/8G/hUYXcb2iiRKtmGfhl4OCw8JkuGQezru2U0wzXd/idSIYnsgdxIUMvwxQW/jL3muX06QaBfplmrvUJft/aLeP9tsK7atDwLFmO9GT8cd8y/xVewVKYNiA8gZwGPuXtCvSe6+CFhUdKtEMlR79XShrxt5XdrqcdbeBrs+iJ6Ou3xuSXudZ1JSXOJSVABx9z9UqiEipYrKQ8S6d3bmcNWGBZHTcWl6KDKAFJuoV1Jc4pK4Wlgi5RBbgjlq9Xjv3nB2c8HDVZkBIdfwnUicFECkxyrHlNeiZVs93o3hKvUwJKn0O4zUhdQaiutfLnDabamyrB6nSbsPSM+jHojUhG4t4EuTeo18s7oKfb8uw0hnN3W7jSK1QgFEakK2YZyCF/JlyDWrq9DXjJy+K1JHNIQlNU2JZJH4qAciNa3k0iLdpMAlktAAYmb7E9TRGgM4QbmU1cA9wFBgPXCuu78fUxMlocqVK4miISuRzpL6e9RPgEfdfRRBIcZVBLshPunuI4Anw8cincwaN4DZRw7IXqOqbdvu+lQi0i2J64GY2b7A3wEXALh7K9BqZlOBSeFltwNPAVdWv4VSCzoNbS2aCWtvCabT9moAv4iGhh+omq1INyUugBBsm7sRuM3MxgFLgMuAj7p7M4C7N5vZQVFPNrMZwAyAIUOGVKfFklxZ9hWfddoM+K/PwcAT4PW7Os5fP3Jt5CZPynmIdGUF1kWsGjObADwHTHT3583sJ8AW4FJ33z/tuvfd/SO5XmvChAm+ePHiyjY4ZtWuUltVLc3w+PFw6jOlV61N732k9GqAvUfAllXB3hxtEVvE9h+kNR1St8xsibtPyHddEnsgTUCTuz8fPl5AkO9428wODnsfBwPvxNbCBCmmSm1cwabk9102d3dJ9FKr1mZbGb5lRepBUXWqRGS3xAUQd/+rmb1hZoe5+2rgFGBl+DUduDb8U7Uh8sisBVXtkuj5Xj/n+2aWRC9174ywF5E1iLVtZVaZyqqL1JvEBZDQpcBvzKwBWAt8hWDG2L1mdiHBjoifj7F9Nac7QaLQzZVSUj2LYp/XSdQOft34kM8axHrvA69qcyeRUiQygLj7UiBq/O2UardFig8+qetLDh5Rie+/3ApvPgKTnyv/B30ZN3cSqSeaWyLJE1USvb0VPny9817i5aJquSIlUQCpcUmdXtqtkiJRiW/CgLL2tvIvBJzmmnElUoJEDmHVs2JnLJWjFlRcQSjr+57+Qufpu+lTcTXcJJIYCiAJU+iMpe4kqJNQ0ylnG9Kn7475buRCwGKT3tlmoSW1BydSCxRAalSu4DH7yAFlrUpb7gKFOT+0M6fv7vqgLFvEFrzOpRyLF0XqhAKI5JXvw7eQYJXe47j+5U2Rz2noBbN2ZEzf3bAgeiHgG/dD86NdP+i7GwDKsXhRpE6oAy9Vl3OYLnO4KrVSfJp3/mo8e/cHfbr0AFCszN6PqvaK5KQeSB0rV2mTfENcReUZChmuamnm+l5X0Hr41cHjTivur2AWN5W2er3MixdFejoFkDLr7odyuZK9hbxOuUqblLWOVtRwVdNDnT/Il82lte/VkU9v7b138E2xASBL1V6tUBfJTgGkzIr9UC50NlVre5BrSAWifAGiZivx7j08d/4i9UE/MjqAdCg2AEQtXlQvRCQnBZCYlVompGYDRD75EthRH/TZFBMAslXtzez9iEgHBZAqKufU2lrW0Mtpbbeux9u2kbf6buQq9SyKCQBaiS5SNAUQqbpZO74Ha/4dRnxj94d7arU55O45pD7ocwXjacnaJE2kp9I03jqWLTFf0dXZUVNlsyWwc0yjjaXtItKJeiBlVu5V25UUSx4laqqse9EJ7B6bAxKpIYnbE72ckrQneqXzHzWxB3pLMywcDm3bdx/r3R/67gfbI3ob2pdcJBa1vCc6AGbWG1gMvOnuZ5nZMOBu4ADgReBL7l5gNjXZiilumC0Q1USvJ9tU2cGf00wnkRqU5BHjy4BVaY9/DNzg7iOA94ELY2lViTRmT+6psiJScxLZAzGzRuAzwDXALDMz4GRgWnjJ7cD3gZsq3ZZylftI/PBSNWg4SqRHSWQAAW4ErgD2CR8fCPzN3XeFj5uAQVFPNLMZwAyAIUOGdLsh5Sr3kVKugCQiErfEDaCY2VnAO+6+JP1wxKWR2X93n+fuE9x9wsCBAyvSxu4od0BKtJZmWHiIqtqK9FCJCyDARGCKma0nSJqfTNAj2d/MUj2mRuCteJoXr5rKpRRaWl2BRqQmJW4Iy93nAHMAzGwS8G13P9/MfgucQxBUpgN1mXnt1jBXNXfby1wwmKuooTZxEqlJSfy9NZsrCRLqawhyIrfE3J7a053Nlkp5r8wFg1G0iZNIzUp0AHH3p9z9rPD7te5+rLsf6u6fd/cd1WhDTQ0Z5VLMB3V3h5SKKU1SaKARkcRJ3BBW0pR7ZlS5NowqWjG77XV3SKnQvTW0iZNITVMpk3qQrYTIlLVdP6jTr01dgxeXO3mgEVre7Ho8szRJqgJv+uLCXg1wyEXKhYjEqOZLmUgZFbPbXrZih8X0SApdMKhNnERqmnog9aDQHsH7S+H3R9FpiU2vfsGf7du79lqqOatLRKpGPRDZrdAewTNfpMv6zPQeQmavRdNvReparc0lkkppaYYtKyJOtIdf7E5yv/8yPPiJ4HtNvxWpW+qBJETsNbKWzQ0S2O2tuxPZ7l2T3N4Gz5wPH26g4/ePfLO6RKRHUgDppnJ98MdaIyvbdNo++0QnuTt6Khk9E02/FakrGsLqph5RHDHbLK09BgK9YMRMmObB16EXE/nfRosAReqOAkhPU8oq8mzTabesoFM1UCOWAAAKfElEQVSOI301eyZtDCVSdzSE1dOUMjMqapZW+iK/9PUgmT0VLfwTqVvqgfQkxRYmzNZbyZYTaXpAW9KKSAf1QBKiLDWyiql3lbo+qreSLScy+HPqaYhIBwWQbipXccS8M7byrfoutjBhrv06VGJERAqgANJNVdvHPF9uo5h6V5nXZ15X6Mp1EalryoHUgkJyG7l6Ddler5D9OkREskhcD8TMBgN3AB8jmC86z91/YmYHAPcAQ4H1wLnu/n5c7ayqQnIbxfQaiu2tiIhESGIPZBfwz+5+OHAccImZjQZmA0+6+wjgyfBxz1eJ3kIxvRURkSwS1wNx92agOfx+q5mtAgYBU4FJ4WW3A08R7JPes1Wit6Ach4iUQRJ7IB3MbChwJPA88NEwuKSCzEFZnjPDzBab2eKNGzdWq6mVo96CiCRU4nogKWa2N3Af8C1332JmBT3P3ecB8yDYUKpyLawS9RZEJKES2QMxs74EweM37n5/ePhtMzs4PH8w8E5c7RMRkQQGEAu6GrcAq9z9+rRTC4Hp4ffTAY3hiIjEKIlDWBOBLwHLzGxpeOw7wLXAvWZ2IbAB+HxM7RMRERIYQNz9v4FsCY9TqtmWkmSWHMlXgkREpEYlbgir5qWXHIl6LCLSQyiAlFNmyZH3Xy6uvLqISA1RACmnzJIjz5zftQSJiEgPoQBSLlElR7asKH/BwlK2rBURqQAFkHKJKjmSKaoXUmxAUE5FRBJCAaRcokqOZIoqQVJMQCh2y1oRkQpK3DTempWv5Ej6dN70Y9l2BYxS7Ja1IiIVpB5ItUT1NKICQjbaBEpEEkYBpBqihp6KDQi5yrqLiMRAAaQaonoaxQYElXUXkYRRDqTSsvU0+u6XPSBE5TVU1l1EEkYBpNKy9TQGf04JcBGpaRrCKlbmuo186zg09CQiPZR6IMVKn011zC+6Ps6koScR6aHUA4mSrVehYokiIh0UQKJkWx2ens9o3wWPHqNiiSJStxRAMmUrF5I5m8p3Bl9a2CcidarmAoiZnW5mq81sjZnNLvsbZFsdXmqxRBGRHqqmAoiZ9QZ+AZwBjAa+YGajy/YGuVaHl1osUUSkh6q1WVjHAmvcfS2Amd0NTAVWluXVc60OT82mammGhcOhbfvua3r3hylrtee5iNSVmuqBAIOAN9IeN4XHOpjZDDNbbGaLN27cWNyrF7JmQzWpRESA2uuBWMQx7/TAfR4wD2DChAkecX12hazZyBVktLJcROpIrQWQJmBw2uNG4K2KvVv6Hh6p4SktDBQRAWpvCOsFYISZDTOzBuA8YGHF3k3bx4qIZFVTAcTddwHfBP4ArALudfcVFXkzbR8rIpJTTQUQAHf/nbuPdPdD3P2air1RMbsFiojUoZoLIFWh7WNFRPJSAImiqboiInkpgETRHh4iInnV2jTe6tBUXRGRvNQDERGRkiiAiIhISRRARESkJAogIiJSEgUQEREpibkXV7C2lpjZRuD1Ep8+ANhUxuZUgtpYHmpjeaiN3ZeU9n3C3Qfmu6hHB5DuMLPF7j4h7nbkojaWh9pYHmpj9yW9fZk0hCUiIiVRABERkZIogGQ3L+4GFEBtLA+1sTzUxu5Levs6UQ5ERERKoh6IiIiURAFERERKogASwcxON7PVZrbGzGbH3R4AMxtsZv9pZqvMbIWZXRYeP8DMHjez18I/PxJzO3ub2Utm9nD4eJiZPR+2755wL/s427e/mS0wsz+H9/JTCbyHl4f/xsvN7C4z6xf3fTSzW83sHTNbnnYs8r5Z4Kfhz88rZnZUjG38t/Df+hUze8DM9k87Nyds42ozmxxXG9POfdvM3MwGhI9juY/FUADJYGa9gV8AZwCjgS+Y2eh4WwXALuCf3f1w4DjgkrBds4En3X0E8GT4OE6XEexXn/Jj4Iawfe8DF8bSqt1+Ajzq7qOAcQRtTcw9NLNBwD8BE9x9DNAbOI/47+N84PSMY9nu2xnAiPBrBnBTjG18HBjj7mOBV4E5AOHPznnAEeFzfhn+7MfRRsxsMHAqsCHtcFz3sWAKIF0dC6xx97Xu3grcDUyNuU24e7O7vxh+v5Xgg28QQdtuDy+7HfiHeFoIZtYIfAb4j/CxAScDC8JL4m7fvsDfAbcAuHuru/+NBN3DUB+gv5n1AfYEmon5Prr708B7GYez3bepwB0eeA7Y38wOjqON7v6Yu+8KHz4HNKa18W533+Hu64A1BD/7VW9j6AbgCiB9VlMs97EYCiBdDQLeSHvcFB5LDDMbChwJPA981N2bIQgywEHxtYwbCX4IUvsBHwj8Le0HOO57ORzYCNwWDrP9h5ntRYLuobu/CVxH8JtoM7AZWEKy7mNKtvuW1J+hrwK/D79PTBvNbArwpru/nHEqMW3MRgGkK4s4lpi5zma2N3Af8C133xJ3e1LM7CzgHXdfkn444tI472Uf4CjgJnc/EviA+If8OgnzCFOBYcDHgb0IhjIyJeb/ZISk/btjZlcRDAP/JnUo4rKqt9HM9gSuAv531OmIY4n6d1cA6aoJGJz2uBF4K6a2dGJmfQmCx2/c/f7w8Nupbm345zsxNW8iMMXM1hMM+51M0CPZPxyKgfjvZRPQ5O7Ph48XEASUpNxDgL8H1rn7RnffCdwPfJpk3ceUbPctUT9DZjYdOAs433cvfEtKGw8h+GXh5fBnpxF40cw+RnLamJUCSFcvACPCWS8NBIm2hTG3KZVPuAVY5e7Xp51aCEwPv58OPFTttgG4+xx3b3T3oQT37P+5+/nAfwLnxN0+AHf/K/CGmR0WHjoFWElC7mFoA3Ccme0Z/pun2piY+5gm231bCHw5nEV0HLA5NdRVbWZ2OnAlMMXdP0w7tRA4z8z2MLNhBInqRdVun7svc/eD3H1o+LPTBBwV/l9NzH3Myt31lfEFnEkwY+MvwFVxtyds0/EE3ddXgKXh15kEeYYngdfCPw9IQFsnAQ+H3w8n+MFcA/wW2CPmto0HFof38UHgI0m7h8APgD8Dy4H/A+wR930E7iLIyewk+JC7MNt9Ixh6+UX487OMYEZZXG1cQ5BHSP3M3Jx2/VVhG1cDZ8TVxozz64EBcd7HYr5UykREREqiISwRESmJAoiIiJREAUREREqiACIiIiVRABERkZIogIiISEkUQEREpCQKICIiUhIFEBERKYkCiEgVmFlDWEL+r2Y2MOPcI2bWkpCNy0QKpgAiUgUebE52PrAf8KvUcTObSVDT7Ap3XxlT80RKolpYIlVkZpcCPwW+BvwX8GL45xmuH0apMQogIlUUlmj/HXACsBY4GBjrSSvTLVIABRCRKgs3C/oLwX7n57j7fTE3SaQkyoGIVN+JBMEDYFycDRHpDvVARKrIzBoJNrNaSTCENQ043t2fi7VhIiVQABGpkjD/8QRwDEHP413gZaANGO/u22JsnkjRNIQlUj2zgJOBy9x9nbtvAb4EDANuiLVlIiVQD0SkCszsk8ALwCPu/j8zzv0ImANMdfeFcbRPpBQKICIiUhINYYmISEkUQEREpCQKICIiUhIFEBERKYkCiIiIlEQBRERESqIAIiIiJVEAERGRkiiAiIhISf4/CRmAxG3BxaIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from random import choice\n",
    "\n",
    "#simulate observed data points\n",
    "N = 150 #number of data points\n",
    "\n",
    "\n",
    "#x values will be a number\n",
    "x = [i for i in range(N)]\n",
    "\n",
    "#x2 values in this case will be a category\n",
    "#represented by 0 or 1\n",
    "x2 = [choice([0,1]) for i in range(N)]\n",
    "\n",
    "x1_slope = 1  #increase in y with each x\n",
    "x2_slope = 35 #difference between cyan and orange\n",
    "intercept = 0 #intercept\n",
    "error_stdev = 5.0 #measure of 'scatter' or 'spread' of points around the model prediction\n",
    "interaction_slope = -0.5\n",
    "\n",
    "observed_y = linear_simulation_with_interaction(x,x2,slope=x1_slope,slope2=x2_slope,\\\n",
    "  interaction_slope = interaction_slope,intercept=intercept,error_stdev=error_stdev)\n",
    "\n",
    "true_y = linear_simulation(x,x2,slope=x1_slope,slope2=x2_slope,error_stdev=0)\n",
    "\n",
    "plot_colored_scatterplot(x=x,y=observed_y,x2=x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how now the orange and blue series now not only have different intercepts, but also different slopes. Specifically, the higher that x1 is, the bigger the product of x2 * x1 * the interaction slope (-0.5) becomes. When x1 is low, this product is small, and so the large x2 slope causes the blue points to be higher. As we move along the x axis, the product of x1*x2*interaction_slope becomes larger for the blue points (but not orange ones since orange is assigned an x2 value of 0 and therefore this product is always 0 for the orange points). This is an interaction effect!\n",
    "\n",
    "**Stop and think**. What do you think would happen if you set the parameter interaction_slope to 0 in the above code? What about +0.5? Jot down your predictions, then try it out for yourself. Can you explain why the graph looks the way it does?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.982\n",
      "Model:                            OLS   Adj. R-squared:                  0.982\n",
      "Method:                 Least Squares   F-statistic:                     2648.\n",
      "Date:                Thu, 15 Jul 2021   Prob (F-statistic):          5.04e-127\n",
      "Time:                        14:54:57   Log-Likelihood:                -444.73\n",
      "No. Observations:                 150   AIC:                             897.5\n",
      "Df Residuals:                     146   BIC:                             909.5\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept      0.9836      1.070      0.919      0.359      -1.131       3.098\n",
      "x              0.9965      0.012     80.405      0.000       0.972       1.021\n",
      "x2            34.3087      1.547     22.177      0.000      31.251      37.366\n",
      "x:x2          -0.5005      0.018    -27.873      0.000      -0.536      -0.465\n",
      "==============================================================================\n",
      "Omnibus:                        0.206   Durbin-Watson:                   2.240\n",
      "Prob(Omnibus):                  0.902   Jarque-Bera (JB):                0.145\n",
      "Skew:                          -0.076   Prob(JB):                        0.930\n",
      "Kurtosis:                       2.982   Cond. No.                         438.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "df = pd.DataFrame({\"x\":x,\"y\":observed_y,\"x2\":x2})\n",
    "df.set_index(\"x\")\n",
    "\n",
    "model = smf.ols(formula='y ~ x + x2 + x * x2', data=df)\n",
    "results = model.fit()\n",
    "summary = results.summary()\n",
    "print(summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Interpretation**. The model is interpreted as above. Pay careful attention to the coef column that indicates each slope, the R<sup>2</sup> value that indicates effect size, and the p-values for each term (which indicate statistical singificance). The main difference from the other examples in this section is that now there is a row in the results labeled `x:x2`. This row represents the slope and significance of the interaction between `x` and `x2`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Is error model error or measurement error? \n",
    "\n",
    "So far we're established that in any real experiment there will be measurement error. That error will cause our *observed* measurements to deviate from expected values, and those differences can be measured with residuals. \n",
    "\n",
    "If measurement error were the only reason data could deviate from our model, then our job would be relatively simple: we could simply measure the average distance our points fell from our line — that is, the average magnitude of our residuals — and from that we'd have a sense of the scale of the measurement error.\n",
    "\n",
    "There is one major problem with applying our discussion above to the real world. The whole time we assumed that the biological rate of DNA replication really was perfectly linear, and it really occurred at the rate we said it did (the rate we put into the `rate_of_dna_replication` variable). Those assumptions were useful to discuss measurement error, but in reality we can rarely or never just assume our model is correct. Instead, we have to test that using data.\n",
    "\n",
    "But we have just established that the presence of measurement error means our data will never perfectly match our model, so how do we tell whether differences between our actual data and our prediction are due to **measurement error** or **model error**?  \n",
    "\n",
    "For instance, the results shown up above may look pretty good with a regression line drawn through them, but what if in truth the `y` value were not changing at all? This might happen, for instance, if replication was not actually happening. If the 'no replication' model were correct, we would have a slope of 0, reflecting the idea that over the course of the experiment no new DNA is being synthesized. Our results, with a 'null model' of no replication drawn on top of them, might look like this:\n",
    "\n",
    "<img src=\"./resources/linear_scatter_high_r2_linear_regression_null.png\" width=\"400\"  description=\"A cartoon of a scatter plot with x and y axes. A series of orange points is plotted. The points have some scatter to them, but roughly form a diagonal line from the lower left to the upper right of the plot. A black regression line runs diagonally through the points. \">\n",
    "\n",
    "To my eye at least, it is suddenly somewhat less clear that these results represent a real linear trend. If someone told us that this uninteresting null model were the truth, do we really have enough data to establish that that is very unlikely?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frequentist testing of null models with p-values\n",
    "It turns out there are several different approaches to the is question in statistics.  **Frequentist** approaches ask how often we would get results as different from our model as we observed, if our null model (e.g. here that the slope &beta;<sub>1</sub> is 0) were in fact true. This frequency is then represented in a `p` value, with lower p values indicating lower changes of getting these data under the null model. Traditionally, p values for a null model that are less than or equal to 0.05 — a 5% or smaller chance of getting errors this big under the null model — are considered statistically significant. It's important to note that rejecting our null model (in which &beta;<sub>1</sub> = 0) doesn't directly prove that some other model is correct, it merely indicates that somethng more interesting than the null model may be going on.\n",
    "\n",
    "Note that in the output from `statsmodels` `OLS` model, we get a p-value for *each* slope. This reflects that we may find that the influence of some variables under a model may not be significant while others are."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Information-theoretic model comparison with the Akaike Information Criterion\n",
    "An alternative way of trying to figure out what model best describes our data is to use **information theory** to compare among several models. To do this, we would first translate several conceptual models for what is going on into mathematical equations. Then, we can compare these various models based on which has the highest likelihood, which is just the chances of getting the data we actually saw if the model were true. Because more complex models tend to look like they fit data better than simple ones — even if they are totally off-base —, model-comparison approaches apply a penalty for each parameter in a model. The [Akaike Information Criterion](https://en.wikipedia.org/wiki/Akaike_information_criterion), often abbreviated AIC is a simple equation that let's us adjust the score of a model based on it's likelihood — which again is just its chances of producing the data we actually saw — based on the number of parameters in the model:\n",
    "\n",
    "$$AIC = 2k - 2ln(\\mathcal{L})$$\n",
    "\n",
    "Where `k` is the number of parameters in the model (so 2 for the equation for a line, since a line has a slope and intercept parameter), and L is the likelihood you get for the best fit of the model to the data (e.g. the chances of getting the data due to measurement error if the data were really produced by a given model). A *lower* AIC indicates a model we should trust more. But notice that each extra parameter (adding 1 to `k`), has to be balanced by a full natural logs worth of likelihood (ln(*L*)) in order to prevent the score of our model from getting worse. In this way, a straight line that explains the data pretty well will be preferred over a more convoluted model, *unless* that more complex model does a much better job of predicting the data and is therefore justified in its complexity.  \n",
    "\n",
    "In the `OLS` data from `statsmodels.py`, we can look at the `AIC` entry for our model. If we then fit *another* model to the data, we can compare the two `AIC` values. The model with the smaller `AIC` value is a better explanation of the data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian model comparison with posterior probabilities\n",
    "Finally, it is also possible to use **Bayesian** statistics to try to figure out a true **posterior probability** that your hypothesis is correct, given your results. To do so however — and this is a big catch — you must specify all the possible mutually-exclusive hypotheses that could explain your data, as well as the **prior probability** that they were true before looking at your data (i.e. based on previous evidence you already had going into your study). In real situations, it is often very difficult to specify all of these, or to find a defensible prior probability for each. Therefore, the calculated posterior probability is an approximation. Nonetheless, Bayesian approaches can be a very powerful approach, and can be adapted to complex situations.\n",
    "\n",
    "In terms of calculation, the Bayesian posterior probability of a hypothesis (let's call it hypothesis 0 or **H<sub>0</sub>**) given some data — **p(H<sub>0</sub>|Data)** — is just the likelihood of your hypothesis based on the data  — p(Data | Hypothesis <sub>0</sub>) — times the prior probability of your hypothesis — p<sub>prior</sub>. That whole quantity is then normalized by dividing the result by the same calculation applied to all the other hypotheses that could explain your data:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$p(H_0|Data) = \\frac{p_{prior}(H_0) * p(Data|H_0)}{\\sum_{i=0}^{n} p_{prior}(H_i) * p(Data|H_i) }$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does all this mean? There are several implications hidden in Bayes' Theorem that are worth spelling out.\n",
    "\n",
    "1. **Equal evidence means equal chances each hypothesis is correct if all your priors are equal**, and the chances of getting your data are equal under all your hypotheses, then the posterior probability of your \"main\" hypothesis **H<sub>0</sub>** is just 1 divided by the number of hypotheses. This means that if you have no prior reason to prefer one hypothesis over another, and your data could be explained by all hypotheses equally, it's basically drawing straws or rolling dice to figure out which hypothesis is correct — they all have equal posterior probabilities.\n",
    "\n",
    "2. **Any hypothesis that couldn't generate the data you saw is ruled out.** If any hypothesis literally could not explain the data you saw — that is, if you could not generate the data you saw using the mathematical model that represents that hypothesis — then p(Data|Hypothesis) is 0.0 for that hypothesis. Since any number times 0 is 0, this means that the whole numerator is 0, and indeed the final probability that the hypothesis is true is 0. Many well-designed experiments attempt to generate clear-cut situations like this in order to rule out some hypotheses, if warrented by the data.\n",
    "\n",
    "3. **Extraordinary claims require extraordinary evidence**. If before doing an analysis, the chances that a hypothesis is true (as represented in the equation by p<sub>prior</sub>) are very low, then you will need a lot of evidence (as represented by the likelihood of that hypothesis given the data) for that hypothesis to become plausible. For instance, even if model<sub>1</sub> is 10 times more likely to generate the observed data than model<sub>2</sub>, if model<sub>1</sub> starts with a 1 in 100 chance of being true (p<sub>prior</sub> = 0.01), then even after seeing that data you would treat model<sub>2</sub> as more likely. This is similar to the idea that even if you have seen tape of an Unidentified Flying Object that would make sense if high-tech aliens had visited the Earth, and that isn't easily explained otherwise, you may not yet be convinced of the presence of actual aliens since your prior probability that aliens are visiting the earth may be very very low."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- A linear model is one that uses the equation for a line to predict some response variable. \n",
    "<br><br>\n",
    "- A linear regression tries to infer the slope and intercept of a linear model based on the data. \n",
    "<br><br>\n",
    "- Regardless of whether you came up with your model *a priori* or inferred it using linear regression, there are likely to be differences between the data and your model. That is, even if your model is basically correct, data you observe generally won't fall perfectly on the line predicted by your linear model due to **measurement error**. Data may also not fall on the line predicted by your model due to **model error** - that is, because your model is wrong.\n",
    "<br><br>\n",
    "- We can simulate measurement error relative to our model, to see what types of data we should expect if the model is correct, for various levels of **measurement error**\n",
    "<br><br>\n",
    "- The **Central Limit Theorem** says that the sum of many independent and identically distributed random variates form a normal distribution - even if the distribution those measurements come from isn't itself normal. Therefore, although we can't typically know all the sources of measurement error comes from, if we assume that the total measurement error is the sum of many smaller independent errors from various sources, then the Central Limit Theorem suggests that a *normal distibution* is a reasonable way to approximate measurement error.\n",
    "<br><br>\n",
    "- We can build a linear model that incorporates error by adding an error term to the equation for a line: `y = mx + b + error`.If we want this equation to look fancy and statistical, we can write it in a more statistically conventional way: y = &beta;<sub>0</sub> + &beta;<sub>1</sub>x<sub>1</sub>  +  &epsilon; . In this version, &beta;<sub>0</sub> is the constant or intercept term, &beta;<sub>1</sub>x<sub>1</sub> is the slope, and &epsilon; is the error. \n",
    "<br><br>\n",
    "- Because error is random, the error term &epsilon; in our stochastic linear model for each point is not the same. Instead, in this notation each &epsilon; is assumed to be drawn at random from a normal distribution.\n",
    "<br><br>\n",
    "- We can use several techniques, including p-values for null models, AIC values for model comparison, or Bayesian posterior probabilities — to get a sense of whether differences between our model predictions and our actual observations are easily explained by **measurement error** or might instead be due to **model error**.\n",
    "<br><br>\n",
    "- In **frequentist statistics**, p-values can be used to describe the probability that the difference between the predictions of any given model and our actual data is due to chance. We can do this either for our main model (to show that deviations between our model and the data could easily be explained by measurement error) or more typically for a **null model** in which the slope of the regression line is 0. If that null model is a very unlikely explanation for the data (typically because p <= 0.05), that lends support to the idea that the actual model is one in which the slope is not 0.\n",
    "<br><br>\n",
    "- In **information-theoretic** model comparison, we try to figure out which of several models best explains the data, while penalizing models for their complexity. This is important because more complex models are not always better - on small datasets especially, more complex models have many knobs to turn and levers to pull that can make it look like they fit the data almost perfectly, even though in reality the results will not generalize at all to other datasets. While simpler models may not fit the data perfectly, they often generalize better to other datasets. Therefore, the Akaike Information Criterion balances the number of parameters against the natural log likelihood of the model (how well it explains the data). This means that more complex models can be justified - but only if they do a much better job of explaining the data than a simpler model.\n",
    "<br><br>\n",
    "- Using **Bayesian statistics**, we list all the hypotheses that could explain the data, and then compare them using a calculation that incorporates both prior information about the chances that a hypothesis is correct (it's 'prior probability'), and the ability of each hypothesis to explain the data at hand. Bayes' theorem allows us to combine this information to — in principle at least — calculate the actual chances that our Hypothesis of interest is true.\n",
    "<br><br>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1**. Sketch what shape of line you expect from the following equation:\n",
    "\n",
    "y = &beta;<sub>0</sub> + &beta;<sub>1</sub>x\n",
    "\n",
    "if &beta;<sub>0</sub> = 3.5 and &beta;<sub>1</sub> = 0\n",
    "\n",
    "You can sketch this in a graphics program, by hand on graph paper or generate a plot using python code.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2** You are studying the transmission rate of an airborne virus in two different hosts: humans and bats. You expect that these hosts differ in transmission rate of the virus. In both cases, you expect that the transmission rate of the virus depends on viral load in an individual. Write the equation for a linear model that predicts transmission rate based on viral load and host, and incorporates an error term. Note that in Jupyter notebooks you can use the text `&beta;` to generate a Beta symbol &beta;, `&epsilon;` to generate an epsilon &epsilon;, and `x<sub>1</sub>` to generate subscripts for variables like x<sub>1</sub>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3**. Examine your above equation. Did you include an interaction between host and viral load? Why or why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4**. Imagine that you put slips of paper with the numbers 1,2,3,4, and 5 into a hat. Then you reached into the hat, drew a number and wrote it down. You then repeated this hundreds of times and summed up your result. What type of statistical distribution would you expect the results to form? Sketch or descibe how you expect the graph to look.\n",
    "\n",
    "*Hint*: you are calculating the sum of a large number of random variates from a uniform distibution\n",
    "*Hint2*: if you're stuck, reread the section on the Central Limit Theorem "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5**. Write python code that simulates the situation in Exercise 4. Do your simulation results match your expectations?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 6**. You have measured data on the armor thickness of a particular type of crustacean at several ages. In each observatin you have also recorded the presence or absence of predators nearby. Each row of your data table represents one individual crustacean that you measured.\n",
    "\n",
    "Here is how your data table looks:\n",
    "\n",
    "\n",
    "|  Armor Thickness (mm) | Age (months)        \t| Predators (1 = present) \t|\n",
    "|--------------------\t|--------------------\t|------------------------\t|\n",
    "| 20.1                  | 0                   \t| 0                   \t    |\n",
    "| 19.5                  | 0                     | 0                         |\n",
    "| 70.3                  | 0                     | 1                         |\n",
    "| 69.0                  | 0                     | 1                         |\n",
    "| 21.9                 \t| 1                  \t| 0                  \t    |\n",
    "| 20.9                  | 1                     | 0                         |\n",
    "| 22.1                 \t| 2                  \t| 0                  \t    |\n",
    "| 73.4                 \t| 3                  \t| 1                  \t    |\n",
    "| 23.5                  | 3                     | 0                         |\n",
    "| 74.7                 \t| 4                  \t| 1                  \t    |\n",
    "| 75.1                 \t| 5                  \t| 1                  \t    |\n",
    "| 25.0                  | 5                     | 0                         |\n",
    "| 25.3                  | 5                     | 0                         |\n",
    "\n",
    "In python, fit a multiple linear regression of Age and the presence of Predators using your data.\n",
    "\n",
    "**Answer the following in your response:**\n",
    "- How much on average does the presence of Predators increase Armor thickness?\n",
    "- How much on average does each month of age increase armor thickness?\n",
    "- Is the effect of predators on armor thickness significant in these data?\n",
    "- Is the effect of age on armor thickness significant in these data?\n",
    "\n",
    "*Hint*: it is probably easiest to build a data table in Excel (by copying the data on this page), save it as a .csv file, then load that into a `pandas` `DataFrame` using the `read_csv` function. From there you should be able to fit an OLS model in `statsmodels` by adapting the code from earlier in the chapter.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 7**. Graph the data from Exercise 6. You can use the code from earlier in the chapter to help.\n",
    "What do you notice in the graph? Does this make sense with the statistical results you saw earlier?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading Responses & Feedback \n",
    "\n",
    "^Make this a hyperlink using Google Forms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References & Further Reading\n",
    "\n",
    "Jonas Kristoffer Lindeløv, \"Common statistical tests are linear models\". This outstanding python notebook inspired the discussion in this chapter. Check it out [here](https://eigenfoo.xyz/tests-as-linear/)!\n",
    "\n",
    "StatTrek, \"Interaction Effects in Statistics\". This page has a nice summary of interaction effects(https://stattrek.com/multiple-regression/interaction.aspx#:~:text=In%20regression%2C%20an%20interaction%20effect,or%20more%20other%20independent%20variables.)\n",
    "\n",
    "Russell A. Poldrack, *Statistical Thinking in the 21st Century*. https://statsthinking21.github.io/statsthinking21-core-site/fitting-models.html \n",
    "\n",
    "The Method of Maximum Likelihood for Simple Linear Regression\n",
    "https://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/06/lecture-06.pdf\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}