History of frmtmb.spline • frmtmb.survival

frmtmb.spline held royston_parmar(), rp_floored() and the curve functions frm_curve(), frm_curve_deriv() and frm_curve_feature() until it was retired. The Royston-Parmar family is in this package now, and the curve functions are in frmtmb. Its NEWS follows as it stood at frmtmb.spline 0.10.2, with every heading two levels down, so that the record of the family stays with it.

frmtmb.spline 0.10.2

Runs on frmtmb 0.68.1 and later. On frmtmb 0.71.0 the standard errors of a curve on an ill-conditioned design come from a factor of the covariance. frmtmb.spline 0.10.1 refuses such a fit on frmtmb 0.71.0, so update the two packages together.

Bug fixes

The grid covariance, the derivative’s and the feature’s standard errors are formed through the factor where frm_lp_basis() returns one. On a raw polynomial of degree 8, C V C' from the whole matrix loses every digit of a band. Through the factor, a curve’s standard errors are within 1e-4 of lm()’s. A first derivative’s are within 1e-3, the accuracy of the central difference.

frmtmb.spline 0.10.1

Needs frmtmb 0.68.1, whose joint covariance is positive semi-definite on a fit that lost a standard error; its tests fail against 0.68.0.

Bug fixes

frm_curve() no longer refuses a fit whose smoothing sds ran to zero with “disagrees with frm_linpred(se.fit = TRUE) by 1 relative”. A standard error is now sqrt(max(q, 0) + max(e, 0)) from the coefficient variance q and the variance e that is not coefficient uncertainty, the order core uses. Clamping q + e as a whole gave 0 where core gave sqrt(e), whenever q came out below -e. The same arithmetic now serves frm_curve_deriv() and the height and slope of frm_curve_feature().
The rows without a standard error are read once from frm_lp_basis(): nonest (a rank-deficient design, NA) as well as se_nonest (a lost direction, NaN). The check against frm_linpred(se.fit = TRUE) also compares WHICH rows have none, and refuses with its own message when the two routes disagree.
A difference curve evaluates a row that both grids share once, so a grid differenced with itself is exactly zero on every BLAS. With OpenBLAS one estimate came out -2.2e-16.

frmtmb.spline 0.10.0

Needs frmtmb 0.68.0, for frm_lp_basis(extra_cov = TRUE) and frm_extra_cov_deriv().

Breaking changes

The "Sigma" attribute of frm_curve(), frm_curve_deriv() and a difference curve is the grid’s whole covariance, A V A' plus the covariance of the part that is not coefficient uncertainty (an exact gp()‘s kriging residual, an unseen grouping level’s draw). It was A V A' alone. Code that read it as the coefficients’ covariance gets a larger matrix wherever such a part exists, and the same one elsewhere.
frm_curve_deriv() leaves out the gp() kernel’s nugget, which is white noise at each position and has no derivative, and takes the rest of that part’s derivative from frmtmb’s frm_extra_cov_deriv() rather than by differencing it. Its .se at an exact gp()’s unseen positions therefore differs from earlier builds.

Bug fixes and changes

