library("gamlss2")
Normal <- function(...) {
fam <- list(
"family" = "Normal",
"names" = c("mu", "sigma"),
"links" = c("mu" = "identity", "sigma" = "log"),
"score" = list(
"mu" = function(y, par, ...) {
(y - par$mu) / (par$sigma^2)
},
"sigma" = function(y, par, ...) {
-1 + (y - par$mu)^2 / (par$sigma^2)
}
),
"hess" = list(
"mu" = function(y, par, ...) {
1 / (par$sigma^2)
},
"sigma" = function(y, par, ...) {
rep(2, length(y))
},
"mu.sigma" = function(y, par, ...) {
rep(0, length(y))
}
),
"loglik" = function(y, par, ...) {
sum(dnorm(y, par$mu, par$sigma, log = TRUE))
},
"mu" = function(par, ...) {
par$mu
},
"pdf" = function(y, par, log = FALSE) {
dnorm(y, mean = par$mu, sd = par$sigma, log = log)
},
"cdf" = function(y, par, ...) {
pnorm(y, mean = par$mu, sd = par$sigma, ...)
},
"random" = function(n, par) {
rnorm(n, mean = par$mu, sd = par$sigma)
},
"quantile" = function(p, par) {
qnorm(p, mean = par$mu, sd = par$sigma)
},
"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)
}
## load the abdominal circumference data
data("abdom", package = "gamlss.data")
## specify the model formula
f <- y ~ s(x) | s(x)
## estimate model
m1 <- gamlss2(f, data = abdom, family = Normal)
## plot estimated effects
plot(m1, which = "effects")
## plot diagnostics
plot(m1, which = "resid")
## predict parameters
par <- predict(m1)
## predict quantiles
pq <- quantile(m1, probs = c(0.05, 0.5, 0.95))
## visualize
plot(y ~ x, data = abdom, pch = 19,
col = rgb(0.1, 0.1, 0.1, alpha = 0.3))
matlines(abdom$x, 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(y, par, log = FALSE, ...) {
shape <- par$sigma
scale <- par$mu/par$sigma
dgamma(y, shape = shape, scale = scale, log = log)
},
"cdf" = function(y, par, 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(p, par, 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)
}
## example using the Munich rent data
data("rent", package = "gamlss.data")
## model formula
f <- R ~ ti(Fl) + ti(A) + ti(Fl, A, bs = "ps") |
ti(Fl) + ti(A) + ti(Fl, A, bs = "ps")
## estimate model
m2 <- gamlss2(f, data = rent, 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, ...)
Arguments
object
|
A fitted “gamlss2” model object.
|
…
|
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(y, par, 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]. -
hess: a named list of functions computing second derivatives (the negative Hessian). For parametersmuandsigma, this includes:hess[“mu”],hess[“sigma”], and optionally cross derivatives likehess[“mu.sigma”]. -
loglik(y, par, …): a function computing the total log-likelihood. -
cdf(y, par, …): cumulative distribution function. -
quantile(p, par, …): quantile function. -
random(n, par, …): 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 quantile residuals are required, a cdf() function must be available.
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