Tweedie distributions

Tweedie distributions are exponential dispersion models, with a mean \(\mu\) and a variance \(\phi \mu^\xi\), for some dispersion parameter \(\phi > 0\) and a power index \(\xi\) (sometimes called \(p\)) that uniquely defines the distribution within the Tweedie family (for all real values of \(\xi\) not between 0 and 1).

Special cases of the Tweedie distributions are:

• the normal distribution, with \(\xi = 0\) (i.e., the variance is \(\phi\) and not related to the mean);
• the Poisson distribution, with \(\xi = 1\) and \(\phi = 1\) (i.e., the variance is the same as the mean);
• the gamma distribution, with \(\xi = 2\); and
• the inverse Gaussian distribution, with \(\xi = 3\).

For all other values of \(\xi\), the probability functions and distribution functions have no closed forms.

For \(\xi < 1\), applications are limited (non-existent so far?), but have support on the entire real line and \(\mu > 0\).

For \(1 < \xi < 2\), Tweedie distributions can be represented as a Poisson sum of gamma distributions. These distributions are continuous for \(Y > 0\) but have a discrete mass at \(Y = 0\).

Examples

Likelihood computations are not need to fit a Tweedie generalized linear model, using the tweedie() family from the statmod package:

library(statmod)

set.seed(96)
N <- 25
# Mean of the Poisson (lambda) and Gamma (shape/scale)
lambda <- 1.5
# Generating Compound Poisson-Gamma data manually
y <- replicate(N, {
  n_events <- rpois(1, lambda = lambda)
  if (n_events == 0) 0 else sum(rgamma(n_events, shape = 2, scale = 1))
})

mod.tw <- glm(y ~ 1, 
              family = statmod::tweedie(var.power = 1.5, link.power = 0) )
              # link.power = 0  means the log-link

However, likelihood computations are necessary in other situations.

tweedie::rtweedie(10, xi = 1.1, mu = 2, phi = 1)
#>  [1] 1.536587 2.520597 4.624115 1.999447 2.532301 4.259932 1.725366 0.000000
#>  [9] 0.000000 2.545023

y <- seq(0, 2, length = 100)
xi <- 1.1
mu <- 0.5
phi <- 0.4

twden <- tweedie::dtweedie(y, xi = xi, mu = mu, phi = phi)  
twdtn <- tweedie::ptweedie(y, xi = xi, mu = mu, phi = phi)

plot( twden[y > 0] ~ y[y > 0], 
      type ="l",
      lwd = 2,
      xlab = expression(italic(y)),
      ylab = "Density function")
points(twden[y==0] ~ y[y == 0],
      lwd = 2,
      pch = 19,
      xlab = expression(italic(y)),
      ylab = "Distribution function")

plot(twdtn ~ y,
     type = "l",
      lwd = 2,
     ylim = c(0, 1),
      xlab = expression(italic(y)),
      ylab = "Distribution function")

tweedie::tweedie_plot(y, xi = xi, mu = mu, phi = phi,
                      ylab = "Density function",
                      xlab = expression(italic(y)),
                      lwd = 2)

tweedie::tweedie_plot(y, xi = xi, mu = mu, phi = phi, 
                      ylab = "Distribution function",
                      xlab = expression(italic(y)),
                      lwd = 2, 
                      ylim = c(0, 1), 
                      type = "cdf")

library(tweedie)

qqnorm( statmod::qresid(mod.tw) )

Estimating \(\xi\)

To use a Tweedie distribution in a glm, the value of \(\xi\) is needed. To find the most suitable value of \(\xi\), tweedie_profile() can be used:

out <- tweedie::tweedie.profile(y ~ 1, 
                                xi.vec = seq(1.2, 1.8, by = 0.05), 
                                do.plot = TRUE)
#> Warning: `tweedie.profile()` was deprecated in tweedie 3.0.5.
#> ℹ Please use `tweedie_profile()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> .............Done.

# The estimated power index:
out$xi.max
#> [1] 1.432653