GenerelizedLinearModels/Generalized Linear Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What are Generalized Linear Models?\n",
    "\n",
    "GLMs were one of those concepts that didn't click until I got lucky to see it from a very particular angle. When I first encountered it, it looked arbitrary, random and unjustified. The subject finally made sense when I understood its original motivation: to generalize linear regression. Someone broke linear regression into pieces, generalized each of them and put them back together, yielding a machine with more nobs to turn.\n",
    "\n",
    "I'll now share that view with you. Hopefully you'll have the same reaction I did."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What are GLMs?\n",
    "\n",
    "As advertised, we'll start in the bubblegum land of linear regression. Here, we assume our dependent variable (a scalar $y$) is normally distributed *given* a linear combination of our independent variables (a vector $\\mathbf{x}$). This can be written as:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(y) & = \\mathcal{N}(y|\\mu,\\sigma^2) \\tag{1}\\\\\n",
    "\\mu & = \\mathbf{x}^T\\boldsymbol{\\beta} \\tag{2}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "I'll call (1) the 'distribution assumption' and (2) the 'expected-$y$ and $\\mathbf{x}$ relation assumption'. The latter has this name because $\\mu$ is the expectation of $y$ given $\\mathbf{x}$ (That is, $\\mu=\\mathbb{E}[y|\\mathbf{x}]$, but I'll use $\\mu$ going forward). \n",
    "\n",
    "A GLM generalizes both of these pieces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalize the distribution assumption\n",
    "\n",
    "Here, we generalize from:\n",
    "\n",
    "\"$y$ is **normally distributed** given $\\mathbf{x}$\"\n",
    "\n",
    "to\n",
    "\n",
    "\"$y$ has **any distribution within the exponential family** given $\\mathbf{x}$\"\n",
    "\n",
    "This is a good idea considering the exponential family is pretty damn awesome (link). To refresh, the *vector* $\\mathbf{y}$ has a distribution from the exponential family if its probability [1] is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(\\mathbf{y})=&\\frac{1}{Z(\\boldsymbol{\\theta})}h(\\mathbf{y})\\exp\\Big\\{{T(\\mathbf{y})\\cdot\\boldsymbol{\\theta}\\Big\\}}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "for some parameter vector $\\boldsymbol{\\theta}$. Our choices of the functions $h(\\cdot)$ and $T(\\cdot)$ will determine which distribution within the exponential family we're using. So we *could* choose them such we get the normal distribution. But a different choice might give us the binomial or poisson or something else. Though we can't choose $Z(\\boldsymbol{\\theta})$ - it's there to ensures $p(\\mathbf{y})$ is a properly normalized distribution. So our choice of $h(\\cdot)$ and $T(\\cdot)$ will force $Z(\\cdot)$ to be something.\n",
    "\n",
    "One thing to call out is that our dependent variable just went from necessarily a scalar to a vector (which could also be a scalar). So we've sneakily generalized to multi-output predictions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Before we do our next generalization...\n",
    "\n",
    "Let's pretend we've chosen $h(\\cdot)$ and $T(\\cdot)$ to recreate linear regression and write what we have so far:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(y) & = \\frac{1}{Z(\\theta)}h(y)\\exp\\Big\\{{T(y)\\cdot\\theta\\Big\\}} \\tag{1}\\\\\n",
    "\\mu & = \\mathbf{x}^T\\boldsymbol{\\beta} \\tag{2}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Notice an issue? How does our $\\mu$ connect to $\\boldsymbol{\\theta}$? We have to build that bridge. We can get that done with a function we'll call $\\Psi(\\cdot)$, which is defined to solve our problem: $\\theta = \\Psi(\\mu)$. This means if we select $h(\\cdot)$ and $T(\\cdot)$ to recreate linear regression, $\\Psi(\\cdot)$ will fall out as required. It is purely a result of the distribution we name and that distribution's connect between its $\\boldsymbol{\\theta}$-values and the resulting expectation of $\\boldsymbol{y}$[2].\n",
    "\n",
    "Now un-pretend and let's get back to generalizing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalize the expected-$y$ and $\\mathbf{x}$ relation assumption\n",
    "\n",
    "Now we generalize from:\n",
    "\n",
    "\"the scalar expectation of $y$ is a linear combination of $\\mathbf{x}$\"\n",
    "\n",
    "to\n",
    "\n",
    "\"**some simple function of our *vector* - **expectation of $y$ is a linear combination of $\\mathbf{x}$\"\n",
    "\n",
    "That 'simple' function is known as a link function (call it $g(\\cdot)$), so this may be written as $g(\\boldsymbol{\\mu})=\\mathbf{x}^T\\boldsymbol{\\beta}$. It's a generalization because $g(\\cdot)$ *could* be the identity function and the length of $\\boldsymbol{\\mu}$ *could* be one, landing us back in bubblegum land. But it could not and we could land elsewhere. Note that if $\\boldsymbol{\\mu}$ is a vector with length greater than 1, then $\\boldsymbol{\\beta}$ is a matrix.\n",
    "\n",
    "Also, the 'simple' requirement on $g(\\cdot)$ means its invertible, so we may write $\\boldsymbol{\\mu}=g^{-1}(\\mathbf{x}^T\\boldsymbol{\\beta})$\n",
    "\n",
    "## Put it back together.\n",
    "\n",
    "Now that we have both generalizations, let's put it all together:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(\\mathbf{y}) & = \\frac{1}{Z(\\boldsymbol{\\theta})}h(\\mathbf{y})\\exp\\Big\\{{T(\\mathbf{y})\\cdot\\boldsymbol{\\theta}\\Big\\}} \\tag{1}\\\\\n",
    "\\boldsymbol{\\theta} &  = \\Psi(\\boldsymbol{\\mu}) \\tag{'bridge'}\\\\\n",
    "\\boldsymbol{\\mu} & = g^{-1}(\\mathbf{x}^T\\boldsymbol{\\beta}) \\tag{2}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "And that's it. So if you are applying a generalized linear model, you have two choices to make:\n",
    "\n",
    "1. Given the independent variables, what distribution do you want $y$ to have? Normal, the poisson, the binomial, something else? You make that choice and that tells you $h(\\cdot)$ and $T(\\cdot)$ which forces $\\Psi(\\cdot)$ and $Z(\\cdot)$.\n",
    "2. What link function, $g(\\cdot)$, do you want? You have a few choices, but typically your answer to the previous question will guide what you should do here. Not all link functions can be combined with all exponential family distributions.\n",
    "\n",
    "With that, you have your model specification and you can start learning the parameters $\\boldsymbol{\\beta}$. Not too bad right?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An example\n",
    "\n",
    "To ground this concept, let's follow those steps in an easy example. Let's say we have a classification problem where $y$ is either 1 or 0. Let's answer those questions:\n",
    "\n",
    "1. Given an $\\mathbf{x}$, what should the distribution of $y$ be? A natural choice is the Bernoulli distribution. If we ask the internet, we discover this implies $h(y)=1$ and $T(y)=y$. These then tell us $\\Psi(\\mu)=\\log(\\frac{\\mu}{1-\\mu})=\\theta$ and $Z(\\theta) = \\exp(\\theta)+1$.\n",
    "2. The expectation of $y$ needs to be between 0 and 1 and $\\mathbf{x}^T\\boldsymbol{\\beta}$ can vary over the whole real line, so we should pick a function that maps from 0 and 1 to the real line. How about the logit function (which is coincidentally our $\\Psi()$ as well)? That is, $g(\\mu)=\\log(\\frac{\\mu}{1-\\mu})=\\mathbf{x}^T\\boldsymbol{\\beta}$.\n",
    "\n",
    "So if we sub in these settings and reduce things a bit, our model specification becomes:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(y) & = \\frac{1}{\\exp(\\mathbf{x}^T\\boldsymbol{\\beta})+1}\\exp\\Big\\{{y\\cdot\\mathbf{x}^T\\boldsymbol{\\beta}\\Big\\}}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "And whadda-ya-know, it's logistic regression! Easy!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A few things worth pointing out...\n",
    "\n",
    "If you want to survive in the wilderness of GLMs, you'll need to keep a few things in mind:\n",
    "\n",
    "1. I've presented the general form of GLMs, but it's easy to make choices of $h(\\cdot)$, $T(\\cdot)$ and $g(\\cdot)$ such that the optimization to determine $\\boldsymbol{\\beta}$ is nearly impossible. Because of this, most GLM software out there will restrict your choices, such that no matter what you pick, it'll be able to optimize that problem.\n",
    "2. At the same time, there are additional generalizations that software will offer that I haven't. They may offer a parameter that allows you to smoothly vary $h(\\cdot)$ or $T(\\cdot)$  or $\\Psi(\\cdot)$. Sometimes a weighted version will be offered. I'm sure there are others I'm not aware of.\n",
    "3. $\\Psi(\\cdot)$ can be pretty weird. Sometime it maps from a single input to multiple outputs. Also, it often holds your exogenously decided parameters - parameters that impact that distribution of $y$ but aren't learned from the data.\n",
    "\n",
    "These two points can make the specific form you see look different from the above. In such a circumstance, we should think carefully about the specific implementation offer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Footnotes\n",
    "\n",
    "[1] 'Probability' isn't correct here. I should say probability density or probability mass, but I'm trying to keep things simple (and wrong I guess?).\n",
    "\n",
    "[2] It's not necessary for my dicussion, but the following is good to know in general. Since picking a distribution within the exponential family fixes $\\Psi(\\cdot)$ and this function relates $\\boldsymbol{\\mu}$ (the expected-$\\mathbf{y}$) and the parameters $\\boldsymbol{\\theta}$, then if you know $\\boldsymbol{\\mu}$, then you know the parameters as well. Because of this, $\\boldsymbol{\\mu}$ is often called the mean *parameters* of the distribution. In the same thread, $\\boldsymbol{\\theta}$ are called the canonical parameters.\n",
    "\n",
    "### Sources\n",
    "\n",
    "[1] Kevin Murphy's book Machine Learning: A Probabilistic Perspective - I *really* love that book.\n",
    "\n",
    "[2] I decided to write this after reading the documentation from the Statsmodels package (link), so I should give them so love.\n"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}