Check variable option after appending prediction

[avehtari]

A simple binomial model and data, and we compute theta=plogis(b_Intercept)

data_bin <- data.frame(N = c(10), y = c(7))
fit_bin <- brm(y | trials(N) ~ 1, family = binomial(), data = data_bin, refresh = 0)
theta <- as_draws_df(fit_bin) |>
  mutate_variables(theta=plogis(b_Intercept)) |>
  subset_draws(variable='theta')

As prediction appends the draws, the following prints sensitivity results also for b_Intercept

powerscale_sensitivity(fit_bin, prediction = \(x) theta)

Adding option variable='theta'

powerscale_sensitivity(fit_bin, prediction = \(x) theta, variable='theta')

causes an error Error: The following variables are missing in the draws object: {'theta'}

Adding option variable='theta'

powerscale_sensitivity(fit_bin, prediction = \(x) theta, variable='b_Intercept')

causes an error Error in prediction(x, ...) : unused argument (variable = "b_Intercept")

It would be nice to be able to focus only on the quantities of interest.

[avehtari]

This

causes an error Error in prediction(x, ...) : unused argument (variable = "b_Intercept")
is solved by

powerscale_sensitivity(fit_bin, prediction = \(x, ...) theta, variable='b_Intercept')

It would be useful to mention in the doc, that the prediction function should have ...

Or maybe allow prediction (or other option name) to be a draws object

[n-kall]

I think the better workflow is to do all mutations to the draws beforehand (i.e. adding the theta) and then the sensitivity analysis.

Two options that might offer a better workflow:

data_bin <- data.frame(N = c(10), y = c(7))
fit_bin <- brm(y | trials(N) ~ 1, family = binomial(), data = data_bin, refresh = 0, backend = "cmdstanr")

# here we manually add the log_liks into the draws object before power-scaling
fit_bin |>
  as_draws_df() |>
  mutate_variables(theta = plogis(b_Intercept)) |>
  bind_draws(log_lik_draws(fit_bin)) |>
  powerscale_sensitivity(variable = "theta")


# here we specify the fit object in the power-scaling call, which calculates the log_liks
fit_bin |>
  as_draws_df() |>
  mutate_variables(theta = plogis(b_Intercept)) |>
  powerscale_sensitivity(fit = fit_bin, variable = "theta")

[n-kall]

I'll need to work on how to do this kind of workflow with moment-matching though