logLik() gives the log-likelihood of the data at the estimates.
logLik_penalized() gives the value that the optimizer maximized: the
log-likelihood plus the log prior term. For a fit with no prior, or
with a prior that only sets bounds, the two are equal.
Arguments
- object
A
frmtmb_fit.- ...
For
AIC()andBIC(), more fitted models, as instats::AIC(). Refused otherwise.- k
The penalty per parameter, as in
stats::AIC().
Value
logLik(): an object of class "logLik" with attributes
df, nobs and REML, and df_why when df is NA.
logLik_penalized(): one number, with attribute log_prior, the
log prior term at the estimates (see "The log prior term"). That is
0 with no prior, and NA when it does not separate from the
integral (see "Random effects and REML"). AIC() and BIC(): as in
stats::AIC().
Details
A prior (frm(prior = )) adds a log prior term to the objective, so
the fit is a maximum a posteriori (MAP) fit. logLik() subtracts that
term at the estimates from the maximized value. The result is the
log-likelihood at the MAP estimates: the data part of the log density,
which brms's log_lik() gives row by row. glmmTMB's logLik() with
priors = gives the penalized value, which is logLik_penalized()
here.
The log prior term
The term, attribute log_prior of logLik_penalized(), is exactly
what the optimizer adds. It is the sum of each prior's log density on
the scale that the optimizer moves, so it includes the log-Jacobian
of every map from that scale to the parameter the prior is about: the
log of a standard deviation (sd, sigma and other positive
parameters), a dpar's link, the chart of an autocorrelation
coefficient, the increments of ordinal thresholds, the chart of a
simplex (simo and other Dirichlet priors) and the Cholesky chart of
an lkj() prior. A prior with a bound has no truncation normalizer:
the bound is a box on the optimizer, not a change in the density.
So the term is not brms's lprior. brms evaluates each density on the
constrained scale with no Jacobian, which Stan adds to lp__
separately, and divides a bounded prior by the probability inside its
bounds. For normal(0, 1) on an sd at the value s, the term here
is dnorm(s, 0, 1, log = TRUE) + log(s); brms's lprior is
dnorm(s, 0, 1, log = TRUE) + log(2). Only the term here is what
separates logLik() from logLik_penalized().
Random effects and REML
With random effects, the objective is an integral over them: the
Laplace approximation, quadrature = TRUE or importance =. A prior
on a parameter outside the integral (a fixed effect, a standard
deviation, a correlation) is a constant in that integral, so the
split is exact. Under REML = TRUE the fixed effects are inside the
integral, so a prior on them (class b, Intercept, or an nlpar)
does not separate from it. Then logLik() stops with an error, and
logLik_penalized() still gives the maximized value.
Information criteria
The MAP estimates do not maximize the log-likelihood, and a parameter
that a prior shrinks is not a full degree of freedom. Thus, under a
prior with a density, logLik() has df = NA, with the reason in
attribute df_why. AIC(), BIC(), extractAIC() and so step()
stop with an error when they dispatch on such a fit, as
anova() and drop1() do. A call that reaches
stats's own methods with the df of logLik(), such as
AIC(lm_fit, map_fit), BIC(lm_fit, map_fit) or
AIC(logLik(map_fit)), gets NA for that fit. A package that counts
the parameters itself is not covered: insight::get_loglikelihood()
replaces the NA with its own count, so performance can still give
an AIC for such a fit. Refit without the prior to compare models by
AIC.
A prior that only sets bounds (set_prior("", class = "sd", lb = 0.1)) adds no density, and these methods work. So does a uniform()
prior on a parameter that the optimizer moves on its own scale (a
coefficient of class b or Intercept, or an nlpar): its log
density is a constant inside its bounds, so the estimates maximize the
likelihood there, as under the bounds alone. A uniform() prior on a
standard deviation, a dpar on its natural scale or a correlation is
read through a log-Jacobian, which is not constant, so it is a
density like any other.
See also
anova.frmtmb_fit() for likelihood-ratio tests between
maximum-likelihood fits.
Examples
set.seed(1)
dd <- data.frame(x = rnorm(60), g = factor(rep(1:12, 5)))
dd$y <- 1 + 0.5 * dd$x + rnorm(12)[dd$g] + rnorm(60)
fit <- frm(y ~ x + (1 | g), data = dd,
prior = set_prior("normal(0, 1)", class = "b"))
# the log-likelihood at the MAP estimates
logLik(fit)
#> 'log Lik.' -85.82969 (df=NA)
# the value the fit maximized, and the log prior it holds
logLik_penalized(fit)
#> [1] -86.88418
#> attr(,"log_prior")
#> [1] -1.054497
# an information criterion needs the maximum-likelihood fit
try(AIC(fit))
#> Error : AIC() counts each parameter as one degree of freedom at the maximum of the log-likelihood, and this fit has a prior that adds its log density to the objective (frm(prior = )). Its estimates maximize the log-likelihood plus the log prior, not the log-likelihood, and a parameter that the prior shrinks is not a full degree of freedom. Refit without the prior to compare models by AIC(); a prior that only sets bounds can stay. To compare the models under the prior, use loo() or bayes_factor() of frmtmb.sample
AIC(frm(y ~ x + (1 | g), data = dd))
#> [1] 179.6546