chapter-03/3.2.3-multiple-regression-from-simple-multivariate-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.2.3 Multiple Regression From Simple Univariate Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we have a *univariate* (p = 1) model with no intercept (3.23):\n",
    "$$Y=X\\beta+\\varepsilon$$\n",
    "\n",
    "The least squares estimate and residuals are (3.24):\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\hat{\\beta} = \\cfrac{\\sum_1^N {x_iy_i}}{\\sum_1^N {x_i^2}} \\\\\n",
    "r_i = y_i - x_i\\hat{\\beta}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "With the inner product:\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\hat{\\beta} = \\cfrac{\\langle \\mathbf{x}, \\mathbf{y} \\rangle}{\\langle \\mathbf{x}, \\mathbf{x}\\rangle}\\\\\n",
    "\\mathbf{r} = \\mathbf{y} - \\mathbf{x}\\hat{\\beta}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Suppose that the columns of the matrix **X** are orthogonal; that is $\\langle \\mathbf{x}_j, \\mathbf{x}_k \\rangle = 0$\n",
    "then it is easy to check that $\\hat{\\beta_j} = \\langle \\mathbf{x}_j, \\mathbf{y} \\rangle / \\langle \\mathbf{x}_j, \\mathbf{x}_j \\rangle$, i.e the inputs have no effect on each other's parameter estimates.\n",
    "\n",
    "Suppose next that we have an intercept and a single input x (3.27):\n",
    "$$\\hat{B}_1 = \\cfrac{\\langle \\mathbf{x} - \\overline{x}\\mathbf{1}, \\mathbf{y} \\rangle}{ \\langle \\mathbf{x} - \\overline{x}\\mathbf{1}, \\mathbf{x} - \\overline{x}\\mathbf{1} \\rangle}$$\n",
    "\n",
    "We can view the estimate as the result of two simple regression:\n",
    "\n",
    "1. Regress **x** on **1** to produce the residual $\\mathbf{z} = \\mathbf{x} - \\overline{x}\\mathbf{1}$\n",
    "\n",
    "2. Regress **y** on the residual **z** to give the coefficient $\\hat{\\beta}_1$.\n",
    "\n",
    "Regress **b** on **a** means $\\hat{\\gamma}=\\langle \\mathbf{a},\\mathbf{b} \\rangle / \\langle \\mathbf{a}, \\mathbf{a}\\rangle$ and the residual vector $\\mathbf{b} - \\hat{\\gamma}\\mathbf{a}$.\n",
    "\n",
    "This recipe generalizes to the case of *p* inputs, as shown in Algorithm 3.1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Algorithm 3.1 Regression by Successive Orthogonalization**\n",
    "1. $\\mathbf{z}_0 = \\mathbf{x}_0 = \\mathbf{1}$\n",
    "\n",
    "2. For $j = 1, 2, \\cdots, p$\n",
    "   \n",
    "   * Regress $\\mathbf{x}_j$ on $\\mathbf{z}_0,...,\\mathbf{z}_{j - 1}$ to produce $\\hat{\\gamma}_{lj}=\\langle \\mathbf{z}_l, \\mathbf{x}_j \\rangle / \\langle \\mathbf{z}_l,\\mathbf{z}_l \\rangle$ $l=0,\\cdots,j-1$, and residualt vector $\\mathbf{z}_j=\\mathbf{x}_j - \\sum_{k=0}^{j-1} \\hat{\\gamma}_{kj}\\mathbf{z}_k$\n",
    "\n",
    "3. Regress $\\mathbf{y}$ on the residual $\\mathbf{z}_p$ to give the estimate $\\hat{\\beta}_p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats, linalg\n",
    "\n",
    "df = pd.read_csv('../data/prostate/prostate.data', delimiter='\\t', index_col=0)\n",
    "mask_train = df.pop('train')\n",
    "df_y = df.pop('lpsa')\n",
    "df = df.apply(stats.zscore)\n",
    "\n",
    "def orthogonalize(X):\n",
    "    p = X.shape[1]\n",
    "    G = np.eye(p)\n",
    "    Z = X.copy()\n",
    "    for j in range(1, p): \n",
    "        for l in range(j):\n",
    "            G[l, j] = np.dot(Z[:, l], X[:, j]) / np.dot(Z[:, l], Z[:, l])\n",
    "        for k in range(j):\n",
    "            Z[:, j] -= G[k, j] * Z[:, k]\n",
    "    return Z, G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of this algorithm is (3.28):\n",
    "\n",
    "$$\\hat{\\beta}_p=\\cfrac{\\langle \\mathbf{z}_p, \\mathbf{y} \\rangle}{\\langle \\mathbf{z}_p,\\mathbf{z}_p \\rangle}$$\n",
    "\n",
    "If $\\mathbf{x}_p$ is highly correlated with some of the other $\\mathbf{x}_k$'s the residual vector $\\mathbf{x}_p$ will be close to zero, and from (3.28) the coefficient $\\hat{\\beta}_p$ will be unstable. \n",
    "\n",
    "From (3.28) we also obtain an alternative formula for the variance estimates, (3.29):\n",
    "\n",
    "$$Var(\\hat{\\beta}_p) = \\cfrac{\\sigma^2}{\\langle \\mathbf{z}_p, \\mathbf{z}_p \\rangle}=\\cfrac{\\sigma^2}{||\\mathbf{z}_p||^2}  $$\n",
    "\n",
    "On other words, the precision with which we can estimate $\\hat{\\beta}_p$ depends on the lengths of the residual vector $\\mathbf{z}_p$;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algorithm 3.1 is known as the *Gram–Schmidt* procedure for multiple regression. We can represent step 2 of Algorithm 3.1 in matrix form (3.30):\n",
    "\n",
    "$$\\mathbf{X}=\\mathbf{Z\\Gamma}$$\n",
    "\n",
    "where $\\mathbf{Z}$ has as columns the $z_j$ (in order), and $\\mathbf{\\Gamma}$ is the upper triangular matrix\n",
    "with entries $\\hat{\\gamma}_{kj}$. Introducing the diagonal matrix $\\mathbf{D}$ with $D_{jj}=||z_j||$, we get (3.31):\n",
    "\n",
    "$$\\mathbf{X}=\\mathbf{Z}\\mathbf{D}^{-1}\\mathbf{D}\\mathbf{\\Gamma}=\\mathbf{QR}$$\n",
    "\n",
    "the so-called QR decomposition of $\\mathbf{X}$. Here $\\mathbf{Q}$ is an N × (p +1) orthogonal\n",
    "matrix, $\\mathbf{Q}^T\\mathbf{Q} = \\mathbf{I}$, and **R** is a (p + 1) × (p + 1) upper triangular matrix.\n",
    "\n",
    "The least squares solution is given by:\n",
    "\n",
    "$$\n",
    "\\hat{\\beta}=\\mathbf{R}^{-1}\\mathbf{Q}^T\\mathbf{y}\n",
    "$$\n",
    "\n",
    "*Proof*:\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf{X}^T\\mathbf{y}=\\mathbf{X}^T\\mathbf{X}\\hat{\\beta}\\\\\n",
    "\\mathbf{R}^T\\mathbf{Q}^T\\mathbf{y}=\\mathbf{R}^T\\mathbf{Q}^T\\mathbf{Q}\\mathbf{R}\\hat{\\beta}\\\\\n",
    "\\mathbf{R}^T\\mathbf{Q}^T\\mathbf{y}=\\mathbf{R}^T\\mathbf{R}\\hat{\\beta}\\\\\n",
    "\\mathbf{Q}^T\\mathbf{y}=\\mathbf{R}\\hat{\\beta}\\\\\n",
    "\\end{equation}\n",
    "$$\n",
    "And the predicted training values:\n",
    "\n",
    "$$\n",
    "\\hat{\\mathbf{y}}=\\mathbf{QQ}^T\\mathbf{y}\n",
    "$$\n",
    "\n",
    "*Proof*:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{\\mathbf{y}}&=\\mathbf{X}\\hat{\\beta}\\\\\n",
    "&=\\mathbf{QR}\\mathbf{R}^{-1}\\mathbf{Q}^T\\mathbf{y}\\\\\n",
    "&=\\mathbf{QQ}^T\\mathbf{y}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can obtain from it not just $\\hat{\\beta}_p$, but also the entire multiple least squares fit.\n",
    "\n",
    "*Proof*:\n",
    "We can easily derive that:\n",
    "$$\n",
    "\\mathbf{R}\\hat{\\beta}=\\mathbf{Q}^T\\mathbf{y}\n",
    "$$\n",
    "\n",
    "which can be expanded into:\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{bmatrix}\n",
    "    R_{0 0} & R_{02}  & \\dots   & R_{0p} \\\\\n",
    "    0       & R_{11}  & \\dots   & R_{1p} \\\\\n",
    "    \\vdots  & \\vdots  & \\ddots  & \\vdots \\\\\n",
    "    0       & 0       & \\dots   & R_{pp}\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "    \\hat{\\beta_0} \\\\\n",
    "    \\hat{\\beta_1} \\\\\n",
    "    \\vdots        \\\\\n",
    "    \\hat{\\beta_p} \n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "    {Q_{0}}^T\\mathbf{y} \\\\\n",
    "    {Q_{1}}^T\\mathbf{y} \\\\\n",
    "    \\vdots        \\\\\n",
    "    {Q_{p}}^T\\mathbf{y} \n",
    "\\end{bmatrix}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Now by applying the backward substitution it is possible to obtain the entire multiple least squares fit. For example to find the $\\hat{\\beta}_p$:\n",
    "$$\n",
    "\\begin{equation}\n",
    "R_{pp}\\hat{\\beta}_p = {Q_{p}}^T\\mathbf{y}\\\\\n",
    "\\hat{\\beta}_p = \\cfrac{\\langle Q_p, \\mathbf{y} \\rangle}{R_{pp}}=\\cfrac{\\langle \\mathbf{z}_p, \\mathbf{y} \\rangle}{\\langle \\mathbf{z}_p,\\mathbf{z}_p \\rangle}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient:  [ 2.46493292  0.67601634  0.26169361 -0.14073374  0.20906052  0.30362332\n",
      " -0.28700184 -0.02119493  0.26557614]\n"
     ]
    }
   ],
   "source": [
    "def least_squares_qr(data_x, data_y):\n",
    "    X = np.c_[np.ones((len(data_x), 1)), data_x]\n",
    "    Z, G = orthogonalize(X)\n",
    "\n",
    "    D = linalg.norm(Z, axis=0)\n",
    "    Q = Z / D\n",
    "    R = np.diag(D) @ G\n",
    "    beta = linalg.solve_triangular(R, Q.T @ data_y)\n",
    "    return beta\n",
    "\n",
    "beta = least_squares_qr(df[mask_train == 'T'].as_matrix(), df_y[mask_train == 'T'].as_matrix())\n",
    "print (\"Coefficient: \", beta)"
   ]
  }
 ],
 "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
}