A fitted curve on a grid, with pointwise and simultaneous bands

Evaluates a fitted linear predictor on a grid of covariate values and returns it with two intervals: the usual pointwise interval, and a SIMULTANEOUS band that covers the whole curve at once.

Usage

frm_curve(
  object,
  newdata,
  contrast = NULL,
  dpar = NULL,
  resp = NULL,
  re_formula = NA,
  allow_new_levels = FALSE,
  level = 0.95,
  simultaneous = TRUE,
  nsim = 10000L,
  transform = FALSE,
  seed = NULL,
  tol = 1e-06
)

Arguments

object: A frmtmb_fit from frm().
newdata: The grid, as a data frame. Every variable the linear predictor reads must be a column, held at the value the curve is wanted at.
contrast: A second grid with the same number of rows, or NULL for an ordinary curve. With it the curve is newdata minus contrast, row by row.
dpar: Distributional parameter to read the curve off. NULL, the default, is the location parameter mu.
resp: Response name, for a multivariate fit.
re_formula: NA (the default) drops the (x | g) group-level terms, the convention mgcv and gratia plot. It KEEPS every smooth, as brms does, a smooth indexed by a grouping factor included, so on a model with s(t, g, bs = "fs") the grid has to carry g at NA too and the curve is that level's. NULL keeps every random effect, so the grid must carry every grouping column. ?frm_linpred has the full rule.
allow_new_levels: Allow a grouping level the fit did not see. FALSE, the default, refuses such a row by name. TRUE predicts it at the population level, and it is the way to read the POPULATION curve off a model with s(t, g, bs = "fs"): put g in the grid at a level the fit did not see, and the factor smooth contributes nothing. Passed to every frm_lp_basis() and frm_linpred() call, so the design, the estimate and the covariance check read the same rows. See the section "The population curve of a factor-smooth model" of frm_curve(). At an unseen level of a term re_formula keeps, such as (1 | g) at NULL, every route carries that level's draw; see the section "An unseen level of a (1 | g) term" of frm_curve().
level: Coverage of both intervals.
simultaneous: Compute the simultaneous band. FALSE returns the pointwise interval alone and skips the simulation.
nsim: Draws in the max-deviation simulation. The default 10000 puts the Monte Carlo error of the critical value near 0.013; 200000 puts it near 0.003.
transform: Return the curve and both bands through the link inverse. The bands are transformed end to end rather than rebuilt, which keeps their coverage under any monotone link.
seed: Seed for the simulation, for a reproducible band.
tol: Largest relative disagreement with frm_linpred(se.fit = TRUE) the assembled covariance may show before the call refuses.

Value

A data frame of class frmtmb_curve: the columns of newdata, .eta where there are several curves, then .estimate, .se, .crit, .lower_ci, .upper_ci, and when simultaneous = TRUE also .crit_sim, .lower_sim and .upper_sim. The grid covariance, A V A' plus the variance that is not coefficient uncertainty, is the "Sigma" attribute, the fit is the "fit" attribute, and "check" carries the covariance agreement and the predict() call count.

print() shows the critical values, the covariance check and the first six rows. Its digits = is passed to print.data.frame() for those rows: significant digits, getOption("digits") when not given. The header lines keep their own precision.

Details

A pointwise interval is the wrong tool for the question a curve usually raises. "Is the velocity above zero at 300 ms" is pointwise; "does this curve have the shape I claim" is a statement about every point at once, and a 95 percent pointwise band covers the whole curve far less than 95 percent of the time. The simultaneous band is the max-deviation simulation of Ruppert, Wand and Carroll (2003, ch. 6): draw the curve's own deviation process from its joint covariance, standardize each draw by the pointwise standard error, take the largest absolute value over the grid, and use the level quantile of those maxima in place of qnorm(0.975). It is the construction gratia::confint(type = "simultaneous") uses on an mgcv fit, and this simulation reproduces gratia's critical value inside its Monte Carlo error.

The route to the covariance

A penalized smooth's wiggly part is a random-effect block in the fitted objective even when the smooth is a population term, so the covariance of a curve needs the joint covariance of the fixed AND random coefficients.