frm_curve(), frm_curve_deriv() and frm_curve_feature() give a NaN standard error and band to a grid row that moves along a direction the fit’s Hessian lost, and warn once per call with the number of such rows (class frmtmb_se_lost_prediction, as core’s predictions). Those rows stay out of the covariance cross-check and out of the simultaneous critical value. Before, a curve of a alone on y ~ a + b with a and b both ~ 1 + x stopped with “the assembled covariance of this grid disagrees with frm_linpred(se.fit = TRUE) by NaN relative … this package is reading the seam wrongly”, which blamed this package for a property of the fit. The curve of the mean, which the data determine, keeps its band. A difference curve takes a row as lost when either grid’s row is, which is conservative. frmtmb 0.68.0 marks those rows (frm_lp_basis()$se_nonest).
frm_curve_deriv(order = 2) past an exact gp()’s data, and at an unseen grouping level, has the right standard error. The design stencil was applied to the extra covariance as well, which divided its rounding by the step to the fourth power. On the review’s construction (y ~ gp(x), 50 points on [0, 4], grid 3.0, 4.2, 4.8, 5.5) .se was 0, 0, 0, 0.875 and is now 0.285, 0.692, 0.830, 0.959, which a 24 times larger step reproduces to 2e-5 relative; the simultaneous critical value over [3.5, 5.5] was 19.34 and is 2.50 (dev/gpby-rev-deriv2b.R, dev/gpby-p1-deriv2b-lane.txt). At an unseen level of y ~ t + (1 + t | g), linear in t, the order-2 .se was 0 to 1.495 and is 0, the true value (dev/gpby-p1-deriv2-lane.txt).
A simultaneous band over an exact gp() past its observed positions is no longer too narrow. The band drew its deviation process from A V A' and divided by a standard error that carried the kriging variance, so it was narrow exactly where it should widen. It now draws from the grid’s whole covariance, A V A' plus the kriging covariance core returns, and its "Sigma" attribute is that whole covariance. On y ~ fac + gp(x), 60 points on [0, 6], a grid on [7, 12]: the critical value was 2.022, the reviewer’s lower bound 2.330, and it is now 2.431 (dev/gpby-crit.R). Over 400 seeds of a GP drawn from the model itself on a grid running past the data, the whole-curve coverage moves from 0.7775 to 0.8100 (binomial mcse 0.0208 and 0.0196), 13 seeds covered by the new band only and none by the old one only (dev/gpby-cov-band.R); the rest of the shortfall from 0.95 is the plug-in hyperparameters’, which the pointwise band shares at 0.887 per point.
A difference curve whose two grids sit at different exact gp() positions is computed, not refused. One frm_lp_basis() call on the stacked grid returns the kriging covariance between the two grids. sp_same_latent() is gone, and with it the refusal of y ~ fac + s(x, by = fac) + gp(x) contrasted across fac.
An unseen grouping level carries its draw into every route. The simultaneous band, frm_curve_deriv() and frm_curve_feature() refused under allow_new_levels = TRUE at a level of a term re_formula keeps, because the seam returned the level’s variance one row at a time; they now read its covariance. A difference between two different unseen levels carries both draws, and one at the same unseen level cancels it.
A peak’s standard error leaves out the gp() kernel’s nugget too, for the same reason as frm_curve_deriv().
A grid with no part outside the coefficients (a smooth, a Hilbert-space gp(k = ), seen levels) allocates no extra n x n matrix: extra_cov arrives as an empty sparse matrix and is skipped. On a 2000-row s(x) grid frm_curve_deriv(nsim = 1000) peaks at the previous release’s 467 MB, and frm_curve() on a 2000-row exact gp(x) grid past its data at 255 MB against 192 MB, the difference being the kriging covariance the band now carries (dev/gpby-p1-mem.R, dev/gpby-p1-mem2.R). ### frmtmb.spline 0.9.0
Needs frmtmb 0.65.0. frm_curve() and its relatives take re_formula from core, which from 0.65.0 keeps every smooth under re_formula = NA, one indexed by a grouping factor included. The help pages say so, and test-curve.R asserts the new rule, so it fails against an older frmtmb.
frm_curve(), frm_curve_deriv() and frm_curve_feature() take allow_new_levels = FALSE and pass it to every frm_lp_basis() and frm_linpred() call they make, so the design, the estimate and the covariance check read the same rows. With TRUE at a subject the fit did not see, a model with s(t, subject, bs = "fs") gives its POPULATION curve: the factor smooth contributes a row of zeros and no variance, so the band is the uncertainty of the population smooth alone. The standard errors are identical() to the seam restricted to the intercept and s(t) columns, and agree with mgcv’s predict.gam(exclude = ) to a ratio between 1.000 and 1.004.
The curve-inference vignette failed R CMD build against frmtmb 0.65.0: it drew the population curve on a grid with no subject column at re_formula = NA, which now keeps the factor smooth. It now names an unseen subject with allow_new_levels = TRUE, and its prose no longer says that re_formula = NA drops the per-subject curves. Every number it prints is unchanged from 0.8.0.
frm_curve(contrast = , allow_new_levels = TRUE) refuses when either grid carries variance that is not coefficient uncertainty. An unseen level of a (1 | g) term loads no design column, so two different unseen levels look like one. Without the refusal their independent draws cancel and the difference standard error is exactly 0.
With allow_new_levels = TRUE and a grid row at an unseen level of a term re_formula keeps, such as (1 | g) at NULL, frm_curve(simultaneous = TRUE), frm_curve_deriv() and frm_curve_feature() refuse and name the grouping factor. One unseen level is one draw shared by every row, frm_lp_basis() returns its variance one number per row, and all three build their answer from A V A' without it. Measured before the refusal: a simultaneous critical value of 0.646 where the full covariance gives 2.146, derivative standard errors 0.40 to 0.50 of the right ones, feature standard errors about a quarter. The pointwise frm_curve() band adds that variance and is answered.
frm_curve_deriv() and frm_curve_feature() on a stored curve reuse its allow_new_levels only when the call does not supply the argument. A supplied value wins, FALSE included.

frmtmb.spline 0.8.0

rp_floored() sees a group with no events. It tested only the EVENT rows for a non-positive d(eta)/d(log t), so a random effect or a covariate on gamma1 could put a group whose rows are all censored on a falling slope, a fitted survival function that RISES with time, and the check returned zero. On dev/frailty/frailty-floor3.R’s design that happened on 4 of 6 seeds, with nlminb code 0 and nothing warning. It now also tests every censored row at its own time (an interval row at both ends) and returns the count as n_nonmonotone_censored, with the rows in the "rows" attribute as nonmonotone_censored. n_nonmonotone keeps its meaning, the event rows whose density is floored.
rp_floored(), frm_curve(), frm_curve_deriv() and frm_curve_feature() WARN on a fit with such a censored row and still answer; the fit warns as it is returned too. The warning gives the largest rise in the fitted survival, reported as the new field max_survival_rise, because the size is what a user judges. An EVENT row on a non-positive slope still refuses, as before. The line is the one flexsurv draws: dsurvspline() makes an event row there -Inf, so its fit avoids it, and never checks a censored row, so a rise there passes silently. brms’s cox() cannot rise at all (an M-spline times non-negative weights).
BREAKING: the report has two new fields, n_nonmonotone_censored and max_survival_rise, between n_nonmonotone and scale.
On ordinary right-censored designs without a cure fraction the new count never fired: 0 of 1600 simulated fits and 0 of 270 fits on nine designs from eight real datasets, at df 1 to 5, proportional hazards and gamma1 ~ x (dev/phase3a-findings.md).
The widened check DOES fire on one ordinary design: a cure fraction with a time-varying effect (gamma1 ~ arm). Past an arm’s last event its spline is fitted to censored rows only and can turn over there. Measured on two arms of 250 with 40 and 55 percent cured, df 3 to 6: 27 of 80 fits fire, every flagged row lies past its own arm’s last event, and the fitted survival rises across them by 5.2e-04 to 0.13. Those fits warn and answer; ?rp_floored says why the check fires.

