stan-dev / stan-dev/bayesplot

ppc_stat_grouped() fails when subset() is used inside mvbf()

Open
#323 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
R
Stars
442
Forks
93
Avg merge
3h 11m
Merged PRs (30d)
1

Description

Summary

  • Thanks for an amazing package!
  • ppc_stat_grouped() works with a multivariate model fit via brms::brm()
  • ppc_stat_grouped() fails with a multivariate model when subset() is used inside the mvbf()
  • Passing pre-filtered newdata is a workaround; other workarounds welcome
  • May apply to ppc_* other than ppc_stat_grouped(); that's all I tested

In the second example below, I get the error message:

Using all posterior draws for ppc type 'stat_grouped' by default.
Error in `validate_group()`:
! length(group) must be equal to the number of observations.
Run `rlang::last_trace()` to see where the error occurred.
> rlang::last_trace(drop = FALSE)
<error/rlang_error>
Error in `validate_group()`:
! length(group) must be equal to the number of observations.
---
Backtrace:
     ▆
  1. ├─bayesplot::pp_check(...)
  2. └─brms:::pp_check.brmsfit(...)
  3.   └─brms::do_call(ppc_fun, ppc_args)
  4.     └─brms:::eval2(call, envir = args, enclos = envir)
  5.       └─base::eval(expr, envir, ...)
  6.         └─base::eval(expr, envir, ...)
  7.           ├─bayesplot (local) .fun(y = .x1, yrep = .x2, group = .x3, stat = .x4)
  8.           │ └─base::eval(ungroup_call("ppc_stat", call), parent.frame())
  9.           │   └─base::eval(ungroup_call("ppc_stat", call), parent.frame())
 10.           └─bayesplot::ppc_stat(...)
 11.             └─bayesplot::ppc_stat_data(...)
 12.               └─bayesplot:::validate_group(group, length(y))
 13.                 └─rlang::abort("length(group) must be equal to the number of observations.")

Reprex

library(tidyverse)
library(brms)
library(bayesplot)

rbernoulli <- function (n, prob = 0.5) {
  sample(c(0L, 1L), size = n, replace = TRUE, prob = c(1 - prob, prob))
}

dat <-
  tibble(
    grp = LETTERS[1:9],
    y0 = map(seq(0.1, 0.9, by = 0.1), ~ rbernoulli(300, prob = .)),
    y1 = map(seq(0.1, 0.9, by = 0.1), ~ rbernoulli(300, prob = . / 2))) %>%
  unnest(
    c(y0, y1))

#' Just a little helper function to simplify the cases below
fit_logistic <- function (...) {
  brm(
    ...,
    data = dat,
    cores = 4, backend = "cmdstanr", # comment out as you see fit
    family = bernoulli())
}

#' Works great
fit <- fit_logistic(mvbf(y0 ~ (1|grp), y1 ~ (1|grp)))
pp_check(fit, resp = "y1", group = "grp", type = "stat_grouped", stat = "mean")

#' Variation: use `subset(...)` inside the `mvbf()`
#' Fails with "length(group) must be equal to the number of observations"
fit_sub <- fit_logistic(mvbf(y0 ~ (1|grp), y1 | subset(y0 == 1) ~ (1|grp)))
pp_check(fit_sub, resp = "y1", group = "grp", type = "stat_grouped", stat = "mean")

#' If we pass in `newdata` that matches the `subset()` criteria --- here, `y0 == 1` --- then success again
newdata <- filter(dat, y0 == 1)
pp_check(fit_sub, newdata = newdata, resp = "y1", group = "grp", type = "stat_grouped", stat = "mean")

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Run the supplied reprex and compare the working and subset cases. Trace the backtrace through pp_check.brmsfit, ppc_stat, ppc_stat_data, and validate_group, focusing on how group and y are built for subsetted multivariate responses. Done means the grouped statistic works without pre-filtered newdata in the reported case.

Written by the indexing model from the issue text.

Assessment

Tech stack
r
Domain
data-visualization
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Stale
Clarity
Mostly clear
Newbie friendliness
42/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.