Okay, got it (I think).
rvar is a little bit stricter than manipulating draws directly, and this is deliberate — it tries to make things that are “dangerous” require breaking out of the abstraction. In your case, if I’m understanding correctly, you want to basically create mixtures of the posterior distributions of route_difficulty[...]. I’m not entirely sure what those mixtures correspond to in your model. I can tell you how to make them, but I’d double check they actually are what you want.
As an example, here is an rvar vector of length 10:
library(posterior)
x = rvar_rng(rnorm, 10, 1:10)
x
## rvar<4000>[10] mean ± sd:
## [1] 0.99 ± 1.01 1.97 ± 1.00 3.00 ± 1.01 4.00 ± 1.00 4.99 ± 1.01
## [6] 6.01 ± 1.01 7.00 ± 1.01 8.00 ± 0.97 9.01 ± 1.00 9.98 ± 1.00
We can overplot all those distributions (kinda illegible)
library(ggdist)
library(ggplot2)
data.frame(x) |>
ggplot(aes(xdist = x)) +
stat_halfeye(alpha = 0.25)
You basically want an operation that merges all k draws in an rvar of length n returning a new rvar of length 1 with k*n draws. To do that, you can extract the draws from an rvar with draws_of(), remove the dimension information with as.vector(), and create a new rvar() from the result:
x_mix = rvar(as.vector(draws_of(x)))
x_mix
## rvar<40000>[1] mean ± sd:
## [1] 5.5 ± 3
Note it now has 40,000 draws, not 4,000. It is a mixture of all the above distributions:
data.frame(x_mix) |>
ggplot(aes(xdist = x_mix)) +
stat_halfeye()
If you want to do this inside a data frame, it looks something like this:
library(dplyr)
draws = tibble(
g = rep(c("a", "b"), 5),
x = rvar_rng(rnorm, 10, 1:10)
)
draws
## # A tibble: 10 × 2
## g x
## <chr> <rvar>
## 1 a 0.97 ± 1.01
## 2 b 2.01 ± 1.02
## 3 a 3.00 ± 1.01
## 4 b 4.01 ± 0.99
## 5 a 4.98 ± 1.01
## 6 b 5.99 ± 0.99
## 7 a 7.00 ± 0.99
## 8 b 7.97 ± 1.01
## 9 a 9.00 ± 1.01
## 10 b 10.00 ± 1.01
draws_by_group = draws |>
group_by(g) |>
summarise(x = rvar(as.vector(draws_of(x))))
draws_by_group
## # A tibble: 2 × 2
## g x
## <chr> <rvar>
## 1 a 5 ± 3
## 2 b 6 ± 3
draws_by_group |>
ggplot(aes(y = g, xdist = x)) +
stat_halfeye()
