The design of a linear predictor over the coefficient vector

frm_linpred(se.fit = TRUE) builds a matrix A with one row per prediction and one column per contributing coefficient, forms A V A' and keeps only its diagonal. Every delta-method quantity over a fitted curve needs the whole thing: a contrast between two grids, an average marginal effect with a correct standard error, a simultaneous band, a derivative, the time of a peak. This returns the pieces so that an extension does not have to rebuild A by perturbation, one predict() call per coefficient.

Usage

frm_lp_basis(
  object,
  newdata = NULL,
  dpar = NULL,
  resp = NULL,
  re_formula = NULL,
  allow_new_levels = FALSE,
  extra_cov = FALSE
)

Arguments

object: A fitted frmtmb_fit.
newdata: Data to build the design at, or NULL for the training data.
dpar, resp: The distributional parameter and response to take the linear predictor of. Both default the way predict() defaults them.
re_formula: NULL keeps every random effect, NA drops the group-level terms and keeps every smooth, a one-sided formula keeps the terms it names. See frm_linpred() for the full rule.
allow_new_levels: Whether a grouping level the fit never saw is allowed.
extra_cov: If TRUE, also return extra_cov, the full covariance over the rows of the variance in extra_var.

Value

A list with

eta: the linear predictor, exactly frm_linpred(type = "link").
A: n x p; d eta / d coef.
coef_pos: length p; the rows of V the columns of A belong to, in V's own order.
V: p x p; frm_joint_cov() subset to coef_pos.
coef_names: length p; the labels of those rows.
extra_var: length n; variance that is NOT coefficient uncertainty, kept separate rather than folded into A V A' because an exact gp()'s kriging variance and a new grouping level's marginal variance are not.
nonest: length n; rows that load on a direction the rank-deficient design could not identify.
extra_cov: only with extra_cov = TRUE: n x n; the covariance between rows of the variance in extra_var, whose diagonal it is, exactly. A sparse Matrix::Matrix() when its blocks cover little of the grid (all zero when nothing contributes), a dense matrix otherwise.

var(eta) is rowSums((A %*% V) * A) + extra_var, and frm_linpred(se.fit = TRUE) is written that way. The covariance of the whole prediction is A V A' + extra_cov.

Beside a latent residual (ar(), ma(), arma(), cosy() or unstr() on a family with no residual), coef_pos leaves out every column that no row of A loads: a latent, a level the rows do not name, a coefficient whose covariate is 0 on every row. Such a column adds nothing to A V A'. A new row loads no latent, so A and V stay small. The fitted rows load every latent, and there A and V are dense, about 16 n^2 bytes for n rows (1 GB at 8000 rows). A V of more than getOption("frmtmb.latres_dense_rows", 20000)^2 numbers is refused. frm_linpred(se.fit = TRUE) reads the fitted rows from the sparse factor of the joint precision instead, at no such cost.

Under a cov = FALSE autocorrelation term eta is the predictor without the term, as brms's emmeans() reads it (incl_autocor = FALSE), and frm_linpred() adds the one-step part from the earlier residuals.

With no dpar on a predictor that frm_linpred() returns as several columns (a cs() term, several location predictors), every row of the grid comes once per column, the column slowest, so eta is as.numeric() of that matrix, and the list also has layers, the column names. Under a cs() term A has the term's coefficients beside the shared design. dpar = "mu" is the shared part.

On a fit whose fixed design is ill-conditioned (a raw polynomial of high degree, a covariate far from zero) the list also has

V_root: p x k; a factor of the coefficient covariance.
V_rest: p x p; the part of V that V_root does not carry, zero without random effects.

with V = V_rest + V_root V_root'. Use rowSums((A %*% V_root)^2) + rowSums((A %*% V_rest) * A) for var(eta), as frm_linpred() does. V is accurate entry by entry, but A V A' formed from it in double precision is not: on a raw polynomial of degree 8 it can be wrong by a factor of 3.

Details

The name follows emmeans::emm_basis(), which is the same idea for the fixed block alone.

Variance that is not coefficient uncertainty

Three sources reach extra_var. A new row of a me() term adds the variance of its latent value's draw from normal(x, sdx); rows at one level of me(gr = ) read one draw and covary fully, and other rows are independent. An exact gp() evaluated at a position the fit did not see adds the field's conditional (kriging) covariance given its values at the fitted positions, k(x, x') - k(x, X) K^-1 k(X, x'); rows at an observed position carry none. Every unseen row of one gp() term loads the SAME draw of the field, so those rows covary, and a simultaneous band, a difference between two grids or a derivative needs the whole matrix rather than its diagonal. A grouping level the fit never saw (allow_new_levels = TRUE) adds its block's covariance: rows at the same unseen level load one draw of its effect and covary fully, and rows at two different unseen levels load independent draws.

For a difference between two grids, call this once on the stacked grid rbind(newdata, newdata2): the off-diagonal block of extra_cov is the covariance between them. For a DERIVATIVE of the curve, do not difference extra_cov: a finite difference of a covariance divides its rounding by the step to the power of twice the order. frm_extra_cov_deriv() returns that covariance from its sources.

A nonlinear body

For a nonlinear predictor A is a JACOBIAN rather than a design, and it is computed by taping the body against the coefficients it reaches through. That includes a ps() block, whose coefficients enter the body through a spline evaluated at an argument the parameters move. frm_linpred(se.fit = TRUE) stays refused for a nonlinear predictor; this is the route.

The Jacobian is exact, and the delta method built on it is still a first-order approximation, which for a warped curve is a stronger assumption than it is for a linear one. allow_new_levels = TRUE is refused there, and so is a contributing exact gp(), because neither variance has a chain rule through the body that has been measured.

A reduced-rank block

An rr() block's loadings live in theta, so a design over (beta, b) alone is INCOMPLETE. A carries the loading columns too, through rr_jacobians(), and coef_pos names their theta rows; a caller therefore gets the whole delta method rather than discovering a missing piece.

Examples

set.seed(1)
dd <- data.frame(x = rnorm(120), g = factor(rep(1:12, each = 10)))
dd$y <- rnorm(120, 1 + 2 * dd$x + rnorm(12, 0, 0.5)[dd$g], 0.4)
fit <- frm(bf(y ~ x + (1 | g)), data = dd)
nd <- data.frame(x = c(-1, 0, 1), g = factor(1, levels = levels(dd$g)))
lb <- frm_lp_basis(fit, newdata = nd, re_formula = NA)
str(lb$A)
#>  num [1:3, 1:2] 1 1 1 -1 0 1
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : chr [1:3] "1" "2" "3"
#>   ..$ : chr [1:2] "(Intercept)" "x"
lb$coef_names
#> [1] "beta.(Intercept)" "beta.x"          

# the covariance of the WHOLE grid, which frm_linpred() reduces to
# diagonal
Sigma <- lb$A %*% lb$V %*% t(lb$A)
all.equal(sqrt(diag(Sigma)),
          frm_linpred(fit, newdata = nd, re_formula = NA,
                      se.fit = TRUE)$se.fit)
#> [1] TRUE