families.Rmd
Family objects are important building blocks in the design of BAMLSS models. They specify the distribution by collecting functions of the density, respective log-likelihood, first-order derivatives of the log-likelihood w.r.t. predictors (the score function), and (optionally) second-order derivatives of the log-likelihood w.r.t. predictors or their expectation (the Hessian).
The bamlss package can be easily extended by constructing families for specific tasks, i.e. problems for which a likelihood can be formulated. However, the most important distributions are shipped with bamlss. Finally, it is also possible to employ gamlss families.
We illustrate how to build a bamlss family by hand along the Gaussian distribution, with density \[ f(y\,|\,\mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma} \cdot \exp \left( \frac{-(y-\mu)^2}{2\sigma^2} \right), \] and log-likelihood function \[ \ell(\mu,\sigma\,|\,y) = - \frac{1}{2} \log(2\pi) - \log(\sigma) - \frac{(y-\mu)^2}{2\sigma^2}, \] for an individual observation. The sum of the log-likelihood function over all observations is the target function of the optimization problem.
In the distributional regression framework the parameters are linked to predictors by link functions, \[ \mu = \eta_\mu, \qquad \log(\sigma) = \eta_\sigma. \] For the Gaussian \(\mu\) and \(\sigma\) are linked to \(\eta_\mu\) and \(\eta_\sigma\) by the identity function and the logarithm, respectively.
The score functions in bamlss are the first derivatives of the log-likelihood w.r.t. the predictors: \[ s_\mu = \frac{\partial\ell}{\partial\eta_\mu} = \frac{\partial\ell}{\partial\mu} \cdot \frac{\partial\mu}{\partial\eta_\mu} = \frac{y-\mu}{\sigma^2}, \] and \[ s_\sigma = \frac{\partial\ell}{\partial\eta_\sigma} = \frac{\partial\ell}{\partial\sigma} \cdot \frac{\partial\sigma}{\partial\eta_\sigma} = -1 + \frac{(y-\mu)^2}{\sigma^2}. \]
For the second derivative of the log-likelihood we are able to obtain the negative expectation, \[ \mathsf{E}(-\partial^2\ell / \partial\eta_{\mu}^{2} ) = \sigma^{-2}, \] and \[ \mathsf{E}(-\partial^2\ell / \partial\eta_{\sigma}^{2} ) = 2. \]
Now we have to write a function that returns a
family.bamlss
object (S3) which encapsulates functions for
density, score and Hessian, and the names of the family, parameter and
link functions:
Name of element | Value |
---|---|
family |
character string with the name of the family |
names |
vector of character strings with the names of the parameters |
links |
vector of character strings with the names of the link functions |
d |
a function returning the density with arguments
y, par, log = FALSE (see below) |
score |
a list with functions (one for each parameter) returning the first derivatives of the log-likelihood w.r.t. predictors |
hess |
a list with functions (one for each parameter) returning the negative second derivatives of the log-likelihood w.r.t. predictors |
Nearly all functions take as first argument the response
y
and as second argument a named list holding the evaluated
parameters par
of the distribution on the scale of the
predictors \(\eta_{\star}\).
The implementation looks like this:
mygauss <- function(...) {
f <- list(
"family" = "mygauss",
"names" = c("mu", "sigma"),
"links" = c(mu = "identity", sigma = "log"),
"d" = function(y, par, log = FALSE) {
dnorm(y, mean = par$mu, sd = par$sigma, log = log)
},
"score" = list(
mu = function(y, par, ...) {
drop((y - par$mu) / (par$sigma^2))
},
sigma = function(y, par, ...) {
drop(-1 + (y - par$mu)^2 / (par$sigma^2))
}
),
"hess" = list(
mu = function(y, par, ...) {
drop(1 / (par$sigma^2))
},
sigma = function(y, par, ...) {
rep(2, length(y))
}
)
)
class(f) <- "family.bamlss"
return(f)
}
library("bamlss")
## Simulate some data.
d <- GAMart()
## Specify formula.
f <- ~ s(x1) + s(x2) + s(x3) + s(lon, lat)
f <- list(mu = update(f, num ~ .), sigma = f)
## Model.
b <- bamlss(f, data = d, family = mygauss)
Optionally, the family.bamlss
object can be extended by
functions for
p(y, par, ...)
,q(p, par)
,r(n, par)
,loglik(y, par)
,mu(par, ...)
,which can help to speed up optimization, or be convenient for predictions and simulations.
There are some families implemented:
Or normal distribution.
Density: \[ f(y\,|\,\mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma} \cdot \exp \left( \frac{-(y-\mu)^2}{2\sigma^2} \right), \] with \(\sigma > 0\).
Links:
Function: gaussian_bamlss(...)
Alternative parameterizations:
gaussian2_bamlss(...)
models inverse sigma
with a log link.Gaussian_bamlss(...)
one parametric family modelling
only mu
.Or skewed-logistic distribution.
Density: \[ f(y\,|\,\mu,\sigma,\alpha)=\frac{\alpha\cdot\exp\left(-\frac{y-\mu}{\sigma}\right)} {\sigma\cdot\left(1+\exp\left(-\frac{y-\mu}{\sigma}\right)\right)^{2\alpha}}, \] with \(\sigma, \alpha > 0\).
Links:
Function: glogis_bamlss(...)
Two-parameter distribution with \(\mu > 0\) for the expectation of the distribution and \(\sigma > 0\) for its shape.
Density: \[ f(y\,|\,\mu,\sigma)=\frac{y^{\sigma-1}\cdot\exp\left(-\frac{\sigma\cdot y}{\mu}\right)} {\left(\frac{\mu}{\sigma}\right)^\sigma\cdot\Gamma(\sigma)}. \]
Links:
Function: gamma_bamlss(...)
Density: \[ f(y\,|\,\xi,\sigma)=\frac{1}{\sigma}\cdot\left(1+\frac{\xi\cdot y}{\sigma} \right)^{-\left(1+1/\xi\right)}, \] with \(\xi, \sigma > 0\).
Links:
Function: gpareto_bamlss(...)
Density: \[ f(y\,|\,\lambda,\alpha)=\frac{\alpha}{\lambda}\cdot\left(\frac{y}{\lambda}\right)^{\alpha-1} \cdot\exp\left(-\left(\frac{y}{\lambda}\right)^\alpha\right), \] with \(\lambda, \alpha > 0\).
Links:
Function: weibull_bamlss(...)
Density: \[ f(y\,|\,\mu,\sigma)= \begin{cases} \frac{1}{\sigma} \phi\left(\frac{y-\mu}{\sigma} \right) & y > 0 \\ \Phi\left(\frac{-\mu}{\sigma} \right) & \mbox{else,} \end{cases} \] where \(\phi\) is the probability density and \(\Phi\) the cumulative distribution function of the standard normal distribution.
Links:
Function: cnorm_bamlss(...)
The beta distribution, with density \[ f_{Beta}(y\,|\,\alpha,\beta)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\cdot\Gamma(\beta)} \cdot y^{\alpha-1} \cdot (1-y)^{\beta-1} \] is implemented using the Balding-Nichols parameterization (Balding and Nichols 1995).
Density: \[ f(y\,|\,\mu,\sigma)=f_{Beta}\left(y\,\middle|\,\alpha = \mu \cdot \frac{1-\sigma}{\sigma}, \beta = (1-\mu) \cdot \frac{1-\sigma}{\sigma} \right), \] with \(0 < \mu < 1\) and \(0 < \sigma < 1\).
Links:
Function: beta_bamlss(...)
Often called binomial, though the Bernoulli is a special case of the binomial distribution.
Density: \[ f(y\,|\,\pi) = \pi^y \cdot (1-\pi)^{1-y}, \] where \(\pi\) is a probability.
Links:
Function:
binomial_bamlss(link = "logit", ...)
Say, \(y\) is a categorial response with \(c+1\) nomial categories, and \(y_{1}^{\star}, \dots, y_{c+1}^{\star}\) are associated dummy variables. \(p_1, \dots, p_{c+1}\) are probabilities for an observation falling in one of the categories. Category \(c+1\) serves as reference, i.e. \(p_{c+1} = 1 - p_1 - \dots - p_{c}\).
Density: \[ f(y\,|\,\pi_1,\dots,\pi_c) = \prod_{s=1}^{c+1} p_{s}^{y_{s}^{\star}}, \quad \text{where} \quad p_r = \frac{\pi_r}{1 + \pi_1 + \dots + \pi_c}, \]
Links:
where \(\eta_r\) specifies the log-odds between category \(r\) and the reference category \(c+1\).
Function: multinomial_bamlss(...)
Density: \[ f(y\,|\,\lambda) = \frac{\lambda^y \cdot \exp(-\lambda)}{y!}, \] with \(\lambda > 0\).
Link:
Function: poisson_bamlss(...)
Density: \[ f(y\, |\, \mu, \theta) = \frac{\Gamma(\theta + y)}{\Gamma({\theta}) \cdot y!} \cdot \frac{\mu^y \cdot \theta^\theta}{(\mu + \theta)^{\theta + y}}, \] with \(\mu,\theta > 0\).
Links:
Function: nbinom_bamlss(...)
Density: \[ f(y\,|\,\mu,\theta)=\frac{f_{\mathrm{NB}}(y\,|\,\mu,\theta)} {1-f_{\mathrm{NB}}(0\,|\,\mu,\theta)}, \] where \(f_{\mathrm{NB}}\) is the density of the negative binomial.
Links:
Function: ztnbinom_bamlss(...)
Say \(y\) is a \(k\) dimensional observation.
Density: \[ f(y\,|\,\mu,\Sigma) = \frac{1}{\sqrt{(2\pi)^{k} |\Sigma|}} \cdot \exp \left(-\frac{1}{2} (y-\mu)^\top \Sigma^{-1} (y-\mu) \right), \] where \(\Sigma\) can be decomposed into \(\Sigma = D \Omega D\) with \[ D= \begin{pmatrix} \sigma_1 & 0 & \cdots & 0 \\ 0 & \sigma_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma_k \end{pmatrix}, \qquad \text{and} \qquad \Omega= \begin{pmatrix} 1 & \rho_{12} & \cdots & \rho_{1k} \\ \rho_{12} & 1 & \cdots & \rho_{2k} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{1k} & \rho_{2k} & \cdots & 1 \end{pmatrix}. \]
Links:
Function: mvnorm_bamlss(k = 2, ...)
Alternative parameterizations:
mvnormAR1_bamlss(k = 2, ...)
with \[
\Omega=
\begin{pmatrix}
1 & \rho & \rho^2 & \cdots & \rho^{k-1} \\
\rho & 1 & \rho & \cdots & \rho^{k-2} \\
\rho^2 & \rho & 1 & \cdots & \rho^{k-3} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\rho^{k-1} & \rho^{k-2} & \rho^{k-3} & \cdots & 1
\end{pmatrix},
\] and \(\mathrm{rhogit}(\rho) =
\eta_{\rho}\).The bamlss package supports gamlss families via a transfer
function that builds the bamlss.family
object from a GAMLSS
distribution funciton (Stasinopoulos and Rigby
2007), e.g. NO()
, internally:
## Load the gamlss families.
library("gamlss.dist")
## Estimate model.
b <- bamlss(f, data = d, family = NO)