venpopov/fitting-ind-vs-group.ipynb

## fitting-ind-vs-group.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The dangers of fitting computational models to individual subjects with few observations\n",
    "\n",
    "I work with computational models. Today I was confused by the behavior of my model, and I finally figured out that in this case *fitting the model to individual subjects was the wrong thing to do*. In fitting computational models, you typically have the choice to either fit the overall data, ignoring by-subject variation, or to fit a separate model to each subject. Many authors recommend the latter, when possible, and for good reasons. However, fitting individual subjects can be problematic, as in this case. \n",
    "\n",
    "I know my model. I know that for this dataset conceptually it should be able to fit a nice logistic function as in the figure below, which the data showed. When I saw the data I was really excited; however, after I fit the model, suddenly it predicted almost a linear function."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image.png](https://i.imgur.com/dXoYbtr.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I was confused. Why is the model behaving like this, when I know it should be able to fit the shape shown above? Turns out, this is a type of overfitting bias when attemtping to fit individual subjects. This case is quite interesting, because overfitting can actually lead to a worse fit of the training data itself. Overfitting problems with generalization are well known, but typically overfitting at least fits well to the data used to train the model. \n",
    "\n",
    "As it turns out, overfitting can lead to worse results when it *overfits selectively*. I'll elaborate.\n",
    "\n",
    "The model in question is a mechanistic memory model with many parameters and complex implementation and is not actually fitting a logistic curve, but for the current purposes I will show a simplified example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that you observe the following logistic growth behavior of *y* as a function of *x*. In my case, y was probability of recall, but it doesn't really matter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAAMFBMVEUAAAAzMzNNTU1oaGh8\nfHyMjIyampqnp6eysrLHx8fQ0NDZ2dnh4eHp6enw8PD///9LiKPpAAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAIaUlEQVR4nO3di1bbOBhFYWFgWiit3/9tJ4RcHJMYydbl6HjvtTowaWr+5a9O\nHCczCiNZF1oPQGUD2DyAzQPYPIDNA9g8gM0D2LxcwPxFEQ1g8wA2D2DzADYPYPMANi8OZjg0\n/X6Y3QawbFEww+Uf1683twEsG8DmJQMPd26L3Q4VKYSFvZ8OfH4Kvt72dAjgZoWwJLzuCB44\ngmUKYVE4/Tn49D3AzQvTHt4pZksAF2nxufPevWfd3P7wj8Vsm4foEi27POBM3c464NlJVvR2\naNLULYHz/pYWfjNqE+erVtMrWFzJ2tZaztQfI7adHbT+UF3108S2Y9yctIovwBV6dKzW8AW4\nZJUehZdnENuORbWeX2MCOGdCsOcATu7B06ma7CmAU/t+JqwpewrgxOq+it0ewIl1I3sK4LT6\nOXRPAZzSJ2xfvgDHN3nrp/EkKQEcV1+qkwCOqVfdEeCIuj14jwG8XN+6I8DL9a47ArxQ9wfv\nMYDv56E7Anw/F90R4DvZHLzHAL7NS3cE+DY33RHgSXYH7zGAv/LUHQH+ylV3BHg0PniP7R3Y\nW3fcO7C77rhrYPuD99hegfehO+4VeC+6466A+/zQ3Nb2A/z1cdd96Y47Au7tA+u5Atg8gM3b\nDXBv/8lJrnYDvMPzq2N7AVafr1g7ARYfr2D7ANaermi7AJYernB7AFaerXj+wLs8d75mDyw7\nWKXcgVXnqpY5sOhYFfMG1pyqatbAkkNVzhlYcabqGQMLjtQgX2C9iZpkCyw3UKNcgdXmaZYp\nsNg4DfME1pqmaZbAUsM0Lm5f3K5TOIqvXag0S/OidsZspdH5Eu/R26mT0CgCrQAetJeX1ZlE\nonTgQXr94J2/vf+9HMBPh0T2q8gYQiUDX3EFT7I0ppAqFfi6EPj1tujtlE5iCLGSgb+SBFaY\nQa41L5NEj2CBEQTzAW4/gWRpV7ImJ1pqJ1nNBxDN5Vp0658vmwkwvo/yAMb3YRbA+D7OARjf\nhQyA8V2qf2B8F+seGN/legfG94f6Bubt/R/rGhjen+sZGN+IOgbGN6Z+gfGNqltgfOPqFRjf\nyDoFxje2PoHxja5LYHzj6xEY34Q6BMY3pf6A8U2qO2B80+oNGN/EOgPGN7W+gPFNridg3t5f\nUUfA8K6pH2B8V9UNML7r6gUY35V1Aozv2vSBP1eFxXd18sA7Xdc5W+rAe12ZPVsAmweweerA\nPAdvTB54xHdT8sDobgtg89SB8d0YwOaJA+O7NYDN0wbGd3MAmwewedLA+G4PYPOUgfHNEMDm\nCQPjmyOAzdMFxjdLAJsnC4xvnmb78fn3nyzb2R7AeZrtxxDC8N/79u1sD+A8zfbjv7fXz08x\nvrz93badzeGbqTs78v3XcDB+nh7HtwtR1licEuBM3duRf38dP4z8crlhtpTsUH55WXxz9X1P\nfrweD98/L+H1fNMt5lBh/WCAczXfk+8vl0fn6wfObzCHyQLRj7eTdypa3fxlUgivH+ffmhyz\n4x3gy3Pw0yGARZu/TPr1cedO89W/j79uD+OsJPjma/4y6e6dJsCXb0s+BwOcr6h9OQX+qigw\nvhlLBb78O8B9BLB5cTvzfNVqcqJV7koWvjkTfLsQ4JzpAeObNYDNkwPGN28Am6cGjG/m1GAA\nzpwaDMCZE4PBN3cAm6cFjG/2ADZPChjf/AFsnhIwvgUC2DwhYHxLBLB5AJunA4xvkQA2TwYY\n3zIBbJ4KML6FAtg8EWB8SwWweRrA+BYLYPMANk8CGN9yAWyeAjC+BQPYPAFgfEsGsHntgfEt\nGsDmNQfGt2wAmwewea2B8S0cwOY1Bsa3dACb1xYY3+IBbF5TYHzLB7B5LYHxrRDA5gFsXkNg\nfGsEsHntgPGtEsDmNQPGt04Am9cKGN9Kxe3o23UKp4uBp21n/f1pZVE7erbS6Bl7w+KU+NZq\nBfBwWQwcYP3SgSerfwOsXw7gp0NpYvhWKxl4mP1K2c7Ke9OGUoFvDt61wPjWKxn4q8vXpO2s\nujNtas3LpK2vg/GtGMDmpV3JmpxoTb9Gbyf1rrS5BteiAa5ZfWB8qwaweQCbVx0Y37oBbF5t\nYHwrB7B5lYHxrR3A5tUFxrd6AJtXFRjf+gFsHsDm1QTGt0EAm1cRGN8WAWxePWB8mwSwedWA\n8W0TwOYBbF4tYHwbBbB5lYDxbRXA5tUBxrdZAJtXBRjfdgFsXg1gfBsGsHkAm1cBGN+WAWxe\neWB8mwawecWB8W0bwOaVBsa3cQCbVxgY39YBbB7A5pUFxrd5AJtXFBjf9gFsXklgfAUC2Lxi\nwCHgq1Ap4BAQlqgQcAgIawSweQCbx3OweZxFm1fxf4RGLYqDuVkq+Lwa6c1tAIsWBTNbXnYY\nvt8GsGgrgIcR4H5KBx5urJO2Q/XLAfx0CGDRkoGHcfMS71SxVOD5wZy0HapfMvBX15tTtkP1\nW/MyafzuC7Bq64FvfM/vLpToqeC2V9XBQInAl6tW0xOtYXZ5q1hPNX5ISl0N1MFDa1f7s0kA\n562rgQBOr6uBOgCmLQFsHsDmAWwewObJA9e6nJKS1kDLe0gd+M5V0eZp/Y37YQ8BnNygNU/n\nwMfUdqjaPEsBnJoecNfPwZ+J7U+1gZYnAjgtvXOC/p+DtXZnvffBI+seWGlnnpIaqXdgqZ15\nSmqmzoHlHhE/05qn7ytZtDGAzQPYPIDNA9g8gM0D2DyAzQPYPIDNA9g8gM0D2DyAr72Gj3H8\nCC+t58gawNf+hedxfPlUNgrgSb/D+1v41XqKvAE8TfDDBVsDeNpbCG+tZ8gcwNMANm94fuYh\n2rjDSdZ7+N16irwBfO34Muk5/Gs9R9YAvna60PHaeo6sAWwewOYBbB7A5gFsHsDmAWwewOYB\nbN7//qsjrPB0gg4AAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(tidyverse, quietly=TRUE)\n",
    "options(repr.plot.width=4, repr.plot.height=3, warn=-1)\n",
    "library(IRdisplay, quietly=TRUE)\n",
    "\n",
    "# prep the data\n",
    "x = c(1,2,3,4,5,6)\n",
    "y = 1/(2 + exp(-x))\n",
    "df = data.frame(x=x,y=y)\n",
    "\n",
    "# plot the data\n",
    "ggplot(df, aes(x,y)) +\n",
    "    geom_point() +\n",
    "    geom_line() +\n",
    "    theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's generate data for some subjects. Each observation is a 1 or a 0 with a probability determined by the underlying model plotted above. One key issue is that we are dealing with relatively little observations per subject (90 overall, 15 per condition), which, however, is not atypical for memory and cognitive experiments.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>x</th><th scope=col>y</th><th scope=col>nObs</th><th scope=col>subject</th><th scope=col>obs</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>1        </td><td>0.4223188</td><td>15       </td><td>1        </td><td>0.3333333</td></tr>\n",
       "\t<tr><td>2        </td><td>0.4683105</td><td>15       </td><td>1        </td><td>0.3333333</td></tr>\n",
       "\t<tr><td>3        </td><td>0.4878556</td><td>15       </td><td>1        </td><td>0.5333333</td></tr>\n",
       "\t<tr><td>4        </td><td>0.4954626</td><td>15       </td><td>1        </td><td>0.2666667</td></tr>\n",
       "\t<tr><td>5        </td><td>0.4983212</td><td>15       </td><td>1        </td><td>0.5333333</td></tr>\n",
       "\t<tr><td>6        </td><td>0.4993811</td><td>15       </td><td>1        </td><td>0.2666667</td></tr>\n",
       "\t<tr><td>1        </td><td>0.4223188</td><td>15       </td><td>2        </td><td>0.4666667</td></tr>\n",
       "\t<tr><td>2        </td><td>0.4683105</td><td>15       </td><td>2        </td><td>0.3333333</td></tr>\n",
       "\t<tr><td>3        </td><td>0.4878556</td><td>15       </td><td>2        </td><td>0.4000000</td></tr>\n",
       "\t<tr><td>4        </td><td>0.4954626</td><td>15       </td><td>2        </td><td>0.5333333</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllll}\n",
       " x & y & nObs & subject & obs\\\\\n",
       "\\hline\n",
       "\t 1         & 0.4223188 & 15        & 1         & 0.3333333\\\\\n",
       "\t 2         & 0.4683105 & 15        & 1         & 0.3333333\\\\\n",
       "\t 3         & 0.4878556 & 15        & 1         & 0.5333333\\\\\n",
       "\t 4         & 0.4954626 & 15        & 1         & 0.2666667\\\\\n",
       "\t 5         & 0.4983212 & 15        & 1         & 0.5333333\\\\\n",
       "\t 6         & 0.4993811 & 15        & 1         & 0.2666667\\\\\n",
       "\t 1         & 0.4223188 & 15        & 2         & 0.4666667\\\\\n",
       "\t 2         & 0.4683105 & 15        & 2         & 0.3333333\\\\\n",
       "\t 3         & 0.4878556 & 15        & 2         & 0.4000000\\\\\n",
       "\t 4         & 0.4954626 & 15        & 2         & 0.5333333\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "x | y | nObs | subject | obs | \n",
       "|---|---|---|---|---|---|---|---|---|---|\n",
       "| 1         | 0.4223188 | 15        | 1         | 0.3333333 | \n",
       "| 2         | 0.4683105 | 15        | 1         | 0.3333333 | \n",
       "| 3         | 0.4878556 | 15        | 1         | 0.5333333 | \n",
       "| 4         | 0.4954626 | 15        | 1         | 0.2666667 | \n",
       "| 5         | 0.4983212 | 15        | 1         | 0.5333333 | \n",
       "| 6         | 0.4993811 | 15        | 1         | 0.2666667 | \n",
       "| 1         | 0.4223188 | 15        | 2         | 0.4666667 | \n",
       "| 2         | 0.4683105 | 15        | 2         | 0.3333333 | \n",
       "| 3         | 0.4878556 | 15        | 2         | 0.4000000 | \n",
       "| 4         | 0.4954626 | 15        | 2         | 0.5333333 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "   x y         nObs subject obs      \n",
       "1  1 0.4223188 15   1       0.3333333\n",
       "2  2 0.4683105 15   1       0.3333333\n",
       "3  3 0.4878556 15   1       0.5333333\n",
       "4  4 0.4954626 15   1       0.2666667\n",
       "5  5 0.4983212 15   1       0.5333333\n",
       "6  6 0.4993811 15   1       0.2666667\n",
       "7  1 0.4223188 15   2       0.4666667\n",
       "8  2 0.4683105 15   2       0.3333333\n",
       "9  3 0.4878556 15   2       0.4000000\n",
       "10 4 0.4954626 15   2       0.5333333"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# function for generating data\n",
    "draw_samples <- function(p, n) {\n",
    "    samples = sample(c(0,1), prob=c(1-p, p), size=n, replace=T)\n",
    "    prop = mean(samples)\n",
    "    return(prop)    \n",
    "}\n",
    "\n",
    "# replicate the dataframe for each subject\n",
    "set.seed(1212276)\n",
    "nSubj = 1000\n",
    "df$nObs = c(15, 15, 15, 15, 15, 15)\n",
    "data <- df[rep(1:nrow(df), nSubj),]\n",
    "data$subject <- rep(1:nSubj, each=nrow(df))\n",
    "\n",
    "# generate data for each subject\n",
    "data <- data  %>% \n",
    "    group_by(x, subject)  %>% \n",
    "    mutate(obs = draw_samples(y, nObs))\n",
    "\n",
    "head(data, n=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so we have plenty of variation in our subjects. Let's plot the generated data against the underlying model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAAOVBMVEUAAAAAv8QzMzNNTU1o\naGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD4dm3///8abOaXAAAACXBI\nWXMAABJ0AAASdAHeZh94AAALw0lEQVR4nO3diXaiSBhA4XLQRDsmCu//sMPiwlJIFdT6c+85\nnYXuoJNvCgHRUhWJTsW+A+Q3gIUHsPAAFh7AwgNYeAALD2DhOQHm/5J0A1h4AAsPYOEBLDyA\nhQew8AAWHsDCA1h4AAsPYOEBLDyAA1RGvG2AAwSw8ADOpZVSAOcSwMIDWHgACw9g4QEsvFVS\nZZ3rO2IewDatgSq7nN8XwwC2aQVTWcYVBti8VUwAZ5M1VDnI3x37mIFNUdf/uhgt2wnwklQ5\nyejHfLdsU7w+vD8Plu0SeIZz6cfCB7BhZpy6HzzUpbwX3ccsNMuMVpJ/qx9MD11e7pRBlsDP\nh+DXsv+a8gO2U2pZVwIfDnGFV4zgQsAINld6k67b1uYF/FywE+DheLWGOgyyvIeuAnjuX4w3\nxktSh0lGP+Y7NtG6v9U90s64GbtvvcsrswcuKtHAc7tRZpzzP2h9Fx1lcSarfwYr9zNZeuDZ\nfeRND6YxeXd7Lnrq+HHkvr5Yh8UT/oEbY+pxx5zrt7UAh21wwmJ+h0q71Pd9c94Ogd+npDS4\nn0dpzKG4sl0CN4ZT3Kg7u97aH3DjW7YfXoviHsj4bX/A1aEbwr1949j3yGe7Ay6fvJKHba+d\nAdePuzvCbdoT8GO3CuA4K/Hdc6eZERxpJV57HhO1rHvy3Qfwe/A+P++FN39go2fuuy/6qBme\nklqZdODX6arhmAU4/ErW9VGqN3hHm2SAw69kXR+kZgbv0o8JSypw+YkX4BgrWdfctTePL3a0\nuzyTQOCFwbuzxAG/n+WFt0kWcAnvOEnAfV14H+UN3LvuhsGrL2vg9+WRveur4B2UM7DmBfds\nm8eJAH4uQHdaxsDj8QuvroyBmwtgu0tgm2/g1Zc18KG7QLLkoXe+nIEP5a6urlpXzsD4GpQx\n8PsC9ig3n0n5Apc7u8B5ZdkCl3u7gn1luQL3X6MQ/tYzKlPg12Xs8C6UJ/DDt9rTxVUryxL4\n7QvwUjkC93wBXipD4L4vwEvlBzzwpaWyA8bXrtyAn8/9hrvFzMsMGF/b8gLG17qsgPG1Lydg\nfFeUETC+a8oHGN9VZQP8eklogNuSVC7A+K4sE2B815YHML6rs5hW5/ntdJlvYHzXZz21XRFv\nYix8V2QLXESY+YwDpA1ZAhcRprZjA72lrcD+J4jGd1PrJqcMuJOF77asgN9zVL6Xma1kdfhu\nzA64KyAwvluzPkwKOoLx3VzSwPhuz+JM1mCa9xA7Wfg6KOFz0e93R/Ky+p2ULjC+TkoWGF83\npQqMr6MSBcbXVckAl9pv8N1aksD4uitF4N57PztY8c5LEBhfXUrpv1MLv/z0gPHVJgYYX5Py\nBd63710d289H9Vf9filVnKuG8K84dZC9ZdWXOt2qB/D9W6nvu3aViQHv27eq1Rq0W+18VW3n\nhvCkvlvI/rKaWhX3B3DRLD5q15gW8N59a8JmfJ7VtR7EP1X11/C1oi1kf9npXp066qq6NF+c\n1T/dGhMB7mZe6B0q7RS4OrZPt7dPxd6ul1OH+doUD5c1A/0B3/yd+tKtMA3g0ewp+/Wt/qnf\n6ldd6q9O3fb4Sdt+1Czr/jyXT0sCeDo/zsYV5tu9frg9q/qx9Vsd/11vQ0zdsiyB9+vbIN7a\nbW23dzzFfC6bbqL1pQA8mQBp+z3Kt996JP5WDdxvdT+Ngd/LTs1Xl275udnJ+qkXaUoBeDSC\nd+3bDMf2eOespo+3/WX9w6R7e5hUHztrSg945771btZP+/m7HqW/o52s/rIv9fXau761f6Fd\nXRLAFQPYW2kAVzwA+yoVYF4E7KnEgPF1XRrAz9OU+DpvYvOv3v+uTto9bvOV2PYAxtd9I5v7\nsT3O6o61167EvnYaQnx9NLL5VufmyGrmrIjhSqwrmcLMWyOb9+nrDSuxjnli/ZUCMPPEeky/\niT6r7w0rsQ3gmVxMCjXeyepOW6vitmEltgE8kwfgqrrU+9HHs/4KPeOVWMVMz3N5AQ6/kvdU\nzy7ui6TkAFfw6vIArNSn63sMV2IX5yhnA1h43jbRt9Nl+0oMw3c+f4/Bd2UlDLCPhleKjzL+\nlc/8w3CbaIBnmrwYYNBG4B9VaJdvvbVp+M40eTHAsLXAr32ss82dAdh9c8CPHeDmKV3V+342\nPXBh5Quww8qZHn+tHh8eyq/vZ4t8ogPfuWZG8Au4//lTACeaIbDVWziofjZ3ZjUwl1LON/MQ\n/PjQ2zSbPwZHAsZXn3YXWgOc8mMwwB/THCLltZPFFvpzumPgjYdJr1coBtlEM4A/5+Fc9Dno\nYzDAn/MAXKi/k7rdT0EufMd3IS/PB1cXda3uQS58B3ghP8DX5g21+ptoXxNEs4sVopHNl/pp\n3rvltwfsbYJoBnCIRjaNbPtuW+8L371NEA1wiMY212P7Vh+9Z5Nm5w+eX4lRbKGDNLLRvC5Y\nD/x8DF4/QTQDOEjjnazjdfwvdBNEF8NhDHC6jWyOShWX4etWPE0QzRZ6ORe/nbHN7Vwo9dU/\nzeFpgmgG8HI+gOt+z0odf17f+pk/mAFskCfg5m1M54+DXQLj+zlfI7g+TDr23h7eywTRABvk\n7TH42+qphjXAbKFN8gDcvPr7n92rv6crMYgBbJIHYPU1OQ62X4lBAC/n5jXT4/focLGS5dhC\nL+foXQ/iXJPFAF7M1fuWRAFmAC+nB1a9j2YBnFiHmbq/zQkY34/NbKIfl0ObFwOYAWwQwNKb\n2cdSlr/taMD4LqXfhc4AmAFsmu53lAswvgZpf0mWYuGBGcDGASy8jIHxNUm3j2ULFhyYAWye\nt0t2fK4EYPPyBcY3WKGBGcCBA1h4gYHxDR3AwosBjG/AwgIzgIMHsPAiAOMbsqDADODwASy8\nkMD4Rghg4QEsvIDA+MYIYOEBLLxwwPhGCWDhBQPGN04ACw9g4YUCxjdSAAsPYOEFAsY3VgAL\nLwzwY4YngMMXEhjfCAEsvCDA+MYLYOEBLLwQwPhGDGDhBQDmIDhmwYDxjZMBsG4G8GowVxbA\n6bYMrJsBvPd5cSVsoaNmC1y8JiK1A8Y3UpbAvZlHAc6ircCLM4CzhY6bHXAx+mOwEgZw3KyA\nB4PXDJgBHDk74McM4K+ZwJdXwgCOnPVhkuVxMMCR8wzMFjp2Fmey+jOAV1bA+MbL87logGPn\nF5gtdPRCAOMbMa/ADOD4ASy8AMD4xswnMAM4gfwD4xs1gIXnEZgtdAp5B8Y3bv6AGcBJBLDw\nfAPjGzlvwAzgNAJYeJ6B8Y2dL2AGcCJ5AC5fHwCOnydgfFMJYOF5BcY3fn6AGcDJBLDwfALj\nm0BegBnA6QSw8JwDl3XdV/imkGvgsqv5EuAUcgxclm9hgFPIHzC+SeQWuAQ4tbyNYHzTCGDh\neduLBjiNfB0H45tIvq7oADiRABaeJ2B8Uwlg4QEsPD/A+CYTwMIDWHheXrqCbzoBLDznwIc6\nF+skN7kGPnS5WCu5yDHw4YBwWgEsPLfAB4BTixEsPICFx1608Cym1Xl+O13GcXC6WU+MVRTT\nZU42A+QlW+CiAjirVs4f3F8GcMptBV6cIJriZgfcnxjadIp3iprd/MHvLwHOJDvg3sTQ/SMn\ngNNt7fzBfV+AE24l8MAX4ISzOJPV39EqBqeylIP+c7GSHG7Nh+KHkhl8/3FrXgJY3K0NA1jc\nrQ0DWNytDUsGmPwEsPAAFh7AwgNYeIkAj677CnCDAW8q9H/boDSANee3Pd9guBsL/t82bJ/A\nRcAbA/hZ0N95YOB4Aez/xngMbgs7pgLeWtibG7dD4MAPijwGd4UcUkXIIxeA20L/BhjBYQv+\nCwA4aGE3mu0tBrwp9qLJXwALD2DhASw8gIUHsPAAFh7AwgNYeAALD2DhASw8gIUnFfhL/VXV\nnzrFvh/Rkwp8V8eqOjXKO08qcHVR1x91jn0v4icWOPLz7MkkF/hHqZ/Y9yGBABaeXODieGQT\nLRi43sm6qkvsexE/qcDtYdJR3WPfj+hJBX6c6PiKfT+iJxWYHgEsPICFB7DwABYewMIDWHgA\nCw9g4f0P2CnKstD43KcAAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data  %>% \n",
    "    gather(variable, value, y, obs)  %>% \n",
    "    ggplot(aes(x, value, color=variable)) +\n",
    "    stat_summary(fun.data=mean_se) +\n",
    "    stat_summary(fun.data=mean_se,geom='line') +\n",
    "    theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's plot a couple of subject just to show their variation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAAPFBMVEUAAAAAv8QaGhozMzNN\nTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD4dm3///+sNJpQAAAA\nCXBIWXMAABJ0AAASdAHeZh94AAAPtUlEQVR4nO2diZaiMBBFY49L27vw//86ijTNEkIqqSSV\n4t1zpkfBImVd2SJE0wLVmNIJgLRAsHIgWDkQrBwIVg4EKweClQPByokUbERSR3Z5gODsxCoj\nGooMF7iJHwsumYcdCI4GgiftRYaLLqHs7DK1FxkuuoSys8vUXmS46BLKzi5Te5HhoksoO7tM\n7UWGiy6h7OwytRcZLrqEsrPL1F5kuOgSys4uU3uR4WzpvnAtKIXglztMi9qrYLYCphD8MvyJ\nZ6eCX0SvwRDMgGTBHRAcBwSnAYI9wUFWJBCcBgj2gy09CI4mheAU2eUBgn3gSw6C40lwHvzC\n15W1W8F8oC960l5kuOgSys4uU3uR4aJLKDu7TO1FhosuoezsMrUXGS66hLKzy9ReZLjoEsrO\nLlN7keGiSyg7u0ztRYaLLqHs7DK1FxkukjqyywMEZydWGdFQ0XCQHAhmpWma0inMgGBOmkac\nYQhmpGnkGYZgRiBYORCsnadfUYohmJPf9VeQYghmZKRVjGIIZmQiVYhiCOZjblSEYghmw6JT\ngGII5sLusrhiCOZizWRhxRDMhEPj76lTEdMQzINbXtf7UcYwBPOw5a5YLyYEs7BpDoKrZlsc\nBFeNhzjsgyvGSxyOoquldF+GEwiOB4J1I9ovBEcj2y8ExyLcr5eh4x3bY89w3SgQfBz+TB97\nhutGul8IjkO83yDBtHDVKBQ87IP/Pdi5YPl+yYKPLTbRAxX4xT44BgjWTQ1+ITicKvxCcDhq\nBA+9V8fRY0K4Uurwi77oUCrxC8GB1OIXggOB4G3CLlISMESCuGEaXJQTHHaZoYBBMASkQKCY\n4LALhQUMcyIgBQoQTEVAChQgmIqAFChgH1xjCgQ0CT4cDuHp+GfgOpB3pJAnuwXlNtHDHwvr\nxWjus1aiDoccNWxc2TlSyJPdklKC+zrZXa0Xo+lmWaMOB9YariyqcWXnSIE5O3/KCSaoOhAI\nTciWwzLrRkZ2BAoJbn7L8mfYpyzNIqqdxwcmtLUwL2mOFHYmuLFU7G/uejH6z4V9y85ZQWt2\nzUZ2u9wHL9/VrHCrK6NlejMscM2wbTLpoHftQzdq06XKIbGM39SCJ8UYitaMZpFUWR5t4jix\nsme3nDVrspCqMNIKtq8PzTBz9MyDZuWxO2i948mxti4t1tO1MSWL4OnExvnUQSLB5LarIqlg\n+8qxGJTVs63G8cwRtSp4ZdX1aLsm0gnuSufhN7HglX3wWnY+TVdFIsFD4bb9+pYvdMXvjqKn\nL17Pzq/pmkgg2L1a5Bc8ezX9GLhmv+yCt8pnLZZPBZevCTEccoZTtV8WwYQt3kqxPIYCDIrq\nGE7Hws5f6/bLIbg/VPEp31qxUgompGdro3K/DIJZzja2ymif71P8UXoBrmrXSx9lp50O42BI\n3QXr5Uou2Pv15AaEQ7y7sO1v8h+F+wt2lctdyvBN+zQ9Whe2Br90wcf5GszTXZBM8DS9yO6R\nCqEKPi420UzdBWH2PQ3TAtrqro11EC6YOsrORrXSCZ5FeIXsVvBskB3KQXjYcdRWZKIj470K\nno/gQBC8Waxsgj1SKfcrOOwQBT8hhT/xKFbYkRS/4aK/Y8UO/TQpaA0O6232iQ3S4AhS4nVA\nkODA3gxWw9r00kfZaYMEh30h6BcbpsT+rZY6vZmuiw67ZsM3lsuwRr07Fry49k+l3jyC/UsX\n9q0+g2GtenkEu84oiKeTYdflBBtufk+JwhZQAwyCXX0C1P6CrIL77DTr5RDs6tWj9/g1jme+\nUd5tqemPdMAm2AVhgRDMjLQ1eHroExTlHwTBXuGc++CsgvV8o+Agw1E0cZGN5RElitSWer8C\nR5vNKXgHyBM8uArftIM/IFg5AgX3rsL33WAEBCtHouDOVcD1kvyJKECm4KDz0x2c1AagTDAM\nz5EoOKwLcRcdj3QgWDkQXAnG2J+ZDQUSBWMfbEGV4Hp/MysXtQveKzdz6v4/me/262LM8do+\nFH4fz0+Ro2ntxZx/2l7w7dWY15t1kRAsiot5SPu5e/40HdeHwrN57USOp91Vm+OtF3x8TD5Z\nlwjBovh82Guv5vO+En+07fdDX2e0Ezmedr6156fqtn17PLiad9sSIVgWp+evcD8e/ny+nZ8y\nh03xdNpjRe/FP+aZi22B1FF28AvgaXk3X+2Xebs/Oj+3x79qu7+Wac9/v9OXBNzh34bcHwz8\nuN13t1dz37e+mtP7589Upm0aBFfGq/nptrXPo+OlzN9py020nYD7gyE4JV/3NfGrfYj7am/n\nueC/aefHo7fn9OvjIOvjPslCuGDqKDvAj9PzfOdqlvvb8bTxadKtO026nztbyHSHP/Dm/XEu\n1D621eb8NTvIGk+7mMtwdP3TzbAuDoKVQxc89gvB4iELnviFYPFQBU/9QrB4Fobe74dn7Xly\nQDYaZWc+EBqQzkzw7dQdhj9PxYACZoJf76fMd8ErJ82gPmaC/3o3y6QDuIFg5dg30VfzWiYd\nwM38IOvZq2mOP2XSAdwsNsVv9+Po09V+AReoD+xrlRN7TZZI6sguD4ujaFoa5crkoo7s8hAt\nmD+lWPyv+i9BYcFPfs5vvuGiSyg7u0ztWafejKdh2SWUnV3cYiJfiE10Kgg5uW6mixT8YTy/\nEpRdQtnZbeG8HTZU8HCMdfUMF11C2dltsHJDe38A/PhK14yer7e3DDf9LYo+yC6h7OzsbIy5\nbfo/veXh+Xp7adMtQd2CB+xr8CB4/L+zPVp6i3DRJZSd3Rb2LXT/Z1hzKXf4B/S3yC6h7Ow2\nsR1jLQRvLVGI4Jc7PEtKJZgrvcicLIJr2Ae/DH8YSCKY7fPHJbiugyzxgl+krMGxp0nDDWwF\n9sGCBb+I2UTT25s+vRY8yILgFMyaO5rvs/m5nX0vfOcroeSDrJdWjeB782/ms735Xvi+izU4\n0RFCHpaCPx/jLWETPeLlCcuySgu+mI/H0B5fmQWLP4pWswY/zHaDMXle+A7BVEpf0fF56kaC\nyP1tEnqyUjFrzjpQiytcdG+v7OwytTd7evqkhYsuoezsMrU3fXoy5vhGuG9FdgllZ5epvdnz\nn+vRmIv3/f2ySyg7u0ztLSd9XY05fXiGiy6h7OwytWeb+IMv/JMhQPDX/TTpZB093BIuuoSy\ns8vU3ux5tw9+JeyDJVJHdnlYHkWf3klH0RKpI7s8zM+DL8TzYMZUQApmhqhDN0CwN/Sf7Toc\nDvHNxl4hFJ/BTqD/8N7hwGEYgvNA/+nMw8Fi2Iz++gHBefAWfFjhOReCpcK0Brf95dD+QHAm\nmPbBECwWpqNoQ6w5BOeC6beNIVgoXL9dDcFCYftxcmLJITgTEKwbNr8QLBO2XTC14hCcBb4V\nmAoEZwGClQPBuinnF4KzAMHKgWDdFPQLwTmAYOVAsG5K+oXgDECwcsQLHv+q++wX3iF4m6J+\nfQwdhz/Tx57huweClVOfYFr43inrly542Af/ewDBm1Qm+NhiE02jMsGjx57hO6ewXwhODQTr\nprRfCE5MDYKH3qvj6DEhfNdUIThduHqK+4XgtECwciDYjuNm6bBZbBHE5SdduhciBTuGOwib\nRW+HBwi24hiwJGwWvR0mINjK2q+ce0FvJ8Vb6BtIt2hvBAvmm0VvhwcIXkHLPhiCV3DVPWyW\n/fXOhoKYjHzUrM/Kh0TBjaMYjhE6qYN3NpP/CPhm16zPykdBwWtvuHEUwz6+39aslWbmD7yy\nc6mapdCsz8pHOcGrb7hZFGNtgM71QTs3GW2arYIJn7EE2TFSTDClTo4or1lLmtUnfNktt9B7\nFLyoU+MuhqNMoX5thsM+Y9MUFsdYKvfBy3e1VriepnUXw1Emut8+YqRiI7vhJdspzD81Zfym\nFjwpxmH5bBHQ/M6Ly8vJ2O+fYZ/s+nivcysJJ8FtasHO9WHdb1qmfinZ9fFeZ89C/OYR7L24\nrH4n2fm37NnBuQvBzl2ZjQxVWfiF4MDwrnRi/S6y827bT7AUv6kED4UT6deaHc2wX0MCSCA4\n+IQvfVH6o+Wo1ruvKDZfrVZwxAlODr+u9Pzab2b/xywrBxyCg7bHC5LVxDs92tmt89WqBPeH\nKpGdEyn9eqbncXLbrD0hLyoX8YLJ50J2UpXkNz2Ozqf5/PVLDzwaywR1lJ12OoyDoXdm2Em7\nAvv5pe9Y115fl+DZHYXzm89YBKeryN1t03j6dedhm7d2eZhncxkgCz7O12CWL8JSXtpIufCK\nvFu1X+Dp21wGqIKPi000x3c/qf0yGKZsjXUIZhxlJ+mlq8SLn8kbXcucigXPBtlhuqpSvmDS\nvlmSX6Lg+QgOPILTFoR8v4PntPW5NQt+QgrfJHU96NfDb09xzq9YcDt7yCFYVDk66Gvk1kWa\nBYHgJQErZLPyuDzUUXZabsGyyvEkYIX0/BIiP6XvTRJWjp6RLu8+ML/vEbMDwTaaxQP/IGHv\nqLBgYdUYCJLV0EPSwyE48CbbxHdfx9H4XJezCEo+bA8dBsGB98mnvr0+jrDsBL6neMGBI10k\nHyAjirDsJL4nCLYCwX/hEBwZlRbsg+1gHzwKDzngFHnEOSYsO3nvie08mPjGpNVBLXwdHSRl\n8JsLxp4sgjT4zQZnV6V/x3xco4AAb180x+XlgBXmLxs4bhABnHB/m0S8vAWkhv3rQuoNXCAt\n/N8HU2/gAklJ8IW/6/4e+M1Niis6arhtdjckuWRH/k2V+yHNNVnib6rcD4kuuhN+y92OSHVV\npexb7nZEsstmJd+QtSfSXRc9HXIorhkQTMoL38Xer7MnqMMokX7iXejtOrsiYAiH1v8GcL+B\nO0FC0gqWeJnhzgi4AZyyBgu8UHhnhAv2GUYJgouTdggHCC5O4jE64Lc0AYORksLhtzABg5ES\nw0FRAsaqnIRv8m/7JdxRlWSXB+IwSvOR7rb5F5RWrijZ2XGQ/AMlu4Sys+MAgrO0Ex4VC46S\nlAPByoFg5UCwciBYORCsnMSCqd0if4E52pKdHQ9pBS++SvYOpAfR25KdHRMyBR8DgvIJzpMd\nExn2wUHVCCwhHdnZxaNJcNBeTnZ28aQXHPZpD4gKCpOdHQMCBQfursLCZGfHQHLBIR928pfO\nQ0PUKNnZcZBacOh7yrOOyM6OhdQdHfkCA0ooOzseEp8HB23OusiwxsgBcrNjAn3RyoFg5UCw\nciBYORCsHAhWDgQrB4KVA8HKgWDlQLByIFg5EKwcCFaORsEX89223+ZcOg8RaBR8M6e2PT8s\nA5WC2zfz+WGupbOQgUrBxS6fEIhOwR/GfJTOQQgQrBydgo+nEzbRT1QKvh9kfZq30lnIQKPg\n7jTpZG6l8xCBRsF9R8eldB4i0CgYjIBg5UCwciBYORCsHAhWDgQrB4KVA8HK+Q+7cG1+xqzL\nsQAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data  %>% \n",
    "    gather(variable, value, y, obs)  %>% \n",
    "    filter(subject <= 4)  %>% \n",
    "    ggplot(aes(x, value, color=variable)) +\n",
    "    geom_point() +\n",
    "    geom_line() +\n",
    "    facet_wrap(~subject) +\n",
    "    theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, since we have few observations, each subject is a poor representation of the overall data, yet, the average reflects it well. Now, for the interesting part. Let us try to fit an logistic growth model to:\n",
    "1. The overall sample\n",
    "2. Each subject individually"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting the overall sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 1.0958580 2.1848688 0.9864729\n"
     ]
    }
   ],
   "source": [
    "# Logistic function\n",
    "myPredict <- function(par, x) {\n",
    "    pred = par[1]/(par[2] + exp(-par[3]*x))\n",
    "    return(pred)\n",
    "}\n",
    "\n",
    "# Calculate error for minimazation\n",
    "errorFunction <- function(par, obs, x) {\n",
    "    predicted = myPredict(par, x)\n",
    "    error = sqrt(mean((obs-predicted) ** 2))\n",
    "    return(error)\n",
    "}\n",
    "\n",
    "# Fit the model and extract predictions\n",
    "overall.fit <- optim(c(1,1,1), errorFunction, obs=data$obs, x=data$x)\n",
    "data$overall.pred <- myPredict(overall.fit$par, data$x)\n",
    "print(overall.fit$par)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alright, we see a pretty good fit and recovery of the original model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAAPFBMVEUAAAAAujgzMzNNTU1h\nnP9oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD4dm3////s2uQOAAAA\nCXBIWXMAABJ0AAASdAHeZh94AAAM3UlEQVR4nO2di3rqKhBGyYn3XVsN7/+uJ+QmUaIhDjAM\n//q+VkVL6V57JoCJozQQjUo9ABAWCBYOBAsHgoUDwcKBYOFAsHAgWDgkgvG/hC8QLBwIFg4E\nCweChQPBwoFg4UCwcCBYOBAsHAgWDgQLB4KFA8H+NKkH4AME+wPBwoHg/PByBsH5AcHhO0kK\nBIfvJCkQHL6TpEBw+E6SAsHhO0mKh7OmJdxAyIHgjvXKmp6AY6EFgjtWC2uaqiUjwxCsvbJu\n59coDjogQiDYI+sOL+wMhx8WDSvc1C32/fqpLXfBb7JuM6dtqUaiD3Mjn93U07fH7awte8FW\n1n01alNZpBjpFiC46cP3o9FealXtWsQKrh1tqzrhzHhU7f09G32i82sURx/mRjwFj4fgqe0/\nA0fB66a505xpbeLdjXw7vlhsiOA6iwheNynuQ9Yj6woXPDYIEPxY7azPujsLkkFGoEjB85Xs\nG2e7Oe9fzJPSUrTjUOv091E5/YjD4C+41hkIdk2Y+sXui/gPRhdeTjPMCHjsZNk7WMx3sp7n\nxP0j52aV92E1K71C96LnOxPdrWuzebDqnXWzeaPBIFvwEMWvcmdCfbMuBCfmg91XmZllXS8k\nCn4Yfk7Mi5GaVVB6IVSwEemQm2pACREouA3czm77bWjJa2FDizjB5sC7M3atKXLqMaVEmOBu\nXjXqLd1thyjB0PuKIMGDXp3hhnFAxAh+6IVgGyGC53qze0cgICIEP+nVxS56HWQkeGm36VXv\nu1eXRv6CnXoheCR3wf1bgd3deVKG4J68BVvv4j8fcyG4J2fBnd7+7uuUCoJ78hX8Vi8YyVUw\n9K4kT8GWXg29b8lR8Ewv/L4nG8HTabDQ60Uugsez6KDXk0wEP86ThF4/MhMMvb7kIXjSO1yN\nAr2ryUPw8Dk4wzwLej3IQ7AJ4Wa4Ch96vchEcLPrr8NvoNeTTATvqgZn4WwiE8GjXwj2JQ/B\nk18I9iUTwTgRditZCK5wpvNmchBsPiYUfjeSgeBquhol6K8RSg6Cdf+uPs6y2gJ/wf3upIbg\nbbAXPPmF4E1wF9wfgLu7ELwF9oLNN8yutsNcMPx+C2/B8Ps1rAX3529A8DfwFwy/X8FZMPwS\nwFgw/FLAVzAOwCSwFdwggElgKxh+aeAqGAmaCKaC4ZcKnoJxACbDo6zO+PC1jVpwgwAmw7u0\nXR2+MBb8EuIruI5Q+Qx+CfEUXEcobYcDMCXfCqYvEI0ETcq24pQBJ1lF+FXK/UjRLzi9BD9q\nVD7a1nWymiL88hXcA8G08BE8PQ4mWJLfu9p3t3v1p3+PStVnbRT+1YdepNWmj+pw04Pg+0mp\n051qFLwES/KrW2tG2q31fFUdZ6PwoE6dSLutVa3q+yC4Ns17qkF47GTNyrwHmWTJ8tsqNPF5\nVtc2iH+0/jP6OqOdSLvtcNeHXrXWF3PnrP4RDYLTXnQjbYtj3xdNN3dv18uhlzml4nmbCfRB\nvHlOHYnGwEqwML/6n/rVv+rS3jv0+XhU2313tPVfYzsJjASL89tOs05tsm2PrSe1/3e9zWW6\n2kQLnj7FjqIzJpzUrcu1/ez4VebY9pqiyWAjWNwB2PDbRuKvNuJ+9f3wLPjRdjD3Ln372Uyy\nftomGtIKbqx7Av2acOzWO2f1ery12+xl0r1bJrVrZxrYCBaYoLWZZv10t6c2Sn+fJll221Ed\np9n1rXuCagRMBAv1ywAegoUmaA6wEAy/4eAgePQLwQFgIRh+w8FAMBJ0SNILRoIOSkrBjSkN\nC79hSSi46YHfoKQT3EBwDFIK7kvlDAW8KcbBmkQf1JdMcO/XfJlH8v0WJ7irZtYMNUUL8Fuk\n4KaC4NAkE1xZgkvwW5xgK4KL8Fui4KaC4OAk3OgoKoC7Xbsl6K9IIu5761ZlQQfgnoVnowr+\nd1RKH/xO+dq4VVnpx+f1y2bctVswHFHwfd+d6Nef7Lm1k5U0pmBdEX4XBQ8nuJt/cGU9puSp\nv5M6m1P7PE/L/UKw8GpXjZvhWTV8GyxPjyl56u5x/cQXnayjqcopV7cQwZNg+5YYCI7BSsER\nrvAfUvRZnb7oZBW7kgQvzKJfBAdQ/DzJ6q+bUPXti07WMdUE3vLD2eGcQjsEhz4Ga31p59H7\ns99HRGx6u7Aswc6drASTrHidFJWhDa41cPxlUrxOqr4mMMWvzwMee9Fq4otO1tCUskk5AcHC\n4SG453a4fN/Je0o5j+MBJ8H6rrwM+wtuIDgSC25Cp+jyMnQq3G5+VO1s9+rkHeUFcCqWJlnn\nLzr5TIEZOhVuwbWX322C4TcKaTY6EMDRSCK4yAxdpfm1thtls7WTNRSZoUsSXGIAcxAcq5Mi\nM3Rpgovzy0XwOUKKLjKAmQg+RzgGl5mhmQiu1d9B3e6HkCe+Q3BMXk+bvairvoc88b1Iv8OV\n7v6sOV363QteBV9NRRc7RRMXiC4ygKuq2qiYWPBR/ZjiAb+WYOoC0SUKrqrNhokFG7NduZfH\nie/UBaIL9LsoePhI//7+7CRL68H4zzueezl7umvzEKyv++6z5q13kxbrBy938paiBFduhmeH\n86BHwc+nSY92R8HK/fT7szOennNcF+wWPB6DvQtEl5ihlyJ4ZunpQgf9eKBeX+2+KMLF8yRr\nf31+hatAdD0PY1/BpfldK1jZKfuRv52C7dd6CN4rVV/m160QF4guMYCXZtGTHTU6fjRbYb0Y\nwdpfsL6da6WO9jYHbYHoIjO0XlgHvxPsOAY/Wr4Q3PJ7Vmr/Mz2krR9cZIbucCyRlPvbw7t7\nkvX8Wp9JVs9NLa+DvxRcaABr91bltOmv7IfTAyuClb1Msn/Ud5nU8tsuk/ZWfWLKAtGlZmhN\nsBfttRhd/KnuGHzyLKwGwSvgIdhc/f3P7+pvr19drl8mgtXxZR1M+asLDmAmbxd6B6+rk2VK\nFpyIbXG/tRP4jU5UwQjg+ECwcKIKht/4xBSMAE4ABMci0R8eU3C3EoTguEQUXHYAQ7B05Asu\nO0MXIrhcv/IFF56h5QsuPEOXIbhcv8k+Ojma4LIzdLoPP48muOgMvUtXviCq4FL9Lgi2T6gM\nRizBhWbonZv+SVGCi87QSynavpwhFBAcA/mCC83QEwtzrPcXJZAQSXDZAayX1sEQLAjXXy9G\ncOkZ2uD884P7jSQYAQzB4pEsGBlaO//+APWgX39HjE4QwFr224UQnI4YgpGhExJDMAI4IRAs\nnAiCGwhOSATB8JsSCBYOBAsnvGD4TQoECweChRNcMPymBYKFA8HCCS0YfhMDwcKBYOEEFtx/\nhi4EpyOGYPhNCAQLJ6xgZOjkRBAMvymBYOGsEOyqAK5ntbKWOkGGTs9nwa4K4Nbtu04QwOnx\nFVxPhUghOAs8BVuVRz8LRoZmwLeC31UARwAzwE9w/fT1vhMIZoCX4FnwfhSMDM0BP8FDBfCp\nEvjbTiCYA97LpPXrYGRoDkCwcDx2suwK4PqzYGRoFoTbi0YAswCChRNMMDI0D8IKht/kQLBw\nQglGhmZCUMHwmx4IFk4gwcjQXIBg4YQUDL8MCCMYAcwGCBZOQMHwy4EgghHAfIBg4YQTDL8s\nCCEYAcwIasGNhmBWQLBwAgiujGT4ZUIIwd0dCOYBBAuHXjAyNCsCCO5uIZgJECwcWsFVCzI0\nK0gFVz0agvlAKbiqJsMQzAVywY0RDL9sIBQ8+YVgRhAKbh4pGoLZEOYYDMFsoBQ8hnADv3wg\nXSZ1hpsGAcwI2o2OpkVjkcSJEG/4wy8jIFg4ECycAILhlxMQLBwIFg69YPhlBQQLh/7SFQhm\nBQQLh1bwroWiQ0AGqeBdD0WXgAhKwbsdDLMDgoVDKHgHwQxBBAsHgoWDWbRwPMrqjA9f27AO\n5ot3Yay6fm0jSQMgCL6Caw3BWbGxfrDdBsGc+VbwuwLRgAF+gu3C0J8LRAMG+NUPftyF4Ezw\nE2wVhrZXThDMl631g22/EMyYjYJnfiGYMR47WfZEq55tZamt/Lf5J7l0HcIJKWlH+B+6Dg0E\ns+qaHghm1TU9EMyqa3r4zxLAV0CwcCBYOBAsHAgWTkrBTyd7Ufceqt+goyYnoWDHpjZl74F6\nDjtqeqQKrkP1DMF+hNMQUnBGQLB3zzgGexAwzEJ1HbDvEIgUHPI4iWOwD8GirA62mIFgD4L+\nOyGCO5JudOTYOwSvJmAe7boP1S9m0YARECwcCBYOBAsHgoUDwcKBYOFAsHAgWDgQLBwIFg4E\nCweChZO34KP60/pPHVKPgzF5C76rvdYHYxkskLdgfVHXH3VOPQrOZC44t7ff45O74B+lflKP\ngTUQLJzcBdf7PVL0OzIX3E6yruqSehScyVtwt0zaq3vqcTAmb8HDRscx9TgYk7dg8BEIFg4E\nCweChQPBwoFg4UCwcCBYOBAsnP8B1ys2kYPzDfAAAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data  %>% \n",
    "    gather(variable, value, y, obs, overall.pred)   %>% \n",
    "    ggplot(aes(x, value, color=variable)) +\n",
    "    stat_summary(fun.data=mean_se) +\n",
    "    stat_summary(fun.data=mean_se, geom='line') +\n",
    "    theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting individual subjects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how fitting the individual subjects will fare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>subject</th><th scope=col>subjectPred</th><th scope=col>x</th><th scope=col>y</th><th scope=col>nObs</th><th scope=col>obs</th><th scope=col>overall.pred</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>1        </td><td>0.3283832</td><td>1        </td><td>0.4223188</td><td>15       </td><td>0.3333333</td><td>0.4284447</td></tr>\n",
       "\t<tr><td>1        </td><td>0.3779846</td><td>2        </td><td>0.4683105</td><td>15       </td><td>0.3333333</td><td>0.4715567</td></tr>\n",
       "\t<tr><td>1        </td><td>0.3884094</td><td>3        </td><td>0.4878556</td><td>15       </td><td>0.5333333</td><td>0.4899402</td></tr>\n",
       "\t<tr><td>1        </td><td>0.3903222</td><td>4        </td><td>0.4954626</td><td>15       </td><td>0.2666667</td><td>0.4971675</td></tr>\n",
       "\t<tr><td>1        </td><td>0.3906641</td><td>5        </td><td>0.4983212</td><td>15       </td><td>0.5333333</td><td>0.4999174</td></tr>\n",
       "\t<tr><td>1        </td><td>0.3907249</td><td>6        </td><td>0.4993811</td><td>15       </td><td>0.2666667</td><td>0.5009506</td></tr>\n",
       "\t<tr><td>2        </td><td>0.3729196</td><td>1        </td><td>0.4223188</td><td>15       </td><td>0.4666667</td><td>0.4284447</td></tr>\n",
       "\t<tr><td>2        </td><td>0.4181209</td><td>2        </td><td>0.4683105</td><td>15       </td><td>0.3333333</td><td>0.4715567</td></tr>\n",
       "\t<tr><td>2        </td><td>0.4718946</td><td>3        </td><td>0.4878556</td><td>15       </td><td>0.4000000</td><td>0.4899402</td></tr>\n",
       "\t<tr><td>2        </td><td>0.5367985</td><td>4        </td><td>0.4954626</td><td>15       </td><td>0.5333333</td><td>0.4971675</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllllll}\n",
       " subject & subjectPred & x & y & nObs & obs & overall.pred\\\\\n",
       "\\hline\n",
       "\t 1         & 0.3283832 & 1         & 0.4223188 & 15        & 0.3333333 & 0.4284447\\\\\n",
       "\t 1         & 0.3779846 & 2         & 0.4683105 & 15        & 0.3333333 & 0.4715567\\\\\n",
       "\t 1         & 0.3884094 & 3         & 0.4878556 & 15        & 0.5333333 & 0.4899402\\\\\n",
       "\t 1         & 0.3903222 & 4         & 0.4954626 & 15        & 0.2666667 & 0.4971675\\\\\n",
       "\t 1         & 0.3906641 & 5         & 0.4983212 & 15        & 0.5333333 & 0.4999174\\\\\n",
       "\t 1         & 0.3907249 & 6         & 0.4993811 & 15        & 0.2666667 & 0.5009506\\\\\n",
       "\t 2         & 0.3729196 & 1         & 0.4223188 & 15        & 0.4666667 & 0.4284447\\\\\n",
       "\t 2         & 0.4181209 & 2         & 0.4683105 & 15        & 0.3333333 & 0.4715567\\\\\n",
       "\t 2         & 0.4718946 & 3         & 0.4878556 & 15        & 0.4000000 & 0.4899402\\\\\n",
       "\t 2         & 0.5367985 & 4         & 0.4954626 & 15        & 0.5333333 & 0.4971675\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "subject | subjectPred | x | y | nObs | obs | overall.pred | \n",
       "|---|---|---|---|---|---|---|---|---|---|\n",
       "| 1         | 0.3283832 | 1         | 0.4223188 | 15        | 0.3333333 | 0.4284447 | \n",
       "| 1         | 0.3779846 | 2         | 0.4683105 | 15        | 0.3333333 | 0.4715567 | \n",
       "| 1         | 0.3884094 | 3         | 0.4878556 | 15        | 0.5333333 | 0.4899402 | \n",
       "| 1         | 0.3903222 | 4         | 0.4954626 | 15        | 0.2666667 | 0.4971675 | \n",
       "| 1         | 0.3906641 | 5         | 0.4983212 | 15        | 0.5333333 | 0.4999174 | \n",
       "| 1         | 0.3907249 | 6         | 0.4993811 | 15        | 0.2666667 | 0.5009506 | \n",
       "| 2         | 0.3729196 | 1         | 0.4223188 | 15        | 0.4666667 | 0.4284447 | \n",
       "| 2         | 0.4181209 | 2         | 0.4683105 | 15        | 0.3333333 | 0.4715567 | \n",
       "| 2         | 0.4718946 | 3         | 0.4878556 | 15        | 0.4000000 | 0.4899402 | \n",
       "| 2         | 0.5367985 | 4         | 0.4954626 | 15        | 0.5333333 | 0.4971675 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "   subject subjectPred x y         nObs obs       overall.pred\n",
       "1  1       0.3283832   1 0.4223188 15   0.3333333 0.4284447   \n",
       "2  1       0.3779846   2 0.4683105 15   0.3333333 0.4715567   \n",
       "3  1       0.3884094   3 0.4878556 15   0.5333333 0.4899402   \n",
       "4  1       0.3903222   4 0.4954626 15   0.2666667 0.4971675   \n",
       "5  1       0.3906641   5 0.4983212 15   0.5333333 0.4999174   \n",
       "6  1       0.3907249   6 0.4993811 15   0.2666667 0.5009506   \n",
       "7  2       0.3729196   1 0.4223188 15   0.4666667 0.4284447   \n",
       "8  2       0.4181209   2 0.4683105 15   0.3333333 0.4715567   \n",
       "9  2       0.4718946   3 0.4878556 15   0.4000000 0.4899402   \n",
       "10 2       0.5367985   4 0.4954626 15   0.5333333 0.4971675   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Fit parameters of the logistic function for every subject and extract them\n",
    "fittedPars <- data  %>% \n",
    "    group_by(subject)  %>% \n",
    "    nest(-subject)  %>%\n",
    "    mutate(fit = map(data, ~ try(optim(par=c(1,1.5,1), fn = errorFunction, obs = .$obs, x=.$x))),\n",
    "           pars = map(fit, ~ .$par))  %>% \n",
    "    select(-fit)\n",
    "\n",
    "# Get final predictions for each subject\n",
    "fitted <- fittedPars  %>% \n",
    "    mutate(subjectPred = map2(data, pars, ~ myPredict(.y, .x$x)))  %>% \n",
    "    unnest(data, subjectPred)\n",
    "\n",
    "head(fitted, n=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've fit a logistic growth model to each subject, let's plot the overall fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAAP1BMVEUAAAAAv8QzMzNNTU1o\naGh8fHx8rgCMjIyampqnp6eysrK9vb3HfP/Hx8fQ0NDZ2dnh4eHp6enw8PD4dm3////9WFhj\nAAAACXBIWXMAABJ0AAASdAHeZh94AAANyklEQVR4nO2di3ajIBRFyZhH06ZpG/n/bx1BRUQ0\nkIDA5ey1xkZMKOmei4CIjAPSsNQFAHGBYOJAMHEgmDgQTBwIJg4EEweCiRNEMP6X5AsEEweC\niQPBxIFg4kAwcSCYOBBMHAgmDgQTB4KJA8HEgWDiQLA/h9QF8AGCfTl0pC6DBxDsyaEndTGc\ngWA/DofCDEOwHxBMmwMEEweCiQPBtCnOLwR7IbwWpReCfShK7AgEu1KkXgh2pky9EOxIoeHL\nIdiNYvVCsAvlhi+HYAdK1gvBTyk6fDkEP6NwvRC8TenhyyF4CwJ6IXgDCnoheBUS4csheA0i\neiHYDpXw5RBshY5eCLZAKHw5BC8hpReCTWiFL4fgOeT0QrBknEdHTy8EC8qbC+sBBNP2C8EQ\nvFcmqSjvfjIvIBgRvFMmyYDgfTJJxaG4VRm8qF3wocD7ybyoWzBdr4qaBVegt2rBNeitWHAV\n4curFVyL3loFV6O3TsH1hC+vUXBVeusTXJne6gTXprcywdWFL69KcI16axJcpd56BIcM3zZc\nVvGpQ/BTvV7OIDg3nkcvBMfPJBoutbOHs7bj9cLsDnnBbidfd2XtuaMgxUQFq5tR3PS6B6X0\nKxS/XLSdcXDTdOivGyMtQ8HjNDq3prNPUJ5H3inenjx306jN9HOWlp9gv4mwbkHZTtxuEJwW\nP8FmULZWzhO3DqKCG0uaUya74nczSiuq5xWVGsO7byMRyx8UT8HjKVil/RNkJtgjgm1CNz9A\nXfCwyTuCnQRPXt2c3QzCljgenufg4XXegrduRrEEqunMVGk/GPtLhIKe4M7r+Sz06hKV1mWP\nd02lHeKC86+iZdgODV39vNq3n/R3rgbpNoX5fUFwwzMWLGtl4XeqhudqLTo9g7IovV4jWfoI\nVp4jWeNZd2g5aWo3JZYWlF5QGosewlfa7UecHIOTrl5KgkX49rXyJDd1mdZgzL7Hwv8lqQiW\nbefzcGEo+6YuBPtyGFvK89Zx6mK5AcE9q1d7ul7veEl3kJqF3wc7yp9H9st/Low1Vy4U/jan\nXqSWxi/s9McHwY8Pxj4eoUpRvGAxqDG0lTWl6fV2XJiQ9td5vjPJVSg8sQ8pUk/rVLPmMQhu\nRPIxVCHKFiyGpg6qI/Ts3XtzF/b4ld27IP7m/Ffok0alSD3t9OCnXjXnn+LFlX0FKkS5gvuB\nx77rm0O8Wjj2Qwfi5d/989TLVFXxPE0E+iBeHGOXQGUoRvD8Qt4wqiz0ZlEb2/liP/yHfXav\nTn19PKqVW0ta/29MD0IpgvVLQGoG1SHX0B14dKfbK+vOrR/s+HX/m8u0pdUreHZFfqirD3nb\nFXywP1nX9q3jpcwxbVlFB6M0wfpVg+i/9W1+ukj84ULcD3+cTMFT2km8+uzTr6KR9d0lhaEM\nwWd1gajfL8Ku4Nj3d65seb7V0/Ru0kN2k7q+cxjKEMzPt1ZNZux+FHOr75foC3FRV7PTj9HI\n0tMu7KJa13/yQKgSFCL4Nps6VYreHChE8HmyW+md+q9ShOBWC1/o9aMEwV3TSgmGXk8KEDw0\nrG4lNa7yIXvBrewZDXoRwN7kLrjv+Qq3t04vItibzAX3Axs34k9GiUnWgsfqmfqjb2KSs+Ch\nehYbAn4TzUDIV3A7Vs9cDl0VrheCTQy/5QPBM87q9EvELwTr6OFLxC8Ea1D0C8ETZ4p+IXhk\nDN+hexQy66RAcA/J6vnJaokxx/JzE0yzeh4XV1s5Wo/gdpxUJ7fk/K4Zrkbw7PRLyO+q4GGC\nO+NsnIwX/AbSjAS3s+qZhl/7upejZzZsBstqPyT5CDaqZxJ+R1YiWAnWfwYmrWDtC1M9/Uoc\nBdO7w199YVU9k/S71opeCI6gOA/B8/Cl53elH2wRTOwcPHxn+n65dSSLeiNr+F/dnmvwax2q\npN1NGp6FocKXtt/6xqKHFbnN6pmq3/oED3oJjm7YqU1wO/qdpk6S9ludYOX3TP/0m5KFm69L\n15A7+S0g8HoVLQWPyzHAbwQMN4+jbKn3S4e8mokbIoDlwhvwGxPDzQe7im635yIvr9Xzt7O+\nZCjd5nNaDDfTalxvZOLGrdUFQ28kEgoe12UQO/AbC3sVfWUfb2TihApgua7KCxmURqLHWZqN\nrH4VLtb8vZGJG9rCOTX4zUQw559dO/p49Vtw/KV+sLZwzgsfL49cBO+USatWzqnEb22Cz/1q\nsRUtu5GHYKZ4IxMHugCWa6oUf1u3O3UJlgEMwTtgdfN3+nw/ky30AIbgqNjdPJiXYW/BcvG6\nqgQ/fXb8Gi7TpbfesHIsbhU981uFYPOp4+7EEfzNGmv6C7/AhhS88Rx2ciyfK+9MYMGqjXX1\nLIQPvV8CKyM5syZ4WNK/fz2bZKntjH/ece7l7LBM8xfcePn1FNwqwV4fK5KzneHoMA96FGxO\nkx7tjoKZ/fD2+TTBQMfotwbBAysRPLNk3OjApx22fLf9pggb+wtuIXjEFMz0Knuqv62C9fe6\nCmY6PuX3enOFftda0coOGx1PyVpYr0Ywz1FwjQHMV/rBW4It5+Ap5UXBL+OTifJbl2DrSBaz\nbybv9kaW+d6sGlnV+rUOVaqqkum7akeLYKZ3k/SPenaT1APXYlXREPwyr8Wi8alr5HPw5BeC\nfQkiuGG/J/b3OEWa+N7WG8CZCO4i95Pd+SPSxPeKAziTy4Wd4Dv7etIwe5bJKm3NghNhuLmw\nb/Eo6p84guF3fww3wuxJtLFiTHyv5U7grDDd3I/yycVRriZpAQzBe2G4ee3B4m6CEcApMBtZ\nx/v7mayAAE6B4ebIWPPpd9/KMhM7egDXKDjRVzbd/F0bxi5ewxxugtvaAzgTwR0/V8aO39N+\n06Edbsa0KdFFcO0BnJHgLoz1sehGbfpdY381kzkzvxC8H7YI7rpJxy+1OxfamMJXMjGpPoBz\nESzPwR/6OXgmtOlfGrOmnwtGAGciWNz9/TVvRVsFq1PwP4GbYOUXgnfE7AdfFv1gXXCj/eMe\njSwEsNMkf2b8NA9vXKlfVWAcsHSBNZnm6fhp7iPVB7DbbTobf0dtUp71Xa6CLeiCG7175C54\nHsAVCna80S65YLXvV0WfxUKrNftdFazPsjOm1WmT6tiwnQR73J70umAt6ZlgsalT8MHOcNSc\nJ6tNjF0c4noEe9ye5NCFHWtlraFljG5tZ2L4rUrwwEoEzxpVc5uzJD4F9TzkZ2+189pMLq9M\nKg7gkfUquv8xbHTB8/uU9AhWH3C6PSm+YAQw32hFLypbPYL5tmDbzvIXvFNsp0wgWLDehLaf\ncMsRDL8DtvAdNmMzya2RxXWnIRpZz9nIRHaRIJjbv7jR0VnvJhldIWb7tJ3ogsVG8wvBDgSR\nEjKv9UwQwAqigsUGASxw/uZhH18YV/DsqWYSCN6ZHQTDr4SiYARwemIKXrSwIHh/ogoWG/hN\nS0TBCOAciClYbGZ+IXh/4glGAGdBRMFigwBW3J6/JQbRBC+7SBCcgriC4XeCmGAEsAktwZYW\nFgQnIZZguUUAa5ASjABeQkuw3ELwhHqW/d5EEYwANtGeZr8zcQTLLQQrpodh705owSJ2rX4h\neCFYu60hHhEE9xU0Alhws9MfLFawfIEA1lipouez2eMQVrB8OJCthQXBJARrj39CAM9YaWNt\nPzAlCCEF93pbYRgBbGBvQhcp+AzBVmxdpLIEqwA+w68Fax84ut+QgtspgiF4SfmCpyp6+V0g\nmILgMYRbBLAFWxsrvt+w3SRpuF36hWBO5HJh22H7KhBMRLAcqoTfnIBg4oQVvHLNE4LTEVTw\nyogr/CYkpOC1eQsQnBAIJk5AweaEhRH4TckOEQzBKYFg4sRvRcNvUsL3g82DEJyU2IuRQnBi\noguG37RAMHEgmDixBcNvYiCYOBBMnMiC4Tc1EEwcCCZOXMHwmxwIJg4EEyeqYPhNDwQTB4KJ\nE1Mw/GYABBMHgokTUTD85gAEEweCiRNPMPxmAQQTB4KJE00w/OYBBBMHgokTSzD8ZgIEEweC\niRNJMPzmAgQTx0Fw06HvLtMgOF+eC27Upt9tlmmLTOA3G3wFNxyCi8JTcMMhuCzeFfxPYGYC\nv/ngJ7jhThEMwfngJVi9hOBi8BPc80ww/GaEdzfJIYIhOCMgmDgeI1laQ2tzJAt+cyLCWDQE\n5wQEEye8YPjNitCCDxCcFxBMnLCCDx0hMgTBCCr40BMiSxCIkIIPBxjODggmTkDBBwjOEEQw\ncSCYOGhFEwf9YOLEf/IZSAoEEweCiQPBxIFg4kAwcSCYOBBMHAgmDgQTB4KJA8HEgWDiQDBx\nwgh+lX8vfzKXrEP8+aKStoT/kHVsIDirrMMDwVllHR4Izirr8OTfSgBvAcHEgWDiQDBxIJg4\nKQUbK42Hzj1WvlFLHZyEgs0V1gLnHinnuKUOD1XBTaycIdiPeBpiCi4ICPbOGedgDyKGWays\nI+YdA5KCY54ncQ72IVqUNdE6MxDsQdS/EyJYknSgo8TcIdiZiPWozD5WvmhFg4yAYOJAMHEg\nmDgQTBwIJg4EEweCiQPBxIFg4kAwcSCYOBBMnLIFX9gv57/slLocGVO24Ac7cn4SlsEKZQvm\nn+z+za6pS5EzhQsu7fL7/pQu+Jux79RlyBoIJk7pgpvjEVX0FoUL7hpZd/aZuhQ5U7Zg2U06\nskfqcmRM2YKHgY5L6nJkTNmCwVMgmDgQTBwIJg4EEweCiQPBxIFg4kAwcf4DnjTA5ur+TSEA\nAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fitted  %>% \n",
    "    gather(variable, value, y, obs, overall.pred, subjectPred)   %>% \n",
    "    ggplot(aes(x, value, color=variable)) +\n",
    "    stat_summary(fun.data=mean_se) +\n",
    "    stat_summary(fun.data=mean_se, geom='line') +\n",
    "    theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Surprise! Fitting the model to each subject makes very different predictions and it overpredicts probabilities at the higher values of X where the underlying logistic model levels off. Why does that happen? To understand it, let's plot the fit to several subjects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAHgCAMAAAC7G6qeAAAAPFBMVEUAAAAAv8QaGhozMzNN\nTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD4dm3///+sNJpQAAAA\nCXBIWXMAABJ0AAASdAHeZh94AAAY2UlEQVR4nO2di3ajOgxFydykr5nOtPD//3oDSRMexk/J\nlsTZa00mTTEIedcxYEw3AGCIrnUAAFACoYEpIDQwBYQGpoDQwBQQGpgCQgNTQGhgikKhO+Fo\niVNfoFKB0DJQF6hUioWmCYOHuSct4wiiMFCpQGgJKAxUKhBaAgoDlQqEloDCQKUCoSWgMFCp\nQGgJKAxUKhBaAgoDlQqEloDCQKUCoSWgMFCpSBL61xXC1fF5gkDlIkjoX48XKpg8QaCCgdDp\nIFDBSBOaFIWeqAlUKqKEVtI1hdCCkST0mHoNnqg5KJwSCqGTiqMPTQHjXx5a6LTiEJoCzvPQ\nEDqpOISmQGGgUoHQ6SBQwQgSWs2xFgIVjCShyVE4REJNoFKB0BJQGKhUILQEFAYqFQgtAYWB\nSgVCS0BhoFKB0BJQGKhUILQEFAYqFQgtAYWBSqVYaNloiVNfoFKB0DJQF6hUSoWmiQIAIiA0\nMAWEBqaA0MAUEBqYAkIDU0BoYAoIDUwBoYEpIDQwBYQGpoDQwBQQGpgCQgNTxBh5vjJ7+/gh\nsjgA9Ygw8vx4mX8QXRyAimQIPXcbQgNZQGhginShf979NwKhgSzyhY4tDkBFIDQwRbLQC58h\nNBAGhAamgNDAFClXCjdXWCA0kAbGcgBTQGhgCggNTHEsofu+dQSAGQgNTHEsoQcIbZ2DCQ2j\nrQOhgSmOJjSMNg6EBqY4nNAw2jYQGpgCQgNTHE9oGG0aCA1McUChYbRlIDQwxRGFhtGGgdDA\nFBAamOKQQsNou0BoYIpjCg2jzZI4g/+wnMAfQgNhJM6ctJ5rRqvQMNoqEBqYIkfotOIy77Qm\nDUrQLjKFImgPAyQL/ehDx87gLzQXlFEJ2kUIHV5kIfR5SO1y9DK/3kmFFrSLPOrJ2b8Q3H3o\nMb8ys0EYlSChWULR0z6zC93PXoVBF1RPurYyONoPKfsWBa/Q/ep/UVAFJemPtn+8EK9TDZxC\nP7+pROYEQketUcaORZM6g3/ClcJ+570UiGKS9C3UL/4jXKMe+MZy9J6fZEASU+941wzyPy4B\n+5QIl9DrbyqJmYHQodU91nM60ayRHyahtxkVUN0bCGIS1a8i/uN6ruR0dKEd+Wxe2w4gdHht\np5Mimwceod0Hxs2re0t5SKIOFEj/uK6VqM3lCQah93Jp0Oh1+ba72O/+kI5Glyfohd7NJITm\nherbQmXD/IBaaN95eHNGyzr0LTyvNEp8d1nbxZQ5xEL7MyEvT0URuQo33MXCr4tnuyyvmhKg\nFTqQCoGZKglJltCbLSeEsuhkCKylBCiFDn9TyctVQUTuos12Mbf/s+oxa+5ujBAKHZMJcdnK\nD0jYyZy8r4vN4Z+4CkqFTuioVMjLV3ZE6oV2ncyQVz+pUAkd+00lLmO5Ae2Xa7KLqf0f55k5\n7d2NERKh+4RMiMtZXkC+Ui12cWebzo/d55lTKlEwRELHlxCXNRNCx/V/rh7vXzSx4XODue3E\n5S0nIGnn2+Muz3qvAIqrlzwg9GGE9l7PNtI+N5l9VFrqMuIJFam9i+HxBoHhGdLqJJ8W0+lK\ny15yPOEClXcx0KMPDjaSViMFQGgDQnu31oeHzkmrkBKaTHguLYGJ8Yi7JOrZ2LVxDvaPpFVH\nETkTnpdPp9sqhTuNVe/75abkNnjnad3ANiN+GctuNu9dDWlHsKwkTjSznk9XndDiINirnWw+\n1u3vkRAEIIk2QktLY0o4scvW20XXlhZ/KuGTIHbgn/DcjaxEslzprLWLri5Q5Bg6W93nieT5\noZMnPHcjLJMcmjYRerq6HRuKsEogodkzVmQlMzoajs5JGevhGs5euTMUWVVARGofevW+4Kyf\nrHSyjH6tsovrBjpiqdsnsiqACgh9Q63Q8214zplsQpGVfjoaPtZNVkpZbripsIvPTSSN1pCV\nfEIyhKZ6kqysnGoXOm24hqzcU8I44XkQUVnluZ7NvosPn6OXHMx2nyeaPrxeVF55Rhxx72LE\nVfvlooOwtFMDoX9QKfS0+sjr5/3qf5s0FVpWbnlG7fPuYp8yvqmfvZqlrdCishuIJTdUzl3s\n04br9ba7zxOthRaUX28s+Tf5c+5i4mA9IzMVeIHQD3o/+aulDHJBcOz+CknZ5qKx0KL6HFww\n7WP4XpQj0lroYxjNsJPTxOT0q1UPhK4C9V5OnedjpC6R5kIfpFoI9/Jx49YxMpcIhK4EXbfD\nxJMj2Ggv9GEqhmg/4+59PSwQuh4kO/o483yYtKUBoStS3u0w82wfNgQIfaSqKdzVxdN9ylZl\nFQhdl6J9hc9hJAh9qMrJ73YsB24cKWcpQOjqZO7tchzSsVKWgAihD1Y9OXtr7nmCXEDoBqR3\nO9bDRA+WsAQgdBPSdtjiAzK5kCH08SooZYd9s0+DFRC6EdHdDu9s6mBN6gz+A91EMwsOWEVx\nuxw/9yIYSZ/bjm6imTlHrKKIfU6YSxRMJAt95mmhD1lJwW7H3lOMySOxQ8YjKbZC5z8X5+Qb\nPEaxWuHsKH3/eOcRR/DZQ77Qsxn8mz51xw1HrhjYFfrK6UR9A/oRaDj76ApU0xItf5PCSBN6\nM/c5hGZB0XeMNBKFPtM8eNMJjH4AnbPJeCQFUwsNoQEBgoSG0aCc1Bn8BwgNJCNkLMcNGA1K\ngdDAFBAamEKU0DAalAKhgSlkCQ2jQSEQGphCmNAwGpQBoYEpIDQwhTShYTQoAkIDU4gTGkaD\nEiA0MIU8oWE0KABCA1NAaGAKgULDaJAPhAamkCg0jAbZQGhgCpFCw2iQC4QGpkidwX81mz+T\n0DAaZJIx++jANLfdHAgN8oDQwBQZc9vVEBpGgzzyhZ7N4M8AhAZZpD80qMpBIYwGeUjtckBo\nS3Sd+6eO3h+xQsNoQ8gVutZZjgFC2wVCAw18d5fp/0v3b/j72nXn92FU9t/55Sbu7LPhtXv5\nGu5Cf7913ds3VRSpM/hXulI4AqOV8dqNkn5dvf7sJt5HZV+6t0nc+WdXtbvz913o8/jxhSoI\nmWM5JiC0Mj5HW4f37vPaSP8Zhn+jrpPBk7jzz16+h5eb2sPwMb55734TBSFYaBitjcvtO3x8\n+/X58XKT99G1WH42NuR30cffda9EMUBoQMbv7u/wt/u4vnu59S9+VJ5eHZ/d/v18ToJooX1G\n5/8ScPF97S6/d9e+8Vt3+f35tZTX9dnhhM6HMy6wx1v3NfUdbmcvtvL+fLbtcpAhWWigjb/X\nlvbvMIr6d/h+WQv9/OxlfPdx+/x9PCj8c/2IBggNCLnczr+9d9v+8vyz+Wm77+m0XfePKAII\nDQj5PZ6bG8a+R/fyd3VQOP/stXt9nP34mn5BFQGEBqaA0MAUEBqYAkIDU2yM/H09Ah1eqI45\nAajLSujvy3Rm5XY2EQB1rIR+697Hcyl057kBqMpK6OcF9jbhAFAGhAamcHc53ru3NuEAUMb6\noPB2Yb07f7UJB4AyNl2Lj0vXXd7J7lkEoCroKwNTlF4pFI6WOPUFKpXNWY60wNslNg4tceoL\nVCrFQtOHRAfrFD2UKAxUKs4Av14+YouL3kGFnqgJVCruAL+7SKNl76BCT9QEWrQWipWkrRtd\njqooDLRoLRQrSVr3n+7s/HxbXHT+FXqiJtCitVCsJG7dj2PC98jiovOv0BM1gYZwTo1yP9sw\njk/uZj9T4hb6HOmz8Pwr9ERNoG780/5095e71Y+fSSMsLC46/wo9URNoVun7y+J/YiC0BBQG\nmlX6/vJomRn2dr7GjCtCsvOv0BM1gWaVvr/Muhqsfei2Qv+6vV4hWyWHJz8BHj3QrNL3l0Xf\n2Wof+p79xwsFDJ78BHj4QLNK31+OcFD4a4DQy3eFSBS6wWm7xxyRrbocy3eFcHVNf63+L4Yx\nULruxojwLv6wEfq9dR96+a4QCD199YnpQ9dgFeC5+/fSfX2/xE40Qy+09GOtWYAqAr29UK1T\nndDXgD+6z+E7dqKZQ7fQwoUe5PWhK7AV+nN8ZBy6HD5+zV5JgNBkrAJ87f6MT3P520poNScP\naI+2jnKWowKrAEeTp+fJRU40c1yhSc8eQGgy1gF+XqaHYTQYbafrSuGvX5SRHuRKYQ1WAabO\nCy17BxUOkVATqFTWB4WXz7TiondQoSdqApXKKsBL150/EuYBk72DCj1RE6hU1gF+vZ+77jV6\n/n7ZO6jQEzWBSsUR4N/3rrv8iSwuegcVeqImUKk4A/zCAP+6KAxUKq4W+u3aQv+OLC56BxV6\noiZQqTj70G8JfWjZaIlTX6BS2Z7luPxOOsshGy1x6gtUKuvz0K+J56EJQwFHo1v9v/61509o\nV7zVL1IfRQGhQT4ee7rnr51LxQqdCoQG+UBooIbTyfHhvQcxWdMNj5tlu+ev7v9Pr0+hu+Hx\n4W0RCA2qcHJz/+1T5Pv/3ePnza+GeQu9+jPwnT6E0KAai4PApb2Lj4Zno71s0heL+rdRGCLI\nwznprF3uDatL6Lu7C6Gf5W6FZ4tAaJkcTGhH52HeQg9+oV0/bDdQGF9ZcXA0ofc6zBDaBsfy\neX1YF3dQOMwdxkGhbI4l9PrE2/5pu9Wpuc5VemcThRGWFT88BxN6B0qLIHRTIPQIhLYCfB58\n/eGstUUsc74ye3uePcMQQhcBocmJMPL8eJl/EF0c7AOhyckQeu42hC4CQpMDoRsCn+lJF/rn\n3X8jELoECE1PvtCxxcEuEJoeCN2Q3j0MHhSQLPTCZwhdQr9zXwcoAEK3YxS6dQzmgNDtgNAM\npFwp3FxhgdBF9PCZHozlaAYaaA4gdDPQQHMAoZvRw2cGIHQzIDQHELoV8JkFCN0KCM0ChG4F\nfGYBQjfihJFJLEDoRqCB5gFCtwENNBMQug0QmgkI3QT4zAWEbsEJN6twAaFbAKHZgNANOOF2\nQjYgdAPQQPMBoeuDBpoRCF2d8ZIKhOYCQlcHQnMCoWsDn1mB0JWZxnBAaDYgdGUgNC8Qui63\nQXYQmo3EGfyH5QT+EDoVNNDMJM6ctJ5rBkKngQaaGwhdk/uofgjNR47Qm+JMT6xW8yDs6ED5\nhc7PmZps+0kW+tGHns3gD6HjlqvQQEPo8CILoc+Ds8txdKPjFqvR48hfuZpseyHqQ7PkoteS\n5D5u/3/ui2VtobPXnl9SFFQHhRzJ6LlWTE2cC4/7vKUKrSLXIcjOcjAko+daMTH9EBVljQY6\nX0s1rUcAwUL/rFF8tyNO6EoNNIQOM5/Bf/9KIXkyesc7iUR+kVQTOmsD/eI/xRCO5aBORu98\nK484oZ8zJUkWWnamY5Ar9GJ1krsdUUJX9BlCUxWnTUbv/VEOcSrUFDpjE1q6d2HUCC0201FC\nz6ZmhNCskI6HpkzGdl0yux1xKlRtoNM30u+8VwjtAH/CZLhWJTHXUULXbaDLhBaZ5XhUCS0w\n11Ft2ylmIQJyOw695ydtEN+CRZaMnRWJ63ZIEjr7PKeW4+8YqO8ppErG7nqEZTtG6PoNdJnQ\n0nKchDqhZWU74sv6dGogdNJ2NsuKSnEi5Hd90yTDtxZJ3Y4ooYOLMMRSJLRmo+mnMSBJhn8l\nYvId8119Ci5BRO6hnZbzSXGoFFpMviOEruYzhJ5gmGiGIBnBVcjodkR8V58Cv6cj98jOvZyI\n/ObAMXNSeTIi1iBB6bDQ66cRqhFardEQOp/wd/Xm6ZoVhY7c1t5S7dObB8vcdsXJiFpB85QH\nha7pM7XQ7dObh0ihC6uiFiGht08/rtlAQ2jC4oXJiCzeOOUZR1NVhY7ammcZnUYzTadblozC\n45lKhIR2PJ5eldA6jZYodHThphkPHU019zlme7HfJyfH3siEa8LzkroruyZQi4DQLgNEC712\n9rSkKLqKsM3gX1B5CUXbGb2/5ek3TgFqCx3c4Kx/tCGmvEByZvCPmh86Pxdl42pq4Rfa3aDx\nRZt5daR/mJxXXiCJMyetJ4j2FM9ORlLBZin3Cl3b5wyh47oS+ozmEzo7GSqE9p4dqN/exV3v\nm/RN61GoM5pkBn83mblILNYo477N7rV61XscW6GT7zYwL/TZNYP/DnnJSC3VJOV+n6uP9wkc\noU5suxc6BhgkwvrQoKxkaBf6lHvKoYCQ0O7OspIRBmmk9qFX72/FvQfJ3uOO7S9jz+Y/f+nI\nOPdpU7/P2WeFc/Ff7ts78tMxwCAVIqHFQZmjLfuV7HkQYROhPanQMcAgFd7nFGbkok6RQna3\n6HssUH2hR5tjutd5a5dJhtApj0ZOTkaVbncpIZ+zhwrl4V7zT9scd0Yvff1CoZvB302dQ7zK\nKc8a9FNT6Hmnq1xoVUazjeW4Y1HocAPtWqaO0FeR193m8nMuEHpGWjKqXF0sZbfL6lumUgPt\nOAgkOImoyGh2oStlrmLKY3zeLlWtgfYv4PsschuyqSB0lS+39kIHZuCoInT8CedUodUYLUvo\ngrTVy/jOSYXAUmzxzS9uh5bwf+bbCISeUeUEUbWMx/m8Xq5CA71/ManmEWprzAhdrY6cV9qD\nIyX4hfZcHIXQxMXrXJOqU0murbi/6/vdHyjpfSHsbd2uz3WErjMOpqHQ4SWZhU688QRCFxav\nM7CrRjU5OxzhZXl7HOGxWLX6P+2pI3RcCovzXKGiHJeZY5ZlbaAjhhZCaOriMTlUIHSKz/Ol\nGYWOGyrb77w3hymh+atqs4HIYzE+oSNHfkNo6uJZI9SS4a6r1foDrWO/eUNN9J0MEJq8eDCN\nCoRe+xy7PFdYp5xT/KZ9tiY0c20t1x5uHXmF9t6PshcKhCYrHkgkUZ45qyvV50cJlqB2by/3\nhmLc55pCB1KpTOiUkwtcPkPoLeaErjNOM6p5fpRhCCljdlDGPy5B1BTam0yyRFcROnqaBCah\n/fe/ekKB0KTFPdmkSzT3OKCkeWwYhH7eAQuhtxgUmqvOnj4nliKOJ3xD934o5n1OnvB8SJuX\nY81uPikTzSx02qxMDELPVp0aCoQeNrMlpc3LsWYvoaSJZqm1PJ/HcsQ9jtmaG4cikWShz0Ut\ndB2hWertts70WfMgdFUyZvAvEnonpcSJpr+n877GjEkgiWMp8FnRra755AsdnvDciWqhcyY1\npY3FOzlTCAg9UTJZowNnUsVnOttnWhIfKHFA0oTeTBWdc9bPVRHiK6cX8TTVRQTic9aERKHP\n0c8p3MdREeLrpmefQz0G+BwmYwb/0hbaURXSKyf2RidmIHSYFkJv60J65fTyfBafs0akTng+\ncAgtvW5E6Ayfo6g7luOH3vujMGToHJrdFNxoI/SqOiRXTtJ9ToxkPTXzgEgQWnLdJN7nxAZ8\njqSR0IsKkVs5p9T7nNiA0JG0EnpeI1IrJ+M2Jy4iH5wCJAgttG4Cj/mrCnyOppnQzzoRWTmP\nUxsSogs9HgA8aSe05Pvqs+9yYiH4eADwpL3QAuumZMwxPfA5hYZCS70N+VQ05pgeCJ1CS6HZ\nZmEp4SRtyHHsozTBRHOhZdXNSdyANvicRlOhGe7xL+C0GfMsILbw87XAgrZCi7oPeTMISUBo\n8c86BjdaC123bjyb246pa6jNT5zu9hk+e4DQE+FHwValv3PqnbQLTD6Nha6OUwbniGcB2sgY\niK2Lowm9baT3bn9tLzR8zuBwQq9MjXkMbCPgcw4HFHrm6v7dVfBZKUcU+sdW382CzYWGz3kc\nUuipI+2997W5zxA6k2MKvXskeKe1zyLmadJJ6gz+q9n8dQp9Ct362l7oxgHoJWP20aFwbrvW\nhO+tau0zyOdYQl9dnt2MUudpAqAqGXPbqRY6ZngohFZMvtCZM/iLQufk68BD+kODDBwUPnF0\nO+Czao7V5XCgbmpf4OXwQuuaCRWEONZZDid6ZkIFYSD0siMNn5WTOoO/iSuFG2RPSwYSOOhY\njjWCpyUDSUDoG0JncQKpQOg7vag5QkAuEPoBbqe2AIR+AJ8tAKGBKSA0MAWEBqaA0MAUEBqY\nAkIDU0BoYIpSoSP4L2YhnpJJcaoJtGGc8huwChH+p6WklkC1xNkECN1sc9kltcTZBAjdbHPZ\nJbXE2QT5nSIAEoDQwBQQGpgCQgNTQGhgCggNTMEu9Greg8TCFbepJVAtcTaCW+jNPGJJhfMK\nZm1TS6Ba4myFZKHPmQWrC10zUC1xtqJKHzo7jQX5r1iyfqBa4myATaFzu3xaAtUSZwNqCJ0v\nV2bJ3KJaAtUSZwvECl3Qb8suqiVQLXG2oILQua3COfdkUVVPGgSqJc4m8AtdkoiqDYqWQLXE\n2Qb+Cyv1C9f3pPr5xVwgdCn533NT6fyN5pTREKiWOFuBsRzAFBAamAJCA1NAaGAKCA1MAaGB\nKSA0MAWEBqaA0MAUEBqYAkIDU0BoYAoIDUwBoYEpIHQyr92/YfjXvbSOA7iA0Ml8d5dheBmt\nBvKA0Ol8dJ9/uvfWUQAnEDoDRTdwHA4IncGfrvvTOgbgBkJnAKHlAqEzOF8u6HIIBUKncz0o\n/Ow+WkcBnEDoZKbTdpfuu3UcwAWETuZ+YeW1dRzABYQGpoDQwBQQGpgCQgNTQGhgCggNTAGh\ngSkgNDAFhAam+B86q+LPqgFCzQAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=4)\n",
    "fitted  %>% \n",
    "    filter(subject %in% c(8, 6, 13, 32, 36, 9))  %>% \n",
    "    gather(variable, value, obs, subjectPred)   %>% \n",
    "    ggplot(aes(x, value, color=variable)) +\n",
    "    geom_point(size=0.2) +\n",
    "    geom_line() +\n",
    "    theme_classic() +\n",
    "    facet_wrap(~subject)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you see it? These subjects represent two classes - those that show generally increasing or decreasing *value* with *x*. We have both due to the random sampling, which introduces noise. The logistic growth model, however, can only fit an increasing function, and it does that for subjects 9, 32 and 36. Subject 36 and 9 show similar to the original pattern; however, Subject 32 shows a much more extreme slope. And while the model can overfit subject like 32, it fails with subjects like 6, 8 and 13 and is in essence biased. For those three subjects, the best thing it can do to minimize the error is to fit a flat function. Indeed, if we calculate for each subject the error at each value of x, we cna see clearly that the error is biased at the extremes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAAM1BMVEUAAAAzMzNNTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////p0TvgAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAGu0lEQVR4nO3dAVPiOBxA8VBw8dxT+f6f9gBPSQGBNin95/X9Zk4ZbieT\n9AmWwG7TTmhp7gloWgaGMzCcgeEMDGdgOAPDGRiuXmB/VEIyMJyB4R7N0u1dud2d/oSBQ3ow\nS/fzpXe7M3B0RYE7H8HhlQTufIqOr0rg1YGBQyoI3O18BMc3PnB+14CR9FwFgb8MHUnPVfYy\nyUdweAaGG7qT1WW3d3cCmzyAKfeiDRyAgeGmC5z2qg2usSYLnL5UG17jTBU4JQuHYGC4iQIn\nAwfhIxjOwHCeRcP5OhjOnSw4A8P5wXc4A8MZGM7AcAaGMzCcgeEMDGdgOAPDGRjOwHAGhjMw\nnIHhDAxnYDgDwxkYzsBwBoYzMJyB4QwMZ2A4A8MZGM7AcAaGMzCcgeEMDGdgOAPDGRjOwHAG\nhjMwnIHhDAxnYDgDwxkYzsBwBoYzMJyB4QwMZ2A4A8MZGM7AcAaGezRLdln30+Xes/sMHNSD\nWbqfL6fb+X2Pj6TnMjBcQeBd/t3AQVUJvDowcEiFgT3Jis7AcGWBs74GjqkocN7XwDGVBO71\nNXBMQ3eyutPtruttZRk4JPei4QwMZ2A4A8MZGM7AcAaGMzCcgeEMDGdgOAPDGRjOwHAGhjMw\nnIHhDAxnYDgDwxkYzsBwBoYzMJyB4QwMZ2A4A8O1ENgfnQIGhjMwnIHhDAxnYDgDwxkYLn7g\ntDfR0EsQPnD6Ms3gCxA9cEoWLnJ24DZ/ao1UiYELnR24bvyBnCRBMnChswP3vtl+1BmpkqqB\nl/hTcrbmgkeMgUOKHrjqWbSBY4x0Nm69X8AGjjHSZCMbeLf73K5TWm8/y0eqx8AFztb80X39\nyuuGn0sbOKSzNf9Jm33aj00avuFh4JAuzqL738ePVJGBCywo8DL3w5bzFL3QLc8WTrKqWOqm\ndgsvk6owcJyRprDYt6Wivx9cjYGPor0fXI+Bj8K9H1zPMvvGf7uwniXmXVTgJqZY3VLOoo8a\nmGJ1izmLPmhgitUt5iz6oIEpVrecs+hdE1OszpMsOAPDPbrm/GruP5d7z6/w3sLRa2CK1T24\n5u7ny+l2ft/jI82pgSlWd7Hmf172T8+b97N7GYGX6CzL5/r4+zelf/v3G7hVFx/Z2R4+j/WW\nNv37bwdeHSTF8Wvgw//6/u/xwNdGUgwGhrv+FL09/1SlgXMtLfX8JOuXT1UaONfSUi/m+nr1\nU5UGzrW01KE7WV1+u7WdrFpaWuqi3vCvpaWlGniElpZq4BFaWqqBR2hpqQYeoaWlGniElpZq\n4MH6m/nRGXioi/drYjPwQFfekQvNwAMZmO3ae+qhGXggA8MZmK6tvgYerqW8Bh6lpaUaeISW\nlmrgEVpaqoFHaGmpBh6hpaUaeISWlmrgEVpaqoFHaGmpBoYzMJyB4QwMZ2A4A8MZGM7AcAaG\nMzCcgeEMPKNnHDIDz8jAcAaGMzCcgeEMDGdgOAPDGRjtOX+JzcBzedJfQzXwTJ71F8kNPBMD\nsz3tH3Mx8EwMDGdgOs+i6XwdjOdOFpyB4QwMZ2A4A8MZGM7AcAZWMQPDGRjOwHAGhjMw3KNZ\n8su5L/4S7y15MEv38+V0O7/v8ZH0XAaGKwi8y78bOKgqgVcHBg6pMLAnWdEZGO5ulq/XQr8E\nzvoaOKaiR3De18AxlQTu9TVwTEN3srrT7a7rbWUZOCT3ouEMDGdgOAPDGRjOwHAGhjMwnIHh\nDAxnYDgDwxkYzsBwBoYzMJyB4QwMcOvQGxjAwHAGhjMwnIHhDAxnYLTb/3K8gVt359oPBm7c\nvau3GLhxBma7ewU1AzfOwHAGpvMsms7XwXjuZMEZGM7AcAaGMzCcgRfMwHAGhjMwnIHhDAxn\nYDgDwxkYzsBwBoYzMJyB4QwMVzHwZFbTDT1OtAldzmeKwNNZzT2Bc9EmdGs+Bh4h2oQMXFm0\nCRm4smgTaj2wChgYzsBwBoYzMFz8wL2LFAcRakK3D1D4wJfXkp9fqJ+4OwfIwMN1oebTeuCj\nSAf0MJlI87kzFwMPFi5w27+DDyIdz+NkIk2ouz0hAw8U7pwA8Ds40uE8vigJ9cKt/cCBDua3\nSFNqPnCkg/kt0pxaDxztGfEo1Hwa38lSGQPDGRjOwHAGhjMwnIHhDAxnYDgDwxkYzsBwBoYz\ncOYlve9272kz9zxqMnDmM613u82hMoeBc6/p71vazj2LqgzcE/DDBYUM3POW0tvcc6jLwD0G\nhuvWa5+iwfYnWX/T69yzqMrAmePLpHX6nHseNRk48/9Gx8vc86jJwHAGhjMwnIHhDAxnYDgD\nwxkYzsBw/wE80awkS52t/AAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=3)\n",
    "fitted  %>% \n",
    "    mutate(error = obs - subjectPred)  %>% \n",
    "    ggplot(aes(x, error)) +\n",
    "    stat_summary(fun.data=mean_se) +\n",
    "    geom_hline(yintercept=0) +\n",
    "    theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The end result of this is that when we average the predictions, the model ends up predicting a much steeper almost linear function, and fails to recover the original data generating function. It has overfit the random noise in a biased way due to its ability to only account for increasing patterns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Would having more observations per subject help?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I've made the point that one cause of this problem is the relatively low number of observations per subject. Let's test this idea and see if the bias will disappear if we had more observations. Here's a function for the entire process, and we'll vary the number of observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modelFunction <- function(nObs, nSubj, sds, a=1, b=1.5, c=1) {\n",
    "    # Prep the data. Generate parameters from each subjects from a normal distribution with sd=sd\n",
    "    data <- expand.grid(x=x, subject=1:nSubj, nObs = nObs)\n",
    "    data$a <- rnorm(n=nrow(data), mean=a, sd=sds[1])\n",
    "    data$b <- rnorm(n=nrow(data), mean=b, sd=sds[2])\n",
    "    data$c <- rnorm(n=nrow(data), mean=c, sd=sds[3])\n",
    "    data$y <- data$a/(data$b+exp(-data$c*data$x))\n",
    "    \n",
    "    # generate data for each subject\n",
    "    data <- data  %>% \n",
    "        group_by(x, subject)  %>% \n",
    "        mutate(obs = draw_samples(y, nObs))\n",
    "    \n",
    "    # get overall fit\n",
    "    overall.fit <- optim(c(a,b,c), errorFunction, obs=data$obs, x=data$x)\n",
    "    data$overallPred <- myPredict(overall.fit$par, data$x)\n",
    "\n",
    "    # Fit parameters of the logistic function for every subject and extract them\n",
    "    fittedPars <- data  %>% \n",
    "        group_by(subject)  %>% \n",
    "        nest(-subject)  %>%\n",
    "        mutate(fit = map(data, ~ try(optim(par=c(1,1.5,1), fn = errorFunction, obs = .$obs, x=.$x))),\n",
    "               pars = map(fit, ~ .$par))  %>% \n",
    "        select(-fit)\n",
    "\n",
    "    # Get final predictions for each subject\n",
    "    fitted <- fittedPars  %>% \n",
    "        mutate(subjectPred = map2(data, pars, ~ myPredict(.y, .x$x)))  %>% \n",
    "        unnest(data, subjectPred)\n",
    "\n",
    "    # Calculate the bias\n",
    "    bias <- fitted  %>%\n",
    "        group_by(x)  %>% \n",
    "        summarise(subjError = mean(obs - subjectPred),\n",
    "                  overallError = mean(obs - overallPred))\n",
    "    return(bias)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate the error for different number of observations\n",
    "set.seed(19882)\n",
    "bias <- data.frame(nObs = seq(15,105,15))\n",
    "bias <- bias  %>% \n",
    "    group_by(nObs)  %>% \n",
    "    do({modelFunction(.$nObs, 2000, c(0,0,0))})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAASFBMVEUAAAAAtusAwJQzMzNN\nTU1TtABoaGh8fHyMjIyampqliv+np6eysrK9vb3EmgDHx8fQ0NDZ2dnh4eHp6enw8PD4dm37\nYdf///+QBu+fAAAACXBIWXMAABJ0AAASdAHeZh94AAAR/ElEQVR4nO2di2KjKBhG0Y1Nx7aZ\ndtrI+7/piqCiooJy8+93dmeamhQcT3/kJjAOSMNSnwAICwQTB4KJA8HEgWDiQDBxIJg4EEwc\nr4Lx25IfEEwcCCYOBBMHgokDwcSBYOJAMHEgmDgQTJxwgp/Pp8+0wTGCCX4+YTgHQgl+PmE4\nCyCYOBBMHNyDiYNaNHGCtoNhOD0QTJywPVkwnBwIJg4EEyfwYAMMpwaCiRN6uBCGEwPBxAk+\n4A/DaYFg4kAwccLPyYLhpEAwcSLMqoThlEAwcWLMi4bhhEAwcRwE31rmr/Vj64nBcDrsBd+G\nv8bX+rGNxCA4HRBMnHOCuf51IzEYToY3wf+1QHAI/t4Zu390L5l7nfi8YItKFgyfoNUrePnm\nEEyRO7t/cf7VfuFpBGt+txKDYRcY+35lt7f21WcnlgvRn53gV3YXsczfb+zlwyYp60zNgnW/\nEOwLxm6iWG4Nv7IveeiTvYrjr+3h2w/nb125bWH4pOCJXwj2RVur+uEf7KYXyuJVd/wuxLcx\nzr/YbSsR9WP2uQ69V+Pr223SlbWVGAw7IPT1Ssdj6vg3e2mvPPvzaZeU1/PaeA+CHZBaTYL7\nV583VbHeTcrreW29CcP2jCKHe/CXvAeP7/57Ybcvi6S8ntfWmxBsjxapYy36r15ECz5sWk0R\n18mCYWu0SB3awS/dkTv/ubN3cQ/+4v88V7LOJgbB1miCv18mPVnTZtK7RVJez2vzXQi2Rr/X\n8o9JX/Qre+3qVm83drPwG3UpQxhOAAQTJ+pipDAcHwgmTtzlhGE4OhBMHAgmTuQV32E4NhBM\nnMiCYTg2EEyc2IJhODIQTBwIJk50wTAcFwgmTnzBMBwVCCZOAsEw7AN5rbuZWRaf85rpLhDs\ngnl7IqV1/4pDcOaYNxhjPGfBMGzP6haBbPhrGwjOkqcZ7RNqYu3eLTiNYBi2xiKC86tkQbAD\na5t8MsMrE2kEw7A9K5t8QjBx0hfRj8dj46MwfJLBbapK1uOxaRiCTzL2ZNl8zmumHY/HtmEI\njkQqwTAcCQgmTqp7MAxHImgtGoLTE7QdvCkYhqMQtqMDIZycwD1Z6OtITULBMByD0H3RCOHE\nBB9sgOG0QDBxwg8XwnBSIJg4EQb80SGdEgfBps0prXZdgeGE2Ate2XUFgvPmpOCb1b5JMOyb\nfiaHzxkdZzbGwsyO45gG5foJWV4n3e0I3tycEiF8GOO4egrB24khhA+yMTMGgi/Mw8zkM1kJ\nhuGDrEWwrFxlJHjdMARvsj63La8IRggfZXVqm2/Bps0pub1ghLBHgtSibfNdBYa9AcHUCdGT\nZZOt/k1VVbO31wxDcECCCa6qhWGEcAJCCa4qB8MQHI6ogmE4PjPB9z+eEoPgTJgJvp2K6J17\nMIfh+MyE/ru/fXtKTPhdGIbg2MwEswEPifGl4DXDEByKsIJhODmhe7JsC2kIDkTwrkoYTsvc\nyc/bC2Mvbz9eEuuwNHxOcF3Xp36eLjMn3zd5B74dqktPEivLUr6wbAyfMVzXMLzCTPAfdm/V\nft/ZoQ4PPbGydDR8QnBdw/Aai1r09OvhxMpyNBw8hCF4nRiCg4dwDcOrhCqiJ4LtDB8XXNZ1\nlxkcLwlWydL9Lgx7DuFWr1QsvsLxhHDNJOHXzfAZwdpvU93JPprWZRim7KSck6Vd5oCG6/LZ\nGtV/uO4DmoBn06Ac72tJ+/pCjQcrVg17FNz6NWxNosrqy3s2D7tyZjOj0vAJf+PBivCFdL26\n94x2PzZ7bprmQIZRMU+c6KfMWtgKNx6sWDPsRXBrtX6Whm2jemZ1rqnmpsnWcGVmeL8X7Dxt\n1vdwIXczbC14UFp36a/5FSzr1SqcmyZjwz3mCGbaH7eJ7x4FF0UhX3g1PA/WuuwPb/ysse3U\nXEKw8R489ZqoFl0Ue4YdBRuLYeVX/jZt/XYMjptBbNM8HvkLNtWipzHoJNhbLboo2kt6xPCy\nqF2/w+p+dww384AVT2W2isvsBa9xrIj2Votu/Yr/5TdrjaWl4YnKDbOS0W9rzmh45lXFsXrk\nVsTvg1/V8OQ+vP0xDW+16LoL4P7WZxvCq5upGun98kLoLYqxtrmI18mJ9dl2H20NX1PxsYfP\nvFWy6rpp/yv6by0NW5pVWYztnS6Em2LwutLKkEL7OK6q9oU4cN1i2oJggtvLvG/4nODhVaHF\n61PekIvl52W5LDV3aruvFReG6SoOVYvuLnjd8OFKOxi2y2zwWxWT8lj9eDHxbFjoohtBbhV3\nJ0E3iIMKboQEZ8N2eY034IJP77eTBIq1eO6nCYhGSGeYqGJdsFYuny6iG0vD21srrTP65aIe\nPKlP6YZV5Jo899NAlOHLVqe3WQjulwY4m5iKqs7wVnP4mOGxgcQNYpThZbk81Sz9iijuzo5m\nMR1MsAwqaVgFsT/Bg19ubuQ8N/enLnr6vhjVciJZTIcTLJGltG/DsoJVGcOX72+bKNGLbdly\nolhMhxbcKpCGuwvpybD0qxfPo1El16KqNr0vi5aTMExNcXjBreExiE2GnQXLbgq9eO7X+9Mi\n19nwo+Jiog+1II4geFJMezA8lAiTgYPFmo52hvlYB3xUVVdU0wriGIInxfRpwyIpEb5aqJmX\nZLWfHqIUP7hoFAvFLueTOVPBzFtX5YxG3Ta50bCL4L7d1Sz8Locu7FOVih/dOZV1TchwJMG9\nj4nhQyFcy0R0vatr7rpM8eoUC8MVLynVtcLOi9ZRhotzIVxPi+fe60q7yGkSn1AshyPEKTZE\nHpGIJ7iXcspwLSNNJmXR2nWbptmm/ZCnVXZVQwqOZ5UsvZD+4zz0v/PbMhTThw3XMnxFOlZ9\nGc4TcYtCGq7EKTZdjevibAhmL2cSM1H3xfTSsJWteghfO7sC16nWqkZdienTIqerh/G6k2/3\nitbuD/TV02OG5fSBxjZ4Fc5PS0jFXRDLMidDxV7Wi36fzbA0bU6pH7O5odfNUNfqsTZciuss\n5rru5jLF2fBDKVbFNE8bxqYnq/rZlO6zKtXjo++Gx0dNm3Lox3az6qjHupbCUnAlZvA9miND\nE+6G+67zcmxvp1JsfHbuuOCNB8D9CB4Ml26Ghd/H42Dr9IBhPg/iMYyjhrP5qbnJnHe3ie/9\nEg6vi096Esz7zt5yHMiRXzYEV0VRi3nqNumbOGRYBnGld5qpCQLhDZdm+rfZ4Y2x1C37x3Dr\n3hG8uTnlhFoFRcltDVfFo14Z+rXkiGHVdy5OTss68novKxGs7LoLfmXy7hswgqVhccksDT/a\nxu9Jv0efLReKyy6Ih0OxF/TxfA/mr7KIXvr1KLgbMFCGlWJp2CS41Vu0FejTncOuhtW5dIqr\naiymo6/Y5LEWvT3Y4FFwVyWVhvmm4bbFK/quSh/DOwcNy97zalLXSt4wzl7wUEwPhg2FtOjP\nqApffo8b7nretCDualo+Tug4J4rodfwK7ovpcug56o4+psNDXdezt8mOhw1LxdO6lp9TOoqX\nnqw5ps0pXXuyNKTh7v4yGlYDu9JyJaux/oZmjxsWv4X6nTjL/ksjoR4+s0D2PE4MT+ZmVLIh\n6nPo/YThXvH4FMU1DCcULK+QMixkVvrkqqr36/U6njEsemYq/TmoSxg2Ovm+v/tLbAPNcBfE\n1ShYFs/9hEx/nDJciofML2bY7OSHHTLsHPbS8Dhdqxr9dsca/9fwjOGu4S4Vy+8vYHjFSZQi\nmvdXaBhCrOSaMip82+AOcAVPGm5axeOzjPkbNjv5y27G44cS20Q33AZx1S21ofQG8XvWcLce\nyIUMr1Wy3jwkZoUy3AdxNfrlYfyeNVzKFX2uUkqbBd8O+T02RVMuCiuDWMSvelpINJFDXbxz\ndelSniW/Rk0r4rTZNTTDw+NgTUi/pwyXwrAYYRoU5204A8GD4WYQ3PVhhrxwxw33A7QVHxRn\n3ak1d/LR1q6+2C1OO7hHGW69dkvWhfd73rDsWi1EfYtnHcQzJx9t+6iblxWnHdwzGO57EUL7\nPWm4O1k5PCIV52t45uSFfbV/Pv5FayYphlJaCg7v92RdWjMsFWdreDkn61M80RCro2NArpEx\nrmgW4YI97RflMjEa7hTnanjm5Ma+/7B/4i7sITE3NMMxApg7LasnmE830WfSqHtxhsycvItG\nsAjgaB0dI3UvOKbfU4ZLbSmZXBXPnbyx22cbyBE7OkaU4Uh+3QUvJwXqhsfncrIih3bwQL89\nnZeT2cVd8K7hWFE8DNnvD91nJbg1XItdCH2ciwUBDHP/ileWUpVXe/+K5yW47OakRtunTNo9\nWUpPzrb2rnh1sVyLGZV2n3DgfASn2QfYraY1fzx5+vtYc6+K15dDtgtgCJa4GF6s6LMwvLJE\nqgOFGe0Tw6xZf9Nm97muYPsgNq3JZTDsayGm1Qhm079WyUtwyudCLBUbF10zGfakeKuEnr8y\nkZngpIslWCk2r6o3q2qpr14U7/m9mOC02Cg2r6pXGQ1zMfgZBDb+DcEu2Bk2LOO0YjiU4lHw\npSpZGWBb27I23Cte7bA4hLrUV6tFZ4FtbWt+oJrciPWqhFAs56rEB4IN2CleGp7UpuuZ4iaN\nYQg2Ylehnh+YGp4FcaLdqCF4Bc+GITg7rJpM8wOrhiE4QywUG1rEK1WtRH4heBMbxdNvxfDw\nqmE/J+UGBG+zr9jacBogeI+LG4bgXXaD2GRYU5zWMARbsKfYYJjnYhiCrdhRvGc4oWIItmRb\n8ay5NDecMIgh2BqXIDYYTjSVAYLt2QniyXdq1/jxQKrJSBDsgppHbRZtMDwGcbLphBDsxtbT\nEFuGIfgqHDQMwVdh84Emo2Gp2K/fIOtFW2TrM7FM2X5ibdpcUg+mDYYPZWj6uX42JWZVBmDn\nicR1w8cwRj4Eh2TnkVO/hs337jCCTZtT6js2/BbBggOGTdvjrFKb6d8OIti0KYcwfSixy+Ns\n2LjB1T7rEcyiCL790gh2MCwxb1Fngbn2LSrPcSIYgk3ohvvG0kHB67XvyILtN6ckgvXokuqW\nPirYRMR78K+NYIfRpdGwr6whOArOhv1l7bcnS7aJIHiBm+E0IIJPYD9VK51hCD6Di+Fuw6D4\nHOjJ0jen/OWCXQxXaQyjL/oce4Z7xVWVyDAEn8RyVjwEXxXLJ1sg+LK4GA5+Mksg+DT2hoOf\nigEIPo/r84dRgeDzOD8lHhMI9oDVUg+mBRAjAME+OLzEZXgg2Av7t2HjIrURgGA/7PZ3QPC1\ngWDqWBmOciZTINgXFoZ9ZievtT6zw3z1IdgbJ3Yy3cK8gBobp2MxvnXlIdgfQQybl0BkHILj\nE0Lw6iKmE8EbFx6CPeLPcGNG+8RU8PrsSgj2if8YdohgVLIiEMiw4fj0Hqx/NXzOFxAc4Da8\nsgwxBKchUFtpCYroRMQyvBCMSlYkIhle9mRtfc5rpr+daKW0FRDsHQimTlaGITgAORmG4BBk\nZBiCg5CPYQgOAgRTJxvDEByIXAxDcCgyMQzBoYBg6uRhGILDkYVhCA5IDoYhOCQZGIbgkEAw\nddIbhuCwJDcMwYFJbRiCQ5PYMASHBoKpk9YwBIcnqWEIjkBKwxAcAQimTkLDDk5Mm1PqxyB4\nnXSG7Z2YNuXQjzkl9utIZhiCI7Gz63Awzgnm+lcI3mJn4/BgeBP86zandOP5TGT4vGBUsmyA\nYOJkLRibU3rgqvdg3S8Eb3HNWvTELwRnyIGeLG1zyttt0pUFwfmBvmjiQDBxIJg4EEwcCCYO\nBBMHgokDwcTxK9gD//lIJP/MfF72bSfRcrLkP2TmFQimmdkABNPMbACCaWY2kJ1g4BcIJg4E\nEweCiQPBxMlM8PRhtgj5xcsp8r+sJy/By4magfOLllfsf9nArxZ8i5cXBI/EvOhxBacAgmPl\nhXuwImpQxcssam4av1dw3Lsi7sE9EWPqFrHpAsGKyNcAERyZ2JcAguMStdTsMoyXE2rRIAQQ\nTBwIJg4EEweCiQPBxIFg4kAwcSCYOBBMHAgmDgQTB4KJQ1zwK/vH+T92T30e6SAu+Ie9cH4X\nln8rxAXzd/b5l72lPouEUBecbKA9F8gL/svY39TnkBIIJg55wbeXFxTRhGkrWZ/sPfVZJIS4\n4K6Z9MJ+Up9HOogLVh0dr6nPIx3EBQMIJg4EEweCiQPBxIFg4kAwcSCYOBBMnP8Bo7Jvtz0j\n8GEAAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bias  %>% \n",
    "    ungroup()  %>% \n",
    "    mutate(nObs = as.factor(nObs))  %>% \n",
    "    ggplot(aes(x, subjError, color=nObs)) +\n",
    "    geom_point() +\n",
    "    geom_line() +\n",
    "    theme_classic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Increasing the number of observations definitely helps, thought it never really eliminates the problem. Increasing the number of observations reduced the noise in the obesrved responses, which reduces the bias problem discussed above. In the absence of more obesrvations, however, I'll have to content with fitting the overall data instead. I've thought about using a Bayesian hierarchical estimation, which can help with cases like this by \"shrinking\" - it imposes a distribution over subject parameters that would minimize the overfitting. Unfortunately, this is not feasible for the model I'm working in, which takes a long time to estimate even as is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortunately, after fitting my model to the averaged data, the problem disappears, as in the simulated example above. I guess that is that. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image.png](https://i.imgur.com/3EB4B1a.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}