Estimating Inverse-Gaussian Dispersion

Why Estimate Dispersion?

In an inverse-Gaussian regression model, the fitted mean describes the response’s location, while the dispersion parameter describes its variability around that mean. Two models can give the same fitted means but very different uncertainty if their dispersion estimates differ. Estimating dispersion is therefore necessary for uncertainty calculations, even when the mean model is the main scientific target.

For independent responses $Y_1,\ldots,Y_n$ with positive means $\mu_1,\ldots,\mu_n$ and a common dispersion $\phi>0$, use the inverse-Gaussian parameterization

f(y_i;\mu_i,\phi) = \frac{1}{\sqrt{2\pi\phi y_i^3}} \exp\left\{-\frac{(y_i-\mu_i)^2}{2\phi\mu_i^2y_i}\right\}, \qquad y_i>0.

In this parameterization, $E[Y_i]=\mu_i$ and $\operatorname{Var}(Y_i)=\phi\mu_i^3$. In this tutorial, the fitted means are treated as given when deriving the conditional likelihood estimate of $\phi$. This distinction matters: the estimate conditional on the means is not the same formula as a residual-degrees-of-freedom adjustment.

Conditional Maximum Likelihood

Definition: Given positive fitted means $\hat\mu_1,\ldots,\hat\mu_n$, define the inverse-Gaussian likelihood residual sum $$ S_{\mathrm{IG}}=\sum_{i=1}^n\frac{(y_i-\hat\mu_i)^2}{\hat\mu_i^2y_i}. $$ The conditional maximum-likelihood estimate of the common dispersion is the value of $\phi$ that maximizes the likelihood with these means held fixed.

The denominator $\hat\mu_i^2y_i$ is specific to the inverse-Gaussian likelihood. It is not the denominator in the Pearson statistic, which comes from the variance function $\mu_i^3$.

Theorem: If $S_{\mathrm{IG}}>0$, the conditional maximum-likelihood estimate is $$ \hat\phi_{\mathrm{ML}}=\frac{S_{\mathrm{IG}}}{n}. $$ It is the unique maximizer over $\phi>0$.

Proof. Independence makes the joint likelihood the product of the individual densities. With the means held fixed, its log-likelihood can be written

\ell(\phi) = C-\frac{n}{2}\log\phi-\frac{S_{\mathrm{IG}}}{2\phi},

where $C$ contains terms that do not depend on $\phi$. Differentiating for $\phi>0$ gives

\ell'(\phi)=-\frac{n}{2\phi}+\frac{S_{\mathrm{IG}}}{2\phi^2} =\frac{S_{\mathrm{IG}}-n\phi}{2\phi^2}.

The derivative is positive when $0<\phi<S_{\mathrm{IG}}/n$ and negative when $\phi>S_{\mathrm{IG}}/n$. Thus the likelihood increases up to $S_{\mathrm{IG}}/n$ and decreases afterwards, so that point is the unique maximizer. Equivalently, $\ell''(S_{\mathrm{IG}}/n)=-n/(2\hat\phi_{\mathrm{ML}}^2)<0$. This proves the result. $\square$

If $S_{\mathrm{IG}}=0$, every observation equals its supplied mean. In that special case the likelihood increases without bound as $\phi$ decreases to zero, so there is no maximizer in the parameter space $\phi>0$. The displayed estimate is an interior maximum only when the residual sum is positive.

When the means come from a regression model, the conditional formula still applies after substituting the fitted means. For the ordinary joint likelihood, maximizing over $\phi$ at each candidate mean-model fit gives this same formula with those candidate means. This does not make the conditional estimate a degrees-of-freedom-adjusted estimate.

Worked Example: Calculating the Conditional Likelihood Estimate

Suppose three positive observations are $y=(1,5,3)$, and the supplied fitted means are $\hat\mu=(2,4,3)$. Calculate the likelihood residual sum and conditional maximum-likelihood estimate.

For the first observation, the contribution is $(1-2)^2/(2^2\cdot1)=1/4$. For the second, it is $(5-4)^2/(4^2\cdot5)=1/80=0.0125$. The third contribution is zero. Therefore,

S_{\mathrm{IG}}=\frac14+\frac1{80}=\frac{21}{80}=0.2625, \qquad \hat\phi_{\mathrm{ML}}=\frac{S_{\mathrm{IG}}}{3} =\frac7{80}=0.0875.

The result is positive, so it is an interior maximum. It estimates the dispersion conditional on these three supplied means; it is not adjusted for how those means may have been estimated.

Pearson-Based Dispersion Estimation

The Pearson statistic divides squared residuals by the model variance function without its dispersion multiplier. For inverse-Gaussian responses, define

X^2_{\mathrm{P}} = \sum_{i=1}^n \frac{(y_i-\hat\mu_i)^2}{\hat\mu_i^3}.

If the means are known rather than estimated, this statistic has expectation $n\phi$. Indeed, each term has expectation

E\left[\frac{(Y_i-\mu_i)^2}{\mu_i^3}\right] = \frac{\operatorname{Var}(Y_i)}{\mu_i^3} = \frac{\phi\mu_i^3}{\mu_i^3} = \phi,

and adding the $n$ expectations gives $E[X^2_{\mathrm{P}}]=n\phi$. Thus, with known means, $X^2_{\mathrm{P}}/n$ is unbiased for $\phi$.

In a regression, the means are estimated from the same observations. A common Pearson-based estimate is $X^2_{\mathrm{P}}/(n-p)$, where $p$ is the number of fitted mean-model coefficients and $n>p$. The residual degrees of freedom account approximately for the fitting of those coefficients; the adjustment is generally not an exact finite-sample unbiasedness result for a nonlinear GLM.

Result: Under a local linear approximation to a regular mean model, the Pearson statistic has expectation approximately $(n-p)\phi$, which motivates the estimate $$ \hat\phi_{\mathrm{P}}=\frac{X^2_{\mathrm{P}}}{n-p}. $$ This is an approximation when the means are estimated, not the conditional maximum-likelihood formula.

To see the degrees-of-freedom calculation, standardize the errors by their standard deviations, so the resulting vector $z$ has mean zero and covariance matrix $I$. Near the fitted values, the standardized fitted residuals are approximately $(I-H)z$, where $H$ is the symmetric idempotent projection matrix onto the $p$-dimensional fitted mean-model space. The Pearson statistic is then approximately $\phi z^\mathsf{T}(I-H)z$. Its expectation is

\phi E\!\left[z^\mathsf{T}(I-H)z\right] = \phi\,\operatorname{tr}\!\left((I-H)E[zz^\mathsf{T}]\right) = \phi\,\operatorname{tr}(I-H) = (n-p)\phi.

Here $E[zz^\mathsf{T}]=I$, and a rank-$p$ projection in $n$ dimensions has $\operatorname{tr}(I-H)=n-p$. The approximation enters through the local linearization of the fitted residuals, not through the trace calculation.

Worked Example: Comparing the Two Residual Sums

Return to the observations and means in the first example. The Pearson contributions are $(1-2)^2/2^3=1/8$, $(5-4)^2/4^3=1/64$, and $0$. Hence

X^2_{\mathrm{P}}=\frac18+\frac1{64}=\frac9{64}=0.140625.

If these means are treated as known, the Pearson estimate is $X^2_{\mathrm{P}}/3=3/64=0.046875$. This differs from the conditional likelihood estimate $0.0875$. If instead these means came from a fitted model with $p=1$, the usual residual-degrees-of-freedom estimate would be $X^2_{\mathrm{P}}/(3-1)=9/128=0.0703125$. These are three different calculations with three different interpretations; the degrees-of-freedom adjustment does not turn the Pearson estimate into the likelihood estimate.

The distinction between the residual sums is also visible algebraically:

\frac{(y_i-\mu_i)^2}{\mu_i^2y_i} \quad\text{versus}\quad \frac{(y_i-\mu_i)^2}{\mu_i^3}.

Their ratio, when $y_i\ne\mu_i$, is $\mu_i/y_i$. The contributions are close when the observation is near its fitted mean, but can differ substantially when it is not. Neither sum should be substituted for the other without changing the estimator being calculated.

Worked Example: How Dispersion Changes Uncertainty

Suppose a response has fitted mean $\mu=4$. Its variance is $\phi\mu^3=64\phi$. Compare two plausible dispersion values, $\phi=0.05$ and $\phi=0.20$.

At $\phi=0.05$, the variance is $64(0.05)=3.2$, giving response standard deviation $\sqrt{3.2}\approx1.789$. At $\phi=0.20$, the variance is $64(0.20)=12.8$, giving standard deviation $\sqrt{12.8}\approx3.578$. Multiplying dispersion by four doubles the standard deviation because standard deviation scales with $\sqrt{\phi}$.

For ten independent responses with this same mean and dispersion, their sample average has variance $\phi\mu^3/10$. Its standard deviation is therefore $\sqrt{3.2/10}\approx0.566$ at the first dispersion value and $\sqrt{12.8/10}\approx1.131$ at the second. The uncertainty changes even though the fitted mean and sample size have not changed.

Dispersion and Regression Uncertainty

The same scaling carries over to coefficient uncertainty. For a regular inverse-Gaussian GLM with design matrix $X$, linear predictors $\eta_i$, and a differentiable link inverse $\mu_i=\mu(\eta_i)$, the Fisher information for the mean-model coefficients is proportional to $1/\phi$. In particular, its usual large-sample covariance approximation has the form

\operatorname{Cov}(\hat\beta) \approx \phi\left[ X^\mathsf{T} \operatorname{diag}\left\{ \frac{(d\mu_i/d\eta_i)^2}{\mu_i^3} \right\} X \right]^{-1},

provided the matrix being inverted is nonsingular and the model is evaluated at the fitted values. Consequently, replacing one positive dispersion estimate by another scales the estimated covariance matrix by their ratio and the coefficient standard errors by the square root of that ratio. The inverse-Gaussian Fisher scoring weights determine the mean-model part of this expression; dispersion multiplies the resulting covariance.

Common Mistakes

Using the Pearson denominator in the likelihood estimate. The likelihood contribution has denominator $\mu_i^2y_i$, whereas the Pearson contribution has denominator $\mu_i^3$. They are not generally equal.
Dividing the conditional MLE by $n-p$. The exact conditional likelihood maximizer is $S_{\mathrm{IG}}/n$. The Pearson residual-degrees-of-freedom estimate is a different estimator.
Calling the Pearson adjustment exactly unbiased after fitting. The $n-p$ argument uses a local linear approximation to fitted residuals. It is a useful rationale, not a universal exact identity for finite-sample nonlinear regression.
Thinking dispersion changes the fitted means. Dispersion controls the variance scale. A change in $\phi$ directly changes uncertainty, but does not by itself change a supplied fitted mean.

Key takeaway: For fixed fitted means, the inverse-Gaussian conditional MLE is the likelihood residual sum divided by $n$. A Pearson estimate instead uses variance-standardized residuals and often divides by $n-p$; its interpretation and residual sum are different.

Check Your Understanding

Use the definitions and distinctions above to answer the following.

Write the likelihood residual contribution for one observation and the Pearson contribution for the same observation.
Why does the conditional maximum-likelihood estimate divide its residual sum by $n$, rather than $n-p$?
For known means, what is $E[X^2_{\mathrm{P}}]$, and which property of the inverse-Gaussian model gives this result?
What approximation motivates dividing the fitted Pearson statistic by $n-p$?
If $\phi$ is multiplied by nine while the fitted mean is unchanged, by what factor does the response standard deviation change?