Profiling Dispersion in Gaussian and Gamma GLMs

Why Profile Dispersion?

A dispersion parameter controls the variability around a model’s mean structure. In a Gaussian GLM it is the error variance; in the usual Gamma GLM parameterization it multiplies the squared mean in the variance function. An interval for dispersion should account for uncertainty in the regression coefficients rather than treating their estimates as known.

Profiling does this by fixing a candidate dispersion value and maximizing the likelihood over the regression coefficients. The resulting profile likelihood is a one-dimensional summary: it records how well the best-fitting regression model can fit the data at each candidate dispersion. Its likelihood-ratio contour gives an interval that need not be symmetric around the maximum-likelihood estimate.

Definition: Write the parameter as $(\beta,\phi)$, where $\beta$ is a vector of regression coefficients and $\phi>0$ is dispersion. The profile log-likelihood for dispersion is $$ \ell_p(\phi)=\sup_{\beta\in\mathcal B}\ell(\beta,\phi), $$ where $\mathcal B$ is the allowed coefficient space. If the joint maximum is attained at $(\hat\beta,\hat\phi)$, the profile likelihood-ratio statistic for a candidate $\phi$ is $$ W(\phi)=2\{\ell_p(\hat\phi)-\ell_p(\phi)\}. $$ A usual approximate $100(1-\alpha)\%$ profile likelihood interval is the set of $\phi>0$ satisfying $W(\phi)\leq \chi^2_{1,1-\alpha}$.

The chi-square cutoff is the Wilks likelihood-ratio cutoff. Its usual justification requires an identifiable model, a true parameter in the interior of the parameter space, a sufficiently smooth likelihood, and nonsingular Fisher information, along with regularity conditions for the likelihood’s local quadratic approximation. The resulting interval is generally an asymptotic interval, not an exact finite-sample interval. If the set crosses the cutoff at two points, those crossings are its endpoints.

Gaussian Regression: The Profile Has a Closed Form

Consider independent observations $Y_i\sim N(x_i^T\beta,\phi)$, for $i=1,\ldots,n$. Suppose the design matrix $X$ has full column rank, and let $p$ be its number of columns. For fixed $\phi$, maximizing the likelihood over $\beta$ is equivalent to minimizing the residual sum of squares. The minimizing coefficient vector is the ordinary least-squares estimate $\hat\beta$, and write $$ \operatorname{RSS}=\sum_{i=1}^n(y_i-x_i^T\hat\beta)^2. $$ The important point for profiling is that this constrained maximizer is the same for every fixed $\phi>0$.

Result: For Gaussian regression with unknown common variance $\phi$, the profile log-likelihood, up to a constant not depending on $\phi$, is $$ \ell_p(\phi)=-\frac n2\log\phi-\frac{\operatorname{RSS}}{2\phi}. $$ If $\operatorname{RSS}>0$, its maximum occurs at $\hat\phi=\operatorname{RSS}/n$, and the profile likelihood-ratio statistic is $$ W(\phi)=n\left\{\log\left(\frac{\phi}{\hat\phi}\right)+\frac{\hat\phi}{\phi}-1\right\}. $$

Derivation. The Gaussian log-likelihood is $$ \ell(\beta,\phi)=-\frac n2\log(2\pi\phi)-\frac{1}{2\phi}\sum_{i=1}^n(y_i-x_i^T\beta)^2. $$ At fixed $\phi$, the first term is constant in $\beta$, and the second is largest when the sum of squared residuals is smallest. Substituting its minimum, $\operatorname{RSS}$, gives the profile formula. Differentiating with respect to $\phi$ gives $$ \ell_p'(\phi)=-\frac{n}{2\phi}+\frac{\operatorname{RSS}}{2\phi^2}. $$ This derivative is zero at $\hat\phi=\operatorname{RSS}/n$; it is positive below that value and negative above it. Thus this point is the maximum. Subtracting the profile log-likelihood at $\phi$ from that at $\hat\phi$, and multiplying by two, gives the displayed expression for $W(\phi)$. $\square$

Worked Example: A Gaussian Dispersion Interval