All of it comes from ONE call to frmtmb's exported seam, frm_lp_basis(), which returns the design A over the coefficient vector, the joint covariance V at exactly the rows A's columns sit at, and the variance that is NOT coefficient uncertainty (a new grouping level's draw, an exact gp()'s kriging residual) as a separate covariance over the grid rows, extra_cov. Sigma is A V A' + extra_cov, and the pointwise standard error is sqrt(diag(Sigma)). The simultaneous band draws from the same Sigma it standardizes by.

Up to frmtmb 0.51.0 there was no such seam. frmtmb.spline, where this function lived until it moved into frmtmb, rebuilt A by unit perturbation, one predict() call per contributing coefficient, and read V out of fit$cache$Vjoint, which was an internal with no precedent to point at. One call replaced the reconstruction and the reach, and what the covariance check verifies has changed from "the reconstruction reproduced core's number" to "the seam is being read correctly".

The check itself stays. Every call compares that pointwise standard error with frm_linpred(se.fit = TRUE), a second route through core, and refuses when the two disagree by more than tol. The measured agreement is in the "check" attribute and is reported by print(); it is at machine precision, 2.2e-16 relative on the curve-inference vignette's model.

The one case with nothing to check against is a nonlinear (nl = TRUE) body: frm_linpred(se.fit = TRUE) is refused there, so frm_lp_basis() is the only route to the number and cov_rel_error is NA. print() says so rather than reporting a check that never ran.

Cost

What this call costs is dominated by ONE thing: the single frm_linpred(se.fit = TRUE) check call, inside which core inverts the fit's joint precision matrix over EVERY coefficient, including the ones this curve does not touch. Measured at re_formula = NA on a 20-point grid, one process each:

s(x, k = 10), 8 random coefficients: 0.29 s.
s(t, k = 8) + (1 + t | subject), 1000 subjects and 2006 random coefficients: 0.98 s.
the same over 4000 subjects, 8006 random coefficients: 6.87 s.

The design rebuild that used to sit beside those figures, and that was a tenth of them at every size, is gone: frm_lp_basis() returns the design core already had, so the predict() call count no longer depends on the number of coefficients at all. The joint-precision solve is now the whole cost, it is paid once because core memoizes it, and it grows with the total number of coefficients in the fit rather than with the grid.

A difference curve

contrast is a second grid of the same height. The curve returned is then the DIFFERENCE of the two linear predictors, row by row, and its covariance is (A1 - A2) V (A1 - A2)', so both bands describe the difference and the simultaneous one answers "is this difference anywhere other than zero" over the whole grid at once. It is the quantity gratia::difference_smooths(group_means = TRUE) reports for a factor-by smooth, and on the same mgcv fit the two agree. Name that argument when you compare: gratia's default, group_means = FALSE, zeroes the intercept and the parametric group columns and reports the smooth-only difference, which is a different quantity. On the fixture test-difference.R uses it is 0.52 away.

A difference is NOT the difference of two calls to this function. The two curves share coefficients, so their covariance is what the difference is made of, and adding two standard errors in quadrature would ignore it.

What the difference path cannot do, and refuses by name:

transform = TRUE. A difference of linear predictors is not a difference of responses under any link but the identity, so there is nothing to map it through.
Two grids that load on different coefficients.
Two grids that load DIFFERENT draws of a latent field whose variance is not coefficient uncertainty. See the next section.
allow_new_levels = TRUE when either grid carries such a variance. An unseen level loads no column of the design, so two different unseen levels look the same as one, and treating them as one would cancel two independent draws. Measured on (1 | subject): the difference standard error came back as exactly 0 where it is at least 0.0515. An unseen level of an fs smooth carries no such variance and is not refused.

The covariance check also means less here, and print() says so. frm_linpred(se.fit = TRUE) returns a marginal standard error per row and never the covariance between the grids, so the check runs on each half and cov_rel_error is the worse of the two: what it licenses is that both designs were read correctly.

An exact `gp()` under a difference

