Introduction
This notes provides the necessary mathematical, statistical and computational details of Generalized Linear Model (GLM).
Three components of GLM
- Random Component: Response variable $Y$ with independent observations $(y_1,\cdots, y_n)$ having probability density or mass function for a distribution
- A linear predictor: For observation $i~~ (i = 1,\cdots, n)$, let $x_{ij}$ denote the value of explanatory variable $x_j~~ (j = 1,\cdots, p)$. Let $x_i$ = ($x_{i1}$, \cdots, $x_{ip}$). Usually, we set $x_{i1}$ = 1, so that it serves as the coefficient of the intercept term in a model.
- A link function: Relates parameter $\eta_i$ related to E($y_i$) to the explanatory variables $x_1$,$x_2$,$x_3$,\cdots,$x_p$ using a linear combination $\eta_i = \sum\beta_jX_{ij},~~~~~~j=1,2,3,\cdots,p,~~~~~~ i = 1,\cdots,n.$
First let us confine to two random components (Bernoulli and Poisson) from an exponential family of distribution and canonical link functions.
Exponential Families
Exponential family of distribution has likelihood of the form
$f(y|\theta,\phi) = \exp(\frac{y\theta – b(\theta)}{a(\phi)} + c(y,\phi))$ where $\theta$ is natural parameter and $\phi$ is scale parameter
Preliminary Results
Result 1:$E[\frac{\partial L}{\partial \theta }] = 0$
Result 2:$E[\frac{\partial^{2} l}{\partial \theta^{2} }] = -E[(\frac{\partial l}{\partial \theta })^{2}]$
Proofs of these two results are available at the end of this notes
From Result 1, it can be observed that for an exponential family of distribution
$E[Y] = \frac{\partial l}{\partial \theta} = b'(\theta)$
Also following is an immediate consequence of Result 2
$E[\frac{\partial^2 l}{\partial \theta^2}] = -\frac{b”(\theta)}{a(\phi)}$
Hence, $\frac{b”(\theta)}{a(\phi)} = \frac{V(Y)}{a(\phi)^2}$
$\Rightarrow V(Y) = a(\phi)b”(\theta)$
Canonical link
If $\eta= \theta$, then the link function $\eta$ is said to be canonical link
Example 1: Bernoulli Distribution
$Y_{1},Y_{2},Y_{3},\cdots,Y_{n} \sim \text{Bernoulli}(p_{i})$
$L(p_i|y) = \prod p_i^{y_i} (1-p_i)^{1-y_i}$ = $= \prod \exp\Big[y_i\log p_i + (1-y_i) \log (1-p_i)\Big]$
$= \exp\sum\Big[y_i\log p_i + (1-y_i) \log (1-p_i)\Big]$
$= \exp\sum\Big[y_i\log \frac{p_i}{1-p_i} + \log (1-p_i)\Big]$
Hence, $l(\theta |y) = \sum\Big[y_i\theta_i – \log (1+\exp~ \theta_i)\Big]$ where $\theta_i = \log \frac{p_i}{1-p_i}$
$\theta_i = \log \frac{p_i}{1-p_i}; b(\theta_i) = \log(1+\exp~\theta_i);~ \text{and}~ a(\phi) =1$
$\Rightarrow b'(\theta_i) = \frac{e^{\theta_i}}{1+e^{\theta_i}} = p_i$ and
$~b”(\theta_i) = \frac{e^{\theta_i}}{(1+e^{\theta_i})^2}=p_i(1-p_i)$
Example 2: Poisson Distribution
$Y_{1},Y_{2},Y_{3},\cdots,Y_{n} \sim \text{Poisson}(\lambda_{i})$
$l(\lambda_{i}|y) = \prod e^{-\lambda_{i}}\frac{\lambda_{i}^{y_{i}}}{y_{i}!} = \prod \exp[-\lambda_{i} \log(\frac{\lambda_{i}^{y_{i}}}{y_{i}!})]$
$= \exp(\sum^{n}_{i=1} ( -\lambda_{i} + y_{i}\log\lambda_{i} – \log y_{i}!))$
$= \exp(\sum^{n}_{i=1} (\frac{y_{i}~\log\lambda_{i}-\lambda_{i}}{1} – \log y_{i}!))$
$\Rightarrow l(\theta|y) = \sum^{n}_{i=1} (\frac{y_{i}~\log\lambda_{i}-\lambda_{i}}{1} – \log y_{i}!)$
$\theta_{i} = \log \lambda_{i}~;~~ b(\theta_{i}) = \lambda_{i} = e^{\theta_{i}}$
$b'(\theta{i}) = b”(\theta_{i}) = e^{\theta_{i}}$
$\Rightarrow b'(\theta_{i}) = b” (\theta_{i}) = \lambda_{i}$
$a(\phi) =$ 1
Remark
These two examples illustrate canonical link functions; logit link for Bernoulli and log link for Poisson
Parameter Estimation
Consider the PDF of an exponential form of distribution
$f(y|\theta,\phi) = \exp(\frac{y\theta – b(\theta)}{a(\phi)} + c(y,\phi))$
The log likelihood is
$l = \sum^{n}_{i=1} (\frac{y_{i}~\theta_{i} – b(\theta_{i})}{a(\phi_{i})} + c(y_{i},\phi_{i}))$
Then the $i^{th}$ component of the likelihood is
$l_i = \frac{y_i~\theta_{i} – b(\theta_{i})}{a(\phi_{i})} + c(y_{i},\phi_{i})$
$\frac{\partial l_{i}}{\partial \theta_{i}} = \frac{y_{i}~-~ b'(\theta_i)}{a(\phi)}$
$\Rightarrow \frac{\partial^{2} l}{\partial \theta_{i}^{2}} = – \frac{b”(\theta)}{a(\phi)}$ is independent of $y_{i}$
$\Rightarrow E[\frac{\partial^{2} l}{\partial \theta^{2} }] = – \frac{b”(\theta)}{a(\phi)}$
Hence,
$\frac{b”(\theta)}{a(\phi)} = \frac{V(y_{i})}{a(\phi)^2}$
$V(y_{i}) = a(\phi) b”(\theta)$
Now using chain relation
$\frac{\partial l_{i}}{\partial \beta_{j}} = \frac{\partial l_{i}}{\partial \theta_{i}} \frac{\partial \theta_{i}}{\partial \mu_{i}} \frac{\partial \mu_{i}}{\partial \eta_{i}}\frac{\partial \eta_{i}}{\partial \beta_{j}}$
$\mu_{i} = b'(\theta_{i})$
i.e ,$\frac{\partial \theta_{i}}{\partial \mu_{i}} = \frac{1}{b”(\theta_{i})}$
Also $\eta_{i} = \sum^{k}_{j=0}x_{ij}\beta_{j}$
$\frac{\partial \eta_{i}}{\partial \beta_{j}} = x_{ij}$
$\Rightarrow\frac{\partial l_{i}}{\partial \beta_{j}} = \frac{y_{i} – \mu_{i}}{a(\phi)}\frac{1}{b”(\theta_{i})}\frac{\partial \mu_{i}}{\partial \eta_{i}}x_{ij}$
=$\frac{y_{i} – \mu_{i}}{a(\phi)}\frac{a(\phi)}{V(y_{i})}\frac{\partial \mu_{i}}{\partial \eta_{j}}x_{ij}$
$\frac{\partial l_{i}}{\partial \beta_{j}} = \frac{y_{i} – \mu_{i}}{V(y_{i})}\frac{\partial \mu_{i}}{\partial \eta_{i}}x_{ij}$
i.e
$\frac{\partial l}{\partial \beta} = \sum^{n}_{i=1} \frac{y_{i} – \mu_{i}}{V(y_{i})}\frac{\partial \mu_{i}}{\partial \eta_{i}}x_{ij}$
The term in the sum is $j^{th}$ element of the matrix $X^{T}DV^{-1}(Y- \mu)$, where D is $\text{Diag}(\frac{\partial\mu_i}{\partial\eta_i})$ and V is $\text{Diag}(V(y_{i}))$
Now, $E[\frac{\partial^{2} l_{i}}{\partial \beta_{h}\partial \beta_{j}}]$
$=~E[\frac{\partial^{2} l_{i}}{\partial \beta_{h}}\frac{\partial^{2} l_{i}}{\partial \beta_{j}}]$
$=~E[\frac{y_{i} – \mu_{i}}{V(y_{i})}\frac{\partial \mu_{i}}{\partial \eta_{i}}x_{ih} \frac{y_{i} – \mu_{i}}{V(y_{i})}\frac{\partial \mu_{i}}{\partial \eta_{i}}x_{ij} ]$
$=~ \frac{x_{ih}x_{ij}}{V(y_{i})} (\frac{\partial \mu_{i}}{\partial \eta_{i} })^{2}$
i.e
$E[\frac{\partial^{2} l_{i}}{\partial \beta_{h}\partial \beta_{j} }] = \sum^{n}_{i=1} \frac{x_{ih}x_{ij}}{V(Y_{i})} (\frac{\partial \mu_{i}}{\partial \eta_{i} })^{2}$
This can be written in matrix form as $X_{P \times n}^{T} W_{n \times n} X_{n \times p}$ where $W = \text{Diag} ( \frac{(\frac{\partial \mu_{i}}{\partial \eta_{i}})^2}{V(y_{i})})$
In Exponential Family of distribution, other than normal distribution a closed form solution is not available to solve p equations arising from equating $\frac{\partial l}{\partial\beta}=0$. Instead, to obtain the maximum likelihood estimator numerically, we must resort to an iterative algorithm such as Newton-Raphson or Fisher scoring method.
That is, at the iteration step m, $\beta^{(m+1)} = \beta^{(m)}-[H^{-1}]^{(m)}q^{(m)}$ where H is the Hessian Matrix. The starting value (m = 0) is assumed to be known
$H = [\frac{\partial^{2}l_{i}}{\partial\beta_{h}\partial\beta_{j}}]_{p \times p}$ and $q = [\frac{\partial l}{\partial \beta}]_{p \times 1}$
In many cases instead of the observed information matrix [-H] , we can use Fisher expected information matrix F=E[-H]. Hence, the iterative formula becomes
$\beta^{(m+1)} = \beta^{(m)} + [F^{-1}]^{(m)}q^{(m)}$
The above procedure can be written as a step-by-step procedure
Example 1: Bernoulli Distribution
$\theta_i = \text{logit}(p_i)$; $b'(\theta_i)= V(y_i)$
$W=\text{Diag}\Big(V(y_i)\Big)_{n\times n}$
The algorithm is,
- Set $\beta_{p \times 1} = \beta_{0}$;
- Find $X_{n \times p} \beta_{p \times 1}$; $\exp(X\beta)$ and $1+\exp(X\beta)$
- Find $p_i = \frac{e^{X\beta}}{1+e^{X\beta}}$
- Find $u_{p \times 1} = X_{p \times n }^{T}(Y- \mu)_{n \times 1}$ and $W_{n \times n} = \text{Diag}(p_{i}(1-p_{i}))$
- Find $F_{p \times p} = X_{p \times n }^{T} W_{n \times n}X_{n \times p}$ and $F^{-1}$
Hence the iterative equation becomes $\beta_{p \times 1}^{(m+1)} = \beta_{p \times 1}^{(m)} – [F^{-1}]_{p \times p}^{(m)}u_{p \times 1}^{(m)}$
Example 2: Poisson Distribution
In the case of Poisson distribution,$\theta_{i} = \log(\lambda_{i}) = \log(\mu_{i})= \sum\beta_{j}x_{ij}$
Poisson model with log link
$\theta_{i} = \log(\lambda_{i}), \mu_{i} = b'(\theta_{i}) = \lambda_{i}$ $a(\phi)=1$
$\Rightarrow V(y_{i}) = b”(\theta_{i})a(\phi) = b”(\theta_{i}) = \lambda_{i}$
$u = X^{T}DV^{-1}(Y- \mu) = X^{T}(Y- \mu)$
$W = \text{Diag}(\frac{(b”(\theta_{i}))^2}{V(y_{i})}) = \text{Diag}(V(y_{i})) = \text{Diag}(\lambda_{i})$
The algorithm is,
- Set $\beta_{p \times 1} = \beta_{0}$
- Find $X_{n \times p} \beta_{p \times 1}$
- Find $\lambda_{n\times 1 } = e^{X\beta}$
- Find $u_{p \times 1} = X_{p \times n }^{T}(Y- \mu)_{n \times 1}$ and $W_{n \times n} = \text{Diag}(\lambda_{i})$
- Find $F_{p \times p} = X_{p \times n }^{T} W_{n \times n}X_{n \times p}$ and $F^{-1}_{p \times p}$
Hence the iterative equation becomes $\beta_{p \times 1}^{(m+1)} = \beta_{p \times 1}^{(m)} – [F^{-1}]_{p \times p}^{(m)}u_{p \times 1}^{(m)}$
Proof of Preliminary Results
Result 1:$E[\frac{\partial L}{\partial \theta }] = 0$
Proof:
Consider l = $\ln(f(x|\theta))$
$E[\frac{\partial l}{\partial \theta }] =\frac{f(x|\theta)} {f'(x|\theta)}$
$E[\frac{\partial l}{\partial \theta }] = \int \frac{f(x|\theta)} {f'(x|\theta)} f(x|\theta) \,dx = \int \frac{\partial }{\partial \theta} f(x|\theta) \,dx = \frac{\partial }{\partial \theta} \int f(x|\theta) \,dx$
$\Rightarrow E[\frac{\partial l}{\partial \theta }] = 0$
Result 2: $E[\frac{\partial^{2} l}{\partial \theta^{2} }] = -E[(\frac{\partial l}{\partial \theta })^{2}]$
Proof:
Now ,$E[\frac{\partial^{2} l}{\partial \theta^{2} }] = \frac{\partial }{\partial \theta} [\frac{\partial l}{\partial \theta } ]$
$E[\frac{\partial^{2} l}{\partial \theta^{2} }] = \int \frac{\partial }{\partial \theta}\frac{\partial l}{\partial \theta }f(x|\theta) \,dx$
Also from
$E[\frac{\partial l}{\partial \theta }] = 0$
$\Rightarrow$ $\frac{\partial }{\partial \theta} E[\frac{\partial l}{\partial \theta }] = 0$
$\Rightarrow$ $\frac{\partial }{\partial \theta} \int \frac{\partial }{\partial \theta} f(x|\theta) \,dx= 0$
$\Rightarrow$ $\int \frac{\partial }{\partial \theta}\frac{\partial l}{\partial \theta } f(x|\theta) \,dx= 0$
$\Rightarrow$ $\int \frac{\partial^{2} l}{\partial \theta^{2} } f(x|\theta) + \frac{\partial }{\partial \theta}\frac{\partial l}{\partial \theta } f(x|\theta) \ dx =0$
$\Rightarrow$$E[\frac{\partial^{2} l}{\partial \theta^{2} }] + \int \frac{\partial }{\partial \theta} \frac{\frac{\partial }{\partial \theta} f(x|\theta)}{f(x|\theta)} f(x|\theta) dx = 0$
$E[\frac{\partial^{2} l}{\partial \theta^{2} }] + \int \frac{\partial }{\partial \theta} \frac{\partial }{\partial \theta} f(x|\theta)\ dx = 0$
$E[\frac{\partial^{2} l}{\partial \theta^{2} }] + E[(\frac{\partial l}{\partial \theta })^{2}] = 0$
$\Rightarrow E[\frac{\partial^{2} l}{\partial \theta^{2} }] = – E[(\frac{\partial l}{\partial \theta })^{2}]$