Gaussian Identity-Link Score Equations

Gaussian Regression With an Identity Link

For a Gaussian generalized linear model with an identity link, the conditional mean is linear in the regression coefficients. This makes the coefficient likelihood equations especially direct: they are the same equations as ordinary least squares. The variance parameter still matters for the likelihood and its curvature, but, when a finite joint maximum exists, it does not change the coefficient fit.

Assume \(Y_1,\ldots,Y_n\) are independent, with fixed design vectors \(\mathbf{x}_i\in\mathbb{R}^p\), and \(Y_i\sim N(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta},\sigma^2)\), where \(\boldsymbol{\beta}\in\mathbb{R}^p\) and \(\sigma^2>0\). Write \(\mathbf{X}\) for the \(n\times p\) design matrix with rows \(\mathbf{x}_i^\mathsf{T}\), and \(\mathbf{y}\) for the vector of responses. The identity link means that the linear predictor equals the mean: \(\mu_i=\eta_i=\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\). The design is treated as fixed throughout.

Coefficient Score and Hessian

For a fixed \(\sigma^2>0\), independence and the Gaussian density give the log-likelihood, up to no omitted terms,

To differentiate with respect to \(\boldsymbol{\beta}\), expand the quadratic term as \(\mathbf{y}^\mathsf{T}\mathbf{y}-2\boldsymbol{\beta}^\mathsf{T}\mathbf{X}^\mathsf{T}\mathbf{y} +\boldsymbol{\beta}^\mathsf{T}\mathbf{X}^\mathsf{T}\mathbf{X}\boldsymbol{\beta}\). Since \(\mathbf{X}^\mathsf{T}\mathbf{X}\) is symmetric, its derivative yields the following score and observed Hessian.

Equivalently, the \(j\)th component of the score is \(\sigma^{-2}\sum_{i=1}^n x_{ij}(y_i-\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})\). Thus each score component is a design-weighted sum of residuals. A coefficient score equation sets each such sum to zero. The Hessian is constant in \(\boldsymbol{\beta}\): it depends on the design and the variance, not on the observed response or the current coefficient values.

Why the Normal Equations Give the Fit

The normal equations do more than identify a stationary point. For any fixed \(\sigma^2>0\), maximizing the coefficient log-likelihood is equivalent to minimizing \(\operatorname{RSS}(\boldsymbol{\beta})\). The next result gives the conditions for the usual unique coefficient solution.

Proof. Full column rank implies that \(\mathbf{X}^\mathsf{T}\mathbf{X}\) is positive definite: for every nonzero \(\mathbf{a}\in\mathbb{R}^p\), \(\mathbf{a}^\mathsf{T}\mathbf{X}^\mathsf{T}\mathbf{X}\mathbf{a}=\|\mathbf{X}\mathbf{a}\|^2>0\). It is therefore invertible, and the normal equations have the stated unique solution. For any \(\boldsymbol{\beta}\), set \(\boldsymbol{\delta}=\boldsymbol{\beta}-\hat{\boldsymbol{\beta}}\). The fitted residual vector is orthogonal to every column of \(\mathbf{X}\), because \(\mathbf{X}^\mathsf{T}(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}})=\mathbf{0}\). Expanding the squared norm gives

The middle term vanishes by orthogonality. The final term is nonnegative and is zero only if \(\mathbf{X}\boldsymbol{\delta}=\mathbf{0}\), which, by full column rank, requires \(\boldsymbol{\delta}=\mathbf{0}\). Thus the solution uniquely minimizes RSS. The log-likelihood differs from \(-\operatorname{RSS}(\boldsymbol{\beta})/(2\sigma^2)\) only by a term independent of \(\boldsymbol{\beta}\), so it has the same unique maximizer. \(\square\)

If \(\mathbf{X}\) is rank deficient, the inverse formula is unavailable and the coefficient vector need not be unique. The normal equations still characterize least-squares solutions, but distinct solutions can give the same fitted values and likelihood. A uniqueness claim therefore requires full column rank.

Likelihood Variance and Unbiased Residual Variance