frmtmb.spline 0.7.0

frm_curve()’s transform check, sp_predict_eta(), the covariance cross-check and the Royston-Parmar monotonicity check read frmtmb::frm_linpred() rather than predict(), which is brms’s predictive summary in frmtmb 0.61.0 and has no type or se.fit. Every refusal and every printed check line names frm_linpred(se.fit = TRUE) for the same reason, including frm_curve_feature()’s. Requires frmtmb 0.61.0.

frmtmb.spline 0.6.0

Requires frmtmb 0.60.0, for frm_stop() and the other condition helpers.
Requires the frmtmb release that exports frm_stop(); the frmtmb (>= 0.59.0) floor must move to it.
BREAKING: every error, warning and message that frmtmb.spline raises is classed. An error has the class c("frmtmb_spline_error", "frmtmb_error", "error", "condition"), and warnings and messages follow the same pattern, so tryCatch(frmtmb_error = ) catches any refusal. The class vector no longer contains simpleError, simpleWarning or simpleMessage. See ?frmtmb::frmtmb-conditions.
BREAKING: a value that matches none of the choices of frm_curve_feature(type =), royston_parmar(scale =) and rp_floored(action =) is refused with a frmtmb_spline_error that names the argument, the value and the choices. The old text was 'arg' should be one of .... A partial value still matches.
A refusal that frmtmb raises about royston_parmar(), such as residuals(type = "deviance"), is a frmtmb_spline_error.
The frmtmb_ps_span_warning that frm_curve() and frm_curve_deriv() raise now also has the classes frmtmb_spline_warning and frmtmb_warning.

frmtmb.spline 0.5.2

Requires frmtmb 0.59.0.
The vignette and tests follow frmtmb 0.59.0’s brms names and varcorr_matrices().

frmtmb.spline 0.5.1

The Royston-Parmar frailty is a MEASUREMENT rather than a smoke test. Item 2.5 of dev/extension-gaps-plan.md; the numbers, the seeds and the scripts are in dev/frailty-findings.md. No density, estimate or fitted value changed; what did change is the compatibility table, in the last bullet below.

