DegreesOfFreedom.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "MathJax.Hub.Config({\n",
       "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
       "});\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "MathJax.Hub.Config({\n",
    "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
    "});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.stats as smstats\n",
    "\n",
    "import matplotlib\n",
    "#matplotlib.use('nbagg')\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In statistics, degrees of freedom are parameters associated with the chi-squared and associated distributions. We need to spcify one such parameter for the chi-squared and t distributions, whereas the F distribution requires two. The t distribution arises from the ratio of a standard normal variable and a chi-squared variable divided by its degrees of freedom, so it inherits the underlying degrees of freedom parameter. The F distribution is a ratio of two chi-squared variables divided by their respective degrees of freedom, and so needs both parameters for its  specification. The interesting question is where the degrees of freedom come from. In order to answer this question, we will have to look at idempotent quadratic forms. \n",
    "\n",
    "Given an $n\\times1$ vector $\\mathbf{x}$, a quadratic form is a scalar value defined as\n",
    "\\begin{equation}\n",
    "\\label{eq-qform}\n",
    "q = \\sum_{i,j} a_{ij} x_i x_j = \\mathbf{x'Ax}\n",
    "\\end{equation}\n",
    "where $\\mathbf{A}$ is an $n\\times n$ symmetric matrix. Now, since $\\mathbf{A}$ is symmetric, its spectral decomposition implies that\n",
    "\\begin{equation}\n",
    "\\mathbf{A = C'\\Lambda C}\n",
    "\\end{equation}\n",
    "where $\\mathbf{\\Lambda}$ is a diagonal matrix containing the eigenvalues of $\\mathbf{A}$ and the columns of $\\mathbf{C}$ are the eigenvectors of $\\mathbf{A}$, and $\\mathbf{C'C = CC' = I}$.\n",
    "\n",
    "Plugging this into equation (\\ref{eq-qform}), we get\n",
    "$$\n",
    "q = \\mathbf{x'C'\\Lambda Cx} = \\mathbf{y'\\Lambda y} = \\sum_i \\lambda_i y_i^2\n",
    "$$\n",
    "where we have defined $\\mathbf{y=Cx}$ and $\\lambda_i$ are the eigenvalues of $\\mathbf{A}$.\n",
    "\n",
    "An idempotent matrix is one whose square is the matrix itself. We have\n",
    "$$\n",
    "\\mathbf{A^2 = A A = C'\\Lambda C C' \\Lambda C = C' \\Lambda^2 C }\n",
    "$$\n",
    "\n",
    "But if $\\mathbf{A}$ is idempotent, $\\mathbf{A^2 = A}$, which implies that $\\mathbf{\\Lambda^2 = \\Lambda}$ or $\\lambda_i^2 = \\lambda_i$ for all $i$. This can only happen if every $\\lambda_i$ is either 0 or 1. This means that\n",
    "$$\n",
    "q = \\sum_{k=1}^K y_k^2\n",
    "$$\n",
    "where $k$ ranges over the non-zero eigenvalues of $\\mathbf{A}$. Note that the total number of non-zero eigenvalues $K$ is also the sum of all the eigenvalues of $\\mathbf{A}$, which is the trace of $\\mathbf{A}$.\n",
    "\n",
    "**Sample standard deviation**\n",
    "\n",
    "The sample standard deviaiton involves the sum of squared deviations from the mean. A vector of deviations from the mean can be expressed in matrix form. To do this, we start with a column vector of ones\n",
    "$$\n",
    "\\mathbf{i} =\n",
    "\\begin{pmatrix}\n",
    "1\\\\\n",
    "1\\\\\n",
    "\\vdots\\\\\n",
    "1\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "This means that $\\sum_i x_i = \\mathbf{i'x}$ and $\\bar{x}=\\frac{1}{n} \\mathbf{i'x}$ and $\\frac{1}{n} \\mathbf{ii'x} $ produces a column vector all of whose entries are $\\bar{x}$. From this, we get\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "x_1-\\bar{x}\\\\\n",
    "x_2-\\bar{x}\\\\\n",
    "\\vdots\\\\\n",
    "x_n-\\bar{x}\n",
    "\\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "x_1\\\\\n",
    "x_2\\\\\n",
    "\\vdots\\\\\n",
    "x_n\n",
    "\\end{pmatrix} -\n",
    "\\begin{pmatrix}\n",
    "\\bar{x}\\\\\n",
    "\\bar{x}\\\\\n",
    "\\vdots\\\\\n",
    "\\bar{x}\n",
    "\\end{pmatrix} =\n",
    "\\left(\\mathbf{I}-\\frac{1}{n}\\mathbf{ii'} \\right)\\;\\mathbf{x} \\equiv \\mathbf{M x}\n",
    "$$\n",
    "\n",
    "$\\mathbf{M}$ is symmetric and idempotent\n",
    "$$\n",
    "\\mathbf{M^2} = \\left(\\mathbf{I}-\\frac{1}{n}\\mathbf{ii'} \\right)\\left(\\mathbf{I}-\\frac{1}{n}\\mathbf{ii'} \\right) =\n",
    "\\mathbf{I} - \\frac{2}{n}\\mathbf{ii'} + \\frac{1}{n^2}\\mathbf{ii'\\; ii'} = \\mathbf{I} - \\frac{1}{n}\\mathbf{ii'}\n",
    "$$\n",
    "since $\\mathbf{ii'\\;ii'}=n\\;\\mathbf{ii'}$.\n",
    "\n",
    "This implies that we can write the sum of squared deviations from the mean as\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i=1}^n (x_i-\\bar{x})^2 = \\mathbf{x'Mx}\n",
    "\\end{equation}\n",
    "\n",
    "Now, if $\\frac{1}{\\sigma}\\mathbf{x}$ is a vector of normal random variables, so is $\\mathbf{y=\\frac{1}{\\sigma}Cx}$ where $\\mathbf{C}$ is the eigenvector matrix for $\\mathbf{M}$. Thus\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^n \\left(\\frac{x_i-\\bar{x}}{\\sigma}\\right)^2 = \\sum_{k=1}^K y_k^2 \\sim \\chi^2_K\n",
    "$$\n",
    "\n",
    "since a sum of squared normal random variables has a chi-squared distribution. All that remains is for us to determine $K$, which equals the trace of $\\mathbf{M}$. Now,\n",
    "$$\n",
    "K = \\text{trace}(\\mathbf{M}) = n(1-\\frac{1}{n}) = n-1\n",
    "$$\n",
    "So, we conclude that\n",
    "$$\n",
    "\\sum_{i=1}^n \\left(\\frac{x_i-\\bar{x}}{\\sigma}\\right)^2 \\sim \\chi^2_{n-1}\n",
    "$$\n",
    "\n",
    "**Regression**\n",
    "\n",
    "In multiple regression, we have a model of the form\n",
    "\\begin{equation}\n",
    "\\mathbf{y = X \\boldsymbol{\\beta + \\epsilon}} \\quad \\boldsymbol{\\epsilon} \\sim N(\\mathbf{0},\\sigma\\mathbf{I})\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathbf{y}$ is an $n\\times 1$ vector, $X$ is an $n \\times k$ design matrix, $\\boldsymbol{\\beta}$ is the vector of $k$ coefficients, and $\\boldsymbol{\\epsilon}$ is the $n \\times 1$ vector of standard normal random errors.\n",
    "\n",
    "We derive an estimate $\\mathbf{b}$ of $\\boldsymbol{\\beta}$ from a sample of observations for which\n",
    "\\begin{equation}\n",
    "\\label{eq-Reg}\n",
    "\\mathbf{y = Xb + e}\n",
    "\\end{equation}\n",
    "by minimizing the squared sum of errors $\\mathbf{e'e}$. The estimate is given by\n",
    "\\begin{equation}\n",
    "\\label{eq-b}\n",
    "\\mathbf{b = \\left(X'X\\right)^{-1} X'y}\n",
    "\\end{equation}\n",
    "\n",
    "Plugging equation (\\ref{eq-b}) into equation (\\ref{eq-Reg}), we get\n",
    "\\begin{equation}\n",
    "\\label{eq-H}\n",
    "\\mathbf{y = X \\left(X'X\\right)^{-1} X'y + e \\equiv Hy + e }\n",
    "\\end{equation}\n",
    "\n",
    "The predicted values $\\mathbf{Hy}$ are often denoted $\\mathbf{\\hat{y}}$, which explains the name \"hat matrix\" for $\\mathbf{H}$.\n",
    "\n",
    "The hat matrix $\\mathbf{H}$ is a projection matrix, since it projects $\\mathbf{y}$ on to the column space of $\\mathbf{X}$. From equation (\\ref{eq-H}), we have\n",
    "\n",
    "$$\n",
    "\\mathbf{e = y - Hy = (I-H)y = (I-H)(X\\boldsymbol{\\beta}+\\boldsymbol{\\epsilon}) = (I-H)\\boldsymbol{\\epsilon}}\n",
    "$$\n",
    "\n",
    "where the last equality follows from the fact that $\\mathbf{(I-H)X\\boldsymbol{\\beta} = \\boldsymbol{0}}$ since $\\mathbf{HX = X}$.\n",
    "\n",
    "It is easy to show that $\\mathbf{H}$ is symmetric and idempotent, which means that so is $\\mathbf{I-H}$. This implies that\n",
    "\n",
    "$$\n",
    "\\frac{\\mathbf{e'e}}{\\sigma^2} = \\frac{\\boldsymbol{\\epsilon}'}{\\sigma} (\\mathbf{I-H}) \\frac{\\boldsymbol{\\epsilon}}{\\sigma}\n",
    "$$\n",
    "\n",
    "is going to have a chi-squared distribution, since $\\boldsymbol{\\epsilon}/\\sigma$ is a vector of standard normal variates. The degrees of freedom for this chi-squared distribution is given by\n",
    "\n",
    "$$\n",
    "\\text{trace}(\\mathbf{I}_{n\\times n}-\\mathbf{H}) = \\text{trace}(\\mathbf{I}_{n\\times n}) - \\text{trace}\\left(\\mathbf{X(X'X)^{-1}X'}\\right) = n - \\text{trace}\\left(\\mathbf{X'X(X'X)^{-1}}\\right) = n - \\text{trace}\\left(\\mathbf{I}_{k\\times k}\\right) = n-k\n",
    "$$\n",
    "\n",
    "Thus, we have established that, in multiple regression\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\mathbf{e'e}}{\\sigma^2} \\sim \\chi^2_{n-k}\n",
    "\\end{equation}\n",
    "\n",
    "This leads to the definition of the **standard error of estimate** $s$:\n",
    "\\begin{equation}\n",
    "s^2 = \\frac{\\mathbf{e'e}}{n-k} = \\frac{\\sum_i e_i^2}{n-k}\n",
    "\\end{equation}\n",
    "\n",
    "so that $(n-k)s^2/\\sigma^2$ has the chi-squared distribution with $n-k$ degrees of freedom. The hypothesis tests  for regression coefficients all follow from this.\n",
    "\n"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}