Exponential family of distributions

• The density function for an exponential family distribution $$f\left(y\right)=\exp \{\frac{y\theta -b\left(\theta \right)}{a\left(\varphi \right)}+c(y,\varphi \left)\right\}$$
  • $a,b,c$: arbitrary functions
  • $\varphi$: an arbitrary scale parameter
  • $\theta$: the canonical parameter; completely depend on the model parameter $\beta$

• Properties about exponential family mean and variance $$\begin{array}{cc}E\left(Y\right)& =b^{\prime }\left(\theta \right)\\ var\left(Y\right)& =b^{\prime \prime }\left(\theta \right)a\left(\varphi \right)\end{array}$$

  • In most practical cases, $a\left(\varphi \right)=\varphi /\omega$ where $\omega$ is a known constant
  • We define a function $$V\left(\mu \right)=b^{\prime \prime }\left(\theta \right)/w,\text{so that}var\left(Y\right)=V\left(\mu \right)\varphi$$

Exponential family examples

Fitting GLMs

• For the GLM model and , assuming $a_i\left(\varphi \right)=\varphi /{\omega }_i$, the log likelihood is $$l\left(\beta \right)=\sum _{i=1}^n{\omega }_i[y_i{\theta }_i-b_i({\theta }_i\left)\right]/\varphi +c_i(\varphi ,y_i)$$

• To optimize, we use the Newton’s method, which is an iterative optimization approach $${\theta }^{(t+1)}={\theta }^{\left(t\right)}-{\left({\nabla }^{2}l\right)}^{-1}\nabla l$$
  • Where both ${\nabla }^{2}l$ and $\nabla l$ are evaluated at the current iteration ${\theta }^{\left(t\right)}$
  • Alternatively, we can use the Fisher scoring variant of the Newton’s method, by replacing the Hessian matrix with its expectation

• Next, we will need to compute the gradient vector and expected Hessian matrix of $l$

Compute the gradient vector and expected Hessian of $l$

• By the chain rule, $$\begin{array}{cc}\frac{\partial {\theta }_i}{\partial {\beta }_j}& =\frac{\partial {\theta }_i}{\partial {\mu }_i}\cdot \frac{\partial {\mu }_i}{\partial {\eta }_i}\cdot \frac{\partial {\eta }_i}{\partial {\beta }_j}\\ & =\frac{1}{b^{\prime \prime }\left({\theta }_i\right)}\cdot \frac{1}{g^{\prime }\left({\mu }_i\right)}\cdot X_{ij}\end{array}$$

• Therefore, the gradient vector of $l$ is $$\begin{array}{c}\frac{\partial l}{\partial {\beta }_j}=\frac{1}{\varphi }\sum _{i=1}^n{\omega }_i[y_i-b_i^{\prime }({\theta }_i\left)\right]\frac{\partial {\theta }_i}{\partial {\beta }_j}=\frac{1}{\varphi }\sum _{i=1}^n\frac{y_i-{\mu }_i}{g^{\prime }\left({\mu }_i\right)V\left({\mu }_i\right)}{X}_{ij}\end{array}$$

• The expected Hessian (expectation taken wrt $Y$) is $$E\left(\frac{{\partial }^{2}l}{\partial {\beta }_j\partial {\beta }_k}\right)=-\frac{1}{\varphi }\sum _{i=1}^n\frac{X_{ij}{X}_{ik}}{g^{\prime }\left({\mu }_i{)}^{2}V\right({\mu }_i)}$$

The Fisher scoring update

• Define the matrices $$\begin{array}{cc}W& =\text{diag}\left\{w_i\right\},w_i=\frac{1}{g^{\prime }\left({\mu }_i{)}^{2}V\right({\mu }_i)}\\ G& =\text{diag}\left\{g^{\prime }\right({\mu }_i\left)\right\}\end{array}$$

• The Fisher scoring update for $\beta$ is $$\begin{array}{cc}{\beta }^{(t+1)}& ={\beta }^{\left(t\right)}+{\left(X^{T}WX\right)}^{-1}{X}^{T}WG(y-\mu )\\ & ={\left(X^{T}WX\right)}^{-1}{X}^{T}W\underset{z}{\underset{}{\left[G\right(y-\mu )+X{\beta }^{\left(t\right)}]}}\end{array}$$

Iteratively re-weighted least square (IRLS) algorithm