(1 | centre) on mu is a shared log-normal frailty, and it is the model rstpm2::stpm2(cluster =, RandDist = "LogN") fits. The two write the same natural cubic spline in different bases, and the change of basis between them is exact to 3.0e-14 of the linear predictor’s own scale. Over 200 replicates at 2000 subjects in 40 centres, 40 percent censored: the treatment coefficient agrees to 3.2e-04 of the run’s own standard error on average and 8.8e-04 at worst, the frailty standard deviation to 1.4e-03 relative, and the two standard errors on that coefficient to 2.6e-06 relative. Recovery on the same run: beta 0.5984 against 0.6 at 93.0 percent coverage and sd 0.4905 against 0.5 at 92.0 percent, against a binomial Monte Carlo error of 0.0154.
The two do NOT report the same log likelihood, and the whole difference is the integration rule. frmtmb takes the Laplace approximation and rstpm2 takes 9-node adaptive Gauss-Hermite quadrature. Against an exact per-cluster integral written in plain R at the same parameters, the Laplace value is 0.058 units low on this design and rstpm2’s is 4.9e-06 low. It costs the ESTIMATE 5.5e-05 units of exact log likelihood, 0.1 percent of the offset itself, though frmtmb’s point was the worse one under the exact criterion on 200 of 200 replicates. The offset grows as the cluster shrinks: 0.0066 units at 200 subjects per centre, 0.058 at 50, and 1.38 at 4, where it also changes sign. Compare logLik() across frmtmb fits, not across packages.
rstpm2’s theta is the frailty’s VARIANCE. sqrt(exp(logtheta)) is what VarCorr()’s standard deviation matches, and the test asserts it as a ratio, so reading theta as a standard deviation fails.
A random effect on gamma1 is NOT a frailty and it works. gamma1 multiplies log t, so a centre deviation gives that centre a cumulative hazard t^(gamma1 + u) and a hazard ratio against another centre of t^(u - u'), which moves with time; at df = 1 it is a per-centre Weibull shape. Over 60 replicates at 2000 subjects in 40 centres with a true sd of 0.2 it recovers 0.1994 at 90 percent coverage (Monte Carlo error 0.028), and the per-centre shape error is 0.089 against 0.163 for the same fit with its own centre deviations dropped, better on 60 of 60.
A block on gamma1 ALONE is ANCHORED at t = 1, and t = 1 is a unit rather than a fact. Rescaling time by c moves log t by log c, so the same model in the new units carries a random intercept of -u log c as well. Measured on one dataset read in years, months and days, with the n_event log c a density owes a rescale: the slope-only fit’s log likelihood moves by 8.19 units and its sd(gamma1 | centre) falls from 0.2201 to 0.0370, while the paired block (1 | c | centre) on mu and gamma1 together moves by 4.3e-08 units and holds sd(gamma1 | centre) at 0.2259 in all three. ?royston_parmar says to pair the block unless the anchor is meant.
FOUND, NOT FIXED, and present in every earlier version: a group with NO EVENTS escapes rp_floored(). Its monotonicity check reads cens == 0 & detadx <= 0, so it tests observed EVENT rows only. With a random effect on gamma1 an all-censored group has no barrier holding its slope positive, since the barrier is in the density, and the score in its deviation is negative for rows past t = 1, so the fit drives that slope down. Measured on 40 centres of 10 with five followed to a common administrative time and no deaths: on 4 of 6 seeds those five come back at slopes of -0.21 to -0.31, the fit converges with no warning and a positive definite Hessian, rp_floored() reports nothing and does not refuse, frm_curve() passes it through, and one centre’s fitted survival RISES from 4.2e-50 to 0.774. ?royston_parmar now says which rows the check tests and gives the two-line check to run instead. Widening the check changes a shipped refusal and needs a measured false-alarm rate first, so it is filed rather than done.
royston_parmar x us and royston_parmar x |ID| join the compatibility table, and rstpm2 joins Suggests. The comparison tests skip without it. This is the one user-visible behaviour change in the release: frm_compat("royston_parmar", "|ID|")$status moves from untested to works, and the us row carries the frailty measurement instead of the blanket covariance-structure note.

frmtmb.spline 0.5.0

frm_curve(object, newdata, contrast = newdata2) returns the difference of two curves with pointwise and simultaneous bands, and frm_curve_feature() on it finds where they meet.

Checked against gratia::difference_smooths(). The estimate gap is provably the fit rather than the arithmetic, at 5.55e-16, and simultaneous coverage is 0.970 over 200 seeds against a binomial Monte Carlo error of 0.0154.
An exact gp() under a difference is answered where the kriging residual is the same random variable and refused where it is not, and sameness is decided on the DESIGN rather than on the numbers. A tolerant test is not merely laxer: on a mirrored grid the two variances agree to 1.1e-16 while being different draws, and differencing them anyway drops a term worth 1.3e-03 of standard error.
FOUND, NOT FIXED, and present in every earlier version: sp_sim_crit() draws from A V A' and standardizes by sqrt(diag(A V A') + extra_var), so a SINGLE grid’s simultaneous band over an exact gp() off the observed positions is at least 17 percent too narrow, exactly where it should widen. simultaneous = FALSE and every pointwise column are correct. The fix waits on a core seam.
frm_curve(object, newdata, contrast = ) returns a DIFFERENCE curve. With a second grid of the same height the reported curve is the difference of the two linear predictors and its covariance is (A1 - A2) V (A1 - A2)', so the pointwise and simultaneous bands are bands on the difference and the simultaneous one answers “is this difference anywhere other than zero” over the whole grid at once. frm_curve_deriv() and frm_curve_feature() carry the second grid through, so frm_curve_feature(type = "crossing", at = 0) on a difference curve locates where the two curves meet, with the implicit-function standard error. Measured against gratia::difference_smooths(group_means = TRUE) on the same mgcv ML fit: the difference agrees to 1.95e-08, which is 1.5 times the gap between the two packages’ own fitted curves and is provably that gap rather than the difference arithmetic, since the difference’s disagreement equals the difference of the two curves’ own to 5.6e-16. The standard errors agree to 1.6 percent, which sits between the two curves’ own 1.5 and 3.1 percent: mgcv’s Vp is conditional on the smoothing parameter and the joint precision this package inverts is not, and nothing measurable cancels. Simultaneous coverage of the true difference over the whole grid is 0.970 over 200 seeds against a pointwise band’s 0.685.

A difference is not two curves compared afterwards, and the argument exists because doing it by hand gets the covariance wrong: the two curves share coefficients, so adding their standard errors in quadrature ignores exactly the term the difference is made of. A grid differenced with itself returns zero with a standard error of identical() zero, where quadrature returns sqrt(2) times a positive number.

Refused by name: transform = TRUE, because a difference of linear predictors is not a difference of responses; two grids that load on different coefficients; a difference whose two grids sit at different exact gp() positions, since the kriging residual arrives per row and the seam returns no covariance between the grids; a new newdata passed to frm_curve_deriv() or frm_curve_feature() on a difference curve without a new contrast beside it, which would otherwise return the FIRST curve’s answer without a word; and, for a derivative or a feature, a contrast whose search variable disagrees with newdata’s, since both are moved together.

A gp() difference at ONE position is computed rather than refused. Both grids then load the same kriging residual of the same field, so it cancels exactly. Sameness is decided on the design rather than 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, since two levels of one grouping block have identical marginal variances by construction and are different draws. Measured on y ~ fac + gp(x) at n = 90 with the grid off the observed positions: the difference across fac comes back flat at 0.0618114535, which is sqrt(vcov(fit)["facB", "facB"]) to 12 digits, and moving the second grid’s x by 1e-10 refuses.

That test is STRICTER than the mathematics needs, and the refusal says so rather than blaming the gp() positions. Only the block carrying the residual has to match for it to cancel, so a contrast across fac on y ~ fac + s(x, by = fac, k = 6) + gp(x) is refused even though its gp() columns are bit-identical, because the by-factor smooth’s own columns differ. It fails closed, so that is an answer you do not get rather than one you should not trust, and the core seam below removes the need for the test altogether.

