Family constructors without a stats::family equivalent, following
brms naming. gaussian(), poisson(), binomial(), and Gamma()
from 'stats' are accepted directly by frm() and +.
Usage
dirichlet(link = "logit", link_phi = "log", refcat = NULL)
logistic_normal(link = "identity", link_sigma = "log", refcat = NULL)
student(link = "identity", link_sigma = "log", link_nu = "logm1")
lognormal(link = "identity", link_sigma = "log")
negbinomial(link = "log", link_shape = "log")
nbinom1(link = "log", link_phi = "log")
Beta(link = "logit", link_phi = "log")
tweedie(link = "log", link_phi = "log")
compois(link = "log", link_nu = "log")
zero_inflated_poisson(link = "log", link_zi = "logit")
zero_inflated_negbinomial(link = "log", link_shape = "log", link_zi = "logit")
hurdle_poisson(link = "log", link_hu = "logit")
hurdle_negbinomial(link = "log", link_shape = "log", link_hu = "logit")
multinomial(link = "logit", refcat = NULL, K = NULL)
dirichlet_multinomial(
link = "logit",
link_phi = "log",
refcat = NULL,
K = NULL
)
cumulative(link = "logit", link_disc = "log", threshold = "flexible")
hurdle_cumulative(
link = "logit",
link_hu = "logit",
link_disc = "log",
threshold = "flexible"
)
beta_binomial(link = "logit", link_phi = "log")
zero_inflated_beta_binomial(
link = "logit",
link_phi = "log",
link_zi = "logit"
)
skew_normal(link = "identity", link_sigma = "log", link_alpha = "identity")
exgaussian(link = "identity", link_sigma = "log", link_beta = "log")
bernoulli(link = "logit")
geometric(link = "log")
exponential(link = "log")
weibull(link = "log", link_shape = "log")
shifted_lognormal(link = "identity", link_sigma = "log", link_ndt = "log")
hurdle_gamma(link = "log", link_shape = "log", link_hu = "logit")
hurdle_lognormal(link = "identity", link_sigma = "log", link_hu = "logit")
zero_inflated_binomial(link = "logit", link_zi = "logit")
zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit")
zero_one_inflated_beta(
link = "logit",
link_phi = "log",
link_zoi = "logit",
link_coi = "logit"
)
xbeta(link = "logit", link_phi = "log", link_kappa = "log")
asym_laplace(link = "identity", link_sigma = "log", link_quantile = "logit")
zero_inflated_asym_laplace(
link = "identity",
link_sigma = "log",
link_quantile = "logit",
link_zi = "logit"
)
huber(link = "identity", k = 1.345, link_sigma = "log")
sratio(link = "logit", link_disc = "log", threshold = "flexible")
cratio(link = "logit", link_disc = "log", threshold = "flexible")
acat(link = "logit", link_disc = "log", threshold = "flexible")
von_mises(link = "tan_half", link_kappa = "log")
categorical(link = "logit", refcat = NULL, levels = NULL, K = NULL)Arguments
- link
Link for
mu. See frmtmb-links.- refcat
For
categorical(),multinomial(),dirichlet_multinomial(),dirichlet()andlogistic_normal(): the name of the reference category, whose predictor is held at zero.NULL(the default) takes the first category: the first level of a factor response, or the first column of a matrix response.NApredicts every category, which every one of them butlogistic_normal()allows; see "Categorical (nominal) responses" for what a fit can then identify.- link_sigma, link_shape, link_phi, link_kappa, link_ndt, link_beta, link_disc
Link for a strictly positive parameter: one of
"log"(the default),"identity","softplus"or"squareplus".- link_nu
Link for
nu.student()'s degrees of freedom take"logm1"(the default) or"identity", which keeps them above one;compois()'s dispersion is an ordinary positive parameter and takes the positive set.- link_zi, link_hu, link_quantile, link_zoi, link_coi
Link for a parameter on the unit interval:
"logit"(the default) or"identity". Inzero_one_inflated_beta(),zoiis the probability of an exact 0 or 1, andcoiis the probability that such a value is 1.- K
For
multinomial(),dirichlet_multinomial()andcategorical(): the number of response categories, for a family built away from the data. The categories are then1..Kand the dparsmu2 ... muK. Leave itNULLandfrm()reads the categories off the response, as brms does: the columns of a multinomial count matrix, the levels of a categorical factor.- threshold
For the ordinal families: the threshold structure, as in brms.
"flexible"(the default) estimates every threshold."equidistant"estimates the first threshold anddelta, the distance between neighboring thresholds."sum_to_zero"holds the thresholds to a sum of zero. See Ordinal thresholds and discrimination.- link_alpha
Link for
skew_normal()'s skewness, which is signed:"identity"(the default),"log","softplus"or"squareplus".- k
For
huber(): Huber's tuning constant, the residual size in units ofsigmawhere the density stops being gaussian and becomes Laplace. It is FIXED, not estimated - the default 1.345 isMASS::rlm()'s, which gives 95% efficiency against a gaussian.- levels
For
categorical(): the response's category labels, in the order that fixes the dpar names. Only needed when the family is built away from the data;frm()reads them off the response.
Details
An ordinal family (cumulative(), sratio(), cratio(), acat())
takes the response's level order as the category order. Supply an
ordered factor, or integer codes 1..K. An unordered factor is
refused, as brms refuses it, because its level order is alphabetical
unless someone set it, and that order is the model.
A bernoulli() response can hold any two values, which are coded 0
and 1 by level order as brms codes them: the levels of a factor, the
sorted values of a number or a character vector (so -2 and -1
become 0 and 1, and so do 1 and 2), and FALSE and TRUE of a
logical. A number with one value is the 1 of 0 and 1 unless the value
is 0. The fit keeps the coding, and codes newdata's response with it.
A third value is refused with brms's message. Two values that both
lie strictly between 0 and 1, such as 0.1 and 0.9, are coded too, as
brms codes them, with a warning: they are usually proportions, which
binomial() with trials() or Beta() fits. simulate() and the
draws return the codes 0 and 1, as brms's posterior_predict() does.
A response with only two outcomes gets brms's message suggesting
bernoulli(): an ordinal or categorical response with two
categories, and a binomial(), beta_binomial() or
zero_inflated_binomial() response whose trials are all one.
Every constructor has brms's fields: $link is the name of the link
for the mean, and $link_<dpar> the link of each other parameter,
as in beta_binomial()$link_phi. See frmtmb_family().
Compositional responses
dirichlet() and logistic_normal() model a row of proportions that
sums to one, brms's simplex families. Write the response as an
n x K matrix: cbind(y1, y2, y3), or a matrix column of the data
(dd$Y <- Y). Every row must sum to one, to within 1.5e-8, and every
proportion must lie strictly between 0 and 1. The column names are
the categories, and the FIRST is the reference category unless
refcat names another, as in brms.
Each non-reference category k gets its own linear predictor
mu<k>, and the main formula applies to every one of them;
bf(Y ~ x, muc ~ z) gives category c a predictor of its own. A
group term (1 | g) gives each category its own effects, and
(1 | ID | g) correlates them across categories, as for
categorical().
dirichlet() is y ~ Dirichlet(phi * softmax(mu)): the category
means are the softmax of the predictors, with the reference
category's held at zero, and phi is the precision, which
bf(phi ~ z) models. refcat = NA predicts every category, which
brms allows; without a prior the predictors are then identified
only up to a common shift, and the fit says so.
logistic_normal() makes the log ratios against the reference
category, log(y_k / y_ref), multivariate normal with means
mu<k>, standard deviations sigma<k> and correlations lncor.
Each sigma<k> takes its own formula; lncor is one correlation
matrix for the model, reported as lncor(b,c) in the summary, and
takes an lkj() prior under class = "lncor".
fitted() on a dirichlet() fit returns brms's n x 4 x K array of
category means, layers named P(Y = k). A logistic-normal mean has
no closed form, so fitted() and conditional_effects() are refused
on it, as in brms. On both, predict() and simulate() draw
compositions (rows on the simplex), frm_linpred(dpar = "mub") gives
one category's predictor, and residuals() is refused, as brms
refuses it for every family whose response is a set of categories.
brms also has dirichlet2, not exported there, reached as
brmsfamily("dirichlet2"): one predictor per category on the scale
of its Dirichlet concentration (log link), and no precision. The
same spelling reaches it here.
The negative binomial in sigma
brmsfamily("negbinomial2") is brms's negbinomial2: the negative
binomial of negbinomial() with the dispersion sigma = 1 / shape,
so the variance is mu + sigma * mu^2, and sigma near 0 is the
Poisson. brms does not export a constructor for it, and neither does
frmtmb. sigma takes the positive links (link_sigma, default
"log") and its own formula, bf(y ~ x, sigma ~ z). The family is
negbinomial() underneath: the same density, simulator, mean and
variance, and the same addition terms, so rate(d) multiplies the
mean and the shape 1 / sigma by d, as in brms. The fit is
negbinomial()'s with sigma in place of 1 / shape: the
log-likelihood is the same and only the parameter is reported
differently.
Hurdle and zero-one-inflated responses
A hurdle family gives every zero to hu, the hurdle probability, and
truncates its count or continuous part at zero. For
hurdle_negbinomial(), P(Y = 0) = hu and a positive count y has
probability (1 - hu) NB(y | mu, shape) / (1 - NB(0 | mu, shape)),
so mu is the mean of the untruncated negative binomial and not the
mean of the positive counts.
zero_one_inflated_beta() takes a response in [0, 1]. zoi is the
probability of an exact 0 or 1, coi is the probability that such a
value is 1, and a value strictly between 0 and 1 follows a beta
distribution with mean mu and precision phi. The expected
response, which fitted() returns, is zoi * coi + (1 - zoi) * mu.
coi is scored only on the values at exactly 0 or 1. A response with
0s but no 1 (or the reverse) sends coi to that boundary, and one
with neither gives coi no data at all, so that no standard error of
the fit is usable; the fit warns in both cases. brms fits such data
through its prior on coi. Here, hold coi at a value with
bf(coi = 0.5).
zero_inflated_beta_binomial() is beta_binomial() with a
zero-inflation probability zi: P(Y = 0) = zi + (1 - zi) BB(0),
and a count y > 0 has probability (1 - zi) BB(y), with the trials
from trials(). fitted() returns (1 - zi) * mu * trials.
hurdle_cumulative() is brms's hurdle over an ordinal response. The
response is 0..K: 0 is the hurdle, with probability hu, and the
categories 1..K follow cumulative() with probabilities scaled by
1 - hu. An ordered factor takes its FIRST level as the category 0.
fitted() returns the K + 1 category probabilities, 0 first. The
expected category that conditional_effects() shows by default is
scored by these codes, so the hurdle scores 0; brms scores each
column by its position and reads one higher. disc and threshold
work as for the other ordinal families (see Ordinal thresholds and
discrimination). thres(x = ) and thres(gr = ) count the
categories above the hurdle, and thres(gr = ) keeps the hurdle in
each group's density, as brms's hurdle_cumulative_*_merged_lpmf
does. cs() takes the category-specific offsets off each row's
thresholds, as in brms; they can make two thresholds cross for some
rows, and the probability of the category between them is then
negative. The density is then NaN, as brms's is, and so are the
row's fitted() probabilities, where brms returns the negative
difference; simulate() gives NA for the row. cs() and
thres(gr = ) together are refused, as in brms.
Ordinal thresholds and discrimination
Every ordinal family (cumulative(), sratio(), cratio(),
acat() and hurdle_cumulative()) has brms's discrimination
parameter disc. The distribution function is read at
disc * (tau_k - eta) (at disc * (eta - tau_k) for cratio() and
acat()), with disc held at 1 unless the formula models it, as in
brms: bf(y ~ x, disc ~ z). Its link is link_disc, "log" by
default. The likelihood cannot tell an intercept in disc apart from
the scale of the thresholds, so write disc ~ 0 + z. brms accepts
disc ~ 1 + z and holds the intercept with its default
normal(0, 1) prior; frmtmb fits it too, and warns unless a prior
holds the intercept, as it does for cell-means columns that add up to
one (disc ~ 0 + h). Held at 1, disc shows nowhere in the output,
as in brms: not on the Links: line, not as a fixed parameter, and
not in fixef(flatten = TRUE) or coef().
threshold selects brms's threshold structure. "flexible"
estimates each threshold. "equidistant" estimates the first
threshold and the distance delta between neighbors, so that
tau_k = tau_1 + (k - 1) * delta; variables() and summary() name
it delta as brms does, and set_prior(class = "delta") gives it a
prior, while class "Intercept" addresses the first threshold only.
delta is positive for cumulative() and hurdle_cumulative(),
whose thresholds are ordered. "sum_to_zero" holds the thresholds
to a sum of zero and does not center the design, as in brms. brms
declares one parameter per threshold and relies on its prior for the
common location, which the likelihood does not see; frmtmb estimates
the K - 2 free directions instead, so class "Intercept" has
nothing to address there and is refused. thres(gr = ) gives each
level its own first threshold and delta, or its own zero sum.
"equidistant" needs two thresholds per vector, since delta has
nothing to measure on one. "sum_to_zero" holds a single threshold
at zero, a vector with no parameter, which brms runs under
thres(gr = ) and frmtmb also fits without it.
Extended-support beta
xbeta() is the extended-support beta of Kosmidis and Zeileis
(2024), brms's xbeta, for a response in [0, 1] with exact 0s and
1s. A latent beta variable with mean mu and precision phi is
stretched by the exceedance kappa to (1 + 2 kappa) Z - kappa and
censored at 0 and 1, so the ends carry the latent mass beyond them.
mu is the mean of the latent variable, not of the response;
fitted() returns the mean of the response, brms's
posterior_epred(). The mass at each end is an incomplete beta
function, which is taken mostly from a continued fraction rather
than from RTMB::pbeta(), whose third derivatives are not finite
everywhere and which therefore failed the Laplace approximation of a
model with random effects.
kappa is placed mostly by the rows at exactly 0 or 1. With none,
only the shape of the interior places it, and it can run to 0,
where the model is Beta(), or grow without bound. When it runs to
0, or a coefficient of a log-link kappa has a standard
error above 10 or none, the fit warns; brms fits such data through its
gamma(0.01, 0.01) prior on kappa. Compare the fit with Beta(),
or hold kappa at a value with bf(kappa = ).
Categorical (nominal) responses
categorical() fits a multinomial logit to an unordered factor. One
category is the reference: its linear predictor is held at zero, and
each other category gets its own predictor named mu<Level>, as in
brms. The reference is the FIRST level, or the level that refcat
names: categorical(refcat = "b") gives the dpars mua and muc.
Changing the reference changes the coefficients, not the model: the
log-likelihood and the category probabilities stay the same. The main
model formula applies to every predictor, and any single category is
overridden by naming it: bf(y ~ x, muc ~ z) gives category "c"
its own predictor and leaves the rest on ~ x. A dpar name is
brms's: dots and underscores in a level are dropped, so the level
"p.q" gives mupq.
refcat = NA gives every category a predictor, as brms allows. The
likelihood then stays the same when one number is added to every
predictor, so a maximum-likelihood fit cannot determine the
coefficients, only the category probabilities, and the fit warns
that their standard errors are not available. A prior on each
predictor (set_prior()) determines them. brms leans on its priors
here too, and with its default priors (flat on the slopes) its
sampler does not converge on this model.
The categories come from the data, so frm() reads them off the
response before it parses the formula. Constructing the family away
from a data set (to inspect it, or to reach the parser through
another entry point) needs them stated: categorical(levels = c("a", "b", "c")), or categorical(K = 3) for a response already coded
1..K, which names the dpars mu2 ... muK. A character or logical
response is coerced to a factor with a message naming the level
order, because that order is the model.
fitted() returns brms's n x 4 x K array of category
probabilities, one layer per category named P(Y = <level>), and
frm_linpred(type = "response") the n x K matrix of the estimates,
rows summing to one: the convention the ordinal families follow.
fitted(scale = "linear") stacks the predictors into an
n x 4 x (K - 1) array, layers eta1, eta2, ... (K layers under
refcat = NA), as brms does; with one predictor (two categories and
a reference) it stays an n x 4 matrix, as in brms. The intervals of
fitted() here and for multinomial() are Wald intervals, the
estimate plus or minus a normal quantile times its delta-method
standard error, held inside the range of the quantity (0 to 1 for a
probability, 0 to the trials for a count), so near the edge of that
range a bound can be exactly 0. frm_linpred(dpar = "mub") gives one
category's predictor, which is where se.fit lives. simulate()
draws factor levels.
multinomial() is the same likelihood on an n x K matrix of counts,
with the number of trials in trials(), as in brms. The columns are
the categories: cbind(y1, y2, y3) gives the dpars muy2 and muy3,
and an unnamed matrix numbers them mu2, mu3. refcat works as
for categorical(). fitted() returns brms's n x 4 x K array of
expected COUNTS, probability times trials, which brms labels
P(Y = <category>) too. multinomial(K = 3) builds the family away
from the data, with the categories 1..K. A one-hot matrix with one
trial a row gives the log-likelihood of categorical().
dirichlet_multinomial() is brms's overdispersed multinomial. The
category probabilities of each row are Dirichlet with mean
softmax(mu) and precision phi, and the counts are multinomial
given them. The expected counts are those of multinomial(), and
the covariance of a row is the multinomial one times
(T + phi) / (1 + phi) at T trials, so a large phi is close to
multinomial(). phi takes a log link and its own formula,
bf(Y | trials(n) ~ x, phi ~ z). Everything else is as for
multinomial(): the response matrix, trials(), refcat,
fitted(), K and simulate(), which draws the probabilities of
each row from the Dirichlet and then the counts. With two categories
it is beta_binomial() with phi the same precision.
Circular responses
von_mises() models an angle in radians on (-pi, pi]. Its mu is
the mean direction and takes the tan_half link, which maps the
whole line onto that interval; kappa is the concentration and takes
a log link, with kappa = 0 the uniform distribution on the circle.
Both are brms's choices. fitted() and frm_linpred(type = "response")
report the mean direction. The normalizing constant needs
log I0(kappa), which RTMB differentiates exactly through its own
besselI method, so nothing here is a series approximation.
Residuals are differences of angles and are NOT wrapped, so read
residuals() on a von Mises fit with that in mind.
Quantile regression inference
asym_laplace() and zero_inflated_asym_laplace() fit quantile
regression through a WORKING likelihood: at a fixed quantile the
point estimates are consistent quantile estimates (they match
quantreg::rq()), but the asymmetric Laplace is not the data's
true density, so Wald standard errors and confint() intervals
computed from it are not calibrated. This is a property of the
asymmetric-Laplace approach, shared by brms. Use
frm_bootstrap() for intervals you can defend. The check
function's kink can also produce a benign false-convergence
warning near the optimum; frm_allfit() confirms the fit when in
doubt.
Degrees of freedom that run off
student() estimates nu on the logm1 link, which holds it above
one but puts no ceiling on it. A student-t reaches gaussian() only
in the limit nu -> Inf, so on data with no heavy tails the
likelihood keeps rising as nu grows and the maximum is never
attained. The fit then reports whatever nu the optimizer last
reached, with a standard error to match: values of 1e9 and above
are ordinary, and two runs of the same model on the same data in a
different row order can differ by orders of magnitude in nu while
agreeing on every coefficient to the last bit.
That is the model saying the data shows no heavy tails, not a
failure to converge. diagnose() names it under "Distributional
parameter at the end of its link". Read the fit as the gaussian one
it has become. If you want a number you can report, refit with
gaussian(), or hold nu somewhere finite with a prior (see
set_prior()); brms does the same with its default
gamma(2, 0.1).
nu runs off the OTHER end too, and there the reading is opposite.
On tails heavier than any identified nu can hold, Cauchy data for
instance, nu goes down to one instead: log(nu - 1) around -20,
a standard error in the thousands, and a natural-scale value that
prints as 1. diagnose() names that under the same heading. The
remedy is not gaussian(), which is the worst fit available for
such data. There the data is the message; a prior is still the way
to hold nu finite.
The likelihood itself stays accurate the whole way. The log density
is formed so that log Gamma((nu + 1) / 2) - log Gamma(nu / 2)
never cancels, which holds it to 7e-15 of a 300-bit reference for
every nu up to 1e50.
Robust regression
huber() fits Huber's least-favorable distribution: gaussian within
k residual standard deviations of mu and Laplace outside, so a
far-out point pulls on the fit with a bounded influence instead of
its squared distance. It is a proper normalized density, not a
penalty, so this is ordinary maximum likelihood and logLik(),
AIC() and the likelihood-ratio machinery all mean what they say.
k is a fixed constant of the family, huber(k = 1.345), not a
distributional parameter. It states where the analyst draws the line
between a residual and an outlier, which is a modelling choice rather
than something the data identifies; MASS::rlm() treats it the same
way. Estimating it would let the likelihood buy fit by widening the
gaussian core, which is the opposite of the point. As k grows the
family collapses to gaussian().
Point estimates track MASS::rlm(psi = psi.huber) closely but not
exactly, and the difference is the scale: rlm() fixes the scale at
a MAD-type estimate and iterates the location, while huber()
estimates sigma by maximum likelihood jointly with mu. The
coefficients agree to about 1e-2 on well-behaved data; the
sigma estimates need not.
The working-likelihood caveat that applies to asym_laplace()
applies here for the same reason. If the data are not actually
Huber-distributed - and the family is chosen precisely because the
error distribution is unknown - then the model is misspecified, the
information matrix is not the variance of the score, and Wald
standard errors and confint() intervals from it are not
calibrated. The point estimates stay consistent for the location.
Use frm_bootstrap() for intervals you can defend.
cens() and trunc() are unavailable: the CDF is a three-piece
function of a parameter-dependent residual and has no branch-free
form for the tape.
rho has a kink at |u| = k, so the objective is only piecewise
smooth and the optimizer often stops with a maximum absolute
gradient around 1e-4 and the accompanying false-convergence
warning. That is the kink, not a bad fit: the same thing happens to
asym_laplace(). The estimates satisfy Huber's own estimating
equations, X' psi(u) = 0 with psi(u) = min(max(u, -k), k), to
the same order. frm_allfit() confirms the fit when in doubt.
Conway-Maxwell-Poisson counts
In compois(), mu is the mean of the count and nu is the
dispersion: nu > 1 gives less variance than the Poisson, and
nu < 1 more. The density is RTMB::dcompois2(), so
exp(Intercept) under the log link is the expected count. This is
not brms's parameterization. brms 2.23.0's com_poisson() fits mu
as the rate lambda of the density lambda^y / (y!)^nu, which is
not the mean unless nu = 1. On counts simulated with lambda = 9
and nu = 2 (n = 400), frmtmb's exp(Intercept) is 2.78, the sample
mean, and brms's is 7.30. The coefficients of the two fits are not
comparable one by one (dev/remeasure-findings.md, R6).
Links
Every constructor takes link for the mean and a link_<dpar> for
each of its other distributional parameters, following brms:
student(link_sigma = "softplus"), zero_inflated_poisson(link_zi = "identity"). frmtmb-links lists the whole roster, what each link
maps, which families take it for the mean, and which set each
parameter admits.
The four families 'stats' owns, gaussian(), poisson(),
binomial() and Gamma(), have no frmtmb constructor to carry
these. Reach their links through brmsfamily():
brmsfamily("gaussian", link_sigma = "softplus").
The link for the mean may be unquoted, student(identity), and is
refused by name when brms does not allow it for the family:
bernoulli("sqrt") is an error.
An ordinal family's link is not a link on a mean. It names the
distribution function the thresholds are read through, so
all five ordinal families take logit, probit, probit_approx,
cloglog and cauchit (and cumulative(), hurdle_cumulative()
and acat() also take softit), and refuse anything else. Off the
logit, acat() reads its categories in brms's second form, a product
of distribution functions times a reversed product of survivals (see
frmtmb-links).
cumulative() and hurdle_cumulative() keep their thresholds
increasing, because their category probabilities are differences of
the distribution function. A cs() term moves each threshold of a
row by its own amount, so it can make one row's thresholds cross;
brms fits that with a warning, and so does this package: such a
row's density and fitted() probabilities are NaN, and
simulate() gives NA there. The
thresholds of sratio(), cratio() and acat() are unconstrained,
as in brms: their category probabilities are positive for any
thresholds, so a fit may have two of them cross.
Examples
set.seed(4)
n <- 120
dd <- data.frame(x = rnorm(n))
# heavier tails than gaussian(), with an estimated df
dd$y <- 1 + 0.8 * dd$x + rt(n, df = 4)
fixef(frm(bf(y ~ x) + student(), data = dd))
#> Estimate Est.Error Q2.5 Q97.5
#> Intercept 0.9572316 0.0937124 0.7735586 1.140905
#> x 0.8342758 0.1046305 0.6292037 1.039348
# counts with more spread than poisson() allows
dd$cnt <- rnbinom(n, mu = exp(0.5 + 0.4 * dd$x), size = 2)
fit <- frm(bf(cnt ~ x) + negbinomial(), data = dd)
fixef_by_dpar(fit)$mu
#> (Intercept) x
#> 0.3690287 0.3371063
# a zero-inflated count: the zi dpar gets its own predictor
dd$zi <- ifelse(runif(n) < 0.3, 0, dd$cnt)
frm(bf(zi ~ x, zi ~ 1) + zero_inflated_poisson(), data = dd)
#> Family: zero_inflated_poisson
#> Links: mu = log; zi = logit
#>
#> Formula: zi ~ x
#> Data: dd (Number of observations: 120)
#> Method: ML logLik: -173.275 AIC: 352.55 BIC: 360.913
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept 0.61 0.12 0.38 0.84 5.20 2e-07
#> zi_Intercept -0.32 0.24 -0.80 0.16 -1.31 0.190
#> x 0.23 0.11 0.02 0.44 2.13 0.033
# an ordered response: level order is the category order
dd$grade <- cut(1 + 0.8 * dd$x + rlogis(n), 3,
labels = c("low", "mid", "high"), ordered_result = TRUE)
frm(bf(grade ~ x) + cumulative(), data = dd)
#> Family: cumulative
#> Links: cdf = logit
#>
#> Formula: grade ~ x
#> Data: dd (Number of observations: 120)
#> Method: ML logLik: -81.675 AIC: 169.35 BIC: 177.713
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept[1] -0.88 0.22 -1.32 -0.45 -4.01 6.0e-05
#> Intercept[2] 3.68 0.51 2.69 4.68 7.26 3.9e-13
#> x 1.11 0.26 0.61 1.61 4.33 1.5e-05
# a proportion in (0, 1)
dd$p <- plogis(0.2 + 0.6 * dd$x + rnorm(n, 0, 0.3))
frm(bf(p ~ x) + Beta(), data = dd)
#> Family: beta
#> Links: mu = logit; phi = log
#>
#> Formula: p ~ x
#> Data: dd (Number of observations: 120)
#> Method: ML logLik: 154.805 AIC: -303.61 BIC: -295.247
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept 0.21 0.03 0.16 0.27 8.07 7.1e-16
#> x 0.58 0.03 0.52 0.64 19.06 < 2e-16
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI
#> phi 50.46 6.45 39.27 64.83
# bounded influence: a few wild points barely move the slope
dd$rob <- 1 + 0.8 * dd$x + rnorm(n)
dd$rob[1:5] <- dd$rob[1:5] + 30
fixef_by_dpar(frm(bf(rob ~ x), family = huber(), data = dd))$mu
#> (Intercept) x
#> 1.1663907 0.8551161
fixef_by_dpar(frm(bf(rob ~ x), family = gaussian(), data = dd))$mu
#> (Intercept) x
#> 2.231838 1.490393
# an unordered factor: one predictor per non-reference category,
# named after the level it belongs to
dd$pick <- factor(sample(c("ale", "stout", "lager"), n, TRUE))
cat_fit <- frm(bf(pick ~ x), family = categorical(), data = dd)
fixef(cat_fit) # mulager and mustout; ale is the
#> Estimate Est.Error Q2.5 Q97.5
#> mulager_Intercept 0.28202930 0.2503092 -0.2085678 0.77262639
#> mustout_Intercept 0.56454273 0.2363012 0.1014010 1.02768450
#> mulager_x -0.08958449 0.2669129 -0.6127242 0.43355523
#> mustout_x -0.43738957 0.2601125 -0.9472007 0.07242154
# reference
head(fitted(cat_fit)) # n x K category probabilities
#> , , P(Y = ale)
#>
#> Estimate Est.Error Q2.5 Q97.5
#> [1,] 0.2564172 0.04058006 0.1851532 0.3435459
#> [2,] 0.2163825 0.04641332 0.1390304 0.3207386
#> [3,] 0.2928210 0.05531228 0.1969847 0.4113962
#> [4,] 0.2768567 0.04647044 0.1954406 0.3763250
#> [5,] 0.3327697 0.08839367 0.1860379 0.5211359
#> [6,] 0.2819019 0.04890836 0.1964583 0.3866265
#>
#> , , P(Y = lager)
#>
#> Estimate Est.Error Q2.5 Q97.5
#> [1,] 0.3334248 0.04389351 0.2535408 0.4241756
#> [2,] 0.3011704 0.05134162 0.2108436 0.4100850
#> [3,] 0.3584389 0.05855308 0.2532796 0.4792380
#> [4,] 0.3479778 0.04992207 0.2574648 0.4509849
#> [5,] 0.3810589 0.08940322 0.2265305 0.5641186
#> [6,] 0.3513701 0.05233584 0.2567179 0.4593533
#>
#> , , P(Y = stout)
#>
#> Estimate Est.Error Q2.5 Q97.5
#> [1,] 0.4101580 0.04632187 0.3233136 0.5029901
#> [2,] 0.4824471 0.05570474 0.3757806 0.5907376
#> [3,] 0.3487401 0.05923469 0.2431006 0.4716782
#> [4,] 0.3751655 0.05205583 0.2798188 0.4812890
#> [5,] 0.2861714 0.07905345 0.1580859 0.4611873
#> [6,] 0.3667279 0.05416057 0.2682803 0.4777160
#>
# one category may take its own predictor
dd$w <- rnorm(n)
frm(bf(pick ~ x, mustout ~ w), family = categorical(), data = dd)
#> Family: categorical
#> Links: mulager = identity; mustout = identity
#>
#> Formula: pick ~ x
#> Data: dd (Number of observations: 120)
#> Method: ML logLik: -128.572 AIC: 265.145 BIC: 276.295
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> mulager_Intercept 0.24 0.24 -0.24 0.72 0.99 0.321
#> mustout_Intercept 0.51 0.23 0.06 0.97 2.22 0.027
#> mulager_x 0.19 0.21 -0.24 0.61 0.86 0.388
#> mustout_w 0.10 0.19 -0.28 0.48 0.52 0.600
# an angle: mu is the mean direction, kappa the concentration
dd$angle <- atan2(sin(0.5 + dd$x), cos(0.5 + dd$x))
vm_fit <- frm(bf(angle ~ x), family = von_mises(), data = dd)
head(fitted(vm_fit)) # the mean direction, in radians
#> Estimate Est.Error Q2.5 Q97.5
#> [1,] 0.82530547 0.01603449 0.79366340 0.85651524
#> [2,] -0.05745388 0.01824794 -0.09319702 -0.02167399
#> [3,] 1.41150621 0.01768376 1.37633078 1.44565802
#> [4,] 1.18172554 0.01719207 1.14764783 1.21504172
#> [5,] 1.83731815 0.01706178 1.80313237 1.87003914
#> [6,] 1.25871540 0.01740169 1.22418332 1.29240058