SBC derived quantities with vector of parameters

General
1

I’d like to include a derived quantity in my SBC run that uses a parameter vector as input. The simple version of the function takes the vector log_lambda as its argument:

tot_abund <- function(log_lambda) {
  sum(rpois(length(log_lambda), exp(log_lambda)))
}

I then declare the derived quantity as

dqs <- derived_quantities(
  tot_abund = tot_abund_pois(log_lambda),
)

After I set up my generator and backend, I run compute_SBC, including dquants = dqs. This gives an error for each fit:

Error in get(name, f_env): object 'log_lambda' not found

It looks like I would need to write a function that takes arguments log_lambda[1], log_lambda[2], etc. and then combine them and do the computation. Is there a better way to do this? Or to compute the derived quantity after the fact and add it to the fit so that the SBC statistics can be computed?

The more complex version of the function depends on the generated data as well. Are these vectors available to derived quantities, or just the posterior parameter samples?

Thank you!

2

I think I managed to construct additional rows for the SBC statistics (sbc$stats) data frame correctly using the following function:

tot_abund_stats <- function(sbc, dataset, backend) {
  # Need the number of SBC replicates to iterate over both the fits and
  # generated variables below
  n_sims <- length(sbc$sims)

  # First calculate posterior samples of the total abundance based on the MCMC
  # draws
  drawmat <- map2(
    sbc$fits,
    dataset$generated,
    \(f, g) {
      as_draws_rvars(f) |>
        pluck("log_lambda") |>
        exp() |>
        # Use `_` as a placeholder argument; only works in recent versions of R
        rvar_rng(rpois, g$N_unm_obs, _) |>
        rvar_sum() |>
        as_draws_matrix() |>
        set_variables("tot_abund")
    }
  )

  # Calculate the rank statistics using `calculate_ranks_draws_matrix`, then
  # construct a data frame (here using `tibble` for convenience) that matches
  # the form of the `sbc$stats` data frame so that `bind_rows` can be used.
  # Iterate over the replicate index so that we can include it as `sim_id` in
  # the stats matrix.
  tot_abund_stats_df <- map(
    seq_len(n_sims),
    \(i) {
      v <- extract_variable(dataset$variables, "tot_abund")[i]
      d <- drawmat[[i]]
      d_thin <- d |>
        thin_draws(SBC_backend_default_thin_ranks(backend))
      # Drop last draw so that number of ranks is even (matches default for other
      # parameters, allowing them to be combined). CAUTION this is not the same
      # logic used in the SBC package but works with the default number of
      # draws and thin rate.
      d_thin2 <- d_thin |> subset_draws(draw = seq_len(ndraws(d_thin) - 1))
      r <- calculate_ranks_draws_matrix(v, d_thin2)
      tibble(
        sim_id = i,
        variable = "tot_abund",
        simulated_value = v,
        rank = r,
        z_score = (simulated_value - mean(d)) / sd(d),
        mean = mean(d),
        sd = sd(d),
        q5 = quantile(d, 0.05),
        median = median(d),
        q95 = quantile(d, 0.95),
        rhat = rhat(d),
        ess_bulk = ess_bulk(d),
        ess_tail = ess_tail(d),
        attributes = "",
        max_rank = attr(r, "max_rank"),
        has_na = anyNA(d),
        all_na = all(is.na(d)),
        all_inf = all(is.infinite(d))
      )
    }
  ) |>
    bind_rows()
}

I did have to work around errors because tot_abund is included in sbc_dataset$variables but not the Stan posterior samples. After using bind_rows to add this to the existing sbc$stats data frame, I am able to generate the diagnostic plots for these quantities.

3

In general, this should work - see e.g. SBC for brms models • SBC and the definition of log_lik_dq_func, could you share a bit more complete example of where it fails for you?

At some point, I’d like to provide more flexibility in building derived quantities, the non-standard evaluation approach is not for everyone anyway.

4

Do you mean you want to use some kind of multivariate “test” rather than a univariate one? I put “test” in quotes here because a lot of the SBC work post Samantha Cook’s original is cagey about formulating it as a hypothesis test.

I haven’t seen this before, but @martinmodrak may know if people have defined this kind of thing for SBC before, and I’d be curious about how they did it for general target densities. My inclination off the cuff would be to define the union of the level sets at the central log density, that is, all the parameter values whose posterior log density is in the central 90% (or some other) interval of log densities. This is close to the idea of a typical set in information theory, but with nominal coverage rather than asymptotic guarantees.