Moment matching for efficient approximate leave-one-out cross-validation
Source:R/loo_moment_match.R
loo_moment_match.brmsfit.RdMoment matching for efficient approximate leave-one-out cross-validation
(LOO-CV). See loo_moment_match
for more details.
Usage
# S3 method for class 'brmsfit'
loo_moment_match(
x,
loo = NULL,
k_threshold = 0.7,
newdata = NULL,
resp = NULL,
check = TRUE,
recompile = FALSE,
...
)
# S3 method for class 'loo'
loo_moment_match(x, fit, ...)Arguments
- x
An R object of class
brmsfitorloodepending on the method.- loo
An R object of class
loo. IfNULL, brms will try to extract a precomputedlooobject from the fitted model, added there viaadd_criterion.- k_threshold
The Pareto \(k\) threshold for which observations moment matching is applied. Defaults to
0.7. Seepareto_k_idsfor more details.- newdata
An optional data.frame for which to evaluate predictions. If
NULL(default), the original data of the model is used.NAvalues within factors (excluding grouping variables) are interpreted as if all dummy variables of this factor are zero. This allows, for instance, to make predictions of the grand mean when using sum coding.NAvalues within grouping variables are treated as a new level.- resp
Optional names of response variables. If specified, predictions are performed only for the specified response variables.
- check
Logical; If
TRUE(the default), some checks check are performed if thelooobject was generated from thebrmsfitobject passed to argumentfit.- recompile
Logical, indicating whether the Stan model should be recompiled. This may be necessary if you are running moment matching on another machine than the one used to fit the model. No recompilation is done by default.
- ...
Further arguments passed to the underlying methods. Additional arguments initially passed to
loo, for example,newdataorrespneed to be passed again toloo_moment_matchin order for the latter to work correctly.- fit
An R object of class
brmsfit.
Details
The moment matching algorithm requires draws of all variables
defined in Stan's parameters block to be saved. Otherwise
loo_moment_match cannot be computed. Thus, please set
save_pars = save_pars(all = TRUE) in the call to brm,
if you are planning to apply loo_moment_match to your models.
References
Paananen, T., Piironen, J., Buerkner, P.-C., Vehtari, A. (2021). Implicitly Adaptive Importance Sampling. Statistics and Computing.
Examples
# \dontrun{
fit1 <- brm(count ~ zAge + zBase * Trt + (1|patient),
data = epilepsy, family = poisson(),
save_pars = save_pars(all = TRUE))
#> Compiling Stan program...
#> Start sampling
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 3.4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 2.176 seconds (Warm-up)
#> Chain 1: 1.564 seconds (Sampling)
#> Chain 1: 3.74 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 2.7e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 2.09 seconds (Warm-up)
#> Chain 2: 1.52 seconds (Sampling)
#> Chain 2: 3.61 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 2.7e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 2.096 seconds (Warm-up)
#> Chain 3: 1.498 seconds (Sampling)
#> Chain 3: 3.594 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 2.8e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 2.219 seconds (Warm-up)
#> Chain 4: 1.662 seconds (Sampling)
#> Chain 4: 3.881 seconds (Total)
#> Chain 4:
# throws warning about some pareto k estimates being too high
(loo1 <- loo(fit1))
#> Warning: Found 9 observations with a pareto_k > 0.7 in model 'fit1'. We recommend to set 'moment_match = TRUE' in order to perform moment matching for problematic observations.
#>
#> Computed from 4000 by 236 log-likelihood matrix.
#>
#> Estimate SE
#> elpd_loo -672.1 36.6
#> p_loo 94.6 14.2
#> looic 1344.2 73.3
#> ------
#> MCSE of elpd_loo is NA.
#> MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
#>
#> Pareto k diagnostic values:
#> Count Pct. Min. ESS
#> (-Inf, 0.7] (good) 227 96.2% 335
#> (0.7, 1] (bad) 8 3.4% <NA>
#> (1, Inf) (very bad) 1 0.4% <NA>
#> See help('pareto-k-diagnostic') for details.
# no more warnings after moment matching
(mmloo1 <- loo_moment_match(fit1, loo = loo1))
#> Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
#>
#> Computed from 4000 by 236 log-likelihood matrix.
#>
#> Estimate SE
#> elpd_loo -672.0 36.6
#> p_loo 94.5 14.2
#> looic 1344.0 73.3
#> ------
#> MCSE of elpd_loo is NA.
#> MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
#>
#> Pareto k diagnostic values:
#> Count Pct. Min. ESS
#> (-Inf, 0.7] (good) 234 99.2% 270
#> (0.7, 1] (bad) 1 0.4% <NA>
#> (1, Inf) (very bad) 1 0.4% <NA>
#> See help('pareto-k-diagnostic') for details.
# }