Suppose a full-rank Gaussian regression has $n=10$ observations and residual sum of squares $\operatorname{RSS}=18$. Then $$ \hat\phi=\frac{18}{10}=1.8. $$ For a 95% likelihood-ratio interval, use the approximate cutoff $\chi^2_{1,0.95}=3.841$. Setting $W(\phi)=3.841$, and writing $r=\phi/\hat\phi$, gives $$ 10\{\log r+1/r-1\}=3.841, \qquad \log r+1/r-1=0.3841. $$ Numerically, the two positive roots are approximately $r=0.46529$ and $r=2.7884$. Therefore the dispersion endpoints are $$ 1.8(0.46529)=0.8375,\qquad 1.8(2.7884)=5.0191. $$ The approximate profile likelihood interval is $(0.8375,5.0191)$. Its unequal distances from $1.8$ are not an error: the profile log-likelihood is not a symmetric quadratic function of $\phi$ over this range.

Gamma GLMs: Profiling the Coefficients

Use the Gamma parameterization $$ Y_i\sim\operatorname{Gamma}\left(\text{shape}=1/\phi,\ \text{scale}=\phi\mu_i\right), \qquad E[Y_i]=\mu_i,\quad \operatorname{Var}(Y_i)=\phi\mu_i^2. $$ Assume independent positive observations $y_i$, and a mean model $\mu_i=\exp(x_i^T\beta)$. The log link ensures $\mu_i>0$. At fixed $\phi$, the Gamma likelihood must still be maximized over $\beta$; the profile is not obtained by simply substituting a coefficient estimate from a different model.

Result: In this Gamma model, maximizing the likelihood over $\beta$ at any fixed $\phi>0$ is equivalent to minimizing $$ Q(\beta)=\sum_{i=1}^n\left\{\log\mu_i+\frac{y_i}{\mu_i}\right\}. $$ Consequently, if this minimization has a finite maximizer, its maximizing coefficient vector is the same for every $\phi>0$. The profile calculation still requires maximizing over $\beta$; the invariance follows from this particular likelihood and mean model.

Derivation. The Gamma density with shape $1/\phi$ and scale $\phi\mu_i$ yields $$ \ell(\beta,\phi)=-n\log\Gamma(1/\phi)-\frac n\phi\log\phi -\frac1\phi\sum_i\log\mu_i +\left(\frac1\phi-1\right)\sum_i\log y_i -\frac1\phi\sum_i\frac{y_i}{\mu_i}. $$ For fixed $\phi$, all terms except those involving $\mu_i$ are constant in $\beta$. The remaining terms are $-Q(\beta)/\phi$, so maximizing the likelihood is equivalent to minimizing $Q$; multiplication by the positive constant $1/\phi$ does not change its minimizer. Under the log link, the score equation for this minimization is $$ \sum_i x_i\left(1-\frac{y_i}{\mu_i}\right)=0. $$ The Hessian of $Q$ is $\sum_i y_i\exp(-x_i^T\beta)x_ix_i^T$, which is positive semidefinite and is positive definite when the weighted design spans the coefficient space. This establishes convexity, though existence of a finite minimizer must still be checked for the data and design. $\square$

A useful expression for the profile follows by measuring fit relative to the saturated means $\mu_i=y_i$. Define $$ D(\beta)=\sum_i\left\{\frac{y_i}{\mu_i}-\log\left(\frac{y_i}{\mu_i}\right)-1\right\}, \qquad \hat D=\min_{\beta\in\mathcal B}D(\beta). $$ Each summand is nonnegative: for $u>0$, $u-\log u-1\geq0$, with equality only at $u=1$. Substitution in the likelihood gives $$ \ell_p(\phi)=-n\log\Gamma(1/\phi)-\frac n\phi\log\phi -\sum_i\log y_i-\frac{n+\hat D}{\phi}. $$ The dependence on $\phi$ is not the Gaussian form. In particular, both the Gamma function and the factor $1/\phi$ contribute to the shape of the profile. A dispersion interval is found by maximizing this expression over $\phi>0$, then solving $2\{\ell_p(\hat\phi)-\ell_p(\phi)\}=\chi^2_{1,1-\alpha}$ at the crossings.

Worked Example: Profiling an Intercept in a Gamma GLM

Take two independent observations $y_1=1$ and $y_2=3$, with a Gamma log-link model containing only an intercept, so both means equal $\mu=\exp(\beta)$. At fixed $\phi$, the coefficient maximizes the likelihood by minimizing $$ Q(\beta)=2\log\mu+\frac{1+3}{\mu}=2\log\mu+\frac4\mu. $$ Differentiating with respect to $\mu>0$ gives $2/\mu-4/\mu^2$, which is zero at $\hat\mu=2$; the derivative is negative for $\mu<2$ and positive for $\mu>2$. Thus this is the minimum of $Q$, and the constrained coefficient estimate is $\hat\beta=\log 2$ for every fixed $\phi$.

