marginaleffects on frmtmb fits

frmtmb registers the methods that marginaleffects needs: get_coef(), set_coef(), get_vcov() and get_predict(). So predictions(), comparisons(), slopes(), hypotheses() and their avg_ forms work on a frmtmb_fit. The standard errors come from the delta method, which marginaleffects applies with a numerical Jacobian.

Value

This page documents the marginaleffects methods and returns nothing.

The parameters marginaleffects perturbs

get_coef() gives every estimated coefficient, then the parameters that a prediction reads outside them: the thresholds and cs() coefficients of an ordinal model, and the free coordinates of each mo() simplex. If a mo() term has a category that no row is in, the data do not determine its simplex. Then get_coef() gives the cumulative weights at the categories that rows are in (zeta<j>_cum<c>), which the data determine. A prediction at an empty category gets a standard error of NaN, with one warning per fit, as fitted() gives it. The other standard errors stay finite. A mo() term on one value is 0 on every row, so no prediction reads its coefficient, and that coefficient gets a variance of 0 in get_vcov().

On a fit with a group-level term, a smooth, gp() or hsgp(), get_coef() then gives the group-level effects at their estimated modes, under the names frm_lp_basis() gives them (b.<term>.<level>.<coef>), and the covariance parameters that a field's coefficients read (theta.<name>, as for an exact gp(), whose field is sdgp L z). A prediction reads the modes, as fitted() does. Under re_formula = NA the group-level effects are not in get_coef(); the coefficients of a smooth stay, because re_formula does not drop a smooth. Inside hypotheses() called on the fit itself, a test of the model's own coefficients, get_coef() and get_vcov() give the population-level and distributional ones alone, as hypotheses() reads those of a brmsfit: "b2 = 0" names the same coefficient as before, and joint = TRUE tests those coefficients and no group-level effect. A function given as hypothesis is not such a test: it gets the model, and what it returns is differentiated with the group-level effects in, as the same call outside hypotheses() is. A latent residual's values are not in get_coef(), though under cov = TRUE a new row continues from them: there is one per fitted row, which would cost one prediction each, and the warning of class frmtmb_interop_se_part names the latent residual on such a fit.

marginaleffects takes its Jacobian with one prediction per parameter, so the time grows with the number of levels. On a poisson (1 + x | g) fit with 40 rows per level, avg_slopes() took 0.30 s at 30 levels and 16.3 s at 500 levels, where it took 0.03 s and 0.09 s with the group-level effects held fixed. vcov = FALSE gives the estimates alone, at the cost of a few predictions, and re_formula = NA drops the group-level effects from get_coef().

The vcov argument

marginaleffects passes its vcov argument to get_vcov():

