Estimating Dispersion at a Fixed Power

Why Dispersion Estimation Matters

In a Tweedie model, the power $p$ determines how the variance changes with the mean, while the dispersion parameter $\phi$ sets its overall scale. Even when $p$ is treated as fixed, the data do not generally reveal $\phi$ directly. It must be estimated, and the estimate affects more than a description of response variability: it also scales the estimated uncertainty of regression coefficients.

There are three common approaches. Maximum likelihood uses the full Tweedie probability model. The Pearson approach uses squared residuals standardized by the fitted variance function. The deviance approach uses the fitted model's unit deviance. The latter two are convenient residual-based estimates, but they are not generally identical to the likelihood estimate or to each other.

Definition: For independent observations $Y_i$ with fitted means $\mu_i$, a Tweedie model with fixed power $p$ has $\operatorname{Var}(Y_i)=\phi V(\mu_i)=\phi\mu_i^p$. Let $r$ be the rank of the regression design matrix and let $\hat\mu_i$ be the fitted means. Define the unscaled Pearson statistic and residual deviance by $$ X_P^2=\sum_{i=1}^n \frac{(y_i-\hat\mu_i)^2}{V(\hat\mu_i)} \qquad\text{and}\qquad D=\sum_{i=1}^n d_p(y_i,\hat\mu_i), $$ where, for $p\ne 1,2$, the Tweedie unit deviance is $$ d_p(y,\mu)=2\left\{\frac{y^{2-p}}{(1-p)(2-p)} -\frac{y\mu^{1-p}}{1-p} +\frac{\mu^{2-p}}{2-p}\right\}. $$ At boundary values of $p$, the unit deviance is defined by its continuous limit. For an ordinary residual-based estimate, the Pearson and deviance estimates are $\hat\phi_P=X_P^2/(n-r)$ and $\hat\phi_D=D/(n-r)$, provided $n>r$.

Three Routes to an Estimate

The likelihood estimate uses the exact probability mass or density at each observed response. For fixed $p$, write the log-likelihood as $\ell(\beta,\phi;p)$, where $\beta$ parameterizes the mean model. The profile log-likelihood for dispersion is $$ \ell_{\mathrm{prof}}(\phi)=\sup_{\beta}\ell(\beta,\phi;p),\qquad \phi>0. $$ The maximum-likelihood estimate $\hat\phi_L$ is a maximizer of this function. Except in special cases, there is no general closed-form formula for it in a Tweedie model. It is found by numerical optimization using the Tweedie likelihood, evaluated exactly or with a stated approximation.

The Pearson estimate instead matches the observed standardized residual variation to its usual residual degrees of freedom. The deviance estimate does the same using a likelihood-based measure of lack of fit. Both are practical scale estimates, but $n-r$ is a degrees-of-freedom adjustment, not a claim that either statistic has an exact chi-squared distribution for every Tweedie model. Their finite-sample behavior can differ from that of the likelihood estimate.

Worked Example: Pearson and Deviance Estimates From Residual Summaries

Suppose a fixed-power Tweedie fit uses $n=8$ observations and a design matrix of rank $r=2$. Its fitted residual summaries are $X_P^2=12.8$ and $D=10.4$. Compute the two residual-based dispersion estimates.

The residual degrees of freedom are $n-r=8-2=6$. Therefore $$ \hat\phi_P=\frac{12.8}{6}=2.1333\ldots,\qquad \hat\phi_D=\frac{10.4}{6}=1.7333\ldots. $$

The Pearson estimate is larger because the observed Pearson statistic is larger than the deviance. This difference is not an arithmetic inconsistency: the two statistics summarize residual departures on different scales. The fit and its power are held fixed in this comparison; only the method of estimating $\phi$ changes.

Likelihood Estimation: A Gamma Special Case

For a general fixed Tweedie power, the likelihood depends on the full Tweedie density, including its dispersion-dependent normalizing terms. Consequently, maximizing only a residual statistic is not equivalent to maximizing likelihood. The Gamma case, $p=2$, illustrates why the likelihood estimate can require solving a nonlinear equation even when the means are held fixed.

Worked Example: Gamma Likelihood Estimate With Fixed Means

