Joint Scores For Regression And Dispersion

Why Treat Regression and Dispersion Together?

In a generalized linear model, the regression coefficients determine the mean structure, while the dispersion parameter controls the scale of variation around that structure. When dispersion is unknown, both sets of parameters belong to the likelihood problem. Their scores are not unrelated: the regression score is scaled by dispersion, and differentiating it with respect to dispersion gives a mixed derivative. Keeping the parameters together makes the joint likelihood equations and curvature explicit.

We work with independent observations \(Y_1,\ldots,Y_n\) from an exponential-dispersion model. The aim is to derive the score vector and the blocks of its derivative, not to solve a particular dispersion equation or to develop an iteration for estimating the regression coefficients.

Model and Joint Log-Likelihood

For observation \(i\), suppose the density or mass function has the form

f(y_i;\theta_i,\phi) = \exp\left\{ \frac{y_i\theta_i-b(\theta_i)}{\phi} +c(y_i,\phi) \right\},

where \(\phi>0\) is a common, unknown dispersion parameter. The support is assumed not to depend on \(\beta\) or \(\phi\), and the displayed expression is assumed to define a normalized distribution on that support. Let \(\beta\in\mathbb{R}^p\), let \(x_i\in\mathbb{R}^p\) be a fixed design vector, and set \(\eta_i=x_i^{\mathsf T}\beta\). A fixed offset could be added to \(\eta_i\) without changing the derivatives below. The link determines \(\theta_i\) through the mean \(\mu_i=b'(\theta_i)\) and the relation between \(\mu_i\) and \(\eta_i\).

Definition: Write \(r_i=y_i-\mu_i\), \(h_i=d\theta_i/d\eta_i\), and \(h_i'=d^2\theta_i/d\eta_i^2\). The contribution \(A_i=y_i\theta_i-b(\theta_i)\) is the part of observation \(i\)'s log-likelihood numerator that depends on its natural parameter. The joint score is the vector of first derivatives of the log-likelihood with respect to \((\beta^{\mathsf T},\phi)^{\mathsf T}\).

Assume \(b\) is twice differentiable over the relevant natural-parameter values, the link is differentiable enough for \(h_i\) and \(h_i'\) to exist, and \(c(y_i,\phi)\) is twice differentiable in \(\phi\). Since the sample is finite, the log-likelihood is a finite sum; differentiating that sum requires no interchange of a limit and an integral.

\ell(\beta,\phi) = \sum_{i=1}^n \left\{ \frac{y_i\theta_i-b(\theta_i)}{\phi} +c(y_i,\phi) \right\} = \sum_{i=1}^n\left\{\frac{A_i}{\phi}+c(y_i,\phi)\right\}.

Deriving the Two Score Components

First differentiate with respect to \(\beta\). By the chain rule,

\frac{\partial A_i}{\partial \beta} = \left\{y_i-b'(\theta_i)\right\} \frac{\partial \theta_i}{\partial\beta} = r_i h_i x_i.

Here \(b'(\theta_i)=\mu_i\), and \(\partial\theta_i/\partial\beta=h_i x_i\) because \(\eta_i=x_i^{\mathsf T}\beta\). Thus the regression score is the sum of the residual contributions, scaled by the common dispersion.

U_\beta(\beta,\phi) = \frac{\partial\ell}{\partial\beta} = \frac{1}{\phi}\sum_{i=1}^n x_i h_i r_i.

For the dispersion score, hold \(\beta\) fixed, so \(A_i\) does not change with \(\phi\). Differentiating both terms in each likelihood contribution gives

U_\phi(\beta,\phi) = \frac{\partial\ell}{\partial\phi} = -\frac{1}{\phi^2}\sum_{i=1}^n A_i +\sum_{i=1}^n c_\phi(y_i,\phi),

where \(c_\phi\) denotes the partial derivative with respect to the second argument. In general, the derivative of \(c\) is essential: it contains the part of the density normalization that depends on dispersion. Dropping it does not give the likelihood score.

Worked Example: Assemble the Score in a Canonical-Link Model

Suppose the link is canonical, so \(\theta_i=\eta_i=x_i^{\mathsf T}\beta\). Then \(h_i=1\). For two observations with \(x_1=(1,0)^{\mathsf T}\), \(x_2=(1,1)^{\mathsf T}\), residuals \(r_1=2\), \(r_2=-1\), and \(\phi=2\), calculate the regression score.

The weighted residual contributions are \(x_1r_1=(2,0)^{\mathsf T}\) and \(x_2r_2=(-1,-1)^{\mathsf T}\). Their sum is \((1,-1)^{\mathsf T}\). Therefore

U_\beta = \frac{1}{2}\{x_1r_1+x_2r_2\} = \frac{1}{2} \begin{pmatrix}1\\-1\end{pmatrix} = \begin{pmatrix}0.5\\-0.5\end{pmatrix}.

For the dispersion component, the data supplied in this calculation do not specify \(A_i\) or \(c_\phi(y_i,\phi)\), so its value cannot be determined from residuals alone. Its correct expression remains \(U_\phi=-\sum_i A_i/\phi^2+\sum_i c_\phi(y_i,\phi)\). This distinction matters: the regression score can be written in terms of residuals, but the dispersion score also depends on the full likelihood normalization.

The Joint Derivative Matrix

The observed Hessian is the matrix of second derivatives of \(\ell\), partitioned according to \(\beta\) and \(\phi\). The regression block follows by differentiating \(r_i h_i\) with respect to \(\eta_i\). Since \(d\mu_i/d\eta_i=b''(\theta_i)h_i\),

\frac{d(r_i h_i)}{d\eta_i} = -\frac{d\mu_i}{d\eta_i}h_i+r_i h_i' = -b''(\theta_i)h_i^2+r_i h_i'.

Consequently, the regression block is

H_{\beta\beta} = \frac{\partial^2\ell}{\partial\beta\,\partial\beta^{\mathsf T}} = \frac{1}{\phi}\sum_{i=1}^n x_i x_i^{\mathsf T} \left\{-b''(\theta_i)h_i^2+r_i h_i'\right\}.

For the mixed block, differentiate the regression score with respect to \(\phi\). At fixed \(\beta\), its residuals and \(h_i\)'s are fixed, and only the factor \(1/\phi\) changes. Differentiating the dispersion score with respect to \(\beta\) provides an independent check: \(\partial A_i/\partial\beta=r_i h_i x_i\), while \(c_\phi(y_i,\phi)\) has no \(\beta\)-dependence.

Result: Under the stated differentiability conditions, the two mixed blocks are transposes of one another:

H_{\beta\phi} = \frac{\partial^2\ell}{\partial\beta\,\partial\phi} = -\frac{1}{\phi^2}\sum_{i=1}^n x_i h_i r_i, \qquad H_{\phi\beta} = \frac{\partial^2\ell}{\partial\phi\,\partial\beta^{\mathsf T}} = -\frac{1}{\phi^2}\sum_{i=1}^n r_i h_i x_i^{\mathsf T}.

Indeed, \(H_{\phi\beta}=H_{\beta\phi}^{\mathsf T}\). Equivalently, the mixed block can be read directly from the regression score:

H_{\beta\phi}=-\frac{1}{\phi}U_\beta.

Finally, differentiating the dispersion score once more gives

H_{\phi\phi} = \frac{\partial^2\ell}{\partial\phi^2} = \frac{2}{\phi^3}\sum_{i=1}^n A_i +\sum_{i=1}^n c_{\phi\phi}(y_i,\phi).

These expressions form the full observed Hessian. The \(h_i'\) term is retained because for a noncanonical link the derivative of \(\theta_i\) with respect to \(\eta_i\) can vary with \(\eta_i\). For a canonical link \(h_i=1\) and \(h_i'=0\), but the mixed regression-dispersion block need not vanish.

Worked Example: Check the Mixed Derivative in a Log-Mean Gamma Model

Consider the gamma exponential-dispersion model with mean \(\mu_i\), variance \(\phi\mu_i^2\), and log link \(\eta_i=\log\mu_i\). Its natural parameter is \(\theta_i=-1/\mu_i\), so \(h_i=d\theta_i/d\eta_i=1/\mu_i\). Take an intercept-only model with two observations \(y_1=1\), \(y_2=2\), \(\mu_1=\mu_2=1\), and \(\phi=1\). Then \(x_1=x_2=1\), \(h_1=h_2=1\), and the residuals are \(r_1=0\), \(r_2=1\).

The regression score is

U_\beta = \frac{1}{1}\{(1)(1)(0)+(1)(1)(1)\}=1.

Differentiating that score with respect to dispersion gives \(H_{\beta\phi}=-U_\beta/\phi=-1\). Directly using the mixed-block formula gives

H_{\beta\phi} = -\frac{1}{1^2}\{(1)(1)(0)+(1)(1)(1)\} =-1.

The agreement checks both the sign and the dispersion power. In particular, the mixed block is not zero merely because the link is log: it is nonzero here because the regression score is nonzero.

Simultaneous Likelihood Equations

Suppose the parameter space for \((\beta,\phi)\) is open at a candidate maximum and the log-likelihood is differentiable there. A necessary condition for an interior local maximum is that every first partial derivative vanish. Thus the simultaneous likelihood equations are

\sum_{i=1}^n x_i h_i(y_i-\mu_i)=0, \qquad -\frac{1}{\phi^2}\sum_{i=1}^n \{y_i\theta_i-b(\theta_i)\} +\sum_{i=1}^n c_\phi(y_i,\phi)=0.

The factor \(1/\phi\) in the regression score has been removed from its equation because \(\phi>0\). Both equations must hold at the same parameter pair: the first uses means determined by \(\beta\), and the second also depends on those means through \(\theta_i\). A solution of the score equations is a stationary point, not automatically a maximum or even a global maximum. Boundary solutions require separate treatment because the interior derivative condition need not apply.

Worked Example: Write the Equations for a Two-Parameter Fit

Let \(\beta=(\beta_0,\beta_1)^{\mathsf T}\) and \(\eta_i=\beta_0+\beta_1 z_i\), with a differentiable link and common unknown \(\phi>0\). For each observation, define \(\mu_i=b'(\theta_i)\), where \(\theta_i\) is determined by \(\eta_i\), and let \(h_i=d\theta_i/d\eta_i\). The two components of the regression equation are obtained from the design vector \(x_i=(1,z_i)^{\mathsf T}\):

\sum_{i=1}^n h_i(y_i-\mu_i)=0, \qquad \sum_{i=1}^n z_i h_i(y_i-\mu_i)=0.

The dispersion equation is

-\frac{1}{\phi^2}\sum_{i=1}^n \{y_i\theta_i-b(\theta_i)\} +\sum_{i=1}^n c_\phi(y_i,\phi)=0.

Together these are three equations for \(\beta_0,\beta_1,\phi\). The first two determine the regression score conditions; the third cannot generally be replaced by a residual sum of squares or another generic residual formula. Its form depends on the density through both \(b\) and \(c\).

Common Pitfall and Practical Check

A frequent error is to differentiate the regression score with respect to \(\phi\) and conclude that the cross-derivative is zero because the fitted means do not change when \(\phi\) changes at fixed \(\beta\). The means indeed stay fixed in that partial derivative, but the score still contains \(1/\phi\). The correct calculation is

\frac{\partial}{\partial\phi} \left\{\frac{1}{\phi}\sum_i x_i h_i r_i\right\} = -\frac{1}{\phi^2}\sum_i x_i h_i r_i,

which is generally nonzero. It vanishes when the regression score itself vanishes, not simply because the mean model is held fixed. Another check is that differentiating the dispersion score with respect to \(\beta\) must yield the transpose of this block; if the signs or powers of \(\phi\) differ, revisit the chain rule.

Key takeaway: The joint score consists of a regression component scaled by \(1/\phi\) and a dispersion component containing both the natural-parameter likelihood contribution and the derivative of its normalizing term. The mixed block is \(-U_\beta/\phi\), and an interior likelihood solution must set both score components to zero simultaneously.

Check Your Understanding

Use the score and derivative formulas to answer the following questions.

Why does the dispersion score contain a derivative of \(c(y_i,\phi)\), even though \(c\) has no dependence on \(\beta\)?
For a canonical link, what are \(h_i\) and \(h_i'\), and which term disappears from \(H_{\beta\beta}\)?
Show from the regression score that \(H_{\beta\phi}=-U_\beta/\phi\).
What assumption about a candidate parameter value is needed to use zero scores as a necessary condition for a local maximum?
Why is a zero mixed block not guaranteed at an arbitrary parameter pair, even when the link is canonical?