loo.Rd
The loo()
methods for arrays, matrices, and functions compute PSISLOO
CV, efficient approximate leaveoneout (LOO) crossvalidation for Bayesian
models using Pareto smoothed importance sampling (PSIS). This is an
implementation of the methods described in Vehtari, Gelman, and Gabry (2017a,
2017b).
The loo_i()
function enables testing loglikelihood
functions for use with the loo.function()
method.
loo(x, ...) # S3 method for array loo( x, ..., r_eff = NULL, save_psis = FALSE, cores = getOption("mc.cores", 1), is_method = c("psis", "tis", "sis") ) # S3 method for matrix loo( x, ..., r_eff = NULL, save_psis = FALSE, cores = getOption("mc.cores", 1), is_method = c("psis", "tis", "sis") ) # S3 method for `function` loo( x, ..., data = NULL, draws = NULL, r_eff = NULL, save_psis = FALSE, cores = getOption("mc.cores", 1), is_method = c("psis", "tis", "sis") ) loo_i( i, llfun, ..., data = NULL, draws = NULL, r_eff = NULL, is_method = "psis" ) is.loo(x) is.psis_loo(x)
x  A loglikelihood array, matrix, or function. The Methods (by class) section, below, has detailed descriptions of how to specify the inputs for each method. 

