Specifying nested random effects in "brms"

Question

I am new to using "brms" and am encountering an issue when specifying the model formula. The error occurs in the specification of the random intercept (three variables/terms denoting the nested, hierarchical group structure). The error message I encounter is shown below:

Setting 'rescor' to FALSE by default for this model
Error in str2lang(x) : <text>:1:17: unexpected ')'
1: ~ var4/var3/var1)
                    ^

Could someone please lend a hand in resolving this issue? Any help would be greatly appreciated. My code is shown below:

result <- brm(
  bf(
    mvbind(var9, var10, var15, var16) ~ var2 + var13 + var14 + var19 + var20 + var34 + var46 + var47 + var50 + var51 + var52 + var56 + var57 + var59 + var60 + var62 + var65 + var66 + var89 + var90 + var91 + var92 + var93 + var94 + cor_ar(var2 + 1 | var4/var3/var1),
           zi ~ var2 + var5 + var114 + cor_ar(var2 + 1 | var4/var3/var1)
    ),
  data = y,
  family = zero_inflated_negbinomial(),
  iter = 1000,
  chains = 4,
  future = TRUE)

Other aspects of the model formula that I am unsure of are:

1. specification of autoregressive covariance structure using cor_ar()
2. setting future = TRUE for parallel computation of chains.

Please let me know if these are specified correctly. Thank you!

Answer 1

The beginning of a simulation setup for testing brm syntax. This is not working right now because, for ar(), values of time within groups must be unique. It's hard to construct the nested factor structure and make sure that time values are unique within groups ... or at least not obvious ...

Nevertheless I thought I'd leave this here temporarily in case it's useful.

library(brms)
library(glmmTMB)
library(dplyr)
nn <- 1000
ng <- 50
set.seed(101)
s <- function() factor(sample(1:ng, size = nn, replace = TRUE))
dd <- data.frame(x1 = rnorm(nn), x2 = rnorm(nn),
     x3 = rnorm(nn), f1 = s(), f2 = s(), f3 = s())
dd <- (dd
   |> group_by(f1, f2, f3)
   |> mutate(time = factor(seq(n())))
)

form <- y ~ x1 + x2 + x3 + 
    ar1(0 + time|f1) + ar1(0 + time|f1:f2) + ar1(0 + time|f1:f2:f3)
dd$y <- simulate_new(
    form[-2], ## formula without LHS
    ziformula = ~1,
    newdata = dd,
    family = gaussian,
    newparams = list(beta = 1:4,
            theta = rep(1, 6),
            betad = 1,
            betazi = 1))[[1]]

brm( y ~ x1 + x2 + x3 + 
    ar(time = time, gr = f1),
    family = gaussian,
    data = dd)