WWS509 - Lecture Notes

B.2 Maximum Likelihood Estimation

An important practical feature of generalized linear models is that they can all be fit to data using the same algorithm, a form of iteratively re-weighted least squares. In this section we describe the algorithm.

Given a trial estimate of the parameters [^(b)], we calculate the estimated linear predictor [^(h_i)] = x_i[^(b)] and use that to obtain the fitted values [^(m_i)] = g^-1([^(h_i)]). Using these quantities, we calculate the working dependent variable

z_i =

^
h_i

+ (y_i -

^
m_i

)

d h_i

d m_i

(B.6)

where the rightmost term is the derivative of the link function evaluated at the trial estimate.

Next we calculate the iterative weights

w_i = p_i/[ b(q_i) (

dh_i

dm_i

)²],

(B.7)

where b(q_i) is the second derivative of b(q_i) evaluated at the trial estimate and we have assumed that a_i(f) has the usual form f/p_i. This weight is inversely proportional to the variance of the working dependent variable z_i given the current estimates of the parameters, with proportionality factor f.

Finally, we obtain an improved estimate of b regressing the working dependent variable z_i on the predictors x_i using the weights w_i, i.e. we calculate the weighted least-squares estimate

^
b

= (XWX)^-1 XWz,

(B.8)

where X is the model matrix, W is a diagonal matrix of weights with entries w_i given by (B.7) and z is a response vector with entries z_i given by (B.6).

The procedure is repeated until successive estimates change by less than a specified small amount. McCullagh and Nelder (1989) prove that this algorithm is equivalent to Fisher scoring and leads to maximum likelihood estimates. These authors consider the case of general a_i(f) and include f in their expression for the iterative weight. In other words, they use w^*_i = fw_i, where w_i is the weight used here. The proportionality factor f cancels out when you calculate the weighted least-squares estimates using (B.8), so the estimator is exactly the same. I prefer to show f explicitly rather than include it in W.

Example: For normal data with identity link h_i = m_i, so the derivative is dh_i/ dm_i = 1 and the working dependent variable is y_i itself. Since in addition b(q_i) = 1 and p_i = 1, the weights are constant and no iteration is required. ^[¯]

In Sections B.4 and B.5 we derive the working dependent variable and the iterative weights required for binomial data with link logit and for Poisson data with link log. In both cases iteration will usually be necessary.

Starting values may be obtained by applying the link to the data, i.e. we take [^(m_i)] = y_i and [^(h_i)] = g([^(m_i)]). Sometimes this requires a few adjustments, for example to avoid taking the log of zero, and we will discuss these at the appropriate time.