Tweedie GLMs: Theory, Methods & Applications

• Tweedie GLMs are regression models defined by a power variance function, uniting Gaussian, Poisson, Gamma, and compound Poisson-Gamma distributions.
• They employ maximum likelihood estimation via IRLS or Newton methods, linking likelihood inference with β-divergence minimization.
• Extensions include double GLMs, smoothing penalties, and boosting to jointly model mean and dispersion for complex semicontinuous data.

A Tweedie generalized linear model (GLM) is a regression approach within the exponential dispersion family characterized by a power variance function $V(\mu) = \mu^p$, $p\in\mathbb{R}\setminus(0,1)$. Tweedie @@@@10@@@@ unite standard models such as Gaussian ($p=0$), Poisson ($p=1$), Gamma ($p=2$), inverse Gaussian ($p=3$), and compound Poisson-Gamma ($1maximum likelihood estimation for Tweedie models is equivalent to minimizing a $\beta$-divergence ($\beta = 2-p$), establishing an explicit connection between likelihood-based inference and divergence minimization (Yilmaz et al., 2012). Estimation leverages iteratively reweighted least squares (IRLS) or Newton-type methods; extensions include double GLMs (DGLMs) for joint mean and dispersion modeling, spatial regularization, Bayesian variable selection, and scalable boosting algorithms.

1. Exponential Dispersion Tweedie Family: Core Characterization

The Tweedie subclass of exponential dispersion models (EDMs) is defined by the density

$$p(y;\,\theta,\varphi) = a(y,\varphi) \exp\left\{\varphi^{-1} \left[y\theta - \psi(\theta)\right]\right\}$$

with canonical parameter $\theta$, cumulant $\psi$, mean $\mu = \psi'(\theta)$, variance $\varphi \psi''(\theta) = \varphi V(\mu)$. For the Tweedie family, $V(\mu) = \mu^p$, $p\in\mathbb{R}\setminus(0,1)$, with canonical link $\theta(\mu) = \frac{\mu^{1-p}-1}{1-p}$ and cumulant $$\psi(\theta) = \frac{[1+(1-p)\theta]^{(2-p)/(1-p)}-1}{2-p}$$ (for $p\neq 2$) (Yilmaz et al., 2012). The power index $p$ interpolates among classical families:

$p$	Model
0	Gaussian
1	Poisson
$1 < p < 2$	Compound Poisson-Gamma
2	Gamma
3	Inverse Gaussian

Formally, $\mathrm{Var}[Y]=\varphi\mu^p$.

2. Mean–Variance Power Law and Compound Poisson-Gamma Representation

The Tweedie mean-variance relationship is governed by $V(\mu)=\mu^p$. For $1$Y=\sum_{m=1}^M C_m$ with $M\sim\mathrm{Poisson}(\lambda)$ and $C_m\sim\mathrm{Gamma}(\alpha,\beta)$, leading to a discrete-continuous distribution with non-trivial mass at zero. The exact formulas for $\alpha$, $\beta$, and $\lambda$ are: $$\alpha = \frac{2-p}{p-1},\quad \beta = \varphi(p-1)\mu^{p-1},\quad \lambda = \frac{\mu^{2-p}}{\varphi(2-p)}$$ and the probability of zero is $\exp(-\lambda)$ (Halder et al., 2019).

3. GLM Formulation, Link Functions, and IRLS Estimation

In GLMs, Tweedie models are specified by $g(\mu)=\eta=X\beta$. Canonical link $g_c(\mu)=\frac{\mu^{1-p}-1}{1-p}$ arises from the exponential family structure; however, the log-link $g(\mu)=\log\mu$ is preferred for practical reasons and is valid across $p$ (Yilmaz et al., 2012, Manna et al., 9 Jul 2025). Parameter estimation proceeds via IRLS or Newton scoring:

$$\beta^{(new)} = \beta^{(old)} + [X^T W X]^{-1} X^T W z$$

where $W$ involves derivatives and variance according to the link function. For log-link, $W_{ii} = 1/(\varphi \mu_i^p)$; for canonical power link, $W_{ii} = \mu_i^p/(1-p)^2$ (Yilmaz et al., 2012). Maximum likelihood estimation for Tweedie GLMs is rigorously equivalent to $\beta$-divergence minimization:

$$d_\beta(y,\mu) = \frac{1}{\beta(\beta-1)} \left[ y^\beta + (\beta-1)\mu^\beta - \beta y \mu^{\beta-1}\right],\quad \beta=2-p$$

so that maximizing likelihood is minimizing $\sum_i d_\beta(y_i,\mu_i)$ (Yilmaz et al., 2012).

4. Extensions: Double GLMs, Smoothing, and Regularization

Double GLMs (DGLMs) extend Tweedie models to jointly model mean and dispersion via separate predictors: $g_1(\mu_i)=x_i^\top\beta$, $g_2(\varphi_i)=z_i^\top\gamma$ (Halder et al., 2019, Halder et al., 2023, Gu, 2024). Smoothing parameter estimation in models with spline or spatial smooths proceeds by penalized likelihood: $\ell(\beta) - \frac{1}{2}\beta^T S_\lambda \beta$ with $S_\lambda = \sum_j \lambda_j S_j$ quadratic penalties. The generalized Fellner–Schall algorithm updates each smoothing parameter $\lambda_j$ by

$$\lambda_j^{new} = \lambda_j \frac{\mathrm{tr}(S_\lambda^{-} S_j) - \mathrm{tr}(V_\lambda S_j)}{\widehat{\beta}^T S_j \widehat{\beta}}$$

yielding efficient convergence for complex models, including tensor-product and adaptive smoothers (Wood et al., 2016).

In the spatial context, location effects $\theta$ are regularized via graph Laplacian penalties (e.g., $S(\theta) = \theta^T L \theta$) in DGLMs (Halder et al., 2020, Halder et al., 2019), often estimated by block coordinate descent. Bayesian DGLM with spike-and-slab priors incorporates spatial processes and variable selection for high-dimensional covariate spaces (Halder et al., 2023).

5. Specializations: Poisson-Tweedie, PET Models, and Dispersion

The Poisson–Tweedie GLM for counts is defined via hierarchical mixture: $Z_i\sim\mathrm{Tw}_p(\mu_i, \varphi)$, $Y_i|Z_i\sim\mathrm{Poisson}(Z_i)$, leading to variance $V(Y_i) = \mu_i + \varphi\mu_i^p$ (Bonat et al., 2016). This structure allows negative $\varphi$ (underdispersion), equidispersion ($\varphi=0$), or overdispersion ($\varphi>0$), with $p$ adapting between zero-inflation, NB-type tails, or extreme heavy tails. Estimation employs quasi-score for $\beta$ and Pearson estimating functions for $(\varphi,p)$, solved via Newton scoring.

Poisson-exponential-Tweedie (PET) models generalize the variance to $V(Y) = m + m^2 + \varphi m^p$, adding quadratic dispersion beyond Tweedie (Abid et al., 2019). Estimation is performed by quasi-moment methods; diagnostic indices for dispersion and zero-inflation are reported, including those based on "unit negative binomial" (zero-shifted geometric) reference (Abid et al., 2019).

6. Zero-Inflated Tweedie and Decision-Tree Boosting Approaches

The zero-inflated Tweedie model combines a Tweedie variate with an extra "perfect zero" state: $$Y = \begin{cases} 0 & \text{w.p. }\pi(x) \ T & \text{w.p. }1-\pi(x) \end{cases}$$ with all model components ($\pi(x)$, $\mu(x)$, $\varphi(x)$) linked via monotone link functions and potentially modeled nonparametrically via gradient-boosted decision trees (Gu, 2024). Parameter estimation uses a generalized expectation-maximization algorithm, iteratively optimizing cross-entropy (for zero probability), Tweedie deviance (mean), and gamma-type loss (dispersion) by boosting. Empirically, explicit modeling of $\varphi(x)$ and $\pi(x)$ yields accurate recovery of mean and dispersion structure under high zero-inflation.

7. Predictive Uncertainty: Conformal Inference and Model Diagnostics

Distribution-free prediction intervals for Tweedie GLMs and LightGBM models with Tweedie loss can be constructed using split-conformal inference (Manna et al., 9 Jul 2025). Non-conformity measures include Pearson, Anscombe, and deviance residuals; locally weighted Pearson residuals (fit via LightGBM) provide tight, valid interval coverage, especially when residual variance is heteroscedastic. Standard diagnostic residuals (Pearson, deviance, randomized quantile) are used to assess model adequacy (Halder et al., 2019).

Residual Type	Coverage (Mean ± SD)	Interval Width
Pearson (GLMNET)	0.9495 ± 0.0061	14.76 ± 0.62
Locally weighted (LGBM)	0.9502 ± 0.0054	13.96 ± 0.51

Locally weighted approaches outperform unweighted ones in terms of interval width, maintaining nominal coverage (Manna et al., 9 Jul 2025).

8. Applications and Empirical Performance

Tweedie GLMs and their extensions are widely deployed for insurance claims modeling, RNA-seq count data, biological assays, and other semicontinuous responses (Halder et al., 2020, Signorelli et al., 2020, Gu, 2024, Halder et al., 2019). Empirical studies confirm unbiased parameter estimation, efficient asymptotic inference via IRLS/Newton or boosting, automatic adaptation to zero-inflation, heavy tails, or overdispersion by estimation of $p$ and $\varphi$, and stability against moderate index parameter misspecification (Halder et al., 2019). In longitudinal count modeling, Poisson–Tweedie mixed-effects models allow for robust inference under high overdispersion and correlation, with superior Type-I error control and RMSE relative to classical negative-binomial approaches (Signorelli et al., 2020).

9. Model Selection, Orthogonality, and Practical Guidance

Tweedie models allow profile-likelihood selection for the index $p$ or cross-validation under $\beta$-divergence loss (Yilmaz et al., 2012). Mean and dispersion parameters exhibit orthogonality in EDMs with power variance; inference for $\beta$ is robust to $\varphi$ and $p$ misspecification. Algorithms are available in R packages (glmnet with Tweedie loss, ptmixed, mcglm, and LightGBM) (Signorelli et al., 2020, Bonat et al., 2016, Gu, 2024, Manna et al., 9 Jul 2025). For spatial data, Laplacian or Gaussian-process penalties address location-referenced effects (Halder et al., 2020, Halder et al., 2023).

10. Summary Table: Tweedie GLM Family Parameters and Model Types

Power $p$	Model Class	Mean–Variance Relation
$p=0$	Gaussian	$\varphi$
$p=1$	Poisson	$\varphi\mu$
$1	Compound Poisson-Gamma	$\varphi\mu^p$
$p=2$	Gamma	$\varphi\mu^2$
$p=3$	Inverse Gaussian	$\varphi\mu^3$
$p$ arbitrary	Poisson–Tweedie (+linear)	$\mu + \varphi\mu^p$
$p$ arbitrary	PET (quad. term added)	$m + m^2 + \varphi m^p$

Comprehensive frameworks for estimation, diagnostics, regularization, and uncertainty quantification position Tweedie GLMs as foundational tools for analyzing complex, semicontinuous datasets in diverse research areas.