08-02 statsmodels.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the statsmodels.api uses numpy array notation\n",
    "# statsmodels.formula.api use formula notation (similar to R's formula notation)\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "bmx = pd.read_sas('BMX_I.XPT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "demo = pd.read_sas('DEMO_I.XPT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged = bmx.merge(right = demo, on = 'SEQN')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "babies = merged.loc[merged.RIDAGEMN <= 12,['RIDAGEMN','RIAGENDR','BMXWT']]  \n",
    "# subset to select those with age in month 12 and under, select age in months, gender, weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "babies.columns = ['age', 'gender', 'weight']  #rename columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      age  gender  weight\n",
      "32    2.0     2.0     5.2\n",
      "63    7.0     2.0     8.2\n",
      "105  10.0     2.0     7.3\n",
      "115   5.0     2.0     7.2\n",
      "121   5.0     1.0     7.5\n"
     ]
    }
   ],
   "source": [
    "print(babies.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 399 entries, 32 to 9533\n",
      "Data columns (total 3 columns):\n",
      "age       399 non-null float64\n",
      "gender    399 non-null float64\n",
      "weight    399 non-null float64\n",
      "dtypes: float64(3)\n",
      "memory usage: 12.5 KB\n"
     ]
    }
   ],
   "source": [
    "babies.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## formula notation\n",
    "use statsmodels.formula.api (often imported as smf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = smf.ols('weight ~ age + gender', data = babies)\n",
    "# gender is currently coded 1 is boys, and 2 is girls\n",
    "# the current ols call will treat gender as a numeric variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<statsmodels.regression.linear_model.RegressionResultsWrapper object at 0x0000019C9C97BD30>\n"
     ]
    }
   ],
   "source": [
    "print(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 weight   R-squared:                       0.707\n",
      "Model:                            OLS   Adj. R-squared:                  0.705\n",
      "Method:                 Least Squares   F-statistic:                     477.0\n",
      "Date:                Wed, 23 May 2018   Prob (F-statistic):          3.47e-106\n",
      "Time:                        14:21:11   Log-Likelihood:                -582.83\n",
      "No. Observations:                 399   AIC:                             1172.\n",
      "Df Residuals:                     396   BIC:                             1184.\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept      6.4509      0.188     34.386      0.000       6.082       6.820\n",
      "age            0.4363      0.014     30.102      0.000       0.408       0.465\n",
      "gender        -0.7327      0.105     -6.987      0.000      -0.939      -0.527\n",
      "==============================================================================\n",
      "Omnibus:                       30.003   Durbin-Watson:                   1.873\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               48.571\n",
      "Skew:                           0.505   Prob(JB):                     2.84e-11\n",
      "Kurtosis:                       4.379   Cond. No.                         28.1\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(results.summary())\n",
    "# the intercept here is artificially too big because boys are 1, and girls are 2\n",
    "# so a baby boy will subtract 1 times 0.7327, and a girl will subtract 2 times 0.7327\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "model2 = smf.ols('weight ~ age + C(gender)', data = babies)  \n",
    "# the C() function makes categories\n",
    "# in our case, it doesn't make a difference because boys are 1 and girls are 2. (There are only 2 categories)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2 = model2.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 weight   R-squared:                       0.707\n",
      "Model:                            OLS   Adj. R-squared:                  0.705\n",
      "Method:                 Least Squares   F-statistic:                     477.0\n",
      "Date:                Wed, 23 May 2018   Prob (F-statistic):          3.47e-106\n",
      "Time:                        14:21:45   Log-Likelihood:                -582.83\n",
      "No. Observations:                 399   AIC:                             1172.\n",
      "Df Residuals:                     396   BIC:                             1184.\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "====================================================================================\n",
      "                       coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------------\n",
      "Intercept            5.7181      0.113     50.484      0.000       5.495       5.941\n",
      "C(gender)[T.2.0]    -0.7327      0.105     -6.987      0.000      -0.939      -0.527\n",
      "age                  0.4363      0.014     30.102      0.000       0.408       0.465\n",
      "==============================================================================\n",
      "Omnibus:                       30.003   Durbin-Watson:                   1.873\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               48.571\n",
      "Skew:                           0.505   Prob(JB):                     2.84e-11\n",
      "Kurtosis:                       4.379   Cond. No.                         17.6\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(results2.summary())\n",
    "\n",
    "# the itercept here is lower...\n",
    "# boys are coded as 0 and will not subtract anything from the intercept\n",
    "# girls are coded as 1 and will subtract 0.7327 from the intercerpt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## numpy array notation\n",
    "similar to sklearn's notation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create arrays for the X and for the y\n",
    "X = babies.loc[:,['age','gender']]  # the feature matrix is a capital X\n",
    "y = babies.weight  # response vector is lowercase y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = sm.add_constant(X)  # need to add a constant for the intercept term.\n",
    "# because we are using the numpy notation, we use sm rather than smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "model3 = sm.OLS(y, X)  # OLS is capitalized in the numpy notation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "results3 = model3.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 weight   R-squared:                       0.707\n",
      "Model:                            OLS   Adj. R-squared:                  0.705\n",
      "Method:                 Least Squares   F-statistic:                     477.0\n",
      "Date:                Wed, 23 May 2018   Prob (F-statistic):          3.47e-106\n",
      "Time:                        14:25:35   Log-Likelihood:                -582.83\n",
      "No. Observations:                 399   AIC:                             1172.\n",
      "Df Residuals:                     396   BIC:                             1184.\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          6.4509      0.188     34.386      0.000       6.082       6.820\n",
      "age            0.4363      0.014     30.102      0.000       0.408       0.465\n",
      "gender        -0.7327      0.105     -6.987      0.000      -0.939      -0.527\n",
      "==============================================================================\n",
      "Omnibus:                       30.003   Durbin-Watson:                   1.873\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               48.571\n",
      "Skew:                           0.505   Prob(JB):                     2.84e-11\n",
      "Kurtosis:                       4.379   Cond. No.                         28.1\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(results3.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to do... categorical variables in numpy notation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## the ANOVA table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = sm.stats.anova_lm(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             df      sum_sq     mean_sq           F         PR(>F)\n",
      "age         1.0  991.502398  991.502398  905.132967  2.449632e-104\n",
      "gender      1.0   53.481241   53.481241   48.822508   1.193519e-11\n",
      "Residual  396.0  433.787039    1.095422         NaN            NaN\n"
     ]
    }
   ],
   "source": [
    "print(table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              df      sum_sq     mean_sq           F         PR(>F)\n",
      "C(gender)    1.0   52.362444   52.362444   47.801169   1.898518e-11\n",
      "age          1.0  992.621194  992.621194  906.154306  2.096787e-104\n",
      "Residual   396.0  433.787039    1.095422         NaN            NaN\n"
     ]
    }
   ],
   "source": [
    "table2 = sm.stats.anova_lm(results2)\n",
    "print(table2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'OLS' object has no attribute 'model'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-39-e78b3ca0d61d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m## todo: fix this\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mtable3\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0manova_lm\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmodel3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtable3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\statsmodels\\stats\\anova.py\u001b[0m in \u001b[0;36manova_lm\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m    324\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    325\u001b[0m         \u001b[0mmodel\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 326\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0manova_single\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    327\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    328\u001b[0m     \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\statsmodels\\stats\\anova.py\u001b[0m in \u001b[0;36manova_single\u001b[1;34m(model, **kwargs)\u001b[0m\n\u001b[0;32m     61\u001b[0m         \u001b[0mrobust\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrobust\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlower\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     62\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 63\u001b[1;33m     \u001b[0mendog\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mendog\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     64\u001b[0m     \u001b[0mexog\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mexog\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     65\u001b[0m     \u001b[0mnobs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mexog\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: 'OLS' object has no attribute 'model'"
     ]
    }
   ],
   "source": [
    "## todo: fix this\n",
    "table3 = sm.stats.anova_lm(model3)\n",
    "print(table3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_new = {'age':4, 'gender':1}  # A dictionary to define the data of a 4 month old boy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'age': 4, 'gender': 1}\n"
     ]
    }
   ],
   "source": [
    "print(x_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    7.463456\n",
       "dtype: float64"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.predict(x_new)  # using the model where gender is coded numerically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    4.985379\n",
       "dtype: float64"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_new = {'age':0, 'gender':2} # a 0 month old girl\n",
    "results.predict(x_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "32       5.858045\n",
       "63       8.039710\n",
       "105      9.348708\n",
       "115      7.167044\n",
       "121      7.899789\n",
       "181     10.954119\n",
       "197      7.167044\n",
       "200      8.772454\n",
       "213      8.476043\n",
       "214      4.985379\n",
       "295      4.985379\n",
       "299      6.154457\n",
       "311      9.785041\n",
       "315      8.476043\n",
       "318      5.858045\n",
       "330      8.912376\n",
       "372     10.221374\n",
       "383      7.899789\n",
       "437      9.348708\n",
       "465     10.954119\n",
       "482      5.858045\n",
       "486      6.730711\n",
       "491      9.208787\n",
       "492      9.208787\n",
       "499      6.590790\n",
       "504      9.785041\n",
       "544      7.603377\n",
       "549     10.081453\n",
       "570     10.517786\n",
       "589      5.421712\n",
       "          ...    \n",
       "8749     8.912376\n",
       "8766     7.603377\n",
       "8790     6.730711\n",
       "8793     8.476043\n",
       "8805     9.645120\n",
       "8820     7.027123\n",
       "8857     9.785041\n",
       "8893    10.954119\n",
       "8930     7.603377\n",
       "8935     7.167044\n",
       "9012     6.154457\n",
       "9050     7.027123\n",
       "9056    10.517786\n",
       "9065     7.603377\n",
       "9079     7.899789\n",
       "9086    10.081453\n",
       "9156     5.718124\n",
       "9187     7.899789\n",
       "9229     9.785041\n",
       "9230     7.463456\n",
       "9249     9.348708\n",
       "9259     6.294378\n",
       "9275     7.899789\n",
       "9293     7.027123\n",
       "9344     6.730711\n",
       "9354     7.603377\n",
       "9388     8.772454\n",
       "9480     6.590790\n",
       "9529    10.954119\n",
       "9533    10.954119\n",
       "Length: 399, dtype: float64"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.predict(X)# running the model on the entire X dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}