The `loo`

methods for arrays, matrices, and functions compute PSIS-LOO
CV, efficient approximate leave-one-out (LOO) cross-validation 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 log-likelihood
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)) # S3 method for matrix loo(x, ..., r_eff = NULL, save_psis = FALSE, cores = getOption("mc.cores", 1)) # S3 method for function loo(x, ..., data = NULL, draws = NULL, r_eff = NULL, save_psis = FALSE, cores = getOption("mc.cores", 1)) loo_i(i, llfun, ..., data = NULL, draws = NULL, r_eff = NULL)

x | A log-likelihood array, matrix, or function. See the |
---|---|

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 |

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 leave-one-out distribution. See the`pareto-k-diagnostic`

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 computing LOO from 3-D pointwise log-likelihood arrays, pointwise
log-likelihood matrices, and log-likelihood functions. The array and matrix
methods are 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 log-likelihood 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 log-likelihood 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 log-likelihood. For each observation`i`

, the`i`

th row of`data`

will be passed to the`data_i`

argument of the log-likelihood function.`draws`

: An object containing the posterior draws for any parameters needed to compute the pointwise log-likelihood. Unlike`data`

, which is indexed by observation, for each observation the entire object`draws`

will be passed to the`draws`

argument of the log-likelihood function.The

`...`

can be used to pass additional arguments to your log-likelihood function. 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 rstanarm is an example of defining a
method that calls `loo.function`

.

Vehtari, A., Gelman, A., and Gabry, J. (2017a). Practical
Bayesian model evaluation using leave-one-out cross-validation and WAIC.
*Statistics and Computing*. 27(5), 1413--1432.
doi:10.1007/s11222-016-9696-4.
(published
version, arXiv preprint).

Vehtari, A., Gelman, A., and Gabry, J. (2017b). Pareto smoothed importance sampling. arXiv preprint: http://arxiv.org/abs/1507.02646/

The loo package vignettes for demonstrations.

`psis`

for the underlying Pareto Smoothed Importance Sampling (PSIS) procedure used in the LOO-CV approximation.pareto-k-diagnostic for convenience functions for looking at diagnostics.

`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 log-likelihood matrix #> #> Estimate SE #> elpd_loo -83.6 4.3 #> p_loo 3.4 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('pareto-k-diagnostic') 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 log-likelihood matrix #> #> Estimate SE #> elpd_loo -83.6 4.3 #> p_loo 3.4 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('pareto-k-diagnostic') for details.# NOT RUN { ### Usage with stanfit objects # see ?extract_log_lik log_lik1 <- extract_log_lik(stanfit1, merge_chains = FALSE) rel_n_eff <- relative_eff(exp(log_lik1)) loo(log_lik1, r_eff = rel_n_eff, cores = 2) # }### Using log-likelihood 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 1fake_data <- data.frame(y,K) dim(fake_data) # N x 2#> [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)#> Warning: Relative effective sample sizes ('r_eff' argument) not specified. PSIS n_eff will not adjusted based on MCMC n_eff.print(loo_3$pointwise[, "elpd_loo"])#> elpd_loo #> -1.267871# Use loo.function method loo_with_fn <- loo(llfun, draws = fake_posterior, data = fake_data)#> Warning: Relative effective sample sizes ('r_eff' argument) not specified. #> For models fit with MCMC, the reported PSIS effective sample sizes and #> MCSE estimates will be over-optimistic.# 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.267871print(loo_3$pointwise[, "elpd_loo"])#> elpd_loo #> -1.267871# 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)#> Warning: Relative effective sample sizes ('r_eff' argument) not specified. #> For models fit with MCMC, the reported PSIS effective sample sizes and #> MCSE estimates will be over-optimistic.all.equal(loo_with_mat$estimates, loo_with_fn$estimates) # should be TRUE!#> [1] TRUE# NOT RUN { ### 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) } # }