An exact gp() evaluated off the observed positions carries a kriging residual that is not coefficient uncertainty. frm_lp_basis() returns its variance one number per ROW and returns no covariance BETWEEN two grids, so var(g1 - g2) has no public route.

It needs none in the ordinary case. A contrast taken across a factor at ONE gp() position leaves both grids loading the same residual, so it cancels exactly and the difference is (A1 - A2) V (A1 - A2)' with nothing left over. That case is computed rather than refused.

Sameness is decided on the design and not on the numbers: the two grids must agree bit for bit on EVERY column of A outside the fixed effects. Equality of the variances would not be enough, because two levels of one grouping block have identical marginal variances by construction and are different draws.

That test is stricter than the mathematics needs, and it is worth knowing where the extra strictness bites. Only the block carrying the kriging residual has to match for the residual to cancel, but the test asks it of every latent column, so a contrast across fac on y ~ fac + s(x, by = fac) + gp(x) is REFUSED even though the gp() columns are identical: the by-factor smooth's own columns differ, which is what a by-factor smooth is for. It fails closed, so the cost is an answer you do not get rather than one you should not trust.

The rest is refused because the SEAM cannot supply it, not because the mathematics is missing. The conditional cross-covariance is k(x1, x2) - Xr1 K Xr2', and core forms every piece of it while predicting, but reduces the result to one variance per row before the seam returns. Hold every latent term equal between the grids and contrast a fixed effect, or read the two curves separately.

The population curve of a factor-smooth model

re_formula = NA keeps a factor smooth such as s(t, subject, bs = "fs"), as brms does, so on v ~ s(t) + s(t, subject, bs = "fs") it gives one subject's curve and the grid must name a subject in every row. To read the POPULATION curve, name a subject the fit did not see and pass allow_new_levels = TRUE. mgcv's fs basis gives an unknown level a row of zeros, so the factor smooth contributes nothing:


pop <- factor("new", levels = c(levels(d$subject), "new"))
frm_curve(fit, newdata = data.frame(t = tt, subject = pop),
          allow_new_levels = TRUE)

The band is then the uncertainty of the intercept and the s(t) coefficients alone. frm_lp_basis() adds no variance for an unseen fs level: its extra_var is exactly zero and the factor smooth's columns of A are exactly zero. Measured on the curve-inference vignette's model, 7200 rows over 20 subjects: the standard errors are identical() to the same seam restricted by hand to the intercept and s(t) columns, and agree with mgcv::predict.gam(exclude = "s(t,subject)", se.fit = TRUE) on an mgcv ML fit of the same model to a ratio between 1.000 and 1.004.

Two unseen fs levels are the population curve twice, so a contrast between them is exactly 0 with a standard error of 0. That is NOT the difference between two new subjects, which this function cannot report.

An unseen level of a `(1 | g)` term

An unseen level of a bar term is different from an unseen fs level. Where re_formula keeps the term, the unseen level is ONE draw of its block's effect shared by every row at that level, so the covariance of the grid is A V A' plus that draw's covariance, which frm_lp_basis(extra_cov = TRUE) returns. Every route reads it: the pointwise and simultaneous bands, frm_curve_deriv(), where a random intercept moves the curve and not its slope and a random slope adds its own variance to the slope's, and frm_curve_feature(). A contrast at the SAME unseen level cancels the draw; one between two different unseen levels carries both. Before frmtmb returned the covariance, the bands drew from A V A' alone, and the simultaneous critical value was 0.646 on (1 | subject) and 0.615 on (1 + t | subject) on the curve-inference vignette's data, against 2.146 and 2.472 from the full covariance, so those routes refused.

re_formula = NA drops the term and gives the population curve on every route.

An exact `gp()` past its positions

An exact gp() evaluated at a position the fit did not see carries the field's kriging residual, one draw of the field shared by every unseen row, and its covariance enters Sigma the same way. The band therefore widens where the grid leaves the observed positions. The 1e-12 nugget frmtmb puts on the kernel's correlation is white noise at each position, which has no derivative: frm_curve_deriv() and a peak's standard error leave it out.

Several predictors in one call

