Fit a model — frm • frmtmb

Fits a model specified with bf() by maximum likelihood, using the Laplace approximation for random effects through RTMB.

Usage

frm(
  formula,
  data,
  family = NULL,
  REML = FALSE,
  start = NULL,
  control = frmtmb_control(),
  se = FALSE,
  na.action = stats::na.omit,
  prior = NULL,
  quadrature = FALSE,
  importance = 0L,
  data2 = list(),
  dry_run = NULL,
  verbose = FALSE,
  drop_unused_levels = TRUE,
  threshold = NULL
)

Arguments

formula

A frmtmb_formula from bf() (with a family attached via +), or a plain formula combined with the family argument. With neither, the family is gaussian(). A . stands for every column of data that is not on the left of the formula, as in brms: y ~ . is expanded with stats::terms() before the formula is read, and each parameter formula is expanded on its own.

data

A data frame. A tibble, a data.table, or a plain named list of equal-length columns is accepted as well, since each reaches stats::model.frame() unchanged. A matrix column is a supported model variable and enters the design as its own block of columns (this is how a functional predictor or a cbind() term is written). A list column is not a model variable and is refused by model.frame() if the formula names one, though it may sit unused in data.

family

A family: a frmtmb_family, a stats::family object, a family constructor with or without its parentheses (for example gaussian(), poisson, cumulative), or a family name as a string. It overrides a family already attached to formula; in a multivariate model it fills only the responses that have none. A multivariate model also takes a list of families, as in brms, one per response in formula order: family = list(gaussian(), poisson()). Each entry fills its response when that response has no family. The default, NULL, means gaussian() - the brms, lme4 and glmmTMB convention - so frm(y ~ x, data = d) is a linear model.

REML

If TRUE, integrate the mu fixed effects out of the likelihood along with the random effects (restricted maximum likelihood). For a gaussian model the Laplace approximation is exact for that integral, so this IS classical Patterson-Thompson REML. For any other family no error-contrast derivation exists; what is computed is the Laplace-approximated integrated likelihood for the variance parameters under a flat prior on the fixed effects, which agrees to first order with the Cox-Reid adjusted profile likelihood (the general-model REML analogue). glmmTMB computes the same quantity. It usually reduces the finite-sample downward bias of variance components, which is REML's purpose, but it is an approximation stacked on an approximation, not an exactness result. The usual restriction carries over: anova() compares REML fits only when their fixed-effect designs agree.

