some custom_family variables are not passed to partial_log_lik_lpmf with cmdstanr backend and threading

[singmann]

Following up on my previous issue #1110, there appears to be a problem when trying to use custom families with additional variables and cmdstanr backend with multithreading. In particular, the additional variables are not passed to partial_log_lik_lpmf.

Here pretty much the same reprex as before:

library("brms")
## Not run: 
## demonstrate how to fit a beta-binomial model
## generate some fake data
phi <- 0.7
n <- 300
z <- rnorm(n, sd = 0.2)
ntrials <- sample(1:10, n, replace = TRUE)
eta <- 1 + z
mu <- exp(eta) / (1 + exp(eta))
a <- mu * phi
b <- (1 - mu) * phi
p <- rbeta(n, a, b)
y <- rbinom(n, ntrials, p)
dat <- data.frame(y, z, ntrials)

beta_binomial2 <- custom_family(
  "beta_binomial2", dpars = c("mu", "phi"),
  links = c("logit", "log"), lb = c(NA, 0),
  type = "int", vars = c("trials[n]", "foo")
)

# define the corresponding Stan density function
stan_funs <- "
  real beta_binomial2_lpmf(int y, real mu, real phi, int N, real foo) {
    return beta_binomial_lpmf(y | N, mu * phi, (1 - mu) * phi);
  }
"
svars <- stanvar(scode = stan_funs, block = "functions")

brm(y | trials(ntrials) ~ z, data = dat, 
    family = beta_binomial2, stanvars =  svars, backend = "cmdstanr", threads = 2)

It fails with

Semantic error in '/tmp/RtmplWXaup/model-7ac84898c50a.stan', line 34, column 68 to column 71:
   -------------------------------------------------
    32:      for (n in 1:N) {
    33:        int nn = n + start - 1;
    34:        ptarget += beta_binomial2_lpmf(Y[nn] | mu[n], phi, trials[n], foo);
                                                                             ^
    35:      }
    36:      return ptarget;
   -------------------------------------------------

Identifier 'foo' not in scope.

And if we look at the code, foo is indeed not passed:

make_stancode(y | trials(ntrials) ~ z, data = dat, 
              family = beta_binomial2, stanvars =  svars, threads = 2)

for which the relevant part from the functions block is:

 real partial_log_lik_lpmf(int[] seq, int start, int end, int[] Y, int[] trials, matrix Xc, vector b, real Intercept, real phi) {

no foo is passed.

[paul-buerkner]

You did not pass foo in any other way as to name it in the custom family. Stan won't be able to find it no matter if you use threading or not. Am Fr., 26. Feb. 2021 um 09:29 Uhr schrieb Henrik Singmann < notifications@github.com>:

…

Following up on my previous issue #1110 <#1110>, there appears to be a problem when trying to use custom families with additional variables and cmdstanr backend with multithreading. In particular, the additional variables are not passed to partial_log_lik_lpmf. Here pretty much the same reprex as before: library("brms")## Not run: ## demonstrate how to fit a beta-binomial model## generate some fake dataphi <- 0.7n <- 300z <- rnorm(n, sd = 0.2)ntrials <- sample(1:10, n, replace = TRUE)eta <- 1 + zmu <- exp(eta) / (1 + exp(eta))a <- mu * phib <- (1 - mu) * phip <- rbeta(n, a, b)y <- rbinom(n, ntrials, p)dat <- data.frame(y, z, ntrials) beta_binomial2 <- custom_family( "beta_binomial2", dpars = c("mu", "phi"), links = c("logit", "log"), lb = c(NA, 0), type = "int", vars = c("trials[n]", "foo") ) # define the corresponding Stan density functionstan_funs <- " real beta_binomial2_lpmf(int y, real mu, real phi, int N, real foo) { return beta_binomial_lpmf(y | N, mu * phi, (1 - mu) * phi); }"svars <- stanvar(scode = stan_funs, block = "functions") brm(y | trials(ntrials) ~ z, data = dat, family = beta_binomial2, stanvars = svars, backend = "cmdstanr", threads = 2) It fails with Semantic error in '/tmp/RtmplWXaup/model-7ac84898c50a.stan', line 34, column 68 to column 71: ------------------------------------------------- 32: for (n in 1:N) { 33: int nn = n + start - 1; 34: ptarget += beta_binomial2_lpmf(Y[nn] | mu[n], phi, trials[n], foo); ^ 35: } 36: return ptarget; ------------------------------------------------- Identifier 'foo' not in scope. And if we look at the code, foo is indeed not passed: make_stancode(y | trials(ntrials) ~ z, data = dat, family = beta_binomial2, stanvars = svars, threads = 2) for which the relevant part from the functions block is: real partial_log_lik_lpmf(int[] seq, int start, int end, int[] Y, int[] trials, matrix Xc, vector b, real Intercept, real phi) { no foo is passed. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1111>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCW2ACW5ERFUK7OHV7CADTTA5LU7ANCNFSM4YIARW5Q> .

[singmann]

Good point, but does not solve the issue. Was just a problem with the reprex (which would not have worked as you correctly noticed).

library("brms")
## Not run: 
## demonstrate how to fit a beta-binomial model
## generate some fake data
phi <- 0.7
n <- 300
z <- rnorm(n, sd = 0.2)
ntrials <- sample(1:10, n, replace = TRUE)
eta <- 1 + z
mu <- exp(eta) / (1 + exp(eta))
a <- mu * phi
b <- (1 - mu) * phi
p <- rbeta(n, a, b)
y <- rbinom(n, ntrials, p)
dat <- data.frame(y, z, ntrials)

beta_binomial2 <- custom_family(
  "beta_binomial2", dpars = c("mu", "phi"),
  links = c("logit", "log"), lb = c(NA, 0),
  type = "int", vars = c("trials[n]", "foo")
)

# define the corresponding Stan density function
stan_funs <- "
  real beta_binomial2_lpmf(int y, real mu, real phi, int N, real foo) {
    return beta_binomial_lpmf(y | N, mu * phi, (1 - mu) * phi);
  }
"
svars <- stanvar(scode = stan_funs, block = "functions") +
  stanvar(scode = "  real foo = 0.5;", block = "tdata")

brm(y | trials(ntrials) ~ z, data = dat, 
    family = beta_binomial2, stanvars =  svars, backend = "cmdstanr", threads = 2)

still fails with:

Semantic error in '/tmp/RtmplWXaup/model-7ac834c49d85.stan', line 34, column 68 to column 71:
   -------------------------------------------------
    32:      for (n in 1:N) {
    33:        int nn = n + start - 1;
    34:        ptarget += beta_binomial2_lpmf(Y[nn] | mu[n], phi, trials[n], foo);
                                                                             ^
    35:      }
    36:      return ptarget;
   -------------------------------------------------

Identifier 'foo' not in scope.

However, the code without threading now works:

brm(y | trials(ntrials) ~ z, data = dat, 
    family = beta_binomial2, stanvars =  svars, backend = "cmdstanr")

Gives as expected:

All 4 chains finished successfully.
Mean chain execution time: 1.0 seconds.
Total execution time: 4.4 seconds.
 Family: beta_binomial2 
  Links: mu = logit; phi = identity 
Formula: y | trials(ntrials) ~ z 
   Data: dat (Number of observations: 300) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.16      0.11     0.95     1.37 1.00     3515     3093
z             0.92      0.48    -0.03     1.89 1.00     3772     3149

Family Specific Parameters: 
    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi     0.79      0.12     0.59     1.06 1.00     3165     2671

Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

[paul-buerkner]

Yes. I see the problem now. I actually have this has a TODO left in my code from some months ago :-D

Let me think about a solution.

[paul-buerkner]

I just pushed a fix to github. You can now specify pll_args in stanvar() to make sure that foo is added to the header of the partial_log_lik functions used in threading:

library("brms")
## Not run: 
## demonstrate how to fit a beta-binomial model
## generate some fake data
phi <- 0.7
n <- 300
z <- rnorm(n, sd = 0.2)
ntrials <- sample(1:10, n, replace = TRUE)
eta <- 1 + z
mu <- exp(eta) / (1 + exp(eta))
a <- mu * phi
b <- (1 - mu) * phi
p <- rbeta(n, a, b)
y <- rbinom(n, ntrials, p)
dat <- data.frame(y, z, ntrials)

beta_binomial2 <- custom_family(
  "beta_binomial2", dpars = c("mu", "phi"),
  links = c("logit", "log"), lb = c(NA, 0),
  type = "int", vars = c("trials[n]", "foo")
)

# define the corresponding Stan density function
stan_funs <- "
  real beta_binomial2_lpmf(int y, real mu, real phi, int N, real foo) {
    return beta_binomial_lpmf(y | N, mu * phi, (1 - mu) * phi);
  }
"
svars <- stanvar(scode = stan_funs, block = "functions") +
  stanvar(scode = "  real foo = 0.5;", block = "tdata",
          pll_args = "real foo")

brm(y | trials(ntrials) ~ z, data = dat, 
    family = beta_binomial2, stanvars =  svars, backend = "cmdstanr")