Maximum-Likelihood Estimation of Gaussian Dispersion

Why Estimate Dispersion After Fitting the Mean?

In a Gaussian regression model, the mean structure describes how the response changes with predictors, while the dispersion describes the remaining variation around that mean. Both are unknown in ordinary regression: the mean coefficients are estimated from the data, and dispersion must be estimated after accounting for the fitted mean. The likelihood provides a direct derivation of this estimate.

The key point is that fitting coefficients changes the residual sum of squares, but it does not change how many response densities contribute the Gaussian likelihood's dispersion term. Under the usual independent-observation model, that term involves the sample size \(n\). This is why the maximum-likelihood denominator is \(n\), not the residual degrees of freedom.

Model and Likelihood

Suppose \(Y_1,\ldots,Y_n\) are independent observations with \(Y_i\sim N(x_i^\mathsf{T}\beta,\sigma^2/w_i)\). Here \(x_i\) is a fixed \(p\)-dimensional covariate vector, \(\beta\) is a \(p\)-dimensional coefficient vector, and each known weight \(w_i\) is positive. Thus \(w_i\) is a relative precision: larger weights mean smaller conditional variances. The common scale \(\sigma^2>0\) is the dispersion to be estimated.

Let \(X\) be the design matrix whose \(i\)th row is \(x_i^\mathsf{T}\), and let \(W=\operatorname{diag}(w_1,\ldots,w_n)\). Assume \(X\) has full column rank \(p\), so the weighted least-squares coefficient estimate is uniquely defined. The density of observation \(i\) is

Independence makes the joint likelihood the product of these densities. Taking its logarithm gives

The weights contribute the term \(\frac12\sum_i\log w_i\), which is constant in both unknown parameters because the weights are known. The coefficient-dependent part of the likelihood is determined by \(Q(\beta)\).

Profiling Out the Regression Coefficients

For any fixed \(\sigma^2>0\), the terms involving \(\beta\) are exactly \(-Q(\beta)/(2\sigma^2)\). Consequently, maximizing the likelihood over \(\beta\) is equivalent to minimizing \(Q(\beta)\). The minimizer is \(\hat\beta=(X^\mathsf{T}WX)^{-1}X^\mathsf{T}WY\), since full column rank of \(X\) and positive weights make \(X^\mathsf{T}WX\) positive definite.

This is a profile likelihood: at each candidate dispersion, the coefficients have first been set to their maximizing values. Importantly, the profile likelihood retains the same \(n\) response-density contributions as the joint likelihood.

The Maximum-Likelihood Estimate

The estimate follows by differentiating the profile log-likelihood. The result below also states the condition needed for a finite, unique maximum.

Proof. Differentiating with respect to \(\sigma^2\), not \(\sigma\), gives

Since \(Q>0\), the denominator is positive throughout the parameter space. The derivative is positive when \(0<\sigma^2<Q/n\), is zero at \(Q/n\), and is negative when \(\sigma^2>Q/n\). Thus the profile log-likelihood increases up to \(Q/n\) and decreases afterward, proving that this point is its unique global maximum. Therefore \(\hat\sigma^2_{\mathrm{ML}}=Q/n\). \(\square\)

The second derivative confirms local curvature at the maximum: \[ \left.\frac{d^2\ell_p}{d(\sigma^2)^2}\right|_{\sigma^2=Q/n} =-\frac{n}{2(Q/n)^2}<0. \] The sign-change argument establishes the global conclusion; the second derivative is an additional check.

Why the Denominator Is \(n\)

The appearance of \(n\) can be read directly from the likelihood. Each of the \(n\) independent normal densities includes a factor proportional to \((\sigma^2)^{-1/2}\). Multiplying these factors gives \((\sigma^2)^{-n/2}\), and taking logarithms gives the term \(-\frac n2\log(\sigma^2)\). The other dispersion-dependent term is \(-Q/(2\sigma^2)\). Setting the derivative of these two terms to zero yields \(Q/n\).

Estimating \(p\) regression coefficients changes the residuals and therefore changes \(Q\). It does not change the fact that the likelihood being maximized is the joint density of \(n\) responses. Replacing \(n\) by \(n-p\) in the score equation would therefore no longer be maximizing this likelihood. A denominator adjusted for a different criterion, such as unbiasedness, is not the maximum-likelihood estimate derived here.

When the Estimate Is Not Finite

The condition \(Q>0\) matters. If the observed response vector lies exactly in the column space of \(X\), then the fitted residuals are all zero and \(Q=0\). In that case the profile log-likelihood is \(C-\frac n2\log(\sigma^2)\), which increases without bound as \(\sigma^2\) approaches zero. Since the model requires \(\sigma^2>0\), there is no finite maximum in the parameter space. Thus one should not mechanically report \(Q/n\) as an interior maximum when \(Q=0\).

For a full-rank design with \(n>p\), exact fitting is a special data configuration, but it remains a mathematical boundary case. If \(n=p\), a full-rank design fits every response vector exactly, so this likelihood cannot produce a finite interior dispersion estimate. A practical calculation should check the residual sum of squares rather than assume it is positive.

Common Mistake

Wrong reasoning: “The model estimates \(p\) coefficients, so \(p\) degrees of freedom have been used. Therefore the maximum-likelihood estimate of \(\sigma^2\) must be \(Q/(n-p)\).”

Correction: That argument substitutes a degrees-of-freedom adjustment into a likelihood score that must be derived from the stated model. The profile log-likelihood contains \(-\frac n2\log(\sigma^2)-Q/(2\sigma^2)\); differentiating it gives \(Q/n\). Estimating coefficients affects \(Q\), but does not change the coefficient \(n\) in the likelihood's logarithmic term. Keep distinct the likelihood maximizer and estimates defined by other criteria.