Distributional-parameter coefficients (a sigma ~ z model's, for example) are deliberately NOT integrated: they belong to the variance-parameter set, exactly as a varExp() coefficient does in nlme::gls(method = "REML"), and the two agree to optimizer precision for gaussian models. Random effects appearing in a distributional parameter's own formula are integrated like any other latent. mgcv's method = "REML" for a location-scale family such as gaulss() integrates the coefficients of every linear predictor, so for a smooth in sigma the two REML criteria are different and give different curves and smoothing parameters.

start

Optional named list of starting values, one entry per parameter component (beta, betad, theta, ...). par_template() returns that list already filled with the defaults, for a fitted model or for a formula and data alone; edit it and pass it back. Its section "The default starting values" tells how the defaults are set; the entries you do not name in start keep those defaults.

An exact gp() on a gaussian response with a constant sigma starts at the maximum of the model's marginal, read from the best cell of six length scales on a coarse grid. Beside ar() in brms's default form (cov = FALSE, any order, a free sigma, no se() or weights()) that marginal holds the ar coefficients too, so they start at its maximum. On step-shaped data the ar coefficient there is often at or just above 1, as brms's unconstrained cov = FALSE likelihood allows. Another family, a non-identity link, sigma ~, ma(), arma(), ar(cov = TRUE), mi(), me() or mo() starts at that marginal of the response on its link's scale, for sdgp and the length scale. Where that marginal has a second maximum at another length scale, the fit also starts there and keeps the better of the two. A field in a parameter of a nonlinear body starts at the marginal of the body made linear at its start, and the fit also descends from the best of a few length scales and keeps the better; a start of your own turns the first off. A mixture and sigma ~ gp() start at the best of a few length scales read from the objective. A theta here replaces all of this, and a beta or betad replaces only itself.

None of this can be certain. A marginal can have more than one maximum. On 744 rows of curves no test had seen, the fit stopped more than 0.001 short of the maximum on 1, the worst 0.56, with no warning, where 0.78.0 read 40 and 62 short on the two tables of 192. A fit from other start values is how to check one that matters.

A component is read by NAME when its vector has names, so start = list(beta = c("(Intercept)" = 5000)) moves that one coefficient and leaves the rest at their defaults; parentheses may be dropped, as everywhere a parameter is named. A vector with no names is positional and must be full length. Naming some entries of one vector and not others is refused.

A nonlinear body can have more than one maximum, and nothing in a fit marks a second one. From the default start, frm() also fits from each sign pattern of the body's coefficients and returns the highest point (see "The default starting values" in par_template()). A start, or options(frmtmb.nl_sign_probe = FALSE), skips that; a located prior places the start and does not skip it. The point a pattern reaches can itself be a second maximum: the fit then converges with no warning and is still below the highest one, as on 1 of 89 such fits measured (1.84 short, where the default start alone stopped 2.29 short with a warning). A start elsewhere, or a located prior, can show a higher maximum; frm_allfit() cannot, as it starts every optimizer from one point.

control

A list from frmtmb_control().

se

If TRUE, run RTMB::sdreport() at fit time. The default (FALSE) defers it until standard errors are first needed (summary, vcov, confint, frm_linpred(se.fit = TRUE)), which cuts roughly a quarter off fit time in fit-and-predict or bootstrap loops. The deferred report is cached, so nothing is computed twice. Either way the fit builds the Hessian at the optimum and warns when a parameter has no standard error; see check_se in frmtmb_control().

na.action

How to handle missing values, as in stats::lm() (default stats::na.omit). Rows dropped for missingness are reported in a message; wrap the call in suppressMessages() to silence it.

prior

Optional set_prior() specification. This makes the fit MAP / regularized ML (glmmTMB's priors= in spirit): useful for stabilizing singular variance components or separating binomials. logLik() is then the log-likelihood at the MAP estimates, and logLik_penalized() the maximized value with the log prior. AIC(), BIC(), anova() and drop1() refuse a prior with a density, as its estimates do not maximize the likelihood. A brmsprior object built by brms's own prior() is translated row by row, so priors copied out of a brms script work whichever package's prior() was in scope. The argument takes brms's spelling, prior; the priors of releases before 0.43 is gone rather than aliased, and a call still using it fails as an unused argument. A named list of prior objects keyed by internal parameter names (list(theta_1 = prior_normal())) is refused, as brms has no such spelling: it gave a parameter a second prior slot beside its class. A plain data frame is refused too, as brms refuses one; as.brmsprior() reads it. set_prior("", class = ...), brms's (flat) on a slot, applies nothing here.

This is also where HARD BOUNDS are written, as lb/ub on a specification that may carry no distribution at all: set_prior("", nlpar = "guess", lb = 0, ub = 1) is a box constraint and nothing else. The lower/upper arguments of releases before 0.49 are gone rather than aliased, and a call still using them fails as an unused argument; every outer parameter they could reach has a set_prior() class: a standard deviation class "sd", on the sd scale, and an internal covariance parameter that no brms class names class "theta" (coef = "theta_12"). Two kinds take no box: a correlation, random-effect, rescor or unstr(), which its class reaches only as a whole matrix (as in brms, it takes no lb/ub), and the internal parameters of a cov = TRUE ar() or ma() of order two or more (brms has no such term). See set_prior()'s Hard bounds section.

The spelling is brms's; the SEMANTICS are not. A prior here is a penalty on the likelihood and the answer is one mode. It is not a posterior, and no interval this fit reports is a credible interval. The two land close where the data dominates: a kidney frailty model given brms's own priors returns sd(patient) 0.38 against brms's posterior mean of 0.40. That is a measurement on one model, though, not a property of the translation. Use frmtmb.sample::frm_sample() when the posterior is the answer.

A prior with a location places a nonlinear start. brms uses its priors to place the sampler, and frm() now reads them the same way, for nonlinear parameters only: where start does not set it, a nonlinear coefficient begins at the location of a normal(), student_t() or cauchy() prior on it, and a message names what was placed. This is what lets the insurance-loss growth curve fit from its brms priors with no start at all. Every other parameter keeps the start it always had.

Without such a prior a nonlinear model still needs start. frm() optimizes rather than samples, so it evaluates the objective AT the starting values, and a nonlinear body is rarely defined at zero. par_template() shows the names and the values in force.

No start identifies a nonlinear model in which a combination of the nonlinear parameters' coefficients changes no fitted value, as when two of them enter only through their sum or difference and their formulas share a term (a + b with a ~ 1 + x and b ~ 1 + x). frm() fits it and then warns, naming the coefficients, when the likelihood at the optimum is flat along such a combination (its curvature below 1e-9 of the largest, on the Hessian of those coefficients scaled to unit diagonal). A curved ridge that does not flatten that far at the optimum (exp(a)^k) is not flagged there; vcov() and summary() then report its standard errors as not finite. A prior that reaches every flat direction identifies the model, as brms's priors do; the split along that direction is then set by the prior alone, so a weak prior gives very wide intervals.

quadrature

If TRUE, marginalize each scalar random effect by adaptive Gauss-Kronrod quadrature instead of the Laplace approximation (the glmer(nAGQ = k) analogue; matches it in tests). Worth it for Bernoulli responses with small clusters, where Laplace biases variance components. Scalar random effects only (one coefficient per level of a block), and not with mi(), trunc(), REML = TRUE, or frmtmb_control(profile = TRUE). A row may read two blocks when one is nested in the other, as in (1 | ga/gb). Such a model is integrated as an iterated integral by Gauss-Hermite rules, frmtmb_control(quad_nodes = ) nodes per level, centered and scaled by the Laplace fit and re-anchored at the estimate. That route needs exactly two blocks on the same predictors, one response, Gaussian blocks and per-row densities, and its memory is limited by frmtmb_control(quad_max_gb = ). Its accuracy depends on the design: against independent references the log-likelihood error at the default nodes was 9e-11 to 7e-6 on seven designs with standard deviations near 1, and 5e-4 to 6e-4 with an inner standard deviation near 3 or an outer one of 3; on a near-separable response it was 1.7e-2, where the rule did not settle and the fit warned. More nodes do not always help on the hard designs. The Laplace error on the same designs was 0.07 to 1.4. A fit whose rule does not settle warns (class frmtmb_quad_unsettled) and gives the rule's own error at the estimate, the log-likelihood at the fit's nodes against 21 x 21 (or 10 more per level); on such a design more nodes are not a sure remedy, since the error need not fall as they rise. Profile and root-search intervals and ranef(condVar = TRUE) are refused on a nested fit (class frmtmb_quad_nested_unsupported); Wald and bootstrap intervals work. Crossed blocks, (1 | g) + (1 | h), are refused, because their integral does not factor into one-dimensional ones; so are (1 | g) in both mu and sigma, where the rule converges slowly, and three or more blocks on a row. The blocks of a by-split gr(g, by = f) each read rows of their own and are allowed. A plain Laplace fit runs first and the quadrature tape is built at its optimum: the Gauss-Kronrod rescaling is fixed when the tape is built, so the starting point decides whether the marginalized objective is finite at all. That fit also supplies the conditional modes, which the marginalized objective no longer carries, so ranef(), fitted() and predict() work as usual.

importance

Number of importance draws per group, or 0 (default) for none. A positive count replaces the Laplace approximation with an importance-sampling correction of it (Skaug and Fournier 2006, the ADMB -is option): each group's marginal likelihood is re-estimated by reweighting draws from the Laplace Gaussian, which is unbiased for the exact integral rather than accurate to O(n^-1). Like quadrature = TRUE it is worth it for binary responses in small clusters, and unlike it, the blocks may have any dimension: this is the correction for a correlated random slope, where quadrature cannot go.

The cost is one tape of N stacked copies of the likelihood, so time and memory grow linearly in the draw count. An odd count is rounded up to the next even one, because the draws are taken in antithetic pairs (measured on this package's own probe design, the pairing is worth two to three times the draws it costs).

Scope, with everything else refused by name: any number of random-effect blocks, of any dimension, provided they share ONE grouping factor and one set of levels, in any covariance structure that is Gaussian within a level and independent between levels; one response; a family with a rowwise density. Not with quadrature, REML = TRUE, frmtmb_control(profile = TRUE), mi(), a residual correlation term, rescor, a nonlinear predictor, or cs(). trunc() and cens() are supported: the draws sit near the conditional mode, which is why the truncation normalizer that underflows at the Gauss-Kronrod nodes does not underflow here.

Distributional regression writes several blocks by construction: (1 | g) in mu and (1 | g) in sigma are two blocks over g, and a level's coefficients from both are drawn together from one joint proposal. Blocks over DIFFERENT grouping factors, crossed or nested, are refused.

The fit records the draw count, the seed, the rounds taken, the per-group effective sample sizes and the Monte Carlo standard error of the corrected log-likelihood in fit$importance; summary() prints them. logLik(), AIC(), vcov() and confint() all come from the corrected objective, with two caveats worth knowing. The reported estimate is the optimum of the previous round's proposal, so the covariance is a Hessian at a point that is stationary only to fit$importance$grad. And confint(method = "profile") walks one tape instead of refitting: its proposal follows the parameters to first order around the estimate, so it warns when a bound lands where the weights no longer cover the integrand.

The Monte Carlo standard error can read high or low. On a binomial model at 200 draws it reads 0.0044 where the error's standard deviation is 0.0017. On a Student response with nu = 5, 15 groups of 6 and 1000 draws, the error's standard deviation is 1.3 to 1.4 times the number, a nominal 95 percent interval covers 0.86, and no warning fires. To judge two corrected fits, refit with another importance_seed or with twice the draws and read how far the log-likelihood moves.

A count that is too small for the model stops the fit. The message gives the worst group's effective sample size and either a larger count at which the same fit settled or the counts it tried. That search can take several times the fit: up to four doublings, to at most 8192 draws and 2e6 rows times draws, measured at 15 times the fit and 3.3 GB at its worst. The seed, the round cap and the effective-sample-size threshold are frmtmb_control() settings.

data2

A named list of objects that are not columns of data: the adjacency matrix of car() and the mesh triple of spde() (both frmtmb.spatial), and the matrices of gr(prec = ), gr(cov = ) and equalto(). This is brms's data2 argument, with one deliberate extension: brms accepts a bare name only, while frmtmb also evaluates compound expressions with data2 in front of the data mask, so gr(g, cov = solve(Q)) finds Q there. Each structural expression resolves from data2 first, then data, then the formula environment; that last step is what a model written before data2 relied on, so old code keeps working. Prefer data2: its objects are stored on the fit, so saveRDS() and a later refit(), influence() or update() in a fresh session do not need the calling environment to still exist.

dry_run

"spec" returns the parsed intermediate representation without touching data; "frame" returns the assembled design matrices and parameter template without fitting; "objective" additionally tapes the objective and returns an UNFITTED object carrying it, which is what frmtmb.sample::frm_sample() samples when it is given a formula rather than a fit. Methods that report a maximum-likelihood quantity refuse on that object.

verbose

Report fit progress; a shortcut for control = frmtmb_control(verbose =), whose value wins when both are given. See frmtmb_control() for the levels and the output.

drop_unused_levels

Drop the factor levels that no row of the data takes, as brms's argument of the same name does. TRUE, the default, is brms's default too. With FALSE an unused level of a predictor gives a design column of zeros, which frmtmb drops with its rank-deficiency message because no data can estimate it (brms keeps the column, and only its prior places the coefficient). A grouping factor loses its unused levels either way, because a level without rows adds nothing to the Laplace likelihood, so ranef() lists the levels the data hold.

threshold

Refused. brms 1.x took the threshold structure of an ordinal model as an argument of brm(), and the JSS overview vignette still passes it there; brms 2.23.0 refuses it too ("passing unknown arguments: threshold"). It is an argument of the ordinal family, as in brms: family = sratio(threshold = "equidistant"). The refusal names that spelling with the family of the call.

Value

An object of class frmtmb_fit. It is a list, and two of its elements are read directly often enough to name here. fit$data2 is the data2 list of known matrices. fit$data is the MODEL FRAME, the object model.frame() returns: one column per term the formula names, under the term's own spelling, so a model with offset(log(Age)) carries a literal offset(log(Age)) column, and a data column no term uses is absent. The variables of an offset() and of an addition term's expression, as Age there, are columns too, because the grids of conditional_effects() and emmeans() hold them. brms keeps its VALIDATED RAW DATA under the same name, so names(fit$data) and ncol(fit$data) differ from brms on any model with a transformed term, and newdata = fit$data is not the brms idiom it looks like. The element exists because brms-shaped code reads fit$data and used to reach the data2 list through $'s partial matching; the model frame is the object this package already had. It is the same object in memory that fit$frame[["data_frame"]] holds, but R's serializer does not deduplicate a shared value, so a SAVED fit carries the frame twice: 5.3% more raw bytes and 14.0% more gzipped on a 20,000-row fit, 0.16% and 1.49% on a 240-row one (dev/adefects-log/p1-dollar.txt). Everything else the fit carries has an accessor, and the accessor is the supported route.

What a nonlinear body sees

The body of a nonlinear formula - the response formula under nl = TRUE, or an nlf() body - is evaluated as a small language for the AD tape, not as ordinary R. A function it CALLS is looked up in RTMB first, so a bare pnorm(), qgamma(), plogis(), besselK() or matrix() in a body is RTMB's tape-capable version. The RTMB:: prefix is no longer needed, and writing it changes nothing: both spellings tape to the same function and give the same fit.

RTMB WINS. It wins over the search path and over a function of the same name that you defined yourself, because the body is a language with its own vocabulary rather than a piece of your session. To reach a different pnorm inside a body, qualify it: stats::pnorm() is the escape hatch, and so is myPkg::pnorm().

Nothing else changes. Only the names RTMB replaces are shadowed, and only where the body calls them, so a helper function of your own, an object the body reads, and every name that is not in RTMB still resolve in the formula's environment exactly as before. A name the body READS is never shadowed either, so a data frame called df or a vector called order keeps its meaning.

The rule stops at the end-user surface. A frmtmb_family() log-density and a frmtmb_structure() log-likelihood are ordinary R functions written by an extension author, and they keep resolving lexically: qualify with RTMB:: there, and see frmtmb_ad_overload() for the rest of that contract.

Uncorrelated group-level terms, ||

(x || g) is brms's (x | gr(g, cor = FALSE)). It is one term, coded with the formula's own contrasts. Each column has its own standard deviation, and there are no correlations. So (1 + f || g), with a factor f of levels a, b and c, has the columns Intercept, fb and fc of (1 + f | g). Its standard deviations are sd_g__Intercept, sd_g__fb and sd_g__fc, and it has no cor_ parameter. (0 + f || g) has fa, fb and fc. lme4 expands || term by term instead. It codes the factor apart from the intercept, with a column for every level, and then correlates them (lme4#818). frmtmb follows brms.

(x | gr(g, cor = FALSE)) and diag(x | g) are the same term as (x || g). (x | mm(g1, g2, cor = FALSE)) is the || term over the membership grouping, (x || mm(g1, g2)). Under ||, a cor argument is ignored, as brms ignores it. Over gr(g, dist = "student"), a || term has one latent scale per level, shared by all its columns, as in brms.

Two combinations are refused. The first is an |ID| key beside ||, as in (x | p || g) or (x | p | gr(g, cor = FALSE)). It asks for correlation across terms and for none at once. brms refuses the first spelling and silently drops the key from the second. The second is a || term over gr(g, cov = A) or gr(g, prec = Q) that has a term with several columns in it (a factor, poly()). A || term over a known matrix is split into one block per term, which is brms's model for numeric columns. The known-covariance block correlates the columns of one term within a level, so a term with several columns is refused.

Monotonic effects

mo(x) fits an ordinal predictor without assuming its categories are equally spaced. x is an ordered factor or an integer vector. As in brms, the categories are coded 0..D from the lowest one: the first level of a factor, or the minimum of the integers. So 1..5 and -2..2 are both coded 0..4. Two categories (D = 1) give a simplex of one weight, which is fixed at 1, so the term is its coefficient times the code, as brms's Jmo = 1 is. One category (D = 0) gives a term that is 0 on every row, as brms's Jmo = 0 is. The data do not determine its coefficient, so in a maximum-likelihood fit its standard error is NaN, and it is not counted in the degrees of freedom of logLik(), as lm() does not count a column it aliases. A group-level term on it, (mo(x) | g), is refused. brms cannot sample such a model, and neither can frmtmb.sample. An integer between the minimum and the maximum that no row has is also a category. For example, 1, 2, 5 gives D = 4. The data do not determine a prediction at such a category, so in a maximum-likelihood fit its standard error is NaN. New data are coded from the minimum of the fitted data. A value that is not a whole number, or that is outside the fitted range, is refused, as brms refuses it. conditional_effects() holds a numeric mo() variable at round(median()), as brms does. The term contributes b * D * sum(zeta[1:k]) at category k, where zeta is a simplex of D non-negative steps summing to one. b is therefore the AVERAGE step, on the coefficient scale of any other predictor, and zeta says how that total is distributed over the categories. b appears in fixef() under the term's label (mox, or mox:z for an interaction); the simplex is held as D - 1 free coordinates in a zeta<j> component of par_template(). summary() reports the simplex itself under "Monotonic Simplex Parameters", as brms does, with delta-method standard errors and Wald intervals held inside [0, 1]; confint(), vcov() and start read the coordinates, which have no reading of their own.

The coordinates are those of stereographic projection onto the unit sphere, whose squared components are the simplex. zeta = 0 is the uniform simplex, and every face of the simplex, where a step is 0, is at a finite coordinate. A softmax puts a face at an infinite coordinate, where the likelihood is flat toward it and an optimizer stops short; this parameterization does not, so a step that is 0 at the maximum is 0 at the estimate and keeps a finite curvature. After the fit, a term whose simplex is on a face is also refitted with its scale coefficient held at the other sign, because each sign of b has its own local maximum, and the better of the two is kept (frmtmb_control(mo_search = FALSE) turns that off). The sampler of frmtmb.sample does not use these coordinates: it samples the softmax under brms's default dirichlet(1) and reports brms's simo_ weights.

Every monotonic TERM gets its own simplex. y ~ mo(x) * z fits two: one shape for the main effect and one for the interaction, because the two terms describe different shapes and there is no reason for the interaction to bend the way the main effect does. The simplexes are numbered in the order brms::brm() enumerates its special terms, which is stats::terms() order and so lists every main effect before any interaction. The j-th monotonic term of a linear predictor is therefore the one brms calls simo_<j> and names <label><1> in the simo rows of its get_prior().

The zeta<j> NUMBER is not that j in general. Simplexes continue the numbering of whatever parameters the family contributed first, so zeta<j> is brms's simo_<j> only for a family that declares none: an ordinal fit spends slot 1 on its thresholds and a cox() fit (frmtmb.survival) on its baseline hazard, and in both the first monotonic simplex is zeta2. Read the terms in order rather than parsing the number.

With no prior, a fit applies no density to a simplex; brms's flat dirichlet(1) contributes only a constant, so the two log-densities agree up to lgamma(D) per simplex. A simo prior row, set_prior("dirichlet(c(2, 1, 1))", class = "simo", coef = "mox1"), makes the fit a MAP fit under that Dirichlet: brms's spelling, with the coef get_prior() lists for each simplex. Like every prior of a MAP fit it is the density of the coordinates, so it carries the chart's Jacobian, prod |u_k| times the projection's conformal factor. In plain terms: the MAP fit under dirichlet(a) has the mode of dirichlet(a + 1/2) on the weights, pulled further toward the uniform simplex by the conformal factor. Thus dirichlet(1) is not the absence of a prior: it moves the mode at least as much as dirichlet(1.5) on the weights does. Only a fit with no simo row applies no density to the simplex. The pull is large where the data say little: with 60 rows, the weights of a dirichlet(1) MAP fit were 0.30, 0.38 and 0.32 where the maximum likelihood weights were 0.05, 0.73 and 0.22; with 400 rows they moved by 0.003 at most. To pull less, give no simo row. A concentration below 1/2 has no maximum in these coordinates and is refused; frmtmb.sample's frm_sample() samples any.

An interaction codes the columns brms's special-term model matrix codes, each with its own simplex and coefficient: mo(x) * f with a factor f of levels a, b, c fits mox, mox:fb and mox:fc, and mo(x):f fits mox:fa, mox:fb and mox:fc. The other factors of a label may be several, mo(x) * z * w, and are named as brms names them, in the order model.matrix() reads the special terms alone: mo(x) * z and z * mo(x) give mox:z, and z + mo(x):z gives z:mox. mo(x):mo(w) is refused. brms 2.23.0 builds the same columns but multiplies every one by the first covariate in its Stan code, and stops on a factor of three or more levels crossed with the main term (brms#1828); frmtmb fits the model brms 2.22.0 fits and brms's Stan data describe.

A group-level term may vary a monotonic coefficient, as in brms: y ~ mo(x) + (mo(x) | g) adds r_g[, mox] to bsp_mox and multiplies the sum by the same transform, (b + r_g) * D * sum(zeta[1:k]), so the group-level term has no simplex of its own. The population-level column it varies is required, and refused without it, where brms 2.23.0 drops the term silently. It is available in a plain (... | g) or (... || g) term. A || term with mo() is one block with no correlations, coded whole like any || term: (mo(x) * f || g) has the columns Intercept, fb, fc, mox, mox:fb and mox:fc.

Category-specific effects, cs()

cs(x) gives a predictor one coefficient per category boundary instead of one for the whole response, in every ordinal family. The term contributes an n by K - 1 matrix of offsets to the thresholds, so the effect of x on the first boundary need not be its effect on the second. On cumulative() and hurdle_cumulative() the offsets can make one row's thresholds cross, and the category between two crossed thresholds then has a negative probability. brms fits these models with a warning that the effects are experimental, and so does this package. Such a row's density is NaN, as brms's is; its fitted() probabilities are NaN too, where brms returns the negative difference; and simulate() gives NA for it. predict() reports NA proportions for a row whose thresholds cross at the estimates, where fitted() is NaN. A row whose thresholds cross only in some of its simulated parameter draws reports the proportions over the draws that do not cross, and the call warns once with the number of draws dropped per row.

A cs() term is expanded by stats::model.matrix() like any other population-level term, which is what brms does too. A numeric predictor is one column; a FACTOR or CHARACTER predictor is its treatment-contrast dummies, one coefficient per dummy per boundary. cs(f) on a three-level factor with levels a, b and c is therefore four coefficients, named fb[1], fb[2], fc[1] and fc[2] in fixef() and bcs_fb[1] and the rest in variables(), with a as the reference level. set_prior(class = "b") covers them all and set_prior(class = "b", coef = "fb") reaches one dummy's pair; default_prior() lists one row per dummy, as brms's get_prior() does. New data are recoded against the levels and contrasts of the FIT, so a one-row newdata whose factor carries only the level it needs still predicts that level, and a level the fit never saw is refused.

A column written on BOTH sides is refused. y ~ x + cs(x) is not identified: adding a constant to the population-level coefficient of x and subtracting the same constant from every one of its category-specific coefficients leaves the likelihood unchanged. Write y ~ cs(x), which fits the same set of distributions. The same holds for a factor on both sides (y ~ f + cs(f)) and for a constant cs() column. brms builds both blocks and samples the ridge; maximum likelihood has no prior to hold it, and what frmtmb reported before the refusal was a coefficient with a standard error of 2.7e5.

The check is a RANK test, so it also catches a column that only a basis or an arithmetic expression spells differently: y ~ I(x) + cs(x), y ~ x1 + x2 + cs(x1 + x2), y ~ poly(x, 2) + cs(x) and y ~ s(x) + cs(x) are all refused, and so are two cs() terms that are the same column under different spellings. The rank comes from qr() at its default tol = 1e-7, over columns scaled to unit norm, so a merely ill-conditioned design is refused once its condition number passes about 1e7. The message names the population-level columns involved and the edit that works, which for poly(x, 2) + cs(x) is I(x^2) + cs(x) (the same maximum, one parameter fewer) and for s(x) + cs(x) is to drop one of the two terms, because there it is the smooth's own unpenalized linear column that is aliased.

The rank test cannot see an mi() or me() term: its design column is a zero placeholder at assembly, filled later with observed-or-latent VALUES rather than with a function of a data column, so there is no basis to test it against. mi(x) + cs(x) and me(x, sdx) + cs(x) are therefore accepted. me(x, sdx) + cs(x) looks genuinely identified, its standard errors being finite, and mi(x) + cs(x) is not settled. mo() has a placeholder column too, but whatever fills it is a monotone function of the predictor's categories, so mo() gets a test of its own against that category basis and y ~ mo(m) + cs(m) is refused.

cs() is not available with thres(gr = ) or under importance, and cs() on the LEFT of a bar ((cs(x) | g)) is the covariance structure of the same name, not a category-specific effect. A MULTIVARIATE fit takes cs() in any response: bf(yo ~ x + cs(f)) + bf(yo2 ~ z + cs(f)) fits, and the coefficients carry the response prefix, yo_fb[1] in fixef() and bcs_yo_fb[1] in variables().

One behavior to know when a cs() level is rare. The model frame drops an unused factor level, as brms does, so a SUBSET in which a level does not occur is fitted with one column fewer. In influence(), deleting the last row of a level therefore gives a refit whose coefficients are a subset of the full fit's: those cells stay NA, cooks.distance() is NA for that unit, and a warning names the unit and the coefficients.

Ordinal thresholds, thres()

The four ordinal families take brms's thres() addition term. y | thres(K) ~ x fits K thresholds, so the response has K + 1 categories even when the top ones are never observed. The data do not place a threshold above the highest observed category: its maximum likelihood estimate runs off toward infinity. Set a prior on class "Intercept" to hold it, as brms's default prior does. Without such a prior the fit warns.

y | thres(gr = g) ~ x gives each level of the factor g a threshold vector of its own and keeps one linear predictor for all rows. The number of thresholds of a level is the highest category observed in that level less one, or n from thres(n, gr = g), where n is a column that is constant within each level. These are brms's counts. A category above the highest one a level takes is not a category of that level, unless thres(n, gr = g) asks for it. An ordered-factor response takes its categories from the levels that occur in the data, as in brms, and thres() may not ask for more. Code the response as integers to fit categories nobody chose.

The thresholds have brms's names: Intercept[a,1] in fixef() and b_Intercept[a,1] in variables() and hypothesis(). set_prior(..., class = "Intercept", group = "a") puts a prior on the thresholds of level a only. As in brms, the design is not centered under grouped thresholds.

A row in level a falls in one of the categories 1..n_a + 1. fitted() and predict() return max(n) + 1 categories, and the probability of a category past a row's own is 0, as in brms's posterior_epred(). newdata must hold the grouping variable, with levels the fit has seen. With gr, cs() is refused, as in brms, and so is residuals(type = "osa").

A refit inside the package keeps the threshold count of the model it came from, and keeps the count of every level of thres(gr = ). frm_bootstrap(), refit(), frm_allfit(), anova(refit = TRUE), confint(method = "profile") and the autoscale pre-fit reuse the assembled design, so the count cannot move. influence() and cooks.distance() rebuild the design from a subset of the data, and carry over the count AND the response's own categories, so a subset that empties ANY category, bottom, interior or top, still fits the fitted model with one threshold that the subset no longer identifies. An ordered-factor response then gives the same table as the same data coded as integers. Two cases cannot be refit and are refused per unit, counted and reported: a groups = deletion that removes a whole thres(gr = ) level, and a data = holding a response category outside the fitted threshold layout or a thres(gr = ) level the fit never saw. A category thres(K) declared but nothing observed is inside that layout, so data = may reach it. See influence().

Calling frm() again on new data is not a refit: it counts the thresholds of that data. That is what update(newdata = ) does, because it re-evaluates the stored call, and what simulate() followed by frm() does. Write thres(K) to pin a count there.

brms's older cat(x) gives the number of CATEGORIES and is read as thres(x - 1), with brms's deprecation warning. cat() and thres() together are refused, because they set one count.

Exposures, rate()

y | rate(d) ~ x fits a count with exposure d, as in brms, for poisson(), negbinomial() and geometric(). The expected count is mu * d. Under the log link that is eta + log(d) on the linear-predictor scale, so on poisson() the model is the same as y ~ x + offset(log(d)). On the negative binomial brms multiplies the shape by d too, so the variance is mu d (1 + mu / shape), and on geometric() the density is the negative binomial of size d. d must be positive.

fitted() and predict() report the expected count, mu * d. fitted(dpar = "mu") and frm_linpred() report mu and its linear predictor without the exposure, as brms's posterior_linpred() does. newdata must hold d. Other families refuse the term, as brms does.

Response subsets, subset() and index()

In a multivariate model, y1 | subset(s) ~ x fits y1 on the rows where the logical variable s is TRUE, so each response can have its own rows. The design matrices and the smooth bases of a response come from its own rows, as brms builds them from the subsetted data. A missing value does not drop a row when every response that uses the variable leaves that row out through its own subset(); this is brms's rule. In a univariate model subset() is a filter on the rows. nobs() counts the rows of the data, as brms's does, in a univariate model too, and nobs(fit, resp = "y1") counts the rows of that response in a multivariate one. The log-likelihood's own count, which BIC() reads, is the rows of the fitted frame: the subset's rows in a univariate model, and in a multivariate one every row kept, including a row that no response's subset uses.

A predictor that reads another response, mi(x), cannot line the rows up when either response uses subset(). As in brms, name each row of x with x | mi() + index(id) and read it with mi(x, idx = ref): each row of the using response takes the value of x on the row whose id equals its ref. A value of ref with no match, a duplicated id, and mi(x) without idx across subsetted responses are refused with brms's messages. The coefficient has brms's name, mixidxEQref.

fitted(), predict(), frm_linpred() and, in frmtmb.sample, log_lik() and the posterior methods need a single resp = on such a model, as brms asks: the responses have different rows. On newdata they answer for the rows where that response's subset is TRUE, and newdata must hold the subset variable. On newdata, mi(x, idx = ) matches idx among the rows of newdata itself. conditional_effects() holds a subset variable at TRUE, as brms does.

Refused: subset() with set_rescor(TRUE), with me() terms (as in brms), with na.action = na.exclude, and, in a multivariate model, with a residual correlation term or a structured family, and vcov_cluster() on such a model.

Three differences from brms remain, all for the same reason. brms keeps every level of the whole data in a subsetted response, and frmtmb keeps the levels that the response's own rows carry, so that it creates no parameter that the data cannot place:

A grouping factor. A level with no rows adds nothing to the marginal likelihood, so the fits agree, but ranef() lists fewer levels, and |ID|-linked terms whose rows carry different levels are refused.
A factor predictor. brms keeps an all-zero column for a level that occurs only outside the subset; frmtmb drops the column with its rank-deficiency message, and the fit is the fit on the subset's rows.
An ordinal response. A category that occurs only outside the subset gives brms one more threshold than frmtmb, because frmtmb counts the categories the response's own rows take.

The Laplace approximation, and how to check it

Random effects are integrated out by the Laplace approximation, which assumes the integrand is close to Gaussian around the conditional mode; the Wald intervals confint() reports assume the log-likelihood is close to quadratic at the optimum. Both degrade in the same places: variance components estimated from few groups, and binary data in small clusters.

Three remedies are in this package. confint(method = "profile") replaces the quadratic assumption with a profile likelihood, frm_bootstrap() with a resampling distribution, and quadrature = TRUE replaces the Laplace approximation itself with adaptive quadrature (the test suite checks that fit against lme4::glmer(nAGQ = 25) and GLMMadaptive in exactly the regime where the Laplace fit is biased).

To MEASURE the violation rather than route around it, install the frmtmb.sample package and call frmtmb.sample::check_laplace(). It runs NUTS on this very objective and reports how far the posterior mean sits from the maximum-likelihood estimate in posterior standard deviations, and the ratio of the posterior standard deviation to the Wald standard error. vignette("diagnostics") works through it.

Examples

if (FALSE) { # \dontrun{
data(sleepstudy, package = "lme4")
fit <- frm(bf(Reaction ~ Days + (Days | Subject)) + gaussian(),
              data = sleepstudy)
summary(fit)
} # }
set.seed(1)
dd <- data.frame(x = rnorm(100), g = factor(rep(1:10, 10)))
dd$y <- rnorm(100, 1 + 0.5 * dd$x + rnorm(10, 0, 0.5)[dd$g], 1)
fit <- frm(bf(y ~ x + (1 | g)) + gaussian(), data = dd)
summary(fit)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + (1 | g) 
#>    Data: dd (Number of observations: 100) 
#>  Method: ML   logLik: -145.746   AIC: 299.492   BIC: 309.913 
#> 
#> Multilevel Hyperparameters:
#> ~g (Number of levels: 10) 
#>               Estimate Est.Error l-95% CI u-95% CI
#> sd(Intercept)     0.56      0.16     0.31     0.99
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.14      0.20     0.74     1.53    5.64  1.7e-08
#> x             0.67      0.11     0.45     0.89    5.92  3.2e-09
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     0.97      0.07     0.83     1.12
fixef(fit)
#>            Estimate Est.Error      Q2.5     Q97.5
#> Intercept 1.1356758 0.2013266 0.7410828 1.5302687
#> x         0.6673571 0.1126869 0.4464948 0.8882195
VarCorr(fit)
#> $g
#> $g$sd
#>            Estimate Est.Error      Q2.5     Q97.5
#> Intercept 0.5572533 0.1627123 0.2383431 0.8761635
#> 
#> 
#> $residual__
#> $residual__$sd
#>   Estimate  Est.Error      Q2.5   Q97.5
#>  0.9658535 0.07199957 0.8247369 1.10697
#> 
#> 

# distributional regression: model sigma too
fit2 <- frm(bf(y ~ x + (1 | g), sigma ~ x) + gaussian(), data = dd)
anova(fit, fit2)
#> Likelihood-ratio tests
#> Each test assumes the smaller model is nested in the larger; see ?anova.frmtmb_fit
#> 
#>                            Df  logLik    AIC  Chisq Chi Df Pr(>Chisq)   
#> y ~ x + (1 | g)             4 -145.75 299.49                            
#> y ~ x + (1 | g), sigma ~ x  5 -141.62 293.23 8.2614      1    0.00405 **
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# a nonlinear model: discover the names, edit, fit. A nonlinear body
# is undefined at the zero start, so `start` is not optional here.
nd <- data.frame(t = rep(seq(0, 10, length.out = 25), 4))
nd$y <- 8 * (1 - exp(-nd$t / 3)) + rnorm(nrow(nd), 0, 0.3)
nf <- bf(y ~ asym * (1 - exp(-t / lrc)), asym ~ 1, lrc ~ 1,
         nl = TRUE) + gaussian()
st <- par_template(nf, data = nd)
st
#> <frmtmb parameter template> starting values
#> $beta  (2)
#> asym_(Intercept)  lrc_(Intercept) 
#>                0                0 
#> $betad  (1)
#> sigma_(Intercept) 
#>         0.7961843 
#> Edit and pass back as frm(start = ) or frm_simulate(newparams = ). Parentheses are optional in names you supply.
st$beta[["asym_(Intercept)"]] <- 5
st$beta[["lrc_(Intercept)"]] <- 1
fixef(frm(nf, data = nd, start = st))
#>                Estimate  Est.Error     Q2.5    Q97.5
#> asym_Intercept 7.955791 0.08151863 7.796018 8.115565
#> lrc_Intercept  2.942713 0.08735254 2.771505 3.113921

# or let a located prior place the same starts, brms-style
fixef(frm(nf, data = nd,
          prior = prior(normal(5, 5), nlpar = "asym") +
            prior(normal(1, 5), nlpar = "lrc")))
#> Nonlinear starting values placed at the prior locations: asym_(Intercept) = 5, lrc_(Intercept) = 1. Give `start` to choose your own.
#>                Estimate  Est.Error     Q2.5    Q97.5
#> asym_Intercept 7.954524 0.08143384 7.794916 8.114131
#> lrc_Intercept  2.941385 0.08725621 2.770366 3.112405