library("gamlss2")
data("faithful", package = "datasets")
## sort for plotting
d <- faithful[order(faithful$eruptions), ]
## construct and estimate a two-normal mixture
fam <- mixture_family(NO, k = 2)
## fit mixture pdf
fit_family(d$waiting, family = fam)
## mixture of experts model
m <- gamlss2(waiting ~ s(eruptions) | . | . | . | ., data = d, family = fam)
## estimated effects
plot(m)
## predict parameters
par <- predict(m)
head(par)
## prior and posterior component probabilities
prior <- family(m)$probabilities(par)
head(prior)
posterior <- family(m)$responsibilities(par, d$waiting)
head(posterior)
## classify observations by their largest
## posterior probability and show estimated quantiles.
fit <- quantile(m)
pc <- max.col(as.matrix(posterior))
ccols <- c(2, 4)
plot(
waiting ~ eruptions,
data = d,
col = ccols[pc],
pch = 19,
xlab = "Eruption duration (minutes)",
ylab = "Waiting time (minutes)",
ylim = range(d$waiting, fit)
)
matlines(d$eruptions, fit, col = 1, lty = 1)
legend(
"topleft", paste("Component", 1:2),
col = ccols, pch = 19, bty = "n"
)
## mixture regression with a numeric and a binary predictor.
data("mtcars", package = "datasets")
cars <- transform(
mtcars,
am = factor(am, levels = 0:1, labels = c("automatic", "manual"))
)
## both component means depend on weight and transmission,
## component scales and the mixing probability are constant
m <- gamlss2(
mpg ~ s(wt) + am | 1 | s(wt) + am | 1 | 1,
data = cars,
family = fam
)
cars_par <- predict(m)
cars_posterior <- family(m)$responsibilities(cars_par, cars$mpg)
cars_component <- max.col(as.matrix(cars_posterior))
table(cars$am, cars_component)
## plot estimated classes and component mean regressions,
## point shape and line type distinguish automatic
## from manual transmission.
wt_grid <- seq(min(cars$wt), max(cars$wt), length.out = 100)
newcars <- expand.grid(wt = wt_grid, am = levels(cars$am))
newcars$am <- factor(newcars$am, levels = levels(cars$am))
cars_fit <- predict(m, newdata = newcars)
plot(
mpg ~ wt,
data = cars,
col = ccols[cars_component],
pch = c(1, 16)[cars$am],
xlab = "Weight (1000 lbs)",
ylab = "Miles per gallon",
ylim = range(cars$mpg, cars_fit)
)
for(j in 1:2) {
for(a in seq_along(levels(cars$am))) {
keep <- newcars$am == levels(cars$am)[a]
lines(
newcars$wt[keep], cars_fit[keep, paste0("c", j, ".mu")],
col = ccols[j], lty = a, lwd = 2
)
}
}
legend(
"topright", paste("Component", 1:2),
col = ccols, lty = 1, lwd = 2, bty = "n"
)
legend(
"bottomleft", levels(cars$am),
col = 1, pch = c(1, 16), lty = 1:2, bty = "n"
)Construct a Finite Mixture Distribution Family
Description
Constructs a “gamlss2.family” representing a finite mixture of two or more component distribution families. Component distribution parameters and the mixing probabilities may all be modeled using separate additive predictors.
Usage
mixture_family(families = NO, k = NULL, reference = 1L,
prefix = "c", initialize = TRUE)
Arguments
families
|
A single distribution family or a list of component families. Each entry may be a family-generating function, a family object, or a character family name accepted by the gamlss2 family infrastructure. If a single family is supplied, it is repeated |
k
|
Number of mixture components. For a single family, |
reference
|
Integer identifying the reference component for the mixing-weight parameterization. Its unnormalized odds are fixed at one. |
prefix
|
Either one non-empty character prefix, from which component labels such as |
initialize
|
Starting-value strategy. |
Details
For component densities \(f_j(y)\) and probabilities \(\pi_j\), the mixture density is
\(f(y) = \sum_{j=1}^k \pi_j f_j(y).\)
It is evaluated using a log-sum-exp calculation to avoid numerical underflow.
Component parameters are ordered component by component. For example, a two-component normal mixture has parameters c1.mu, c1.sigma, c2.mu, c2.sigma, and pi2. Consequently, a gamlss2 formula for this family has five right-hand sides in that order.
Mixing probabilities use a baseline-logit parameterization analogous to MN. For every non-reference component \(j\), the parameter pij has a log link and represents positive odds \(w_j\) relative to the reference component. With the reference odds fixed at one,
\(\pi_j = \frac{w_j}{1 + \sum_{\ell \ne r} w_\ell}, \qquad \pi_r = \frac{1}{1 + \sum_{\ell \ne r} w_\ell}.\)
The probabilities() family function converts these odds to normalized probabilities.
Analytical first derivatives are assembled from the component-family scores. If
\(r_{ij} = \frac{\pi_{ij}f_j(y_i)}{\sum_\ell \pi_{i\ell}f_\ell(y_i)}\)
denotes a posterior responsibility, the score for a component parameter is its component score multiplied by \(r_{ij}\). The mixing-predictor score is \(r_{ij} - \pi_{ij}\). Missing Hessians are supplied numerically by the standard family completion infrastructure.
The family supplies a fused working-response update for estimation. It uses the same finite-difference diagonal Hessian as the standard family completion infrastructure, but shares the component densities at the current parameters. For a component-parameter perturbation, only that component’s density is recomputed; mixing-parameter perturbations reuse every component density. Repeated identity and log links are mapped directly, and a single family repeated k times is completed only once. These optimizations do not change the numerical updates.
The default initialize = “separated” strategy places recognized component location parameters at ordered empirical midpoint quantiles. Other parameters retain their component-family initializers and mixing probabilities start equally. This deliberately broad initialization is robust for flexible mixture regressions, where response-only clusters do not account for the model’s covariates and can select a poor local mode.
With initialize = “cluster”, empirical midpoint quantiles first define ordered centers for a finite numeric response. One-dimensional Lloyd updates assign observations to their nearest center; a rank partition is used if tied centers would otherwise produce an empty cluster. Each component family’s initialization functions are then evaluated on the response observations assigned to that component, and the resulting values are reduced to component-level starts. Mixing probabilities are initialized from cluster proportions and converted to odds relative to the reference component. This can be useful for intercept-only or otherwise simple mixtures, but is not uniformly better for mixture regression.
With initialize = “modal”, sufficiently separated modes of a kernel density estimate initialize continuous component locations and scales, with mixing probabilities based on nearest-mode assignments. If suitable modes cannot be found, this falls back to separated starts. When modes are found, log-linked sigma parameters use a shifted log link with lower bound 0.8 times the density bandwidth, and log-linked JSU tau parameters have lower bound one. These bounds remain in effect during fitting: this option changes the admissible parameter space, not just the starting values. Other strategies retain the component-family links without these bounds.
Every proposed natural-scale value is checked after applying its link. Values on a link boundary are shrunk toward the inverse-link value at zero using one pseudo-observation. If a cluster-specific value remains invalid, initialization falls back to the corresponding full-response component initializer and finally to the inverse-link value at zero. Non-numeric responses or responses that cannot be partitioned use these fallback values and equal mixing probabilities.
The returned family includes cdf(), quantile(), random(), mean(), and variance() whenever all component families provide the corresponding operation. All components must use the same distribution type, for example all continuous or all discrete. Mixture quantiles are obtained by numerically inverting the mixture CDF.
As with finite mixtures generally, the likelihood can have local optima and equivalent solutions obtained by exchanging component labels. Separated initialization helps but does not impose an ordering constraint. Multiple fits, including fits using the optional “cluster” strategy or explicitly supplied starting values, can still be useful for difficult data.
Value
An object of class “gamlss2.family”. In addition to the standard family elements, it contains:
-
probabilities(par, n = NULL): normalized prior mixing probabilities. -
responsibilities(par, y): posterior component probabilities conditional on observations. -
score_matrix(par, y): observation-wise scores for all mixture parameters, evaluated with shared responsibilities. -
components: the completed component families. -
component_map: mappings from original component parameter names to mixture parameter names.
See Also
gamlss2, gamlss2.family