The coefficient equations and the variance estimation problem are distinct. With full column rank, let \(p\) be the number of columns of \(\mathbf{X}\), including an intercept if present, and assume \(n>p\). The residual-variance estimate \(\operatorname{RSS}/(n-p)\) is unbiased under the Gaussian model. By contrast, when the joint likelihood has a maximum over \(\sigma^2>0\), its variance estimate is \(\operatorname{RSS}/n\).

For clarity, the likelihood distinction can be checked by holding the least-squares coefficients fixed. The part of the log-likelihood that depends on \(\sigma^2\) is \(-\frac{n}{2}\log(\sigma^2)-\frac{\operatorname{RSS}}{2\sigma^2}\). If \(\operatorname{RSS}>0\), its derivative is \(-n/(2\sigma^2)+\operatorname{RSS}/(2\sigma^4)\), which is zero at \(\hat{\sigma}^2_{\mathrm{ML}}=\operatorname{RSS}/n\). The second derivative at that point is negative, and the log-likelihood tends downwards at both ends of the positive variance range, so this is the maximum. If instead \(\operatorname{RSS}=0\), the log-likelihood tends to \(+\infty\) as \(\sigma^2\) decreases to zero. Since zero is not in the parameter space, no MLE for \(\sigma^2>0\) is attained.

Derivation of unbiasedness. Let \(\mathbf{H}=\mathbf{X}(\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}\) and \(\mathbf{M}=\mathbf{I}-\mathbf{H}\). The fitted residual vector is \(\mathbf{M}\mathbf{y}\). The matrix \(\mathbf{M}\) is symmetric and idempotent, and has rank \(n-p\), because \(\mathbf{H}\) is the projection onto the \(p\)-dimensional column space of \(\mathbf{X}\). Under the model, \(\mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}\), where \(E[\boldsymbol{\varepsilon}]=\mathbf{0}\) and \(\operatorname{Var}(\boldsymbol{\varepsilon})=\sigma^2\mathbf{I}\). Since \(\mathbf{M}\mathbf{X}=\mathbf{0}\), \(\operatorname{RSS}=\boldsymbol{\varepsilon}^\mathsf{T}\mathbf{M}\boldsymbol{\varepsilon}\). Taking expectation gives \(E[\operatorname{RSS}]=\operatorname{tr}(\mathbf{M}\operatorname{Var}(\boldsymbol{\varepsilon}))=\sigma^2\operatorname{tr}(\mathbf{M})=\sigma^2(n-p)\). Dividing by \(n-p\) proves the claim. \(\square\)

In the two-coefficient example above, \(n=3\), \(p=2\), and \(\operatorname{RSS}=1/6>0\). Thus the variance MLE is \((1/6)/3=1/18\), while the unbiased residual-variance estimate is \((1/6)/(3-2)=1/6\). The two denominators answer different statistical questions; the unbiased estimate is not obtained by solving the Gaussian likelihood equation.

Common Mistakes and Practical Checks

The identity-link Gaussian case is simple enough that errors can be easy to overlook. Keep the following distinctions in view when deriving or checking a fit.

• Dropping the variance factor in the score. The coefficient score includes \(1/\sigma^2\). This factor does not change the zero of the score when \(\sigma^2>0\), but it belongs in the score and Hessian.
• Calling any solution unique. The inverse formula and uniqueness require full column rank. With redundant columns, coefficient vectors can differ while producing identical fitted means.
• Using the unbiased denominator for the MLE. \(\operatorname{RSS}/(n-p)\) is unbiased under the stated conditions; it is not the Gaussian variance MLE, which is \(\operatorname{RSS}/n\) when \(\operatorname{RSS}>0\).
• Assuming \(n>p\) alone guarantees a variance MLE. It does not exclude an exact fit. If \(\operatorname{RSS}=0\), the likelihood has no maximum over strictly positive variances.
• Checking only the intercept residual sum. In a model with several design columns, the score requires the residual vector to be orthogonal to every column of \(\mathbf{X}\), not merely to sum to zero.

For a numerical fit, a useful check is to calculate \(\mathbf{X}^\mathsf{T}(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}})\). It should be the zero vector up to rounding error. Then report the variance estimate with its purpose and conditions stated: likelihood-based when \(\operatorname{RSS}>0\), or unbiased when full column rank and \(n>p\).