Take three positive observations whose ratios to their fixed means are $y_i/\mu_i=1/2,1,3/2$. In the Gamma model with variance $\phi\mu_i^2$, the shape parameter is $\alpha=1/\phi$. Up to terms not involving $\alpha$, the log-likelihood as a function of $\alpha$ is $$ \ell(\alpha)=\sum_{i=1}^3\left\{\alpha\log(y_i/\mu_i)-\alpha(y_i/\mu_i)+\alpha\log\alpha-\log\Gamma(\alpha)\right\}+C, $$ where $C$ does not depend on $\alpha$.

Differentiating and writing $\psi(\alpha)=d\log\Gamma(\alpha)/d\alpha$ gives the likelihood equation $$ \ell'(\alpha)=\sum_{i=1}^3\log(y_i/\mu_i)-\sum_{i=1}^3(y_i/\mu_i) +3\{\log\alpha+1-\psi(\alpha)\}=0. $$

Here $\sum_i\log(y_i/\mu_i)=\log(3/4)$ and $\sum_i y_i/\mu_i=3$. Thus the equation is $$ \log(3/4)-3+3\{\log\alpha+1-\psi(\alpha)\}=0. $$

Equivalently, $\log\alpha+1-\psi(\alpha)=1-\log(3/4)/3$, approximately $1.0959$. Numerical solution gives $\hat\alpha\approx5.38$, so the likelihood dispersion estimate is $\hat\phi_L=1/\hat\alpha\approx0.186$. This example uses fixed means to isolate the dispersion likelihood; in a regression fit, the likelihood is profiled over the mean-model coefficients as well.

How Dispersion Scales Coefficient Uncertainty