1. Initialization: $${\hat{\mu }}_i=y_i+{\delta }_i,{\hat{\eta }}_i=g\left({\hat{\mu }}_i\right)$$

   • ${\delta }_i$ is usually zero, but may be a small constant ensuring ${\hat{\eta }}_i$ is finite

2. Compute pseudo data $z_i$ and iterative weights $w_i$: $$\begin{array}{cc}{z}_i& =g^{\prime }\left({\hat{\mu }}_i\right)(y_i-{\hat{\mu }}_i)+{\hat{\eta }}_i\\ w_i& =\frac{1}{g^{\prime }{\left({\hat{\mu }}_i\right)}^{2}V\left({\hat{\mu }}_i\right)}\end{array}$$

3. Find $\hat{\beta }$ by minimizing the weighted least squares objective $$\sum _{i=1}^n{w}_i{(z_i-X_i\beta )}^2$$ then update $$\hat{\eta }=X\hat{\beta },{\hat{\mu }}_i=g^{-1}\left({\hat{\eta }}_i\right)$$

• Repeat Step 2-3 until the change in deviance is near zero

IRLS example 1: logistic regression

• For logistic regression, $$\begin{array}{cc}g\left(\mu \right)& =\log \left(\frac{\mu }{1-\mu }\right),g^{\prime }\left(\mu \right)=\frac{1}{\mu (1-\mu )}\\ V\left(\mu \right)& =\mu (1-\mu ),\varphi =1\end{array}$$

• Therefore, in Step 2 of IRLS, $$\begin{array}{cc}{z}_i& =\frac{y_i-{\hat{\mu }}_i}{{\hat{\mu }}_i(1-{\hat{\mu }}_i)}+{\hat{\eta }}_i\\ w_i& ={\hat{\mu }}_i(1-{\hat{\mu }}_i)\end{array}$$

IRLS example 2: GLM with independent normal priors

• Assume that the vector $\beta$ has independent normal priors $$\beta \sim \text{N}(0,\frac{\varphi }{\lambda }{I}_p)$$

• Log posterior density (we still call it $l$, with some abuse of notation) $$l\left(\beta \right)=\frac{1}{\varphi }\sum _{i=1}^n{\omega }_i[y_i{\theta }_i-b_i({\theta }_i\left)\right]-\frac{\lambda }{2\varphi }{\beta }^T\beta +\text{const}$$

• Gradient vector and expected Hessian matrix (wrt $\beta$) $$\begin{array}{cc}\nabla l& =\frac{1}{\varphi }\left[X^{T}WG\right(y-\mu )-\lambda \beta ]\\ E\left({\nabla }^{2}l\right)& =-\frac{1}{\varphi }(X^{T}WX+\lambda I_p)\end{array}$$

  • Here, $W$ and $G$ are the same as in Equation and

• IRLS for GLM with independent normal priors $$\begin{array}{cc}{\beta }^{(t+1)}& ={\beta }^{\left(t\right)}+{(X^{T}WX+\lambda I_p)}^{-1}\left[X^{T}WG\right(y-\mu )-\lambda {\beta }^{\left(t\right)}]\\ & ={(X^{T}WX+\lambda I_p)}^{-1}{X}^{T}W\underset{z}{\underset{}{\left[G\right(y-\mu )+X{\beta }^{\left(t\right)}]}}\end{array}$$

Large sample distribution of $\hat{\beta }$

• Hessian of the negative log likelihood (also called observed information) $$\hat{I}=\frac{X^{T}WX}{\varphi }$$

• Fisher information, also called expected information $$I=E\left(\hat{I}\right)$$

• Asymptotic normality the MLE $\hat{\beta }$ $$\hat{\beta }\sim \text{N}(\beta ,I^{-1})\text{or}\hat{\beta }\sim \text{N}(\beta ,{\hat{I}}^{-1})$$

Deviance

• Deviance is the GLM counterpart of the residual sum of squares in normal linear regression $$\begin{array}{cc}D& =2\varphi [l\left({\hat{\beta }}_{max}\right)-l\left(\hat{\beta }\right)]\\ & =\sum _{i=1}^{n}2{\omega }_i[y_i({\tilde{\theta }}_i-{\hat{\theta }}_i)-b\left({\tilde{\theta }}_i\right)+b\left({\hat{\theta }}_i\right)]\end{array}$$

  • Here, $l\left({\hat{\beta }}_{max}\right)$ is the maximized likelihood of the saturated model: the model with one parameter per data point. For exponential family distribution, it is computed by simply setting $\hat{\mu }=y$.
  • $\tilde{\theta }$ and $\hat{\theta }$ are the maximum likelihood estimates of the canonical parameters for the saturated model and the model of interest, respectively

