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
variance-power index \(\xi\) (Tweedie 1984), that is sometimes called \(p\), which 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.
| \(\xi < 0\) |
Extreme stable |
\((-\infty,
\infty)\) |
| \(\xi = 0\) |
Normal |
\((-\infty,
\infty)\) |
| \(0 < \xi <
1\) |
Do not exist |
|
| \(\xi = 1\) (with
\(\phi = 1\)) |
Poisson |
Discrete: \(\{0, 1, 2,
\dots\}\) |
| \(1 < \xi <
2\) |
Poisson–gamma |
\([ 0,
\infty)\) |
| \(\xi = 2\) |
Gamma |
\((0, \infty)\) |
| \(\xi > 2\) |
Skewed right |
\((0, \infty)\) |
| \(\xi = 3\) |
inverse Gaussian |
\((0, \infty)\) |
For \(\xi < 0\), applications are
limited, 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\).
Use cases
Likelihood computations are not needed to fit a Tweedie generalized
linear model, using the tweedie() family from the statmod
package Dunn and Smyth (2018):
library(tweedie)
library(statmod) # For the tweedie() family function, for use in glm()
set.seed(96)
N <- 25
# Poisson mean (lambda): the gamma summands have shape 2 and scale 1
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.
Since likelihoods cannot be computed analytically apart from special
cases, the tweedie package computes Tweedie density
functions and distribution functions using numerical methods Dunn and Smyth (2008).
Computing the quantile residuals
Quantile residuals (Dunn and Smyth
1996) can be computed to assess the fit of a glm. The function
qresid() is in the statmod
package Dunn and Smyth (2018):
qqnorm( qr <- statmod::qresid(mod.tw),
las = 1)
qqline(qr,
col = "grey")
The quantile residuals suggest the model may not be adequate.
Generating random numbers
Random numbers from a Tweedie distribution are generated using
rtweedie():
tweedie::rtweedie(10, xi = 1.4, mu = 2, phi = 1)
#> [1] 1.15415709 0.00000000 0.68100935 3.78377438 2.45493862 1.04947222
#> [7] 1.57083901 0.09651273 7.47423368 1.89911729
Plotting density and probability functions
Accurate density functions are generated using
dtweedie():
y_grid <- seq(0, 2,
length = 100)
xi <- 1.05
mu <- 0.5
phi <- 0.4
twden <- tweedie::dtweedie(y_grid, xi = xi, mu = mu, phi = phi)
twdtn <- tweedie::ptweedie(y_grid, xi = xi, mu = mu, phi = phi)
plot( twden[y_grid > 0] ~ y_grid[y_grid > 0],
xlab = expression(italic(y)),
ylab = "Density function",
type ="l",
lwd = 2,
las = 1)
points(twden[y_grid==0] ~ y_grid[y_grid == 0],
lwd = 2,
pch = 19)
The solid point at \(y = 0\) shows
the probability mass \(\Pr(Y = 0)\),
not a density. Accurate probability functions are generated using
ptweedie():
plot(twdtn ~ y_grid,
xlab = expression(italic(y)),
ylab = "Distribution function",
type = "l",
las = 1,
lwd = 2,
ylim = c(0, 1) )
points(twdtn[y_grid==0] ~ y_grid[y_grid == 0],
lwd = 2,
pch = 19)
However, the function tweedie_plot() is sometimes
more convenient (especially for \(1 < \xi
< 2\), when a probability mass at \(Y = 0\) is present):
tweedie::tweedie_plot(y_grid, xi = xi, mu = mu, phi = phi,
plot_args = list(
ylab = "Density function",
xlab = expression(italic(y)),
las = 1,
lwd = 2
)
)
tweedie::tweedie_plot(y_grid, xi = xi, mu = mu, phi = phi,
type = "cdf",
plot_args = list(
ylab = "Distribution function",
xlab = expression(italic(y)),
las = 1,
lwd = 2,
ylim = c(0, 1))
)
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 (Dunn and Smyth 2018):
y <- rtweedie(50, xi = 1.5, mu = 1, phi = 1)
out <- tweedie_profile(y ~ 1,
xi.vec = seq(1.3, 1.8, length = 5),
do.plot = TRUE)
#> ================================================================================================================================================================================================================================================
# The estimated power index:
out$xi.max
#> [1] 1.581986
# The 95% CI for the power index:
out$ci
#> [1] 1.400437 1.748554
Example: \(p > 2\)
Consider the poison data in the GLMsData
package (Dunn and Smyth 2022), which
records the survival times of animals under various treatments and
poisons (used in (Dunn and Smyth 2018),
Sect. 12.4.2). For convenience, a copy is included with this
package:
poison <- read.csv(system.file("extdata",
"poison.csv",
package = "tweedie"))
# Convert to factors:
poison$Psn <- factor(poison$Psn)
poison$Trmt <- factor(poison$Trmt)
head(poison)
#> Psn Trmt Time
#> 1 I A 0.31
#> 2 I A 0.45
#> 3 I A 0.46
#> 4 I A 0.43
#> 5 I B 0.82
#> 6 I B 1.10
The survival time is likely related to the poison type and type of
treatment:
plot(Time ~ Psn,
data = poison,
xlab = "Poison type",
ylab = "Survival time (tens of hours)",
las = 1)
plot(Time ~ Trmt,
data = poison,
xlab = "Treatment type",
ylab = "Survival time (tens of hours)",
las = 1)
In both cases, the variance of the survival times seems to increase
with the mean survival time. Perhaps a Tweedie glm would be appropriate.
To fit a Tweedie glm, the value of \(\xi\) is needed; since no exact zeros are
present (or possible) for the survival times, very likely \(\xi \ge 2\):
PSNPrf <- tweedie_profile(Time ~ Trmt + Psn,
data = poison,
do.plot = TRUE,
xi.vec = seq(2.5, 5.5, length = 7) )
#> ================================================================================================================================================================================================================================================================================================================================
Then, \(\xi\approx{}\) 4.1:
PSNPrf$xi.max
#> [1] 4.117514
So then the Tweedie glm can be fitted:
PSNMod <- glm(Time ~ Trmt + Psn,
data = poison,
family = statmod::tweedie(var.power = round(PSNPrf$xi.max, 1),
link.power = 0))
anova(PSNMod,
test = "F")
#> Analysis of Deviance Table
#>
#> Model: Tweedie, link: mu^0
#>
#> Response: Time
#>
#> Terms added sequentially (first to last)
#>
#>
#> Df Deviance Resid. Df Resid. Dev F Pr(>F)
#> NULL 47 68.294
#> Trmt 3 21.275 44 47.019 24.531 2.486e-09 ***
#> Psn 2 35.672 42 11.346 61.698 3.157e-13 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Quantile residuals can be used to assess the model:
qr <- statmod::qresid(PSNMod)
qqnorm(qr,
las = 1)
qqline(qr,
col = "grey")
The quantile residuals suggest that this model looks fine.
Example: \(1 < p < 2\)
Consider the quilpie data in the GLMsData
package (Dunn and Smyth 2022), which
records the July rainfalls at Quilpie (used in (Dunn and Smyth 2018)). For convenience, a copy
is included with this package:
quilpie <- read.csv(system.file("extdata",
"quilpie.csv",
package = "tweedie"))
# Convert to factors:
quilpie$Phase <- factor(quilpie$Phase)
head(quilpie)
#> X Year Rain SOI Phase Exceed y
#> 1 1 1921 38.4 2.7 2 Yes 1
#> 2 2 1922 0.0 2.0 5 No 0
#> 3 3 1923 0.0 -10.7 3 No 0
#> 4 4 1924 24.4 6.9 2 Yes 1
#> 5 5 1925 0.0 -12.5 3 No 0
#> 6 6 1926 9.1 -1.0 4 No 0
The rainfall may be related to the SOI Phase:
plot(Rain ~ Phase,
data = quilpie,
xlab = "SOI Phase",
ylab = "July rainfall (mm)",
las = 1)
The variance of the survival times seems to increase with the mean
rainfall: Also, some rainfall totals are zero:
cat("Percentage of Julys with zero rain:",
round( sum(quilpie$Rain == 0) / length(quilpie$Rain) * 100, 1),
"\n")
#> Percentage of Julys with zero rain: 29.4
Perhaps a Tweedie glm would be appropriate. To fit a Tweedie glm, the
value of \(\xi\) is needed, so \(1 < \xi < 2\):
QRainPrf <- tweedie_profile(Rain ~ Phase,
data = quilpie,
do.plot = TRUE,
xi.vec = seq(1.2, 1.7, length = 7) )
#> ================================================================================================================================================================================================================================================================================================================================
Then, \(\xi\approx{}\) 1.4:
QRainPrf$xi.max
#> [1] 1.372874
So then the Tweedie glm can be fitted:
QRainMod <- glm(Rain ~ Phase,
data = quilpie,
family = statmod::tweedie(var.power = round(QRainPrf$xi.max, 1),
link.power = 0))
anova(QRainMod,
test = "F")
#> Analysis of Deviance Table
#>
#> Model: Tweedie, link: mu^0
#>
#> Response: Rain
#>
#> Terms added sequentially (first to last)
#>
#>
#> Df Deviance Resid. Df Resid. Dev F Pr(>F)
#> NULL 67 617.44
#> Phase 4 225.83 63 391.61 8.8803 9.439e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Quantile residuals can be used to assess the model:
qr <- statmod::qresid(QRainMod)
qqnorm(qr,
las = 1)
qqline(qr,
col = "grey")
The quantile residuals suggest that this model looks fine.
Since \(1 < p < 2\), the glm
parameters can be converted to the underlying parameters:
out <- tweedie_convert(xi = QRainPrf$xi.max,
mu = fitted(QRainMod),
phi = summary(QRainMod)$dispersion )
out2 <- cbind(
Year = quilpie$Year,
lambda = round(out$poisson.lambda, 2), # Poisson mean number of rain events
shape = round(out$gamma.shape, 2), # Gamma shape
p0 = round(out$p0, 2), # Pr(Rainfall = 0) = exp( -lambda)
scale = round(out$gamma.scale, 2) # Gamma scale
)
head(out2, 10)
#> Year lambda shape p0 scale
#> 1 1921 2.29 1.68 0.10 8.82
#> 2 1922 1.19 1.68 0.31 5.97
#> 3 1923 0.58 1.68 0.56 3.92
#> 4 1924 2.29 1.68 0.10 8.82
#> 5 1925 0.58 1.68 0.56 3.92
#> 6 1926 1.51 1.68 0.22 6.88
#> 7 1927 1.19 1.68 0.31 5.97
#> 8 1928 1.51 1.68 0.22 6.88
#> 9 1929 1.19 1.68 0.31 5.97
#> 10 1930 1.19 1.68 0.31 5.97
References
Dunn, Peter K., and Gordon K. Smyth. 1996. “Randomized Quantile
Residuals.” Journal of Computational and Graphical
Statistics 5 (3): 236–44.
Dunn, Peter K., and Gordon K. Smyth. 2005. “Series Evaluation of
Tweedie Exponential Dispersion Model Densities.”
Statistics and Computing 15 (4): 267–80.
Dunn, Peter K., and Gordon K. Smyth. 2008. “Evaluation of
Tweedie Exponential Dispersion Model Densities by
Fourier Inversion.” Statistics and
Computing 18 (1): 73–86.
Dunn, Peter K., and Gordon K. Smyth. 2018. Generalized Linear Models
with Examples in R. Springer.
Dunn, Peter K., and Gordon K. Smyth. 2022.
GLMsData: Generalized
Linear Model Data Sets.
https://doi.org/10.32614/CRAN.package.GLMsData.
Smyth, Gordon K., Lizhong Chen, Yifang Hu, Peter K. Dunn, Belinda
Phipson, and Yunshun Chen. 2025. Statmod: Statistical
Modeling.
Tweedie, M. C. K. 1984. “An Index Which Distinguishes Between Some
Important Exponential Families.” Statistics: Applications and
New Directions, 579–604.