The built-in covariance check means less on this path and print() says so. predict(se.fit = TRUE) returns a marginal standard error per row and never the covariance between two grids, so the check runs on each half and cov_rel_error is the worse of the two.
KNOWN DEFECT, not introduced here and not fixed here: a SIMULTANEOUS band over an exact gp() evaluated off the observed positions is too narrow. frm_curve(simultaneous = TRUE) draws its max-deviation process from A V A' and standardizes it by sqrt(diag(A V A') + extra_var), and off the observed positions those are not the same object: an exact gp()’s kriging variance is in the divisor and not in the draw. The standardized maxima are then too small and the level quantile of them is too small with them.

This is a ONE-GRID defect. It has nothing to do with a difference and it affects any fit with an exact gp() term whose grid leaves the observed positions. simultaneous = FALSE is unaffected, and so is every pointwise column, including .se.

Bounded with public quantities, since extra_var comes from frm_lp_basis() and the divisor from predict(se.fit = TRUE). On y ~ fac + gp(x), 60 observations on [0, 6], noise 0.2, nsim = 20000: the critical value is 1.00009 of the correctly scaled one on a grid inside the data, 1.00007 at the edge, and 1.17485 on a grid that extrapolates past the observed positions. So the band is at least 17 percent too narrow exactly where it should be widening. That ratio is a LOWER bound: rescaling fixes the marginal scale and keeps A V A'’s correlation, and a kriging residual that decorrelates faster than the mean function pushes the critical value higher still.

The fix needs the conditional cross-covariance of the gp(), which no public seam returns; it is filed as a core seam in dev/diffcurve-findings.md and as a backlog entry in dev/feature-gaps.md. Until then, read a band over an extrapolating gp() grid as a lower bound on its own width.
frm_curve_feature()’s second span gate asks its question of each grid of a difference. The gate is unchanged for one grid: it fires when a column other than the search variable varies DOWN a grid. It is not asked whether the two grids of a difference differ from one another, which they always do. An out-of-span row in the contrast grid refuses exactly as one in newdata does, and every span message says which grid raised it.

frmtmb.spline 0.4.0

FIX: frm_curve_feature()’s second span check reads the grid you passed rather than the five-point stencil at the located roots. The stencil reaches eps past each root, so the old gate was true when a root happened to land within that distance of a knot and false otherwise, whatever the grid held. On a fit with two ps() terms, one row whose z left the second term’s span was refused when the crossing landed at the end of the bracket and returned in silence, with a bit-identical estimate, when it landed in the middle; and the refusal it did raise quoted z’s span while the gate that fired was t’s stencil fringe. The check now runs when, and only when, a column other than the search variable varies down the grid, which is the one thing the grid scan cannot see, and it refuses on what that question answers. A grid whose other columns are pinned to row 1, which is every grid in this package’s own examples, is unaffected and costs no extra prediction.

frmtmb.spline 0.3.0

A curve past the fitted knot span warns once under the function that drew it, the feature search refuses, and frm_curve_deriv() no longer dies on its default simultaneous band. Requires frmtmb 0.53.0 for the classed span warning.

The three curve functions say when the grid leaves a ps() knot span. They are the doors a user actually draws a curve through, and they were the ones that said nothing useful: frm_curve() and frm_curve_deriv() read the curve through frm_lp_basis(), which did not raise it at all, while frm_curve_feature() reached it through predict(), which did, ELEVEN times on one call. The first two now warn once, with the span in the message and under the name of the function that was called. The third refuses instead, below.
The span surfacing needs the frmtmb this round ships, not the 0.52.0 the Depends: floor names today. It is driven by the frmtmb_ps_span_warning condition class, which is new in core this round; against a released frmtmb 0.52.0 the class is never raised, so frm_curve() and frm_curve_deriv() fall silent again and frm_curve_feature() stops refusing, with nothing to say so. The floor rises with core’s version at consolidation.
FIX: frm_curve_deriv() past a ps() knot span works on its own DEFAULT arguments. It warned and then died in quantile.default() with “missing values and NaN’s not allowed”, naming neither the span nor the function. Past the outer knot the derivative design is exactly zero, so those rows carry a standard error of exactly zero, and the max-deviation simulation standardized their (exactly zero) deviation by it. Such a point is covered with probability one and cannot be the argmax, so it now leaves the maximization; the band is computed over the rows that carry uncertainty and the zero rows get a zero-width one. A grid on which EVERY row is past the outer knot refuses by name instead. Dropping nothing is bit-identical to the old arithmetic, so an ordinary grid is unaffected.
FIX: a grid laid exactly on the knot span is inside it. The span was checked on the widened difference stencil, which reaches a millionth of the grid’s range past both ends, so a grid ending on a knot drew a warning from frm_curve_deriv() and a REFUSAL from frm_curve_feature(), and narrowing the grid to the span, which is what the refusal tells the user to do, reproduced the refusal. Both now re-ask on the grid that was passed, and only when the widened call reported something, so an ordinary call pays nothing. The frm_curve_deriv() warning counts grid rows rather than stencil rows as a result.
BEHAVIOR CHANGE: frm_curve_feature() REFUSES a search whose bracket leaves the span, where it used to warn after the fact. A band drawn past the span is visible on the page; a peak or a crossing located past it leaves the function as a number with a standard error beside it and nothing to say which curve it came off, because the decaying partial sum has peaks and crossings of its own. The bracket is checked at the grid scan, before any root is refined, and again on the difference stencil at the located roots. Narrow newdata to the span, or refit with a larger pad =.
BEHAVIOR CHANGE: rp_floored() returns n_censored_deep where it returned n_censored_floored. Nothing about the count changed; the name did. Since 0.2.0 the censored term is scored exactly through lccdf and no row past the threshold is floored, so a field that says “floored” reports the opposite of what the slot bought. On the review’s 600-subject design the count is 1, at -log S of 55.73, and that row’s contribution is -55.7302136 exactly where the old probability-scale form gave -Inf and the floor gave -35.127363.
The package documentation no longer describes seams core has since supplied. ?frmtmb.spline-package said fit$cache$Vjoint was an internal this package reaches into, which has not been true since 0.2.0, and named “an lccdf slot, or a post-fit family hook” as fixes core owed, both of which landed and both of which this package uses. What is still missing is named for what it is: a per-row log-likelihood slot, which loo() and waic() need, and a per-group one, which frm(importance =) corrects with. The vignette("royston-parmar") paragraphs that predated the same seams are corrected with them.