Suppose the mean is related to covariates through a differentiable link $g(\mu_i)=x_i^\mathsf{T}\beta$. Let $X$ be the design matrix, and define the diagonal matrix $W$ by $$ W_{ii}=\frac{1}{[g'(\mu_i)]^2V(\mu_i)}. $$ For independent observations from a regular exponential-dispersion model, the expected information for $\beta$, conditional on the means, is $X^\mathsf{T}WX/\phi$. The following derivation shows both why this is the information and why dispersion multiplies the coefficient covariance.

Theorem: Suppose the observations are independent, have an exponential-dispersion model with differentiable cumulant function and variance $\phi V(\mu_i)$, and have means linked by $g(\mu_i)=x_i^\mathsf{T}\beta$. Assume the usual regularity conditions for differentiating the likelihood and using its information matrix, and assume $X$ has full column rank. Then the expected information for $\beta$ is $X^\mathsf{T}WX/\phi$. Under the usual large-sample likelihood approximation, the covariance of $\hat\beta$ is estimated by $\phi(X^\mathsf{T}WX)^{-1}$, with $W$ evaluated at the fitted means.

Write the mean-parameter score.
For an exponential-dispersion model, differentiating the log-likelihood with respect to $\mu_i$ gives $(Y_i-\mu_i)/\{\phi V(\mu_i)\}$. Since $\partial\mu_i/\partial\beta=x_i/g'(\mu_i)$, the score contribution is $x_i(Y_i-\mu_i)/\{\phi V(\mu_i)g'(\mu_i)\}$.

Compute its covariance.
The observations are independent and $E[Y_i-\mu_i]=0$, so cross-observation covariances vanish. The $i$th contribution has covariance $x_ix_i^\mathsf{T}\operatorname{Var}(Y_i)/\{\phi^2 V(\mu_i)^2[g'(\mu_i)]^2\}=x_ix_i^\mathsf{T}/\{\phi V(\mu_i)[g'(\mu_i)]^2\}$.

Sum and invert.
Summing the contributions yields $X^\mathsf{T}WX/\phi$. Under the stated likelihood regularity conditions, the inverse expected information is the first-order large-sample covariance, namely $\phi(X^\mathsf{T}WX)^{-1}$. $\square$

In practice, $\phi$ is unknown, so a fitted covariance matrix substitutes an estimate of dispersion. This gives a useful direct comparison: using $\hat\phi_P$, $\hat\phi_D$, or $\hat\phi_L$ multiplies the same estimated inverse information by a different scalar, provided the fitted means and power are held fixed. Standard errors therefore scale as the square root of the chosen dispersion estimate. This usual calculation treats the estimated dispersion as a plug-in quantity; it does not separately account for uncertainty in that estimate.

Worked Example: Recomputing Standard Errors With Correct Information Weights

Consider a log-link model with fixed power $p=1.5$, design rows $(1,0),(1,1),(1,2)$, and fitted means $1,2,4$. For the log link, $g'(\mu)=1/\mu$, so $W_{ii}=\mu_i^{2-p}=\sqrt{\mu_i}$. The weights are therefore $1,\sqrt{2},2$.

The weighted information matrix without the factor $1/\phi$ is $$ X^\mathsf{T}WX= \begin{pmatrix} 1+\sqrt{2}+2 & \sqrt{2}+4\\ \sqrt{2}+4 & \sqrt{2}+8 \end{pmatrix} = \begin{pmatrix} 3+\sqrt{2} & 4+\sqrt{2}\\ 4+\sqrt{2} & 8+\sqrt{2} \end{pmatrix}. $$ Its determinant is $$ (3+\sqrt{2})(8+\sqrt{2})-(4+\sqrt{2})^2=8+3\sqrt{2}. $$ Thus $$ (X^\mathsf{T}WX)^{-1} =\frac{1}{8+3\sqrt{2}} \begin{pmatrix} 8+\sqrt{2} & -(4+\sqrt{2})\\ -(4+\sqrt{2}) & 3+\sqrt{2} \end{pmatrix}. $$

The diagonal entries are $(8+\sqrt{2})/(8+3\sqrt{2})\approx0.7690$ and $(3+\sqrt{2})/(8+3\sqrt{2})\approx0.3606$. Using the Pearson estimate from the first example, $\hat\phi_P=2.1333$, gives coefficient variances $2.1333(0.7689)\approx1.640$ and $2.1333(0.3605)\approx0.769$. The corresponding standard errors are approximately $\sqrt{1.640}=1.281$ and $\sqrt{0.769}=0.877$.

Using instead $\hat\phi_D=1.7333$, the variances are $1.7333(0.7689)\approx1.333$ and $1.7333(0.3605)\approx0.625$, giving standard errors approximately $1.155$ and $0.791$. The means, design, and information weights have not changed; the difference in standard errors comes entirely from the dispersion estimate.

Choosing and Interpreting the Estimate

Likelihood estimation is the direct choice when the Tweedie likelihood is available and the model is being fit by maximum likelihood. It uses more of the distributional information than either residual statistic, but can involve numerical optimization and depends on accurate likelihood evaluation. Pearson and deviance estimates are straightforward to calculate after fitting and are often useful as scale summaries or checks. A substantial gap between them can signal that the fit's residual behavior deserves closer inspection; it does not, on its own, identify the cause.

For a correctly specified model, the residual-statistic estimates are commonly motivated by the fact that fitted residuals have fewer effective degrees of freedom than the original observations. In a Gaussian linear model, the Pearson statistic equals the residual sum of squares divided by the variance, and its expectation is exactly $n-r$. For non-Gaussian Tweedie models, the analogous degrees-of-freedom adjustment is a useful estimation convention, not a universal exact finite-sample identity.

Common pitfall: It is tempting to call $X_P^2/(n-r)$ the “likelihood estimate” because it is a familiar estimate of residual scale. That is incorrect in general: it is the Pearson estimate. The likelihood estimate maximizes the full dispersion-dependent likelihood, while the deviance estimate divides a different residual statistic by the residual degrees of freedom. Label the method used, and remember that whichever estimate is used in the coefficient covariance scales its standard errors by its square root.

The power and the dispersion play distinct roles in the variance function: changing $p$ changes its dependence on $\mu_i$, whereas changing $\phi$ rescales it. This tutorial holds the power fixed, so the comparison concerns only how that common variance scale is estimated and propagated to coefficient uncertainty.

Check Your Understanding

Use these questions to check both the calculations and their interpretation.

A model has $n=20$, design rank $r=4$, and Pearson statistic $X_P^2=24$. What is the Pearson dispersion estimate?
Why does a likelihood estimate of Tweedie dispersion generally require more than dividing a residual statistic by $n-r$?
If the inverse weighted information has a diagonal entry of $0.5$, what coefficient standard error results from using $\hat\phi=2$?
For fixed fitted means and power, what changes in the coefficient covariance when the deviance estimate is substituted for the Pearson estimate?
Is $D/(n-r)$ guaranteed to equal the maximum-likelihood estimate in a non-Gaussian Tweedie model? Explain briefly.