frechet(), gen_extreme_value() and discrete_weibull() are
brms's families of the same names, with brms's parameters, links and
defaults.
Usage
frechet(link = "log", link_nu = "logm1")
gen_extreme_value(link = "identity", link_sigma = "log", link_xi = "log1p")
discrete_weibull(link = "logit", link_shape = "log")Arguments
- link
Link for
mu. See frmtmb-links.- link_nu
Link for
nu:"logm1"(the default) or"identity".- link_sigma
Link for the GEV's
sigma:"log"(the default),"identity","softplus"or"squareplus".- link_xi
Link for
xi:"log1p"(the default) or"identity".- link_shape
Link for the discrete Weibull's
shape:"log"(the default),"identity","softplus"or"squareplus".
Frechet
frechet() is for a strictly positive response with a heavy right
tail. Its parameters are mu and nu. As in brms, mu is the MEAN,
not the scale: the density is Stan's frechet(nu, s) with scale
s = mu / gamma(1 - 1 / nu), which needs nu > 1. The default link
of nu is logm1, log(nu - 1), which keeps it there. Under
link_nu = "identity" a nu at or below 1 has no mean, and the
density is NaN there, as brms's is, in the fitted objective as well
as on doubles.
The variance is finite only for nu > 2; residuals(type = "pearson") divides by Inf where it is not.
Generalized extreme value
gen_extreme_value() has a location mu, a scale sigma and a
shape xi. The density is
exp(-t^(-1/xi)) t^(-1 - 1/xi) / sigma with
t = 1 + xi (y - mu) / sigma, and its limit at xi = 0 is the
Gumbel density, which the family evaluates exactly at xi = 0 and
to full precision near it. So bf(y ~ x, xi = 0) fits a Gumbel
model, as in brms.
The support depends on the parameters: y > mu - sigma / xi for
xi > 0 and y < mu - sigma / xi for xi < 0. Outside the
support the density is 0. Its log falls by 1e8 per unit of
t below 0 rather than jumping to -Inf, so a row a millionth
outside has a log density below -100 and the optimizer meets a slope
that turns it back, not a wall. The start, xi = 0, is inside the
support for every data set. For xi > -1 the density goes to 0 at
an end of the support, so a maximum has every row inside it. That is
the region brms's scale_xi() maps its tmp_xi onto.
In a mixture() a component that does not reach a row adds nothing
to that row's density, which is what a mixture density is. brms's
density is NaN there, which rejects the parameter value, so a brms
mixture holds each component's support over every row.
The default link of xi is log1p, log(1 + xi), brms's link
for a predicted xi, which keeps xi > -1. link_xi = "identity"
is the other choice brms allows. Under it, xi <= -1 has no
likelihood, so every link holds xi above -1. Below -1 the
likelihood is unbounded and has no maximum.
brms parameterizes an xi without a formula differently. It samples
an unbounded tmp_xi and scales it into the support the current mu
and sigma leave, and it reports xi from it after the fit. That is
a change of variable, so the maximum of the likelihood is the same
point. A fit does not use it, because its range reaches below
xi = -1. frmtmb.sample::frm_sample() does, because there the
support's edge is a wall the sampler meets; it holds the range above
-1 (see frmtmb-gev-xi-chart).
A fit whose xi stops within 1e-3 of -1 with the likelihood not
stationary in xi there warns: its log-likelihood is a supremum, not a
maximum. A maximum inside that band gets a message. For
-1 < xi < -0.5 the maximum exists but is not regular (Smith, 1985),
so its standard errors do not have their usual meaning; outside the
band no message says so.
The mean, which fitted() returns, is
mu + sigma (gamma(1 - xi) - 1) / xi, and mu + sigma * 0.5772...
(Euler's constant) at xi = 0. It is infinite for xi >= 1, and
fitted() returns Inf there. brms's posterior_epred() uses the
same formula and returns NaN at xi = 0 and xi = 1 and a finite,
wrong number for xi > 1. The variance is finite only for
xi < 1 / 2.
brms 2.23.0 marks gen_extreme_value deprecated and warns when it
builds the Stan code. frmtmb does not repeat the warning: it is about
brms's plans for its own code, not about the model.
Discrete Weibull
discrete_weibull() is for a count 0, 1, 2, .... Its parameters
are mu in (0, 1) and a positive shape, with
P(Y >= y) = mu^(y^shape), so P(Y = y) = mu^(y^shape) - mu^((y + 1)^shape). mu is not the mean: it is the probability of a
count of at least 1. The mean, which fitted() returns, is the sum
of mu^(y^shape) over y >= 1. frmtmb sums it to a relative
precision of about 1e-12, the tail beyond 1000 terms through its
integral. brms's posterior_epred() stops the sum at the first term
below 0.001, so the two differ where the tail is long.
cens() and trunc() follow brms's discrete conventions (see
frmtmb_family(), Censoring a discrete response): a right-censored
count at k is P(Y > k), and trunc(lb = ) reads its bound as
inclusive.
Addition terms
All three take weights(), cens() and trunc(). frechet and
gen_extreme_value also take mi(), as in brms.
References
Smith, R. L. (1985). Maximum likelihood estimation in a class of nonregular cases. Biometrika, 72(1), 67-90.
Examples
set.seed(1)
n <- 200
dd <- data.frame(x = rnorm(n))
# Frechet: mu is the mean
nu <- 4
mu <- exp(1 + 0.3 * dd$x)
dd$yf <- mu / gamma(1 - 1 / nu) * rexp(n)^(-1 / nu)
fit_f <- frm(yf ~ x, family = frechet(), data = dd)
fixef(fit_f)
#> Estimate Est.Error Q2.5 Q97.5
#> Intercept 1.0225652 0.02844083 0.9668222 1.0783082
#> x 0.2980685 0.02175693 0.2554257 0.3407113
# generalized extreme value
xi <- 0.2
dd$yg <- 1 + 0.5 * dd$x + 0.8 * (rexp(n)^(-xi) - 1) / xi
fit_g <- frm(bf(yg ~ x, sigma ~ 1), family = gen_extreme_value(),
data = dd)
fixef(fit_g)
#> Estimate Est.Error Q2.5 Q97.5
#> Intercept 1.0019245 0.06764939 0.8693341 1.13451483
#> sigma_Intercept -0.1671922 0.06467891 -0.2939605 -0.04042382
#> x 0.4751598 0.06185628 0.3539237 0.59639586
# discrete Weibull counts
q <- plogis(1 + 0.4 * dd$x)
dd$yd <- ceiling((log1p(-runif(n)) / log(q))^(1 / 1.3) - 1)
fit_d <- frm(yd ~ x, family = discrete_weibull(), data = dd)
head(fitted(fit_d))
#> Estimate Est.Error Q2.5 Q97.5
#> [1,] 1.392365 0.1285444 1.161902 1.668541
#> [2,] 1.748125 0.1257542 1.518238 2.012820
#> [3,] 1.312478 0.1357565 1.071637 1.607446
#> [4,] 2.586694 0.3434126 1.994062 3.355457
#> [5,] 1.820845 0.1338922 1.576455 2.103122
#> [6,] 1.318117 0.1352124 1.078047 1.611648