frmtmb.spline 0.2.0

The three things this package had to work around are seams in frmtmb 0.52.0 now, so it stops working around them. Requires frmtmb (>= 0.52.0).

The one internal this package reached into is gone. frm_curve() and its two companions read fit$cache$Vjoint directly, which was written by an @noRd function and documented nowhere; they now call the exported frmtmb::frm_lp_basis(). ?frm_curve’s section on it changes from “The one internal this reaches into” to “The route to the covariance”.
The design rebuild is gone with it. frm_curve() used to rebuild the grid design by unit perturbation, one predict() call per contributing coefficient plus one probe per block of 24 that contributed nothing; frm_lp_basis() returns the design core already had. The predict() call count no longer depends on the number of coefficients at all, and the linearity probe that guarded the perturbation is no longer needed.
frm_curve() now works on a NONLINEAR (nl = TRUE) linear predictor, which it used to refuse. A is a Jacobian there rather than a design, and predict(se.fit = TRUE) is refused for a nonlinear predictor, so there is no second route to check against: cov_rel_error is NA and print() says the check did not run rather than reporting a passing one that never did.
frm_curve() now works on an rr() block at re.form = NULL, which it used to refuse through the covariance check. A reduced-rank block’s loadings live in theta; the perturbation could not see the derivative with respect to them and the assembled standard errors came out 27 percent away from predict(se.fit = TRUE)’s. frm_lp_basis() carries the loading columns, and the two now agree exactly.
An exact gp() term works too. Its kriging variance is not coefficient uncertainty, and frm_lp_basis() returns it separately as extra_var rather than folding it into A V A'.
royston_parmar() declares frmtmb’s new lccdf slot, so a right-censored row is scored from log S directly, in closed form on all three scales: -exp(eta) on "hazard", -log1p(exp(eta)) on "odds" and pnorm(eta, lower.tail = FALSE, log.p = TRUE) on "normal". The floor at -35.127363 is gone, and with it the region past -log S = 30 where the term was flat and its gradient exactly zero. cens() moves from conditional to works in the compatibility table.
rp_floored()’s censored count is a DIAGNOSTIC rather than a refusal. It still reports censored rows whose fitted -log S passes 19.2, because such a row is one the data barely constrain, but it no longer stops: the arithmetic is exact there now. The MONOTONICITY floor is unchanged and still refuses, because a non-positive d(eta)/d(log t) means no hazard exists and the reported likelihood is a pseudo-likelihood rather than a density.
royston_parmar() declares post$fit_check, frmtmb’s new fit-end hook, so a fit with a non-monotone row warns as it is returned rather than only when someone calls rp_floored().

frmtmb.spline 0.1.0

First release. Two things a spline needs that a fit does not give you: inference about the CURVE rather than about its coefficients, and a survival family whose parameter is a spline.

Curve inference

