The GENMOD Procedure

Tweedie Distribution For Generalized Linear Models

The Tweedie (1984) distribution has nonnegative support and can have a discrete mass at zero, making it useful to model responses that are a mixture of zeros and positive values. The Tweedie distribution belongs to the exponential family, so it conveniently fits in the generalized linear models framework. According to such parameterization, the mean and variance for the Tweedie random variable are $\mr{E}(Y)=\mu $ and $\mr{Var}(Y)=\phi \mu ^ p$, respectively, where $\phi $ is the dispersion parameter and p is an extra parameter that controls the variance of the distribution.

The Tweedie family of distributions includes several important distributions for generalized linear models. When $p=0$, the Tweedie distribution degenerates to the normal distribution; when $p=1$, it becomes a Poisson distribution; when $p=2$, it becomes a gamma distribution; when $p=3$, it is an inverse Gaussian distribution.

Except for these special cases, the probability density function for the Tweedie distribution does not have a closed form and can at best be expressed in terms of series. Numerical approximations are needed to evaluate the density function. Dunn and Smyth (2005) propose using a finite series and provide a formula to determine its lower and upper indices in order to achieve a desired accuracy. Alternatively, you can apply the Fourier transformation on the characteristic function (Dunn and Smyth, 2008). These approximations tend to be expensive when a high level of accuracy is demanded or the data volume becomes large. PROC GENMOD uses the series method unless it becomes complicated to do so. In this case, the method that is based on the Fourier transformation is used. The accuracy of approximation is controlled by the EPSILON= option, whose default value is $10^{-5}$.

The Tweedie distribution is not defined when p is between 0 and 1. In practice, the most interesting range is from 1 to 2 in which the Tweedie distribution gradually loses its mass at 0 as it shifts from a Poisson distribution to a gamma distribution. In this case, the Tweedie random variable $Y$ can be generated from a compound Poisson distribution (Smyth, 1996) as

\begin{eqnarray*} Y & =& \Sigma _{i=1}^{T} X_ i \\ T & \sim & \mr{Poisson}(\lambda ) \\ X_ i & \sim & \mr{gamma}(\alpha ,\gamma ) \end{eqnarray*}

where $Y=0$ if $T=0$, $T$ and $X_ i$ are statistically independent, and $\mr{gamma}(\alpha , \gamma )$ denotes a gamma random variable that has mean $\alpha \gamma $ and variance $\alpha \gamma ^2$. These parameters are determined by the Tweedie parameters as follows:

\begin{eqnarray*} \lambda & =& \frac{\mu ^{2-p}}{\phi (2-p)} \\ \alpha & =& \frac{2-p}{p-1} \\ \gamma & =& \phi (p-1) \mu ^{p-1} \end{eqnarray*}

Inversely, given the Tweedie distributional parameters, the parameters of the compound Poisson distribution are determined as follows:

\begin{eqnarray*} \mu & =& \lambda \alpha \gamma \\ p & =& \frac{\alpha +2}{\alpha +1} \\ \phi & =& \frac{\lambda ^{1-p} (\alpha \gamma )^{2-p}}{2-p} \end{eqnarray*}

In terms of generalized linear models parameterizations, the canonical parameter $\theta $ for the Tweedie density can be expressed as

\begin{eqnarray*} \theta & = & \left\{ \begin{array}{ll} \frac{\mu ^{1-p}}{1-p} & p \neq 1 \\ \log \mu & p = 1 \\ \end{array} \right. \\ \end{eqnarray*}

and the function $b(\theta )$ is

\begin{eqnarray*} b(\theta ) & = & \left\{ \begin{array}{ll} \frac{\mu ^{2-p}}{2-p} & p \neq 2 \\ \log \mu & p = 2 \\ \end{array} \right. \\ \end{eqnarray*}

Because of the intractability of differentiating the gradient functions with respect to the variance parameters, PROC GENMOD uses a quasi-Newton approach to maximize the likelihood function, where the Hessian matrix is approximated by taking finite differences of the gradient functions. Convergence is determined by a union of two criteria: the relative gradient convergence criterion is set to $10^{-9}$, and the relative function convergence criterion is set to $2\times 10^{-9}$. Convergence is declared when at least one of the criteria is attained during the quasi-Newton iteration.

Before PROC GENMOD maximizes the approximate likelihood, it first maximizes the following extended log quasi-likelihood which is constructed according to the definition of McCullagh and Nelder (1989, Chapter 9) as

\[ Q_ p(\mb{y}, \bmu , \phi , p) = \sum _ i q(y_ i, \mu _ i, \phi , p) \]

where the contribution from an observation is

\[ q(y_ i, \mu _ i, \phi , p) = -0.5 \log (2 \pi \phi y_ i^ p / w_ i) - w_ i \left( \frac{y_ i^{2-p}-(2-p)y_ i \mu _ i^{1-p} + (1-p) \mu _ i^{2-p}}{(1-p)(1-p)} \right) / \phi \]

and $w_ i$ is the weight for the observation from the WEIGHT statement.

The range of parameter p for the quasi-likelihood is from 1 to 2. For a specified P= value outside this range, PROC GENMOD skips optimization of the quasi-likelihood. To maintain numerical stability, PROC GENMOD imposes a lower bound of 1.1 and a upper bound of 1.99 for computation with the quasi-likelihood. The estimates that are obtained from optimizing the quasi-likelihood are usually near the full-likelihood solution so that fewer iterations are needed for maximizing the more expensive full likelihood.