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,
criterion = "aicc"
)
## estimated effects
plot(m)
## predict parameters
par <- predict(m, type = "parameter")
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 quantiles of the fitted distributions
pc <- max.col(as.matrix(posterior))
ccols <- c(2, 4)
fq <- quantile(m)
plot(
waiting ~ eruptions,
data = d,
col = ccols[pc],
pch = 19,
xlab = "Eruption duration (minutes)",
ylab = "Waiting time (minutes)",
ylim = range(d$waiting, fq)
)
matlines(d$eruptions, fq, 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, type = "parameter")
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))
fit <- predict(m, newdata = newcars, type = "parameter")
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, 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], 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.
Component scores are weighted by posterior responsibilities; missing Hessians are computed numerically. Working-response updates reuse component densities.
The default initialize = “separated” starts recognized locations at ordered empirical midpoint quantiles, uses component initializers for other parameters, and starts with equal mixing probabilities. “cluster” starts from response-based groups and their proportions; “component” uses unmodified component initializers and equal probabilities. “modal” uses separated density modes and their nearest observations, falling back to separated starts if no suitable modes are found.
When modal starts are found, log-linked sigma has a lower bound of 0.8 times the density bandwidth and log-linked JSU tau has a lower bound of one. These bounds remain in effect during fitting. Invalid or unusable response-based starts fall back to full-response component initializers and then to the inverse-link value at zero.
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