With no dpar, an ordinal family with a cs() term has one curve per threshold, eta + cs_k, and a family with several location predictors has one per predictor: the predictors frm_linpred() returns as columns, which is what brms's posterior_linpred() returns as layers. The result holds every row of the grid once per curve, the column .eta naming it, in the order as.numeric() reads that matrix. The simultaneous band covers all the curves at once, and "Sigma" is over all the rows, so the covariance between two thresholds is there. frm_curve_deriv() has the same layout, and its derivative in a cs() variable differs by threshold. frm_curve_feature() searches each curve on its own and adds .eta to its table. dpar = "mu" reads the shared part, one curve, as in brms.

A me() term

newdata must hold the measurement SD. The curve reads the observed value, and the band carries the variance of a new row's latent value, as frm_linpred() does. So .se is not coefficient uncertainty alone. That variance is white noise between rows, so a derivative and the location of a feature leave it out. Under me(gr = ) a level holds one value, so its noisy variable has no derivative and no feature within a level, and both refuse.

With contrast, two rows at one x and one sdx, to all.equal()'s tolerance, read one draw, and it cancels, as a kriging residual at one position does: a contrast in z has the standard error of z's coefficient, 0.051 at a coefficient standard error of 0.051. brms has no such curve. A difference of its draws of the two grids reads a draw per row and has a standard deviation of 0.96 on that fit. Under me(gr = ) one level is one draw in one grid in both packages.

Two readers answer two questions, and each is right for its own:

frm_curve(contrast = ) answers for one unit with one unknown x*, reported as x, compared at two values of z. x* is not in that difference.
emmeans() contrasts answer for two new rows, each with its own measurement error, as brms's draws do. They carry both draws: 0.97 on that fit.

For the effect of z, read frm_curve(contrast = ) or the coefficient. For the effect of x, read frm_curve_deriv(), the slope's 0.053 on that fit: a contrast in x itself has two values of x, so it reads two draws, 0.998.

Past a `ps()` knot span

A ps() basis is a partition of unity only between its frozen outer knots. Past them it is a partial sum that decays to zero, so a curve drawn there bends smoothly to whatever the rest of the body gives, which is exactly the shape a reader does not question. predict(newdata = ) says so, and so does frm_lp_basis(), the seam this function reads. It is surfaced again here, ONCE per ps() term per call and carrying the span, so that the sentence names the function you called: this one reads the seam on the grid, but frm_curve_deriv() reads it on a three-point difference stencil and frm_curve_feature() on a five-point one, and core counts the rows it was handed.

frm_curve_feature() REFUSES instead of warning. A band past the span is visibly wrong on the page; a peak located past it is a number with a standard error beside it and nothing to give it away.

References

Ruppert, D., Wand, M. P. and Carroll, R. J. (2003) Semiparametric Regression. Cambridge University Press, ch. 6.

Examples

set.seed(1)
dd <- data.frame(x = sort(runif(200)))
dd$y <- 2 * sin(pi * dd$x) + rnorm(200, 0, 0.4)
fit <- frm(bf(y ~ s(x, k = 8)),
           family = stats::gaussian(), data = dd)
cv <- frm_curve(fit, newdata = data.frame(x = seq(0, 1, length.out = 25)),
                nsim = 2000)
head(cv[, c("x", ".estimate", ".se", ".lower_ci", ".lower_sim")])
#> <frmtmb curve> , 6 grid points, level 
#>   critical value: pointwise NULL
#>   covariance NOT checked: frm_linpred(se.fit = TRUE) is refused for a nonlinear predictor, so there is no second route to compare against
#>            x .estimate        .se  .lower_ci  .lower_sim
#> 1 0.00000000 0.1320144 0.15875289 -0.1791356 -0.32995857
#> 2 0.04166667 0.3413065 0.11074390  0.1242525  0.01904038
#> 3 0.08333333 0.5478148 0.07848158  0.3939937  0.31943241
#> 4 0.12500000 0.7472178 0.06874628  0.6124776  0.54716530
#> 5 0.16666667 0.9386515 0.06827916  0.8048268  0.73995831
#> 6 0.20833333 1.1261967 0.06628100  0.9962883  0.93331811