mixture(fam1, fam2, ...) builds a K-component mixture: each
component keeps its own distributional parameters, suffixed by the
component index (mu1, sigma1, mu2, ...), and the mixing
proportions come from theta1 ... theta{K-1} (multinomial logit
against the last component, each with its own linear predictor, so
mixing weights may depend on covariates). The main model formula
applies to every component mean; override per component with
bf(y ~ x, mu2 ~ 1).
Arguments
- ...
Two or more component families.
- groups
Optional one-sided formula naming the latent-class grouping factor.
- order
brms's argument.
NULL,"none"orFALSEleaves the components unordered, which is what a fit here does."mu"orTRUE, brms's ordering of themuintercepts, is refused, and any other value is refused as invalid, in brms's words. With ordinal components,"mu"orTRUEgives the components one shared threshold vector, as in brms (see "Ordinal components").
Details
A formula for a mixing weight follows brms's rule. Write one for all
of theta1 ... thetaK but one, and the one left out is the reference
component, whose linear predictor is 0: bf(y ~ 1, theta2 ~ x) on
two components makes component 1 the reference, so theta2's
coefficients are the log odds of component 2 against component 1.
A formula for fewer of them is refused, as brms refuses it.
A mixing weight's RESPONSE scale is the softmax over the component
predictors, so frm_linpred(type = "response", dpar = "theta1") is a
probability while type = "link" stays the predictor the density
works on. Under se.fit = TRUE that probability's standard error is
the delta method through its OWN predictor, p (1 - p) times the
predictor's standard error. For two components that is exact. For
three or more the softmax also moves with the other components'
predictors, and those terms are dropped, so the standard error is
CONSERVATIVE: measured 5.5% to 26.1% wider than the joint delta
method on a three-component fit, never narrower.
The likelihood is a parameter-branch-free logsumexp, so Laplace
machinery is untouched; the usual finite-mixture ML caveats apply
instead: the likelihood is invariant to component relabeling, so the
component means are initialized on spread-out response quantiles,
and multimodality is real (compare starts, or order the intercepts
with bounds: set_prior("", class = "Intercept", dpar = ..., lb = )
per component).
Ordinal components
The components can be ordinal, as in brms: cumulative(),
sratio(), cratio() and acat() in any combination, each with its
own link and threshold structure, or hurdle_cumulative() for every
component. Component k has the latent predictor mu<k> and the
discrimination disc<k> (held at 1 unless its formula models it;
hu<k> too for the hurdle family). The category probabilities of
the mixture are the theta-weighted sum of the components'
category probabilities, and fitted(), predict() and simulate()
use that sum.
order sets the thresholds, as in brms. "none", the default,
gives each component its own thresholds, reported as
b_mu1_Intercept[k], b_mu2_Intercept[k], and addressed by
set_prior(class = "Intercept", dpar = "mu1"). "mu" (or TRUE)
gives all components one threshold vector, brms's
fixed_Intercept: it is ordered if one component is
cumulative() or hurdle_cumulative(), a sum-to-zero component
centers it, and set_prior(class = "Intercept") with no dpar
addresses it. "mu" cannot be combined with
threshold = "equidistant", as in brms. thres(x = ),
thres(gr = ) and cs() work as for one ordinal family, per
component; a cumulative() or hurdle_cumulative() component with
cs() warns once, as brms warns.
The likelihood does not change when two components change places, so
the labels follow the start values. A component has no mu
intercept to start apart, so component k starts with its
thresholds shifted by -F^-1(k / (K + 1)) on its link's scale, which
gives component 1 the lowest latent location; under order = "mu"
the slopes of mu<k> start apart instead. groups = is not
available with ordinal components. Three designs have a flat
direction in the likelihood, and the fit warns about each while the
parameters no prior holds outnumber what the data identify (a
threshold prior holds none of the hurdle's hu<k> or mixing
weights): hurdle components whose
hu<k> and mixing weights have no predictors; a mixture with no
predictor in any distributional parameter, where every row has one
category distribution; and order = "mu" with no
predictor in any mu<k>, when the components are the same family
with the same link and the same held disc, so that they are one
distribution.
Maximum likelihood for an ordinal mixture often has its supremum at
a degenerate boundary, where a component turns into a step function
of its predictors: its slopes and threshold distances run off, and
their standard errors mean nothing. In 160 fits of a two-component
mixture to simulated data with two well-separated classes (300 or
500 rows), 22 ended there, and in 16 of the 22 the degenerate point
had a higher log-likelihood than the best non-degenerate optimum
found from other starts, so a better start does not move such a fit
off it. The fit warns when a component ends at such a point: when
doubling all its latent distances (or, with cs(), one threshold's)
costs the log-likelihood less than 0.1, or when no row lies within a
latent distance of 50 of one of its thresholds. Strong predictors
alone do not set it off: none of 167 sound fits with slopes up to 15
warned, under the logit, probit, cloglog and cauchit. Compare
several starts (frm(start = )), the more so with three or more
components, and hold the components with priors on class = "Intercept" and class = "b" with dpar = "mu<k>", as a brms fit
does; or use fewer components. An ordinal mixture holds each
interior category's probability at exp(-692) or above, so that a
component which stops using a category does not stop the optimizer
with an undefined gradient. This floor binds only where two
thresholds of a component are the same double, and the fit warns
when it ends there.
With groups = ~g the mixture moves to the group level (latent
classes): every observation of a group shares one class draw, and
the marginal likelihood sums the class assignment per group.
Continuous random effects, smooths, and gp() terms are allowed in
the component formulas - the class sum happens conditional on the
latent effects, so one Laplace approximation integrates them
(growth-mixture models). Random effects written in a component
formula are class-specific by construction; the Laplace
approximation of the class-mixture integrand is not exact even for
gaussian responses (a fraction of a log-likelihood unit in typical
well-separated problems). quadrature = TRUE makes the integral
numerically exact when the per-group integrand is univariate (one
scalar random intercept, in one class); with class-specific
intercepts in several classes the coordinates couple and quadrature
remains approximate - use frmtmb.sample::check_laplace() to judge.
Mixing-weight predictors are evaluated at each group's first row
(use group-constant covariates). mixture_probs() returns the
posterior class probabilities per group (or per observation for
ordinary mixtures), conditional on the random-effect modes.
A mixture of multivariate normal components over a matrix response,
with one class per row (model-based clustering in the manner of
mclust), is mixture_mvn() in the extension package frmtmb.latent:
run library(frmtmb.latent) first.
Examples
# two well-separated gaussian components
set.seed(3)
dd <- data.frame(y = c(rnorm(80, 0, 1), rnorm(80, 5, 1)),
x = rnorm(160))
fit <- frm(bf(y ~ 1) + mixture(gaussian(), gaussian()), data = dd)
# one mu and sigma per component, plus the mixing weight theta1
fixef(fit)
#> Estimate Est.Error Q2.5 Q97.5
#> mu1_Intercept -0.002346141 0.09962412 -0.1976058 0.1929135
#> mu2_Intercept 4.999986512 0.11450977 4.7755515 5.2244215
# posterior class probability per observation
head(mixture_probs(fit))
#> class1 class2
#> [1,] 1.0000000 4.288987e-08
#> [2,] 0.9999990 1.000118e-06
#> [3,] 0.9999851 1.488386e-05
#> [4,] 1.0000000 1.798825e-08
#> [5,] 0.9999891 1.087896e-05
#> [6,] 0.9999952 4.800314e-06
# the mixing weight can take its own predictor
frm(bf(y ~ 1, theta1 ~ x) + mixture(gaussian(), gaussian()), data = dd)
#> Family: mixture(gaussian, gaussian)
#> Links: mu1 = identity; sigma1 = log; mu2 = identity; sigma2 = log;
#> theta1 = identity
#>
#> Formula: y ~ 1
#> Data: dd (Number of observations: 160)
#> Method: ML logLik: -327.164 AIC: 666.328 BIC: 684.779
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> mu1_Intercept 0.00 0.10 -0.20 0.19 -0.03 0.98
#> mu2_Intercept 5.00 0.11 4.78 5.22 43.65 <2e-16
#> theta1_Intercept 0.01 0.16 -0.30 0.32 0.06 0.95
#> theta1_x -0.07 0.15 -0.36 0.22 -0.47 0.64
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI
#> sigma1 0.88 0.07 0.75 1.03
#> sigma2 1.01 0.08 0.86 1.19
# ordinal components, each with its own thresholds
set.seed(5)
do <- data.frame(x = rnorm(300))
cls <- rbinom(300, 1, 0.4)
lat <- ifelse(cls == 1, 1.5 * do$x + 1, -do$x - 1) + rlogis(300)
do$y <- findInterval(lat, c(-1.5, 0, 1.5)) + 1
fo <- frm(bf(y ~ x), family = mixture(cumulative(), cumulative()),
data = do)
fixef(fo)
#> Estimate Est.Error Q2.5 Q97.5
#> mu1_Intercept[1] -0.7001189 0.1949475 -1.082208963 -0.3180289
#> mu1_Intercept[2] 0.5564518 0.2793007 0.009032619 1.1038711
#> mu1_Intercept[3] 2.4378557 0.4870638 1.483228307 3.3924832
#> mu2_Intercept[1] -4.1652485 2.7658915 -9.586296181 1.2557991
#> mu2_Intercept[2] -1.0120801 1.0328723 -3.036472624 1.0123125
#> mu2_Intercept[3] -0.1683368 0.9262030 -1.983661228 1.6469876
#> mu1_x -0.7099426 0.2144915 -1.130338235 -0.2895469
#> mu2_x 2.9076463 2.0421557 -1.094905383 6.9101980
head(fitted(fo)[, "Estimate", ])
#> P(Y = 1) P(Y = 2) P(Y = 3) P(Y = 4)
#> [1,] 0.1930635 0.4059203 0.27914618 0.1218700
#> [2,] 0.3743300 0.1682677 0.09804975 0.3593526
#> [3,] 0.2396042 0.3548756 0.28002949 0.1254908
#> [4,] 0.2293582 0.2739112 0.24290635 0.2538242
#> [5,] 0.4108600 0.1509122 0.07995397 0.3582738
#> [6,] 0.1887376 0.3932082 0.28202289 0.1360313
# \donttest{
# latent classes: every observation of a group shares one class
set.seed(4)
n_g <- 40
cls <- rep(c(1, 2), each = n_g / 2)
dg <- data.frame(g = factor(rep(seq_len(n_g), each = 5)))
dg$y <- rnorm(nrow(dg), c(0, 4)[cls[as.integer(dg$g)]], 1)
fg <- frm(bf(y ~ 1) + mixture(gaussian(), gaussian(), groups = ~g),
data = dg)
head(mixture_probs(fg)) # one row per group, not per observation
#> class1 class2
#> 1 1 1.262122e-12
#> 2 1 2.820791e-12
#> 3 1 8.436732e-16
#> 4 1 7.884398e-16
#> 5 1 4.031519e-09
#> 6 1 1.439555e-13
# }