r_eff  Vector of relative effective sample size estimates for the
likelihood ( 
save_psis  Should the 
cores  The number of cores to use for parallelization. This defaults to
the option

is_method  The importance sampling method to use. The following methods are implemented: 
data, draws, ...  For the 
i  For 
llfun  For 
The loo()
methods return a named list with class
c("psis_loo", "loo")
and components:
estimates
A matrix with two columns (Estimate
, SE
) and four rows
(elpd_loo
, mcse_elpd_loo
, p_loo
, looic
). This
contains point estimates and standard errors of the expected log pointwise
predictive density (elpd_loo
), the Monte Carlo standard error of
elpd_loo
(mcse_elpd_loo
), the effective number of parameters
(p_loo
) and the LOO information criterion looic
(which is
just 2 * elpd_loo
, i.e., converted to deviance scale).
pointwise
A matrix with four columns (and number of rows equal to the number of
observations) containing the pointwise contributions of each of the above
measures (elpd_loo
, mcse_elpd_loo
, p_loo
, looic
).
diagnostics
A named list containing two vectors:
pareto_k
: Estimates of the shape parameter \(k\) of the
generalized Pareto fit to the importance ratios for each leaveoneout
distribution. See the paretokdiagnostic page for details.
n_eff
: PSIS effective sample size estimates.
psis_object
This component will be NULL
unless the save_psis
argument is set to
TRUE
when calling loo()
. In that case psis_object
will be the object
of class "psis"
that is created when the loo()
function calls psis()
internally to do the PSIS procedure.
The loo_i()
function returns a named list with components
pointwise
and diagnostics
. These components have the same
structure as the pointwise
and diagnostics
components of the
object returned by loo()
except they contain results for only a single
observation.
The loo()
function is an S3 generic and methods are provided for
3D pointwise loglikelihood arrays, pointwise loglikelihood matrices, and
loglikelihood functions. The array and matrix methods are the most
convenient, but for models fit to very large datasets the loo.function()
method is more memory efficient and may be preferable.
array
: An \(I\) by \(C\) by \(N\) array, where \(I\)
is the number of MCMC iterations per chain, \(C\) is the number of
chains, and \(N\) is the number of data points.
matrix
: An \(S\) by \(N\) matrix, where \(S\) is the size
of the posterior sample (with all chains merged) and \(N\) is the number
of data points.
function
: A function f()
that takes arguments data_i
and draws
and returns a
vector containing the loglikelihood for a single observation i
evaluated
at each posterior draw. The function should be written such that, for each
observation i
in 1:N
, evaluating
f(data_i = data[i,, drop=FALSE], draws = draws)
results in a vector of length S
(size of posterior sample). The
loglikelihood function can also have additional arguments but data_i
and
draws
are required.
If using the function method then the arguments data
and draws
must also
be specified in the call to loo()
:
data
: A data frame or matrix containing the data (e.g.
observed outcome and predictors) needed to compute the pointwise
loglikelihood. For each observation i
, the i
th row of
data
will be passed to the data_i
argument of the
loglikelihood function.
draws
: An object containing the posterior draws for any
parameters needed to compute the pointwise loglikelihood. Unlike
data
, which is indexed by observation, for each observation the
entire object draws
will be passed to the draws
argument of
the loglikelihood function.
The ...
can be used if your loglikelihood function takes additional
arguments. These arguments are used like the draws
argument in that they
are recycled for each observation.
loo()
methods in a packagePackage developers can define
loo()
methods for fitted models objects. See the example loo.stanfit()
method in the Examples section below for an example of defining a
method that calls loo.array()
. The loo.stanreg()
method in the
rstanarm package is an example of defining a method that calls
loo.function()
.
Vehtari, A., Gelman, A., and Gabry, J. (2017a). Practical Bayesian model evaluation using leaveoneout crossvalidation and WAIC. Statistics and Computing. 27(5), 14131432. doi:10.1007/s1122201696964 (journal version, preprint arXiv:1507.04544).
Vehtari, A., Gelman, A., and Gabry, J. (2017b). Pareto smoothed importance sampling. preprint arXiv:1507.02646
The loo package vignettes for demonstrations.
psis()
for the underlying Pareto Smoothed Importance Sampling (PSIS)
procedure used in the LOOCV approximation.
paretokdiagnostic for convenience functions for looking at diagnostics.
loo_compare()
for model comparison.
### Array and matrix methods (using example objects included with loo package) # Array method LLarr < example_loglik_array() rel_n_eff < relative_eff(exp(LLarr)) loo(LLarr, r_eff = rel_n_eff, cores = 2)#> #> Computed from 1000 by 32 loglikelihood matrix #> #> Estimate SE #> elpd_loo 83.6 4.3 #> p_loo 3.3 1.2 #> looic 167.2 8.6 #>  #> Monte Carlo SE of elpd_loo is 0.1. #> #> All Pareto k estimates are good (k < 0.5). #> See help('paretokdiagnostic') for details.# Matrix method LLmat < example_loglik_matrix() rel_n_eff < relative_eff(exp(LLmat), chain_id = rep(1:2, each = 500)) loo(LLmat, r_eff = rel_n_eff, cores = 2)#> #> Computed from 1000 by 32 loglikelihood matrix #> #> Estimate SE #> elpd_loo 83.6 4.3 #> p_loo 3.3 1.2 #> looic 167.2 8.6 #>  #> Monte Carlo SE of elpd_loo is 0.1. #> #> All Pareto k estimates are good (k < 0.5). #> See help('paretokdiagnostic') for details.# \dontrun{ ### Usage with stanfit objects # see ?extract_log_lik log_lik1 < extract_log_lik(stanfit1, merge_chains = FALSE)#> Error in extract_log_lik(stanfit1, merge_chains = FALSE): object 'stanfit1' not found#> Error in relative_eff(exp(log_lik1)): object 'log_lik1' not foundloo(log_lik1, r_eff = rel_n_eff, cores = 2)#> Error in loo(log_lik1, r_eff = rel_n_eff, cores = 2): object 'log_lik1' not found# } ### Using loglikelihood function instead of array or matrix set.seed(124) # Simulate data and draw from posterior N < 50; K < 10; S < 100; a0 < 3; b0 < 2 p < rbeta(1, a0, b0) y < rbinom(N, size = K, prob = p) a < a0 + sum(y); b < b0 + N * K  sum(y) fake_posterior < as.matrix(rbeta(S, a, b)) dim(fake_posterior) # S x 1#> [1] 100 1#> [1] 50 2llfun < function(data_i, draws) { # each time called internally within loo the arguments will be equal to: # data_i: ith row of fake_data (fake_data[i,, drop=FALSE]) # draws: entire fake_posterior matrix dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE) } # Use the loo_i function to check that llfun works on a single observation # before running on all obs. For example, using the 3rd obs in the data: loo_3 < loo_i(i = 3, llfun = llfun, data = fake_data, draws = fake_posterior, r_eff = NA) print(loo_3$pointwise[, "elpd_loo"])#> elpd_loo #> 1.267103# Use loo.function method (setting r_eff=NA since this posterior not obtained via MCMC) loo_with_fn < loo(llfun, draws = fake_posterior, data = fake_data, r_eff = NA) # If we look at the elpd_loo contribution from the 3rd obs it should be the # same as what we got above with the loo_i function and i=3: print(loo_with_fn$pointwise[3, "elpd_loo"])#> elpd_loo #> 1.267103#> elpd_loo #> 1.267103# Check that the loo.matrix method gives same answer as loo.function method log_lik_matrix < sapply(1:N, function(i) { llfun(data_i = fake_data[i,, drop=FALSE], draws = fake_posterior) }) loo_with_mat < loo(log_lik_matrix, r_eff = NA) all.equal(loo_with_mat$estimates, loo_with_fn$estimates) # should be TRUE!#> [1] TRUE# \dontrun{ ### For package developers: defining loo methods # An example of a possible loo method for 'stanfit' objects (rstan package). # A similar method is planned for a future release of rstan (or is already # released, depending on when you are reading this). In order for users # to be able to call loo(stanfit) instead of loo.stanfit(stanfit) the # NAMESPACE needs to be handled appropriately (roxygen2 and devtools packages # are good for that). # loo.stanfit < function(x, pars = "log_lik", ..., save_psis = FALSE, cores = getOption("mc.cores", 1)) { stopifnot(length(pars) == 1L) LLarray < loo::extract_log_lik(stanfit = x, parameter_name = pars, merge_chains = FALSE) r_eff < loo::relative_eff(x = exp(LLarray), cores = cores) loo::loo.array(LLarray, r_eff = r_eff, cores = cores, save_psis = save_psis) } # }