library("gamlss2")
Normal <- function(...) {
fam <- list(
"family" = "Normal",
"names" = c("mu", "sigma"),
"links" = c("mu" = "identity", "sigma" = "log"),
"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))
}
),
"log_likelihood" = function(par, y, ...) {
sum(dnorm(y, par$mu, par$sigma, log = TRUE))
},
"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, ...)
},
"random" = function(par, n) {
rnorm(n, mean = par$mu, sd = par$sigma)
},
"quantile" = function(par, p) {
qnorm(p, mean = par$mu, sd = par$sigma)
},
"support" = c(-Inf, Inf),
"initialize" = list(
"mu" = function(y, ...) { (y + mean(y)) / 2 },
"sigma" = function(y, ...) { rep(sd(y), length(y)) }
),
"mean" = function(par) par$mu,
"variance" = function(par) par$sigma^2,
"valid.response" = function(x) {
if(is.factor(x) | is.character(x))
stop("the response should be numeric!")
return(TRUE)
}
)
class(fam) <- "gamlss2.family"
return(fam)
}
## specify the model formula
f <- Ozone ~ s(Wind) | s(Wind)
## estimate model
m1 <- gamlss2(f, data = airquality, family = Normal)
## plot estimated effects
plot(m1, which = "effects")
## plot diagnostics
plot(m1, which = "resid")
## predict distributions
nd <- model.frame(m1)
nd <- nd[order(nd$Wind), ]
d <- predict(m1, newdata = nd)
## predict quantiles
pq <- quantile(d, probs = c(0.05, 0.5, 0.95), elementwise = FALSE)
## visualize
plot(Ozone ~ Wind, data = nd,
pch = 19, col = rgb(0.1, 0.1, 0.1, alpha = 0.3),
ylim = range(nd$Ozone, pq))
matlines(nd$Wind, pq, lwd = 2,
lty = 1, col = 4)
## another example that defines only the density
## function; in this case all derivatives are
## approximated numerically. For residual diagnostics,
## the $cdf() and $quantile() functions are needed as well
Gamma <- function(...) {
fam <- list(
"names" = c("mu", "sigma"),
"links" = c("mu" = "log", "sigma" = "log"),
"pdf" = function(par, y, log = FALSE, ...) {
shape <- par$sigma
scale <- par$mu/par$sigma
dgamma(y, shape = shape, scale = scale, log = log)
},
"cdf" = function(par, y, lower.tail = TRUE, log.p = FALSE) {
shape <- par$sigma
scale <- par$mu/par$sigma
pgamma(y, shape = shape, scale = scale,
lower.tail = lower.tail, log.p = log.p)
},
"quantile" = function(par, p, lower.tail = TRUE, log.p = FALSE) {
shape <- par$sigma
scale <- par$mu/par$sigma
qgamma(p, shape = shape, scale = scale,
lower.tail = lower.tail, log.p = log.p)
}
)
class(fam) <- "gamlss2.family"
return(fam)
}
## model formula
f <- Ozone ~ s(Temp) + s(Wind) | .
## estimate model
m2 <- gamlss2(f, data = airquality, family = Gamma)
## visualize estimated effects
plot(m2, which = "effects")
## diagnostics, needs the $cdf() and $quantile() function!
plot(m2, which = "resid")Family Objects in gamlss2
Description
A “gamlss2.family” object describes a response distribution for use with gamlss2. It tells the fitting algorithm which distributional parameters exist, which link functions they use, and which distributional functions are available, such as densities, quantiles, random-number generators, or score functions.
Usage
## S3 method for class 'gamlss2'
family(object, ...)
## S3 method for class 'gamlss2.family'
print(x, full = TRUE, ...)
Arguments
object
|
A fitted “gamlss2” model object.
|
x
|
A “gamlss2.family” object.
|
full
|
Logical; print available family-specific and derivative functions? |
…
|
Not used currently. |
Details
A gamlss2 family object is a list of class “gamlss2.family”. It defines the response distribution and the functions needed to fit and use that distribution in a model.
The minimum required elements are:
-
family: the name of the distribution (character string). -
names: a character vector of parameter names (e.g.,c(“mu”, “sigma”)). -
links: a named character vector specifying the link function for each parameter (e.g.,c(mu = “identity”, sigma = “log”)), or a list of link functions, seesoftplus. -
pdf(par, y, log = FALSE, …): a function to evaluate the (log-)density.
The pdf() function must accept the response y, a named list par of evaluated parameter values (e.g., par$mu, par$sigma), a logical log, and optional additional arguments.
Further elements can be added depending on how fully the family should support fitting, prediction, diagnostics, and simulation. Common optional elements are:
-
score: a named list of functions (one per parameter), each computing the first derivative of the log-likelihood with respect to the linear predictor:score[param]. -
hessian: a named list of functions computing second derivatives (the negative Hessian). For parametersmuandsigma, this includes:hessian[“mu”],hessian[“sigma”], and optionally cross derivatives likehessian[“mu:sigma”]. -
log_likelihood(par, y, …): a function computing the total log-likelihood. -
cdf(par, y, …): cumulative distribution function. -
quantile(par, p, …): quantile function. -
support(par, …): lower and upper response support. It must return a numeric matrix with columnsminandmax, with one row per parameter combination. A fixed numeric vectorc(min, max)can be supplied instead. If support is omitted but a quantile function is available, it is inferred from the quantiles at probabilities zero and one. Without either, a discrete family is assumed to be a count distribution on the nonnegative integers. For a continuous family with a density, a parameter-dependent support is inferred heuristically by probing the log-density. Explicit support is recommended for shifted or truncated counts and for unusual, disconnected, or numerically extreme continuous supports. -
random(par, n, …): random number generator. -
mean(par, …): mean function. -
variance(par, …): variance function. -
skewness(par, …),kurtosis(par, …): optional higher-order moment functions. -
initialize: a named list of initialization functions, one for each parameter (e.g.,initialize$mu(y, …)), used to generate starting values. -
valid.response(x): a function that checks whether the response is valid (e.g., numeric, non-factor). -
optimizer(): an optional function defining a custom optimization method for use withgamlss2, see alsoRS.
If analytical score or Hessian functions are not provided, they are approximated numerically. If cdf() is omitted for a family with an available density and support, it is constructed by numerical integration for a continuous family or by summing the probability mass over the integer support for a discrete family. Analytical CDFs are retained unchanged. Numerical CDFs accept lower.tail, log.p, rel.tol, and abs.tol; continuous CDFs additionally accept subdivisions, while discrete CDFs accept max.terms.
If quantile() is omitted but a CDF and support are available, it is constructed by numerical inversion. Continuous CDFs are inverted by bracketing and bisection. For discrete families, the generated function returns the smallest integer whose CDF is at least the requested probability. Numerical quantiles accept lower.tail, log.p, tol, and maxiter; quantiles generated directly from a count PMF additionally accept max.terms. Analytical quantile functions are retained unchanged.
If mean() or variance() is omitted for a family with a density and support, the missing functions are constructed numerically. Continuous moments use adaptive integration; the variance is integrated as a centered second moment. Discrete moments are obtained by batched summation over the integer support. Numerical moments accept rel.tol and abs.tol; continuous moments additionally accept subdivisions, while discrete moments accept max.terms. A failure to establish a finite moment is reported as an error. Analytical moment functions are retained unchanged.
If random() is omitted but a quantile function is available, random values are generated by inverse transform sampling. Its … arguments are forwarded to the quantile function. Existing random-number generators are retained unchanged.
Value
For family.gamlss2(), the family object stored in a fitted “gamlss2” model is returned. A “gamlss2.family” object itself is a list describing the response distribution, its parameters, their link functions, and the distribution-specific methods used during fitting, prediction, diagnostics, and simulation.
See Also
gamlss2