• From the second equality, we can see that deviance is independent of $\varphi$

• For normal linear regression, deviance equals the residual sum of squares

Scaled deviance

• Scaled deviance does depend on $\varphi$ $$D^{\ast }=\frac{D}{\varphi }$$

• If the model is specified correctly, then approximately $$D^{\ast }\sim {\chi }_{n-p}^2$$

• To compare two nested models,
  • If $\varphi$ is known, then under $H_0$, we can use $$D_0^{\ast }-D_1^{\ast }\sim {\chi }_{p_1-p_0}^2$$
  • If $\varphi$ is unknown, then under $H_0$, we can use $$F=\frac{(D_0-D_1)/(p_1-p_0)}{D_1/(n-p_1)}\sim F_{p_1-p_0,n-p_1}$$

Canonical link functions

• The canonical link $g_c$ is the link function such that $$g_c\left({\mu }_i\right)={\theta }_i={\eta }_i$$ where ${\theta }_i$ is the canonical parameter of the distribution

• Under canonical links, the observed information $\hat{I}$ and the expected information $I$ matrices are the same

• Under canonical links, since $\frac{\partial {\theta }_i}{\partial {\beta }_j}=X_{ij}$, the system of equations that the MLE satisfies becomes $$\frac{\partial l}{\partial {\beta }_j}=\sum _{i=1}^n{\omega }_i(y_i-{\mu }_i)X_{ij}=0$$ Thus, if ${\omega }_i=1$, we have $$X^{T}y=X^T\hat{\mu }$$

  • For any GLM with an intercept term and canonical link, the residuals sum to zero, i.e., ${\sum }_i{y}_i={\sum }_i{\hat{\mu }}_i$

GLM residuals

• Model checking is perhaps the most important part of applied statistical modeling

• It is usual to standardize GLM residuals so that if the model assumptions are correct,

  • the standardized residuals should have approximately equal variance, and
  • behave like residuals from an ordinary linear model

• Pearson residuals $${\widehat{ϵ}}_i^p=\frac{y_i-{\hat{\mu }}_i}{\sqrt{V\left({\mu }_i\right)}}$$

  • In practice, the distribution of the Pearson residuals can be quite asymmetric around zero. So the deviance residuals (introduced next) are often preferred.

Deviance residuals

• Denote $d_i$ as the $i$th component in the deviance definition , so that the deviance is $D={\sum }_{i=1}^n{d}_i$

• By analogy with the ordinary linear model,we define the deviance residual $${\widehat{ϵ}}_i^d=\text{sign}(y_i-\widehat{{\mu }_i})\sqrt{d_i}$$

  • The sum of squares of the deviance residuals gives the deviance itself

Quasi-likelihood

• Consider an observation $y_i$, of a random variable with mean ${\mu }_i$ and known variance function $V\left({\mu }_i\right)$

  • Getting the distribution of $Y_i$ exactly right is rather unimportant, as long as the mean-variance relationship $V(\cdot )$ is correct

• Then the log quasi-likelihood for ${\mu }_i$, given $y_i$, is $$q_i\left({\mu }_i\right)={\int }_{y_i}^{{\mu }_i}\frac{y_i-z}{\varphi V\left(z\right)}dz$$

  • The log quasi-likelihood for the mean vector $\mu$ of all the response data is $q\left(\mu \right)={\sum }_{i=1}^n{q}_i\left({\mu }_i\right)$

• To obtain the maximum quasi-likelihood estimation of $\beta$, we can differentiate $q$ wrt ${\beta }_j$, for $\forall j$ $$0=\frac{\partial q}{\partial {\beta }_j}=\sum _{i=1}^n\frac{y_i-{\mu }_i}{\varphi V\left({\mu }_i\right)}\frac{\partial {\mu }_i}{\partial {\beta }_j}⟹\sum _{i=1}^n\frac{y_i-{\mu }_i}{V\left({\mu }_i\right)g^{\prime }\left({\mu }_i\right)}{X}_{ij}=0$$ this is exactly the GLM maximum likelihood solution, which can be obtained through IRLS