normal_family <- family(
Normal(),
links = c(mu = "identity", sigma = "log")
)
normal_familyFamily: Normal(identity(mu), log(sigma))
---
Derivative functions:
..$ score
.. ..$ mu
.. ..$ sigmaFamily Objects
Family objects connect a response distribution with the predictors and optimization machinery used by gamlss2. There are two complementary ways to define a family:
- bridge an existing
distributions3distribution, or - construct a
gamlss2.familylist directly for a distribution that needs custom behavior.
The first route should usually be preferred. It gives the fitted model the same distribution interface used by distributions3 for densities, probabilities, quantiles, moments, and simulation.
For the fitting workflow, start with First Steps. To modify an existing distribution for an observation mechanism such as truncation or censoring, see Family Extensions.
1 The family interface
A gamlss2.family object contains the distribution-specific information needed by the fitting algorithm. The required components are:
| Component | Purpose |
|---|---|
family |
Distribution name |
names |
Names of the distribution parameters |
links |
One link function for each parameter |
pdf(par, y, log = FALSE, ...) |
Log-density or density evaluation |
Useful additional components are cdf(), quantile(), random(), support(), mean(), variance(), score(), hessian(), update(), initialize(), and valid.response(). Missing probability, quantile, moment, and random-number methods can often be constructed numerically from the density and support, but analytical methods are preferable when they are available.
The two distribution interfaces use slightly different conventions. A family callback receives a named parameter list, for example par$mu and par$sigma, while a distributions3 method receives an object inheriting from distribution. The bridge described below translates between these forms.
Callback arguments have consistent meanings throughout a family: par is a named list of distribution parameters on the response scale, y contains the observed response values, and eta contains linear-predictor values used by an update method. The p argument denotes probabilities for quantiles, n is the requested number of random draws, and log = TRUE asks pdf() for log-densities. Additional arguments are passed through ....
2 Bridging a distributions3 distribution
The family() method turns a distributions3 object into a gamlss2.family. The names in links select the parameters to model and the values specify their links.
The resulting object contains the family metadata as well as wrappers for pdf(), cdf(), quantile(), random(), the moments, and the log-likelihood. The distribution constructor is retained in create_distribution, which is used later when a fitted model is converted back to a distributions3 object.
names(normal_family) [1] "family" "names" "links"
[4] "log_likelihood" "mu" "pdf"
[7] "cdf" "random" "quantile"
[10] "crps" "mean" "variance"
[13] "skewness" "kurtosis" "create_distribution"
[16] "support" "score" "valid.response" normal_family$names[1] "mu" "sigma"vapply(normal_family$links, function(x) x$name, character(1)) mu sigma
"identity" "log" The bridge result can now be used like any other family. This example uses the cars data from base R and allows both the normal mean and standard deviation to vary with speed.
m <- gamlss2(
dist ~ s(speed) | s(speed),
data = cars,
family = normal_family,
trace = FALSE
)
summary(m)Call:
gamlss2(formula = dist ~ s(speed) | s(speed), data = cars, family = normal_family,
trace = FALSE)
---
Family: Normal(identity(mu), log(sigma))
*--------
Coefficients:
Estimate Std. Error t value Pr(>|t|)
mu.(Intercept) 42.425 1.782 23.81 <2e-16 ***
sigma.(Intercept) 2.638 0.100 26.38 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
edf
mu.s(speed) 1.2175
sigma.s(speed) 1.0000
*--------
n = 50 df = 4.22 res.df = 45.78
Deviance = 405.6622 Null Dev. Red. = 12.92%
AIC = 414.0972 elapsed = 0.04secBy default, predict() performs the reverse operation: it extracts the fitted parameters and returns a vector of distributions3 objects. The explicit prodist() interface is equivalent.
nd <- data.frame(speed = c(10, 20, 30))
d <- predict(m, newdata = nd)
d 1 2
"Normal(mu = 23.07, sigma = 10.03)" "Normal(mu = 58.88, sigma = 18.55)"
3
"Normal(mu = 96.16, sigma = 34.31)" All ordinary distributions3 methods can then be used. The elementwise argument controls whether a vector of probabilities is paired with the distribution vector or expanded to all distribution-probability combinations.
mean(d) 1 2 3
23.06591 58.88142 96.16345 variance(d) 1 2 3
100.6198 344.1271 1176.9385 cdf(d, 50) 1 2 3
0.99637459 0.31605261 0.08921306 quantile(d, c(0.1, 0.5, 0.9), elementwise = FALSE) q_0.1 q_0.5 q_0.9
1 10.21074 23.06591 35.92109
2 35.10779 58.88142 82.65505
3 52.19786 96.16345 140.12905quantile(d, c(0.1, 0.5, 0.9), elementwise = TRUE) 1 2 3
10.21074 58.88142 140.12905 random(d, 2) r_1 r_2
1 34.22014 11.21147
2 69.41028 28.83345
3 76.39498 125.71896This gives the complete workflow:
distributions3 object -> family() -> gamlss2.family
-> gamlss2() -> fitted model -> predict() -> distributions3 object3 Constructor arguments and distributions3_family()
The lower-level distributions3_family() function is useful when a distribution constructor has fixed arguments that are not model parameters. For example, a binomial distribution can have a fixed number of trials while its probability is modeled.
binomial_family <- distributions3_family(
"Binomial",
links = c(p = "logit"),
create_distribution = "do.call",
size = 10
)
binomial_distribution <- binomial_family$create_distribution(
list(p = rep(0.5, 3))
)
binomial_distribution[1] "Binomial(size = 10, p = 0.5)" "Binomial(size = 10, p = 0.5)"
[3] "Binomial(size = 10, p = 0.5)"distributions3::pdf(binomial_distribution, 1:3)[1] 0.009765625 0.043945312 0.117187500With create_distribution = "do.call", the fixed size argument is passed to the constructor every time a distribution object is created. The default "structure" strategy is faster, but is appropriate only when adding the class to the parameter object is sufficient.
The workhorse also accepts a custom create_distribution function when a distribution needs a specialized construction step. In all cases, the returned object must inherit from gamlss2.family and its parameter names must match the supplied links.
4 Scores, Hessians, and updates
For a gamlss2.family written by hand, score() and hessian() are functions on the linear-predictor scale. A distributions3 score or Hessian is generally provided with respect to a distribution parameter. The bridge applies the link-function chain rule when score = TRUE or hessian = TRUE.
normal_family <- family(
Normal(),
links = c(mu = "identity", sigma = "log"),
score = TRUE,
hessian = TRUE,
update = TRUE
)
names(normal_family[c("score", "hessian", "update")])[1] "score" "hessian" "update" Analytical derivatives can improve speed and numerical stability. If they are not available, gamlss2 approximates the required derivatives numerically. Cross-derivatives can be supplied by the underlying distribution when a second-order optimizer such as CG can benefit from them.
Flexible links are represented by link-gamlss2 objects. When analytical score and Hessian callbacks are requested, the bridge uses derivatives from the link object. For example, the softplus link can be passed to the bridge after construction with make.link2(). In this example the derivative flags are disabled explicitly, so numerical derivatives are used.
normal_family <- family(
Normal(),
links = list(mu = make.link2(softplus), sigma = "log"),
score = FALSE, hessian = FALSE
)
## The softplus link does not provide all second-link derivatives needed by
## the analytic update path, so the bridge uses numerical derivatives here.
m_softplus <- gamlss2(
dist ~ s(speed) | s(speed),
data = cars,
family = normal_family,
trace = FALSE
)
AIC(m, m_softplus) AIC df
m_softplus 414.0374 4.228447
m 414.0972 4.2175055 Numerical completion and support
Families do not need to implement every distributional method. Given a density and a valid support, gamlss2 can construct a missing CDF by integration or summation, invert a CDF to obtain quantiles, calculate moments, and generate random values by inverse transform sampling.
These fallbacks are convenient for prototypes, but they may be expensive and need a trustworthy support. In particular, provide support() explicitly for shifted or truncated distributions, unusual continuous supports, and discrete families that are not counts on the nonnegative integers. Analytical CDFs and quantiles should be retained whenever possible.
The type component should be either "continuous" or "discrete", and valid.response() should reject invalid response values. The type and support are used by residual diagnostics, numerical completion, and the distributions3 methods.
6 A reparameterized distribution from scratch
Suppose we want a Gamma distribution with the GAMLSS parameterization (mu, sigma), where mu is the mean and sigma controls dispersion. The ordinary distributions3::Gamma() uses shape and rate instead. We can define a small distributions3 class that translates between the two parameterizations:
\[ \text{shape}=1/\sigma^2, \qquad \text{rate}=1/(\mu\sigma^2). \]
The following example implements the distribution methods needed by the bridge, including an expected Hessian. The Hessian is supplied on the (mu, sigma) scale; family() then applies the link-function transformation.
Gamma2 <- function(mu = 1, sigma = 1) {
stopifnot(all(mu > 0), all(sigma > 0))
structure(list(mu = mu, sigma = sigma),
class = c("Gamma2", "distribution"))
}
gamma2_base <- function(d) {
Gamma(shape = 1 / d$sigma^2, rate = 1 / (d$mu * d$sigma^2))
}
pdf.Gamma2 <- function(d, x, log = FALSE, elementwise = FALSE, ...) {
dgamma(x, shape = 1 / d$sigma^2,
rate = 1 / (d$mu * d$sigma^2), log = log)
}
log_pdf.Gamma2 <- function(d, x, elementwise = FALSE, ...) {
dgamma(x, shape = 1 / d$sigma^2,
rate = 1 / (d$mu * d$sigma^2), log = TRUE)
}
cdf.Gamma2 <- function(d, x, ...) {
cdf(gamma2_base(d), x, ...)
}
quantile.Gamma2 <- function(x, probs, ...) {
quantile(gamma2_base(x), probs, ...)
}
random.Gamma2 <- function(x, n, ...) {
random(gamma2_base(x), n, ...)
}
mean.Gamma2 <- function(x, ...) x$mu
variance.Gamma2 <- function(x, ...) x$mu^2 * x$sigma^2
support.Gamma2 <- function(d, drop = TRUE, ...) {
min <- rep(0, length(d))
max <- rep(Inf, length(d))
make_support(min, max, d, drop = drop)
}
score.Gamma2 <- function(d, y, which, ...) {
mu <- d$mu
sigma <- d$sigma
if(which == "mu")
(y - mu) / (mu^2 * sigma^2)
else
2 / sigma^3 * (y / mu - 1 - log(y / (mu * sigma^2)) +
digamma(1 / sigma^2))
}
hessian.Gamma2 <- function(d, y, which, expected = TRUE, ...) {
mu <- d$mu
sigma <- d$sigma
if(which == "mu")
-1 / (mu^2 * sigma^2)
else if(which == "sigma")
4 / sigma^4 - 4 / sigma^6 * trigamma(1 / sigma^2)
else
rep(0, length(y))
}
fGA2 <- family(
Gamma2(),
links = c(mu = "log", sigma = "log"),
score = TRUE, hessian = TRUE, expected = TRUE,
initialize = list(
mu = function(y, ...) (y + mean(y)) / 2,
sigma = function(y, ...) rep(1, length(y))
)
)
print(fGA2)Family: Gamma2(log(mu), log(sigma))
---
Derivative functions:
..$ score
.. ..$ mu
.. ..$ sigma
..$ hessian
.. ..$ mu
.. ..$ sigmaair <- subset(airquality, !is.na(Ozone))
m <- gamlss2(
Ozone ~ s(Temp) + s(Wind) | s(Temp),
data = air,
family = fGA2,
trace = FALSE
)
summary(m)Call:
gamlss2(formula = Ozone ~ s(Temp) + s(Wind) | s(Temp), data = air,
family = fGA2, trace = FALSE)
---
Family: Gamma2(log(mu), log(sigma))
*--------
Coefficients:
Estimate Std. Error t value Pr(>|t|)
mu.(Intercept) 3.55382 0.02682 132.49 <2e-16 ***
sigma.(Intercept) -0.87025 0.06331 -13.75 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
edf
mu.s(Temp) 4.4762
mu.s(Wind) 1.0718
sigma.s(Temp) 3.9513
*--------
n = 116 df = 11.5 res.df = 104.5
Deviance = 929.2643 Null Dev. Red. = 14.2%
AIC = 952.263 elapsed = 0.18secThis pattern is useful when an existing distribution has the right likelihood but the wrong parameterization. Prediction and Uncertainty uses this same Gamma family to construct intervals for means, quantiles, and future observations. A user-defined family can also be assembled directly as a gamlss2.family, as shown next.
7 Writing a family directly
When no suitable distributions3 distribution exists, a family can be written as a classed list. The following is a compact normal family with analytical probability and moment methods. The random() signature is random(par, n), matching the gamlss2.family interface.
my_normal_family <- function(mu.link = "identity", sigma.link = "log") {
fam <- list(
family = "my_normal_family",
names = c("mu", "sigma"),
links = list(mu = mu.link, sigma = sigma.link),
pdf = function(par, y, log = FALSE, ...) {
dnorm(y, mean = par$mu, sd = par$sigma, log = log)
},
cdf = function(par, y, ...) {
pnorm(y, mean = par$mu, sd = par$sigma, ...)
},
quantile = function(par, p, ...) {
qnorm(p, mean = par$mu, sd = par$sigma)
},
random = function(par, n, ...) {
rnorm(n, mean = par$mu, sd = par$sigma)
},
mean = function(par, ...) par$mu,
variance = function(par, ...) par$sigma^2,
support = c(-Inf, Inf),
type = "continuous",
valid.response = function(y) is.numeric(y)
)
class(fam) <- "gamlss2.family"
fam
}
my_family <- my_normal_family()
m <- gamlss2(
dist ~ s(speed) | s(speed),
data = cars,
family = my_family,
trace = FALSE
)
summary(m)Call:
gamlss2(formula = dist ~ s(speed) | s(speed), data = cars, family = my_family,
trace = FALSE)
---
Family: my_normal_family(identity(mu), log(sigma))
*--------
Coefficients:
Estimate Std. Error t value Pr(>|t|)
mu.(Intercept) 42.425 1.782 23.81 <2e-16 ***
sigma.(Intercept) 2.638 0.100 26.38 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
edf
mu.s(speed) 1.2175
sigma.s(speed) 1.0000
*--------
n = 50 df = 4.22 res.df = 45.78
Deviance = 405.6622 Null Dev. Red. = 12.92%
AIC = 414.0972 elapsed = 0.02secFor a manually defined family, the minimum required interface is family, names, links, and pdf(). Derivative callbacks have the following shape:
score = list(
mu = function(par, y, ...) {
(y - par$mu) / par$sigma^2
},
sigma = function(par, y, ...) {
-1 + (y - par$mu)^2 / par$sigma^2
}
)
hessian = list(
mu = function(par, y, ...) 1 / par$sigma^2,
sigma = function(par, y, ...) rep(2, length(y)),
"mu:sigma" = function(par, y, ...) rep(0, length(y))
)The Hessian entries are the negative of the expected second derivatives on the predictor scale. An optional update(par, y, eta, which) function can provide a specialized Newton/Fisher-scoring update returning list(eta, weights).
8 Optional scoring rules
The distributions3 bridge also exposes a crps() callback when the optional scoringRules package is available. This is useful for probabilistic forecast evaluation, but it is not needed for fitting, prediction, or the basic family interface. Forecasting and Assessment demonstrates scoring rules and graphical diagnostics.
9 Summary
Use the distributions3 bridge when the response distribution already has a distributions3 implementation. It provides a consistent probability and moment interface and handles link-scale derivative transformations. Use a manual gamlss2.family when a distribution needs custom density, derivative, support, initialization, or optimization logic. In either case, a good family should provide reliable support information and analytical distributional methods whenever practical.