Stasinopoulos 2007 - GAMLSS Model Framework

Summary

A GAMLSS assumes independent observations $y_i\sim f(y_i\mid {\mu }_i,{\sigma }_i,{\nu }_i,{\tau }_i)$ and gives each distribution parameter its own monotone link and an additive predictor of linear terms plus smooth or random-effect terms. It strictly generalises GLM and GAM by dropping the exponential-family restriction and by letting scale and shape depend on covariates.

Overview

GLMs model $E(y)$ through a link and treat dispersion as constant or mean-tied. GAMs add smooth terms to the mean. GAMLSS (Rigby and Stasinopoulos 2001, 2005; Akantziliotou et al. 2002) are “semi-parametric”: a parametric distribution assumption for $y$ (“parametric”) with non-parametric smooths in the parameter predictors (“semi”). The framework suits leptokurtic/platykurtic, skewed, or overdispersed responses and heterogeneity where scale or shape moves with covariates.

Main Content

GAMLSS model (Eq. 1)

Independent $y_i$ ($i=1,\dots ,n$) have density $f(y_i\mid {\mathbf{\theta }}_i)$ with ${\mathbf{\theta }}_i=({\theta }_{1i},{\theta }_{2i},{\theta }_{3i},{\theta }_{4i})=({\mu }_i,{\sigma }_i,{\nu }_i,{\tau }_i)$. For $k=1,\dots ,4$, with monotone known links $g_k(\cdot )$,

$$g_k({\mathbf{\theta }}_k)={\mathbf{\eta }}_k={\mathbf{X}}_k{\mathbf{\beta }}_k+\sum _{j=1}^{J_k}{\mathbf{Z}}_{jk}{\mathbf{\gamma }}_{jk}$$

where ${\mathbf{X}}_k$ is a known $n\times J_k^{\prime }$ design matrix, ${\mathbf{\beta }}_k$ the linear coefficients, ${\mathbf{Z}}_{jk}$ known $n\times q_{jk}$ design matrices, and ${\mathbf{\gamma }}_{jk}\sim N_{q_{jk}}(0,{\mathbf{G}}_{jk}^{-1})$ with ${\mathbf{G}}_{jk}={\mathbf{G}}_{jk}({\mathbf{\lambda }}_{jk})$ a (generalised-inverse) symmetric precision matrix depending on hyperparameters ${\mathbf{\lambda }}_{jk}$. If ${\mathbf{G}}_{jk}$ is singular, ${\mathbf{\gamma }}_{jk}$ has improper density $\propto \exp (-\frac{1}{2}{\mathbf{\gamma }}_{jk}^{\top }{\mathbf{G}}_{jk}{\mathbf{\gamma }}_{jk})$.

The first two parameters are conventionally location ($\mu$) and scale ($\sigma$); remaining parameters are shape (skewness $\nu$, kurtosis $\tau$). Seldom do all parameters need covariates.

Sub-models (Eqs. 2-5)

• Semi-parametric additive (${\mathbf{Z}}_{jk}={\mathbf{I}}_n$, ${\mathbf{\gamma }}_{jk}=h_{jk}(x_{jk})$): $g_k({\mathbf{\theta }}_k)={\mathbf{X}}_k{\mathbf{\beta }}_k+{\sum }_j{h}_{jk}(x_{jk})$ (Eq. 2).
• Parametric linear (no additive terms): $g_k({\mathbf{\theta }}_k)={\mathbf{X}}_k{\mathbf{\beta }}_k$ (Eq. 3).
• Non-linear semi-parametric : $g_k({\mathbf{\theta }}_k)=h_k({\mathbf{X}}_k,{\mathbf{\beta }}_k)+{\sum }_j{h}_{jk}(x_{jk})$ (Eq. 4).
• Non-linear parametric: $g_k({\mathbf{\theta }}_k)=h_k({\mathbf{X}}_k,{\mathbf{\beta }}_k)$ (Eq. 5).

Penalized log-likelihood (Eq. 6)

For fixed hyperparameters ${\lambda }_{jk}$, the ${\mathbf{\beta }}_k,{\mathbf{\gamma }}_{jk}$ maximise

$${\ell }_p=\ell -\frac{1}{2}\sum _{k=1}^p\sum _{j=1}^{J_k}{\lambda }_{jk}{\mathbf{\gamma }}_{jk}^{\top }{\mathbf{G}}_{jk}{\mathbf{\gamma }}_{jk},\ell =\sum _{i=1}^n\log f(y_i\mid {\mathbf{\theta }}_i).$$

For parametric models ${\ell }_p=\ell$.

Restrictions in the R implementation. Only $\log f$ and its first (optionally expected second and cross) derivatives wrt each parameter must be computable; numerical derivatives are allowed. Information-orthogonal parameters (cross derivatives zero) make fitting easier; the authors prefer parameterisations in which $\mu$ is a location (mean, median or mode).

Interval-censored data (Eqs. 7-8): the likelihood $L(\theta )=\prod f(y_i\mid \theta )$ becomes $L(\theta )=\prod [F(y_{2i}\mid \theta )-F(y_{1i}\mid \theta )]$ for observed intervals $(y_{1i},y_{2i}]$. Finite mixtures (Eq. 9): $f_Y(y\mid \psi )={\sum }_{k=1}^K{\pi }_k{f}_k(y\mid {\theta }_k)$, fitted by EM; mixing probabilities may themselves depend on covariates.

Examples

Reading the formula for the CD4 model

For CD4 counts vs age the paper fits cd4 ~ cs(age), sigma.fo = ~cs(age) with the normal family. Here $k=1$: $g_1=\mathrm{identity}$, ${\eta }_1={\beta }_{01}+h_{11}(\text{age})$; $k=2$: $g_2=\log$, $\log \sigma ={\beta }_{02}+h_{12}(\text{age})$. This is a smooth mean and smooth log-standard-deviation. See Stasinopoulos 2007 - Lice, CD4 and Claims Examples.

Connections

• Strictly contains GLMs (Generalized Linear Models) and GAMs.
• Smooth terms are penalised (P-splines, cubic smoothing splines) — see Stasinopoulos 2007 - Distributions and Additive Terms. The Bayesian reading: ${\mathbf{\gamma }}_{jk}\sim N(0,{\lambda }^{-1}{\mathbf{K}}^{-})$ is a (partially improper) Gaussian prior and Eq. 6 gives the posterior mode; see Stasinopoulos 2007 - Bayesian P-Splines and Distributional Regression in PyMC.
• Hierarchical analogue: random-effect terms $\mathbf{Z}\mathbf{\gamma }$ are partial pooling (Hierarchical Models).

See Also

Stasinopoulos 2007 - GAMLSS in R - Overview, Stasinopoulos 2007 - Fitting Algorithms and Model Selection