{ "metadata": { "name": "regression_with_factor_vars" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Regression with factor variables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Key ideas:\n", "\n", "* Outcome is still quantitative\n", "* Covariate(s) are factor variables (`Categorical` in pandas)\n", "* Fitting lines = fitting means\n", "* Want to evaluate contribution of all factor levels at once" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Example: movie ratings\n", "\n", "import pandas as pd\n", "\n", "movies = pd.read_csv('http://www.rossmanchance.com/iscam2/data/movies03RT.txt', sep='\\t', \n", " names=['X', 'score', 'rating', 'genre', 'box_office', 'running_time'], skiprows=1)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "movies.head()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Xscoreratinggenrebox_officerunning_time
0 2 Fast 2 Furious 48.9 PG-13 action/adventure 127.146 107
1 28 Days Later 78.2 R horror 45.065 113
2 A Guy Thing 39.5 PG-13 rom comedy 15.545 101
3 A Man Apart 42.9 R action/adventure 26.248 110
4 A Mighty Wind 79.9 PG-13 comedy 17.781 91
\n", "
" ], "output_type": "pyout", "prompt_number": 2, "text": [ " X score rating genre box_office running_time\n", "0 2 Fast 2 Furious 48.9 PG-13 action/adventure 127.146 107\n", "1 28 Days Later 78.2 R horror 45.065 113\n", "2 A Guy Thing 39.5 PG-13 rom comedy 15.545 101\n", "3 A Man Apart 42.9 R action/adventure 26.248 110\n", "4 A Mighty Wind 79.9 PG-13 comedy 17.781 91" ] } ], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "def jitter(series, factor):\n", " z = float(series.max()) - float(series.min())\n", " a = float(factor) * z / 50.\n", " return map(lambda x: x + np.random.uniform(-a, a), series)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rotten tomatoes score vs rating (average scores per rating are marked with red squares):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "r = pd.Categorical.from_array(movies['rating'])\n", "\n", "plot(jitter(r.labels, 2.), movies['score'], 'ob')\n", "xticks(np.unique(r.labels), r.levels)\n", "xlabel('rating')\n", "ylabel('score')\n", "\n", "meanRatings = movies.groupby('rating')['score'].mean()\n", "\n", "plot(range(4), meanRatings, 'sr', markersize=10);" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEMCAYAAAAoB2Y1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0U1XaP/BvoIWWQi9cmlJuxXKTgtMCIvobmyhvGt9V\n5S3jayuKE0RAFIFxZgkObSWKvARdXlqGWTrC0ryvOshao0ynVSiobXEUuYiUFhSttNyaKG25t2DK\n+f0RkjbtyT0nSZPvZ62uRU/OOdnNIc/ZZ+9n7y0TBEEAERGFhV6BLgAREfkPgz4RURhh0CciCiMM\n+kREYYRBn4gojDDoExGFEZ8F/fnz50Mul2Py5MnWbc3NzVCpVBg3bhyysrJw7tw562vr1q3D2LFj\nMWHCBJSXl/uqGERE5IDPgv6jjz6K7du322zT6XRQqVQ4duwYZs6cCZ1OBwA4cuQIPvjgAxw5cgTb\nt2/Hk08+ievXr/uqKEREZIfPgv6dd96JhIQEm20lJSXQaDQAAI1Gg23btgEA/vnPf2LOnDmIjIxE\nSkoKxowZg7179/qqKEREZEeElCc3Go2Qy+UAALlcDqPRCAA4c+YMZsyYYd1v+PDhOH36dLfjZTKZ\nlMUjIgpJjiZa8FtHrkwmcxjE7b0mCILDn9WrVzvdhz+B++H1Cf4fXqPg/nH3+jgjadCXy+UwGAwA\ngMbGRiQmJgIAhg0bhpMnT1r3O3XqFIYNGyZlUYiICBIH/VmzZkGv1wMA9Ho9cnJyrNu3bNmCa9eu\n4fjx4/jhhx8wffp0KYtCRETwYZv+nDlzUFlZibNnz2LEiBF44YUX8OyzzyI3NxebN29GSkoKtm7d\nCgCYOHEicnNzMXHiREREROCvf/2rx+33SqXSV38CSYDXJ/jxGgU3X18fmeBKI1CAyGQyl9qoiIjI\nzFnc5IhcIqIwwqBPRBRGGPSJiMIIgz4RURhh0CciCiMM+kREYUTSuXeIiMg7ZWVVKC4ux9WrEejb\n14Rly7KQnZ3p8fkY9ImIglRZWRWWL9+Burq11m11dfkA4HHgZ/MOEVGQKi4utwn4AFBXtxYbNuz0\n+JwM+kREQerqVfHGmLa23h6fk0GfiChI9e1rEt0eFdXu8TkZ9ImIgtSyZVlITc232ZaaugpLl6o8\nPicnXCMKQb7O+KDAKSurwoYNO9HW1htRUe1YulTl8Fo6i5sM+kQhRizjIzU1H0VFagb+MMBZNonC\njBQZHxQ6gj7oK5VaqNUFKCurCnRRiHoEKTI+KHQE/eCsykotAO8HJBCFCykyPih0BH1N34KPp0Su\nkSLjg0JH0Nf0O+PjKZFzlqfhDRsKO2V83MOnZALgp6BfVFSETZs2QRAELFy4EMuXL0dzczPy8vLQ\n0NBgXTQ9Pj7e4Xn4eErkmuzsTAb5EOSLVFzJm3dqamqwadMm7Nu3D4cOHUJpaSnq6uqg0+mgUqlw\n7NgxzJw5EzqdzuF5+HhKROHMkopbXv4iKiu1KC9/EcuX73A7yUXyoP/dd9/htttuQ1RUFHr37g2F\nQoF//OMfKCkpgUajAQBoNBps27ZN9HiFQgu1uhBFRXw8JaLw5atUXMmbdyZNmoT8/Hw0NzcjKioK\nH3/8MaZNmwaj0Qi5XA4AkMvlMBqNoscrlQDQG/v2fYaYmOtQmjcQEYUVe6m4BsNJaLVal88jedCf\nMGECVq5ciaysLMTExCA9PR29e9t2yMpkMshkMtHj3fljiIhClb1U3KSkETZx8vnnn3d4Hr+kbM6f\nPx/79+9HZWUlEhISMG7cOMjlchgMBgBAY2MjEhMT/VEUIqIeyVepuH7J3vn555+RmJiIEydO4MMP\nP8SePXtw/Phx6PV6rFy5Enq9Hjk5Of4oClFY4IRrocdXqbh+mXAtMzMTTU1NiIyMxGuvvYa77roL\nzc3NyM3NxYkTJ+ymbHLCNSL3ccK18MZZNonCjFpdgPLyF0W2F2L79jUBKBH5E2fZJAoznHCNHOlR\n0zAQkXNSTLjGPoLQwaBPFGKWLctCXV1+lzb9VVi69B6PzifWR8BZb3sutukThSB3l9hzhH0EPYuz\nuMmaPlEI8uWEa+wjCC3syCUih7goS2hh0Ccih7goS2hhmz4ROeXLPgIS56sMKQ7OIiIKcr4cRc3B\nWUREQc5Xc+W7gkGfiCjA/JkhxaBPRBRg/syQYtAnIgowf2ZIsSOXiCgI+CpDitk7RERhhNMwEJFT\nnEUzfDDoE4U5zqIZXtiRSxTm/JkjToHnl6C/bt06pKWlYfLkyXjooYdw9epVNDc3Q6VSYdy4ccjK\nysK5c+f8URQi6oKzaIYXyYN+fX093nrrLXzzzTc4fPgw2tvbsWXLFuh0OqhUKhw7dgwzZ86ETqeT\nuihEJMJejnhNzVGUlVX5uTQkNcmDfmxsLCIjI3HlyhWYTCZcuXIFycnJKCkpgUajAQBoNBps27ZN\n6qIQkQixHHFgFZqalmD58h0M/CFG8o7cgQMH4k9/+hNGjhyJ6OhoqNVqqFQqGI1GyOVyAIBcLofR\naBQ9XqvVWv+tVCqhVCqlLjJRWLF01mo0eWhquhlAO4B7AGSiri4TGzYUskM3iFVUVKCiosLl/SXP\n06+rq8N9992H3bt3Iy4uDg888ADuv/9+LF26FC0tLdb9Bg4ciObmZtvCMU+fyG+USi0qK7XdtisU\nWlRUdN9OwSngefr79+/HHXfcgUGDBgEAfve73+Grr75CUlISDAYDkpKS0NjYiMTERKmLQkQOcIUs\n6QXDeAjJg/6ECROwZs0atLa2IioqCrt27cL06dMRExMDvV6PlStXQq/XIycnR+qiEIUFTwPLsmVZ\nqKvL7zKn+yosXXqPlMUNG96Oh/DZDUPwg/Xr1wsTJ04UJk2aJPz+978Xrl27JjQ1NQkzZ84Uxo4d\nK6hUKqGlpaXbcX4qHlHIKC2tFFJTVwmAYP1JTV0llJZWuny8Wl0gKBSrBbW6wOXjyLmsrHyb62L5\nUasLnB7rznV1FjeDOqoy6BO5x5vAQtJSKFaLXhuFYrXTY925rs7iJkfkEoUQDrQKXt70mfjyujLo\nE4UQdsYGL2/mzPfldeWEa0QhhJ2xwcvS6bphQ2GnOfPv8XsnO+fTJwoxvlqMgwKra7bO7bcnY8+e\nRqfXlYuoEIWxYMgLJ/eJpXempuajqEjt9PoFfHAWEQUG58nvuexPd+39lBjsyCWfKyurglpdAKVS\nC7W6gBN2BQjnye+5pMzCYk2ffIq1y+DB9M2eS8osLNb0yadYuwweTN/subxJ73SGNX3yKdYugwfT\nN3sub9I7nWHQJ59i7TJ4SBk4SHrZ2ZmSXCumbJJPiaearUJREYMNkT8wT5/8joODeibm9IcGBn0i\ncsqbwUAUXBj0KeBYgwx+anUBystfFNleiO3b1wSgROQpjsilgGLefs/ArKvg5etKE4M+SUrK4eTk\nO8y68g93A7gUlSYGfZIUa5A9A3P6pedJAJei0sSgT5KyV4O8ePEXqNUFbOcPEszpl54nAdyVSlPX\npwdnJA/633//PR588EHr7z/99BPWrFmDuXPnIi8vDw0NDUhJScHWrVsRHx8vdXHIz8RqkElJT+PM\nmTZ8801HxyHb+QNPqsFAZObJU685iFcBKIc5XJsAZFmb3cSeHoC13c7TmeRz74wfPx4HDx7EwYMH\nceDAAfTr1w+zZ8+GTqeDSqXCsWPHMHPmTOh0OqmLQgGQnZ2JoiI11OpCKBRaqNWFGDr0IgyGzTb7\ncX4eCnWe9JvcfnsyIiLeB/AiAC2AFxER8T5mzBgKQPzpwRm/Nu/s2rULY8aMwYgRI1BSUoLKykoA\ngEajgVKpZOAPUV1rkEqlVnQ/tvNTKOt46lXDUnOPjj6KGTMUdo/56qszMJnesNlmMr2BPXsKAdh/\nenDEr0F/y5YtmDNnDgDAaDRCLpcDAORyOYxGo+gxWq3W+m+lUgmlUil1MUlizBTpeTjWwnvZ2ZnY\nt68GL730PlpbzYG8tRV499183Hprlejn6axJyPxdqrjx4yLBT65evSoMHjxY+PnnnwVBEIT4+Hib\n1xMSErod48fikR+VllYKqamrBECw/qSm/lkoLa0MdNFIhPj1WsXr5YGsrHybz9Hyo1YXeLS/2LVx\nFjf9VtP/5JNPMHXqVAwZMgSAuXZvMBiQlJSExsZGJCYmOj2Hdt48oL7e+ZulpED7zjtelZekw0yR\nnoVjLXzH3c5cZ6m0Yt+lHTscl8FvQf/vf/+7tWkHAGbNmgW9Xo+VK1dCr9cjJyfH+Unq66G90Q/g\niNaLcpJ/MFPEv7xpnuFYC98oK6tCTc1R0dfsNW26UkHq+l2SybpPp9GZX4L+5cuXsWvXLrz11lvW\nbc8++yxyc3OxefNma8omhR62BQeet6M62QfjPcs1aGpaAiAfndMqnQ2C83UFyS9BPyYmBmfPnrXZ\nNnDgQOzatcsfb08Bwnl3gkNh4Qeoq9tos82d5hmxJobo6McxY8ZvfF7WUNW9iawQQG8MHvwdioqe\n9Ov3gWvkkmS4Xm7glZVV4ejRS6Kvudo8k52diblzhyE6Og/mxtNCtLY+jHffPY2ysiqflTWU2TaR\nZQJYA0CLtLQJTufeUasLoFRqoVYX+OTz5jQMJBm2BQdecXE52tpGiL7WeVSnsya4r746g9bWD2y2\n1dVlsjPXRZ40kUn1pMygT5JhW3DgmW+8d6NrO3JU1GIsXfqQy4GFN3DveDKhnVRZUwz6JBnO3Bh4\n5huvJUCY25GBdkyc2I7s7Eyo1QV2AwsA6xOAOeukqtO5zHgDd40nacpS3WgZ9EkyzMcPPNsbr/lz\nT01dhRdeeARlZVXYt++k6HFnzlzq9gQQEbEYJhNszsMbuOvczcKR6kmZQZ+84qw9mPn4gWXvxgsA\ny5fvQEuLeHt/fX09Ll78yGabyfQGBg9+EGlpn/EG7gdSPSn3rKCfkuLawKuUFGnLQQCYktlTiN14\nO5p1qtC1vR9YhevXfxU9V1raBFRUaCUqKXVm74a9b18NNJqNMJmiERHRiqeeUkCrfdLl8/aooM+p\nFYILh+f3XB3txV3b+78H8AT69DmNy5e7H8c2fP/qesPWav+KtWurYTJ1ZFKtXbsYwF9dDvzM0yeP\nMaOj57JtL+7IGwfGAMjE6NH9kZqab3OMuWlB5bcyUnd/+Uul6FTLGze6nr8f9DV9pVLL4ftBiimZ\nPY+lD+b06V8QHb3YOsWv2SoA99zo6M0DwE74YGMyRYtu//XXKJfPEfRBv7JSC4BtxcGIKZnBqXPn\n+oULpwD0RWzsEFy48AsaG9s6rVpWhejoPCQmxuHy5QuQy+MwfPhOm+DO75t0PJmXKiKiVXR7ZGSb\ny+/rMOhfuXIFJ0+exPjx410+oVTYVhx8mJIZfDo619UAPgDwK4AkmAdoZcLcaWvJt89Ea2smJkwo\nxPbtfwtYmcORp0kQTz2lwNq1i22aeCIiHseSJW585+xNtP/Pf/5TGDdunDBq1ChBEAThm2++Ee67\n7z4Plw7wDACbxQEUitV+fX+insa86EalAHRdWGPVje2CABS49b0qLa0UsrLyBYVitZCVlc/FU3zA\n3cVUOlu9eqMweHCeEBenEQYPzhNWr95o87qDsC4IgoNFVLRaLb7++mvcddddAICMjAz89NNPrt9N\nJMC2YiLHzJ3r5bBNwcSN3wthruHbdrQ7+l5ptX/FSy9V27T9s6nVM52bc6qrf4TYCGdLEoSjph+t\n9km3UjS7shv0IyMjER8fb7OtV6/AJfuwrZjIOXPnur2vtSXYdwR5R9+rsrIqvPRSpchEa2xqdVf3\n5pwqABsBfAQgBkAWgExERbVLPv7FbtBPS0vDe++9B5PJhB9++AHFxcW44447vH5DdykUWrYVE7lo\n2bIs7N69Ea2i/X3tAJ5GTMwPmDbN+fequLgcra03i77GtFz32I5pqQKwA+Y+F4t8JCW9g6VL50k+\n/sVu1f0vf/kLamtr0bdvX8yZMwexsbF4/fXXvX5Dd1VUaLF9+xoGfCIXZGdnYsUKBaKjF3d5ZTEA\nA4CxiIqKAQCYm3/tMzcVMS3XF2zHtIg3vyUnRyM7O1Py8S+iZzeZTMjOzsbnn3+O//mf//HJGxGR\nf2i1T+LWW6vw3HNLcOTIJbS1jQTwEAAgIuJ9NDV9AMtS046aDcxNRVnoOk1DdPTjWLr0YWn/iBBj\nO6ZFPKgPGDBEZN8OvrrRir57REQEevXqhXPnznVr1/c3tbqAA7OIXNC1869jgNVOtLV9hpqao2hq\n+gDm5oVyABGoq5Phuef+T/T7ZR6HYUn/NE/TEB19FCtWKPh9dJPtmBbHQV3q8S922/RjYmIwefJk\nqFQqxMSYHwdlMhmKi4vdfpNz585hwYIFqK2thUwmw9tvv42xY8ciLy8PDQ0N1oXRxW4w5eUvMluA\nyAl7nX9FRWps374GgHl0e2WlpT25Y78jR55AWVlVt+9XxziMnZ3GYSzh99ADnce0nDr1C376yXY0\ndOegbtn3ueeW4PjxSwCuITY2xmdlkQl2GvbeuTG5mUwmA2Bu/5PJZNBoNG6/iUajgUKhwPz582Ey\nmXD58mWsXbsWgwcPxooVK7B+/Xq0tLRAp9PZFk4mA2AunlpdaP3PS0S21OoClJe/KLK943tj3gcA\nHO9H0isrq+pyM1XZ3EzFbuKpqeabuLObrkwmc9xf4yiJv62tTaiurhaqq6uFa9euOR00IObcuXPC\n6NGju20fP368YDAYBEEQhMbGRmH8+PHd9kGnwVkcmBU6ONjH9xSK1aKDfTp/b0pLK4U+feaK7jd5\n8h8CV3jqxpvBW07Cuv3BWRUVFdBoNBg1ahQA4MSJE9Dr9VAoFK7dym44fvw4hgwZgkcffRSHDh3C\n1KlT8frrr8NoNEIulwMA5HI5jEajnTNoAQCnTn2OigollEqlW+9PwYVz8EvDlc6/7OxM9O37Kq5d\n675fY2OjVEUjD7iTwVNRUYGKigrXT27vbpCRkSF899131t+///57ISMjw417ldm+ffuEiIgIYe/e\nvYIgCMLy5cuFgoICIT4+3ma/hISEbsfiRk0/NfXPrA2GCG9qMGRfaWmlkJpqO/WC2PcmLW2RyBQN\nfxbS0hYFqOQkxt73ZNCgXKdPyA7CuiAIDmr6JpPJZqK1cePGwWQSr004Mnz4cAwfPhy33norAOC/\n//u/sW7dOiQlJcFgMCApKQmNjY1ITEwUPV6tLuTArBDCOfil4erkd8OGDUFtbRY6L5IO3IPhw3f6\nu8jkgFgGT0TE42hqWoLKSvM19fQJ2W7Qnzp1KhYsWIC5c+dCEAS89957mDZtmtuFT0pKwogRI3Ds\n2DGMGzcOu3btQlpaGtLS0qDX67Fy5Uro9Xrk5OSIHs/OpZ5LbP4Qe80Q33xzWDSDhFznynrEHWmY\nnA47mHW9iZvTbZeg81w9no7StZu909bWho0bN+Lf//43AODOO+/Ek08+ib59+7r9Bxw6dAgLFizA\ntWvXkJqairfffhvt7e3Izc3FiRMn7KZsOu2FpqBlL/tg7txhePPNH2EwvNppb/PiHUlJ27BpUw4D\nv8QcZY54Msc7Sc+cbqvttl2h0HZbs9hZ3LQb9C9fvoyoqCj07m1+7G5vb8fVq1fRr18/z0vuJgb9\nnkEsUBQXl9tNIfz550YcPHgRwM0wNy+oYKnBMHUwcLxJEyTplJVVQaPZiKamm2Ee2GWenA0Q/744\ni5t2m3fuvvtufPrpp+jfvz8A84IqarUaX375pdd/BIUOe9k4UVFnRfdva+uN2NjhN37Tir5OgSE+\n0ZcaGs1GTJr0GWv+AWD5fplHUluY2/JTU7d71CxnN+hfvXrVGvABYMCAAbhy5Yrbb0Chzd6MgIMG\n5YnuHxXVbqcWYp4aoLr6JKfekJB23jygvl70tehv66HAFwCAeqSgAfMB7HB5rh7yXtfrU119HMNb\nRmM4lDb7Gfq+j1eK9B5dB7tBv1+/fjhw4ACmTp0KANi/fz+io8UX5aXwZS8bJzk5GfHx9ucPqa7W\nw2CwTOTVMTVASwtQXs7gIpn6emgtEVxUAwBACaBBZDZIzqUvMdHrc6LbbjlRozy+BnaDflFREXJz\nczF06FAAgMFgwJYtWzx6EwpdnbNxRmEeUlAPABh4qh7DhifgdMJ7uH4d6NULGNYvAfte/hJIScGm\nTfPx3HP/h+PH5+DixSswmf4EoADm/5Im1NWpsWHDTgYXL4gtkD7yp3o3zsD02mDlzXpWdoP+8ePH\ncfDgQTQ0NODDDz/E3r17A7pyFgWnzvnEKahHBW7UUloAtDTY7nzjdy1s0wsnTXoctbW2k4AB+Th1\nSrxfgJyzt0B6LBJgqc07Eh9Xj0ERR9HU1P01zqUfeMOGJ3h8rN0ovmbNGsTGxuL8+fP4/PPP8cQT\nT+CJJ57w+I2kUlZWBbW6AEqlFmp1AcrKqgJdpLCSnZ2JoiI11OpCxMfVu3RMU/N5m98NhnMQW1TC\naDwP8oy5r0UNc7PZRgB6mCdaa3bp+PT0FOj1S5Camm+z3dxEp/JxacldgwbGeXys3Zq+JVWztLQU\nCxcuxL333ovCwkKP30gKnMclOFhq7VrlbqDSeS3y9KkWm9+Tk5NFa5SWpkVyn/0F0kdDrI24q+rq\n46h4+TPExhoxZcoSDBgwhMuWeilYxkDYDfrDhg3DokWLsHPnTjz77LNoa2vD9evX/Vk2p6ReS5Kk\nceFCm02GztChMTh8uPt+ycn9u28klzheIN255pbR1sFAqan5eOGFu/md8oKzCqrlhhD9bb3kZbHb\nvLN161ao1WqUl5cjPj4eLS0tePnllyUvkDs4j0vPZGqPQnn5i1i+fAfKyqqwbFkWmxF8bNmyLERH\nH/XJucwVKc7N4w37FdSd1htCefmLOHc+xeP3sDR1O+Nw5az777/f+vvQoUOD7nFb6rUkSSrmTijL\nU5llRKGzycLIdeYF0mvw0ku2KzT16nUMcOmBvR7mbCrz6E9WpLzjqILa0f9SANzIfnOX7ZNE1yY9\nW54//wUBqdeSJKlYOqGqsHfvD1AqtRztKQHLAumdb6Zjz8iBw67MnS/AHB42AqhhRcpL3SuoHYMR\nr127AKARwGYAX8CV7KquxJ4k7OnRQd/V6WQpGJkHZLW0bOFoTwl1nXlTO+8ktHYyP6qrj6O5ZSCA\nZtTjLnRMk/EYZsyYKnVRQ5ptBdV2MKJZPoAq1CPFOvY2MvII/t8dE8VPmJJi86u9JwlRvpv23/eC\nvHjUxWqFovuqDyI/CigEgIupBBvzkovi1yUh4UEubeml0tJKQa0uEBIS8ux8NQo8XsLSdtEVDxdR\nIXJbSorNFGrmmuPoTlvOA2hBY59TSIg52amW04Ftx4HjKOOnpWU8ysu12L17MW666T0MGzaEzXFu\nsjx1madJFtvD9v++O9lrYk3d9jDok89o33nH5ne7c4DfrkVKXxPKy7ufg23HgbNsWRZ2796I1lax\nV83XpbX1DdTWFqK2dg2b4zxkLwHF8hkD7vdNdm7q3rHD8b6cV4Ek4yi7immawSc7OxM5OSno1eux\nLq+sgnnNAwtzjZSpnJ4R+7+flPQ0pkwxQKHQQq0uRFGR+32T2dmZLq1FwZo+ScZRdhU74YNPWVkV\n9u6NwPXrGnSsoXsUgAKdl+nrXCNlc5z7xP/vz/bb/327K2cFA66c1fM5WpqPgoej1ZnMNwBLDdK8\ntCVXOvM9d6ZpcLSvxytnEfmCK4t1kzRcDSKOVmcCMpGQcALDhz+NH388g9bWjsW5OSbGd9yZR8zb\nOcf8EvRTUlIQGxuL3r17IzIyEnv37kVzczPy8vLQ0NBgd2F0InJd1/nzGxtjbRagtxcYbAf2mAcN\nAZEwD8wCpk8fie3b13R6avuMzXE+5s48Yt7OOeaXoC+TyVBRUYGBAwdat+l0OqhUKqxYsQLr16+H\nTqeDTqfzR3GIQk732l8BzFMpd7AXGDoG9nQMGrKIiFiMGTNuAcCnNl/q+hR2+vQvovuJ9Zl4O+eY\n35p3urYxlZSUoPJGsqpGo4FSqWTQJ/JQ99qf64GhI8uq+1TMJtMb2LMnuKZU7+nEmmeioxfDfNO1\nvamKpTB7O+eY32r6//Ef/4HevXvj8ccfx8KFC2E0GiGXywEAcrkcRqNR9FitVmv9t1KphFKp9EOJ\niXqW7rU/1wNDR5ZVpOgxu3fXcbF6HxJrnmltfQPR0Xlobe34fO31mSxbloXq6sdgMAyFOYTXoX//\nGgwalG4TL+3xS9D/97//jaFDh+KXX36BSqXChAkTbF6XyWSQyWSix7ryRxCFu+61vyyYO2OdT0Zo\nCeTm7J3u575ypR3l5Vmoq9thsz95xl7zzJgxyUhOdjWFOQ6dm+/69/8jHnooB9nZmXj++ecdvr9f\ngr5lSuYhQ4Zg9uzZ2Lt3L+RyOQwGA5KSktDY2IjExER/FIUoJHUfE5GJpKR3kJzsfNUrS/tyUlI8\nWloew/Xrmzu9ugrAEgA7uUCRj9hrnklO7m83/bVzH0BNzVE0NS2xed1geDV4OnKvXLmC9vZ2DBgw\nAJcvX0Z5eTlWr16NWbNmQa/XY+XKldDr9cjJyZG6KEQhS3zAzzynQaB7+3IVgDwAyQD6oyMn/yMA\nHIzlC+5OCS/WB9A5pdYiaDpyjUYjZs+eDQAwmUx4+OGHkZWVhWnTpiE3NxebN2+2pmwSkec8ya7p\n3r6ceeOn84AsADgDoIpzI/mAu6PRxefKXwvzNeo4Jmg6ckePHo1vv/222/aBAwdi165dUr89Udhy\nZXCW/XnYO9cazU080dEbsXTpEjv7kzvcuUE7v0ZViI7eiDNnkr1bLpGIei5XR23an/HxO5gXUWmH\npYlnzJiP2J7vI+5MuWDvGg0e/B3k8sfx008ytLZ+gMOHgcOHAWfLJQb1KiVBXjyioGW7qIb9RWpK\nSyuF1NRVNvtERy8SgEoucCMRsc88NXWV3QVqxPf/s1BaWmnnOnMRFaKwI94kYH9N4s7tyzNm/Abv\nvrsDdXXOc8bJfWLTXtTVRUKj2Qi9vntKrKM+gJdf/szt92fQJwpB4gtx21+TuGug6bqgOufZ8R17\n0140NQFVTcF9AAART0lEQVTLl3dvguvaFLR0acfN2n7znANSPL74SpAXjzxkeSxVKFZz3VUf6vy5\nZmQ8KSQlze/0yM81iYNFR5OM82virCmotLRSiI5+nM07FLy8nRaWxIl9rklJf8SUKQsxYMAwVFdz\nTeJg4Wzai87XxNmMmtnZmbjppvdQW2tZ9MZ52iaXSyS/sv+fmMvueUPsczUYXsXx45dw4cIpmEwX\nRY8Ty+0uK6uCWl0ApVILtboAZWVVkpQ5XGVnZ6KoSI1Bg46Kvt75mrgyo+awYUNgHlOhhe3YCnGs\n6ZNfeTstLImz97m2tMSipWUwgN/Dlbl4+CTmH9nZmdDrzW34trNtPo5Tp2Cd4M6VGTXFRvg6wqBP\nfuXttLAkzn6H3jkAb3b63dwMMHjwdygqetLnC3SQ6zpn5Zw5c8m6MlltbSZqa80327lzhzmdsqFr\ndo+lo94erpFLfiVWk0xNXYWiImaHeEN8fpZVAFoBvNZtf4VCi4oKbbftkyY9jtraITDXBzvWyrW3\nP/mGWl2A8vIXu20fNCgPSUnxMBjOITk5GcnJ/Z2uM801cimouDvvCLmm8+e6d+8JtLSMhHkkbbno\n/vba8n/44TpsV9z6o939yXfsNc81Nd2MpiYtACA+Pt9pwHcFa/pEIca21t99CUR7T1ZTpizBwYMb\nu52vT58cfPjhH3ljloAlB3/fvh/R0rKly6tVMK9TfDMsT11q9U670y9bsKZPFGa6Pk1dvPgzAOfz\n6tfXXxI9X58+MQz4Eigrq8KCBfobK2DFAVgM4I0br1YBeB/AB52OyMepU2e9fl8GfaIQ5GgWR0tK\nZvfJvq6K7h8Zed16nKuThJFzhYX/C4MhCR3Naea1DPr374W2tkswmf7V5Yi1MBof9Pp9GfSJwoij\nlMyUlP5oabFN6wRWYfTo/kzllID5yWpTpy3mtQwiI+dg9OgxN2bMtGVZhdAbHJxFFEYcDY5bs+b3\niI8/AXNapxZAIeLjT+CFFx7hoDpJ9LWzvQ+GDo0RfSU5ub/X78qgTxRGnA2Oi4rqY7M9KqqvS8eR\n+1JSxAP46NH9sWxZFlJT8222m/PzVV6/L5t3iMKIo8FxxcXlMBg222w3GMwdwn37imeDMJXTc2vW\n5GHBgj/CYHjVui0p6Wm88EKeS6nNnvaxMOiTy9iR1/M5WpTb3tzsbW298cwzd7u1mDc5l52diU2b\nugb22dbvlLPOeHt9LM74Jei3t7dj2rRpGD58OP71r3+hubkZeXl5aGhosC6KHh8f74+ikIfYkRca\nHNUgi4vtD+TioDppeLKYPeB4ugxn/BL0i4qKMHHiRFy8aJ7pT6fTQaVSYcWKFVi/fj10Oh10Op0/\nikIe4pwsPZsrT2mOngIAzwMU+Z43fSySB/1Tp07h448/Rn5+Pl591dx2VVJSgsobswJpNBoolUoG\n/SDHjryey9WnNNbmew5vJi6UPOg//fTTePnll3HhwgXrNqPRCLlcDgCQy+UwGo12j9dqtdZ/K5VK\nKJVKqYpKDnB2zJ7Lnac01uZ7BtunsgoAFUhI+BSDBo1xeqykQb+0tBSJiYnIyMhARUWF6D4ymQwy\nmczuOToHfQocZ4/+FHwsTTpff30KQAEsM2Za8Cmt5xJ/KluL7OxMvP/+Ow6PlTTof/nllygpKcHH\nH3+MtrY2XLhwAY888gjkcjkMBgOSkpLQ2NiIxMREKYtBPsBH/55FfKplS3aH+ZrxKS14eJIZ5/FT\nmSQr/4qoqKgQ7r33XkEQBOGZZ54RdDqdIAiCsG7dOmHlypWix/ixeOQFLnQePCzXIiEhT3TRbaDg\nxuLaf+Z1ChLOFj93l7O46degf9999wmCIAhNTU3CzJkzhbFjxwoqlUpoaWkRLxyDftDz9X9Y8pzt\ntVgtGvTj4jSCWl3A6xNEsrLyRa+VWl3g0fmcxU2/Dc5SKBRQKBQAgIEDB2LXrl3+emuSEFM5g4ft\ntRDveJ8xY4TT+djJv/ydGce5d8grTOUMHrbXIgsdbfhmvpq7hXzL35lxnIaBvMJUzuBhey0sT1mF\nSEg4genTR7LjPUi5mxnn7XQoDPrkFaZyBo/u1yITqanbUVT0GIN9ELNcm+eeW4jjxy8C6IvYWPEZ\nOH0xHQqDPnmFqZzBg9eiZzt/PhEtLW8BAFpagOXLuwdzX/ShcWF0IqIAU6sLUF7+osj2QpuOd6VS\ni8pKbbf9FAotKirM27kwOhGJsrQNnz79CwyGc0hOTsbQoTGcMjsAXE2I8EUfGoM++RXn5A8OHW3D\nagA7ALyJpibg8GFOmR0IrgZzX/ShMeiT33BO/uDR0TZcANuF0DnOIhBcDea+6Ldh0Ce/4UCu4NHR\nnMBxFsHAnWDu7UyoDPrkNxzIFTw6mhM4ziJY+Gtaa47IJb/hQK7gsWxZFlJT8+Fs5G5ZWRXU6gIo\nlVqo1QUoK6vyf2HJp1jTJ7/hQK7g0dGcsBOnTp2F0fgghg4diuTk/tZmBfbBhCbm6ZNflZVVYcOG\nnZ3aLVUMIEHK1dxx8g9XM9+Yp09Bhcvx9RzsgwkevnzqYps+EYliH0zwsJ/5ttPtczHoE5Gojs7e\nDpyeOTB8+dTF5h2iEOfpKGhO4BY8fPnUxaBPFMK8bQtmH0xw8GXmG7N3SFKcayewmIETOlzNfAt4\n9k5bWxsUCgWuXr2Ka9eu4b/+67+wbt06NDc3Iy8vDw0NDUhJScHWrVsRHx8vdXHIj5jnHXjMwAkd\nlu9McXE52toiUFxcbrPdVZIH/aioKHz++efo168fTCYTfvvb3+KLL75ASUkJVCoVVqxYgfXr10On\n00Gn00ldHPIjzrUTGJ2frmpqjgKoQsfyiWbMwOl5fFWJ8kubfr9+/QAA165dQ3t7OxISElBSUoLK\nykoAgEajgVKpZNAPMaxl+p9YYIiIWAyTCbAEfo6C7hm6No3+8ksL6uo22uxTV7cWzz23JPiC/vXr\n1zFlyhTU1dXhiSeeQFpaGoxGI+RyOQBALpfDaDSKHqvVaq3/ViqVUCqVfigx+QLzvP1P7OnKZHoD\ngwc/iLS0z5iB00OI3byjojSi+9bU/IiHH34UY8eOcuncfgn6vXr1wrfffovz589DrVbj888/t3ld\nJpNBJpOJHts56FPPwrl2/M/e01Va2gTrcnoU/MRu3m1tI0T3vXZtOpqagPfe0wIAnn/+eYfn9mvK\nZlxcHLKzs3HgwAHI5XIYDAYkJSWhsbERiYmJ/iwK+QHzvP2PT1ehQfzmnQWZbAEEYVOnbasA3IO2\nts9cPrfkQf/s2bOIiIhAfHw8WltbsXPnTqxevRqzZs2CXq/HypUrodfrkZOTI3VRKACY5+1ffLoK\nDeI370zExOhw6VIhgN4A2gHcAyATUVGuT8cgedBvbGyERqPB9evXcf36dTzyyCOYOXMmMjIykJub\ni82bN1tTNonIO3y6Cg32bt5z596Ld989jbq6NTbb3bmpc3AWEVEQsjcYy9kgLWdxk0GfiCiEOIub\nnGWTiCiMMOgTEYURBn0iojDCoE9EFEYY9ImIwgiDPhFRGGHQJyIKIwz6RERhhEGfiCiMMOgTEYUR\nBn0iojDCoE9EFEYY9ImIwgiDPhFRGGHQJyIKIwz6RERhhEGfiCiM9PigX1FREegikAO8PsGP1yi4\n+fr6SB70T548ibvuugtpaWmYNGkSiouLAQDNzc1QqVQYN24csrKycO7cOY/Oz/+wwY3XJ/jxGgW3\nHhf0IyMj8dprr6G2thZ79uzBxo0bcfToUeh0OqhUKhw7dgwzZ86ETqeTuihERGFP8qCflJSE9PR0\nAED//v1x88034/Tp0ygpKYFGowEAaDQabNu2TeqiEBGFPZngaNl0H6uvr4dCoUBNTQ1GjhyJlpYW\nAIAgCBg4cKD1d2vhZDJ/FY2IKGQ4CusR/irEpUuXcP/996OoqAgDBgyweU0mk4kGeD/ej4iIwoJf\nsnd+/fVX3H///XjkkUeQk5MDAJDL5TAYDACAxsZGJCYm+qMoRERhTfKgLwgCHnvsMUycOBF/+MMf\nrNtnzZoFvV4PANDr9dabARERSUfyNv0vvvgCmZmZuOWWW6xNOOvWrcP06dORm5uLEydOICUlBVu3\nbkV8fLyURSEiIqGHMhgMwpw5c4SbbrpJmDp1qnD77bcLH330UaCLFdZ69eolpKenC5MmTRIeeOAB\n4cqVK4Ig8Fr5iy8+/61btwoTJ04UevXqJRw4cMC6/euvvxbS09OF9PR0YfLkycKWLVv88jeFM8v1\nnDx5sjB79mzh4sWLvjlvoG86nhAEATk5OVAqlairq8P+/fuxZcsWnDp1KtBFC2v9+vXDwYMHcfjw\nYfTp0wdvvPEGAPBa+YkvPv/Jkyfjo48+QmZmZrftBw4cwMGDB1FeXo4lS5agvb1d8r8pnFmuZ3V1\nNWJjY/Hmm2/65Lx+y97xpc8++wx9+/bFokWLrNtGjhyJp556KoClos7uvPNOVFdX81oFiKef/4QJ\nE0S3R0dHW//d2tqKuLg49O7d27eFJrtuv/12HDp0yCfn6pE1/draWkyZMiXQxSA7TCYTPvnkE9xy\nyy28VgEg1ee/d+9epKWlIS0tDa+++qpPzknOtbe3o7y8HJMmTfLJ+Xpk0O+a0//UU08hPT0d06dP\nD1CJCDDXADMyMnDrrbdi1KhRmD9/frd9eK2kI/XnP336dNTW1uKbb77B8uXLcf78eV8Um+ywXM+h\nQ4fi5MmTWLx4sW9O7JOeAT/79NNPBYVCYbPt7NmzQkpKSmAKRIIgCEL//v27beO18h9PPv958+YJ\n6enpQnZ2ts0+SqXSpiO3q7vvvlvYv3+/94UmuyzX88qVK8Kdd94pfPjhhz45b4+s6d99991oa2uz\ndlQBwOXLlwNYIrKH1yqwnH3+b7/9Ng4ePIjS0tJuxwqdsrnr6+thMpkAAA0NDfjhhx8wduxYCUtO\nFtHR0SguLkZ+fr5PZinokUEfALZt24bKykrcdNNNuO222zBv3jy89NJLgS5WWLM3VxKvlX/44vP/\n6KOPMGLECOzZswfZ2dn4z//8TwDA7t27kZ6ejoyMDDzwwAP429/+htjYWMn+FrK9nunp6RgzZgy2\nbt3q/XkFX9w6iIioR+ixNX0iInIfgz4RURhh0CciCiMM+kREYYRBn8gFr7/+OlpbW62/Z2dn48KF\nCwEsEZFnmL1DdIPlqyCW+jh69Gjs378fgwYN8nexiHyKNX0Ka/X19Rg/fjw0Gg0mT56Mxx57DLfe\neismTZoErVYLACguLsaZM2dw1113YebMmQCAlJQUNDc3o76+HjfffDMWLVqESZMmQa1Wo62tDQCw\nb98+3HLLLcjIyMAzzzyDyZMnB+rPJLJi0Kew9+OPP2LJkiWoqanBK6+8gn379uHQoUOorKxETU0N\nli1bhuTkZFRUVODTTz8FYPs08OOPP+Kpp55CTU0N4uPj8Y9//AMA8Oijj+Ktt97CwYMHERERYXfw\nFJE/MehT2Bs1apR1ArIPPvgAU6dOxZQpU1BbW4sjR444PX706NG45ZZbAABTp05FfX09zp8/j0uX\nLuG2224DADz00EM+GUJP5K0eOZ8+kS/FxMQAAI4fP45XXnkF+/fvR1xcHB599FFrU40jffv2tf67\nd+/eNh2+Fgz4FCxY0ye64cKFC4iJiUFsbCyMRiM++eQT62sDBgxwK1snLi4OAwYMwN69ewEAW7Zs\n8Xl5iTzBmj6FPUtb+29+8xtkZGRgwoQJGDFiBH77299a91m0aBHuueceDBs2zNqu3/X4rr9v3rwZ\nCxcuRK9evaBQKBAXFyfxX0LkHFM2iSRy+fJla9ORTqeD0WjEa6+9FuBSUbhjTZ9IImVlZVi3bh1M\nJhNSUlLwzjvvBLpIRKzpExGFE3bkEhGFEQZ9IqIwwqBPRBRGGPSJiMIIgz4RURhh0CciCiP/H6lr\nSG8xl4VxAAAAAElFTkSuQmCC\n" } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is written:\n", "\n", "$$S_i = b_0 + b_1 1(Ra_i = \"PG\") + b_2 1(Ra_i = \"PG - 13\") + b_3 1(Ra_i = \"R\") + e_i$$\n", "\n", "$1(Ra_i = \"PG\")$ is a logical value that is one if the movie rating is \"PG\" and zero otherwise.\n", "\n", "Average values\n", "\n", "$b_0$ = average of G movies\n", "\n", "$b_0 + b_1$ = average of PG movies\n", "\n", "$b_0 + b_2$ = average of PG-13 movies\n", "\n", "$b_0 + b_3$ = average of R movies" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To do regression with factor variables in Python we use `statsmodels`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from statsmodels.formula.api import ols\n", "\n", "lm1 = ols('score ~ rating', movies).fit()\n", "\n", "lm1.summary()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: score R-squared: 0.020
Model: OLS Adj. R-squared: -0.002
Method: Least Squares F-statistic: 0.9182
Date: Sat, 06 Apr 2013 Prob (F-statistic): 0.434
Time: 17:57:22 Log-Likelihood: -569.90
No. Observations: 140 AIC: 1148.
Df Residuals: 136 BIC: 1160.
Df Model: 3
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [95.0% Conf. Int.]
Intercept 67.6500 7.193 9.405 0.000 53.425 81.875
rating[T.PG] -12.5929 7.849 -1.604 0.111 -28.114 2.928
rating[T.PG-13] -11.8146 7.411 -1.594 0.113 -26.471 2.842
rating[T.R] -12.0200 7.476 -1.608 0.110 -26.803 2.763
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 3.313 Durbin-Watson: 2.100
Prob(Omnibus): 0.191 Jarque-Bera (JB): 3.327
Skew: 0.367 Prob(JB): 0.189
Kurtosis: 2.821 Cond. No. 14.0
" ], "output_type": "pyout", "prompt_number": 5, "text": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: score R-squared: 0.020\n", "Model: OLS Adj. R-squared: -0.002\n", "Method: Least Squares F-statistic: 0.9182\n", "Date: Sat, 06 Apr 2013 Prob (F-statistic): 0.434\n", "Time: 17:57:22 Log-Likelihood: -569.90\n", "No. Observations: 140 AIC: 1148.\n", "Df Residuals: 136 BIC: 1160.\n", "Df Model: 3 \n", "===================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "-----------------------------------------------------------------------------------\n", "Intercept 67.6500 7.193 9.405 0.000 53.425 81.875\n", "rating[T.PG] -12.5929 7.849 -1.604 0.111 -28.114 2.928\n", "rating[T.PG-13] -11.8146 7.411 -1.594 0.113 -26.471 2.842\n", "rating[T.R] -12.0200 7.476 -1.608 0.110 -26.803 2.763\n", "==============================================================================\n", "Omnibus: 3.313 Durbin-Watson: 2.100\n", "Prob(Omnibus): 0.191 Jarque-Bera (JB): 3.327\n", "Skew: 0.367 Prob(JB): 0.189\n", "Kurtosis: 2.821 Cond. No. 14.0\n", "==============================================================================\n", "\"\"\"" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we plot the fitted values:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plot(jitter(r.labels, 2.), movies['score'], 'ob')\n", "xticks(np.unique(r.labels), r.levels)\n", "xlabel('rating')\n", "ylabel('score')\n", "\n", "plot(range(4), lm1.params[0] + np.append(0, lm1.params[1:]), 'sr', markersize=10);" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEMCAYAAAAoB2Y1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4U1W6P/BvoEjLpbQUmnDTarmXQgG5OOfYRrGN5+kM\nU8exiKMTRUCR2zjnATwUplFEqpxRW0SdGfhpHGcG8Bn1dNoRimJT5lEsYIerChbKrU2Btlx7gZT9\n+yMkJM3OpeneO0nz/TxPH8vO3slqlnmz9lrvWkslCIIAIiIKC10CXQAiIlIOgz4RURhh0CciCiMM\n+kREYYRBn4gojDDoExGFEcmC/qxZs6BWq5GcnGw/Vl9fj/T0dAwfPhwZGRm4cOGC/bE1a9Zg2LBh\nGDlyJEpKSqQqBhEReSBZ0H/qqaewdetWp2N5eXlIT0/HkSNHMG3aNOTl5QEADh8+jM2bN+Pw4cPY\nunUrnnvuOdy4cUOqohARkRuSBf17770XsbGxTscKCwuh1+sBAHq9Hp9++ikA4P/+7/8wc+ZMdOvW\nDQkJCRg6dCjKy8ulKgoREbkRIeeT19bWQq1WAwDUajVqa2sBANXV1Zg6dar9vMGDB+PMmTMu16tU\nKjmLR0TUKXlaaEGxgVyVSuUxiLt7TBAEn35yc3N9Ppc/yv6wboL7h/UT3D/trR9vZA36arUaZrMZ\nAFBTU4P4+HgAwKBBg3Dq1Cn7eadPn8agQYPkLAoREUHmoD99+nQYjUYAgNFoRFZWlv34pk2bcO3a\nNRw/fhxHjx7F5MmT5SwKERFBwj79mTNnwmQy4fz58xgyZAheeuklvPDCC8jOzsbGjRuRkJCALVu2\nAABGjx6N7OxsjB49GhEREXj77bc73H+v1Wol+CtIDqyb4Mb6CW5S149K8KUTKEBUKpVPfVRERGTl\nLW5yRi4RURhh0CciCiMM+kREYYRBn4gojDDoExGFEQZ9IqIwIuvaO0RE1H7FxWUoKChBS0sEune3\nYNGiDGRmpkry3Az6RERBpLi4DIsXb0Nl5Wr7scrKHACQJPCze4eIKIgUFJQ4BXwAqKxcjXXrtkvy\n/Az6RERBpKVFvAOmubmrJM/PoE9EFES6d7eIHo+MbJXk+Rn0iYiCyKJFGUhMzHE6lpi4HAsXpkvy\n/FxwjYh8ImdGCTkrLi7DunXb0dzcFZGRrVi4MN3n99pb3GTQJyKvxDJKEhNzkJ+vY+APMlxlk4g6\nTO6MElJOWOTp87aUqGPkzigh5XT6oC/3RAeicCB3Rgkpp9N37/C2lKjj5M4oIeV0+pY+b0uJOs52\nV7xu3UqHjJIHebccghQJ+vn5+diwYQMEQcCcOXOwePFi1NfXY8aMGThx4oR90/SYmBjJX5u3pUTS\nyMxMZZAPUo7jlt7I3r1z8OBBbNiwAbt378a+fftQVFSEyspK5OXlIT09HUeOHMG0adOQl5cny+vz\ntpSIOjPbuGVJycswmQxez5e9pf/9999jypQpiIyMBACkpaXh73//OwoLC2EymQAAer0eWq1WlsDP\n21Ii6szExi09kX1y1vfff4+f//zn+PrrrxEZGYkHHngAd999N/785z+joaEBACAIAvr27Wv/t71w\nKhVyc3Pt/9ZqtdBqtXIWl4gopKSkPIl9+xIcjrzocXKW7C39kSNHYtmyZcjIyEDPnj2RkpKCrl2d\nB1FVKhVUKpXo9QaDQe4iEhGFLLV6MACDw5EXPZ6vSMrmrFmzsGfPHphMJsTGxmL48OFQq9Uwm80A\ngJqaGsTHxytRFCKiTkVs3NITRbJ3zp49i/j4eJw8eRIff/wxdu3ahePHj8NoNGLZsmUwGo3IyspS\noihE5AFnr4eetuOWN4dK3VJkwbXU1FTU1dWhW7dueOONN3Dfffehvr4e2dnZOHnypNuUTS64RqQc\nLqrWOXCVTSLyiU63AiUlL4scX4mtW1cFoETkD66ySUQ+4ez18NDpl2EgIvcc+/APHvwOQBkA564c\nzl7vXBj0icKUWB9+RMSzsFgAW+C3zl5/MDAFJFmwT58oTLnrw+/X71EkJY1s9zZ9FBy8xU229InC\nlLs+/KSkkSgtNShbGFIMB3KJwhRXoA1PDPpEYYor0IYn9ukThbHi4jKsW7fdYQVa9uEHE39mSHNy\nFhFRCPJ3hjQnZxERhSC59vdm0CciCkJyzZBm0CciCkJyZVcx6BMRBSG5sqs4kEtEFKT8ya5i9g4R\nURjhMgxEBIC7YpEVgz5RGBDL+a6stPYXM/CHFw7kEoUBuXK+KfQoEvTXrFmDpKQkJCcn47HHHkNL\nSwvq6+uRnp6O4cOHIyMjAxcuXFCiKERhibtikY3sQb+qqgp/+tOf8O233+LAgQNobW3Fpk2bkJeX\nh/T0dBw5cgTTpk1DXl6e3EUhClv+5nwXF5dBp1sBrdYAnW4FiovL5CgeKUj2Pv3o6Gh069YNjY2N\n6Nq1KxobGzFw4ECsWbMGJpMJAKDX66HVahn4iWSyaFEGKitz2qzj4nlXLI4DdE6yB/2+ffviv//7\nv3H77bcjKioKOp0O6enpqK2thVqtBgCo1WrU1taKXm8wGOy/a7VaaLVauYtM1OnYgvS6dSsdcr4f\n9Bi83Y8DrGTQDyKlpaUoLS31+XzZ8/QrKyvxs5/9DDt37kSfPn3wyCOP4OGHH8bChQvR0NBgP69v\n376or693Lhzz9IkCRqs1wGQyuBxPSzNwZ60gFvA8/T179uAnP/kJ4uLiAAC/+MUv8PXXX0Oj0cBs\nNkOj0aCmpgbx8fFyF4WI2oE7awUfKeZayB70R44ciVWrVqGpqQmRkZH4/PPPMXnyZPTs2RNGoxHL\nli2D0WhEVlaW3EUhCnvtCRr+jANQx7mrI6nGWGQP+uPGjcOvf/1r3H333ejSpQsmTJiAuXPn4vLl\ny8jOzsbGjRuRkJCALVu2yF0UorDW3qDhzzgAdYynOpJqjIVr7xCFCZ1uBUpKXhY5vhJbt64KQImo\nLU911Nzc1acxFu6cRUQAOEErFHiqI6nGWBj0icIEB2aDn6c6kmp9fS64RhQmODAb/DzVkVRjLOzT\nJwoj/mzKQcpyrKPLl89BEFoQHT3Y5xRNbqJCRB3CdfgDQyyTJzExB/n5Oo/vf8AnZxFR6OL6O4Ej\n1zIYHMglxXDFxtDDdfgDR65sK7b0SRFsMYYmpnkGjlzZVmzpkyLYYgxNTPMMHKlSNNtiS58UwRZj\naGKaZ+DItQwGgz4pgi3G0MT1dwIrMzNV8veaKZukCPH0s+XIz2cAIZIS8/QpaHBikHKYWx++GPSJ\nwoy/k3qoc2DQp6DCFqj8uIRyeOOMXAoazNVXBjOlOge5GkgM+qQYuaaVkzNmSoU+g+FtvPaaCU1N\nowBYAGSgsnIbgI43kDg5ixTDFqgy5JrUQ8ooLi7Da6/tR1PTZgAGAC8D2IbKSp0kkxnZ0ifJubst\nZQtUGcytDx1in5WCghI0Nb3b5szVAFZK0kCSPej/8MMPePTRR+3/PnbsGFatWoXHH38cM2bMwIkT\nJ+wbo8fExMhdHJKZp357zu5UjhyTekha7j4rPXo0urmiq72B1JH+fkWzd27cuIFBgwahvLwc69at\nQ79+/bB06VK8+uqraGhoQF5ennPhmL0TcrxljjBXn8jK3WclLm4G6uo2uxyPipqBjz6aDwAeU3KD\nKnvn888/x9ChQzFkyBAUFhbCZDIBAPR6PbRarUvQp9Djrd+eLVAiK3efFY0mBjExznfEUVHPYOnS\nNGRmpkKnWyGaEPG7381BQUGJ19dVNOhv2rQJM2fOBADU1tZCrVYDANRqNWpra0WvMRgM9t+1Wi20\nWq3cxaQOYL99aOL8CeW5+6wMHhyPhQvT24zJ/MpeH65fFqUAjNi//yAslkzvLywopKWlRejXr59w\n9uxZQRAEISYmxunx2NhYl2sULB5JpKjIJCQmLhcAwf6TmPg/QlGRKdBFIzfE62w560xm/n5WMjJy\nnK6x/jge8xw3FWvpf/bZZ5g4cSL69+8PwNq6N5vN0Gg0qKmpQXx8vM/PZXjySaCqyvuJCQkwvP++\nX+Ul/zBzJPRw/kRg+PtZEUuIiIw8heZm315XsaD/t7/9zd61AwDTp0+H0WjEsmXLYDQakZWV5fuT\nVVXBcHM8wBODH+WkjmO/fWhx7S4oA1CCXbtOQ6dbwa4eiXW0K03sy+Ls2V6oqPDxCaS8XXHnypUr\nQlxcnHDp0iX7sbq6OmHatGnCsGHDhPT0dKGhocHlOnfFy01La3tvI/qTm5Ym019E/ioqMgkZGTlC\nWlqukJGRwy6EIODcXWASAHb1yEWurjTn5/Uc1oO605xBv3Nh37H8cnPXC3Fx2UKfPnohLi5byM1d\n7/Ua53oR6y8WBJ1uhQKl7/zE++OleX+LikyCTrfCa9DnMgykGO6TKy+D4W2sXr0fdXWbcfHi+6ir\n24zVq/fDYHjb43WZmanIz9dBp1uJPn1Oi57DpTKkIZ6mWYby8qPQag3Q6VaguLjMr+fOzEz1aRVV\nBn1SDNfekddbb5lgsThP37dY3sX69d6DiC1gTJkyWOTRMhw8+F2HgxKJpWmWAdiGhoZNMJkMKCl5\nGYsXb5P1PWbQJ8Uwh19eFkuU6PHr1yN9fg7XxdrKEBHxV9TVbVYsKHVmru9vCazr6twi990vF1wj\nxXDtHXlFRDSJHu/WzcdcPrhmhnz77X5cvvyx0zlM5/Rf2/d3//5TaGhwPU/Ou18GfVIMc/jltWBB\nGlavfvZmF4817VKlOobo6NtQXFzm8/tsS7ktLi7DL37xg+g57JLzn2NKs3X9Hddz5Lz7ZfcOKU64\nuRiU7b8kDYPhOeTkjEXv3plQqT4A8DIE4a84duwDv7pkVq78ANeu3RB9jF1y0gjE3gchuUcuZ+SG\nJm7YrQyp9sjt2/dRNDQ8B2AbHPudVaqn8Y9/6FlnEnFcefby5XMQhBZcu9YVZvMFDBw4EAMG9Gzf\n0snBtMqmVBjIQxOn+ytDuiyp7gBs9bISQFcArejR4zzrS0KO3WnWRpEO1i/aP6CuDjhwQNq9pNm9\nQ4phyqYypMqSSkjodfO3VACrYF3YZBVGjBBL66SOutUokjejp9MG/eLiMuh0K5hbHESYsqkMqfqJ\nV62aAY3mt07HbrttNlpaLPxMyeBWo0jexlFIdu9442nLPt6WBg5TNpXhKUvKtthXTc1VVFdXQ6OJ\nwaBB/UX7jDMzU7Fhg/V5qquv4Mcfq9HUNB+HDqXi0CF+pqR26dK5m7/J2zjyOJDb2NiIU6dOYcSI\nEZK8WHv5u12iVANZJD1ul6gsxxUdL106h5qaZpjNGx3OyAGgQ2LiNo8D6vxMyau4uAyzZxthNmsA\n2Pr0nRtH+fm+pTf7PZBbWFiIJUuWoKWlBVVVVaioqEBubi4KCwvb87cEBPuOgxeXXVaO2B2vNciX\n4dYA7WoAK+0D6gBEl/09c+YcgBWwhgwLgAwAqfxMSaSgoOTml3EZgO0AzgN4FL16qfAf/zFU0vks\nboO+wWDAN998g/vuuw8AMH78eBw7dkySF5Ub+46JxLOlbEH+VtAHrFk5QHX1lTZfEmXYuXM94uP/\niFOnugNwbOlbu3b4mZLGrYZqKhzrZuJEA7ZuNdjHKD2twW+7q/PGbdDv1q0bYmJinI516RIa477s\nOyZyf8drC/K3WAN3dXU16uo23zxmXQisqWkzTpxYAeeADwCrERU1AwsXzpesvOFCbBMVTw1VX8Yo\nnc9p+0XvzG3QT0pKwl/+8hdYLBYcPXoUBQUF+MlPftLevy8gON2fyP0dry3IWy0H8CASE5cjMjIG\ndXW2445pg+JhYujQgfxMtZO7AP7444PcNlR9md8iflcnzm3Qf+utt/Dyyy+je/fumDlzJnQ6HVau\nXNmuPzCQ2HdM4U7sjlejeR4DB17B9evPo6amBmp1HwwevN0eXA4dsp1pCw1lAL4Tff6BA3uJHif3\n3AXwXbtWIj9fJ9pQXbt2h+hzOY6nuL+rcyV6psViQWZmJr788ku88sorPj8ZEQUP8Tvehzw2hm59\nSVhg6+IB5sPah8/u0o7ylGTirqHqyxil+7s6V6IliIiIQJcuXXDhwgWXfn2lcWNmIv+5CyTuNufe\nvfsg3nprBpqauqKx8VUAy9B2KYZ+/b5Hfv5z/Ez6wZ8kE1/GKMXOccftPUHPnj2RnJyM9PR09OzZ\nE4A1/7OgoMDrk7Z14cIFzJ49G4cOHYJKpcJ7772HYcOGYcaMGThx4gQSEhKwZcsW0S+YkpKXOQmE\nSELu+pV37z6IDz884zCYC9iydByzSpKSDPws+sldAJ86dbDb7Bxfxigdz9m2zXMZ3E7Oev/momYq\nlQqAdRlclUoFvV7f7j9Ur9cjLS0Ns2bNgsViwdWrV7F69Wr069cPS5cuxauvvoqGhgbk5eU5F06l\nAmAtHieBEEnD3USruLgZbQK+zUpY196xXc/PYke0naA4deoAfPjhGclWn/U6qdXTrunNzc3C/v37\nhf379wvXrl3za4f2CxcuCHfeeafL8REjRghms1kQBEGoqakRRowY4XIOAPtu8WlpuX69PgW/oiKT\nkJGRI6Sl5QoZGTlCUZEp0EXq1JKTf2P/XDn+9Oz5K9HjQK7998TE/2H9SCwjI0f0fdfpVvj1fF7C\nuuC2e6e0tBR6vR533HEHAODkyZMwGo1IS0tr17fO8ePH0b9/fzz11FPYt28fJk6ciDfffBO1tbVQ\nq9UAALVajdraWjfPYAAAnD79JUpLtdBqte16fQpuXCdJedXV1aLHW1ouih7v1+97JCUZmPosk46u\nIFBaWorS0lLfX9Ddt8H48eOF77//3v7vH374QRg/fny7v3V2794tRERECOXl5YIgCMLixYuFFStW\nCDExMU7nxcbGulyLmy19ti46L6lbOeRdUtJcAVje5j3/H+H2238pJCY6H+dnT37uPgNxcdl+vfce\nwrogCB5a+haLxWmhteHDh8Ni8T0tyGbw4MEYPHgwJk2aBAD45S9/iTVr1kCj0cBsNkOj0aCmpgbx\n8fGi1+t0K9m66MS4TpLyBg3qj0OHMuC4MQrwIEaN6oqFC9M5qVFh4pk3y1FXNx+LF1tHZaWsA7dB\nf+LEiZg9ezYef/xxCIKAv/zlL7j77rvb/QIajQZDhgzBkSNHMHz4cHz++edISkpCUlISjEYjli1b\nBqPRiKysLNHrOWAUXNyl+nl7zN3jXCdJedYg03bbyuX2AM8gryzb+63Xz0Bd3SjYvoSBVFRWpkq+\ns5zboP/OO+9g/fr19hTNe++9F88995xfL7Ju3Tr86le/wrVr15CYmIj33nsPra2tyM7OxsaNG+0p\nmxTcPPW/A/DYN+/P9HOSh7cUQG9f3iS9zMxUjBmzAyaTweUxqe963aZsXr16FZGRkeja1fqCra2t\naGlpQY8ePSQtgMfC+bmePsnD05rqgiB4XG/d+doyWNd2iUBc3HdYsCANu3bVcI39IMDN6wNnwoTZ\nqKjQoO3y1e1NkfV7Pf37778fX3zxBXr1sq6v0djYCJ1Oh6+++srnF6fOxZ/+d9tjt661Te23BpW6\nOuDDDxlUpOZva52b1wdGcXEZamqi0Xb5ao3mfSxc+KSkr+U26Le0tNgDPgD07t0bjY2Nkr44hRZP\n/e/uWha2vnnrtWUA1gNwngDEoCItT91wuz/6f0BVldtro/5dhTT8C1VIwAm8bz/OgXV5WTdReR13\n4EkkoMp+vPeVo9i99hh2r21zQUICDDcn0LaX26Dfo0cP7N27FxMnTgQA7NmzB1FRUX69CHUO3tYA\n8fTYPfcMxI4df4XFMkr0uRlUpOOptT61uQoGk8nLM5yAFsAJhyMcWJeX7U44AVUohUP9XAFgcp1X\nYejAa7kN+vn5+cjOzsaAAQMAAGazGZs2berAS1Go82UNkFfmjkP/qxfRpQswqEcsdq/9CrvXAvv3\nH8d/WPoCOAOg1KUlyaAiHeeuNOvYCWDB6dPngH6+PsutiVocWJdfe1bJ7Ci3Qf/48eOoqKjAiRMn\n8PHHH6O8vDxkds4i+XhK6cvMTMXuYbEwmPZbDzScaHPGyZv/NTm1JBlUpHWrK815c+1jx55F1RXx\n2bhtde1Si6RRz2Dw4Hjm6ivAdheNSvlfy20UX7VqFaKjo3Hx4kV8+eWXmDdvHubNmyd/iXxk2zNS\nqzVAp1uB4uKyQBeJANTVi0/ld1UFwIBu3X6G/HwGFSktWpSBqKj1aLttXlPTYzh16opPz9F6YzgG\nD47H1q2rWDcKyMxMRX6+Dn1jq2R/LbctfVuqZlFREebMmYOf/vSnQbNzFtdrCV5nTtf7eGYCAANG\njnyedSaxzMxUDB36CQ4caPtICVpvDAdQ49PzcJxFWZmZqdg9NgEwtb1DFudvhpbboD9o0CDMnTsX\n27dvxwsvvIDm5mbcuHHD5z9ATkwrC143bqjadb5tyz1OCJLWgAE9RYK+71vqARxnCWZ19RddGr5l\nZfMwatQHXq91+3/Bli1bsHXrVixZsgQxMTGoqanB2rVt84YCg+u1BK8uXXyfTGfry+edm/TEMq2i\nor4Dmny7PiryOBYsfEmm0lFHnTldj8oG54Zvc/NMVFSs93qtx52zHn74Yfu/BwwYYM/kCTSu1xK8\nBg3uCzSc9Hpe39gqLMl/CZmZqdDpVvDOTWJimVZTp6bhm1dfBZq9Xz90WF++9wpxvMvtu/+4T9e4\n3lHbBu43A/C8pE377veChC97RlJgxPXt49N5Y8cm2IMK79yk462b7Mn333ROwHfD13qk9mlbP/fc\nM9Bp16w0lOJWlpt7rnfUJWg7cO9OSAZ9X/LFKXS4u3O7fPmcwiUJbb50kx29HgUtXDdC6trlMEYn\nDboV7BMSZC9vuBGrn507n0VT02P2f1chAVpY74THjk1w+1x9I6KQWOXY8G1HKG/3Cv0KCvLikYjc\ntDSx/fZcfnLT0uzXFBWZBI3meZdNPTSaWdzAox28bUhTVGQSYmN/LXpOcvJvAlz6zs9d/QArXI75\nsj1sUZFJmDDhOSEy8tcCkO1wvZ+bqBD5JSEBBgD//ncVLlxMcHk4pk8VUlISnFqSmZmpGDDgA5jN\nzpt6mM3SryXemXnqJrO1MhsahoieY8uiAphJJRd39WP9f96ZL+OTtomSxcVl+N3v/ozDh+ehufkd\nr9cx6JOkbItA6XQrYBJbannqShhElomNjh4MsRVF2K/vO08JDrfSnMsA5MCx/9dxPIyZVPJxVz9R\nUd+hySGrqr3jk47Bf926ldi2zfP5XFeBZLFoUQYSE3Ocjln/Z04XPZ8ZWR3n6T2/1cpMBaCDdatE\nA3r2/AWio89i7dod0OlWYOXKD9xkUm1X4k/o1NzVz9KladDpViItzQCdbqXfM9QzM1N9WnefLX2S\nRXsH25mR1XGe3vOCghKHM1Nv/pThxg0zKiretT8SGTkP1rsB53riHVfHBUsCituds4IBd84KL9bb\n0+3cQUsGYt02UVEz0NS0WeTslQCcW4y23ZvY399x7XkP/Xm//d45i0hp3JRbOmLBIj9f59TKrK4e\nKLJUAxAZeRLNDhO4OHNaOu15D2V7vyXIRPLqjjvuEJKTk4WUlBRh0qRJgiAIQl1dnfDAAw8Iw4YN\nE9LT04WGhgaX6xQqHlGnUlRkEhITlzulACYmLndJf3WXQnjXXTOEuLhsoU8fvRAXly3MnLlUyMjI\nuZnumSMAJtGUUPLOW1qtv+c68hY3FRnIValUKC0tRUVFBcrLywEAeXl5SE9Px5EjRzBt2jTk5eUp\nURSiTs/9goTOg7FiA4sazSw0Ng5EXd1mXLz4Purq5uOjjy6ipORlNDQYYd3DdRus/f5W7O/3nXja\nZhnKy4+6LBMv10x1xbp3hDZ9TIWFhTDd3LZNr9dDq9Uy8BNJwNdgITawePZsFCoqXnc4qwQWy7tw\nthrWfn/r9cyw8p1rlpp1zZyGhk2w7WJp68KRK6NNkaCvUqnwwAMPoGvXrnjmmWcwZ84c1NbWQq1W\nAwDUajVqa2tFrzUYDPbftVottFqtAiUmCl3ul7U4A51uhcugoC34FxeX4YknNra5yvOEImZY+a64\nuAznzjUgMlKP5uYhADIgtmaObbFBXzPaSktLUVpa6ntBOtI/5avq6mpBEATh7Nmzwrhx44SysjIh\nJibG6ZzY2FiX6xQqHlGnItanr9E85bLUhWM//61r2vYji/crx8Y+Kuh0K7hMho/E6iQy8lmhV69H\nRd9f2zIMRUUmQadbIaSl5fr8fnuLm4q09G1LMvfv3x8PPfQQysvLoVarYTabodFoUFNTg/j4eCWK\nQtTpee+2sW6YXlnZDXr9ehiN8DBjNwMREc86dfEkJi5Hfv48Zuy0g9g4S3PzO4iLm4ErIjtY2rpw\nvGW0iWVpeSN70G9sbERrayt69+6Nq1evoqSkBLm5uZg+fTqMRiOWLVsGo9GIrKwsuYtCFDYcp+YX\nFJTg2LGrNx9x3jC9rg5YvDgHkZHnbz5uCzDWdZBiY3/AokX3YtcurmjbEe7GWTSaGMTE+Dcp0VNK\np0dS3Lp4cuzYMWHcuHHCuHHjhKSkJOGVV14RBMGasjlt2jSmbBLJxLlLIcdjd01cXLbocaZjSsNT\n+qU/XTientNb3OSMXKJOou2t/rlzDQ7b59la+N0gtrBdcvLzaGzs4dLi9HcdGHIm1irv6Ps7duzz\nOHDgDZFHOCOXqNMTCyqRkXqHM2yBRXwP1YEDe2HhwvSArwvTWUm97k5xcRl+/LEatvEZayi3wJoR\n5Blb+kSdgE63AiUuS1mvgHUylaMyREX9FU1NbQdmGeBDibW+BwLYD8BxHsUcABvY0ifq7MQHCjPQ\nrdsMXL8+DLaWoEZTjWeeuZsDsyHOWt/VsAZ8x9a+2uu1DPpEnYC7CVkqVTScW/u/xaRJY2AwPKdI\nuch/nlbYtNZ3BNpmY1l53iCdQZ9kwSV4lSU2ezMqar3L0slm8+s+bUHJ+gssbytsLlqUgZ0716Op\nyQJvQb4tBn2SHJfgVZ7YQKG7pZO9LdjF+gs894vmrbTPwVi69CBefHFnu5+bA7kkOfFBxVsbcZD8\niovLoNdZQx7xAAAQt0lEQVSvR12d6yYpEybMR79+sfZW/D33DMTXX1c7pHqaUVGxweU61p9ytFoD\nTCaDy/G0NANKS28dj47OwuXLn7Y5iymbpDC5loQl39ha6nV189F2E3SN5nlUVzfj229tX8pl2LHj\nr07LLHDLxMBw7FI7ePA7iNVB2xU2f/vbDKxe/azISqjuMeiT5LjJeWC5dg1Yl1To1+97DBjQCxUV\njitpui6d3Nz8DhyXTrZh/clHrEvNuuYRYKsHseUZrAPyb2P9+kdx/XokunVrxvnz8IhBnyTHTc4D\n68yZc7Dm6DtO2ElFUpJB5GzxEOBuy0SSh1gfvsXyLiIifoaePf8fIiKa8PjjaaJjKgbDc07ZWCqV\n2L7HtzDok+Sknn1IvisuLsOxYyo4p2laB2EjI1tF+npPQ+wLYvToXujfn/WnFHddohbLRFy8aAAA\nfPhhDiZNKutwPXAgl6gTcTeIHhU1Ax99NB8AHLoRygB8CsBxp6wcaDQ12LDhSQZ5BVnrzbapiuMX\n8HYAtwbPJ0yYj717xZfSsPEWN9nSJ+pE3LUYhw4d6BTE161bifLyo2ho2NTmzNUYOHA+A74MPM19\nuOeegS4D6sCzAMY6Pcfhw1dQXNyx1j6DPlEn4m4QfeDAXvbfbXne1rRA13OvX79NdFtF8p+3uQ9f\nf10tkoHzLqwD6rc0N9+Odeu2M+gTkZW1xeicwhcR8QymTh3ncq67L4gff6zGgQO3BgM5MavjvE22\ncneHZtuL2Go5gAfR3LyjQ2Vh0CfqRKwtxsdgS9MEWmGx/Aq7dm13OVd86YZn0NQ0H46LeFVWqvC7\n3/2ZQb8DvM1dcfcFDHwP6/4HrQAeBJCKyEjXumwPBn2iTsQaXFLhOrHKtXUolmV1+jRw6BDQdhGv\nw4fndbgvOZx5m7si9gWs0TwPoCfMZoP9mBSpswz6RJ1IeyfGtd14W6dbgUOHStB2Ea/m5nd8WqiN\nxHmbuyKe5vyQ07HLl89AEG7D2rU7UFBQ4vdYC4M+BQxXcpReRyfGLVqUgbKyjU4Ts2y4DIP/fJm7\n0vYL2PG4pIvg+bHHb7tZLBYhJSVF+OlPfyoIgnVT9AceeMDjpug35w8oUTwKAOdNu60/iYnLfd4U\nmtzzd6Ntm/Hjn+Mm6UHG08bqbXmLm13a/ZXlh/z8fIwePRoqlQoAkJeXh/T0dBw5cgTTpk1DXl6e\nEsWgIOI+m6Fjg1Thznb31NzcFd27W7BwYXq7W4KrVs1AYmKO0zHr3UK6lEWldpByEUPZu3dOnz6N\nf/7zn8jJycHrr1tn/hUWFsJ0M0FYr9dDq9Uy8IcZrsQpPam6ALiMRvCRchFD2YP+888/j7Vr1+LS\npUv2Y7W1tVCrrXs5qtVq1NbWur3eYDDYf9dqtdBqtXIVlRTElTil5y0XvD3c9S9TYHgaqyktLUVp\naanPzyVr0C8qKkJ8fDzGjx/vtlAqlcre7SPGMehT58GVOKUn5d0TB9mDi7e7L8fG8IsvvujxuWQN\n+l999RUKCwvxz3/+E83Nzbh06RKeeOIJqNVqmM1maDQa1NTUID4+Xs5iUBBiF4L0pLp74naJwcnx\n7sv2pbx27Y52fykrtsqmyWTC//7v/+If//gHli5diri4OCxbtgx5eXm4cOGCaJ8+V9kMPmwBBi+x\nYJ2YuByPPz7YaTtEb3XG7S6Dm3g95yA/X4fMzNTgWmXT1o3zwgsvIDs7Gxs3bkRCQgK2bNmiZDHI\nT2wBBjexu6epUwfjww/PtKvOOMge3Do8diNHTqlUgrx4Yac9ucIUHPypM9ZzcEtLyxWtn7S0XEEQ\ngiRPnzoHtgBDjz91tmhRBvP0g1hHx264DAP5jGmWocefOuMge/AQG0PraOYbt0skn7kbKMzPZ0AI\nVqyz0OVpwBYA1q3b7vClbL0LKygoQUnJao9xk0Gf2qW4uMzlfzYGj+DGOgtN7cmicv6CCKLsHQp9\nnKkZelhnoak94zFiGT3ucCCXiFwUF5dhwoT56NtXj759H8WECbNRXFwW6GKFlfaMx7jfbtEVW/ok\nOU7gCm3FxWWYPftTmM3r7ccaGnIwe7YRGzZwToZS2jNg6367RVcM+iQpTuAKfQUFJTCbX29zdDXM\n5pVYt24761Eh7cmiEvuCcIdBnyQl5UqPFBjuuwq6iu6oRfLxdTzG8Qti2zbP57JPnyTFCVyhz31X\nQSvnZASxzMxUn9ZGYtAnSXECV+hbtCgDGs1vb/6rDMAKAHrcdtsBTJ06IIAlIykw6JOkOIU/9GVm\npmLDhizcddejUKk+APAyACOuXfsUH354hlk8IY6Ts0hynAzUOXCJZekpkdkWVEsrU3jgZKDOgeMz\n0gqWzDZ27xCRKI7PSMt9Ztt2RcvBoE9Eojg+I61guXNi9w4RuRUdXYvY2JkAbsOdd/bCSy/NYNed\nn4LlzoktfSJyYet/rqjYgIaGv6GhwYiLF2MCXayQFix3TszeoYDhGj3Bi5k78lAisy3g2TvNzc1I\nS0tDS0sLrl27hp///OdYs2YN6uvrMWPGDJw4ccK+OXpMDFsS4SJYMhlIXLD0P3c2tsw2W4Nn7dod\nKCgoUbbBI8O+vS6uXr0qCIIgXL9+XZgyZYqwc+dOYcmSJcKrr74qCIIg5OXlCcuWLXO5TqHiUQBw\n8+3g5lw/JgHIEYBcIS4uWygqMgW6eCGtqMgkJCYud/r/PjFxuWTvq7e4qchAbo8ePQAA165dQ2tr\nK2JjY1FYWAiTyQQA0Ov10Gq1yMvLU6I4FATYkgxut1Zt1AHYBsB6R1ZXByxezDuyjpByUULHLtJL\nl84BaPF6jSJB/8aNG5gwYQIqKysxb948JCUloba2Fmq1GgCgVqtRW1sreq3BYLD/rtVqodVqFSgx\nyS1YMhlInC346PXrUVe32ekxrpraMVI1eG51kaYDKAXQH8AOr9cpEvS7dOmCf//737h48SJ0Oh2+\n/PJLp8dVKhVUKpXotY5BnzqP9mwQQYGRmZmKMWN24OYNuRPekflPqgaP8x2D9uZ/DQDEY6mNonn6\nffr0QWZmJvbu3Qu1Wg2z2QyNRoOamhrEx8crWRQKsPZsEEGBwzsy6UnV4GnPFomOZA/658+fR0RE\nBGJiYtDU1ITt27cjNzcX06dPh9FoxLJly2A0GpGVlSV3USjIcI2e4Mc7MulJ1eBpzxaJjmTP0z9w\n4AD0ej1u3LiBGzdu4IknnsCSJUtQX1+P7OxsnDx50m3KJvP0iQKPq6YGJ7G0Z2A5gDUe4yYnZxER\nhSjHL+TLl88BuIZvv93AoE9EFC68xU2uvUNEFEYY9ImIwgiDPhFRGGHQJyIKIwz6RERhhEGfiCiM\nMOgTEYURBn0iojDCoE9EFEYY9ImIwgiDPhFRGGHQJyIKIwz6RERhhEGfiCiMMOgTEYURBn0iojDC\noE9EFEY6TdAvLS0NdBHIDdZNcGP9BDep60f2oH/q1Cncd999SEpKwpgxY1BQUAAAqK+vR3p6OoYP\nH46MjAxcuHChQ6/D/3GDF+smuLF+glvIBf1u3brhjTfewKFDh7Br1y6sX78e3333HfLy8pCeno4j\nR45g2rRpyMvLk7soRERhT/agr9FokJKSAgDo1asXRo0ahTNnzqCwsBB6vR4AoNfr8emnn8pdFCKi\nsKcSPG2bLrGqqiqkpaXh4MGDuP3229HQ0AAAEAQBffv2tf/bXjiVSqmiERF1Gp7CeoRShbhy5Qoe\nfvhh5Ofno3fv3k6PqVQq0QCv4PcREVFYUCR75/r163j44YfxxBNPICsrCwCgVqthNpsBADU1NYiP\nj1eiKEREYU32oC8IAp5++mmMHj0av/nNb+zHp0+fDqPRCAAwGo32LwMiIpKP7H36//rXv5Camoqx\nY8fau3DWrFmDyZMnIzs7GydPnkRCQgK2bNmCmJgYOYtCRERCiDObzcLMmTOFu+66S5g4caJwzz33\nCJ988kmgixW2unTpIqSkpAhjxowRHnnkEaGxsVEQBNaTUqR4/7ds2SKMHj1a6NKli7B371778W++\n+UZISUkRUlJShOTkZGHTpk2K/E3hylaXycnJwkMPPSRcvnxZmucN9JdORwiCgKysLGi1WlRWVmLP\nnj3YtGkTTp8+Heiiha0ePXqgoqICBw4cwG233YZ3330XAFhPCpHi/U9OTsYnn3yC1NRUl+N79+5F\nRUUFSkpKMH/+fLS2tsr+N4UrW13u378f0dHR+MMf/iDJ8yqWvSOHHTt2oHv37pg7d6792O23344F\nCxYEsFRkc++992L//v2spwDx9/0fOXKk6PGoqCj7701NTejTpw+6du0qbaFJ1D333IN9+/ZJ8lwh\n3dI/dOgQJkyYEOhikAiLxYLPPvsMY8eOZT0FgFzvf3l5OZKSkpCUlITXX39dkuckz1pbW1FSUoIx\nY8ZI8nwhHfTb5vYvWLAAKSkpmDx5coBKRE1NTRg/fjwmTZqEO+64A7NmzXI5h/UkH7nf/8mTJ+PQ\noUP49ttvsXjxYly8eFGKYpMIW10OGDAAp06dwrPPPivNE0syMhAgX3zxhZCWluZ07Pz580JCQkJg\nCkRCr169XI6xnpTjz/v/5JNPCikpKUJmZqbTOVqt1mkgt637779f2LNnT8cLTaJsddnY2Cjce++9\nwscffyzJ84Z0S//+++9Hc3OzfbAKAK5evRrAEpEY1lNgeXv/33vvPVRUVKCoqMjlWsEho7uqqgoW\niwUAcOLECRw9ehTDhg2TseQEWMdSCgoKkJOTI8kqBSEd9AHg008/hclkwl133YUpU6bgySefxGuv\nvRboYoUtd+slsZ6UIcX7/8knn2DIkCHYtWsXMjMz8V//9V8AgJ07dyIlJQXjx4/HI488gj/+8Y+I\njo6W7W8Jd451mZKSgqFDh2LLli0df15Biq8OIiIKCSHf0iciIt8x6BMRhREGfSKiMMKgT0QURhj0\niXzw5ptvoqmpyf7vzMxMXLp0KYAlIvIPs3eIbrJ9FMTSHu+8807s2bMHcXFxSheLSFJs6VNYq6qq\nwogRI6DX65GcnIynn34akyZNwpgxY2AwGAAABQUFqK6uxn333Ydp06YBABISElBfX4+qqiqMGjUK\nc+fOxZgxY6DT6dDc3AwA2L17N8aOHYvx48djyZIlSE5ODtSfSWTHoE9h78cff8T8+fNx8OBB/P73\nv8fu3buxb98+mEwmHDx4EIsWLcLAgQNRWlqKL774AoDz3cCPP/6IBQsW4ODBg4iJicHf//53AMBT\nTz2FP/3pT6ioqEBERITbiVNESmLQp7B3xx132Bcf27x5MyZOnIgJEybg0KFDOHz4sNfr77zzTowd\nOxYAMHHiRFRVVeHixYu4cuUKpkyZAgB47LHHJJlCT9RRIb2ePpEUevbsCQA4fvw4fv/732PPnj3o\n06cPnnrqKXtXjSfdu3e3/961a1enAV8bBnwKFmzpE9106dIl9OzZE9HR0aitrcVnn31mf6x3797t\nytbp06cPevfujfLycgDApk2bJC8vkT/Y0qewZ+trHzduHMaPH4+RI0diyJAh+M///E/7OXPnzsWD\nDz6IQYMG2fv1217f9t8bN27EnDlz0KVLF6SlpaFPnz4y/yVE3jFlk0gmV69etXcd5eXloba2Fm+8\n8UaAS0Xhji19IpkUFxdjzZo1sFgsSEhIwPvvvx/oIhGxpU9EFE44kEtEFEYY9ImIwgiDPhFRGGHQ\nJyIKIwz6RERhhEGfiCiM/H9bXrSmTgFG1wAAAABJRU5ErkJggg==\n" } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Question 1: What is the difference in average rating between G and R movies?\n", "\n", "Answer: $b_0 + b_3 - b_0 = b3$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "lm1.summary()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: score R-squared: 0.020
Model: OLS Adj. R-squared: -0.002
Method: Least Squares F-statistic: 0.9182
Date: Sat, 06 Apr 2013 Prob (F-statistic): 0.434
Time: 18:06:24 Log-Likelihood: -569.90
No. Observations: 140 AIC: 1148.
Df Residuals: 136 BIC: 1160.
Df Model: 3
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [95.0% Conf. Int.]
Intercept 67.6500 7.193 9.405 0.000 53.425 81.875
rating[T.PG] -12.5929 7.849 -1.604 0.111 -28.114 2.928
rating[T.PG-13] -11.8146 7.411 -1.594 0.113 -26.471 2.842
rating[T.R] -12.0200 7.476 -1.608 0.110 -26.803 2.763
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 3.313 Durbin-Watson: 2.100
Prob(Omnibus): 0.191 Jarque-Bera (JB): 3.327
Skew: 0.367 Prob(JB): 0.189
Kurtosis: 2.821 Cond. No. 14.0
" ], "output_type": "pyout", "prompt_number": 10, "text": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: score R-squared: 0.020\n", "Model: OLS Adj. R-squared: -0.002\n", "Method: Least Squares F-statistic: 0.9182\n", "Date: Sat, 06 Apr 2013 Prob (F-statistic): 0.434\n", "Time: 18:06:24 Log-Likelihood: -569.90\n", "No. Observations: 140 AIC: 1148.\n", "Df Residuals: 136 BIC: 1160.\n", "Df Model: 3 \n", "===================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "-----------------------------------------------------------------------------------\n", "Intercept 67.6500 7.193 9.405 0.000 53.425 81.875\n", "rating[T.PG] -12.5929 7.849 -1.604 0.111 -28.114 2.928\n", "rating[T.PG-13] -11.8146 7.411 -1.594 0.113 -26.471 2.842\n", "rating[T.R] -12.0200 7.476 -1.608 0.110 -26.803 2.763\n", "==============================================================================\n", "Omnibus: 3.313 Durbin-Watson: 2.100\n", "Prob(Omnibus): 0.191 Jarque-Bera (JB): 3.327\n", "Skew: 0.367 Prob(JB): 0.189\n", "Kurtosis: 2.821 Cond. No. 14.0\n", "==============================================================================\n", "\"\"\"" ] } ], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "lm1.conf_int()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
Intercept 53.424778 81.875222
rating[T.PG]-28.113847 2.928133
rating[T.PG-13]-26.471002 2.841772
rating[T.R]-26.803285 2.763285
\n", "
" ], "output_type": "pyout", "prompt_number": 11, "text": [ " 0 1\n", "Intercept 53.424778 81.875222\n", "rating[T.PG] -28.113847 2.928133\n", "rating[T.PG-13] -26.471002 2.841772\n", "rating[T.R] -26.803285 2.763285" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Question 2: What is the difference in average rating between PG-13 and R movies?\n", "\n", "Answer: $b_0 + b_2 - (b_0 + b_3) = b_2 - b_3$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could rewrite our model\n", "\n", "$$S_i = b_0 + b_1 1(Ra_i = \"G\") + b_2 1(Ra_i = \"PG\") + b_3 1(Ra_i = \"PG - 13\") + e_i$$\n", "\n", "Average values\n", "\n", "$b_0$ = average of R movies\n", "\n", "$b_0 + b_1$ = average of G movies\n", "\n", "$b_0 + b_2$ = average of PG movies\n", "\n", "$b_0 + b_3$ = average of PG-13 movies\n", "\n", "What is the difference in average rating between PG-13 and R movies?\n", "\n", "$b_0 + b_3 - b_0 = b3$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from patsy.contrasts import Treatment\n", "\n", "lm2 = ols('score ~ C(rating, Treatment(reference=\"R\"))', movies).fit()\n", "lm2.summary()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: score R-squared: 0.020
Model: OLS Adj. R-squared: -0.002
Method: Least Squares F-statistic: 0.9182
Date: Sat, 06 Apr 2013 Prob (F-statistic): 0.434
Time: 19:23:37 Log-Likelihood: -569.90
No. Observations: 140 AIC: 1148.
Df Residuals: 136 BIC: 1160.
Df Model: 3
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [95.0% Conf. Int.]
Intercept 55.6300 2.035 27.342 0.000 51.606 59.654
C(rating, Treatment(reference=\"R\"))[T.G] 12.0200 7.476 1.608 0.110 -2.763 26.803
C(rating, Treatment(reference=\"R\"))[T.PG] -0.5729 3.741 -0.153 0.879 -7.971 6.825
C(rating, Treatment(reference=\"R\"))[T.PG-13] 0.2054 2.706 0.076 0.940 -5.146 5.557
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 3.313 Durbin-Watson: 2.100
Prob(Omnibus): 0.191 Jarque-Bera (JB): 3.327
Skew: 0.367 Prob(JB): 0.189
Kurtosis: 2.821 Cond. No. 7.05
" ], "output_type": "pyout", "prompt_number": 18, "text": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: score R-squared: 0.020\n", "Model: OLS Adj. R-squared: -0.002\n", "Method: Least Squares F-statistic: 0.9182\n", "Date: Sat, 06 Apr 2013 Prob (F-statistic): 0.434\n", "Time: 19:23:37 Log-Likelihood: -569.90\n", "No. Observations: 140 AIC: 1148.\n", "Df Residuals: 136 BIC: 1160.\n", "Df Model: 3 \n", "================================================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "----------------------------------------------------------------------------------------------------------------\n", "Intercept 55.6300 2.035 27.342 0.000 51.606 59.654\n", "C(rating, Treatment(reference=\"R\"))[T.G] 12.0200 7.476 1.608 0.110 -2.763 26.803\n", "C(rating, Treatment(reference=\"R\"))[T.PG] -0.5729 3.741 -0.153 0.879 -7.971 6.825\n", "C(rating, Treatment(reference=\"R\"))[T.PG-13] 0.2054 2.706 0.076 0.940 -5.146 5.557\n", "==============================================================================\n", "Omnibus: 3.313 Durbin-Watson: 2.100\n", "Prob(Omnibus): 0.191 Jarque-Bera (JB): 3.327\n", "Skew: 0.367 Prob(JB): 0.189\n", "Kurtosis: 2.821 Cond. No. 7.05\n", "==============================================================================\n", "\"\"\"" ] } ], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "lm2.conf_int()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
Intercept 51.606500 59.653500
C(rating, Treatment(reference=\"R\"))[T.G] -2.763285 26.803285
C(rating, Treatment(reference=\"R\"))[T.PG] -7.971015 6.825300
C(rating, Treatment(reference=\"R\"))[T.PG-13] -5.146371 5.557140
\n", "
" ], "output_type": "pyout", "prompt_number": 20, "text": [ " 0 1\n", "Intercept 51.606500 59.653500\n", "C(rating, Treatment(reference=\"R\"))[T.G] -2.763285 26.803285\n", "C(rating, Treatment(reference=\"R\"))[T.PG] -7.971015 6.825300\n", "C(rating, Treatment(reference=\"R\"))[T.PG-13] -5.146371 5.557140" ] } ], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the model:\n", "\n", "$$S_i = b_0 + b_1 1(Ra_i = \"PG\") + b_2 1(Ra_i = \"PG - 13\") + b_3 1(Ra_i = \"R\") + e_i$$\n", "\n", "Question 3: Is there any difference in score between any of the movie ratings?\n", "\n", "To answer this we'll use ANOVA:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from statsmodels.stats.anova import anova_lm\n", "\n", "anova_lm(lm1)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
dfsum_sqmean_sqFPR(>F)
rating 3 570.123813 190.041271 0.918184 0.433975
Residual 136 28148.635044 206.975258 NaN NaN
\n", "
" ], "output_type": "pyout", "prompt_number": 21, "text": [ " df sum_sq mean_sq F PR(>F)\n", "rating 3 570.123813 190.041271 0.918184 0.433975\n", "Residual 136 28148.635044 206.975258 NaN NaN" ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our p-value is the probability $PR(>F)$. As before, in hypothesis-testing parlance, we'll take our null hypothesis, $H_0$, to be that there is no difference in score between any of the movie ratings, and the alternate hypothesis, $H_a$, to be the opposite. Here our p-value is greater than some threshold, e.g. 0.05, thus we can't reject the null hypothesis (i.e. we do not have significant evidence against the null hypothesis at the significance level of 0.05)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way is to use Tukey's HSD test:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from statsmodels.stats.multicomp import tukeyhsd\n", "import statsmodels.sandbox.stats.multicomp as multi\n", "\n", "mc = multi.MultiComparison(movies['score'], r.labels)\n", "res = mc.tukeyhsd() \n", "\n", "print res[0]\n", "print zip(list(r.levels), range(4))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Multiple Comparison of Means - Tukey HSD, FWER=0.05\n", "==============================================\n", "group1 group2 meandiff lower upper reject\n", "----------------------------------------------\n", " 0 1 -12.5929 -33.0089 7.8232 False \n", " 0 2 -11.8146 -31.0934 7.4642 False \n", " 0 3 -12.02 -31.4657 7.4257 False \n", " 1 2 0.7782 -8.6152 10.1717 False \n", " 1 3 0.5729 -9.1586 10.3043 False \n", " 2 3 -0.2054 -7.245 6.8342 False \n", "----------------------------------------------\n", "[('G', 0), ('PG', 1), ('PG-13', 2), ('R', 3)]\n" ] } ], "prompt_number": 22 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }