5.2 Algorithm

The parameter vector $\beta$ of the GLM is estimated using a solution to th log-likelihood function as follows. In the canonical form, for independent observations $y_1,\mathrm{\dots },y_n$, the likelihood is given by:

$$\mathrm{\ell }({\beta }_1,\mathrm{\dots },{\beta }_p)=\sum _{i=1}^n(y_i{\theta }_i-\kappa ({\theta }_i))+\sum _{i=1}^n\log q(y_i).$$

Next we take the first derivative of $\mathrm{\ell }$ with respect to ${\beta }_r$, $1\le r\le p$. Note here that ${\theta }_i$ is a one-to-one function of ${\mu }_i$ with $\frac{d{\mu }_i}{d{\theta }_i}={\kappa }^{\prime \prime }({\theta }_i)$ and ${\eta }_i$ is a one-to-one function of ${\mu }_i$ through the link function and finally ${\eta }_i={𝐱}_i^{\prime }𝜷$. Hence, by the chain rule:

$$\frac{d\mathrm{\ell }}{d{\beta }_r}=\sum _{i=1}^n\frac{d{\mathrm{\ell }}_i}{d{\theta }_i}\frac{d{\theta }_i}{d{\mu }_i}\frac{d{\mu }_i}{d{\eta }_i}\frac{d{\eta }_i}{d{\beta }_r}$$

where

${\displaystyle \frac{d{\mathrm{\ell }}_i}{d{\theta }_i}}=y_i-{\kappa }^{\prime }({\theta }_i)=y_i-{\mu }_i\mathit{\hspace{1em}\hspace{1em}}{\displaystyle \frac{d{\mu }_i}{d{\theta }_i}}={\kappa }^{\prime \prime }({\theta }_i)=\mathrm{var}(Y_i)\mathit{\hspace{1em}\hspace{1em}}{\displaystyle \frac{d{\eta }_i}{d{\beta }_r}}=x_{i,r}.$

This leads to the likelihood equations:

$$\frac{d\mathrm{\ell }}{d{\beta }_r}=\sum _{i=1}^n\frac{(y-{\mu }_i)x_{i,r}}{\mathrm{var}(Y_i)}\frac{d{\mu }_i}{d{\eta }_i},\mathrm{for}1\le r\le p.$$

We denote the above likelihood vector form by:

$$𝐮=\sum _{i=1}^n(y_i-{\mu }_i){𝐱}_i\frac{1}{\mathrm{var}(Y_i)}{\left(\frac{d{\mu }_i}{d{\eta }_i}\right)}^2\frac{d{\eta }_i}{d{\mu }_i}=\sum _{i=1}^n{W}_i(y_i-{\mu }_i)\frac{d{\eta }_i}{d{\mu }_i}{𝐱}_i,$$

where

$W_i={\displaystyle \frac{1}{\mathrm{var}(Y_i)}}{\left({\displaystyle \frac{d{\mu }_i}{d{\eta }_i}}\right)}^2={\displaystyle \frac{1}{\mathrm{var}(Y_i){\{g^{\prime }({\mu }_i)\}}^2}}.$

(5.2)

In the sequel, all sums are over $i$ from 1 to n, unless otherwise and the subscript $i$ is omitted from the summands. The $r$-the component, $1\le r\le p$ of $𝐮$ is:

$u_r={\displaystyle \sum W(y-\mu )\frac{d\eta }{d\mu }{x}_r}$

(5.3)

and let the expectation of the negative Hessian matrix be:

$$𝐈=𝔼\left[-\frac{d{}^2\mathrm{\ell }}{d{\beta }_{r}d{\beta }_s}\right],$$

where both $𝐮$ and $𝐈$ are evaluated at the current estimate $𝐛$ of $𝜷$. Then from (5.3), for $1\le r,s\le p$,

	$-{\displaystyle \frac{d{u}_r}{d{\beta }_s}}$	$=-{\displaystyle \sum \left[(y-\mu )\frac{d}{d{\beta }_s}\left\{W\frac{d\eta }{d\mu }\right\}{x}_r+W\frac{d\eta }{d\mu }{x}_r\frac{d}{d{\beta }_s}(y-\mu )\right]}$
		$=-{\displaystyle \sum \left[(y-\mu )\frac{d}{d{\beta }_s}\left\{W\frac{d\eta }{d\mu }\right\}{x}_r-W{x}_r\frac{d\eta }{d\mu }\frac{d\mu }{d{\beta }_s}\right]}$
		$=-{\displaystyle \sum \left[(y-\mu )\frac{d}{d{\beta }_s}\left\{W\frac{d\eta }{d\mu }\right\}{x}_r\right]}+{\displaystyle \sum W{x}_r\frac{d\eta }{d{\beta }_s}}$
		$=-{\displaystyle \sum \left[(y-\mu )\frac{d}{d{\beta }_s}\left\{W\frac{d\eta }{d\mu }\right\}{x}_r\right]}+{\displaystyle \sum W{x}_r{x}_s}.$

Hence:

$I_{r,s}=-𝔼\left[{\displaystyle \frac{d{u}_r}{d{\beta }_s}}\right]=-{\displaystyle \sum \left[𝔼[Y-\mu ]\frac{d}{d{\beta }_s}\left\{W\frac{d\eta }{d\mu }\right\}{x}_r\right]}+{\displaystyle \sum W{x}_r{x}_s}={\displaystyle \sum W{x}_r{x}_s}.$

(5.4)

To apply Fisher’s scoring method, note that the $r$-th component of $\mathrm{𝐈𝐛}$ is:

${\displaystyle \sum _s}{I}_{r,s}{b}_s={\displaystyle \sum _s}{\displaystyle \sum _{i=1}^n}{W}_i{x}_{i,r}{x}_{i,s}{b}_s={\displaystyle \sum _{i=1}^n}{W}_i{x}_{i,r}{\displaystyle \sum _s}{x}_{i,s}{b}_s={\displaystyle \sum _{i=1}^n}{W}_i{x}_{i,r}{\eta }_i$

(5.5)

where ${\eta }_i={𝐱}_i^{\prime }𝐛$ is the $i$-th linear predictor evaluated estimate. Hence from (5.3) and (5.5):

$\mathrm{𝐈𝐛}+𝐮={\displaystyle \sum W\left\{\eta +(y-\mu )\frac{d\eta }{d\mu }\right\}𝐱}={\displaystyle \sum Wz𝐱},$

(5.6)

where

$z_i=z_i(𝐛)={\eta }_i+(y_i-{\mu }_i){\displaystyle \frac{d{\eta }_i}{d{\mu }_i}}={𝐱}_i^{\prime }𝐛+(y_i-{\mu }_i)g^{\prime }({\mu }_i)$

(5.7)

with all quantities (${\mu }_i$ and ${\eta }_i$) evaluated at the current estimate $𝐛$. Consequently from (5.1), (5.5) and (5.6):

$${𝐛}^{*}={𝐈}^{-1}(\mathrm{𝐈𝐛}+𝐮)={\left(\sum _{i=1}^n{W}_i{𝐱}_i^{\prime }{𝐱}_i\right)}^{-1}\left(\sum _{i=1}^n{W}_i{𝐱}_i^{\prime }{𝐳}_i\right)={({𝐗}^{\prime }\mathrm{𝐖𝐗})}^{-1}{𝐗}^{\prime }\mathrm{𝐖𝐳}$$

where $𝐖$ is a diagonal matrix with $i$-th diagonal entry $W_i$ of (5.2).

Remark 1 on implementation: For implementing IRWLS, start with initial $𝐛$ and first compute the linear predictor ${\eta }_i={𝐱}_i^{\prime }𝐛$. Then calculate ${\mu }_i=g^{-1}({\eta }_i)$; however, often the initial ${\mu }_i$’s are taken as $Y_i$’s and one evaluates ${\eta }_i=g({\mu }_i)$ (There are obvious problems with such choice, for example, when some $y_i$’s are zero and one has to take the logarithm as in the Poisson case). Finally, $z_i(b)$ of (5.7) is evaluated and the iteration continues.

Remark 2: Consider a multiple linear regression model with observed $z_i$’s of (5.7) defined as:

	$z_i={𝐱}_i^{\prime }𝜷+(y_i-{\mu }_i){\displaystyle \frac{d{\eta }_i}{d{\mu }_i}}$	$={𝐱}_i^{\prime }𝜷+{\displaystyle \frac{(y_i-{\mu }_i)}{{\{\mathrm{var}(Y_i)\}}^{1/2}}}\left[{\displaystyle \frac{{\{\mathrm{var}(Y_i)\}}^{1/2}}{\frac{d{\mu }_i}{d{\eta }_i}}}\right]$
		$={𝐱}_i^{\prime }𝜷+{\displaystyle \frac{(y_i-{\mu }_i)}{{\{\mathrm{var}(Y_i)\}}^{1/2}}}{W}_i^{-1/2}$		(5.8)

Since $𝔼[(Y_i-{\mu }_i)/\mathrm{var}{(Y_i)}^{1/2}]=0$, the updated estimate ${𝐛}^{*}$ is nothing but the weighted least squares estimate of $𝜷$ with weights given by $W_i$’s, where these weights and $z_i$’s are calculated using the current value of $𝐛$ since that is the best approximation at the current stage. The hypothetical model (5.8) can be motivated from a one-step Taylor approximation:

$$g(y_i)\approx g({\mu }_i)+(y_i-{\mu }_i)g^{\prime }({\mu }_i)={𝐱}_i^{\prime }𝜷+(y_i-\mu )\frac{d{\eta }_i}{d{\mu }_i}$$

Remark 3: From (5.4), $-𝔼[\frac{d{u}_r}{d{\beta }_s}]=-\frac{d{u}_r}{d{\beta }_s}$ if $W_i\frac{d{\eta }_i}{d{\mu }_i}$ is a constant function of $𝜷$. This happens under the canonical link function and consequently Fisher’s scoring method and Newton-Raphson method for finding $\widehat{𝜷}$ coincide resulting in fast convergence. This is because:

$$W_i\frac{d{\eta }_i}{d{\mu }_i}=\frac{1}{\mathrm{var}(Y_i)g^{\prime }({\mu }_i)}=\frac{1}{v({\mu }_i)g^{\prime }({\mu }_i)}$$

and for this to be free from $𝜷$, $g^{\prime }(\mu )v(\mu )$ is a constant or:

$$g(\mu )=c\int \frac{1}{v(\mu )}𝑑\mu .$$

In particular:

• •

  Simple linear regression – $v(\mu )=1$, $g(\mu )=\mu$

• •

  Poisson regression – $v(\mu )=\mu$, $g(\mu )=\log (\mu )$

• •

  Logistic regression – $v(\mu )=\mu (1-\mu )$, $g(\mu )=\log \{\mu /(1-\mu )\}$