frm_curve() evaluates a fitted linear predictor on a grid and returns it with two intervals. The pointwise one is the usual delta-method interval. The simultaneous one covers the whole curve at once, by the max-deviation simulation of Ruppert, Wand and Carroll (2003, ch. 6): draw the curve’s own deviation process from its joint covariance, standardize by the pointwise standard error, take the largest absolute value over the grid, and use that distribution’s quantile in place of 1.96. This is the construction gratia::confint(type = "simultaneous") uses, and this package’s simulation reproduces gratia’s critical value on the same mgcv fit to inside two Monte Carlo standard errors at 10000, 100000 and 500000 draws.
The critical value is reported WITH its Monte Carlo standard error, which is the quantile standard error sqrt(p(1-p)/n) / f(q). A critical value without one invites a comparison that simulation noise alone would fail: two packages agreeing to 0.007 at nsim = 10000 have agreed to half a standard error, and that is the honest way to say it.
The comparison itself needed care. gratia standardizes a SMOOTH-ONLY deviation by smooth_estimates()’s .se, which is the FULL linear predictor’s standard error, intercept column included; the two differ by up to 13.5 percent on the test model and the critical value that comes back is 8 percent smaller than the self-standardized one. Both bands are exact, because the same divisor calibrates and scales, but the critical VALUES are not comparable unless the divisor is. The package’s own bands standardize by their own standard error.
frm_curve_deriv() gives the first or second derivative of the same curve with delta-method standard errors and the same two intervals. The derivative is taken of the DESIGN, not of the fitted values, so the estimate and its standard error describe one function.
The step size is not gratia’s. Swept against the exact derivative of a known function, a central difference bottoms out near eps = 1e-6 at 4.0e-10 for the first derivative, and near eps = 1e-4 at 2.2e-07 for the second, because a second difference divides by eps^2 and turns cancellation into the dominant error. gratia’s fixed eps = 1e-7 is 14 times off the optimum at order 1, which is negligible, and eight orders of magnitude off at order 2, which is not: on the test fit gratia’s second derivative differs from mgcv’s own same-stencil answer by 1.76, 52 percent at the grid edge, and this package’s agrees with it to 8.6e-05. So eps is scaled by the covariate’s range and set per order.
frm_curve_feature() locates a peak, a trough or a level crossing and gives its POSITION a standard error, by the implicit-function delta method: var(t*) = var(f'(t*)) / f''(t*)^2 for a stationary point and var(f(t*)) / f'(t*)^2 for a crossing. These are the numbers a movement paper reports, the time of peak velocity and movement onset, and they are usually reported without an interval because nothing hands one over.
At a stationary point the curve’s own value has a simplification worth knowing and the function reports it: the derivative of f(t*) with respect to the coefficients loses its f'(t*) dt*/dc term because f'(t*) is zero there, so the standard error of the PEAK HEIGHT is the ordinary pointwise standard error at t*, not inflated by the uncertainty in where the peak is.
Every root the grid brackets is returned, so a curve with two peaks gives two rows and a curve with none gives a zero-row answer rather than an error. “This curve does not peak in this window” is an answer.

The covariance these three needed, and the core seam

All three need the same object, the covariance of the whole grid prediction, and frmtmb exports no route to it. A penalized smooth’s wiggly part is a random-effect block even when the smooth is a population term, so the curve covariance needs the joint covariance of the fixed AND random coefficients. vcov(full = TRUE) excludes b under both of its branches, by construction: its documented invariant is that its row names are confint()’s. predict(se.fit = TRUE) forms the grid covariance internally and returns its diagonal.
So the package rebuilds it, out of two facts. The linear predictor is LINEAR in the coefficients, so the difference between a prediction and the same prediction with one coefficient raised by one is that coefficient’s design column, exactly and with nothing to tune. And the joint covariance is the inverse of the fit’s own joint precision.
Neither piece came from an exported function, so neither is trusted. Every call recomputes sqrt(diag(Sigma)) and compares it with predict(se.fit = TRUE), and refuses when they disagree by more than tol. Measured agreement on the package’s own models is 3e-15 to 9e-15 relative, which is machine precision, and it is reported by print() rather than hidden.
The covariance is core’s own, read from its cache rather than recomputed. predict(se.fit = TRUE), which this package calls anyway for the check, memoizes the inverted joint precision on the fit; the curve then subsets it. Recomputing it instead, which the first draft did, paid for a second sdreport() per call, took a Schur complement over the coefficients the curve does not touch on a densified copy of a sparse matrix, and bypassed autoscale_sdreport(). Measured at 8006 random coefficients: 114 s and 2.1 GB before, 8.2 s and 1.2 GB after, with the covariance subset itself at 0.00 s. An autoscaled fit now works rather than being refused, at 3.7e-16 against predict(se.fit = TRUE).
The cost is the joint-precision solve, not the call count. Measured at re.form = NA on a 20-point grid: design rebuild 0.01 s against predict(se.fit = TRUE) 0.29 s at 8 random coefficients, 0.07 s against 0.98 s at 2006, and 0.28 s against 6.87 s at 8006. The term the call count measures is a tenth of the cost at every size. Size a job from the number of coefficients in the FIT, not from the grid.
The design rebuild is one predict() call per contributing coefficient plus one probe per block of 24 that contributes nothing. Under the default re.form = NA a per-subject grouping block contributes nothing and is skipped in blocks rather than one coefficient at a time: measured, the vignette’s 20-subject factor-smooth model costs 32 calls against 110 random coefficients, and the count does not grow with the number of LEVELS. It is not the ideal 13 either, for two stated reasons: a chunk that straddles the boundary between a live block and a dead one is expanded whole, and a component short enough to skip the probe is expanded whether or not it contributes.
dev/spline-seam-proposal.md says what core would have to export to make the rebuild unnecessary, and what a penalized coefficient block handed to a NONLINEAR body would need beyond that.

The Royston-Parmar family

royston_parmar() writes the log cumulative hazard as a natural cubic spline in log time, which is the flexible parametric survival model of Royston and Parmar (2002). scale = "hazard" is the proportional hazards version, "odds" the proportional odds version and "normal" the probit one; all three are here, because the second and third are one switch arm each once the first is written.
It is parameterized exactly as flexsurv::flexsurvspline() parameterizes it, and that is testable rather than claimed. Taking flexsurv’s fitted coefficients on its own bc data (686 rows, 299 events) and evaluating frmtmb’s objective at them reproduces flexsurv’s log likelihood to between 1.4e-16 and 1.6e-15 relative across all three scales and 0, 1 and 3 interior knots. That is an identity, not an agreement. The two optima then match as well, and frmtmb’s is never the worse of the two.
mu is gamma0 and the rest are gamma1, gamma2, and so on, all with identity links. A formula on mu is proportional hazards; a formula on any other coefficient is a TIME-VARYING effect, which is what flexsurv spells anc =. Each is an ordinary distributional parameter, so s(), a random effect and a prior all reach them.
Knots go at equally spaced quantiles of the log UNCENSORED times, as Royston and Parmar place them and as flexsurv defaults. The quantiles need the response, so they are taken at frame assembly through family_finalize() and the family the fit carries has them baked in; knots = and bknots = pin them at construction instead. df counts interior knots plus one, so flexsurv’s k is df - 1.
The family declares an lcdf as well as an lpdf, so cens() and trunc() both reach the likelihood, and right, left and interval censoring all work. Both are registered CONDITIONAL rather than working, for the reason below. Note the coding: frmtmb’s cens() reads 0 as an observed event and 1 as right censored, which is the OPPOSITE of a Surv() status column.
rp_floored() is the check this family cannot do without, and it refuses by default. Two things in royston_parmar() are floors rather than answers and both are silent in the fitted object. This function recomputes them at the fitted parameters and stops, naming the row count, the maximum, the threshold, the reason and the remedy. frm_curve() and its two companions call it, so the documented way to inspect this family will not draw a curve off a fit whose likelihood is a floor artifact.
The first floor is the one that made the family unshippable without this check. frmtmb forms a right-censored contribution as log(1 - F(y)) on the PROBABILITY scale, and core offers a family no complementary log-CDF slot to hand back log S directly, so the scored log S carries absolute error about .Machine$double.eps / S whatever the family does internally: exact to 1.3e-13 at -log S of 10, wrong by 1.7e-04 at 30, floored at -35.127363 past 36. Past 30 the term is FLAT, its gradient exactly zero, so the optimizer prices such a row at a constant and fits the rest as if it were free. Measured: a 600-subject fit with one subject censored far beyond every event time converges without a warning and puts its treatment coefficient tens of percent out, on data flexsurv declines to fit at all. How far the reported log likelihood is from the model’s depends on the data, not on the family: a floored row contributes -35.127363 instead of its own -log S, so the shortfall is about -log S - 35 per row, and two runs of that design differing only in seed give 2.4e+03 and 2.166e+04. The threshold is one quantity for all three scales, because the error depends only on S: -log S is H on "hazard", log(1 + exp(eta)) on "odds" and -log(Phi(-eta)) on "normal", and 19.2 is where eps / S passes 1e-8.
The refusal is POST-FIT and cannot be otherwise. logLik() reads object$opt$objective directly and the family protocol has no hook that runs when a fit finishes, so nothing in this package can make logLik() or AIC() refuse on their own, and the optimizer may have walked through or stopped inside the flat region before the check is called. The real fix is a core one, an lccdf slot, and it is five sites rather than one: see dev/spline-seam-proposal.md.
Monotonicity of the cumulative hazard is NOT enforced, and flexsurv does not enforce it either. What this family does instead is refuse to answer NaN for it: where the spline’s derivative in log time goes non-positive, the log density is a large finite negative number, since NaN stops the optimizer and one NaN component poisons a mixture’s log-sum-exp. That floor is NOT inert when it is used: it turns logLik() and AIC() into a pseudo-likelihood, measured at 6 floored rows and 3952 units on a cure-fraction dataset. rp_floored() counts those rows too and refuses on them. The floor is a smooth positive part rather than a branch, because RTMB refuses a comparison on an AD type outright and has no CondExp. It changes the log density by 2.4e-15 at a derivative of 1.
That floor had to be written twice, and the measurement is what caught it. The textbook spelling, 0.5 * (u + sqrt(u^2 + eps2)), is right on paper and wrong in double precision on the side it exists for: at u = -35.75 the square root rounds to |u|, the sum cancels to exactly zero and log() returns -Inf after all. Measured, 647 of 686 rows reached -Inf that way. Recovering the small branch through (u + s)(s - u) = eps2 and ADDING it back does not help either, because adding 1.4e-16 to 71.5 loses it. The form that works SELECTS between the two cancellation-free branches with sign(), which is available on an advector and whose zero derivative is the correct one because the branches agree at u = 0.
The spline basis is written branch-free, with the truncated power spelled 0.5 * (e + abs(e)), so that it tapes if the response is ever promoted to a parameter. It is a function of the response, which is data, so nothing normally needs that; it costs one line to have it anyway.

Deliberate omissions

No post$mean_fn, so fitted() and predict(type = "response") are refused. The mean of a Royston-Parmar survival time has no closed form, and core’s cox() refuses for the same reason. Read the fitted log cumulative hazard with frm_curve() instead.
No exact basis derivative in frm_curve_deriv(). The design is rebuilt through predict(), which evaluates a basis and never differentiates one, so an exact derivative would mean reading mgcv smooth objects out of the fitted frame. That is a deeper reach into core than anything else here makes, for an error already at the tenth significant figure.
Derivatives of order 3 and higher are refused rather than offered: a third central difference divides by eps^3 and there is no step size at which it is accurate.
No accurate log survival past -log S of about 19.2, and none at all past 36. See the rp_floored() entries above: the squeeze that keeps the term finite rather than -Inf cannot also keep it accurate, because the number core asks a family for is F and the complement of a double near 1 is not representable. flexsurv computes log S directly and has no such region. The package’s answer is refusal by name rather than a silent floor, and the real fix is a core lccdf slot.