{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Chapter 15: One Factor Models\n",
"\n",
"Load the packages"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import statsmodels.api as sm\n",
"import statsmodels.formula.api as smf\n",
"import seaborn as sns\n",
"from scipy import stats"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Coagulation example"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" coag \n",
" diet \n",
" \n",
" \n",
" \n",
" \n",
" 1 \n",
" 62 \n",
" A \n",
" \n",
" \n",
" 2 \n",
" 60 \n",
" A \n",
" \n",
" \n",
" 3 \n",
" 63 \n",
" A \n",
" \n",
" \n",
" 4 \n",
" 59 \n",
" A \n",
" \n",
" \n",
" 5 \n",
" 63 \n",
" B \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" coag diet\n",
"1 62 A\n",
"2 60 A\n",
"3 63 A\n",
"4 59 A\n",
"5 63 B"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coagulation = pd.read_csv(\"data/coagulation.csv\", index_col=0)\n",
"coagulation.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`pandas` version of boxplot"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEVCAYAAAAIK+VbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFlZJREFUeJzt3X+w5XV93/HnKyxEfmMELqQaNkwUdtgIZS+pCJi9rto0qVUzNM0mmSLudKfpZEfT6QTSbatOsx20ts1W27GOq8EMLlqUaKGuELhXSxtIdsk6Ra7UgYAYRNAqehFdd+fdP+7Zcl1hz9nd77lnz/k8HzN37jnfH5/7vp89+zrf+znf7+ebqkKS1JafGHUBkqTlZ/hLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4qxlJXpLkk0meTPLNJO9L8hNJ/kWSR5I8keQjSU5dss9/TfJ4kqeSfD7JBUvWvSjJf0vynSR/keQPktw1mt9OOjSGv5qQ5BjgFuARYCXwN4AbgTf3vmaAc4GTgPct2fUzwEuBM4F7gRuWrPtPwNPAWcBVvS9pLMS5fdSCJJcCnwbOrqq9S5bfAXyiqv5z7/l5wH3A8Uu36607DfgWcBqwAHwfWF1VD/TW/wGwtqouX4ZfSToiHvmrFS8BHjkw0IGfZvGvgf0eAVYAU0mOSXJdkgeTfAd4uLfN6cAZve0eXbLv0sfSUc3wVyseBX4myYoDlj8GnLPk+c8Ae4GvA78BvAF4DXAqi8NFAAGe7G334iX7vqTzqqUhMfzVij8HvgZcl+TEJC9IchmwHfjdJD+b5CTg3wAf6/2FcDLwA+CbwAm9dQBU1T7gk8A7kpyQ5HzgHy7vryQdPsNfTeiF9euBnwO+AnwV+AfAh4A/Bj4P/BWL4/ibert9hMVhoL8G7gfuPqDZ32HxL4LHe21sZ/HNQjrq+YGv1JEk7wLOqirP+tFRzyN/6TAlOT/Jy7PoF4ANwM2jrksaxIEffkka3MksDvX8NPAE8O+AT420ImlADvtIUoMc9pGkBhn+ktSgZR3zP/3002vlypXL+SMP29NPP82JJ5446jImkn07HPbrcIxbv+7atesbVXVGv+2WNfxXrlzJzp07l/NHHra5uTnWrl076jImkn07HPbrcIxbvyZ5pP9WDvtIUpMMf0lqkOEvSQ0y/CWpQYa/JDXI8Ney2b59O6tXr2bdunWsXr2a7du3j7qkiWC/6nA4t4+Wxfbt29m8eTPbtm1j3759HHPMMWzYsAGA9evXj7i68WW/6nB55K9lsWXLFrZt28bMzAwrVqxgZmaGbdu2sWXLllGXNtbsVx0uw1/LYn5+nssv/9H7ml9++eXMz8+PqKLJYL/qcBn+WharVq3irrvu+pFld911F6tWrRpRRZPBftXhMvy1LDZv3syGDRuYnZ1l7969zM7OsmHDBjZv3jzq0saa/arD5Qe+Whb7P3zctGkT8/PzrFq1ii1btvih5BGyX3W4DH8tm/Xr17N+/fqxmyjraGe/6nA47CNJDTL8JalBhr8kNcjwl6QG9Q3/JOcl2b3k6ztJ3pbkp5LcnuTLve8vXI6CJUlHrm/4V9UDVXVRVV0ErAG+B9wMXAvcUVUvBe7oPZckjYFDHfZZBzxYVY8AbwCu7y2/Hnhjl4VJkobnUMP/14H988VOVdXXAHrfz+yyMEnS8KSqBtswOQ54DLigqr6e5NtVddqS9d+qqh8b90+yEdgIMDU1tebGG2/spvIhW1hY4KSTThp1GRPJvh0O+3U4xq1fZ2ZmdlXVdL/tDuUK378D3FtVX+89/3qSs6vqa0nOBp54rp2q6gPABwCmp6drXK5A9GrJ4bFvh8N+HY5J7ddDGfZZz7NDPgCfBq7qPb4K+FRXRUmShmug8E9yAvBa4JNLFl8HvDbJl3vrruu+PEnSMAw07FNV3wNedMCyb7J49o8kacx4ha8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1aNDbOJ6W5KYkX0oyn+TSJBcluTvJ7iQ7k/zCsIuVJHVjoNs4AluBHVV1ZZLjgBOAjwPvrKrPJPll4N3A2uGUKUnqUt/wT3IK8CrgzQBVtQfYk6SAU3qbnQo8NqQaJUkdG+TI/1zgSeDDSS4EdgFvBd4GfDbJe1gcPnrl0KqUJHUqVXXwDZJp4G7gsqq6J8lW4DssHu1/rqo+keTXgI1V9Zrn2H8jsBFgampqzY033tj17zAUCwsLnHTSSaMuYyLZt8/a9MimUZfwvN57zntHXcJRYdxerzMzM7uqarrvhlV10C/gLODhJc+vAG4FnuLZN48A3+nX1po1a2pczM7OjrqEiWXfPuuca27prK0u+7XLusbduL1egZ3VJ4urqv/ZPlX1OPBokvN6i9YB97M4xv+LvWWvBr48+HuTJGmUBj3bZxNwQ+9Mn4eAq4FPAVuTrAC+T29oR5J09Bso/KtqN3DgGNJdwJrOK5IkDZ1X+EpSgwx/SWqQ4S9JDTL8JalBhr8kNWjQUz0lDcnKa2/trrEd3bR16vHHdtKOjl6GvzRCD1/3K521tfLaWzttT5PNYR9JapBH/upMkk7bqz6TDrbiUPo17xpsO/tWHvmrM4NMJlVVnHPNLQNtp0WD9uvs7OzA20qGvyQ1yGEf9XXhO2/jqWd+2GmbXZ3hcurxx/KFt7+uk7aklhj+6uupZ37Y6Vkkc3NzrF27tpO2Oj1NUmqIwz6S1CDDX5IaZPhLUoMMf0lq0EDhn+S0JDcl+VKS+SSX9pZvSvJAki8mefdwS5UkdWXQI/+twI6qOh+4EJhPMgO8AXh5VV0AvGdINUo6iO3bt7N69WrWrVvH6tWr2b59+6hL0hjoe6pnklOAVwFvBqiqPcCeJL8NXFdVP+gtf2KIdUp6Dtu3b2fz5s1s27aNffv2ccwxx7BhwwYA1q9fP+LqdDQb5Mj/XOBJ4MNJ/jLJB5OcCLwMuCLJPUk+l+SSoVYq6cds2bKFbdu2MTMzw4oVK5iZmWHbtm1s2bJl1KXpKDfIRV4rgIuBTVV1T5KtwLW95S8EXgFcAnw8ybl1wMQhSTYCGwGmpqaYm5vrsPzhWVhYGJtal0OXfdF137b87zQ/P8++ffuYm5v7//26b98+5ufnm+6XLk1sFgwwAdRZwMNLnl8B3ArsANYuWf4gcMbB2lqzZk2Ni9nZ2VGXcNQ455pbOm2vy77turZxc8EFF9Sdd95ZVc/265133lkXXHDBCKuaLOOWBcDOGmByv77DPlX1OPBokvN6i9YB9wN/ArwaIMnLgOOAb3T4viSpj82bN7NhwwZmZ2fZu3cvs7OzbNiwgc2bN4+6NB3lBp3bZxNwQ5LjgIeAq4GngQ8luQ/YA1zVe9eRtEz2f6i7adMm5ufnWbVqFVu2bPHDXvU1UPhX1W5g+jlW/Va35Ug6VOvXr2f9+vWdTpinyecVvpLUIMNfkhrkfP7q6+RV1/Lz11/bbaPXd9PMyasAurvXgNQKw199fXf+Om/mIk0Yh30kqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMGmtI5yWnAB4HVQAFvqao/6637Z8C/Bc6oKm/gPqE6nzp5RzftnXr8sZ20I7Vm0Pn8twI7qurK3k3cTwBI8hLgtcBXhlSfjgJdzuUPi28kXbcp6dD0HfZJcgrwKmAbQFXtqapv91b/B+D3WPxrQJI0JgY58j8XeBL4cJILgV3AW4F1wF9X1ReSPO/OSTYCGwGmpqaYm5s70pqXxcLCwtjUOo7s2+75mj10MzMznbY3OzvbaXvDlKqDH7QnmQbuBi6rqnuSbAX2sPjXwOuq6qkkDwPT/cb8p6ena+fOnd1UPmRd3mpQP8phn+HwNTsc4/Z6TbKrqqb7bTfI2T5fBb5aVff0nt8EXAz8LPCFXvC/GLg3yVmHWa8kaRn1Df+qehx4NMl5vUXrgHur6syqWllVK1l8g7i4t60k6Sg36Nk+m4Abemf6PARcPbySJEnDNlD4V9Vu4HnHkHpH/5KkMeEVvpLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNWjQi7wmxsEmoTtc/eZHasWh9G3e1X8b+1WH68J33sZTz/yws/a6vJ/Fqccfyxfe/rrO2jtczYX/oIEybpM5HQ0G7VsnINOwPfXMDzv7/9v167XzGyMdJod9JKlBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0a6CKvJKcBHwRWAwW8BfhV4PXAHuBB4Oqq+vaQ6hyIV/VJ0mAGvcJ3K7Cjqq7s3cf3BOB24Peram+SdwG/D1wzpDoH4lV9kjSYvsM+SU4BXgVsA6iqPVX17aq6rar29ja7G3jx8MqUJHVpkCP/c4EngQ8nuRDYBby1qp5ess1bgI89185JNgIbAaamppibmzuigvvpqv2FhYXOax327z4uhtG3sl8PZBb0UVUH/QKmgb3A3+o93wr86yXrNwM3A+nX1po1a2qYzrnmls7amp2d7aytqm5rG3dd960W2a/PajkLgJ3VJ4uraqCzfb4KfLWq7uk9vwm4GCDJVcDfBX6z90MlSWOgb/hX1ePAo0nO6y1aB9yf5JdY/ID371XV94ZYoySpY4Oe7bMJuKF3ps9DwNXAXwA/Cdzeu4nH3VX1j4dSpSSpUwOFf1XtZnHsf6mf674cSdJy8ApfSWpQc7dxlDT5Tl51LT9//bXdNXh9d02dvApg9LeINfwlTZzvzl/n1f59OOwjSQ0y/CWpQYa/JDXI8JekBhn+ktSgiTrbx9O7JGkwExX+nt4lSYNx2EeSGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYNFP5JTktyU5IvJZlPcmmSn0pye5Iv976/cNjFSpK6MeiR/1ZgR1WdD1wIzAPXAndU1UuBO3rPJUljoG/4JzkFeBWwDaCq9lTVt4E38OwECNcDbxxWkZKkbg1y5H8u8CTw4SR/meSDSU4EpqrqawC972cOsU5JUocGmdtnBXAxsKmq7kmylUMY4kmyEdgIMDU1xdzc3OHUObCu2l9YWOi81mH/7uNiGH0r+/VAZkEfVXXQL+As4OElz68AbgUeAM7uLTsbeKBfW2vWrKlhOueaWzpra3Z2trO2qrqtbdx13bdaZL8+q+UsAHZWnyyuqv7DPlX1OPBokvN6i9YB9wOfBq7qLbsK+FRn70iSpKEadErnTcANSY4DHgKuZvHzgo8n2QB8Bfj7wynx0HQ6dfKO7to69fhjO2tLUn9mwcENFP5VtRuYfo5V67ot58h0NZc/LL5wumxP0vIxC/rzCl9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwa9wndiJBl823cNtt3idBqSND6aO/IfZMKjqmJ2dnbgbSVp3DQX/pIkw1+SmmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1aKCLvJI8DHwX2AfsrarpJBcB7wdeAOwF/klV/fmwCpUkdedQrvCdqapvLHn+buCdVfWZJL/ce762y+IkScNxJMM+BZzSe3wq8NiRlyNJWg6DHvkXcFuSAv5LVX0AeBvw2STvYfFN5JVDqlGS1LFBw/+yqnosyZnA7Um+BFwJ/G5VfSLJrwHbgNccuGOSjcBGgKmpKebm5rqpfMgWFhbGptZxY98Oh/06PJPYrznUicmSvANYAP4lcFpVVRanynyqqk452L7T09O1c+fOw611Wc3NzbF27dpRlzGR7NvhsF+HY+W1t/Lwdb8y6jIGlmRXVU33267vmH+SE5OcvP8x8DrgPhbH+H+xt9mrgS8ffrmSpOU0yLDPFHBzbx78FcBHq2pHkgVga5IVwPfpDe1Iko5+fcO/qh4CLnyO5XcBa4ZRlCRpuLzCV5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkho0UPgneTjJ/06yO8nOJcs3JXkgyReTvHt4ZUqSujTIPXz3m6mqb+x/kmQGeAPw8qr6QZIzO69OkjQURzLs89vAdVX1A4CqeqKbkiRJwzZo+BdwW5JdSTb2lr0MuCLJPUk+l+SS4ZQoSeraoMM+l1XVY72hnduTfKm37wuBVwCXAB9Pcm5V1dIde28WGwGmpqaYm5vrrPhhWlhYGJtax419Oxz266GbmZkZaLu8a7D2Zmdnj6Ca5ZUDsrr/Dsk7gAXgNSwO+8z1lj8IvKKqnny+faenp2vnzp3Pt/qoMjc3x9q1a0ddxkSyb4fDfh2OcevXJLuqarrfdn2HfZKcmOTk/Y+B1wH3AX8CvLq3/GXAccA3nq8dSdLRY5Bhnyng5iT7t/9oVe1IchzwoST3AXuAqw4c8pEkHZ36hn9VPQRc+BzL9wC/NYyiJEnD5RW+ktQgw1+SGmT4S1KDDH9JapDhL0kNOuSLvI7ohyVPAo8s2w88MqfjdQvDYt8Oh/06HOPWr+dU1Rn9NlrW8B8nSXYOcpWcDp19Oxz263BMar867CNJDTL8JalBhv/z+8CoC5hg9u1w2K/DMZH96pi/JDXII39JapDh/xySvClJJTl/1LVMiiT7kuxO8oUk9yZ55ahrmhRJzkpyY5IHk9yf5L/3plnXYVryev1i7zX7T5NMVF467PMcknwcOBu4o6reMeJyJkKShao6qff4bwP/vKp+ccRljb0szrX+v4Drq+r9vWUXASdX1f8YaXFj7IDX65nAR4H/WVVvH21l3Zmod7IuJDkJuAzYAPz6iMuZVKcA3xp1ERNiBvjh/uAHqKrdBn93quoJFm9F+zu9N9uJMOg9fFvyRmBHVf2fJP83ycVVde+oi5oAxyfZDbyAxb+qXj3ieibFamDXqIuYdFX1UG/Y50zg66Oupwse+f+49cCNvcc39p7ryD1TVRdV1fnALwEfmaSjKDVhol6vHvkvkeRFLB6Rrk5SwDFAJfk9b1HZnar6sySnA2cAT4y6njH3ReDKURcx6ZKcC+xjgl6vHvn/qCuBj1TVOVW1sqpeAvwVcPmI65oovbOojgG+OepaJsCdwE8m+Uf7FyS5JIkfpnckyRnA+4H3TdJBoEf+P2o9cN0Byz4B/AbgB2hHZv+YPyz++XxVVe0bZUGToKoqyZuAP0xyLfB94GHgbSMtbPztf70eC+wF/hj496MtqVue6ilJDXLYR5IaZPhLUoMMf0lqkOEvSQ0y/CWpQZ7qKfUkeQewwOLcQ5+vqj89yLZvBm6rqseWpzqpW4a/dICq+lcDbPZm4D7A8NdYcthHTUuyOckDSf4UOK+37I+SXNl7vCbJ55LsSvLZJGf31k0DN/TmfD9+hL+CdFgMfzUryRoWp+3+m8CvApccsP5Y4L3AlVW1BvgQsKWqbgJ2Ar/Zm6zumeWtXDpyDvuoZVcAN1fV9wCSfPqA9eexOGXy7b0JSI8BvrasFUpDYvirdQeb3yTAF6vq0uUqRlouDvuoZZ8H3pTk+CQnA68/YP0DwBlJLoXFYaAkF/TWfRc4eflKlbrlkb+aVVX3JvkYsBt4hANmbq2qPb0Pd/9jklNZ/P/yhyzOof9HwPuTPANc6ri/xo2zekpSgxz2kaQGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXo/wG8btw1K/8laQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"coagulation.boxplot('coag',by='diet')\n",
"plt.suptitle(\"\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`seaborn` version of boxplot"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFgNJREFUeJzt3X+QXWd93/H3R7ZsywZjsNbgYVFcV7aZwoCxlxRDMNjGLkoJv8alIWFmwZlokk6lQpsmTtNpyTR0XELbeJ0pROMfiI4JMQYDzVRGCp2QpgSHlRGOMXYkXAFby/ZKxomFZEu2vv3j3g220I8r7T17dve8XzOau/fcc+756szd/dznOec8T6oKSVJ3LWm7AElSuwwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjTmy7gEEsX768zjnnnLbLkKQFZfPmzTurauRo6y2IIDjnnHOYnJxsuwxJWlCSfG+Q9ewakqSOMwgkqeMMAknqOINAkjrOINC8tHPnTtasWcOuXbvaLmXB81jqaAwCzUvr16/nnnvuYf369W2XsuB5LHU0BoHmnZ07d7Jhwwaqig0bNvhNdhY8lhqEQaB5Z/369cxMoXrgwAG/yc6Cx1KDMAg072zatIn9+/cDsH//fjZu3NhyRQuXx1KDMAg071x55ZUsXboUgKVLl3LVVVe1XNHC5bHUIAwCzTvj4+MkAWDJkiWMj4+3XNHC5bHUIAwCzTvLly9n1apVJGHVqlWceeaZbZe0YHksNYgFMeicumd8fJzt27f7DXYIPJY6msxcUTCfjY2NlaOPStKxSbK5qsaOtp5dQ5LUcQaBJHWcQSBJHWcQSFLHNRYESS5IsuVZ//42yQeTvCjJpiRb+48vbKoGSdLRNRYEVfVAVV1YVRcCFwN7gDuAa4GvVNV5wFf6zyVJLZmrrqErgO9W1feAdwAzI1+tB945RzVIkg5hroLg54E/7P/84qraAdB/POtQGyRZnWQyyeT09PQclSlJ3dN4ECQ5CXg78Nlj2a6q1lXVWFWNjYyMNFOcJGlOWgSrgLur6pH+80eSnA3Qf3x0DmqQJB3GXATBe/lxtxDAl4CZQU/GgS/OQQ2SpMNoNAiSnApcCXz+WYuvA65MsrX/2nVN1iBJOrJGRx+tqj3AmQct20XvKiJJ0jzgncWS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxTU9VeUaS25Pcn+Q7SS5JcmGSryfZkmQyyU83WYMk6cganaoSuB64s6quTnIScCpwG/DbVbUhyc8CHwXe3HAdkqTDaCwIkpwOXAq8H6Cq9gH7khRwen+1FwAPNVWDJOnommwRnAtMA7ckeTWwGfgXwAeBLyf5GL2uqdc3WIPUqmuuuYYdO3Yc9/ZPPfUUBw4cGGJFx2/JkiWcfPLJs3qPs88+m5tvvnlIFWlYmgyCE4GLgDVVdVeS64Fr6bUCPlRVn0vyHuAm4C0Hb5xkNbAaYMWKFQ2WKTXn8ccfZ/ee3cf/m3YAqGFWdPyeOfAM+/ftP/43eLp3PDT/NBkEU8BUVd3Vf347vSD4GXotA4DPAjceauOqWgesAxgbG5snvwrSsRkdHWU60xx48/z4Vt+mJX+6hNGXjrZdhg6hsauGquph4AdJLugvugK4j945gTf1l10ObG2qBknS0TV91dAa4Nb+FUMPAh8Avghcn+RE4En63T+SpHY0GgRVtQUYO2jxnwMXN7lfSdLgvLNYkjrOIJCkjjMIJKnjDAJJ6rimrxqS9HjvGvrW7O4/Pq+9EgB4HHhpyzXokAwCqUErV65suwS2bu3dqnPeS89rt5CXzo/joZ9kEEgNWrt2bdsl/F0NExMTLVei+cog0NBNTEywbdu2Wb3H1NQU0Bui4XitXLlyXvwhno1hHMuZFsFsj8ViOJ46NINA89LevXvbLmHRWLZsWdslaJ5L1fwfz21sbKwmJyfbLkNzyO4MafaSbK6qg0d3+Am2CPQcw+iKGIZhdWfMlt0h6gKDQM+xbds2/vreu1nxvGdareOk/b3LLZ/c/o3Wavj+7hNa27c0lwwC/YQVz3uGfzu2++grLnK/M9n2hffS3PDOYknqOINAkjrOIJCkjjMIJKnjGg2CJGckuT3J/Um+k+SS/vI1SR5I8u0kH22yBqnrdu7cyZo1a9i1a1fbpWiearpFcD1wZ1W9HHg18J0klwHvAF5VVa8APtZwDVKnrV+/nnvuuYf169e3XYrmqcaCIMnpwKXATQBVta+qHgd+Fbiuqp7qL3+0qRqkrtu5cycbNmygqtiwYYOtAh1Sky2Cc4Fp4JYk30xyY5LTgPOBNya5K8lXk7y2wRqkTlu/fj0zw8gcOHDAVoEOqckgOBG4CPh4Vb0G+BFwbX/5C4HXAf8auC1JDt44yeokk0kmp6enGyxTWrw2bdrE/v37Adi/fz8bN25suSLNR00GwRQwVVV39Z/fTi8YpoDPV89fAgeA5QdvXFXrqmqsqsZGRkYaLFNavK688kqWLl0KwNKlS7nqqqtarkjzUWNBUFUPAz9IckF/0RXAfcAXgMsBkpwPnATsbKoOqcvGx8eZaXAvWbKE8fHxlivSfNT0VUNrgFuT3ANcCPxH4Gbg3CT3Ap8BxmshjIUtLUDLly9n1apVJGHVqlWceeaZbZekeajRQeeqagtwqLGw39fkfiX92Pj4ONu3b7c1oMNy9FFpkVu+fDk33HBD22VoHjMI9BxTU1P86IkTHIIZ+N4TJ3Baf+5kaTFzrCFJ6jhbBHqO0dFRnnx6hxPT0JuY5pTR0bbLkBpni0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs7RR/UTvr+7/fkIHtnT+47y4lMPtFbD93efwPmt7V2aOwMFQZIXHWLxE1W1/yjbnQHcCLwSKOCaqvqL/mu/BvwuMFJVTl4/T6xcubLtEgDYt3UrAKecc15rNZzP/DkeUpMGbRHcDbwM+CEQ4AxgR5JHgV+uqs2H2e564M6qujrJScCpAEleBlwJfH82xWv41q5d23YJwI/rmJiYaLkSafEbNAjuBO6oqi8DJLkKeCtwG/DfgH948AZJTgcuBd4PUFX7gH39l/8r8OvAF2dRuyQdk4mJCbZt2zar95jqT186OstJi1auXDlvvngNerJ4bCYEAKpqI3BpVX0dOPkw25wLTAO3JPlmkhuTnJbk7cD/q6pvHWmHSVYnmUwyOT09PWCZktSsvXv3snfv3rbLGKpBWwSPJfkN4DP95/8U+GGSE4DDnc07EbgIWFNVdyW5HvgwvVbCVUfbYVWtA9YBjI2N1YB1StJhDeMb+GLsthy0RfALwCjwBXrdOSv6y04A3nOYbaaAqaq6q//8dnrB8PeAbyXZ3n/Pu5O85LiqlyTN2kAtgv5VPWsO8/IhO9yq6uEkP0hyQVU9AFwB3F1VV8ys0w+DMa8akqT2DHr56Ai9k7uvAE6ZWV5Vlx9l0zXArf0rhh4EPnCcdUqSGjLoOYJbgT8C3gb8CjBO70TwEVXVFmDsCK+fM+D+JUkNGfQcwZlVdROwv6q+WlXXAK9rsC5J0hwZtEUwcwfxjiT/GHiI3oleSdICN2gQ/E6SFwD/CrgBOB34UGNVSZLmzKBXDf1x/8e/AS5rrhxJ0lwb6BxBktEkdySZTvJIks8lsWtIkhaBQU8W3wJ8CTgbeCnwP/rLJEkL3KBBMFJVt1TV0/1/nwRGGqxLkjRHBj1ZvDPJ+4A/7D9/L7CrmZLmniMSDtcwjufW/nwEszkWi+FY6seG8bkahmF8NodhmJ/vQYPgGuD36Q0fXcDX8C7h51hsoxG2bdmyZW2XoHlm27Zt3L9lC20PTDbTjfL4li2t1fDwkN9v0CD4D8B4Vf0Q/m7Gso/RC4gFzxEJh6vtb0pavF4C/BJpu4zW3cRwB2Qe9BzBq2ZCAKCqHgNeM9RKJEmtGDQIliR54cyTfovAie8laREY9I/5fwa+luR2eucI3gN8pLGqJElzZtA7iz+VZBK4nN7k9e+uqvsarUySNCcG7t7p/+H3j78kLTKDniOQJC1SBoEkdVyjV/4kOQO4EXglvZPM1wDvBn4O2Ad8F/hAVT1+vPvwbsPn8m5aSceq6UtArwfurKqr+/MWnwpsAn6zqp5O8p+A3wR+43h3sG3bNr75V/dx4NQXDafi45R9vRs8Nn932Pf8DW7Jnsda27ekhauxIEhyOnAp8H6AqtpHrxWw8VmrfR24erb7OnDqi3jyH7xttm+z4J1y3x8ffSVJOkiT5wjOpTfB/S1JvpnkxiSnHbTONcCGBmuQJB1Fk0FwInAR8PGqeg3wI+DamReT/BbwNHDroTZOsjrJZJLJ6enpBsuUpG5rMgimgKmquqv//HZ6wUCSceBtwC9W1SFHT6qqdVU1VlVjIyNOfSBJTWksCKrqYeAHSS7oL7oCuC/JW+mdHH57Ve1pav+SpME0fdXQGuDW/hVDD9Kbw+AbwMnApiQAX6+qX2m4DknSYTQaBFW1BRg7aPHKJvcpSTo2DiUtaUGYmpriCYY/KctCtAPY3Z8edxgcYkKSOs4WgaQFYXR0lMd37nSqSnqtojNGR4f2frYIJKnjDAJJ6jiDQJI6ziCQpI5b8CeLp6amWLLnbxx5E1iyZxdTU0+3XYakBcYWgSR13IJvEYyOjvLIUyc6HwG9+QhGR1/SdhmSFhhbBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxzUaBEnOSHJ7kvuTfCfJJUlelGRTkq39xxc2WYMk6ciabhFcD9xZVS8HXg18B7gW+EpVnQd8pf9cktSSxoIgyenApcBNAFW1r6oeB94BrO+vth54Z1M1SJKOrskWwbnANHBLkm8muTHJacCLq2oHQP/xrENtnGR1kskkk9PT0w2WKUnd1mQQnAhcBHy8ql4D/Ihj6AaqqnVVNVZVYyMjI03VKEmd1+Too1PAVFXd1X9+O70geCTJ2VW1I8nZwKOz3dGSPY+1Ph9BnvxbAOqU01urYcmexwBHH9Xi9TC9idvbtKv/eGaLNTwMnDHE92ssCKrq4SQ/SHJBVT0AXAHc1/83DlzXf/zibPazcuXKWdc6DFu3PgHAeX+/zT/EL5k3x0Matvny2Z7euhWAM847r7UazmC4xyNVzaVrkguBG4GTgAeBD9DrjroNWAF8H/gnVfXYkd5nbGysJicnG6tzGNauXQvAxMREy5VIatJC+l1Psrmqxo62XqMT01TVFuBQRVzR5H4lSYPzzmJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOq7Ry0cXiomJCbZt2zar99jav8lk5hrj47Vy5cpZv4ckHQuDYEiWLVvWdgmSdFwMAmb/LV6SFjLPEUhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHdfoDWVJtgNPAM8AT1fVWH/6yk8ApwBPA/+sqv6yyTokSYc3F3cWX1ZVO5/1/KPAb1fVhiQ/23/+5jmoQ5J0CG10DRVwev/nFwAPtVCDJKmv6RZBARuTFPAHVbUO+CDw5SQfoxdEr2+4BknSETQdBG+oqoeSnAVsSnI/cDXwoar6XJL3ADcBbzl4wySrgdUAK1asaLhMSequRruGquqh/uOjwB3ATwPjwOf7q3y2v+xQ266rqrGqGhsZGWmyTEnqtMaCIMlpSZ4/8zNwFXAvvXMCb+qvdjmwtakaJElH12TX0IuBO5LM7OfTVXVnkt3A9UlOBJ6k3/0jSWpHY0FQVQ8Crz7E8j8HLm5qv5KkY+OdxZLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HGNBkGS7Un+KsmWJJPPWr4myQNJvp3ko03WIEk6sibnLJ5xWVXtnHmS5DLgHcCrquqpJGfNQQ2SpMNoo2voV4HrquopgKp6tIUaJEl9TbcICtiYpIA/qKp1wPnAG5N8BHgS+LWq+kbDdUgSExMTbNu2bVbvsXXrVgDWrl07q/dZuXLlrN9jWJoOgjdU1UP97p9NSe7v7/OFwOuA1wK3JTm3qurZGyZZDawGWLFiRcNlStJgli1b1nYJQ5eD/v42t6Pkw8Bu4C30uob+tL/8u8Drqmr6cNuOjY3V5OTk4V6WJB1Cks1VNXa09Ro7R5DktCTPn/kZuAq4F/gCcHl/+fnAScDOw72PJKlZTXYNvRi4I8nMfj5dVXcmOQm4Ocm9wD5g/OBuIUnS3GksCKrqQeDVh1i+D3hfU/uVJB0b7yyWpI4zCCSp4wwCSeo4g0CSOs4gkKSOm7MbymYjyTTwvbbrGMByvCdimDyew+OxHK6Fcjx/qqpGjrbSggiChSLJ5CB38WkwHs/h8VgO12I7nnYNSVLHGQSS1HEGwXCta7uARcbjOTwey+FaVMfTcwSS1HG2CCSp4wyCIUjyriSV5OVt17LQJXkmyZYk30pyd5LXt13TQpbkJUk+k+S7Se5L8j/7w7/rGD3rs/nt/ufzXyZZFH9D7RoagiS3AWcDX6mqD7dczoKWZHdVPa//8z8C/k1Vvanlshak9MaA/xqwvqo+0V92IfD8qvrfrRa3AB302TwL+DTwf6rq37db2ewtijRrU5LnAW8Afgn4+ZbLWWxOB37YdhEL2GXA/pkQAKiqLYbA7FXVo/Sm0v3n/cBd0Jqes7gL3gncWVV/neSxJBdV1d1tF7WALUuyBTiFXivr8pbrWcheCWxuu4jFqqoe7HcNnQU80nY9s2GLYPbeC3ym//Nn+s91/PZW1YVV9XLgrcCnFsM3Li1ai+KzaYtgFpKcSe8b6yuTFHACUEl+3ek3Z6+q/iLJcmAEeLTtehagbwNXt13EYpXkXOAZFsFn0xbB7FwNfKqqfqqqzqmqlwH/F/iZlutaFPpXYZ0A7Gq7lgXqfwEnJ/nlmQVJXpvEk++zlGQE+ATw+4vhS58tgtl5L3DdQcs+B/wC4Am54zNzjgB6ze7xqnqmzYIWqqqqJO8Cfi/JtcCTwHbgg60WtnDNfDaXAk8D/x34L+2WNBxePipJHWfXkCR1nEEgSR1nEEhSxxkEktRxBoEkdZyXj0oDSvJhYDe9MZD+rKr+5Ajrvh/YWFUPzU110vEzCKRjVFX/boDV3g/cCxgEmvfsGpKOIMlvJXkgyZ8AF/SXfTLJ1f2fL07y1SSbk3w5ydn918aAW/vj1y9r8b8gHZVBIB1GkovpDS3+GuDdwGsPen0pcANwdVVdDNwMfKSqbgcmgV/sD6C3d24rl46NXUPS4b0RuKOq9gAk+dJBr19Ab6jnTf0BUk8AdsxphdIQGATSkR1pDJYA366qS+aqGKkJdg1Jh/dnwLuSLEvyfODnDnr9AWAkySXQ6ypK8or+a08Az5+7UqXjZ4tAOoyqujvJHwFbgO9x0IiyVbWvf2J4IskL6P0+/R69eQA+CXwiyV7gEs8TaD5z9FFJ6ji7hiSp4wwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjvv/RrBdIwvNslIAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sns.boxplot(x=\"diet\", y=\"coag\", data=coagulation)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Maybe seaborn plot is preferable (although without the colors). Pandas version is annoying over the title and failure to label the y-axis.\n",
"\n",
"`swarmplot` is good when there are replicated points:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAG2hJREFUeJzt3XucVXW9//HXe27AjNwGB0SFAE09at4YNCxRsdIMNTtkp7LQSk43fyc76bF6nA6d6uij0820o/HT0Hp4CVFP0AVFy7RfhQ0XTZSLggoCMqAg98vM5/fH3qAMMGycvWbtPev9fDx47P1de6/Zb+cxznvWd+313YoIzMwsuyrSDmBmZulyEZiZZZyLwMws41wEZmYZ5yIwM8s4F4GZWca5CMzMMs5FYGaWcS4CM7OMq0o7QCEOPvjgGDJkSNoxzMzKyqxZs1ZHRMP+nlcWRTBkyBCamprSjmFmVlYkvVjI8zw1ZGaWcS4CM7OMcxGYmWWci8DMLONcBGZmGeciMDPLOBeBmVnGuQis9Gx6FZ7/Q+7WOmzJuiX8beXf2NG6I+0oVqLK4oIyy5D5v4Upl8OOLVDVA8b+DI45P+1UZeu6mddx1/y7ABjUcxCTzp3EgLoBKaeyUuMjAistD341VwIAOzbDg19LN08ZW7x28a4SAFi6fim3z7s9vUBWslwEVlrWv9JmvDKdHF3Aqs2r9tjWvLk5hSRW6lwEVlpOuKT9sRVseP/hHFp36G7bxgwbk1IaK2U+R2Cl5fzvQf1QeGkmDD4N3vmFtBOVrerKaiadN4lJT09izZY1jBk2hrMGnZV2LCtBioi0M+xXY2NjePVRM7MDI2lWRDTu73meGjIzyzgXgZlZxrkIzMwyzkVgZpZxLgIzs4xLrAgkHS1p7pv+vS7pS5LqJc2QtCh/2zepDGZmtn+JFUFELIiIkyLiJGA4sAl4ALgWeCQi3g48kh+bmVlKOmtq6Bzg+Yh4EbgIuCO//Q7gg52UwczM9qKziuCfgLvz9wdExAqA/G3/TspgZmZ7kXgRSKoBLgTuPcD9xktqktTU3OyFsszMktIZRwTvB2ZHxM5lJV+RNBAgf7vnEolAREyMiMaIaGxoaOiEmGZm2dQZRfBR3pgWApgKjMvfHwf8qhMymJnZPiRaBJJqgfcC979p8/XAeyUtyj92fZIZzMysfYkuQx0Rm4B+bbatIfcuIjMzKwG+stjMLONcBGZmGeciMDPLOBeBmVnGuQjMzDLORWBmlnEuAjOzjHMRmJllnIvAzCzjXARmZhnnIjAzyzgXgZlZxrkIzMwyzkVgZpZxLgIzs4xzEZiZZZyLwMws45L+qMo+kqZImi/pWUkjJZ0k6a+S5kpqknRqkhnMzKx9iX5UJXADMD0ixkqqAWqBycA3I+J3ks4HvguclXAOsy5v5caVTHt+GtUV1Vx45IXUd69n0/ZNTHt+Gqu3rOa8IedxRJ8jiAgeeekR5q2Zx4gBIzj9sNPTjm4pS6wIJPUCRgGXAUTENmCbpAB65Z/WG1ieVAazrFi5cSVjp41l3dZ1ANw1/y6mXDCFzz38OZ5a/RQAt/39Nm4/73YefvFhJs2bBMCtf7+Va0ZcwyeO/URq2S19SU4NDQOagUmS5ki6VVId8CXgvyUtBb4HfDXBDGaZMPX5qbtKAGDFxhVMmjdpVwkAbG/dzt3z7+bu+Xfvtu8vnvlFp+W00pRkEVQBpwA3R8TJwEbgWuBzwFURMQi4CrhtbztLGp8/h9DU3NycYEyz8ldVsefBfU1FzR7bqiuqqayo3O++li1JFsEyYFlEzMyPp5ArhnHA/flt9wJ7PVkcERMjojEiGhsaGhKMaVb+LjriIgbUDtg1PqL3EVx2/GWcfugb8/911XVceuylfPr4T+/aJsQV77iiU7Na6UnsT4GIWClpqaSjI2IBcA7wDLkpozOBR4HRwKKkMphlRb8e/bj/ovuZ8cIMaiprOGfwOfSo6sFPzvkJjy59lObNzZwz+Bz61/bnqL5HMXzAcOatmUfjgEb+od8/pB3fUqaISO6LSycBtwI1wGLgcuA4cu8mqgK2AJ+PiFntfZ3GxsZoampKLKeZWVckaVZENO7veYlODkbEXKBtiD8Bw5N8XTMzK5yvLDYzyzgXgZlZxrkIzMwyzkVgZpZxLgIzs4zzJYVWWlpbYe6dsHQmDDoNTvo4VPjvlbfqtS2vcc/8e1i9eTUfGPYBThlwStqRrAS5CKy0zPh3+MtNuftzfgHN8+Hc76SbqUy1tLZw+fTLeX7d8wBMWTSFie+dyGkDT0s5mZUa/6llpWXW7e2PrWCzV83eVQIArdHKfYvuSzGRlSoXgZWW7r3bH1vBetX0KmibmYvASss53wDlfyxVkRvbW3J0/dG8f8j7d43ru9fzyWM/mWIiK1WJrjVULF5rKGNeXQLLmuDwRqgfmnaasjfrlVms3ryadx/2buqq69KOY52oJNYaMntL6oe6AIpo+AAv7WXt89SQmVnGuQjMzDLORWBmlnEuAjOzjHMRmJllXKJFIKmPpCmS5kt6VtLI/PYrJS2QNE/Sd5PMYGZm7Uv67aM3ANMjYqykGqBW0tnARcAJEbFVUv+EM1g52bEV/vITWPoEDDoVRn4RqmrSTlW2lm9YzqSnJ7FmyxrGDBvD6MGj045kJSixIpDUCxgFXAYQEduAbZI+B1wfEVvz21cllcHK0G+vhtl35O4v/B289gJc+ONUI5Wr7S3buXz65SzfuByAGS/O4MbRN3LWoLPSDWYlJ8mpoWFAMzBJ0hxJt0qqA44CzpA0U9IfJY1IMIOVm6cmtz+2gs1aNWtXCez068W/TimNlbIki6AKOAW4OSJOBjYC1+a39wXeCVwNTJaktjtLGi+pSVJTc3NzgjGtpPQc0GZ8SDo5uoD+PfacdW3o0ZBCEit1SRbBMmBZRMzMj6eQK4ZlwP2R8wTQChzcdueImBgRjRHR2NDgH97MOPc6qOqRu1/VA879r3TzlLFhfYbxsWM+tms8qOcgLjvusvQCWclKdNE5SY8Dn4mIBZImAHXA88ChEfENSUcBjwCDo50gXnQuYza9CiuehIEnQm192mnK3pJ1S1i9eTUn9z+ZqgovL5YlpbLo3JXAnfl3DC0GLic3RfQzSU8D24Bx7ZWAZVBtPRxxdtopuoyhvYcytLcX8bN9S7QIImIusLc2ujTJ1zUzs8L5ymIzs4xzEZiZZZyLwMws41wEZmYZ5yIwM8s4v6nYSsvW9fCH62DpTBh0Gpz9VejWM+1UZl2ai8BKy7R/gafvy91/uQk2vAJjb0s3k1kX56khKy3PTmsznppODrMMcRFYaenb5grY+mHp5DDLEBeBlZYPfB9q++Xu1/aD87+Xbh6zDPA5AistQ8+Aq56BNc9BvyOhunvaicy6PBeBlZ7q7nDI8WmnMMsMTw2ZmWWci8DMLONcBGZmGeciMDPLOBeBmVnGFVQEkur38q+6gP36SJoiab6kZyWNfNNjX5EUkvb44HozM+s8hR4RzAaagYXAovz9JZJmSxrezn43ANMj4hjgROBZAEmDgPcCL73V4NZFbVwN910BPzohd7txddqJzHZ59a67WHzBBSz5yEfY8NhjaccpmkKvI5gOPBARDwJIeh9wHjAZ+B/gtLY7SOoFjAIuA4iIbeQ+rB7gh8A1wK86kN26ol99ARZOz91f+yJsfR0+9st0M5kB63//B175z2/tGi/9whc5cvrvqD7ssBRTFUehRwSNO0sAICIeAkZFxF+BbvvYZxi5I4dJkuZIulVSnaQLgZcj4skOJbeu6bmH24wfSSeHWRsb//Sn3Tds387GmU+kE6bICi2CVyX9m6S35f9dA7wmqRJo3cc+VcApwM0RcTKwEZgAfB34xv5eUNJ4SU2SmpqbmwuMaWWv/7G7jwccu/fnmXWybkcdVdC2clRoEXwMOBz4X3LTOYPz2yqBS/axzzJgWUTMzI+nkCuGocCTkl7If83Zkg5pu3NETIyIxohobGhoKDCmlb0Lb3xjBdK+Q+GCH6ebxyyvz4cuptcFF0BFBerenYOv/CI9jj8u7VhFoYhI7otLjwOfiYgFkiYAdRFx9Zsef4HctFO7ZwQbGxujqakpsZxWYiJgwyo4qD9Iaacx203LunWoupqK2tq0o+yXpFkR0bi/5xV0slhSA7mTu8cBu5aDjIjR+9n1SuBOSTXAYuDyQl7PMk6CngPSTmG2V5W9e6cdoegKfdfQncAvgTHAZ4Fx5E4Etysi5gL7bKOIGFLg65uZWUIKPUfQLyJuA7ZHxB8j4lPAOxPMZWZmnaTQI4Lt+dsVkj4ALCd3otfMzMpcoUXwbUm9gX8FbgR6AVcllsrMzDpNQUUQEb/O310HnJ1cHDMz62yFLjp3uKQHJDVLekXSfZI8NWRm1gUUerJ4EjAVGAgcBkzLb8u8HS2tTJg6jxO/+RDv+cEf+cP8VWlHKm9rX4KffxC+c2judu3StBNZBrWsX8/LX/5XFpwynCWXfITNT88DoPnGm1h4+rt47j3vZd3UqSmnLJ6CLiiTNDciTtrftqSU8gVltz6+mG//5tld4+7VFfzl2nPoW1eTYqoydvsYeOHxN8ZDR8G4aenlsUxaMWECa+95Y7HD6kMPpeHLV7H8K1e/8aSKCob95td0Gzo0hYSFKfSCskKPCFZLulRSZf7fpcCajkXsGp5Y8upu4y3bW3ly2dqU0nQBL/65/bFZJ9jU5g/P7cuX77nsdGsrm2fP6cRUySm0CD5Fbk2hlcAKYCy+ShiAkwf33W1cU1nBcYd2vSsPO83hI3YfH7bfP2bMiq7HCSfuNq7q35/a0/ZYbZ8eJ7yjsyIlqtAi+BYwLiIaIqI/uWKYkFiqMvKpdw9h7PDDqa4Uh/Tqzg8+ciINPfe1Mrft10U/eeOX/2GNubFZJ+t/9VeoO3MUSNQMG8ZhP/wBfS6+mPpx41D37lT27csh//ENur397WlHLYpCzxHMyS8l3e62pJTyOYKdWluDigovkFY0rS1QUZl2Csu4aG1FFRX73Vaqin2OoELSrjkQSfUUfjFaJrgEiswlYCVgb7/wy6UEDkShv8y/D/xZ0hQgyJ0v+E5iqczMrNMUemXxzyU1AaMBAR+KiGcSTWZmZp2i4Omd/C9+//I3M+tiut5kl5mZHRAXgZlZxiVaBJL6SJoiab6kZyWNlPTf+fFT+YXs+iSZobNs2d5Ca2tyn/9sZqUhtm0jduxIO0ZRJX1EcAMwPSKOAU4EngVmAMdHxAnAQuCrCWdI1JbtLVx59xyO+48HOfW/HuaBOcvSjmRmCYiWFlZ885ssGN7IwpGn8+odd6QdqWgSKwJJvYBRwG0AEbEtItZGxEMRsbNO/0qZf9LZbX9awrQnl9PSGqzesI2r732KVa9vSTuWmRXZuv/9FWvvvofYvp3W9et55brr2bJgQdqxiiLJI4Jh5D7gfpKkOZJulVTX5jmfAn63t50ljZfUJKmpubk5wZgd81SbBeZ2tAbPrHg9pTRmlpQt857ec1t+eepyl2QRVAGnADfnl6LYCFy780FJXwd2AHfubeeImBgRjRHR2NDQkGDMjhk5rN9u4x7VlZw8qO8+nm1m5ar21DaLzlVWUjuiayyKmGQRLAOWRcTM/HgKuWJA0jhgDPDxKGSxoxL2iZFD+Oczh9HQsxvHHdqLiZ8cTu/a6rRjmVmR9TrvXBq+/GWqBg6k5ogjOOz736Nm8OC0YxVFQYvOveUvLj0OfCYiFkiaANQBjwA/AM6MiILmfMph0Tkzs1JT6KJzSS8cdyVwp6QaYDG5zzD4G9ANmCEJ4K8R8dmEc5iZ2T4kWgQRMRdo20ZHJvmaZmZ2YHxlsZlZxrkIzMwyzkVgZpZxLgIzs4xzEZiZHYAdq1fT8nrXWj3ARWBmVoDWbdtY9n/+hUVnjGLRu97Nqh/+KO1IReMiMDMrwNopU1j/0EMQQWzfzpqf/pTNTz2VdqyicBGYmRVg23PP77Ft6162lSMXgZlZAQ4668zdxqqpoe70kSmlKS4XgZlZAQ4aNYqB3/4W3Y89ltrGRgbdcjPVhxySdqyiSHqtITOzLqPP2LH0GTs27RhF5yMCM7OMcxGYmWWci8DMLONcBGZmGeciMDPLOBeBmVnGJVoEkvpImiJpvqRnJY2UVC9phqRF+du+SWboDBHBvOXrWLluS9pRzCxhWxYuZNuyl9OOUVRJX0dwAzA9IsbmP7e4Fvga8EhEXC/pWuBa4N8SzpGY1Ru2cumtM5m/cj0Vgs+fdSRfOffotGOZWZG1btzI0n/+LJuamgDo8+GxDPzWt1JOVRyJHRFI6gWMAm4DiIhtEbEWuAi4I/+0O4APJpWhM0x8bDHzV64HoDXgpj88x5LVG1NOZWbF9trke3eVAMDae6ew6W9/SzFR8SQ5NTQMaAYmSZoj6VZJdcCAiFgBkL/tv7edJY2X1CSpqbm5OcGYHbPstU0FbTOz8rb95T2ng7btZVs5SrIIqoBTgJsj4mRgI7lpoIJExMSIaIyIxoaGhqQydtj57xi427h/z26MGFKfUhozS0qvc98H0q5xRV0dB51xRoqJiifJcwTLgGURMTM/nkKuCF6RNDAiVkgaCKxKMEPixpxwKNt2tHL/7Jdp6NmNL44+ku7VlWnHMrMiqx0xgsNvupHX7vklFbW19LviCqr69Us7VlEoIpL74tLjwGciYoGkCUBd/qE1bzpZXB8R17T3dRobG6PpTXNzZma2f5JmRUTj/p6X9LuGrgTuzL9jaDFwObnpqMmSPg28BHw44QxmZtaORIsgIuYCe2ujc5J8XTMzK5yvLDYzyzgXgZlZxrkIzMwyzkVgZpZxLoIi2LqjhccWNvPsitfTjmJmCYqWFjb+dSabn3oq7ShF5Q+v76CX127mklv+wstrNwPw0VMHc92H3pFyKjMrtpa1a3nxE59g66LnADho9GgO/8lN6E1XG5crHxF00P99bPGuEgC4+4mXWPjK+hQTmVkSXpt8764SANjw+9+z8c9/TjFR8bgIOmj1hq17blu/5zYzK28ta1bvue3VV1NIUnwugg76x1MO3208uL6WEUO96JxZV9NrzAVQ9cZsemXfvhx05pkpJioenyPooLOP6c/tl4/ggTkv03BQNz5zxjCqK92vZl1Nj3ccz9vuuJ21k++loq6W+k9+kspevdKOVRSJLjpXLF50zszswBW66Jz/dDUzyzgXgZlZxrkIzMwyzkVgZpZxLgIzs4xL9O2jkl4A1gMtwI6IaJR0EnAL0B3YAXw+Ip5IMoeZme1bZ1xHcHZEvPmSvO8C34yI30k6Pz8+qxNymJnZXqQxNRTAzqswegPLU8hgZmZ5SR8RBPCQpAB+GhETgS8BD0r6HrkiOj3hDGZm1o6ki+BdEbFcUn9ghqT5wFjgqoi4T9IlwG3Ae9ruKGk8MB5g8ODBCcc0M8uuTltiQtIEYAPw70CfiAjlFvJeFxHtLtjhJSbMzA5c6ktMSKqT1HPnfeB9wNPkzgnsXLJvNLAoqQxmZrZ/SU4NDQAeyH96TxVwV0RMl7QBuEFSFbCF/PSPmZmlI7EiiIjFwIl72f4nYHhSr2tmZgfGVxabmWWci8DMLONcBGZmGeciMDPLOBeBmVnGuQjMzDLORWBmlnEuAjOzjHMRmJllnIvAzCzjXARmZhnnIjAzyzgXgZlZxrkIzMwyzkVgZpZxLgIzs4xzEZiZZVyiRSDpBUl/lzRXUtObtl8paYGkeZK+m2QGMzNrX5KfWbzT2RGxeudA0tnARcAJEbFVUv9OyGBm1mFblyxh3dSpVPSopc/Yf6Sqvj7tSEXRGUXQ1ueA6yNiK0BErEohg5nZAdm6eDFLxn6Y2LQJgLWTJzNs2lQqevRIOVnHJX2OIICHJM2SND6/7SjgDEkzJf1R0oiEM5iZddi6++/fVQIA25ctY8Ojj6YXqIiSPiJ4V0Qsz0//zJA0P/+afYF3AiOAyZKGRUS8ecd8cYwHGDx4cMIxzczap5pue27r1j2FJMWX6BFBRCzP364CHgBOBZYB90fOE0ArcPBe9p0YEY0R0djQ0JBkTDOz/erzkUuo6v/GKc3uJ57AQaPOSDFR8SR2RCCpDqiIiPX5++8D/hPYAIwGHpV0FFADrN73VzIzS1/1gAEM+82vWf/wI1TU1tLz7LNQVRqnWYsvyf+KAcADkna+zl0RMV1SDfAzSU8D24BxbaeFzMxKUWXPnvS5+INpxyi6xIogIhYDJ+5l+zbg0qRe18zMDoyvLDYzyzgXgZlZxrkIzMwyzkVgZpZxLgIzs4xzEZiZZZyLwMws41QO13JJagZeTDtHAQ7GV0kXk7+fxePvZXGVy/fzbRGx3zV6yqIIyoWkpohoTDtHV+HvZ/H4e1lcXe376akhM7OMcxGYmWWci6C4JqYdoIvx97N4/L0sri71/fQ5AjOzjPMRgZlZxrkIikDSxZJC0jFpZyl3klokzZX0pKTZkk5PO1M5k3SIpHskPS/pGUm/zX8glB2gN/1szsv/fH5ZUpf4HeqpoSKQNBkYCDwSERNSjlPWJG2IiIPy988FvhYRZ6Ycqywp96lQfwbuiIhb8ttOAnpGxOOphitDbX42+wN3Af8vIv4j3WQd1yXaLE2SDgLeBXwa+KeU43Q1vYDX0g5Rxs4Gtu8sAYCImOsS6Lj857CPB76YL9yy1jU+cDNdHwSmR8RCSa9KOiUiZqcdqoz1kDQX6E7uKGt0ynnK2fHArLRDdFURsTg/NdQfeCXtPB3hI4KO+yhwT/7+PfmxvXWbI+KkiDgGOA/4eVf4i8u6rC7xs+kjgg6Q1I/cX6zHSwqgEghJ14RPvnRYRPxF0sFAA7Aq7TxlaB4wNu0QXZWkYUALXeBn00cEHTMW+HlEvC0ihkTEIGAJ8O6Uc3UJ+XdhVQJr0s5Spn4PdJN0xc4NkkZI8sn3DpLUANwC3NQV/ujzEUHHfBS4vs22+4CPAT4h99bsPEcAucPucRHRkmagchURIeli4EeSrgW2AC8AX0o1WPna+bNZDewAfgH8IN1IxeG3j5qZZZynhszMMs5FYGaWcS4CM7OMcxGYmWWci8DMLOP89lGzAkmaAGwgtwbSYxHxcDvPvQx4KCKWd046s7fORWB2gCLiGwU87TLgacBFYCXPU0Nm7ZD0dUkLJD0MHJ3fdruksfn7wyX9UdIsSQ9KGph/rBG4M79+fY8U/xPM9stFYLYPkoaTW1r8ZOBDwIg2j1cDNwJjI2I48DPgOxExBWgCPp5fQG9z5yY3OzCeGjLbtzOAByJiE4CkqW0eP5rcUs8z8gukVgIrOjWhWRG4CMza194aLALmRcTIzgpjlgRPDZnt22PAxZJ6SOoJXNDm8QVAg6SRkJsqknRc/rH1QM/Oi2r21vmIwGwfImK2pF8Cc4EXabOibERsy58Y/rGk3uT+f/oRuc8BuB24RdJmYKTPE1gp8+qjZmYZ56khM7OMcxGYmWWci8DMLONcBGZmGeciMDPLOBeBmVnGuQjMzDLORWBmlnH/H8MQ23siYljzAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sns.swarmplot(x=\"diet\", y=\"coag\", data=coagulation)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fit the model:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"OLS Regression Results \n",
"\n",
" Dep. Variable: coag R-squared: 0.671 \n",
" \n",
"\n",
" Model: OLS Adj. R-squared: 0.621 \n",
" \n",
"\n",
" Method: Least Squares F-statistic: 13.57 \n",
" \n",
"\n",
" Date: Wed, 26 Sep 2018 Prob (F-statistic): 4.66e-05 \n",
" \n",
"\n",
" Time: 10:18:05 Log-Likelihood: -52.540 \n",
" \n",
"\n",
" No. Observations: 24 AIC: 113.1 \n",
" \n",
"\n",
" Df Residuals: 20 BIC: 117.8 \n",
" \n",
"\n",
" Df Model: 3 \n",
" \n",
"\n",
" Covariance Type: nonrobust \n",
" \n",
"
\n",
"\n",
"\n",
" coef std err t P>|t| [0.025 0.975] \n",
" \n",
"\n",
" Intercept 61.0000 1.183 51.554 0.000 58.532 63.468 \n",
" \n",
"\n",
" diet[T.B] 5.0000 1.528 3.273 0.004 1.814 8.186 \n",
" \n",
"\n",
" diet[T.C] 7.0000 1.528 4.583 0.000 3.814 10.186 \n",
" \n",
"\n",
" diet[T.D] -1.776e-14 1.449 -1.23e-14 1.000 -3.023 3.023 \n",
" \n",
"
\n",
"\n",
"\n",
" Omnibus: 0.476 Durbin-Watson: 2.205 \n",
" \n",
"\n",
" Prob(Omnibus): 0.788 Jarque-Bera (JB): 0.029 \n",
" \n",
"\n",
" Skew: 0.074 Prob(JB): 0.985 \n",
" \n",
"\n",
" Kurtosis: 3.084 Cond. No. 5.79 \n",
" \n",
"
Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: coag R-squared: 0.671\n",
"Model: OLS Adj. R-squared: 0.621\n",
"Method: Least Squares F-statistic: 13.57\n",
"Date: Wed, 26 Sep 2018 Prob (F-statistic): 4.66e-05\n",
"Time: 10:18:05 Log-Likelihood: -52.540\n",
"No. Observations: 24 AIC: 113.1\n",
"Df Residuals: 20 BIC: 117.8\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 61.0000 1.183 51.554 0.000 58.532 63.468\n",
"diet[T.B] 5.0000 1.528 3.273 0.004 1.814 8.186\n",
"diet[T.C] 7.0000 1.528 4.583 0.000 3.814 10.186\n",
"diet[T.D] -1.776e-14 1.449 -1.23e-14 1.000 -3.023 3.023\n",
"==============================================================================\n",
"Omnibus: 0.476 Durbin-Watson: 2.205\n",
"Prob(Omnibus): 0.788 Jarque-Bera (JB): 0.029\n",
"Skew: 0.074 Prob(JB): 0.985\n",
"Kurtosis: 3.084 Cond. No. 5.79\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod = smf.ols(\"coag ~ diet\", coagulation).fit()\n",
"lmod.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Examine the coefficients:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 61.0\n",
"diet[T.B] 5.0\n",
"diet[T.C] 7.0\n",
"diet[T.D] -0.0\n",
"dtype: float64"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod.params.round(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Look at the design matrix:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"DesignMatrix with shape (24, 4)\n",
" Intercept diet[T.B] diet[T.C] diet[T.D]\n",
" 1 0 0 0\n",
" 1 0 0 0\n",
" 1 0 0 0\n",
" 1 0 0 0\n",
" 1 1 0 0\n",
" 1 1 0 0\n",
" 1 1 0 0\n",
" 1 1 0 0\n",
" 1 1 0 0\n",
" 1 1 0 0\n",
" 1 0 1 0\n",
" 1 0 1 0\n",
" 1 0 1 0\n",
" 1 0 1 0\n",
" 1 0 1 0\n",
" 1 0 1 0\n",
" 1 0 0 1\n",
" 1 0 0 1\n",
" 1 0 0 1\n",
" 1 0 0 1\n",
" 1 0 0 1\n",
" 1 0 0 1\n",
" 1 0 0 1\n",
" 1 0 0 1\n",
" Terms:\n",
" 'Intercept' (column 0)\n",
" 'diet' (columns 1:4)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import patsy\n",
"patsy.dmatrix('~ diet', coagulation)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Test the significance of diet."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" df \n",
" sum_sq \n",
" mean_sq \n",
" F \n",
" PR(>F) \n",
" \n",
" \n",
" \n",
" \n",
" diet \n",
" 3.0 \n",
" 228.0 \n",
" 76.0 \n",
" 13.571429 \n",
" 0.000047 \n",
" \n",
" \n",
" Residual \n",
" 20.0 \n",
" 112.0 \n",
" 5.6 \n",
" NaN \n",
" NaN \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" df sum_sq mean_sq F PR(>F)\n",
"diet 3.0 228.0 76.0 13.571429 0.000047\n",
"Residual 20.0 112.0 5.6 NaN NaN"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sm.stats.anova_lm(lmod)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fit without an intercept term:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"OLS Regression Results \n",
"\n",
" Dep. Variable: coag R-squared: 0.671 \n",
" \n",
"\n",
" Model: OLS Adj. R-squared: 0.621 \n",
" \n",
"\n",
" Method: Least Squares F-statistic: 13.57 \n",
" \n",
"\n",
" Date: Wed, 26 Sep 2018 Prob (F-statistic): 4.66e-05 \n",
" \n",
"\n",
" Time: 10:18:05 Log-Likelihood: -52.540 \n",
" \n",
"\n",
" No. Observations: 24 AIC: 113.1 \n",
" \n",
"\n",
" Df Residuals: 20 BIC: 117.8 \n",
" \n",
"\n",
" Df Model: 3 \n",
" \n",
"\n",
" Covariance Type: nonrobust \n",
" \n",
"
\n",
"\n",
"\n",
" coef std err t P>|t| [0.025 0.975] \n",
" \n",
"\n",
" diet[A] 61.0000 1.183 51.554 0.000 58.532 63.468 \n",
" \n",
"\n",
" diet[B] 66.0000 0.966 68.316 0.000 63.985 68.015 \n",
" \n",
"\n",
" diet[C] 68.0000 0.966 70.387 0.000 65.985 70.015 \n",
" \n",
"\n",
" diet[D] 61.0000 0.837 72.909 0.000 59.255 62.745 \n",
" \n",
"
\n",
"\n",
"\n",
" Omnibus: 0.476 Durbin-Watson: 2.205 \n",
" \n",
"\n",
" Prob(Omnibus): 0.788 Jarque-Bera (JB): 0.029 \n",
" \n",
"\n",
" Skew: 0.074 Prob(JB): 0.985 \n",
" \n",
"\n",
" Kurtosis: 3.084 Cond. No. 1.41 \n",
" \n",
"
Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: coag R-squared: 0.671\n",
"Model: OLS Adj. R-squared: 0.621\n",
"Method: Least Squares F-statistic: 13.57\n",
"Date: Wed, 26 Sep 2018 Prob (F-statistic): 4.66e-05\n",
"Time: 10:18:05 Log-Likelihood: -52.540\n",
"No. Observations: 24 AIC: 113.1\n",
"Df Residuals: 20 BIC: 117.8\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"diet[A] 61.0000 1.183 51.554 0.000 58.532 63.468\n",
"diet[B] 66.0000 0.966 68.316 0.000 63.985 68.015\n",
"diet[C] 68.0000 0.966 70.387 0.000 65.985 70.015\n",
"diet[D] 61.0000 0.837 72.909 0.000 59.255 62.745\n",
"==============================================================================\n",
"Omnibus: 0.476 Durbin-Watson: 2.205\n",
"Prob(Omnibus): 0.788 Jarque-Bera (JB): 0.029\n",
"Skew: 0.074 Prob(JB): 0.985\n",
"Kurtosis: 3.084 Cond. No. 1.41\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmodi = smf.ols(\"coag ~ diet-1\", coagulation).fit()\n",
"lmodi.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Do the ANOVA explicitly:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" df_resid \n",
" ssr \n",
" df_diff \n",
" ss_diff \n",
" F \n",
" Pr(>F) \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 23.0 \n",
" 340.0 \n",
" 0.0 \n",
" NaN \n",
" NaN \n",
" NaN \n",
" \n",
" \n",
" 1 \n",
" 20.0 \n",
" 112.0 \n",
" 3.0 \n",
" 228.0 \n",
" 13.571429 \n",
" 0.000047 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" df_resid ssr df_diff ss_diff F Pr(>F)\n",
"0 23.0 340.0 0.0 NaN NaN NaN\n",
"1 20.0 112.0 3.0 228.0 13.571429 0.000047"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%%capture --no-display\n",
"lmodnull = smf.ols(\"coag ~ 1\", coagulation).fit()\n",
"sm.stats.anova_lm(lmodnull, lmod)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Do it with sum contrasts:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"OLS Regression Results \n",
"\n",
" Dep. Variable: coag R-squared: 0.671 \n",
" \n",
"\n",
" Model: OLS Adj. R-squared: 0.621 \n",
" \n",
"\n",
" Method: Least Squares F-statistic: 13.57 \n",
" \n",
"\n",
" Date: Wed, 26 Sep 2018 Prob (F-statistic): 4.66e-05 \n",
" \n",
"\n",
" Time: 10:18:05 Log-Likelihood: -52.540 \n",
" \n",
"\n",
" No. Observations: 24 AIC: 113.1 \n",
" \n",
"\n",
" Df Residuals: 20 BIC: 117.8 \n",
" \n",
"\n",
" Df Model: 3 \n",
" \n",
"\n",
" Covariance Type: nonrobust \n",
" \n",
"
\n",
"\n",
"\n",
" coef std err t P>|t| [0.025 0.975] \n",
" \n",
"\n",
" Intercept 64.0000 0.498 128.537 0.000 62.961 65.039 \n",
" \n",
"\n",
" C(diet, Sum)[S.A] -3.0000 0.974 -3.081 0.006 -5.031 -0.969 \n",
" \n",
"\n",
" C(diet, Sum)[S.B] 2.0000 0.845 2.366 0.028 0.237 3.763 \n",
" \n",
"\n",
" C(diet, Sum)[S.C] 4.0000 0.845 4.732 0.000 2.237 5.763 \n",
" \n",
"
\n",
"\n",
"\n",
" Omnibus: 0.476 Durbin-Watson: 2.205 \n",
" \n",
"\n",
" Prob(Omnibus): 0.788 Jarque-Bera (JB): 0.029 \n",
" \n",
"\n",
" Skew: 0.074 Prob(JB): 0.985 \n",
" \n",
"\n",
" Kurtosis: 3.084 Cond. No. 2.68 \n",
" \n",
"
Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: coag R-squared: 0.671\n",
"Model: OLS Adj. R-squared: 0.621\n",
"Method: Least Squares F-statistic: 13.57\n",
"Date: Wed, 26 Sep 2018 Prob (F-statistic): 4.66e-05\n",
"Time: 10:18:05 Log-Likelihood: -52.540\n",
"No. Observations: 24 AIC: 113.1\n",
"Df Residuals: 20 BIC: 117.8\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"=====================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"-------------------------------------------------------------------------------------\n",
"Intercept 64.0000 0.498 128.537 0.000 62.961 65.039\n",
"C(diet, Sum)[S.A] -3.0000 0.974 -3.081 0.006 -5.031 -0.969\n",
"C(diet, Sum)[S.B] 2.0000 0.845 2.366 0.028 0.237 3.763\n",
"C(diet, Sum)[S.C] 4.0000 0.845 4.732 0.000 2.237 5.763\n",
"==============================================================================\n",
"Omnibus: 0.476 Durbin-Watson: 2.205\n",
"Prob(Omnibus): 0.788 Jarque-Bera (JB): 0.029\n",
"Skew: 0.074 Prob(JB): 0.985\n",
"Kurtosis: 3.084 Cond. No. 2.68\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from patsy.contrasts import Sum\n",
"lmods = smf.ols(\"coag ~ C(diet,Sum)\", coagulation).fit()\n",
"lmods.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Diagnostics"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGa9JREFUeJzt3X+U3XV95/HneyYTkoFgsofMohIMsdVs2bVIs1Yrxoj4Y9UNslstHt1Ta9kbS93YbEWbsrCAh9YsRX6UFcKJlhY5sOIWDRY1YEyzlsIxgYCNBsQQTJA40CYamAmZH+/94/vNZBIyk5ubufOdH8/HOXPufL/3zv2+zjcz95Xvz09kJpIktVQdQJI0NlgIkiTAQpAklSwESRJgIUiSShaCJAmwECRJJQtBkgRYCJKk0pSqAxyNk046KefOnVt1DEkaVzZu3PhcZs4+0uvGVSHMnTuXDRs2VB1DksaViHiqnte5y0iSBFgIkqSShSBJAiwESVLJQpAkAePsLCNJY8u6LZ2sXL+V7bu6mDOrnSUL57FofkfVsdQgtxAkNWTdlk4uXb2Zzj17mTm9jc49e7l09WbWbemsOpoaZCFIasjK9Vtpaw3ap04honhsaw1Wrt9adTQ1yEKQ1JDtu7qY3tZ60Lzpba3s2NVVUSIdKwtBUkPmzGqnu6fvoHndPX2cMqu9okQ6VhaCpIYsWTiPnr6ka18vmcVjT1+yZOG8qqOpQRaCpIYsmt/BFYtPp2PGNH7R3UPHjGlcsfh0zzIaxzztVFLDFs3vsAAmELcQJEmAhSBJKlkIkiTAQpAklSwESRJgIUiSShaCJAmwECRJJQtBkgRYCJKkkoUgSQLGQCFERGtEPBwR36g6iyRNZpUXAvBJ4EdVh5Ckya7SQoiIU4D3AquqzCFJqn4L4Vrg00B/xTkkadKrrBAi4n1AZ2ZuPMLrahGxISI2PPvss6OUTpImnyq3EN4MLI6IbcAdwNkR8eVDX5SZN2fmgsxcMHv27NHOKEmTRmWFkJnLM/OUzJwLnA+szcyPVJVHkia7qo8hSJLGiDExpnJmrgPWVRxDkiY1txAkSYCFIEkqWQiSJMBCkCSVLARJEmAhSJJKFoIkCbAQJEklC0GSBFgIkqSShSBJAiwESVJpTNzcbrSs29LJyvVb2b6rizmz2lmycB6L5ndUHUuSDlLVZ9Wk2UJYt6WTS1dvpnPPXmZOb6Nzz14uXb2ZdVs6q44mSQOq/KyaNIWwcv1W2lqD9qlTiCge21qDleu3Vh1NkgZU+Vk1aQph+64upre1HjRvelsrO3Z1VZRIkl6qys+qSVMIc2a1093Td9C87p4+TpnVXlEiSXqpKj+rJk0hLFk4j56+pGtfL5nFY09fsmThvKqjSdKAKj+rJk0hLJrfwRWLT6djxjR+0d1Dx4xpXLH4dM8ykjSmVPlZFZnZ9IWMlAULFuSGDRuqjiFJ40pEbMzMBUd63aTZQpAkDc9CkCQBFoIkqWQhSJIAC0GSVLIQJEmAhSBJKlkIkiTAQpAklSwESRJgIUiSSpUVQkTMiYjvRsSPImJzRHyyqiySpGrHVO4F/jgzH4qIGcDGiLg3M3/YrAUuu+MhVj+6k77+pLUlWPy6k7nm/DObtThJGlcq20LIzGcy86Hy+z3Aj4BXNmt5y+54iLs2PUNff3F3177+5K5Nz7DsjoeatUhJGlfGxDGEiJgLvB54sFnLWP3oznJZB74Gz5ekya7yQoiIE4D/C/xRZv7yMM/XImJDRGx49tlnG17O/i2DeudL0mRTaSFERBtFGdyWmX97uNdk5s2ZuSAzF8yePbvhZbW2xFHNl6TJpsqzjAL4IvCjzPx8s5e3+HUnA5B54GvwfEma7KrcQngz8F+AsyNiU/n1nmYt7Jrzz+S8M14+sEXQ2hKcd8bLPctIkkqOqSxJE5xjKkuSjoqFIEkCLARJUslCkCQBFoIkqWQhSJIAC0GSVLIQJEmAhSBJKlkIkiTAQpAklaocQnPUXX/f46z63pO8sK+P46e2csFZp7H0nNdUHUsat9Zt6WTl+q1s39XFnFntLFk4j0XzO6qOpQbVtYUQEW+OiOPL7z8SEZ+PiFc1N9rIuv6+x7lu7RN09/QxpQW6e/q4bu0TXH/f41VHk8aldVs6uXT1Zjr37GXm9DY69+zl0tWbWbels+poalC9u4xuBLoi4teBTwNPAX/TtFRNsOp7T9ISMKWlhZZoKR+L+ZKO3sr1W2lrDdqnTiGieGxrDVau31p1NDWo3kLozeI+2ecC12XmdcCM5sUaeS/s6+PQwdFaopgv6eht39XF9LbWg+ZNb2tlx66uihLpWNVbCHsiYjnwEeDvIqIVaGterJF3/NRWDh0+uT+L+ZKO3pxZ7XT3HPwfqu6ePk6Z1V5RIh2regvhd4AXgd/PzJ3AK4GrmpaqCS446zT6E3r7++nP/vKxmC/p6C1ZOI+evqRrXy+ZxWNPX7Jk4byqo6lBdZ1lVJbA5wdN/5Rxdgxh/9lEnmUkjYxF8zu4guJYwo5dXZziWUbj3rBDaEbEHuBwLwggM/PEZgU7HIfQlKSjV+8QmsNuIWTmuDpwLElq3FFdmBYRHcC0/dPlriNJ0gRQ74VpiyPix8CTwN8D24BvNjGXJGmU1XuW0WeBNwKPZ+ZpwNuBf2haKknSqKu3EHoy85+BlohoyczvAmc0MZckaZTVewxhd0ScAKwHbouITqC3ebEkSaOt3i2Ec4FuYBnwLeAnwH9sVihJ0uir98K0FwZN/nWTskiSKlRXIRxygdpUivsYvTDaF6ZJkpqn3i2Egy5Qi4j3A29oSiJJUiUaGkIzM78GnD3CWSRJFap3l9F/GjTZAizg8Pc4OioR8W7gOqAVWJWZnzvW95QkNabe004Hn1HUS3Gl8rnHsuByTIX/DbwD2AF8PyJWZ+YPj+V9h+P4r9LI8m+qOZbd8RCrH91JX3/S2hIsft3JXHP+mU1fbr3HEH6vCct+A/BEZm4FiIg7KEqmKYWwf/zXttY4aPzXK8BfYKkB/k01x7I7HuKuTc8MTPf1Zzn9UNNLYdhjCBHxlxFx/VBfx7jsVwLbB03vKOcN6amnnuLuu+8GoLe3l1qtxj333APA3r17qdVqrFmzBoDnn3+eWq3G2rVrAfjCmkd5cd2NtPx8CxHBcX1dvLjuRj5/2zcA2LlzJ7VajQcffLAIs2MHtVqNjRs3ArBt2zZqtRqPPPIIAE888QS1Wo3NmzcD8Nhjj1Gr1XjssccA2Lx5M7VajSeeeAKARx55hFqtxrZt2wDYuHEjtVqNHTt2APDggw9Sq9XYuXMnAPfffz+1Wo3nnnsOgPXr11Or1di9ezcAa9eupVar8fzzzwOwZs0aarUae/fuBeCee+6hVqvR21tcP3j33XdTq9UG1uVdd93FhRdeODB95513snTp0oHp22+/nWXLlg1M33rrrVx00UUD07fccgvLly8fmF61ahWXXHLJwPRNN93E5ZdfPjB9ww03cOWVVw5MX3vttaxYsWJg+uqrr+bqq68emF6xYgXXXnvtwPSVV17JDTfcMDB9+eWXc9NNNw1MX3LJJaxatWpgevny5dxyyy0D0xdddBG33nrrwPSyZcu4/fbbB6aXLl3KnXfeOTB94YUXctdddw1M12q1hn/3du/eTa1WY/369QA899xz1Go17r//fmD8/u594b6iDOLpR/mXb9/A9JZ+2lqDq754h797x/C7t/ZLn+O4pzcQARHwso1/xXE/e5jVj+5s+HevXkc6qLwB2Ehxh9MzgR+XX2cAxzoYcRxm3kuOS0RELSI2RMSGnp6ehhf29O5uWuLgRbZE8NzzLzb8ntJk9vTuw4+p/C8v7Kso0cQw1MHZvkPHAG6CYQfIGXhRxHeBd2ZmTzndBqzJzLc1vOCINwGXZea7yunlAJn550P9zLEMkPOhmx+gc89e2qce2EvWta+XjhnTuL32xobeU5rM/Jtqjlf/6T309SeD//+aCa0twU/+7D0NvWe9A+TUe9rpK4DB1yKcUM47Ft8HfjUiTouIqcD5wOpjfM8hOf6rNLL8m2qOxa87GShKYP/X4PnNVO9ZRp8DHi63FADeClx2LAvOzN6I+ATwbYrTTr+UmZuP5T2H4/iv0sjyb6o5igPH1ZxlVNcuI4CIOBn4zXLywczc2bRUQ3BMZUk6eiOyyygi5pePZ1LsItpefr2inCdJmiCOtMvovwM14OrDPJd4+wpJmjCGLYTMrJWPDZ9NJEkaH+o6yygiPhARM8rv/0dE/G1EvL650SRJo6ne004vycw9EXEW8C6KQXJuOsLPSJLGkXoLYf9Vye8FbszMr1MMlCNJmiDqLYSnI2Il8EHgnog47ih+VpI0DtT7of5BigvI3p2Zu4F/BVw0/I9IksaTugohM7uATuCsclYvxU3uJEkTRL1nGf1P4DPA/nvOtgFfblYoSdLoq3eX0XnAYuAFgMz8GQff7E6SNM7VWwj7srjpUQJExPHNiyRJqkK9dzv9SnmW0cyI+K/Ax4BVR/iZMef6+x5n1fee5IV9fRw/tZULzjqNpee8pupY0rjlmMrNUdV6PZq7nb4DeCfFSGffzsx7mxnscI7lbqfX3/c41619gpaAloD+LL4+efavWApSAwaPqTy9rZXunj56+pIrFp9uKRyDZqzXkR4gh8y8NzMvysxPAWsj4sMNJavIqu89SUvAlJYWWqKlfCzmSzp6K9dvpa01aJ86hYjisa01WLl+a9XRxrUq1+uRbn99YkQsj4gbIuKdUfgEsJXi2oRx44V9fbQcMopzSxTzJR297bsOP6byjl1dFSWaGKpcr0faQrgVeC3wA+ACYA3wAeDczDy3ydlG1PFTWzl0jOr+LOZLOnpzZrXT3XPwf6i6e/o4ZVZ7RYkmhirX65EKYV5mfjQzVwIfAhYA78vMTU1PNsIuOOs0+hN6+/vpz/7ysZgv6eg5pnJzVLlej1QIPfu/ycw+4MnM3NPcSM2x9JzX8Mmzf4Xpba309hebYB5Qlhq3aH4HVyw+nY4Z0/hFdw8dM6Z5QHkEVLlehz3LKCL6KC9Gozi7aDrQVX6fmXli0xMO4pjKknT06j3L6EgjprmDXZImCW9hLUkCLARJUslCkCQBFoIkqWQhSJIAC0GSVLIQJEmAhSBJKlkIkiTAQpAklSophIi4KiK2RMSjEXFXRMysIock6YB6x1QeafcCyzOzNyJWAMuBzzR7oY6pLI0sx1RujqrWayVbCJm5JjN7y8kHgFOavcz9Yyp39/QxpaUYcOK6tU9w/X2PN3vR0oS0f+zfzj17mTm9jc49e7l09WbWbemsOtq4VuV6HQvHED4GfLPZC3FMZWlkOaZyc1S5Xpu2yygi7gNOPsxTF2fm18vXXAz0ArcN8z41oAZw6qmnNpznhX3FlsFgjqksNW77ri5mTm87aJ5jKh+7Ktdr0wohM88Z7vmI+F3gfcDbc5hRejLzZuBmKAbIaTTP8VNb6e7poyUOzHNMZalxc2a107lnL+1TD3yMOKbysatyvVZ1ltG7KQ4iL87MUfnvhGMqSyPLMZWbYyyPqdwsNwAzgHsjYlNE3NTsBTqmsjSyHFO5OcbsmMpjjWMqS9LRq3dM5bFwlpEkaQywECRJgIUgSSpZCJIkwEKQJJUsBEkSYCFIkkoWgiQJsBAkSSULQZIEWAiSpFJVQ2hKmgAcQnNicQtBUkMcQnPisRAkNcQhNCceC0FSQ7bv6mJ628EjDjqE5vhmIUhqyJxZ7XT3HDwmuUNojm8WgqSGOITmxGMhSGqIQ2hOPJ52Kqlhi+Z3WAATiFsIkiTAQpAklSwESRJgIUiSShaCJAmwECRJJQtBkgRYCJKkkoUgSQIsBElSyUKQJAEVF0JEfCoiMiJOqjKHJKnCQoiIOcA7gJ9WlUGSdECVWwjXAJ8GssIMkqRSJYUQEYuBpzPzkSqWL0l6qaaNhxAR9wEnH+api4E/Bd5Z5/vUgBrAqaeeOmL5JEkHi8zR3WMTEf8O+A6wfyTuU4CfAW/IzJ3D/eyCBQtyw4YNTU4oSRNLRGzMzAVHet2oj5iWmT8ABoZYiohtwILMfG60s0iSDvA6BEkSMAbGVM7MuVVnkCS5hSBJKlkIkiTAQpAklSwESRJgIUiSShaCJAmwECRJJQtBkgRYCJKkkoUgSQIsBElSyUKQJAFj4OZ2o2ndlk5Wrt/K9l1dzJnVzpKF81g0v+PIPyhJk8Ck2UJYt6WTS1dvpnPPXmZOb6Nzz14uXb2ZdVs6q44mSWPCpCmEleu30tYatE+dQkTx2NYarFy/tepokjQmTJpC2L6ri+ltrQfNm97Wyo5dXUP8hCRNLpOmEObMaqe7p++ged09fZwyq72iRJI0tkyaQliycB49fUnXvl4yi8eevmTJwnlVR5OkMWHSFMKi+R1csfh0OmZM4xfdPXTMmMYVi0/3LCNJKk2q004Xze+wACRpCJNmC0GSNDwLQZIEWAiSpJKFIEkCLARJUikys+oMdYuIZ4GnRuCtTgKeG4H3GQ3jJet4yQlmbYbxkhMmZ9ZXZebsI71oXBXCSImIDZm5oOoc9RgvWcdLTjBrM4yXnGDW4bjLSJIEWAiSpNJkLYSbqw5wFMZL1vGSE8zaDOMlJ5h1SJPyGIIk6aUm6xaCJOkQE74QImJmRHw1IrZExI8i4k0R8YGI2BwR/RExZs42GCLrVeX0oxFxV0TMrDonDJn1s2XOTRGxJiJeMRZzDnruUxGREXFSlRn3G2KdXhYRT5frdFNEvKfqnDD0eo2I/xYRj5V/X/+r6pww5Hr9P4PW6baI2DRGc54REQ+UOTdExBuaGiIzJ/QX8NfABeX3U4GZwL8BXgusAxZUnfEIWd8JTCnnrQBWVJ1zmKwnDnp+KXDTWMxZfj8H+DbFdS0nVZ1zmHV6GfCpqrPVmfVtwH3AceX8jqpzDvc7MOj5q4FLx2JOYA3wH8p57wHWNTPDhL79dUScCCwEPgqQmfuAfcDu8vnKsh1qmKxrBr3sAeC3Rz3cIYbJOtjxQKUHqI6Q8xrg08DXKwl3iKGyjqXf0f2GyfoHwOcy88VyfmdlIUtH+l2NYgV/EDi7inyDcgy1ThM4sXzZy4CfNTPHRN9lNA94FviriHg4IlZFxPFVhxpCPVk/Bnxz9KO9xJBZI+LKiNgOfBi4tMqQDJEzIhYDT2fmIxXnG2y4f/9PlLvivhQRsyrMuN9QWV8DvCUiHoyIv4+If19tTODIf1dvAX6emT+uJt6AoXL+EXBV+Tf1F8DyZoaY6IUwBTgTuDEzXw+8APxJtZGGNGzWiLgY6AVuqybeQYbMmpkXZ+YcipyfqC4icPiclwEXU31ZHWqodXoj8GrgDOAZit0bVRsq6xRgFvBG4CLgK1H9Js6RPgM+BNxeRbBDDJXzD4Bl5d/UMuCLTU1R9X6zJu+TOxnYNmj6LcDfDZpexxg5hjBcVuB3gX8E2qvOWc96Lee9CvinMZjzO0AnsK386gV+Cpw8BrMeuk7nVr1Oh8sKfAtYNGj+T4DZYzFr+f0U4OfAKWN4nf6CA5cHBPDLZuaY0FsImbkT2B4Rry1nvR34YYWRhjRU1oh4N/AZYHFmdlUWcJBhsv7qoJctBraMerhBhsj5UGZ2ZObczJwL7ADOLF9bmWHW6csHvew84J9GPdwhhvm7+hrlvviIeA3FgdFKbyJ3hM+Ac4AtmbmjknCDDJPzZ8Bby3lnA03dtTXhL0yLiDOAVRS/nFuB3wMWAX8JzKY4wLwpM99VVcb9hsj6feA44J/Llz2QmR+vJuEBQ2RdRXH2Vj/F2Tsfz8ynKwvJ4XNm5q5Bz2+j2Eqs/O6XQ6zT6yl2FyXFFs2SzHymqoz7DZH1BeBLFHn3UZwdtbaykKWhfgci4haKv6ebqsy33xDr9HTgOoqtmb3AhZm5sWkZJnohSJLqM6F3GUmS6mchSJIAC0GSVLIQJEmAhSBJKlkImjAiom/QHSw3RcTciFgQEdeXzy+KiN8a9Pr3R8SvNbCc50co74i8jzRSJvTN7TTpdGfmGYfM2wZsKL9fBDwP3F9Ovx/4BmP0YkVptLmFoAmt3Cr4RkTMBT4OLCu3Ht5KcTX1VeX0q8uvb0XExoj4fxExv3yP0yLiHyPi+xHx2SGWsyIiLhw0fVlE/HFEnBAR34mIhyLiBxFx7lAZB03fEBEfLb//jfJGcRsj4tv7r1yOiKUR8cPypnd3jNgK06TmFoImkulxYKCTJzPzvP1PZOa2iLgJeD4z/wIgIlYD38jMr5bT36G4uvrHEfGbwBcobhdwHcVNx/4mIv5wiGXfAVxb/gwUt1R+N8XVpedl5i+jGIjngYhYnXVcERoRbRRX1J+bmc9GxO8AV1Lc9fZPgNMy88UYI4MmafyzEDSRHG6XUV0i4gTgt4A7B92g87jy8c3Afy6/v5VioKKDZObDEdERxShxs4FdmfnT8kP9zyJiIcUtPV4J/GugnnsnvRb4t8C9ZaZWijueAjwK3BYRX6O4h5B0zCwEqdAC7B6mUOq5x8tXKQYwOpliiwGKcSFmA7+RmT3lvZOmHfJzvRy8+3b/8wFszsw38VLvpRhQZTFwSUScnpm9dWSUhuQxBE0me4AZh5vOzF8CT0bEB6AYSSsifr183T8A55fff3iY97+jfN1vU5QDFKNcdZZl8DaK24If6ing1yLiuIh4GcWdLgEeA2bHgfGK2yLi9IhoAeZk5ncpRn2bCZxQ1xqQhmEhaDK5GzivPIj8FooP8IvKEapeTfFh//sR8QiwGdh/APiTwB9GxPcpPuAPKzM3UxTM04PuSHobsCAiNpTv/5JbgmfmduArlLuBgIfL+fsoymVFmWkTxW6tVuDLEfGD8rXXZObuRleKtJ93O5UkAW4hSJJKFoIkCbAQJEklC0GSBFgIkqSShSBJAiwESVLJQpAkAfD/AVsEOOTRT/jKAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sns.residplot(lmod.fittedvalues, lmod.resid)\n",
"plt.xlabel(\"Fitted values\")\n",
"plt.ylabel(\"Residuals\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/anaconda/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYVOWVx/HvAVRoRNCACmrTigtqFBBc47iiccWoGBVEUEOPmcRlJgnqkEmcKHGNUWM0KQQBaQFRxyBxN24xxgACIuOAxLAJRgRBtBFZzvzx3pam6a663V1Vt6r793meeqrq1q17T5dYp859732PuTsiIiItkg5AREQKgxKCiIgASggiIhJRQhAREUAJQUREIkoIIiICKCGIiEhECUFERAAlBBERibRKOoD66Nixo5eVlSUdhohIUZkxY8Yn7t4p03pFlRDKysqYPn160mGIiBQVM1sUZz0dMhIREUAJQUREIkoIIiICKCGIiEhECUFERAAlBBGRolZRAWVl0KJFuK+oaPi2iuq0UxER2aKiAsrLobIyPF+0KDwHGDiw/ttThSAiUqSGD9+SDKpUVoblDaGEICJSpBYvrt/yTJQQRESKVGlp/ZZnooQgIlKkRoyAkpKtl5WUhOUNoYQgIlKkBg6EVAq6dgWzcJ9KNWxAGXSWkYhIURs4sOEJoCZVCCIiAighiIhIRAlBREQAJQQREYkoIYiICKCEICIiESUEEREBlBBERCSihCAiIoASgoiIRJQQREQEKICEYGYtzWymmU1NOhYRkeYs8YQAXAO8l3QQIiLNXaIJwcz2BM4EHkwyDhERSb5CuBsYBmxOOA4RkWYvsYRgZmcBH7v7jAzrlZvZdDObvmLFijxFJyLS/CRZIXwL6GdmC4GJwElmNr7mSu6ecvc+7t6nU6dO+Y5RRKTZSCwhuPsN7r6nu5cBFwF/cvdLkopHRKS5S3oMQURECkRB9FR291eAVxIOQ0SkWVOFICIigBKCiIhElBBERARQQhARkYgSgoiIAEoIIiISUUIQERFACUFERCJKCCIiTYF7ozehhCAiUszmz4dhw+CQQ2D9+kZtSglBRKTYrF8PEyfCSSfBAQfAr38d7j/9tFGbLYi5jEREJIZ582DkSBgzBlauhH32gVtugSFDYPfdG715VQgi0mxUVEBZGbRoEe4rKgp7uwB8+SU88giccAJ07w733AMnngjPPw/vvw/XX5+VZACqEESkmaiogPJyqKwMzxctCs8BBg4svO3y3nuhGhg7Flatgm7d4NZbQzWw226N2HDdzLMwMp0vffr08enTpycdhogUobKy8GVdU9eusHBhgWx33Tp4/HFIpeD112G77eDcc0OGOfHEUII0gJnNcPc+mdZThSAizcLixfVbntftzp0bqoFx48LA8L77wu23w+DBsOuujQuwHpQQRKRZKC2t/Zd8aWlC2123DiZPDtXAG2/A9tvDeeeFauCEE8CscYE1gAaVRaRZGDECSkq2XlZSEpbndbtz5sDVV0OXLqECWLEC7rwTli6FCRPCoaEEkgGoQhCRZqJqgHf48HA4p7Q0fGk3auA37nYrK+HRR0M18OaboRro3z9UA8cdl1gCqEmDyiIiuTJ7dhgbePhh+OyzcNpoeTkMGgQdO+YtDA0qi4gk4YsvYNKkUA289RbssANccEFIBMceWzDVQG2UEEREsmHWrJAExo+HtWvhwAPh7rtDNbDLLklHF4sSgohIQ33+eZhTKJWCadOgdWv47ndDNXDMMQVdDdRGCUFEpL7efjskgYqKkBQOPjhMKTFoEOy8c9LRNZgSgohIHGvXhtNCUymYMQPatIELLwzVwFFHFV01UBslBBGRuriHL/9UKkww98UXoe/AffeF80o7dEg6wqxSQhARqemzz0ICSKVg5sxwpdlFF4Vq4IgjmkQ1UBslBBERCNXAtGkhCUyYEC4m69ED7r8fBgyA9u2TjjDnlBBEpHlbsyYMDqdS4UKytm1DAigvhz59mmw1UBslBBFpftzDRWOpVDhtdN066NkTHnggJIOddko6wkQoIYhI87F6dbhwLJUKk8y1bQuXXBKqgd69m1U1UJuMCcHMugFL3X29mZ0AHAqMc/fVjdmxme0FjAN2BzYDKXe/pzHbFBHZhnuYUC6VClNKfPklHHYY/P73cPHF0K5d0hEWjDjTXz8ObDKzfYFRwN7AI1nY90bgR+5+IHAU8AMzOygL2xWRPGpIP+GG9iCu1/tWrYJ77w2niX7rW6ET2eDBMH16OJW0vFzJoCZ3T3sD3o7ufwJcFT2emel99b0BfwBOSbdO7969XUQKx/jx7iUl7uFneLiVlITl2XxP7Pdt3uz++uvugwa5t24dVjr8cPeRI93Xrs3K31yMgOke43s44/TXZvYWcDcwHDjb3f9hZu+6+zezlZTMrAx4Dfimu39W13qa/lqksDSkn3BDexCnfd+MlWGK6VQqNKdv1y6MDQwdCr16Zfw7mrpsTn99GXAlMCJKBnsD4xsbYBUz25FwWOra2pKBmZUD5QClje11JyJZ1ZB+wg3tQbzt686/8Drli1Kwx2Owfn2YQmL06DDBXNu26Tco24jVIMfM2gCl7j4vqzs32w6YCjzn7ndlWl8VgkhhSaJC+AafcCnjKCdFd+bxmbVnpx8MCtXAoYfW+29oDuJWCBkHlc3sbGAW8Gz0vKeZTclCgEYYpH4vTjIQkcLTkD7FDept7M7Iga8wqeUAPmQP7uJHrOQblG8/hqcfXAa/+Y2SQTZkGmQAZgDtqTaQDMyJM0CRYbvHAg68Q0g4s4Az0r1Hg8oihWf8ePeuXd3Nwn2mweF6vefjj93vuMN9v/3cwb8s6eCj213th/BO7H1JlgeV3f1IM5vp7r2iZe+4e97TsQ4ZiTQDmzfDK6+EAeInnoANG0LryfLy0Ji+TZukIyw62RxUftfMBgAtzWw/4GrgL40NUERkKx9/DGPGhKb0CxaERjM/+EEYGzhIlyjlQ5wL064CDgbWAxOAz4BrcxmUiDQTmzfDCy+EJvR77gnXXQdduoTpJZYtg1//WskgjzJWCO5eSbgGYXjuwxGRZuGjj7ZUAx98EJrQX3UVfO97oTm9JKLOhGBmTxEGfWvl7v1yEpGINE1V1UAqBVOmwMaNcPzxcPPNcO65oUG9JCpdhXBn3qIQkaZr+fJwsdiDD4YLDTp2hGuvDdXAAQckHZ1UU2dCcPdX8xmIiDQhmzbB88+HauCpp8Lzk06CW2+F73wHdtgh6QilFukOGT3q7t81sznUcugoidNORaTAffjhlmpg8WLo1Al+9KNQDey3X9LRSQbpDhldE92flY9ARKRIbdoEzz4bBoinTg3P+/aFO++Ec86B7bdPOkKJKd0ho+XRw39z9+uqv2ZmtwHXbfsuEWk2li6FUaPCbckS2HVX+MlPQjXQrVvS0UkDxLkO4ZRalp2e7UBEpAhs3BjGBM4+O8xGd+ON4TTRxx4LSeGWW5QMili6MYTvA/8G7GNm71R7qR3wRq4DE5ECsnjxlmrgww9h993h+uvhiitgn32Sjk6yJF2F8AhwNjAluq+69Xb3S/IQm4jEELetZH3aT1ZUQLeuGznHpvBSyVlsLtsbbroptKN84omQIEaMUDJoYtKNIawB1gAXm1lLYLdo/R3NbEd3z9DOQkRyraIizPlWWRmeL1oUngMMHFj/9QCevGcRi388itc2jmIPlrFsXWdub/WfdL/jCr5zbVlO/x5JVpzZTn8I3Aj8E9gcLXbNdiqSvLjNZjKut2ED/PGPkEqx+ZlnAXiG0xnJUKZyFptolbGBjRSubM52ei1wgLuvbHxYIpJNcdtR1rVei0X/gJ+OCtcOLF8Oe+zBzfwXo7icxXSNtQ1pOuKcZbSEcOhIRApMXW3Gay6v/rwVGziXJ3iG01hAt3BmUO/eYX6hhQsZ3fW/t0kG6fYlTUechPAB8IqZ3WBm/1F1y3VgIpJZ3HaUI0bAQa0/YAT/yWJKeYLz+abN5d1zfxaOA1WdStqqVcNaXEqTEOeQ0eLotn10E5ECUTUgPHx4OKRTWhq+uL8eKP7qK5gyhYFjUwz88gU20YI/ciZPdiqn752nMeDSbb8CMm5TmqyMg8qFRIPKIjEtWBDmE3roodCJrLQ0XEF82WWhEY00K1kbVDazTsAwQte0rycsd/eTGhWhiGTX+vXw5JNhTqGXXoKWLcNhoPJyOPXU8FwkjTiHjCqASYRJ7q4EBgMrchmUiNTD/PmhGhgzBlasCOeS3nxzqAa6dEk6OikicRLCN9x9lJldE/VIeNXM1CtBJEnr14crhkeOhJdfDr/++/WDf/3XMNOoqgFpgDgJYUN0v9zMzgSWAToIKZKEefNC05mxY2HlSth7b/jlL2HIEOjcOenopMjFSQg3m1l74EfAb4CdgH/PaVQissWXX4ZqIJWCV1+FVq1C17Hycjj55DA5kUgWZEwI7j41ergGODG34YjI1957LxwSGjsWVq0KE8ndemuoBnbbLenopAmKc5bRQ9TeQvPynEQk0pytWxd6C4wcCa+/Dtttt6UaOOkkVQOSU3EOGU2t9rg1cC5hHEFEsmXu3JAExo2DTz+FffeF22+HwYNDJzKRPIhzyOjx6s/NbALwYs4iEmkuKitDNZBKwRtvhGrg/PNDNXD88aoGJO/iVAg17QdomiuRhpozJ1QDDz8Mq1fD/vuHhvSXXgqdOiUdnTRjccYQ1hLGECy6/wi4LsdxiTQtlZXw6KOhGnjzTdh+e+jfH4YODdWAWdIRisQ6ZNQuVzs3s9OAe4CWwIPufmuu9iWSiHfeCUlg/HhYswa6d4e77oJBg6Bjx6SjE9lK2oOUZtbGzL5nZndFtwFmlpUZT6O2nL8FTgcOIrTqPCgb2xapr2z2JZ446gt+0nE0f7WjoEcPNqUeDHMKvfYa/O//wr//OxXPdcx6H2SRRnP3Wm/AIcBCYCxwNXBN9Hga0AG4ua73xrkBRwPPVXt+A3BDuvf07t3bRbJt/Hj3khJ32HIrKQnL67XezJk+7+Tv+xrauYO/y0F+NXf7Hm1WbrWtrO1PJCZgusf5Xq7zBXgZOKWW5X2BpdW/zBtyA/oTDhNVPR8E3JfuPUoIkgtdu279pVt169o183ptWevDdhnpfvjh7uDrrLWP4VI/hj87bK51W43ZX23riWQSNyGkG0Po7O4v1FJRvGhmGwjXIzRGbaNo21wAZ2blQDlAqXr4SQ40pC/xYcygnBQDeIR2qz6HLt+Ee++ly9WX8Ck7p91WY/sgq7ex5Eq6MYQWZrZDzYVm1hrY4O6Vjdz3UmCvas/3pJYL3tw95e593L1PJ52SJzkQty/xgXuuZSgpptObGfRhEA/zGP05b/e/hMHjq65ip67bJoOa22pIH+Q4y0UaK11CGAc8bmZlVQuix48CD2dh39OA/cxs72ig+iJgSha2K1IvaXsIu8P06VBezqyPO5PiX9mODfyQ39CFZfyw5CHOv/Por08bjdOPuD59kNXbWPIq3fEk4IeEfsqfRLdFwFVxjkXFuQFnAPOBvwPDM62vMQTJlfHjw7F5s3A/MbXG/YEH3Hv12jKae/nl/syNf/WupZu/Xq+2Ad6a22roOvVZTyQdYo4hxOqpbGbtouSxNhdJKS71VJaccodp08J1AxMmhIvJevQITWcGDID27ZOOUKRBstZTGZJPBCI5tWZNOME/lYLZs6Ft25AAysuhTx9dRSzNRkPmMhIpfu7w1lshCUycGKad7tULHnggJIOddko6QpG8U0KQ5mX16jCNRCoVJpnbcccwjUR5OfTunXR0IomKM7ldCaF9Zqm7DzWz/YADfEsnNZHC5h4mlEulYNKk0JKyd2/4/e/h4ouhXc6m6xIpKnEqhIeAGYSpJiBcPzCZrRvniBSeVau2VANz54ZqYPDgMMOoqgGRbcRJCN3c/UIzuxjA3deZaZRNCpR7aDaTSsHkyaEaOPzw0H/gootCUhCRWsVJCF+ZWRuiaSXMrBuwPqdRidTXypWh4UwqFZrTt2sHl10WqoFevZKOTqQoxEkIPweeBfYyswrgW8CQXAYlEot7aESfSoVWlOvXw5FHwqhRcOGF4fRREYktToOcF8zsbeAowoR017j7JzmPTKQun3wSmtGnUjBvXjhF9HvfC9VAjx5JRydStOpMCGZ2WI1Fy6P7UjMrdfe3cxeWSA3u8OqrIQk8/jh89RUcfTQ89BBccIGqAZEsSFch/CrNaw6clOVYRLa1YgWMHRsSwfvvQ4cOYSqJoUPhkEOSjk6kSakzIbj7ifkMRORrmzfDyy+HM4OeeAI2bIBjj4X/+q/QmL5Nm6QjFGmS0vZUhtD/wMz+w8yeMLPHzezaqCeCSL1k7A/88cdw222s7bw/9O3LqknPM6r1D3jqtrlh8HjQoK2SQTb7IIsI6ae/jmZCfRQYBZwY3VLA5DhTqWb7pumvi1ed/YHHbXJ/4QX3Cy5w3247d/DXWhznAxjvO7Cu0f2G1ZdYJIvTX5vZbHfvkWlZPmj66+JVVgaLFm15vhsfMYQxfL/VSLpu/AB22QWGDOHkiUP507Lu27y/a1dYuLDu7TV2PZGmLO701xkPGQEzzeyoahs+EnijMcFJ87N4MRibOYXnmUx/lrAXt3IDH2wshUcegQ8/hF/9ipeXb5sMqt6f7nlj1xOReAnhSOAvZrbQzBYCbwLHm9kcM3snp9FJ07B8Obe2/yV/pxvP822O51Xu5loO4P+4rOvLYYK51mFYKtv9htWXWCS+OAnhNGBv4Pjotjeh9eVZwNm5C02K2qZN8OyzcN55sNdeDFs9nEUt9uFCJrInSxnGHSwtOaDBfYTVl1gkB+IMNAA7A4cCh1Xd4rwv2zcNKheBDz90v+mm0AAY3Dt1ch82zH3+/Kz3EVZfYpF4yOKg8k2EuYv+TjTBXcgjnvcL0zSoXKA2bYLnngsXj02dGp737RuazpxzDmy/fdIRijRr2eyp/F3CFNhfNT4saVKWLoXRo+HBB2HJEth1V/jxj8NVxN26JR2diNRTnITwLtAB+DjHsUgx2LgxjA2kUvDHP4arik85Be66C/r1UzUgUsTiJIRbCKeevku1Pgju3i9nUUnhWbIkTCs9alSoDHbbDa67Lswyus8+SUcnIlkQJyGMBW4D5gCbcxuOFJSNG+Hpp0M18Mwz4ULfb38b7rkHzj4bttsu6QhFJIviJIRP3P3enEcihWPRoi3VwLJl0Lkz3HBDqAbKypKOTkRyJE5CmGFmtwBT2PqQkfohNCUbNoQxgVQqjBEAnH46/Pa3cOaZqgZEmoE4CaGqIe1R1ZapH0JT8Y9/hEpg9GhYvhy6dIGf/hSuuCJM+CMizUacFprqi9DUbNgATz0VqoHnnwczOOOMcN3A6adDqzi/E0SkqYn1f76ZnQkcDHzdB8Hdf5GroCRHPvggXDMwejT885+w557w85/D5ZfDXnslHZ2IJCxjQjCz3wElhF4IDwL9gb/lOC7Jlq++gilTQjXwwguhS8xZZ4Vq4LTToGXLpCMUkQIRp0I4xt0PNbN33P2/zexXwBO5DkwaacGCUA089FDoRFZaCr/4BVx2WagMRERqiJMQ1kX3lWbWBVhJmPG0wczsDsJMqV8R5ki6zN1XN2abQqgGnnwyVAMvvRR+/Z99dqgGTj1V1YCIpBVn+uupZtYBuAN4G1gITGjkfl8AvunuhwLzgRsaub0mo0H9f+fPh2HDwi//Cy/k89kLuKPDzeyxaTFlM/+HilWnx0oG9d23ehWLNDFxpkStugE7AO3r854Y2zwXqIizblOf/rpe/X+//NJ9wgT3E08MK7Zs6X7eef7SsGd9xzYb691DuL69h9WrWKR40Njpr83scGCJu38UPb8UOB9YBNzo7quykZDM7ClgkruPz7RuU5/+Olb/33nzYORIGDMGVq6EvfcOs4sOGQKdOze4h3B936dexSLFI+701+kSwttAX3dfZWbHAROBq4CewIHu3j9DAC8Cu9fy0nB3/0O0znCgD3Ce1xGImZUD5QClpaW9F9X2LdREtGgRfmvX1JovWVfxRBgbePXVcJ3Ad74TxgZOPjm8McM2zMLEpPXdd13va+h+RCT/stEPoWW1KuBCIOXujwOPm9msTBt2974ZAhxMaMN5cl3JINpOCkhBqBAy7beYlZZu/au7O+8xlJFc1mIsDFwVZhW99dZQDey2W6xtVF9en31nel9D9yMihSvdoHJLM6tKGCcDf6r2WqMuZTWz04DrgH7uXtmYbTUlI0bALm3WMZDxvMpxvMdB/JD7+OzwvvDii/D++2HK6TqSQdU2GtJDuL7vU69ikSaorsEFYDjwBvAHYCZbDi/tC7wRZ4AizbYXAEuAWdHtd3He16QHld991/2aa/zLtju7g89nXx/R4Xaf/Nt/1ntTDe0hXN/3qVexSHEgGz2VzewooDPwvLt/ES3bH9jRE5jttMkNKq9bB5Mnh7GBN94IM4qef34YGzj++K3GBkREGiorPZXd/a+1LJvfmMAEmDMnnCn08MOwejXsvz/ceSdceil06pR0dCLSTGlay3yprIRHHw3VwJtvht7D/fuHauC448LpOSIiCVJCyLXZs0M1MH48rFkD3buHhvSDBkHHjklHJyLyNSWEXPjiC5g0KVQDb70FO+wAF1wQqoFjj1U1ICIFSQkhm2bNCklg/HhYuxYOPBDuvjtUA7vsknR0IiJpKSE01uefw8SJIRFMmwatW8N3vxuqgWOOUTUgIkVDCaGh3n47JIGKipAUDj4Y7r0XLrkEdt456ehEROpNCaE+1q6FCRNCIpgxA9q02VINHH20qgERKWpKCJm4hy//VAoeeSQMGB9yCNx3HwwcCB06JB2hiEhWKCHU5bPPQgJIpWDmzDBRz0UXhWrgiCNUDYhIk6OEUJ17GBhOpcKhocpK6NED7r8fBgyA9u2TjlBEJGc0WQ6EC8buvx969YIjjwxnDQ0YAH/7W6gOvv/9WpOBWkiKSFPSfCsE93DRWCoVLiKrrAwJ4Xe/g4svhp12Svv2iopw9Kgymrx70aLwHMLQgohIsUk722mhycpsp6tXhwvHUqkwydyOO4ZqoLwceveOvRm1kBSRYpGV2U6bDPcwoVwqFSaYW7cO+vQJzy+6CNq1q/cmFy+u33IRkULXtBPCp5+GKaZTKZg7N3zxDx4cmtIfdlijNq0WkiLS1DS9QWV3+POfQ2+BLl3gmmugbVt48EFYtgweeKDRyQDUQlJEmp6mUyGsWgXjxoVq4L33wqDw5ZeHaqBnz6zvrmrgePjwcJiotDQkAw0oi0ixKv6EMHs23HEHPPYYrF8PRx0Fo0eHKSXats3prgcOVAIQkaaj+BPCwoUwdWqoBIYOhUMPTToiEZGiVPwJ4cwzw9hAzQP6IiJSL8WfEFq1CjcREWmUpneWkYiINIgSgoiIAEoIIiISUUIQERFACUFERCJKCCIiAighiIhIRAlBRESAhBOCmf3YzNzMOiYZh4iIJJgQzGwv4BRALWVERApAkhXCr4FhQPH08BQRacISSQhm1g/40N1nJ7F/ERHZVs5mhTOzF4Hda3lpOPCfwKkxt1MOlAOUqj+liEjOmHt+j9iY2SHAS0BltGhPYBlwhLt/lO69ffr08enTp+c4QhGRpsXMZrh7n0zr5X3eaHefA+xa9dzMFgJ93P2TfMciIiJb6DoEEREBCqBBjruXJR2DiIioQhARkYgSgoiIAEoIIiISUUIQERFACUFERCJKCCIiAighiIhIRAlBREQAJQQREYkoIYiICKCEICIiESUEEREBijwhVFRAWRm0aBHuKyqSjkhEpHglPttpQ1VUQHk5VEZtdhYtCs8BBg5MLi4RkWJVtBXC8OFbkkGVysqwXERE6q9oE8LixfVbLiIi6RVtQigtrd9yERFJr2gTwogRUFKy9bKSkrBcRETqr2gTwsCBkEpB165gFu5TKQ0oi4g0VNGeZQThy18JQEQkO4q2QhARkexSQhAREUAJQUREIkoIIiICKCGIiEjE3D3pGGIzsxXAojzvtiPwSZ73mQ2KO78Ud34p7vrp6u6dMq1UVAkhCWY23d37JB1HfSnu/FLc+aW4c0OHjEREBFBCEBGRiBJCZqmkA2ggxZ1fiju/FHcOaAxBREQAVQgiIhJRQqjBzC4ws7lmttnM6jwbwMwWmtkcM5tlZtPzGWMd8cSN+zQzm2dmC8zs+nzGWEc8u5jZC2b2fnS/cx3rbYo+61lmNiXfcVaLI+3nZ2Y7mNmk6PW3zKws/1FuK0bcQ8xsRbXP+HtJxFkjptFm9rGZvVvH62Zm90Z/0ztmdli+Y6xNjLhPMLM11T7rn+U7xjq5u27VbsCBwAHAK0CfNOstBDomHW994gZaAn8H9gG2B2YDByUc9+3A9dHj64Hb6ljv8wL4jDN+fsC/Ab+LHl8ETCqSuIcA9yUda42YjgMOA96t4/UzgGcAA44C3ko65phxnwBMTTrO2m6qEGpw9/fcfV7ScdRXzLiPABa4+wfu/hUwETgn99GldQ4wNno8FvhOgrFkEufzq/73PAacbGaWxxhrU4j/3TNy99eAVWlWOQcY58FfgQ5m1jk/0dUtRtwFSwmh4Rx43sxmmFl50sHEtAewpNrzpdGyJO3m7ssBovtd61ivtZlNN7O/mllSSSPO5/f1Ou6+EVgDfCMv0dUt7n/386NDL4+Z2V75Ca1RCvHfc1xHm9lsM3vGzA5OOpgqRd0gp6HM7EVg91peGu7uf4i5mW+5+zIz2xV4wcz+L/plkDNZiLu2X6o5P80sXdz12Exp9HnvA/zJzOa4+9+zE2FscT6/RD7jDOLE9BQwwd3Xm9mVhCrnpJxH1jiF+FnH8TZhKonPzewM4Elgv4RjApppQnD3vlnYxrLo/mMz+x9CWZ7ThJCFuJcC1X/57Qksa+Q2M0oXt5n908w6u/vyqNz/uI5tVH3eH5jZK0AvwnHxfIrz+VWts9TMWgHtSf7wQca43X1ltacjgdvyEFdjJfLvubHc/bNqj582s/vNrKO7Jz43kw4ZNYCZtTWzdlWPgVOBWs8oKDDTgP3MbG8z254w6JnYGTuRKcDg6PFgYJtKx8zW579MAAAE+UlEQVR2NrMdoscdgW8B/5u3CLeI8/lV/3v6A3/yaCQxQRnjrnHsvR/wXh7ja6gpwKXR2UZHAWuqDj8WMjPbvWpcycyOIHwPr0z/rjxJelS70G7AuYRfHuuBfwLPRcu7AE9Hj/chnKkxG5hLOGRT8HFHz88A5hN+XRdC3N8AXgLej+53iZb3AR6MHh8DzIk+7znAFQnGu83nB/wC6Bc9bg1MBhYAfwP2Sfozjhn3LdG/5dnAy0D3Aoh5ArAc2BD9274CuBK4MnrdgN9Gf9Mc0pwVWGBx/7DaZ/1X4JikY6666UplEREBdMhIREQiSggiIgIoIYiISEQJQUREACUEERGJKCFI3pnZN6rN9PiRmX0YPV5tZnm9vsDMekZXi1Y979fQWWCjGXA71rK8vZmNM7O/R7eKumZ1bYx0f4uZ3WhmP872PqVpUUKQvHP3le7e0917Ar8Dfh097glszvb+oiuG69KTcI5+VWxT3P3WLIcwCvjA3bu5ezfCNQpjsrwPyM/fIk2YEoIUmpZmNtJCb4fnzawNgJl1M7Nno8kEXzez7tHyrmb2UjQp20tmVhotH2Nmd5nZy8Bt0dXlo81smpnNNLNzoqt2fwFcGFUoF1roC3BftI3dzOx/oknIZpvZMdHyJ6M45maa2NDM9gV6AzdVW/wLoIeZHRDNjT+12vr3mdmQ6PHPonjfNbNUtatbXzGz28zsb2Y238z+JdPfUiOmuj7LC6J9zTaznE7DIoVJCUEKzX7Ab939YGA1cH60PAVc5e69gR8D90fL7yNMgXwoUAHcW21b+wN93f1HhIn0/uTuhwMnAncA2wE/I/Qs6Onuk2rEci/wqrv3IMxvPzdafnkURx/gajNLN5vpQcAsd99UtSB6PJPQwyKd+9z9cHf/JtAGOKvaa63c/QjgWuDnHqa1Tve3VFfXZ/kz4NvR39svQ2zSBDXLye2koP3D3WdFj2cAZWa2I2H6ism2pbXADtH90cB50eOHCQ13qkyu9kV8KtCv2nH01kBphlhOAi6Fr7/E10TLrzazc6PHexGSWF1z0Ri1z8AZp0fCiWY2DCgBdiEkpKei156I7mcAZTG2FXaa/rN8AxhjZo9W2740I0oIUmjWV3u8ifDLuAWwOhpnyKT6l+8X1R4bcL7XaCJkZkfWJzgzOwHoCxzt7pUWZl5tneYtc4FeZtbC3TdH22gBHEqYBrmUrSv11tE6rQm/3Pu4+xIzu7HGfqo+p03U7//jOj9Ld78y+jzOBGaZWU/fehZUaeJ0yEgKnofpgv9hZhfA1710e0Qv/4UweyfAQODPdWzmOeCqasfhe0XL1wLt6njPS8D3o/VbmtlOhOmsP42SQXdC68Z0sS8gHB76abXFPwVecvfFwCLgIAu9mNsDJ0frVH35fxL9qu+fbj8x/paqeOr8LM2sm7u/5e4/Az5h66mlpRlQQpBiMRC4wsyqZpitagF5NXCZmb0DDAKuqeP9NxHGDN6x0Py8apD3ZcIX8iwzu7DGe64hHLaZQzg0czDwLNAq2t9NhNkqM7mcMP30AjNbQUgiVwK4+xLgUeAdwhjIzGj5akJfgjmEBirTYuwn3d9SXV2f5R1mNif6fF4jzMYpzYhmOxXJIzM7AHiaMKj7dNLxiFSnhCAiIoAOGYmISEQJQUREACUEERGJKCGIiAighCAiIhElBBERAZQQREQk8v88TRJ8boebQwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = sm.qqplot(lmod.resid, line=\"q\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Do the Levene test:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" coag \n",
" diet \n",
" meds \n",
" mads \n",
" \n",
" \n",
" \n",
" \n",
" 1 \n",
" 62 \n",
" A \n",
" 61.0 \n",
" 1.0 \n",
" \n",
" \n",
" 2 \n",
" 60 \n",
" A \n",
" 61.0 \n",
" 1.0 \n",
" \n",
" \n",
" 3 \n",
" 63 \n",
" A \n",
" 61.0 \n",
" 2.0 \n",
" \n",
" \n",
" 4 \n",
" 59 \n",
" A \n",
" 61.0 \n",
" 2.0 \n",
" \n",
" \n",
" 5 \n",
" 63 \n",
" B \n",
" 65.5 \n",
" 2.5 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" coag diet meds mads\n",
"1 62 A 61.0 1.0\n",
"2 60 A 61.0 1.0\n",
"3 63 A 61.0 2.0\n",
"4 59 A 61.0 2.0\n",
"5 63 B 65.5 2.5"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coagulation['meds'] = coagulation.groupby('diet').transform(np.median)\n",
"coagulation['mads'] = abs(coagulation.coag - coagulation.meds)\n",
"coagulation.head()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" df \n",
" sum_sq \n",
" mean_sq \n",
" F \n",
" PR(>F) \n",
" \n",
" \n",
" \n",
" \n",
" diet \n",
" 3.0 \n",
" 4.333333 \n",
" 1.444444 \n",
" 0.649189 \n",
" 0.592646 \n",
" \n",
" \n",
" Residual \n",
" 20.0 \n",
" 44.500000 \n",
" 2.225000 \n",
" NaN \n",
" NaN \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" df sum_sq mean_sq F PR(>F)\n",
"diet 3.0 4.333333 1.444444 0.649189 0.592646\n",
"Residual 20.0 44.500000 2.225000 NaN NaN"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmodb = smf.ols('mads ~ diet', coagulation).fit()\n",
"sm.stats.anova_lm(lmodb)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1 62\n",
"2 60\n",
"3 63\n",
"4 59\n",
"Name: coag, dtype: int64"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"coagulation.coag[coagulation.diet == \"A\"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Do the Bartlett test:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BartlettResult(statistic=1.6679561088149575, pvalue=0.6440812243179765)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"stats.bartlett(coagulation.coag[coagulation.diet == \"A\"],\n",
" coagulation.coag[coagulation.diet == \"B\"],\n",
" coagulation.coag[coagulation.diet == \"C\"],\n",
" coagulation.coag[coagulation.diet == \"D\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Pairwise comparisons"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assemble the pieces necessary for pairwise comparisons:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 1.183216\n",
"diet[T.B] 1.527525\n",
"diet[T.C] 1.527525\n",
"diet[T.D] 1.449138\n",
"dtype: float64"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod.bse"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 6.100000e+01\n",
"diet[T.B] 5.000000e+00\n",
"diet[T.C] 7.000000e+00\n",
"diet[T.D] -1.776357e-14\n",
"dtype: float64"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod.params"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.0859634472658364"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"stats.t.ppf(0.975,20)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1.8136382, 8.1863618])"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod.params[1] + np.array([-1, 1]) * stats.t.ppf(0.975,20) * lmod.bse[1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Can do the Tukey pairwise:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"Multiple Comparison of Means - Tukey HSD,FWER=0.05 \n",
"\n",
" group1 group2 meandiff lower upper reject \n",
" \n",
"\n",
" A B 5.0 0.7244 9.2756 True \n",
" \n",
"\n",
" A C 7.0 2.7244 11.2756 True \n",
" \n",
"\n",
" A D 0.0 -4.0562 4.0562 False \n",
" \n",
"\n",
" B C 2.0 -1.8242 5.8242 False \n",
" \n",
"\n",
" B D -5.0 -8.5773 -1.4227 True \n",
" \n",
"\n",
" C D -7.0 -10.5773 -3.4227 True \n",
" \n",
"
"
],
"text/plain": [
""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"thsd = sm.stats.multicomp.pairwise_tukeyhsd(coagulation.coag, coagulation.diet)\n",
"thsd.summary()"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAF1CAYAAAAeDGUGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHU9JREFUeJzt3XuUZVddJ/DvrwkkFCEEEyLPVIkIEZQJ0DxEB1p8IG/wAcRGgjaUOOgMLBEVXAqjLb6WgIMMdKMYoBQQeSMMMkzHBzbQjSCCyMsuEyIhCQkhdGgg7PnjnArVlarq3enq3Kruz2etXlX33HPO/d19z63zPXvve7taawEA4NA2TboAAICNQnACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITG1pVtaq60yr3f7SqtnTua19V/eCaFbdBVNU7qurcSddBUlUz4zF9wnh7V1U9+Sg8Tvf7omNfJ1bVx6rq1muwr1dX1XPXoKzF+7zNWN9N1nK/HL8EJyZiDClfrarTlyz/0HjimLke+/yzqvqtxctaa3drre06omL7HvuUqnphVf1HVV1VVZ8ab59+6K0nq7X2kNbaeZOuI7n2uLh6bMPLq+rtVXWHzm23VNWFR7vGtTDW2qrqWUewj+dW1dfGtrqiqt5bVd/Ts+0avy9mk/xta+1zYwi/avz3tfE9vnD7pWv0eIeltfafSf4+ybZJPD7HHsGJSfr3JOcs3Kiq705y08mVc/2MV7L/N8ndkvxIklOS3D/JZUnuM8HSVlWD9fg34BGttZOT3CbJxUn+14TrORrOTfKF8eeReO3YVrfKEA7eUFV1JDtc6O06DD+b5FXJtSH85LGmuSS/t3C7tfbUI6nrCM2NdcIRW49/NDl+vCrJExfdPjfJKxevsHSooqqeVFV/v3RHVTWbZGuSZ41Xt28dl187/DZeob++ql5bVV+qqg9W1X9ZrrCq2lRVv1JVn66qy6rqdVX1LSs8jycmOTPJY1prH2utfaO19vnW2m+21v563N93js/linGY5JGLHuvPquoli67W/6Gqbj32WF1eVR+vqnssWn9fVf3qOPxweVW9oqpOGu+7ZVW9raouGe97W1Xdfkl7bq+qf0iyP8kdF7dxVd2pqs6vqi9W1aVV9dpF296/qj4w3veBqrr/kv3+5lj7l6rqXQu9bVV10jgEc9n4/D9QVd+6Qlteq7X2lSSvT3LXRY9zYlX9wdizd3FVvbSqblpVN0vyjiS3XdTDcdux92qhjl+rqq9X1Snj7d+qqheutt9Fj/vwGnpDF3p27r7k9XhmVf3z2DavXXg9llNVU0l+PMnTknxHVW0+VFt0tNXXkpyX5NZJTquqb6+q94xtfmlVzVXVqUtqXvq+eHVVXZnkSVV1n6raU1VXju3xhys8lzOTfHuS9/XUWVVPrqpdi26fUCv0MNfQi/u3VfWCGpxUVX9YVReMNb1k0XH/8ap6yKJtTxyP/+8aF/1jkrOq6nY9dcJqBCcmaXeSU2oIFTdK8rgkr74+O2qt7cjBV7iPWGHVRyX5yyTfkuTPk7ypqm68zHr/PcmjkzwwyW2TXJ7kj1fY5w8meWdr7arl7hz3/9Yk70pyRpJfSDJXVXdZtNpjk/xaktOTHMjwh/6D4+3XJ1l64tqa5MEZTlp3HrdNhvf0K5JMZwhzVyd58ZJtfyrD8MrNk8wvue83xzpvmeT2GXt7agiNb0/yR0lOG+t5e1Wdtmjbn0zy0+NzvEmSZ47Lz01yiyR3GLd96ljXqsaA8bgMx8mC3x2f79lJ7pTkdkl+vbX25SQPSXLRoh6Oi5J8IMNrmCQPGJ/v9y66ff5q+x3ruGeSP83QY3FakpcleUtVnbiorsdm6G38tiR3T/KkVZ7ajyW5KsNx+H9y8MXD9TLW8qQkF7bWLk1SSZ6f4dj9zgxt/9xVdvGoDMfZqRneRy9K8qLW2ikZjrHXrbDddyf5TGvt60f6HBYbw+57kryntfaMNvzfYH+Qb7bvdySZSfKccZNXJnnCol08PMm+1tq/JElr7atJPpNk2QslOByCE5O20Ov0Q0k+nuSzR/nx9rbWXj9eof9hkpOS3G+Z9X42yXNaaxe21g5kOOn8eC0/jHFakv9c5THvl+TkJL/TWvtqa+09Sd6WRcOUSd7YWts79rK8MclXWmuvbK1dk+S1Se6xZJ8vbq1d0Fr7QpLtC/tqrV3WWvur1tr+1tqXxvseuGTbP2utfbS19vWxHRb7WobQddvW2ldaawu9ew9L8snW2qvG7f4iw+u1OKC+orX2idba1RlOtGcv2udpSe7UWrtmfJ5XrtJeb6qqK5JcmeG4+P1kGFpM8pQkz2itfWF8fr+d5PGr7Ov8JA8cX7e7Zwh+Dxx7Ku6d5O869vuUJC9rrb1vrP+8DOF28XHzR621i8bX462Lnvtyzs0wxHZNhvB+zgrhvcdjx7a6IMm9MoT9tNY+1Vr7m9bagdbaJRmO9aXHwWL/2Fp709hbenWG1+xOVXV6a+2q1truFbY7NcmXrmftK7ldhtdtrrX23GToAU7y5CRPb61dPh4/z883X6NXJXlEVZ083v6pcdliXxrrhSMiODFpr8rQU/GkLBmmO0ouWPiltfaNJBdmuCpfajrJG8ehmSuS/GuSa5IsN8R0WYb5OCu5bZILxsdbMJ/hBLHg4kW/X73M7ZNzsAsW/T6/8ByqaqqqXlZV8+Owy98mOXXs0Vtu26WelaG34v01DCn+zKLnsLR3aulz+Nyi3/cvqvlVGXpWXlNVF1XV7x0iKDy6tXZqkhOT/HyS82v4xNatkkwl2bvodXnnuHwl5yfZkuSeST6S5G8yBIj7JfnU2DtzqP1OJ/nFhfvG+++Qg4+blZ77QWqY6P79GXp1kuTNGcL7w1Z5Dqt5XWvt1NbaGa21B7XW9o6Pc0ZVvaaqPjseB6/O0Hu5kqXHxLYMPXAfr2Fo9eErbHd5hp7LtfTIJDdOsnPRsltnOB4+vOg1eFuG3s201i5I8v4kjxl7R384Qyhd7OZJrljjWjkOCU5MVGttPsMk8YcmecMyq3w5w0ltwWofeW4dD3ntJ7TGq9jbJ7lomfUuSPKQ8aS08O+k1tpyPWLvTvLgGubZLOeiJHeogydin5kj611b/EmzM/PN5/CLSe6S5L7jMMsDxuWLJwyv2E6ttc+11p7SWrtthl63l9TwdQ8XZQgQi3U9h9ba11prz2ut3TXDpPmHp2N4auzdeUOGwPp9SS7NECLvtug1ucU4EXml5/XeDO3xmCTnt9Y+Ntb9sHxzmO5Q+70gyfYlx8LU2Ot2uH4qw9/dt1bV5zIMH52UNRiuW+L5Gdrj7uNx8IQcfAwsdVDbtdY+2Vo7J0Mw+d0kr1/h+P7nDPPkeieU97yfX5rk/2UYCl5Y9+IkX01ylyWv0S0WbXdehuf5uIyf8lu4o4YPcNwxyYc764QVCU6sB9uSPGicp7LUh5L86NiTcqes/pHiizP8cVzNvarqR8c/9E/PMOSy3DDES5Nsr6rpJKmqW1XVo1bY56synFz/qqrOqmFi+WlV9eyqemiGibNfzjBx/cY1fH/OI5K85hC1ruZpVXX78er62RmG85LhqvrqJFeM9/3G4ey0qn6ivjmZ/PIMJ9Rrkvx1kjtX1U+OE3ofl2HS9ts69vn9VfXdY6/XlRmGga7p2K7GNr9lkn8de+x2JnlBVZ0xrnO7qnrwuMnFGSZGX3syba3tT7I3w0TshaD03gyh8PxxnUPtd2eSp1bVfceablZVD6uq69PT8sQkz8swlLfw78eSPGzJfLEjdfMM86iuGCdE/9LhbFxVT6iqW41ts9BLc53XrLV2YZJPpv/Tox9OcvfxeLhplj8+W4Z5cJ/JMJfspHFY8+VJXji+F2s8/n940XZvSHLfDL2US3uv75fkEytc+MBhEZyYuNbap1tre1a4+wUZrjQvznBFObfCeknyJ0nuOnblv2mFdd6c4Yr08gxX/z+6zDyfZJgc+5Yk76qqL2UIV/ddof4DGSaIfzzDUNCVGYYNTk/yvnFi6iMzTF6+NMlLkjyxtfbxVZ7Lofx5hkncnxn/LXx/1QszfKXDpWPN7zzM/d47yfuq6qoMz/9/tNb+vbV2WYaeol/MMDT5rCQPH4e6DuXWGSYeX5lhyPP8rP4hgLeOj39lhjla57bWPjre98tJPpVk9zgE9e4MPUoZ2/MvknxmPAYWhtLOzzD08/5Ft2+eYRgzHfvdk2Ge04szHDefyuqTv5dVVffLMKH5j8eevYV/bxn3ec6qOzg8z8swPPnFDJP6l+vNXc2PJPno+Dq8KMnjx/l3y3lZhvfSIY09fr+dZFeSf8vBr8Hi9VqGi6TPZxgyPzHDsTef4XX8Yobj/zsWbfPlJG/K0KO49P2/NcPFEByxGo5POPbV8I3Ed2qtPeFQ665nVbUvyZNba++edC0whpp/SvIDbfiyyUnW8j+TnNlae9KiZbfJ8D1rZ48XMXBEDveLzgDgWmOP610PueJRNg51/nSGHuVrjWFu4vVx7DBUB8CGVlU/l+Q/kry5tfbeSdfDsc1QHQBAJz1OAACdBCcAgE5rOjm8hv9odTZJbnazm93rrLPOWsvdAwAcFXv37r20tbba/0SQ5CjOcdq8eXPbs2elr+YBAFg/qmpva23zodYzVAcA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQ64VArVNU1ST6S5MZJvp7kvCQvbK194yjXBgCwrvT0OF3dWju7tXa3JD+U5KFJfuPolsXxam5uLjMzM9m0aVNmZmYyNzc36ZIA4FqHNVTXWvt8ktkkP19VdXRK4ng1NzeX2dnZzM/Pp7WW+fn5zM7OCk8ArBuHHKpbqrX2maralOSMJBevfUnHny1btky6hHVh9+7dOXDgwEHL9u/fn23btmXnzp0Tqmr92LVr16RLADjuXd/J4cv2NlXVbFXtqao9l1xyyRGUxfFoaWg61HIAuKFVa231Faquaq2dvOj2HZN8IMnpbZWNN2/e3Pbs2bNmhXLsm5mZyfz8/HWWT09PZ9++fTd8QQAcN6pqb2tt86HWO6wep6q6VZKXJnnxaqEJro/t27dnamrqoGVTU1PZvn37hCoCgIP1BKebVtWHquqjSd6d5F1Jnnd0y+J4tHXr1uzYsSPT09OpqkxPT2fHjh3ZunXrpEsDgCQdQ3XXl6E6AGCjOCpDdQAAxzPBCQCgk+AEANBJcAIA6CQ4AQB0EpwAADoJTgAAnQQnAIBOghMAQCfBCQCgk+AEANBJcAIA6CQ4AQB0EpwAADoJTgAAnQQnAIBOghMAQCfBCQCgk+AEANBJcAIA6CQ4AQB0EpwAADoJTgAAnQQnAIBOghMAQCfBCQCgk+AEANBJcAIA6CQ4AQB0EpwAADoJTgAAnQQnAIBOghMAQCfBCQCgk+AEANBJcAIA6CQ4AQB0EpwAADoJTgAAnQQnAIBOghMAQCfBCQCgk+AEANBJcAIA6CQ4AQB0EpwAADoJTgAAnQQnAIBOghMAQCfBCQCgk+AEANBJcAIA6CQ4AQB0EpwAADoJTgAAnQQnAIBOghMAQKfu4FRVt66q11TVp6vqY1X111V156NZHADckObm5jIzM5NNmzZlZmYmc3Nzky6JdeaEnpWqqpK8Mcl5rbXHj8vOTvKtST5x9MoDgBvG3NxcZmdns3///iTJ/Px8ZmdnkyRbt26dZGmsI13BKcn3J/laa+2lCwtaax86OiUBsBa2bNky6RI2lN27d+fAgQMHLdu/f3+2bduWnTt3TqiqjWXXrl2TLuGo6x2q+64kew+1UlXNVtWeqtpzySWXHFllAHADWhqaDrWc41Nvj1OX1tqOJDuSZPPmzW0t9w3A4Tkerv7X0szMTObn56+zfHp6Wltyrd4ep48mudfRLAQAJmn79u2Zmpo6aNnU1FS2b98+oYpYj3qD03uSnFhVT1lYUFX3rqoHHp2yAOCGtXXr1uzYsSPT09OpqkxPT2fHjh0mhnOQaq1vRK2qbpvkhRl6nr6SZF+Sp7fWPrnc+ps3b2579uxZozIBAI6eqtrbWtt8qPW65zi11i5K8tgjqgoAYAPzzeEAAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6dQWnqrqmqj5UVR+uqg9W1f2PdmEA9Jmbm8vMzEw2bdqUmZmZzM3NTbokOGad0Lne1a21s5Okqh6c5PlJHnjUqgKgy9zcXGZnZ7N///4kyfz8fGZnZ5MkW7dunWRpcEzqDU6LnZLk8rUuBGDBli1bJl3ChrF79+4cOHDgoGX79+/Ptm3bsnPnzglVtfHs2rVr0iWwQfQGp5tW1YeSnJTkNkketNxKVTWbZDZJzjzzzDUpEICVLQ1Nh1oOHJlqrR16paqrWmsnj79/T5KXJ/mutsrGmzdvbnv27FmzQgG4rpmZmczPz19n+fT0dPbt23fDFwQbVFXtba1tPtR6h/2putbaPyY5Pcmtrk9hAKyd7du3Z2pq6qBlU1NT2b59+4QqgmPbYQenqjoryY2SXLb25QBwOLZu3ZodO3Zkeno6VZXp6ens2LHDxHA4SnqH6q5J8pGFm0me3Vp7+2rbGKoDADaK3qG6rsnhrbUbHXlJAAAbm28OBwDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACdBCcAgE6CEwBAJ8EJAKCT4AQA0Kk7OFXVY6qqVdVZR7MgAID16nB6nM5J8vdJHn+UaoHMzc1lZmYmmzZtyszMTObm5iZdEgBcqys4VdXJSb43ybYIThwlc3NzmZ2dzfz8fFprmZ+fz+zsrPAEwLpxQud6j07yztbaJ6rqC1V1z9baB49mYZOyZcuWSZdw3Nq9e3cOHDhw0LL9+/dn27Zt2blz54SqOr7t2rVr0iUArCu9Q3XnJHnN+PtrxtvXUVWzVbWnqvZccskla1Efx5GloelQywHghlattdVXqDotyYVJPp+kJbnR+HO6rbLx5s2b2549e9awVI51MzMzmZ+fv87y6enp7Nu374YvCIDjRlXtba1tPtR6PT1OP57kla216dbaTGvtDkn+Pcn3HWmRsNj27dszNTV10LKpqals3759QhUBwMF6gtM5Sd64ZNlfJfnJtS+H49nWrVuzY8eOTE9Pp6oyPT2dHTt2ZOvWrZMuDQCSdAzVXV+G6gCAjWIth+oAAIjgBADQTXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdBKcAAA6CU4AAJ0EJwCAToITAEAnwQkAoJPgBADQSXACAOgkOAEAdDphLXdWVbNJZsebV1XVv63l/pOcnuTSNd7n8Ug7rg3tuHa05drQjmtHW66NjdSO0z0rVWvtaBeyZqpqT2tt86Tr2Oi049rQjmtHW64N7bh2tOXaOBbb0VAdAEAnwQkAoNNGC047Jl3AMUI7rg3tuHa05drQjmtHW66NY64dN9QcJwCASdpoPU4AABOzboNTVe2rqo9U1Yeqas+47Oyq2r2wrKruM+k617uqOrWqXl9VH6+qf62q76mqb6mqv6mqT44/bznpOjeCFdry98fb/1xVb6yqUydd53q3XDsuuu+ZVdWq6vRJ1rhRrNSWVfULVfVvVfXRqvq9Sde53q3w3na+OUxVdZexvRb+XVlVTz/WzjnrdqiuqvYl2dxau3TRsncleUFr7R1V9dAkz2qtbZlQiRtCVZ2X5O9aay+vqpskmUry7CRfaK39TlX9SpJbttZ+eaKFbgArtOV9kryntfb1qvrdJNGWq1uuHVtrV1TVHZK8PMlZSe61+L3P8lY4Ju+R5DlJHtZaO1BVZ7TWPj/RQte5FdrxdXG+ud6q6kZJPpvkvkmelmPonLNue5xW0JKcMv5+iyQXTbCWda+qTknygCR/kiStta+21q5I8qgk542rnZfk0ZOpcONYqS1ba+9qrX19XG13kttPqsaNYJVjMklekORZGd7nHMIqbflzSX6ntXZgXC40rWKVdnS+OTI/kOTTrbX5HGPnnPUcnFqSd1XV3vEbyZPk6Ul+v6ouSPIHSX51YtVtDHdMckmSV1TVP1XVy6vqZkm+tbX2n0ky/jxjkkVuECu15WI/k+QdN3xpG8qy7VhVj0zy2dbahydc30ay0jF55yT/tareV1XnV9W9J1vmurdSOzrfHJnHJ/mL8fdj6pyznoPT97bW7pnkIUmeVlUPyHAl9YzW2h2SPCPjFQIrOiHJPZP879baPZJ8OcmvTLakDWvVtqyq5yT5epK5yZS3YSzXjs/NMLT06xOsayNa6Zg8Icktk9wvyS8leV1V1cSqXP9Wakfnm+tpHO58ZJK/nHQtR8O6DU6ttYvGn59P8sYMc0nOTfKGcZW/HJexsguTXNhae994+/UZ/kBcXFW3SZLxp678Q1upLVNV5yZ5eJKtbb1OGlw/VmrHb0vy4XFu4+2TfLCqbj2ZEjeMldrywiRvaIP3J/lGhv8vjOWt1I7ON9ffQ5J8sLV28Xj7mDrnrMvgNHbd33zh9yQ/nORfMowxP3Bc7UFJPjmZCjeG1trnklxQVXcZF/1Ako8leUuGPwoZf755AuVtKCu1ZVX9SJJfTvLI1tr+iRW4QazQjh9srZ3RWptprc1kOJHdc1yXFazy/n5Thr+Pqao7J7lJNs5/snqDW6UdnW+uv3PyzWG65Bg756zLT9VV1R0z9DIlQzfqn7fWtlfV9yV50bjsK0n+W2tt74TK3BCq6uwMn1S6SZLPJPnpDIH5dUnOTPIfSX6itfaFiRW5QazQlh9IcmKSy8bVdrfWnjqZCjeG5dqxtXb5ovv3ZcknalneCsfkl5P8aZKzk3w1yTNba++ZWJEbwArteLc43xy2qppKckGSO7bWvjguOy3H0DlnXQYnAID1aF0O1QEArEeCEwBAJ8EJAKCT4AQA0ElwAgDoJDgBAHQSnAAAOglOAACd/j9jh5E6M8vTnAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = thsd.plot_simultaneous()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## False discovery rate\n",
"\n",
"Read in the data:"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" school \n",
" class \n",
" gender \n",
" social \n",
" raven \n",
" id \n",
" english \n",
" math \n",
" year \n",
" \n",
" \n",
" \n",
" \n",
" 1 \n",
" 1 \n",
" 1 \n",
" girl \n",
" 9 \n",
" 23 \n",
" 1 \n",
" 72 \n",
" 23 \n",
" 0 \n",
" \n",
" \n",
" 2 \n",
" 1 \n",
" 1 \n",
" girl \n",
" 9 \n",
" 23 \n",
" 1 \n",
" 80 \n",
" 24 \n",
" 1 \n",
" \n",
" \n",
" 3 \n",
" 1 \n",
" 1 \n",
" girl \n",
" 9 \n",
" 23 \n",
" 1 \n",
" 39 \n",
" 23 \n",
" 2 \n",
" \n",
" \n",
" 4 \n",
" 1 \n",
" 1 \n",
" boy \n",
" 2 \n",
" 15 \n",
" 2 \n",
" 7 \n",
" 14 \n",
" 0 \n",
" \n",
" \n",
" 5 \n",
" 1 \n",
" 1 \n",
" boy \n",
" 2 \n",
" 15 \n",
" 2 \n",
" 17 \n",
" 11 \n",
" 1 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" school class gender social raven id english math year\n",
"1 1 1 girl 9 23 1 72 23 0\n",
"2 1 1 girl 9 23 1 80 24 1\n",
"3 1 1 girl 9 23 1 39 23 2\n",
"4 1 1 boy 2 15 2 7 14 0\n",
"5 1 1 boy 2 15 2 17 11 1"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"jsp = pd.read_csv(\"data/jsp.csv\", index_col=0)\n",
"jsp.head()"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"jsp['mathcent'] = jsp.math - np.mean(jsp.math)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Try two different plotting methods"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEVCAYAAAALsCk2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xu8HWV97/HPLxcSTCKI4BZEQ1S0uRjhbEBtEbINkMjFSD0W06JAadKAhNNX1LOB9NWibQ4JWHxBkMTmRC6lCfRoUQmQC8naCqVUQhEISUGQoMi1ct2BbEnyO3/MrJ3Z67r3WmvWmpn1fb9e67X3mt/MPM88M2ue9TzPzCxzd0RERKKGtToDIiKSPKocRESkiCoHEREpospBRESKqHIQEZEiqhxERKSIKgeRQTKz7WZ2QqvzIdIMqhxESjCz683s7xOQjx4z+4tW50PajyoHEREpospBMiXs+vmGmT1sZjvMbKWZdZjZnWb2hpndZWbvCuf9f2b2vJm9ZmY/M7PJ4fS5wJ8B/9vMes3stkgSR4Trfs3MbjGz0ZG0Z5nZL8zsdTN70sxmhtP3C/PxnJn91sz+3syGh7GzzeweM/u2mb1iZk+Z2WfD2CLg08A1YT6uaUohiqDKQbLpC8CJwEeA04A7gUuAAwmO+QvD+e4EDgfeA/wn8M8A7v6P4f+Xu/tYdz8tsu4/AWYCE4CpwNkAZnYMcCPwDWB/4Dhge7jMDcAu4MPAkcBJQLSr6BPAY2H+LgdWmpm5+0LgbuCCMB8X1FcsIoM3otUZEInBUnd/AcDM7gZedPcHw/e3AtMB3P37+QXM7FLgFTPbz91fq7Duq9392XCZ24AjwunnAt939w3h+9+G83QAnwX2d/e3gB1m9h1gLvC9cN6n3X1FOP8NwLVAB/B87UUgUh9VDpJFL0T+f6vE+7Fht84i4IvAQcCeMH4gUKlyiJ6w3wQOCf9/P3BHifnHAyOB58wsP20Y8JtS63T3N8P5xlbIg0jsVDlIu/pTYBZwAkH3z37AK0D+DD7UxxX/BvhQmel9wIHuvquGfOqxydISGnOQdjWO4KT9O+AdwP8piL8AfHAI61sJnGNm081smJm9z8z+wN2fA9YD/2Bm7wxjHzKz4we53qHmQ6QhVDlIu7oReJpgbGArcF9BfCUwycxeNbMfVVuZu/8cOAf4DkG31E8JupQAvgLsE6bzCvAD4OBB5vMq4H+GVzJdPchlROpm+rEfEREppJaDiIgUUeUgIiJFVDmIiEiRhlQOZvZ9M3vRzLZEpl0aPirgF+Hr5EakJSIi8WtUy+F6gkcKFPqOux8RvkrdICQiIgnUkJvg3P1nZnZYves58MAD/bDD9q5mx44djBkzpuS8ijU2lrT8KJaMWNLyo1j9sQceeOC/3f2gkjNHxH2H9AVm9hVgM/A1d3+lcIbwCZhzATo6Ovj2t7/dH+vt7WXs2NJPEVCssbGk5UexZMSSlh/F6o91dXU9XXLGQu7ekBdwGLAl8r4DGE7QdbWI4KFkFdfR2dnpUblczstRrLGxpOVHsWTEkpYfxeqPAZt9EOf02K5WcvcX3H23u+8BVgDHxJWWiIg0VmyVg5lFHw9wOrCl3LwiIpIsDRlzMLPVwDTgQDN7BvhbYJqZHUHwVMntwF82Ii0REYlfo65Wml1i8spGrFtERJpPd0iLiEgRVQ4iIlJEvwQXivyEY/5S3EHFJP1q3fdpOi66urqAZOQzTeXWaGk61tRyCLk747vXlCz4SjFJv1r3fZqOiyTlM03l1mhpOtZUOYiISBFVDiIiUkSVg4iIFEnsgHStAzDNHrhJWj418Fhako6LpA08Vjpm8vkZaiwt4ijvJB1rUPs5IbEth1oHYJo9cJO0fCZpoC9JA49JOi6SNvBYKb18foYaS4s4yjtJxxrUfk5IbOUgIiKto8pBRESKqHIQEZEiiR2QjkPSBopqVesAU63LDWZQcqjrTdOAbVroYoTGysI21KOtWg5JGyiqVa3rrHW5wQxKNnJAPmkDtmmRpHLJwn7KwjbUo60qBxERGRxVDiIiUqQhlYOZfd/MXjSzLZFpB5jZBjP7Zfj3XY1IS0RE4teoAenrgWuAGyPTLgI2uvtiM7sofN/diMTafaAoDlm421XSLy0XKrTDOaghLQd3/xnwcsHkWcAN4f83AJ9vRFphem09UBSHLNztKumXlgsV2uEcFOelrB3u/hyAuz9nZu8pNZOZzQXmAnR0dNDT0zMgXvi+XWK9vb1Nz0tc6816uWV9++Jab9bLLS2xsty9IS/gMGBL5P2rBfFXqq2js7PTo8Z3r/Fysh7L5XJNTS+u9Wa93NKwfVMvXefju9f4+O41PvXSdUNKr1o8LbGsHxdDiQGbfRDn9DivVnrBzA4GCP++GGNaIlLGa2+9zfUzx7B98Sm89tbbrc6OpESclcNPgLPC/88CfhxjWiIi0kANGXMws9XANOBAM3sG+FtgMfAvZnYu8Gvgi41IS6ScrF9BkvXti0vWy+3j31wPwGEX3c5++47kob89qSHrbdTVSrPd/WB3H+nuh7r7Snf/nbtPd/fDw7+FVzOJNJRn/AqSrG9fXLJebnF1G7bVg/ca7ePfXN+/Mxpda4uItJIen1GH1956m+2LT9Fgn4hkjioHEREpom6lFqk0iJT1AbS0ULdh40XL9OPfXK/yTDC1HFqk0iBS1gfQ0kLdho0XLVOVZ7Kp5SAi0gRpa4mq5SAi0gRpa4mqchARkSLqViJ9zb1m0gCilFLpuNAxkw1qOZC+5l4z1TqA+PFvruewi27v/1+ypdJxoUHnbFDlILHQCUIk3TLVraTuIRGJW7ucZzLVclD3kIjErV3OM5lqOUj6peVbWaV8Jm0bxk28iPlPAzfAuIkAp7QsL1JatWMmrsdyV5KploOkX1q+lVXKZ9K24Y1ti1k6fimPnPUIb2xb3NK8SGnVjplW/JqfKgcRESkSe7eSmW0H3gB2A7vc/ai40xQRkfo0a8yhy93/u0lpiYhInVI3IJ20wb5yWpHPVgxapV1ajqdqat33lZZr9p3Ozd4XabqooJK4LjhoRuXgwHozc+B77v6P0aCZzQXmAnR0dNDT0zNg4cL3+YGZ3t5exo4dy9lrdwyYp6enh97e3v5phctXmtbI5arlE6gpvUqx1956m2uO9bLpVVtnpXz09vYCVnb5OJZrZNmUi9VzPNUaa/Q25Lej0r4vl5dKy0XL5oJ73i752Sy3f2uJDeYzU277a4lVSq8V55lq6ywXe2Pb4ro/9yW5e6wv4JDw73uAh4Djys3b2dnpUeO713ih/LRcLlc0T6VYreuMKy/ju9cMOb1qsXrWWW2ZXC5XMd04liucZ+ql63x895r+19RL1w1qnbWmF0es1nxWi1Xa9/XE3Evvw7hjhflJy76PI738+0adS4DNPohzd+xXK7n7s+HfF4FbgWPiTlOyqdZLRKPPeTrsotv1rKcUStrlwe0g1srBzMaY2bj8/8BJwJY40xQppBOLyNDFPebQAdwa/ibyCGCVu6+NOc0hGzfxIj52w0XBm4JBnUoDU5WWS5okDVanZbAvTfs3SaLlVlhmlWJZkJXzBcTccnD3X7n7x8PXZHdfFGd6tXpj22IeOeuRkneRVvrWWWm5pGnFHZaV8pKGb/JJ27/RCj7JXWPRcisss0qxLMjK+QJ0h7RIaiSpgpfsU+WQQGn5higi2aXKIYH0DVFEWi11d0hXkqYBn1rvakzqwLJ+K7i9xDGwXOvnt9kXOMRxnonr3FXP+SJTLYc0DfjU+hjlJLUq9FOg7SuOgeVaP7/NvsAhjvNMXOeues4XmWo5SG2S1AJIU+tPJMsy1XKoRHfJlpekFkCaWn8iWdY2lUNarq0XEUmCRHYrpeUOWmlflbq/ar3jvhWy/vvSSSvvWrViPyWy5aBv+ZJ0lbq/ar3jvhWy/vvSSSvvWrViPyWychARkdZS5SAiIkUSOeYgMlRZ6VsWSYrUVQ5puQ4+LfmMS7MfzZzvW+7p6WHatGn9ly1nSaVByTgGLLP+eG2pLHXdSmm5Dj4t+YxL1h/N3AqVBiXjGLDUPmxvqWs5NFu7twBEpD3F3nIws5lm9piZPWFmF8WdXqO1ewtARNpT3L8hPRz4LvBZYBIw28wmxZmmiIjUL+5upWOAJ9z9VwBmdjMwC9gac7ptKS2DkuqqKy9Jj2TPulrvcm8XcXcrvQ/4TeT9M+E0iUFaBiXVVVdekh7JnnW13uXeLuJuOViJaT5gBrO5wFyAjo4Oenp6AOjp6aG3t3fA+7wsxIBExfJ5BWtqbCjlVurbXE/PmLrWmbQY6LhIyr5IUgzi2b8VuXtsL+BTwLrI+4uBi8vN39nZ6e7u47vXuLt7Lpcb8D4rsfz7JMXyeW12rDA/7RzLv09SLJ9XHRfZOS6AzT6I83fc3Ur3A4eb2QQz2wf4EvCTmNMUEZE6xdqt5O67zOwCYB0wHPi+uz8aZ5pZV+ugc5oezdx/d/PaYFA2Thp4TAcNHjdf7DfBufsdwB1xp5PXzBNLK7yxbTHXzxwz5EdE1Lpcs21fHHyoD7vo9v7/4/TGtsWZf+xGFlTaT9qHcPbaHQ0/52XqDulmn1jikqZv+SIyeHF8ed2++JRYznmZqhyyIi3f8kVk8NL25TV1D94TEZH4tU3LodqgVdbHKirJ+qOZ6xmwzMpxEUeftGRb21QOlQat0tbca7Ro2Zy9dkers9NwtQ5YZuW4iKtPWrKtbSoHEUm+rLTUskCVg4gkQlZaalmhAWkRESmilkMLaZCwcdJyl2xa8pmX7+YpdYzWGsuypF38UM89U6ocWkSDhI2Vlrtk05JPGNjNU/i7ErXGsi5pFz/Uc89UIiuHNF12mqS8SHNp36ef9mF5iawc0nLZaZLyIs2lfZ9+2oeVJbJyEClH3/QkCdrhOFTlMAjtcCDUotl3Vsf1TS9N+7fSRQy6wKE5WtXiaPb+VeVQhZqe5WXhzuo07d9KFzHoAodsa8X+VeUQo7i+kablG2K7Xs4oyZG2S4eTJLbKwcwuBeYAL4WTLgl/+KctxPWNNC3fENv5ckZJjjRdOpw0cbccvuPu3445DZGq0jSuIJIE6lbKmLR0OTVTK8YVVBmVF0d3o8q78eKuHC4ws68Am4GvufsrMaeXGbWc5NPS5ZR1aRrkbrY4uhtV3vGoq3Iws7uA95YILQSWAX8HePj3H4A/L7GOucBcgI6ODnp6egDo6emht7d3wPtCpaZVisWxzjhi188cw9lrd3D9zDFl56klvUrbni8bsIbFKuWnnfdvtXw2exuqpRdHmrXEat2/aVkurthg9m9J7h77CzgM2FJtvs7OTnd3H9+9xt3dc7ncgPdRpaZVisWxzjTFxnevKbvt0bJpVKxSfsZ3rxnwmnrpukFvR6V15vNSbZ7BrjOOWLV8Trl+yoBX3PmsdFzElWYtsVr3b1qWiytWav8Cm30Q5+04r1Y62N2fC9+eDmyJKy1JF3UDlFfPg9JEGinOMYfLzewIgm6l7cBfxpiWiNQoLRcxNHvQud0HuWOrHNz9y3GtW0QaIy0XMTS7tanWrX4JTkRESlDlICIiRVQ5iIhIEd0h3UZqHXjUA/RE2k9bVQ7tfPVBrQOPeoCeSHtqm8pBVx+IiAyexhxERKSIKgcRESnSNt1KIknRzmNfkh6qHESaSGNfkhbqVhIRkSKqHEREpIgqBxERKZLYMQcN2jVXku6CNrPg7xLyPxY1KDpmRIrV+mSERFYOGrRrrqTdBe3u9PT0MG3atEEvo2NGpFg9j2RXt5KIiBRR5SAiIkXqqhzM7Itm9qiZ7TGzowpiF5vZE2b2mJnNqC+bIiLSTPWOOWwB/hj4XnSimU0CvgRMBg4B7jKzj7j77jrTExGRJqir5eDu29z9sRKhWcDN7t7n7k8BTwDH1JOWiIg0T1xXK70PuC/y/plwWhEzmwvMBejo6KCnp2dAvPC9YvHFalm2p6eH3t5ewMouW0t+ent7E1U2zYxFt71cmTcrL61Ks5Gx/DHarDJtdnr1xspy94ov4C6C7qPC16zIPD3AUZH33wXOjLxfCXyhWlqdnZ0eNb57jZejWGNjtSybn5bL5couW2t+crlcTculPTa+e03/tlcq82bkpVVpNjIWPUarzZPG9GqJAZu9yrnY3au3HNz9hKFXOTwDvD/y/lDg2RrWIyIiLRBXt9JPgFVmdiXBgPThwM9jSit2td6x2w6SdGe1iDROvZeynm5mzwCfAm43s3UA7v4o8C/AVmAt8FVP8ZVK7k4ul1PFUGD74lP677xMwp3VItI49V6tdKu7H+ruo9y9w91nRGKL3P1D7v5Rd7+z/qxKO5g/fz6jR4+mq6uL0aNHM3/+/FZnSaQtJfLZStKe5s+fz/Lly1myZAmTJk1i69atdHd3A7B06dIW506kvejxGZIYK1asYMmSJSxYsIDRo0ezYMEClixZwooVK1qdNZG2k7mWQxYGj6ttw9NLTk319pXT19fHvHnzBkybN28eX/va11qUI2m2OD6/epR7bTLXcsjC4HG1bUj79pUzatQoli9fPmDa8uXLGTVqVItyJM3W6M9v9KKJ7YtP0YUTQ5C5loOk15w5c/rHGCZNmsSVV15Jd3d3UWtCROKnykESIz/ofMkll9DX18eoUaOYN2+eBqNFWiBz3UqSbkuXLmXnzp3kcjl27typikGkRVLZcsjCoLPUph32fbnf/G2HbZfkSGXLIQuDzlKbrO/7SoOnWd92SZZUVg6SXbpDWiQZUtmtJNmkO6RFkkMtB0kM3SEtkhyqHGJkZsHdzOFAolRW7g7pvr6+FuVIpLEqnROSdr5Q5RAjDSAOje6QlqyrdE5I2vlCYw6SGLpDWiQ5VDlIYugOaZHkqPeX4L5oZo+a2R4zOyoy/TAze8vMfhG+lldaj0ie7pAWSYZ6Ww5bgD8Gvlci9qS7H1Hn+ptGd59mWxz7tx2Omf7B0cXF21cplhZZeER4XMdhvT8Tus3dH2tUZlopaYNB0lhx7N92OGby2zjUWFpk4RHhcR2HcY45TDCzB4HXgb9297tLzWRmc4G5AB0dHfT09AyIF77P6+3trSlWap0bN27kpptu4te//jUf+MAHOPPMM5k+fXpd64xzuVrXWS1Wz7KNjsWx/XGUW1z7Ignpxfm5qBTLcplWS68Vx3ZZ7l7xBdxF0H1U+JoVmacHOCryfhTw7vD/TuA3wDurpdXZ2elR47vXeDm5XK6mWOE6V61a5RMmTPBNmzb5hg0bfNOmTT5hwgRftWpVzeuMe7la11kpVs+yccTi2P44yi2OdSYhvbg/F0nYxiSm14xjG9jsVc7F7l69W8ndT3D3KSVeP66wTJ+7/y78/wHgSeAjQ6+64rdo0SJWrlxJV1cXI0aMoKuri5UrV7Jo0aJWZ02kZfS5kFi6lczsIOBld99tZh8EDgd+FUda9dq2bRvHHnvsgGnHHnss27ZtizXddhrMTPo2VspnrduQlm2H0r9J3qrPRa3SVN5pUe+lrKeb2TPAp4DbzWxdGDoOeNjMHgJ+AMxz95fry2o8Jk6cyD333DNg2j333MPEiRNjTdfbaDAz6dtYKZ+1bkNath1K/yZ5qz4XtUpTeadFXS0Hd78VuLXE9B8CP6xn3c2ycOFCzj33XFauXMnu3bvJ5XKce+65aj5LW9PnQtr+DunZs2cDweOit23bxsSJE1m0aFH/dJF2pM+FtH3lAMEHYfbs2fT09DBt2rRWZ0ckEfS5aG96KmsbSdLjgJP2eOIsqFSeWShrHTPNpcqhjSRpwE4DiI1XqTyzUNY6ZporU5XD6tWrmTJlCtOnT2fKlCmsXr261VmShNMxI1JaZsYcVq9ezcKFC/uvrhg+fDjnnnsugAbRpCQdMyLlZabloDs6Zah0zIiUl5mWQ9ru6IxLqbtdpTQdM/GJ467zZkvSNrSizDLTckjbHZ1x0YDd4OmYiU8cd503W5K2oRVllpnKIX9HZy6XY9euXf13dC5cuLDVWZOEysoxk7RB9aTlJ8tmzJjBsGHD6OrqYtiwYcyYMaNh685Mt5Lu6JShysIxk7RB9aTlJ8tmzJjB+vXrOe+88zj55JO54447WLZsGTNmzGDdunXVV1BFZloOEBx8W7ZsYePGjWzZskUHo1SV9mMmaYPqSctPlm3YsIHzzjuPa6+9lrFjx3Lttddy3nnnsWHDhoasP7EthywMWklzteO+SNqgetLyk2XuzmWXXTZg2mWXXcayZcsasv7EthyyMGglzdWO+yJpg+pJy0+WmRkXX3zxgGkXX3xxwx4vktjKQUSqS9qgetLyk2Unnngiy5Yt4/zzz6e3t5fzzz+fZcuWceKJJzZk/YntVhKR6pI2qJ60/GTB6tWrWbRoUX95Lly4kNmzZ7Nu3TpmzJjB8uXLWbZsGWbGSSed1JDBaKizcjCzK4DTgN8T/E70Oe7+ahi7GDgX2A1c6O6NybGIDJC0R2snLT9pVu3qr3xFEEdZ19uttAGY4u5TgceBiwHMbBLwJWAyMBO41syG15mWiEhbaeXVX3VVDu6+3t13hW/vAw4N/58F3Ozufe7+FPAEcEw9aTWCngcvSaDjUAarEVd/1XqsWaOu7DCz24Bb3P0mM7sGuM/dbwpjK4E73f0HJZabC8wF6Ojo6Lz55pv7Y729vYwdO7Zkeoo1NgZw9todXD9zTENjSdrGtMSSVJ6tSDMLZdqo9M455xwuvPBCjjzyyP7Ygw8+yNVXX811111X0zq7uroecPejSs4c5e4VX8BdwJYSr1mReRYCt7K3svkucGYkvhL4QrW0Ojs7PSqXy3k5jYytWrXKJ0+e7MOGDfPJkyf7qlWrYk0viTF39/HdaxoeS9I2pmXfJ6k8W5FmFsq0UemtWrXKx40b5yNHjnTAR44c6ePGjSs6ToeyTmCzVzkXu3v1AWl3P6FS3MzOAk4FpocJAzwDvD8y26HAs1VrqhbQ7f7tS/teku7ee+9lx44dHHTQQbz44osccMABvPTSS9x7772xH6N1jTmY2UygG/icu78ZCf0E+JKZjTKzCcDhwM/rSSsuut2/fWnfS9KtWLGCK664gueff55Nmzbx/PPPc8UVV7BixYrY0673aqVrgHHABjP7hZktB3D3R4F/AbYCa4GvuvvuOtOKhW73b1/a95J0fX19zJs3b8C0efPm0dfXF3va9V6t9GF3f7+7HxG+5kVii9z9Q+7+UXe/s/6sxkO3+7cv7XtJulGjRrF8+fIB05YvX86oUaNiT7vtH5+xcOFCzjjjDCZMmMD06dOZMGECZ5xxhm73T6Bafydg/vz5jB49mq6uLkaPHs38+fMBPepBkm/OnDl8/etfx8zo6urCzPj617/OnDlzYk9bj8+I2DueLklT6+Dx/PnzWb58OUuWLGHSpEls3bqV7u5uAJYuXdo/jx71IEn0+OOP4+4MGzaMPXv29P99/PHHY0+77VsOixYt4pZbbuGpp55i06ZNPPXUU9xyyy0alEyYWgePV6xYwZIlS1iwYAGjR49mwYIFLFmypH9AL+2/5yDZlv/Nht27d5PL5di9e3dDf7OhkravHDQoWV3+jt5WqnU/tXJAT5IhzXeke5nfbGhGL0fbVw4alKzOw99JaKVa91MrB/QkGfLHbxq7jeP+zYZK2mbModxjb/MD0mPGjOHpp59m/Pjx7Nixg6uuuqrVWZaI/OBxfswhP3hcrVtpzpw5/WMMkyZN4sorr6S7u7uoNSGSRPnfbAA4+eST+3+z4aSTToo97baoHCoNZkalsdnZLmr9nYD8oPMll1xCX18fo0aNYt68ef3TRZIs7t9sqKQtupUqDWZGB6Q3btyoAekEq3XweOnSpezcuZNcLsfOnTtVMUiqrFu3jj179pDL5dizZ09TKgZok5ZDtcFMDUiLSKPkeyBsSbovj2+LlkOlwUwNSItII6V5ADyqLVoO1QYzaxnoFBHJsraoHAYzmKm7ZEVE9mqLygEq/+i5fhBdRGSgtqkcpLL+y3gXD62fNCuDbyIyUFsMSEt1td4FnZXBNxEZSJUDtT8KWtJB+1dk6OrqVjKzK4DTgN8DTwLnuPurZnYYsA14LJz1vugPASWJfkc427R/RWpTb8thAzDF3acCjwPRJ0Q9WeoX4pJGvyOcbdq/IrWp92dC17v7rvDtfcCh9WepufTI7mzT/hWpjTVqINHMbgNucfebwm6lRwlaE68Df+3ud5dZbi4wF6Cjo6Pz5ptv7o/19vYyduzYkuk1KnbOOedw4YUXcuSRR/bHHnzwQa6++mquu+66puallbGk5aed9+/Za3dw/cwxichLK9JULN5YV1fXA+5+VMmZo9y94gu4C9hS4jUrMs9C4Fb2VjajgHeH/3cCvwHeWS2tzs5Oj8rlcl5Oo2KrVq3yCRMm+KZNm3zDhg2+adMmnzBhgq9atarpeWllLGn5aef9O757TWLy0oo0FYs3Bmz2Kudid68+IO3uJ1SKm9lZwKnA9DBh3L0P6Av/f8DMngQ+AmyuWls1Wa2PgpZ00P4VqU29VyvNBLqB4939zcj0g4CX3X23mX0QOBz4VV05jZHukM427V+Roav3DulrCLqQNoR3yuYvWT0O+JaZ7QJ2A/Pc/eU60xLJPN1xLklRV+Xg7h8uM/2HwA/rWbdIO3J3tXAkEXSHtIiIFFHlICIiRVQ5iIhIEVUOIiJSRJWDiIgUUeUgIiJFVDmIiEgRVQ4iIlKkYU9lbQQzewl4OjLpQOC/y8yuWGNjScuPYsmIJS0/itUfG+/uB5WZd6/BPJ2vVS8qPD1QscbGkpYfxZIRS1p+FGv8/i33UreSiIgUUeUgIiJFkl45/KNiTYu1Ik3Fkh9rRZqKNS9WVqIGpEVEJBmS3nIQEZEWUOUgIiJFEl85mNnR4d/JZvYHBbFPhH87zew9ZjbczGaZ2UkV1vfVEtOmmNmXImkdHP41M/u8mV0cxkeY2efM7B1l1j3SzE4zsz8M359pZl81s/3D9x8zs7lm1m1mZ5vZIbWViohIvBIz5mBmpSoqA9YCDwMdwC7g3cCfu/tLZrYJeCqcrw84CHgWeB14DzAR8Mi6ACYDW4A33X2mmf3FHtk/AAAJ6klEQVQVMB24Hfgj4LfAMe7+GTO7CngL2AQcARwFHEtwo94LwK3AT9z9lXAbbgXuB/YHOoE7CG4++VPgQWBf4CGgC9hJ8BOq97r7jeHyncAngXcBrxL87OrmCmV2tLvfb2aTgd3u/l+R2Cfc/T/Cdf4G+B1wKvCWu68vs76vuvt3C6ZNAaYAT4ZpHezuz1nwe5azwjJ+CvgBcDJwl0d+TzyynpHATOB37n6vmZ0J7Af8s7u/amYfAz4VbvsLwHp3f9bMhgOfLywX4EfuvqvMdpwGrAnzsztc154wNovgmDoV+GWY9z8n2M83uvvOEuv7lrv/jZnt7+6vhtNOzZdLuO1HuPuDZrYvMA/4g3Ddy4GvAGvcveh31M3sAODPCPbPvwLfAN4JXAtsD/MZLZfb3f3+cNkhHS/hMpk/Zmopm8GUS4q2fcifmZJlkqDK4U2CDTAGntCnAlvc/fhwvqnA1QQfoiXA8EjsEXf/WPh/DrgtXP56d+8Jp9/p7p81s01hBfBToCty8rgH2OnuJ5jZXe5+QiSPOQB37zKzCcAfA6cRVEw/Br7o7l3hvFvcfUpkuT3uPj2yrg3ufmI+DTP7DsHvcd8FvEZwgjiB4OT2V6WKjDaoOIETw23cWFAuHwf+pky5XB/mY3tYLtOBv3D3x8JyeR34T2B4mN6PwmkzgE8Avwb2lCiXXeG2XxZux4/DcjmU4K7Tz5jZDcC/R8rlbIKK4iHgveE++1d3fyQsl/VhfvcnqFQuJTgpfxN4Ang0XPYzwDjgZYL92UGZ48XdL2znL1vAkeXKhto/SyPTsO3ufqOZ/RNlPjPufmaJ7S+tljvn4ngBDwD7lZi+Afg3YJ/ItHeFO+AF4N8i00+L/N8T/t0HOB+4GfgccGc4/XngRuAZYN/IcpuBLwP/F7gOuAmYAywFrgByJfLYAcwl+Lb618Dfhfn+GnAOwYG3GugGPgssBq4Kl82Ff39Wplx+BrwZHlS58G/+/98BP43MOxXoAY4O54nGHon8nwMWEJyUpkWm58tmU/j3p8CwSPwegm865P8WrDO/LRPCbe8B1oXln4vMu6VguY2F+zyfBnB3mXK5m+Ck9f1wP0Vfz+T3fzjvIcB6gm9umwrycn/k/43AF4BV4X4bUa5cCvLy03BZC7fXCmL5cnlHuP6bCI6zywvyubVUeUbzF/lMlD1ewr/tfMzE8VlKxbbnPxvlPjOlppd7JanlcDBBE+r3BdNHAP8D2O7uL0amDwe+CDwC/Je7747E9gFmuvtPCtbzZeCj7n6RmY2PJPOsu79tZmOBT7v7neF4wAyCE/9rBLXyQ2Y2w93XldmGfQmagk8SdFmcRXDCWAX0AqcDHwQeA25z9z1mdogH3SdXEpw87iI46b2T4FtIH/Bp4DPu/lpBehvCZbry5WZm7yI4+RwFPOHufxROP83dbwv/73H3aWE5/QVwXJjH8zxoVT1PcDL9DHC4u78VLrcZuAo4nuBb90iCD8RUgm9FR3vYcorksYPgpPw5gpbhKILm7lqCb8FnAK8AvyD4tnM8QWX9v8IW1x3htJ5IuRxP8EE/HZjl7i8VpHkLcDBwqru/Hk7bh+B679MJPmz5cjna93bT9Lj7tPD/kwm6g+4FPudB6+7VMI+TgA970LwfRvANrxv4S4IWx3sJTgwTw/L5fIlyGRGW75fDfezAaIIugJfDNJ4Ly/hhYBrByfrvw3J5kDLHi7v/lZk9QLKOmeOAEQw8ZnYCRw3imBlN0KJbS3Cs/AkDj5lpwOjIMVO2bAg+S9M97B6sUi4HAP8U5rGjyrZvIPgWP9jPSyO3fV8PWos5D3o1vkHpz8zd7n45g5SYykHAzI4k6Efcn+Ak8e8EH6hnKF1xHk1Q+RRWnCOAS4AfUrrivMjdv1Uw/5eBj4bL9K8LeM7dfx9WnAvc/VulKk6CFtoBpSrOMJ9bKF1xbibo4ilVcR7n7j8zs2OBj4Vl8hrByfiDBK3Nj1Kij5igVfm+sNyisU+G6R8CvJ2PhX28n3T3u6P9zmbWRdCH/HMP+uOnhLFt4XLvAKa6+31h7NMEH958Pj9EcCLYVSafPw+XcYIKZUZYLq+7+7+F294F/NjdHw6X+0MP+qGPITjpjQjX7+6+OJznEIJvpL8kMnZA8M3zSIITVTS2i6BP+2GCCjw65rCHYOW32d7xiJeBCwm6Uj9lZocRPOBtGMHJ+eTINq81s+MIulqGE3SB7Hb3y81sJvASJcY4wrK9gOBz8BDBN+pXCXoBegm6iPYFtoZ53EFQgT4XfpY+yd7P0oHu/ndm9l7g5ehnyYKLVJYTdN8UfpYuAF5y91si0/Kfl8vd/aAyXzS/SjDusDbyeXlPeFy8290XmdlJXjCWY+E4RuSL5hPh66xwn90Qbufp4f59jKD7/HyC7srnwvUcCBwTbtMTBJX+/QyBKoeEqNJHPGOIMQiap0NdLomxhwk+VLsp7gd+kPJ9xHHFKuWl1nzWss6nwnL6PQVjA+4+18xWUn7sYHhMsZL5qRIbEcYamZf82Eh+bACCltijFI9pwt6xg1KxUsvl45WWK4zVk5d8vHC56DonE7SIjzOztb53DOQEgu7uPwJ+6+4XMVhD6YPSK74Xe/tCo698X2ilftLBxnI1xsqts9Z8bqqQXql1DrZ/vJ1jA8YGwr+Vxg6yHqs0PpDpWPh/2TGQIZ2T4j7p6TXIHVF5QL6dYxUvRmjnWGR6qQsxysazHgv/lrwQpU1iZS+2KXf+KfVq+UlRr/4dd3D0RBCZPqLNY8cQdJVEpw8HvtTmsckEl3FHY/sQDJ5TKZ71WIlj6BxgcZnjK3MxYHzkNTKcNhb4bOHylV4acxARkSKJf3yGiIg0nyoHEREpospBZJAseFjiNQ1a1/bwWnSRRFLlICIiRVQ5SNszszFmdruZPWRmW8zsDDM72szuDaf93MzGhbMfYmZrzeyXZnZ5ZB2zzeyRcPkl1aaLJN2I6rOIZN5MgscenAJgZvsR3J18hgePXX4nwaMnIHii5pEEd+Y+ZmZLCe5gXkLwqIJXgPVm9nmCR2MUTXf3HzVv00Rqo5aDSPDwxhPMbImZfRr4AMEzpe4HcPfXfe9z8De6+2se/O7DVoJryY8muPnqpXC+fyZ4OFu56SKJp8pB2p67P07w7f4R4DKCh5qVuwGoL/L/boLWt5WZt9x0kcRT5SBtL3xq5pvufhPwbYKneR5ie382dlz4JM5y/gM43swOtOBR8rMJnmtTbrpI4mnMQSR4HPgVZrYHeBs4j+Bb/9Lw0clvETzdsiQPHhF9McFD3wy4w91/DFBuukjS6fEZIiJSRN1KIiJSRJWDiIgUUeUgIiJFVDmIiEgRVQ4iIlJElYOIiBRR5SAiIkVUOYiISJH/DxPIkol6J3t+AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"jsp.boxplot('mathcent',by='school')\n",
"plt.suptitle(\"\")\n",
"plt.xticks(fontsize=8, rotation=90)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztvXt8HNWV4P891qPxU7Zl+dHGMsHCQEKYMDjkCQGHAOEXcMJvk5D1TLLJ7pCQ9ybOEmJmkgk2hpDHZpkdzzizyUwmzjKTGUB2AsaAcTAhEOwhCGMbY5FIttvGUttq2ZbcUqvP/lEluVvurm61+lHdfb6fT3+kurfq1qlbt+rUveeec0VVMQzDMIx0TCi1AIZhGIa/MUVhGIZheGKKwjAMw/DEFIVhGIbhiSkKwzAMwxNTFIZhGIYnpigMwzAMT0qqKETkxyJyRER2JqR9S0QOisjv3d/1pZTRMAyj2il1j+IfgetSpP9AVd/i/h4uskyGYRhGArWlPLmqPiUi54y3nFmzZuk554y7GMMwjKpix44d3aralGm/kioKDz4vIh8HtgNfVdVjXjufc845bN++vTiSGYZhVAgi0pHNfqUeekrFWmAR8BbgEPC9VDuJyC0isl1Etnd1dRVTPsMwjKrCd4pCVV9X1SFVjQM/Ai5Ls986VV2iqkuamjL2nAzDMIwc8Z2iEJF5CZsfAnam29cwDMMoPCW1UYjI/wWuBGaJyAHgm8CVIvIWQIE/Ap8umYCGYRhGyWc9fSxF8v8puiCGYRhGWnw39GQYhmH4C1MUhmEYhid+9aPwPWvXrmXbtm0AXH755dx6661Z5RnlRy73eu3atbS3txMKhVIe5wf8JGM1PjOJ1xwMBlm0aFFS+8kmr1h1ZYpiHPT39+eUZ5Qfud7rcmgHfpHRL3IUk3JpV6KqRT1hIViyZImWwjN7xYoVAHz3u98dU55RfuR6r8uhHfhFRr/IUUxK3a5EZIeqLsm0n9koDMMwDE9MURiGYRieVIyNIlfDnJfRqJJl9Ivx0I9ylLId5GrgrEYZC0G+22OxJwwUanJFxSiKYXIx8hTbMOQXGf1iPKxWOdKdz09GzHKQMd+U+7NWiHtTMYpiWDsOG3my1Za33nor7e3tQOGNaH6SMbHMUn7x+VGOUrYDLzlMxsKT7/aY6zM/nvOlkn+8cpiNwjAMw/DEFIVhGIbhScUMPRUCvxgPx1PmWA1zhTC451qPuRjm/ORt7Bf8OGGgHO9Lucs/HkxRZMAvxsNil5dvg3sh5Khkg2q+8Ut9+EWOXCl3+XPFFIUHfjEejqfMsRrmCmFwz7UeczHMFdt4WA74ccJAOd6Xcpd/PJiNwjAMw/DEFIVhGIbhSamXQv0x8AHgiKpe5KbNBP4FOAdnKdSPqOqx8ZzHDJz5werRKDa5Tlwo94kofqPUPYp/BK4blfZ14AlVPQ94wt3OC/39/VVrjMonVo9GsfFqc+nyit1OK/m5KPWa2U+JyDmjkpcBV7r//xOwFbhtPOcxA2d+sHo0ik2uExfKfSKK3/DjrKc5qnoIQFUPicjsVDuJyC3ALQDNzc1FFM8wyp/h4RJg5O+KFSvKLoifURz8qCiyQlXXAevAWbioxOIYRlnR3t5O25691DTOJa7OCHTbnr0llsrwK35UFK+LyDy3NzEPOFJqgYzKoNxDYHuRi9dwTeNcJt34qZHtvg0/Lph8fqIS20G6HiKQl+vzo6LYAHwCuNv921pacYxKolKNjVDZ15ZvKq2u2tvb2bdrD80NTdQPOWkDB8N0RrryUn6pp8f+XxzD9SwROQB8E0dB/KuI/FegE/hw6SQ0KolyD4HtRTV7DY+VSm0HzQ1N3PGu5Nflqt/8Ii9ll3rW08fSZL23qIIUATMeGoZRrvhx6KkiaW9vZ+/uNuY1CDVDju197+62EktlGIaRGVMURWReg/AXV9SPbP/oqYGMxwz3RNIZqCrRMFcOFNp4WC1YPZYHpih8Tnt7O6/sbqPOvVPHDrVxpCd5n0ozzJUDzvTS3UjjdFTjALzUdQgN92Q40kikvb2d3Xv2MaNxIXF1PqIOdw1yLNxRYsmMRExRlAGzp8PNV52+Vfc/GRv5v1INc+WANE6n9salSWmxDVtKJE35MqNxIe9bdkdS2mOtq0okjZEKUxSGYRhFJtOQst8wRWEYhlFkHL+HvQRqnOG2gQPH6ew9VGKp0mOKYoyY8e1Myu3ryMg/Xs/FyZMnmTx5csq8UCjEhLqmYovrC5qnzWPl2z89sr362b8voTTemKIYI+3t7by66yUWNNRQN+S4QJ46uIv9kaESS1Y62tvb2bWnjQlua+o80saJo6WVySgu7e3tvLznVSbOamYA5yv5te4o/d2dnFU3geigMm1WMzE370D3AL3dnQTqhMl1pZTcyAZTFDmwoKGG//HOaUlp33mmt0TS+IMpM+HS62Vke8fD3nEarWdWeUyc1cziZcnLx+xtvRuNHGDarGbetuwbSXnPtd5FNLK/mCIaOWKKwigJ7e3t7NzTRn0jDLo6ZW9XGwPh0splGMaZVLWisLH10lLfCGcvS15k8UBrvETSGMbYqZbQPFWtKNrb29m3exeBmhoABkL76YwcL7FUhhfl8GCGQqGRj47RMiamlVL+UCjEUO/xpNDiQ+FDhAZPFE2GSsCZvfQqzVMXUB9zjC37dr2ac3np2g5knhQwi0DO581EVSsKgOaGqay8/K0j26u3PV9CaYxMOB7RO2FWAHBCoLTt2VlaoUbR39/vem3PRNUZV2vbs3skv23PnlF5e0oip5Efmqcu4BtLvjayfdf2e3Muq7+/n327XqF52lzqY05ve+BAhM7ew0igFo3GaG6YTf2QYw8cOHiMzsgRJFALZ5miMIzTzApQ88EFI5tDD/nPICqNM6m74ZqR7cGNm0flXZ+Q93BRZQMnLtixut4zFi4KNk3zOMooBs3T5rLyHZ9MSlv925+wP9pNc8Ns7njn8qS8Vc+sZ/+pwk4znJB5F8MwDKOaMUVhGIZheOLboScR+SNwHBgCYqq6pLQSGYZhVCe+VRQuV6lq93gK8Jol4xeKKaM5umWPV12FQiGok7TH+oHM8k8Z03GJs7ZSllffmFf5vWaP5budjicEySxm5lWWXIhGo3REu85Y+rQj0sVkiY67fL8rinHjTIHdzcKGBurdkBv7du/OcFRxGV5zYu50YULcmQnzSoFWv2tvb2f37jZmzIC467Jw+HAbx44V5HRlzek1JxpQddrOS10hNBxhcl091E0qsYTeOPK/woTGOcTVGWXe2dVDPPw6k+tq0yoK57i91DQGiaszdbxtz96R/Jf2vAq1zgybXV39DIUPMKluAlKfsric6e/vZ8+efcxqXAjuWhV79uzL70lc2tvbeXX3PuY3NFM75JyrLzTAwUgnE+qFeFRZMK2ZupiTd+rAAPt7O5kQEMjzdfsRPysKBTaLiAJ/r6rrEjNF5BbgFoDm5mbPghY2NHDHFVeObK96amueRR0/c6cLH7/ydNCbn24dLNi5ZsyAa65O/hre/Lh3yI1yoBA+FtLYQO2yy5PSYq3boLc8Foua0DiHiTf+WVJa/4afQa+3C3xNY5ApN352ZPvEhr9NyDubqcu+OrJ9vPV70BvKk8TJzGpcyLIb/nJku3XjnUBh7vX8hma+cPnKpLT7tq3mUP9+Fkxr5mtvSw5Bcu9zd3Ew6o8Zd4FAgAVnTeeOd304KX3Vb35BfXD8PT0/K4p3qWpIRGYDj4nIHlV9ajjTVRzrAJYsWVL+bzlj3Dhfwi/DrEmc9rF42fOYzMMzhZPXyJ329nZe2bOPppkLEbe38UqBehuGjxWFqobcv0dE5EHgMuAp76OMqmfWJGo/eOHIZuwh72FGR7nsQhqnoeqsHPhS1wE03MvkugDUTSyouEbuNM1cyM3Xn+5t3P/wnSWUprLxpaIQkcnABFU97v5/DfDtEosFOF+ZJyOxM6LF7o/E0JN707rf7927F2LKj54aGDnmUI9ynBDBYLBI0ntTTONhLnIM5+V7TFgap1G77O1JabHWZ6E3vREwGo2i4YEzlj7VcA/RuAKpx/+rkWg0SrS7g+da70pK7+3uYGjwFAODHWcsfXos3IHGByhnvIbHQqEQJ3tPJK1B0dF7iEFiFDASR874UlEAc4AHRQQcGX+uqptKK1JmYrEYr+5qY37DBGqHHEtx38GdHIzEiVHr28oepr+/n92725ieYOjeXSCjeiY52va0wSzBMVVBW/dL0K1MrptcFcbDbMk8dOZvg3slczoO1HzqY87TP544UKXEl+8uVX0N+JNSy5GKYDDIKe1JuR7FwVMB5p0V5wvvTv4kuO/pKIdOBZg9cYi/uOL0W+5HTw0w1Se9iWGmz4Arrzlt6N66uUTmn1lCzU3JBoKhBwYhUhpxRhMIBBiaNonaG5cmpcc2bCHQ20dfkeQ4PbOpibg6921n11Hi4S53ZlPpFUUgECDQsCDtehSTpzXzvmV3JOU91rqKk72dxRSzIDRPnc83LvvCyPZdv7sPcN4jA/HjZ6xwtz96pOgyZoMvFYWRHRYm3QCY0NhE4IaPJqVFN/4L9NqcZyM/mKIoY9rb29mzu4069y6GD7XR1VNamQzDqDyqWlGEQiFO9hxPCi3e0XOcyRRmTnguhEIhjkfg/idjI2lHeqBfHSN403S46b01I3kPPDGU0Rg8oUgRvjIapW3qacUQjUaJd3eyt/XupPS+7k4mxAdzss9Go1G6wx0jvhMA3eEOBgYneU4AKWZP2/GI3p8UWrzj+H4mhybnNEnFKe8wq3/7k6T0jt7DDDIEZ41dxkzro2RDVSuKSqW/v589u9tonAHqGqW7DrcRPgZ19ZNxoxEURY6X97QxqREGXFPHH7ra6AvDWXWTTVEYBcGJdLCP+jpHPUVeH+Tw0Y4SS1U6nDUu9tDc0Ej9kPMg7ts1tjVQqlpRBINBBhg6Y+Gieh8ZmIPBIMekm5uvOn2r7n8yxox53jI2zoAPvDe56/DLJ+L0nsxNjlzjUU1qhAs/kCzH7l/GifemOSAD0WgUuuPJa1B0RwkN+GeacTUSCASQhrNZvOzrSel7W+9GIwdyLrNhWvMZntmzmjJ/YcyduZCPX3PaQP7Tzas89s6dQCDAgvr5ZyxcVB/M4dN/uLzArLTrUeRKc0Mjd7x72cj2qqdbx3R8VSsKL9J1X0OhkA9CgBWf9vZ2du1uo2EmuB8l7CrB1FnDMIqPKYo0OHOgXyZQ6y5HePA1OiP9SOCsnMYJK4GGmfDO605PnX1mU/GnzgYCAWIN8TNWuAvOst6EYRSKilAUXV1d4zbWpKK5YSK3X754ZHvNtr3sP2VhpQxvhr22E5c/1fBRQoNOBFrtjSQtf6rhMKHBmG+mO4dCIYZ6TyQFAhwKhwgNTinr4b1MkyumM6tksmWDY+g+wqpn1ield0SO5GzozpaKUBTRaNQNJT7ztLHGZ6HEDSMTjvPcHqh1HsuXurrQsHeUVyN7+vv72bt7H/OmN1MTdxxfjx8a4FBPJzX1wnQL65WWilAUAAsbZnLH5deObK/a9mgJpTGqGcdrewp1N1wzkja4cTPBpjkAhOtqqbvh+oS8hwk2zQZAGhupv+HGkbyBjRuKJPVpgsEgx+pOnBFmPNhU/vGr5k1v5palyV7g67as4kifP8KFe+GEEp/JHe9cnpS+6pn17D91tKDntjWzDcMwDE8qpkdhGF54zWIzfw7D8KbiFUUoFKKvpydpVbuOnh4mgW8Mc473tSatane4RzmplesbEI1GGQrDgdZ4cnoYauLjX+N3NMPrTlDndKLbujuh+4S75oQP4zpnSSgUIt7b68R2SiAePkLf4CA6GHNWtEvKe51ofCjn88V6Tzqr2rnEwgcgPlitkwGrgopXFIYxwqwp1C7705HNWOt/QKRwS84aRqVQ8YoiGAwyCGesmV1XgC/1aDTKgWic+55O/iI+EIkTIwppZlUEg0EiEj5jzeyGDN7X5Yxj8I1x9rJkM9mB1jg1vQFiRQvUXd4Eg0GO1p2VMnrsxN5j9E9rTLlmdqA3zKkcz9dT13/GmtmBAq2ZbfgDM2YbhmEYnvi2RyEi1wE/BGqAf1DVuzMcUnICAe+Fi7CvZMMwyhBfKgoRqQH+N/A+4ADwvIhsUNVdpZWs/IlGowwMwObHkz3Mjx0D1SiTcogsOxwKfcfDp8s8HgbRaLrRtowy0q3OinaJdCvReJRqm6aUOVy7mZHHizNs3MF921YnpR/o6SA6dIr90Q7ufS55ze/9vR0MMlgVS/P6UlEAlwH73CVREZH7gWWAKQqj6nDWEN+DNDahri5+qSuMhruYXFdnisIoOH5VFPOBRFfJA8DbSiRLRREIBJg8OcY1V0tS+ubHlZMnA0As9YEeBINBYrXdXHr96TJ3PKwMHs+tPCfw31DKNbMDkQAx4imPC4VC0NtH7KGE8C3dfRURglwamwjccFNSWnTjA9BrSxrmg0AgwLyJC/jC5SuT0u/btppD/fuZH1jA196WvOb3vc/dxcGo/z2684FfjdmSIi1prEREbhGR7SKyva/Pxv4NwzAKhV97FAeABQnbZ0Py+qSqug5YBzB37lwL6WoQDAbpro9S+8ELR9JiD+22EOSGMU78qiieB84TkTcAB4Gbgf9cWpHKh2g0SnjAWdEukfAxiGs07VKow4burZtP692eY6Bxnw3ddEedFe4iA872oFKICNFOuPAosdZtSeka7iEaB5iU/5MaWREKheiN9HH/w6fX0z4S7uBUzLknxyN9SavaHT7awZAO2C3LEV8qClWNicjngUdxpsf+WFVfLsS5OiPH+eIjvwZgzpRJdEaO0+Kjd6KRzMSJE0fWGWmPODOAFl2wiEWLFo3MCDKMaqQz0sWq3/yC1086dqs5k6fTGelCAuOfJehLRQGgqg8DD2fccRyMvHDcF0x9cAEtQcrmpeN8VcEDT5yO29Pl2jbTr5md3sAcCASYNDnGldecNhFt3azMm+sfzRkMBvnud78LnJ4iOno7Xzje4xOpXXZ5UnqsdRuB3n7ziikhwWCQs2oHufn60+tp3//wncyc7bwUIzWDZ6yZ3X2is+hyetHZe4gvPuFMuZ0zuZHO3kNIoCansiZOnEjQfZ8NtDsvgfr5jbTMb3QmeYxzcN63iqIYDK8aNvqFk5hmGIaRb874SD17Ki1MdV7qqSf1eZLpA2rgYNe45K1qRVHuBINBAtLNTe89/RXywBND9PTlNi3VMIzikO4jdcWKFQwciJRMrnSYokhDKBTiZKSPNdv2jqR1RPoYZAL7o0N855lejpx0hnxmT65hf2SICeUbrTojoVCISASe2XS6Dxs5CjLkM0P3GDltsH42KV3Dva7BOr1vuYZ7iG3YgkZOACANU9BwD9SVh6vuUPgwfRt+TDziLLeqsUFomlZiqQw/YopijNTW1nLe4jcCMOh2G8+av4jz5rsOX9pfSvGMIpFkVO91jepN86BpnrMGSimFywJHfkfBt/c6wxKLLlhcNvY5o7hkpShE5B5VvS1TWiURDAYZ0FPcfvnikbQ12/ZSP/9cz7HAvoOFXbu2VASDQbSmm3ded9rQ/cwmJTinfHsTMGywDlC77O1J6bHWZwn0RtO+8DONCYe7Xi+YzPmgmJMCjPIn2x7F+4DRSuH9KdIMwzDyRne4g9aNdxKJHAZgMBZlVlNLiaWqPjxDeIjIrSLyEnC+iLQl/P4AtBVHRMMwqpGJEydywQUtzGqqAxkAGeCCC1pGhvyM4pGpR/Fz4BFgDfD1hPTjqlqZYyw+5EgPrN3ozGKaMcXZnjGvxEKVGU7AwOPO8qfDdB8nGhfAH7MQQqEQ2tvLwMYNI2kaDhON5zBfsoAMhQ/Q81PndVDT0MRQ+ADUTaC/u5O9rXcTjRwBINAwm/7uTs6qm0BvdyfPtd7FyYgzJDe5YQ693Z0E6lKFdXMo9+GxUCjEyeMnuOt3942kdRw/wOTQFM/jOnsPs/q3P+H1k84rds7kmXT2HkYCuZmUnYk5vax6uvW0HJEwkyX7ZYA9z6yqESACfMxdI2KOe8wUEZmiqv7yYKlARs+3njFvETPmlY9ToFFZjG6Pi5omQtN5nDx5kslubJh2N7TKubMCMCt13tmz6mFWi6PAjRESe0sD7d0A1J/dQAsNeXGcy5VsjdmfB74FvM5pdxAFLi6MWMYw5hSYH5yAgTFql/3pSFqs9T8IRAZ9M0MpGAwSrquj/oYbR9IGNm4g0NvrGxm92uMwY8lbsWIFh7uy/7ItJ4LBIANDfXzjsi+MpN31u/uoD6YPODVcv5DOce5YbnJoHXe8e9lI2qqnW6kPNmVdRrZ9mS8D56tqeGwiFofBwUE6eo6yatujI2kdPUeZxBDBYJCOSIRVT23l8Alnvnt0aIiWEsz9PxRRfvTUAOETzmfBwBBMLe9JQ0YaNHyUwY2b0chxJyE2CE1zSiuU4cnBSCf3bVtN1wlneKxpyhwORjqZUJ9+eKxayFZR7McZgio7krtyTne5ZfHiohvEEmOxHHHlWLzYgtnlnW534aLIKWd7MF6QyLJeJLat9l7n42TRBRfavfYxzvPpfLXF2p3hsUnBes4LtmQMq9F5fD93bb+X1/sc20x0KEoL5xVc5mKSraJ4DdgqIr8CosOJqvr9gkg1Rurq6lg4fSZ3XH7tSNqqbY9SF5zj2ZUrBAcjce57OkrXSadlNU2ewMFInPPeWN6GuVxwvJ5h9y+Tn7K+MEg8mndvz6QXdJ4iy2o4Qqx12yjv6wg0pe8KZho+yDfxcBfRjf9CPOIEg5vQMJ14uAvqaomHX6d/w8+IR465eTOIh1+Hpul5l6OcyWQ4P3VgIOVxiQpmoN0ZQms5/7yCfRR0Ro6w6pn1vH7SuZ9zJs+gM3KElvkz8n6uRLJ9VjvdXz1VsZR4biS+qGJuI5nkem3blL7Ck+8XdGrv6yA0BX1zP5M9rJ2Xx6KmmdA0M9mI3Bt286ZD03TfyF/uFHNmVvLoiDMjqn7+DFrmzyj4/cxKUajqXwOIyGRVPVlQicqYYvde/E4gEGDCtBgXfiDZXWf3L+PEewMM+TxwoddLwC9kK6Nf5Teyp5Tvl2xnPb0D+D/AFKBZRP4E+LSqfraQwlUTh3uUn24d5OiwoTsGDePwlQgfc9afGLalNkx10vwUr24gDAda4wy61q+6Bidt4vjXWTGMM4hGoxwa6GDdllVJ6Yd6OjiuaZZ9NIDsh57+J3AtsAFAVV8UkSsKJlWVkdht7HaHrM6/cFHO3cnEIZPek055TXMX0TTXdTyj9J3C1MM6i6DJdRAq4oRQDfcSa30WjTj1Ig2T0XAvZD970DAqmqztiaq6XyRpmthQun3Hg4h8C/gLYHiljW+4q91VLPnuUmYaN929u5vNjyvH3d7G1Klw7BjUF7G3kUnG7u7izMROrbDOhqbC2ZWcqbMPo5FeJyEWg6bZBTlXOXAs3MFjras47sZzmtowl2PhDuaOI6bT4aMd/M9//zwAM6fO5fDRDgKBALMnLeCWpXck7btuyyqmzvNRV9uHZD09VkTeCaiI1ANfBHYXTix+oKo2mJoFXT3w41ZHZ0+f4mw3egxZJb78Trq9jblzFzHXR72NYlJsO0Ty1FlHUy+64IKqnTqb1B57nZlFc5vqmNuUe0yn0d7jDXPqaJhjXuDjIVtF8Rngh8B84ACwGfhcoYQysmP0A9E4bxGN87y/hDPNDDp0uLtQ4hoUf+qs3ymEgdZr9bjjh1JPczW8yXbWUzewvMCyJPJ5Efk4sB34qqqe4bcuIrcAtwBMnTq1iKL5h2zCKZQt3crQA4MQcYPbNAh0a9Gd5/yA44/SRXTjA0npGu4iNBhNc5Rh5I9sZz39E/AlVe1xt2cA31PVT+VyUhF5HJibImslsBa4EyeW1J3A94AzzqOq64B1AHPnzi1RqCyjEKR0nJu1CGaZP4phlIJsh54uHlYSAKp6TEQuyfWkqnp1NvuJyI+AX+Z6HqM8MX+UZJxV+KYTuOGmpPToxgcINjUW5JxD4RAnNvwt8YgzFKmxAWhanOGo4tJ1tIP7H76Tnt7TixrNnF2YRY3293Zy73N3ceSkEwdq9uQ57O/t5DyqYxGlbBXFBBGZMTwEJCIzx3DsmBCReap6yN38ELCzEOcxDCM1yQZ358U4vJ62X0iU5dhxx+5wfoEWNUosc9CNA3XW2fWcR/UsopTty/57wDMi8m/u9oeB1YURie+IyFtwhp7+CHy6QOcxDCMF5dCjK6aM5VAfhSZbY/ZPRWQHcBUgwE2quqsQAqnqnxeiXMMwDCM3xjJ8tAc4NnyMiDTbCneGYRiVT7aznr4AfBNnhbshnF6FrXBXgfQcg62blRPD6+3EYF6q+WmjOHEUtt3vTD6bNM3ZDljMJsOoCLLtUXwJH69wZ+SHJCOm67V94YWZ13MY7fjXPHsRzHY8vaNV5ultGJVIxa9wNx46I/18edNLAMyZHKAz0k/L/BILVUBy9Rr28oT9Q5d5ehtGueOpKETkK+6/vl7hrhCM/kqun38uLbYAkWEYVUimHsVwbIxUK9xVtDd0RYfHMAzDGAOeiiJhZbsPq+ovEvNE5MOFFMwwqhENh4n+7J8BkIYGNByGOpsVYJSWCZl3AeD2LNMMw8iRRYsWcfEFFzBZhMkivLmpiYsvuICJEyeWWjSjyslko3g/cD0wX0T+V0LWNPD5gseGUWZ4TQoId9mEQ6N0ZLJRhHBCfd8I7EhIPw7890IJZRiGYfiHTDaKF4EXReTnqjpYJJkMwzAMH5GtH8U5IrIGeCNw1nCiqp5bEKkMwzAM35CtovgJTgiPH+AEBvwkThgP39AROcqqbY9y2I09ER0apCU4p8RSVRaRo/DMJuXkcHiPQZhfgipeu3Yt7e3tIz4ua9euTXIWTEv3CWI/fdr5v2ESdJ+oyhXzjOqjMxJm1dOtvH7S8ZuOxmK0zG/K+vhsFcVEVX1CRERVO4Bvicg2HOVRcgKBAC0XXgjAQPsJAFoWX2jOcXkkKbzHCecFvdgN71EqxjIbaLQD5aJZzbZinlEVJLbxgfZeAFoWXzCmtp+tojglIhOAV0Xk88BBYHbWZylb9uxVAAAY60lEQVQwTU1NZ4SaMAe5/OKnmPxZ9R7SHFNq2Q2j2Hg9u5/97GezKiNbP4ovA5OALwKXAn8GfDxbQQ3DMIzyJVtFocA/AxuAJcBi4Ee5nlREPiwiL4tIXESWjMq7XUT2icgrInJtrucwDMMw8kO2Q0/rga8BLwHxPJx3J3AT8PeJiSLyRuBm4E1AEHhcRBar6lAezmkYhmHkQLaKoktVN+TrpKq6G0DkjIlTy4D7VTUK/EFE9gGXAb/N17kNwzCMsZGtovimiPwD8ATJYcYfyLM884FnE7YPuGlnICK3ALcANDc351kMwzAMY5hsFcUngQuAOk4PPSmQVlGIyONAqkU0V6pqa7rDUqSlDGeuquuAdQBLliyp6JDnhmEYpSRbRfEnqvrmsRSsqlfnIM8BYEHC9tk48aYMwzCMEpHtrKdnXUNzodkA3CwiARF5A3Ae8LsinDevJHoOr127ttTilJS+MOz+ZZwX1ju/3b+M02eBUA2jrMi2R/Fu4BMi8gccG4UAqqoX53JSEfkQcB/QBPxKRH6vqteq6ssi8q/ALpww5p8r1xlPtobAKG/uXscj+g1Ni6DJPKINo5zIVlFcl8+TquqDwINp8lYDq/N5vmJz66235uQ9XGkU2pt73759rFixgu9///uce67FpzSMQpHV0JOqdqT6FVo4w/Dinnvuoa+vjzVr1pRaFMOoaLK1URiGr9i3bx8dHc63SkdHB6+99lqJJTKMyiXboaeKZdjwDM7wyKJFi8pq2MhL/sS8rENxlwn33HNP0vaaNWv40Y9yjipj+JScQ8qn4VBPJ+u2rCJ84nUAGqfM4VBPJ1PnteRF3kql6hUFlL/h2Uv+cr+2dAz3JtJtG5VFPtpx4gSKI+0DAEydV8/UeS02uSIDVa8oyt3w7CV/uV+bFwsXLkxSDgsXLiyhNEahyGf79VOo/HLDbBRGWXLbbbclbd9+++0lksQwKh9TFEZZ0tLSMtKLWLhwoU2PNYwCUhWKYu3ataxYscK8pSuM2267jUmTJiX1Jir1Xmu4i+jGBzj1sx9z6mc/JrrxATTcNZJv0QCMQlJVNopKNexWKy0tLTz00EMp8yrpXid7uPc4aU2N0NSYlFdJ12z4i6pQFJVq0K12UnlmV+K9zsYIW8kTF4zSUxVDT0ZlYp7ZhlEcTFEYZYl5ZhtG8TBF4RMq1QhbKFJ5ZhtGOZE4AWHFihVJz7xXXimoChtFOWEGyewwz2yjEiiXqAqmKHyCGSLHhnlmG+VOOUVVsKEnoywxz2zDKB4lURQi8mEReVlE4iKyJCH9HBHpF5Hfu7+/K4V8hv8xz2zDKB6l6lHsBG4CnkqR166qb3F/nymyXFljnrClJ5Vndq7k23hYyZMTyskIO1by/Vx3Hj/IXb+7jy9t/Su+tPWv6Dx+MA9SpqdQ76WS2ChUdTeAiJTi9HnDT8amasTLMzsXCnE/K7WNlIsRNhfyJX+i1/xAewyAlvPPK3hI80LUvx+N2W8QkReAXuAOVd1WaoFSceutt/KRj3yEu+66i49+9KOlFscYJ/k2HvrJEJlvEtv+N77xDWbOnJmUV87Xnk/5SxHWvFD1X7ChJxF5XER2pvgt8zjsENCsqpcAXwF+LiLT0pR/i4hsF5HtXV1dqXYpOOvXr2fnzp2sX7++JOc3jFJhbb+6KJiiUNWrVfWiFL9Wj2Oiqhp2/98BtAOL0+y7TlWXqOqSpqamwlyEB+FwmM2bN6OqPProoxw9erToMhhGKbC2X334anqsiDSJSI37/7nAeYAvYzOsX7+eeDwOQDweL+iXVbkbCL0oh0kBXjLmIn85GLq9ZCxm28+VSn5mSkGppsd+SEQOAO8AfiUij7pZVwBtIvIi8G/AZ1TVl58rW7ZsIRZzDFSxWIwnnniioOebOHFi2RsJ01EO1+YlY67yl+t1F7vt50o51G+5UKpZTw8CD6ZI/3fg34sv0dhZunQpmzZtIhaLUVtby3vf+96CnavcDYRelMO15duD1u/XC94yFrPt50o5tKtywldDT+XE8uXLmTDBqb4JEyawfPnyEktkGMXB2n71YYoiRxobG7nmmmsQEa699tqkKYKGUclY268+TFGMg+XLl3PRRRdV1BeVnwyt5WDoLgfS1eN47rUf2r6f2mql40eHu7KhsbGR733ve6UWoyD4xQjoFznKnXx7Uvup7VsbKTymKHDmhd91112sXLmy6rvRfjIA+skgWc5tJF09+qVuc6Xc5S8nbOgJ8zI1MmNtxKhmql5RmJepkQlrI0a1U/WKohy8TAuBGYqzp1rbSL7Jt4d7sSnEpIB8ylFIql5RlIuXaSEwz9XsqOY2km8K4eFeTPwif7HrquqN2eXgZVoI/GQo9juV0Eb8YIz3Cs1fDu3RL5MCMtXVjh07WLlyJWvWrOGSSy7JyzmrvkdhXqZGJiqhjfjFGO8XOSqZ1atXE4/HufPOO/NWZtUrCvMyNTJR7m3EL8Z4v8hRyezYsYMTJ04AcOLECV544YW8lFsxiqIcvEzNk7S0jCf0tB88kXPFL8Z4v8hRyaxevTppO1+9iopRFMPkYuQZ9jIt1pdiORjtKpVc677YbSSf+MUY7xc5Kpnh3kS67VypGGO23w1hUB4yVjLlYDAtBH4xxvtFjkog3eSEKVOmJCmHKVOm5OV8FdejMAwjGb8Y4/0iRyWQblLAypUrk7b/8i//Mi/nM0VhGBWOX4zxfpGj3PGaFHDppZeO9CKmTJlS3tNjReReEdkjIm0i8qCITE/Iu11E9onIKyJybSnkG8bW3TWKTaHanF+M8X6Ro5zJNClg5cqVTJgwIak3Md6JNKXqUTwGXKSqFwN7gdsBROSNwM3Am4DrgL8VkZoSyQiY4dkoPoVoc34xxvtFjnIm06SASy+9lE2bNqXsTeTatkq1ZvbmhM1ngf/k/r8MuF9Vo8AfRGQfcBnw2yKLCGQ2fvrB29UoDsW619VqcDeyZ+nSpTzyyCMMDQ1RU1OT1aSA8bYpP9goPgU84v4/H9ifkHfATTsDEblFRLaLyPaurq4Ci5ga8zKtHuxeG35h+fLlDA0NATA0NFSUYbyCKQoReVxEdqb4LUvYZyUQA4afPklRlKYqX1XXqeoSVV3S1NSU/wvIgHmZVg92rw0/cezYsaTtnp6egp+zYIpCVa9W1YtS/FoBROQTwAeA5ao6rAwOAAsSijkbCBVKxvFgXqbVg91rw0/cc889Sdtr1qwp+DlLNevpOuA24EZV7UvI2gDcLCIBEXkDcB7wu1LImAnzMq0e7F4bfqKjo8NzuxCUykbxN8BU4DER+b2I/B2Aqr4M/CuwC9gEfE5Vh0okoydLly6lpsaZkJWtQckoHuFwmK9+9atjGibat28fH/zgB3nttdeS0pcuXUptrTPvwzyKjVKzcOFCz+1CUBJFoaotqrpAVd/i/j6TkLdaVRep6vmq+ohXOaVk+fLlDI+YqarNC/cZuRif77nnHvr6+s7oyptHseEnbrvttqTt22+/veDn9MOsJ8PIK7kYn/ft2zfShe/o6EjqVZhHseEnIpGI53YhMEWRI+vXr0/6yqx2A6efvNhzMT5nMhCaR3Hl4Ke2mguFCiXuhSmKHDED55n4xYs9l3uTyUBoHsWVhV/aai4UKpS4FxUTZrxQpPPItZDJyfjJoziXe7Nw4cIk5VAMA6FRGvzUVnOhUKHEvbAeRQbSGUWXL1+eNLxhQxL+IRfjcykMhIaRC4UKJe6FKQoPzCO3PMnF+NzS0jLSi1i4cCHnnntuocU0jJwoVChxL0xReOBlFF2/fj0iTsQREal6Y7bfyMX4fNtttzFp0iTrTRh5Y7zhvdORKpR4ITFF4YGXUXTLli1JgbnMmO0vcjE+t7S08NBDD1lvwsg7+Taee4USLwRmzPbAyyhqxmzDMDJRzkbzRKxH4YGXUdS8dQ3DqBZMUXjgZRQ1b13DMKoFG3rKwPLly+no6EjZY/DKMwzDqBRMUWRg2Cg61rxKJDH0wfB2NmOww8cBrFixgkWLFlXM2K1hVAOmKIwxk8vsjXINl2AYhimKcZEuvEelkmsvoFxCJlTb/TSMbDFj9jjIZc0Dw7/Y/TSM1JRqKdR7RWSPiLSJyIMiMt1NP0dE+t1V70ZWvvMjFt6jsrD7aRjpKVWP4jHgIlW9GNgLJMZMaE+18p3fyGXNA8O/2P00jPSUainUzaoaczefBc4uhRzjwdajqCzsfhpGevxgo/gUkLg29htE5AUR+bWIXJ7uIBG5RUS2i8j2rq6uwks5iqVLl1Jb68wFsBAe5Y/dT8NIT8EUhYg8LiI7U/yWJeyzEogBw/38Q0Czql4CfAX4uYhMS1W+qq5T1SWquqSpqalQl5EWC+FRWdj9NIz0FExRqOrVqnpRil8rgIh8AvgAsFxV1T0mqqph9/8dQDuwuFAyjgcL4VFZ2P00jPSUxI9CRK4DbgPeo6p9CelNwFFVHRKRc4HzgNdKIWM2WAiPysLv99M83I1SUSqHu78BAsBj7uI/z7oznK4Avi0iMWAI+Iyq+naeYrWF8Kh0yuF+moe7UQrEHfUpa5YsWaLbt28vtRiGYRhlhYjsUNUlmfbzw6wnwzAMw8eYojAMwzA8MUVhGIZheGKKwjAMw/DEFIVhGIbhiSkKwzAMwxNTFIZhGIYnpigMwzAMTyrC4U5EuoCOhKRZQHea3dPl5XJMsfP8IofJaDL6SQ6TMfe8haqaOaqqqlbcD9g+1rxcjil2nl/kMBlNRj/JYTLmLy/dz4aeDMMwDE9MURiGYRieVKqiWJdDXi7HFDvPL3J45flFDq88v8jhlecXObzy/CKHV55f5PDK84scaakIY7ZhGIZROCq1R2EYhmHkCVMUhmEYhielWuGu4IjIW1X1eRF5EzCkqnsS8t4GxID9QBhn7e5+Vd2cpqzPqer/HpV2EXAR0O6eZ56qHhJnyb5lwIXAH4B/A64HHteEZV8TyqkDrgPCqvqMiPwZ0ACsBxYA7wBmAK8Dm1U1NI5qMQzDGDNlb6MQkVS9IgE2AW3AHByl0Ah8SlW7RCTk5keBJiAE9AKzcV7wmlAOwJuAnUCfql4nIl8G3gv8CngXcBC4TFWXisgPgX5gC/AWYAnwbhyHwNeBB4ENqnrMlf9B4HlgOnAp8DCOM8w3gY3Ai8BVwCmc5WGfUdWfusdeCrwdR5H04Cwpm3Kpv0yKU1Wfc8srhvI8haP0xqQ4VbVHRN5MCuUpIjXAB0fXB/CQqsbSXMMNQNyt182qGnfTl+G0jw8Ar7oyfwrnvv5UVU+lKOvbqvpXIjJdVXvctA8M1wfOB8NbVPUFEZkIfAa4wC3774CPA79U1ddGlTsTWI5zTx4AvgZMA/4W+KMrY2J9/EpVn3ePzbp9uPtX9MdVvutDVZ/LdN3FuOZ8Pxcp66ICFEUfzoULyS/4i4Gdqvoed7+Lgf+F86BtVtUZbvpLqvpm9/8ncV7OFwP/qKpb3fRHVPX9IrLFVQa/Bq5KeLE8DZxS1atF5HFVvTpBvicBVPUqEXkDcBNwA46SagU+rKpXufvuVNWL3P+PqurMhHIeU9X3DZcvIj/AWXf8cSCC8/K4Guel9+XR1YS34tyC01iF4ijPvwZ+z9gU538GXgAmkkJ5Au9zr++JUfXxJ8BfcSYC/BZnBkjMlf2/qeorbn30Av8B1LjneshNuxZ4G9CJo2RG10fMveY17jW0uvVxNo4X7FIR+Sf33MP18V9wlMaLwFz3Xj2gqi+JyGbgH92yPgN8C+dF/dfAPuBl97ilwFTgqHsP55B9+xiuj4r9uAIuyXN9bAHq0lw3qjqzSNec83Ohqn+W4rpTM1YPPb/9gB1AQ4r0x4DfAPUJaTPcmzaQkHZDwv9b3b/1wGeB+4EbgUfc9MPAT4EDwMSE47YDfw78A/AT4GfAXwD3AfcCT6aQbw5wC/BL4A7gTlfmrwKfxHkQbwPeD9wN/NA97kn371Np6mPIbXxPun+H/w8Dv07Y72JgK/BWd5/EvJcS/n8S+ArOy+rKhPThOtni/v01MCEh/2mcryOG/ybkHXP/vsG93q3Ao26dP5mw385Rcjwx+h4Plw9sS1Mf23BeaD92703iL5qwXxDYjPOVt2WUHM8n/P8E8P8DP3fvU226+hglx6/dY8W9VhmVN3xfJ7nl/8xtV50J++0aVR9PjjrHEwltP137eAroy6GNHCti+xi5tjG2kaMe7SPf9bEl3XUX+Zpzfi5Spaf7VUKPYh5Od2xgVHot8KfAH1X1SEJ6DfDfgR+o6lBCej1wnapuGFXGnwPnq+rXRWRhwilCqjooIlOAy1X1EREJ4nxxzsHR3s+o6osicq2qPppG/ok4Xcp2nGGOT+C8TO7H+Ro5F3gF2KiqcREJqtOl/D7OS+VxnBfhNHf/m4DFqhoZdZ7H3P2vGq4rEZmB80JaAuxT1Xe56Teo6kb3/62qeqVbP/8NuALnJXmrOr2swzgv2KXAeara7x63Hfgh8B6cr/I6nAfnYmCZqjaPkm8Ozkv6RpweYgCnu7wJ5yv5o8AxnJ5Im1vuRFX9kttre9hN25pQH+/BeRF8yD1n16hzHgFaVLXX3a7H6WF8COdhHK6Pt+rp4Zytqnql+//1OENGzwA3qtPT63Hle6Nbdo87PPo8juL/NE5PZC7Oi+NCt14+qG7PMkG+Wvf6/4jz5XoWztDBUbf8Q269tgFX4rzAV7n18QKp20cUuBxYOsY28j5VrXfT8tk+rsCxlSa2j1PAkhT1MbqNnIXTu9uE0zZWuWUO18dZCe0jU328V90hwwz1MRP4Z1e+OamuG6f38hjO1322z0Su1/wRkp+LK3Geiy9meC62qep3yJKyVxTVjIhcgjMuOR3nBfJbHCW1JYXifCuOAhqtOGuBbwD/DuxJoTy/rqrfHrX/nwPnu8eMlAUcUtUBV3l+RVW/nUp5Aleo6n0pruetOMMVqRTndpyhoA9xpvK8QlWfEpF3A2926yKC83I+F6fXeT5njjN/EOjC+dBITH+7e+4gMDic544Tv11VtyWOW4vIVTjj0L9Tx9ZzkZu32z1uEnCxqj7r5l2O83APy7gI50URSyHj23CUpuIolmvd+uhV1d+413wV0Kqqbe4x71RnHPsynJdhrVu2qurd7j15g3uNI/YGnK/TS3BeYIl5MZwhlSdItlHEcQrdKKftW0eBL+IMqb5DRM7BCUI3AedlfX3CdW4SkStwhmVqcHrDQ6r6HRG5zr03Z9hE3Pr8PE57fxHni7sH+AVOL28isMuV7ySOAj3kPi9v5/TzMktV7xSRuTi9kZFnRkQ+h2M7upQzn5nPA12q+i8JacPPxXdwPryGGf6g/ByOnWJTwjMx220Djaq6WkSu0VE2n2GbR8IH5T739wn3Pv2Te40fcu/pKzjDb5/FGb48JCKzgMvca9mH81H4PGPAFEWZIt5G/GuzTB/m0TR5XsflO2+85bXhPHhDnDmW/AKpx5k7cXocqcaf0x2TTZ6XHLnI+OscyvuDWz8DnGlTqCG9PSpd3n/CUfJjOWY4z0uOdHnDMzLzJeOwHWXYlgBOr+xlzrRvwmkbS6q8VMcN57/Jle/FDOWNR47h/NHHJZaZykZ0Nc5Q97uAg6r6dbJlLONU9vPPj9Pjqom/J90Gk2681WssNjHvyRzz0pWZi4xbPM41uryx2F/Sjbtne0yx8/Ih42ibQi55PXkurxB5XjJ62VHymlfMc2WRl9ZeMqb3TbFfcPbLz4/0RvzeNOmPeRxT7LxCyJhu4sLrHnkDORxT7LycZUxIT5qwkWNeJM/lFSIvrYzu35STVAqRV8xzZZAj7QScbN4zI/uPZWf7+ecHzEt8SSSkn50mvdbjmGLnFULGy4DZo9JrgJs98lbmcEyx83KV8U1Azai8evclkkve5/JcXiHy0sqYor18Erg7TVvKW14xz5UqD1iY8Ktz06YA7x99rNfPbBSGYRiGJxbryTAMw/DEFIVhGIbhiSkKw8gTIvJfRORv8lTWH93574ZRckxRGIZhGJ6YojCMDIjIZBH5lYi8KCI7ReSjIvJWEXnGTfudiEx1dw+KyCYReVVEvpNQxsdE5CX3+HsypRuGn6jY9SgMI49chxOK4f8DEJEGHK/oj6oTOnoaTggMcCKBXoLjHfyKiNyH4zl9D04IhWPAZjd8yO9SpavqQ8W7NMPIjPUoDCMzLwFXi8g9InI50IwT1+p5AFXt1dOx/Z9Q1Yg6a1bswpm//lYcp68ud7/1OEHk0qUbhq8wRWEYGVDVvThf/S8Ba3ACsKVzQIom/D+E02uXNPumSzcMX2GKwjAy4Eb77FPVnwHfxYlAGnSj3SIiU93ooel4DniPiMwSJ8z9x3Bi76RLNwxfYTYKw8jMm4F7RSQODOKsOSDAfW74536cyJwpUSfU8+04wekEeFhVWwHSpRuGn7AQHoZhGIYnNvRkGIZheGKKwjAMw/DEFIVhGIbhiSkKwzAMwxNTFIZhGIYnpigMwzAMT0xRGIZhGJ6YojAMwzA8+X8rsQbp9xshjgAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sns.boxplot(x=\"school\", y=\"mathcent\", data=jsp)\n",
"plt.xticks(fontsize=8, rotation=90)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"OLS Regression Results \n",
"\n",
" Dep. Variable: mathcent R-squared: 0.082 \n",
" \n",
"\n",
" Model: OLS Adj. R-squared: 0.068 \n",
" \n",
"\n",
" Method: Least Squares F-statistic: 5.935 \n",
" \n",
"\n",
" Date: Wed, 26 Sep 2018 Prob (F-statistic): 1.46e-33 \n",
" \n",
"\n",
" Time: 10:18:07 Log-Likelihood: -11032. \n",
" \n",
"\n",
" No. Observations: 3236 AIC: 2.216e+04 \n",
" \n",
"\n",
" Df Residuals: 3187 BIC: 2.246e+04 \n",
" \n",
"\n",
" Df Model: 48 \n",
" \n",
"\n",
" Covariance Type: nonrobust \n",
" \n",
"
\n",
"\n",
"\n",
" coef std err t P>|t| [0.025 0.975] \n",
" \n",
"\n",
" C(school)[1] -3.3685 0.769 -4.383 0.000 -4.875 -1.861 \n",
" \n",
"\n",
" C(school)[2] 0.6714 1.229 0.546 0.585 -1.738 3.081 \n",
" \n",
"\n",
" C(school)[3] 2.2964 1.064 2.158 0.031 0.210 4.383 \n",
" \n",
"\n",
" C(school)[4] -2.6619 0.869 -3.064 0.002 -4.365 -0.958 \n",
" \n",
"\n",
" C(school)[5] 1.8320 0.809 2.264 0.024 0.245 3.419 \n",
" \n",
"\n",
" C(school)[6] 0.2381 0.952 0.250 0.802 -1.628 2.104 \n",
" \n",
"\n",
" C(school)[7] 1.7227 1.181 1.459 0.145 -0.592 4.037 \n",
" \n",
"\n",
" C(school)[8] 0.5909 0.790 0.748 0.455 -0.959 2.141 \n",
" \n",
"\n",
" C(school)[9] 2.2458 0.914 2.456 0.014 0.453 4.039 \n",
" \n",
"\n",
" C(school)[10] -1.2453 1.505 -0.827 0.408 -4.196 1.705 \n",
" \n",
"\n",
" C(school)[11] -0.4619 1.246 -0.371 0.711 -2.905 1.981 \n",
" \n",
"\n",
" C(school)[12] 0.2082 0.840 0.248 0.804 -1.439 1.855 \n",
" \n",
"\n",
" C(school)[13] -1.6474 0.888 -1.856 0.064 -3.388 0.093 \n",
" \n",
"\n",
" C(school)[14] -3.2898 1.124 -2.926 0.003 -5.494 -1.086 \n",
" \n",
"\n",
" C(school)[15] -0.4921 1.013 -0.486 0.627 -2.478 1.493 \n",
" \n",
"\n",
" C(school)[16] -3.7374 1.013 -3.691 0.000 -5.723 -1.752 \n",
" \n",
"\n",
" C(school)[17] 1.0881 1.648 0.660 0.509 -2.144 4.320 \n",
" \n",
"\n",
" C(school)[18] -0.8619 0.994 -0.867 0.386 -2.811 1.087 \n",
" \n",
"\n",
" C(school)[19] 1.6108 1.111 1.449 0.147 -0.568 3.790 \n",
" \n",
"\n",
" C(school)[20] -0.2095 1.138 -0.184 0.854 -2.440 2.021 \n",
" \n",
"\n",
" C(school)[21] -3.6185 0.769 -4.708 0.000 -5.125 -2.111 \n",
" \n",
"\n",
" C(school)[22] -1.7271 1.087 -1.589 0.112 -3.858 0.404 \n",
" \n",
"\n",
" C(school)[23] -0.8171 0.968 -0.844 0.399 -2.715 1.081 \n",
" \n",
"\n",
" C(school)[24] 2.6238 0.929 2.825 0.005 0.803 4.445 \n",
" \n",
"\n",
" C(school)[25] 1.4232 1.075 1.323 0.186 -0.685 3.532 \n",
" \n",
"\n",
" C(school)[26] 0.4919 0.835 0.589 0.556 -1.145 2.129 \n",
" \n",
"\n",
" C(school)[27] -2.3058 0.863 -2.672 0.008 -3.998 -0.614 \n",
" \n",
"\n",
" C(school)[28] -5.8286 1.064 -5.478 0.000 -7.915 -3.742 \n",
" \n",
"\n",
" C(school)[29] 0.2090 0.936 0.223 0.823 -1.627 2.045 \n",
" \n",
"\n",
" C(school)[30] 0.2392 0.773 0.309 0.757 -1.276 1.754 \n",
" \n",
"\n",
" C(school)[31] 3.8241 0.713 5.366 0.000 2.427 5.221 \n",
" \n",
"\n",
" C(school)[32] -0.0268 0.857 -0.031 0.975 -1.707 1.654 \n",
" \n",
"\n",
" C(school)[33] 1.0633 0.644 1.651 0.099 -0.200 2.326 \n",
" \n",
"\n",
" C(school)[34] 1.9359 0.769 2.519 0.012 0.429 3.443 \n",
" \n",
"\n",
" C(school)[35] 1.9041 1.013 1.880 0.060 -0.081 3.890 \n",
" \n",
"\n",
" C(school)[36] 2.3605 0.781 3.021 0.003 0.828 3.893 \n",
" \n",
"\n",
" C(school)[37] 1.7618 0.960 1.836 0.067 -0.120 3.644 \n",
" \n",
"\n",
" C(school)[38] -1.6106 1.181 -1.364 0.173 -3.925 0.704 \n",
" \n",
"\n",
" C(school)[39] -1.4842 1.099 -1.350 0.177 -3.639 0.671 \n",
" \n",
"\n",
" C(school)[40] -4.9855 1.264 -3.943 0.000 -7.464 -2.506 \n",
" \n",
"\n",
" C(school)[41] 0.8809 1.246 0.707 0.480 -1.562 3.324 \n",
" \n",
"\n",
" C(school)[42] -0.3261 0.644 -0.506 0.613 -1.589 0.937 \n",
" \n",
"\n",
" C(school)[44] 0.5881 1.843 0.319 0.750 -3.026 4.202 \n",
" \n",
"\n",
" C(school)[45] -4.6392 1.111 -4.174 0.000 -6.818 -2.460 \n",
" \n",
"\n",
" C(school)[46] 3.0636 1.032 2.968 0.003 1.039 5.088 \n",
" \n",
"\n",
" C(school)[47] 2.3969 0.730 3.284 0.001 0.966 3.828 \n",
" \n",
"\n",
" C(school)[48] 0.2119 0.514 0.412 0.680 -0.795 1.219 \n",
" \n",
"\n",
" C(school)[49] 2.7964 0.869 3.219 0.001 1.093 4.500 \n",
" \n",
"\n",
" C(school)[50] -2.6520 0.734 -3.615 0.000 -4.090 -1.214 \n",
" \n",
"
\n",
"\n",
"\n",
" Omnibus: 136.092 Durbin-Watson: 1.213 \n",
" \n",
"\n",
" Prob(Omnibus): 0.000 Jarque-Bera (JB): 151.071 \n",
" \n",
"\n",
" Skew: -0.521 Prob(JB): 1.57e-33 \n",
" \n",
"\n",
" Kurtosis: 2.814 Cond. No. 3.59 \n",
" \n",
"
Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: mathcent R-squared: 0.082\n",
"Model: OLS Adj. R-squared: 0.068\n",
"Method: Least Squares F-statistic: 5.935\n",
"Date: Wed, 26 Sep 2018 Prob (F-statistic): 1.46e-33\n",
"Time: 10:18:07 Log-Likelihood: -11032.\n",
"No. Observations: 3236 AIC: 2.216e+04\n",
"Df Residuals: 3187 BIC: 2.246e+04\n",
"Df Model: 48 \n",
"Covariance Type: nonrobust \n",
"=================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"---------------------------------------------------------------------------------\n",
"C(school)[1] -3.3685 0.769 -4.383 0.000 -4.875 -1.861\n",
"C(school)[2] 0.6714 1.229 0.546 0.585 -1.738 3.081\n",
"C(school)[3] 2.2964 1.064 2.158 0.031 0.210 4.383\n",
"C(school)[4] -2.6619 0.869 -3.064 0.002 -4.365 -0.958\n",
"C(school)[5] 1.8320 0.809 2.264 0.024 0.245 3.419\n",
"C(school)[6] 0.2381 0.952 0.250 0.802 -1.628 2.104\n",
"C(school)[7] 1.7227 1.181 1.459 0.145 -0.592 4.037\n",
"C(school)[8] 0.5909 0.790 0.748 0.455 -0.959 2.141\n",
"C(school)[9] 2.2458 0.914 2.456 0.014 0.453 4.039\n",
"C(school)[10] -1.2453 1.505 -0.827 0.408 -4.196 1.705\n",
"C(school)[11] -0.4619 1.246 -0.371 0.711 -2.905 1.981\n",
"C(school)[12] 0.2082 0.840 0.248 0.804 -1.439 1.855\n",
"C(school)[13] -1.6474 0.888 -1.856 0.064 -3.388 0.093\n",
"C(school)[14] -3.2898 1.124 -2.926 0.003 -5.494 -1.086\n",
"C(school)[15] -0.4921 1.013 -0.486 0.627 -2.478 1.493\n",
"C(school)[16] -3.7374 1.013 -3.691 0.000 -5.723 -1.752\n",
"C(school)[17] 1.0881 1.648 0.660 0.509 -2.144 4.320\n",
"C(school)[18] -0.8619 0.994 -0.867 0.386 -2.811 1.087\n",
"C(school)[19] 1.6108 1.111 1.449 0.147 -0.568 3.790\n",
"C(school)[20] -0.2095 1.138 -0.184 0.854 -2.440 2.021\n",
"C(school)[21] -3.6185 0.769 -4.708 0.000 -5.125 -2.111\n",
"C(school)[22] -1.7271 1.087 -1.589 0.112 -3.858 0.404\n",
"C(school)[23] -0.8171 0.968 -0.844 0.399 -2.715 1.081\n",
"C(school)[24] 2.6238 0.929 2.825 0.005 0.803 4.445\n",
"C(school)[25] 1.4232 1.075 1.323 0.186 -0.685 3.532\n",
"C(school)[26] 0.4919 0.835 0.589 0.556 -1.145 2.129\n",
"C(school)[27] -2.3058 0.863 -2.672 0.008 -3.998 -0.614\n",
"C(school)[28] -5.8286 1.064 -5.478 0.000 -7.915 -3.742\n",
"C(school)[29] 0.2090 0.936 0.223 0.823 -1.627 2.045\n",
"C(school)[30] 0.2392 0.773 0.309 0.757 -1.276 1.754\n",
"C(school)[31] 3.8241 0.713 5.366 0.000 2.427 5.221\n",
"C(school)[32] -0.0268 0.857 -0.031 0.975 -1.707 1.654\n",
"C(school)[33] 1.0633 0.644 1.651 0.099 -0.200 2.326\n",
"C(school)[34] 1.9359 0.769 2.519 0.012 0.429 3.443\n",
"C(school)[35] 1.9041 1.013 1.880 0.060 -0.081 3.890\n",
"C(school)[36] 2.3605 0.781 3.021 0.003 0.828 3.893\n",
"C(school)[37] 1.7618 0.960 1.836 0.067 -0.120 3.644\n",
"C(school)[38] -1.6106 1.181 -1.364 0.173 -3.925 0.704\n",
"C(school)[39] -1.4842 1.099 -1.350 0.177 -3.639 0.671\n",
"C(school)[40] -4.9855 1.264 -3.943 0.000 -7.464 -2.506\n",
"C(school)[41] 0.8809 1.246 0.707 0.480 -1.562 3.324\n",
"C(school)[42] -0.3261 0.644 -0.506 0.613 -1.589 0.937\n",
"C(school)[44] 0.5881 1.843 0.319 0.750 -3.026 4.202\n",
"C(school)[45] -4.6392 1.111 -4.174 0.000 -6.818 -2.460\n",
"C(school)[46] 3.0636 1.032 2.968 0.003 1.039 5.088\n",
"C(school)[47] 2.3969 0.730 3.284 0.001 0.966 3.828\n",
"C(school)[48] 0.2119 0.514 0.412 0.680 -0.795 1.219\n",
"C(school)[49] 2.7964 0.869 3.219 0.001 1.093 4.500\n",
"C(school)[50] -2.6520 0.734 -3.615 0.000 -4.090 -1.214\n",
"==============================================================================\n",
"Omnibus: 136.092 Durbin-Watson: 1.213\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 151.071\n",
"Skew: -0.521 Prob(JB): 1.57e-33\n",
"Kurtosis: 2.814 Cond. No. 3.59\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod = smf.ols(\"mathcent ~ C(school) - 1\", jsp).fit()\n",
"lmod.summary()"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" df \n",
" sum_sq \n",
" mean_sq \n",
" F \n",
" PR(>F) \n",
" \n",
" \n",
" \n",
" \n",
" C(school) \n",
" 48.0 \n",
" 15483.816175 \n",
" 322.579504 \n",
" 5.935264 \n",
" 1.458040e-33 \n",
" \n",
" \n",
" Residual \n",
" 3187.0 \n",
" 173212.333393 \n",
" 54.349650 \n",
" NaN \n",
" NaN \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" df sum_sq mean_sq F PR(>F)\n",
"C(school) 48.0 15483.816175 322.579504 5.935264 1.458040e-33\n",
"Residual 3187.0 173212.333393 54.349650 NaN NaN"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sm.stats.anova_lm(smf.ols(\"mathcent ~ C(school)\", jsp).fit())"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"C(school)[1] 1.210916e-05\n",
"C(school)[2] 5.848061e-01\n",
"C(school)[3] 3.099485e-02\n",
"C(school)[4] 2.203529e-03\n",
"C(school)[5] 2.364073e-02\n",
"C(school)[6] 8.024944e-01\n",
"C(school)[7] 1.445850e-01\n",
"C(school)[8] 4.547159e-01\n",
"C(school)[9] 1.410363e-02\n",
"C(school)[10] 4.080158e-01\n",
"C(school)[11] 7.108935e-01\n",
"C(school)[12] 8.042920e-01\n",
"C(school)[13] 6.351130e-02\n",
"C(school)[14] 3.455110e-03\n",
"C(school)[15] 6.270231e-01\n",
"C(school)[16] 2.273774e-04\n",
"C(school)[17] 5.092721e-01\n",
"C(school)[18] 3.859690e-01\n",
"C(school)[19] 1.473420e-01\n",
"C(school)[20] 8.538620e-01\n",
"C(school)[21] 2.610820e-06\n",
"C(school)[22] 1.121715e-01\n",
"C(school)[23] 3.986798e-01\n",
"C(school)[24] 4.759221e-03\n",
"C(school)[25] 1.857784e-01\n",
"C(school)[26] 5.556972e-01\n",
"C(school)[27] 7.572699e-03\n",
"C(school)[28] 4.647573e-08\n",
"C(school)[29] 8.233416e-01\n",
"C(school)[30] 7.569780e-01\n",
"C(school)[31] 8.645474e-08\n",
"C(school)[32] 9.750611e-01\n",
"C(school)[33] 9.889090e-02\n",
"C(school)[34] 1.182717e-02\n",
"C(school)[35] 6.015578e-02\n",
"C(school)[36] 2.541872e-03\n",
"C(school)[37] 6.650601e-02\n",
"C(school)[38] 1.725456e-01\n",
"C(school)[39] 1.769603e-01\n",
"C(school)[40] 8.213962e-05\n",
"C(school)[41] 4.796616e-01\n",
"C(school)[42] 6.127513e-01\n",
"C(school)[44] 7.496910e-01\n",
"C(school)[45] 3.070621e-05\n",
"C(school)[46] 3.023142e-03\n",
"C(school)[47] 1.035959e-03\n",
"C(school)[48] 6.800311e-01\n",
"C(school)[49] 1.301069e-03\n",
"C(school)[50] 3.046839e-04\n",
"dtype: float64"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod.pvalues"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"from statsmodels.sandbox.stats.multicomp import multipletests\n",
"reject, padj, _, _ = multipletests(lmod.pvalues, method=\"bonferroni\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These schools are rejected by Bonferroni."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"C(school)[1] -3.368450\n",
"C(school)[16] -3.737400\n",
"C(school)[21] -3.618450\n",
"C(school)[28] -5.828595\n",
"C(school)[31] 3.824053\n",
"C(school)[40] -4.985458\n",
"C(school)[45] -4.639201\n",
"C(school)[50] -2.652027\n",
"dtype: float64"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lmod.params[reject]"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Index(['C(school)[28]', 'C(school)[31]', 'C(school)[21]', 'C(school)[1]',\n",
" 'C(school)[45]', 'C(school)[40]', 'C(school)[16]', 'C(school)[50]',\n",
" 'C(school)[47]', 'C(school)[49]', 'C(school)[4]', 'C(school)[36]',\n",
" 'C(school)[46]', 'C(school)[14]', 'C(school)[24]', 'C(school)[27]',\n",
" 'C(school)[34]', 'C(school)[9]'],\n",
" dtype='object')"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"selsch = np.argsort(lmod.pvalues)[np.sort(lmod.pvalues) < np.arange(1,50)*0.05/49]\n",
"lmod.params.index[selsch]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Schools identified using false discovery rate."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"C(school)[1] -3.368450\n",
"C(school)[4] -2.661928\n",
"C(school)[9] 2.245764\n",
"C(school)[14] -3.289835\n",
"C(school)[16] -3.737400\n",
"C(school)[21] -3.618450\n",
"C(school)[24] 2.623786\n",
"C(school)[27] -2.305764\n",
"C(school)[28] -5.828595\n",
"C(school)[31] 3.824053\n",
"C(school)[34] 1.935898\n",
"C(school)[36] 2.360544\n",
"C(school)[40] -4.985458\n",
"C(school)[45] -4.639201\n",
"C(school)[46] 3.063562\n",
"C(school)[47] 2.396895\n",
"C(school)[49] 2.796405\n",
"C(school)[50] -2.652027\n",
"dtype: float64"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"reject, padj, _, _ = multipletests(lmod.pvalues, method=\"fdr_bh\")\n",
"lmod.params[reject]"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"application/json": {
"Software versions": [
{
"module": "Python",
"version": "3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]"
},
{
"module": "IPython",
"version": "6.5.0"
},
{
"module": "OS",
"version": "Darwin 17.7.0 x86_64 i386 64bit"
},
{
"module": "pandas",
"version": "0.23.4"
},
{
"module": "numpy",
"version": "1.15.1"
},
{
"module": "matplotlib",
"version": "2.2.3"
},
{
"module": "seaborn",
"version": "0.9.0"
},
{
"module": "scipy",
"version": "1.1.0"
},
{
"module": "patsy",
"version": "0.5.0"
},
{
"module": "statsmodels",
"version": "0.9.0"
}
]
},
"text/html": [
"Software Version Python 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)] IPython 6.5.0 OS Darwin 17.7.0 x86_64 i386 64bit pandas 0.23.4 numpy 1.15.1 matplotlib 2.2.3 seaborn 0.9.0 scipy 1.1.0 patsy 0.5.0 statsmodels 0.9.0 Wed Sep 26 10:18:07 2018 BST
"
],
"text/latex": [
"\\begin{tabular}{|l|l|}\\hline\n",
"{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
"Python & 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE\\_401/final)] \\\\ \\hline\n",
"IPython & 6.5.0 \\\\ \\hline\n",
"OS & Darwin 17.7.0 x86\\_64 i386 64bit \\\\ \\hline\n",
"pandas & 0.23.4 \\\\ \\hline\n",
"numpy & 1.15.1 \\\\ \\hline\n",
"matplotlib & 2.2.3 \\\\ \\hline\n",
"seaborn & 0.9.0 \\\\ \\hline\n",
"scipy & 1.1.0 \\\\ \\hline\n",
"patsy & 0.5.0 \\\\ \\hline\n",
"statsmodels & 0.9.0 \\\\ \\hline\n",
"\\hline \\multicolumn{2}{|l|}{Wed Sep 26 10:18:07 2018 BST} \\\\ \\hline\n",
"\\end{tabular}\n"
],
"text/plain": [
"Software versions\n",
"Python 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]\n",
"IPython 6.5.0\n",
"OS Darwin 17.7.0 x86_64 i386 64bit\n",
"pandas 0.23.4\n",
"numpy 1.15.1\n",
"matplotlib 2.2.3\n",
"seaborn 0.9.0\n",
"scipy 1.1.0\n",
"patsy 0.5.0\n",
"statsmodels 0.9.0\n",
"Wed Sep 26 10:18:07 2018 BST"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%load_ext version_information\n",
"%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels"
]
}
],
"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.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}