At this fitted mean, the relative deviance quantity is $$ \hat D=\left\{\frac12-\log(1/2)-1\right\} +\left\{\frac32-\log(3/2)-1\right\} =\log 2-\log(3/2)=\log(4/3)\approx0.287682. $$ For these data, $\sum_i\log y_i=\log3$, so the profile log-likelihood is $$ \ell_p(\phi)=-2\log\Gamma(1/\phi)-\frac2\phi\log\phi-\log3 -\frac{2+\log(4/3)}{\phi}. $$ This formula fully specifies the one-dimensional profile. Its maximum and the two 95% endpoints can be obtained by numerical optimization and root-finding; at each candidate $\phi$, no new coefficient optimization is needed here because the intercept estimate is invariant to $\phi$. The example also shows why the Gamma profile should be evaluated from its full likelihood, rather than by applying the Gaussian RSS formula.

Endpoint Interpretation and a Shape-Scale Conversion

A profile interval is defined by a likelihood-ratio threshold, not by equal distances on either side of the estimate. The lower and upper endpoints solve the same equation, but the profile may rise or fall at different rates on the two sides. Consequently, a numerical interval can be strongly asymmetric even when it is connected and contains the maximum. Reporting the endpoints directly preserves this information; replacing them by an estimate plus or minus a symmetric standard error can obscure it.

Worked Example: Converting Gamma Shape Endpoints

Some Gamma software reports shape $\kappa=1/\phi$ rather than dispersion. Suppose a profile calculation has produced a shape interval $(0.50,11.7)$. Since $\phi=1/\kappa$ is strictly decreasing for positive $\kappa$, transform both endpoints and reverse their order: $$ \left(\frac1{11.7},\frac1{0.50}\right) =(0.085470085\ldots,2.00)\approx(0.08547,2.00). $$ The lower dispersion endpoint is the reciprocal of the upper shape endpoint. The order reversal is essential: taking reciprocals without reversing the endpoints would produce an incorrectly ordered interval.

The profile interval is useful when dispersion is itself scientifically relevant, such as when assessing the degree of extra variation relative to the mean-variance relationship. It also makes clear how uncertainty in the regression coefficients enters: the likelihood at each candidate dispersion is maximized over those coefficients before comparison with the global maximum. In the Gaussian and Gamma cases developed here, that constrained maximizer has a convenient dispersion-invariance property. The property should not be assumed for other likelihoods without checking the likelihood’s dependence on the regression parameters.

Common Mistakes

Holding the regression coefficients fixed. A dispersion profile uses the best coefficient fit at each candidate dispersion. In these two models the maximizing coefficient vector happens to be invariant to $\phi$, but this follows from their likelihoods; it is not the definition of profiling.
Using the Gaussian formula for Gamma data. The Gaussian profile depends on the residual sum of squares. The Gamma profile contains $\log\Gamma(1/\phi)$ and a different mean-fit term, so its endpoints must come from its own likelihood.
Assuming the interval is symmetric. The cutoff is applied to the profile likelihood-ratio statistic. There is no requirement that the endpoint distances from $\hat\phi$ match.
Reversing the Gamma shape-to-dispersion conversion incorrectly. Because $\phi=1/\kappa$ decreases with $\kappa$, reciprocal endpoints must be reordered.
Calling a likelihood-ratio interval exact. The chi-square cutoff is generally justified asymptotically under regularity conditions. A finite-sample exact interval requires a separate exact pivot or calibration argument.

Key takeaway: For a dispersion profile, maximize over regression coefficients at each fixed dispersion, then retain values whose likelihood-ratio statistic is below the chosen cutoff. Gaussian regression gives a closed-form profile; the Gamma profile has a distinct shape, and its endpoints can be notably asymmetric.

Check Your Understanding

Use the likelihood and parameterization appropriate to each model.

In Gaussian regression, why does the coefficient estimate that maximizes the likelihood at fixed $\phi$ not depend on $\phi$?
For a Gaussian profile, write the equation in $r=\phi/\hat\phi$ that determines the likelihood-ratio interval endpoints.
In the Gamma log-link model, identify the objective in $\beta$ that must be minimized at fixed dispersion.
Why does the Gamma coefficient maximizer remain the same as $\phi$ varies in the model considered here?
If a Gamma shape interval is $(a,b)$, where $0<a<b$, what is the correctly ordered dispersion interval?