TRUE or NULL: The model-based covariance (the default).
FALSE: No standard errors.
A matrix: A covariance over get_coef(), matched by its names. An unnamed matrix is read in the order of get_coef(), with marginaleffects' warning. A matrix over the whole outer parameter vector, as vcov(fit, full = TRUE) and vcov_cluster(full = TRUE) give it, is cut down to the rows of get_coef(). A REML fit has no fixed effects in that vector, so such a matrix is refused there. If the matrix has no finite value in a row of get_coef(), as vcov(fit, full = TRUE) has for a parameter the fit does not determine, a warning names the row: every standard error is then NaN. On a fit whose get_coef() holds group-level effects, a matrix over the outer vector is carried into them: their block is their covariance given the outer parameters plus the given matrix carried through the modes' derivative in those parameters, as the default carries the model's own. A matrix over the rows of get_vcov(fit) is read as given. A matrix over the coefficients alone holds the effects at their modes, with a warning of class frmtmb_interop_se_part.
A function: A function of the fit that returns one of these matrices.
"HC0": The sandwich with each row as its own cluster (vcov_cluster(type = "CR0")).
"HC1" or "stata": HC0 times n / (n - k), as for lm(), with k the number of population-level coefficients (those of every linear predictor but the distributional parameters') plus the steps of each mo() simplex that the data determine; a mo() term on one value counts nothing, as lm() does not count a column it aliases. So k is the column count of lm() on factor() of the mo() variable. marginaleffects maps "stata" to HC1 for lm().
A one-sided formula, ~ g: The cluster sandwich over g, vcov_cluster(type = "CR1"), with the factor G / (G - 1).

On a fit with a group-level term, "HC0" and "HC1" are refused (each row is its own cluster, which every grouping level crosses), so the sandwich there is a cluster formula over a grouping factor the effects nest in, ~ g. It covers the outer parameters, and it is carried into the group-level effects as a matrix over the outer vector is: only the part of their variance that comes from the outer parameters becomes robust, and their conditional variance given those parameters stays model-based. The bread of these sandwiches is the inverse of the observed information. sandwich::vcovHC() and vcovCL() use the expected information of a glm. The two are equal for a gaussian model and for a glm with its canonical link, so there "HC0" equals vcovHC(type = "HC0"), and ~ g equals vcovCL() for a glm (for lm(), vcovCL() also multiplies by (n - 1) / (n - k)). For another link they differ: on 300 rows of a Gamma(link = "log") model, the standard error of the slope is 1.128 times vcovHC()'s. The sandwiches are formed on the directions the data determine, so they stay finite where the default covariance does: on a fit with an empty mo() category, a mo() term on one value or a variance at its boundary.

The other strings that marginaleffects knows are refused, each with its reason. "HC2" to "HC5" and "robust" (which marginaleffects maps to HC3) weight each row by its leverage, the diagonal of a hat matrix, which frmtmb does not compute. "HAC", "NeweyWest" and "kernHAC" need a time order of the rows. "satterthwaite" and "kenward-roger" are denominator degrees of freedom for a linear mixed model fitted by REML. vcov_cluster() refuses a fit whose random effects cross the clusters, and so "HC0" and "HC1" refuse such a fit too. It also refuses a residual correlation term (ar(), ma() or arma(), either form, and the other residual covariance terms), whose density is joint over each group: under one, "HC0", "HC1", "stata" and a cluster formula are refused, vcov = TRUE gives the model-based standard errors, and frm_bootstrap() is the route to a robust one.

What the standard error holds

marginaleffects forms a standard error from the parameters get_coef() holds, by the delta method. With the group-level effects and the field parameters among them, get_vcov() gives the joint covariance of the fixed effects, the variance parameters and the predicted effects, which is the covariance fitted.frmtmb_fit() and frm_linpred(se.fit = TRUE) read. So the standard error of a level from predictions() is fitted()'s Est.Error, and an average, a slope or a contrast holds the same uncertainty. On 168 fits with (1 | g) or (1 + x | g) on four families, the standard errors of 2016 levels were within 1e-5 of fitted()'s, and those of avg_predictions(), avg_slopes() and avg_comparisons() were within 25 percent of the standard deviation of 150 draws from the joint covariance on 164 of the fits. On the other 4, a contrast of a group-level covariate at 8 groups, the delta method itself is that far from the draws.

Through 0.79.1 get_coef() held the population-level and distributional coefficients alone, and the standard error held the group-level effects and a smooth's coefficients fixed. On the same 168 fits, avg_predictions() was more than 25 percent too wide on 77 percent of them (median 2.03 times) and avg_slopes() on 52 percent (median 1.30 times); on a gp() field every level read one value.

What no parameter holds is left out, with a warning of class frmtmb_interop_se_part that names it:

A me() term. Every row reads as a new row, whose latent value brms draws from N(x, sdx). The standard error of a level reads 0.07 where frm_linpred() reads 1.11. A slope, a contrast or a test of a coefficient on an identity link needs no such draw and is right: -0.491, standard error 0.0485. get_vcov() warns.
A latent residual (an autocorrelation term on a family with no residual). A new row's latent value is in no parameter, and the fitted latents it continues from under cov = TRUE are left out of get_coef() (see above). get_vcov() warns.
A fit whose objective integrates the group-level effects out (quadrature = TRUE or importance =). Its covariance has no rows for them, so they are held at their modes, as fitted() holds them, and says. get_vcov() warns.
A grouping level the fit did not see (allow_new_levels = TRUE), and an exact gp() at a position the fit did not see. fitted() adds the block's variance and the field's kriging variance there. get_predict() warns once per call, from the first prediction of the Jacobian, which marginaleffects does not muffle, and counts the rows of each kind. A kriging variance below 1e-4 of the field's own variance is not counted; at the data and between data points it was below 1e-7 in the cases measured, and the level's standard error equals fitted()'s to 4 digits there.
A matrix passed as vcov that covers the coefficients alone (see above).

On the response scale of a me() fit, the estimate is the mean over the draw, as fitted.frmtmb_fit() reads it. vcov = FALSE gives the estimates alone, with no standard error and no warning.

marginaleffects' Jacobian is a finite difference; fitted() and frm_lp_basis() differentiate exactly. Far from its data an exact gp() is strongly curved in its length scale, and the two part: on a gp(x) + (1 | g) fit with data on 0 to 6, the forward difference read 0.9996, 1.0217, 1.0749, 1.0232 and 0.6265 times fitted()'s standard error at x = 7, 8, 9, 10 and 12, and numderiv = "fdcenter" read 1.0006, 1.0081, 1.2061, 1.1244 and 0.6441, no closer. Read a level there from fitted(). A slope of an exact gp() at a small step eps reads the field's nugget, which grows as 1 / eps: on one fit with x on 0 to 6 its standard error was 0.0308 at the default step and at eps = 1e-3, 0.0328 at 1e-4 and 0.1124 at 1e-5, for the same estimate.

brms's prediction arguments

re_formula, allow_new_levels, dpar (also under brms's name nlpar) and resp go to frm_linpred(), as they go to brms's prediction for a brmsfit. re_formula = NA or ~ 0 gives the population-level prediction, and NULL (the default) keeps every group-level effect. The arguments that marginaleffects passes for other model classes and frmtmb does not take are refused by name: lme4's re.form and allow.new.levels, brms's sample_new_levels, ndraws, draw_ids and incl_autocor, and glmmTMB's zitype. marginaleffects 1.0.0 keeps a fixed list of the arguments that it knows for each model class, and a package cannot add to that list. So marginaleffects warns that re_formula is "not known to be supported" for a frmtmb_fit, and then passes it on. The argument has its effect.

Examples

# \donttest{
if (requireNamespace("marginaleffects", quietly = TRUE)) {
  set.seed(1)
  d <- data.frame(x = rnorm(120), g = factor(rep(1:12, each = 10)))
  d$y <- 1 + 0.5 * d$x + rnorm(12)[d$g] + rnorm(120)
  fit <- frm(y ~ x + (1 | g), data = d)
  nd <- data.frame(x = c(-1, 1), g = factor(1, levels = 1:12))
  # the standard error holds the group-level effects, as fitted() does
  marginaleffects::predictions(fit, newdata = nd)
  fitted(fit, newdata = nd)
  # population-level predictions; marginaleffects warns that it does not
  # know re_formula for this class, and passes it on
  suppressWarnings(marginaleffects::predictions(fit, newdata = nd,
                                                re_formula = NA))
  # a cluster-robust covariance
  fit0 <- frm(y ~ x, data = d)
  marginaleffects::avg_slopes(fit0, variables = "x", vcov = ~ g)
}
#> 
#>  Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
#>     0.553      0.113 4.89   <0.001 20.0 0.332  0.775
#> 
#> Term: x
#> Type: response
#> Comparison: dY/dX
#> 
# }