Parametric bootstrap

Simulates nsim response vectors from the fitted model (by default with new group effects each draw; see re_formula), refits the model to each through refit() (warm-started, no re-parsing), and collects FUN of every refit. Draws whose refit fails are kept as NA rows; draws whose optimizer does not report convergence are kept but flagged.

Usage

frm_bootstrap(
  fit,
  FUN = function(f) fixef(f, flatten = TRUE),
  nsim = 500,
  seed = NULL,
  re_formula = NA
)

Arguments

fit

A frmtmb_fit for a univariate model.

FUN

Function of a frmtmb_fit returning a numeric vector. Default: fixef(f, flatten = TRUE), the fixed effects named as vcov() names its rows, so a bootstrap standard error lines up with the Wald one by name.

nsim

Number of bootstrap draws.

seed

Optional seed.

re_formula

Which group-level terms stay at their fitted values; the rest are redrawn in every replicate, from the fitted covariance. It is read exactly as predict() and simulate.frmtmb_fit() read it. The default NA (or ~0, ~1) redraws every group-level effect, in the manner of lme4's bootMer(use.u = FALSE); NULL holds all of them and redraws the observation noise alone; a one-sided formula holds the terms it names and redraws the others.

A smooth, gp() or hsgp() term is a term of the formula, not a group effect, and it is never redrawn, under any re_formula: a smooth indexed by a grouping factor (s(g, bs = "re"), s(x, g, bs = "fs"), a t2() with an re margin) included. Every replicate holds it at its fitted coefficients and the refit estimates it again, so the spread is the uncertainty of the fitted curve.

Value

A frmtmb_boot object: t0 (FUN at the original fit), t (nsim x length(t0) matrix), and converged. confint() gives percentile intervals.

Details

There is no standard bootstrap generic to implement (boot::boot and lme4::bootMer are plain functions), hence the frm_ prefix.

Ordinal thresholds

Every refit reuses the assembled design, so the number of thresholds is the fitted model's and FUN receives the same parameters in the same order in every replicate. A replicate whose response never reaches the top category keeps that threshold, which the replicate's own data no longer identify: its estimate is large, and the replicate's log-likelihood is the one a fit with thres(K) pinned on that replicate gives, to optimizer tolerance. With thres(gr = ) the per-level counts are kept the same way.

simulate() followed by frm(), written out by hand, is not a refit: frm() counts the thresholds of whatever data it is given, so a replicate that never reaches the top category fits one threshold fewer, at the same maximized log-likelihood. Write thres(K) there to pin it.

Examples

set.seed(3)
dd <- data.frame(x = rnorm(80), g = factor(rep(1:8, 10)))
dd$y <- rnorm(80, 1 + 0.5 * dd$x + rnorm(8, 0, 0.5)[dd$g], 1)
fit <- frm(bf(y ~ x + (1 | g)) + gaussian(), data = dd)
bs <- frm_bootstrap(fit, nsim = 20, seed = 1)
bs
#> Parametric bootstrap: 20 refits, 0 failed or not converged
#> 
#>                   estimate      bias       se       lwr     upr
#> (Intercept)       1.235000 -0.017843 0.137530  0.943890 1.40490
#> x                 0.570230  0.028611 0.126600  0.423170 0.88196
#> sigma_(Intercept) 0.055309  0.017652 0.086742 -0.071832 0.24671
confint(bs)
#>                           lwr       upr       est
#> (Intercept)        0.94389449 1.4049436 1.2349783
#> x                  0.42317126 0.8819600 0.5702318
#> sigma_(Intercept) -0.07183161 0.2467122 0.0553095