Brms custom family with vint and two vreals

Interfaces brms
1

I am tyring to create a custom family in brms using an lpdf that works fine in rstan. When I try to run the code below, I get an error that says “could not find function vreal”. I’m guessing that the problem is in the custom_family function below or in the specification in the brm() model, but I cannot find the solution to it anywhere else online or on this forum.

paretocounts <- custom_family(
  "paretocounts", dpars = c("mu"),
  links = c("identity"),
  lb = -Inf, ub = Inf,
  type = "real", vars = c("vint1[counts]", 
                          "vreal1[xmin]",
                          "vreal2[xmax]")
)

stan_funs <- "
real paretocounts_lpdf(real x, real mu, real xmin, real xmax, real counts){
    if(mu != -1)
    return(counts*(log((mu+1) / ( xmax^(mu+1) - xmin^(mu+1))) + mu*log(x)));
    else
    return(counts*(log(log(xmin) - log(xmax)) + mu*log(x)));
}
"
stanvars <- stanvar(scode = stan_funs, block = "functions")

#### load attached data called dw_short.csv  ####

fit2 <- brm(
  dw | vint(counts) | vreal(xmin, xmax) ~ mat_s + (1|sample), 
  data = dw_short,
  family = paretocounts,
  stanvars = stanvars
)

Data attached

Please also provide the following information in addition to your question:

  • Operating System: Windows 10 Enterprise
  • brms Version: 2.19.0

dw_short.csv (494.4 KB)

2

just don’t use vint and vreal in the vars argument of custom_family(). brms handles that itself.

3

Thanks for the quick reply. How should the addition terms be added to the brm formula? Currently I have them as follows:

fit2 ← brm(
dw | vint(counts) | vreal(xmin, xmax) ~ mat_s + (1|sample),
data = dw_short,
family = paretocounts,
stanvars = stanvars
)

That didn’t work to begin with. I’ve tried it by removing the addition terms and the error message appears to be looking for the addition terms somewhere:

Error in stanc(file = file, model_code = model_code, model_name = model_name, :
0

Semantic error in ‘string’, line 58, column 16 to column 47:

Ill-typed arguments supplied to function ‘paretocounts_lpdf’. Available signatures:
(real, real, real, real, real) => real
Instead supplied arguments of incompatible type: real, real.

4

they should be there but don’t use the vint part of the vars argument. also combine multiple addition terms with a +. please see the brms_customfamilies vignette for more details

5

Thank you. It’s getting closer but still gives an error and I couldn’t find the solution in the vignette. I think I’m confused about how to enter multiple addition terms in the formula. The current error is this:
Error: The following addition terms are invalid:
‘counts’, ‘xmin’, ‘xmax’

Based on this code:

paretocounts ← custom_family(
“paretocounts”, dpars = c(“mu”),
links = c(“identity”),
lb = -Inf, ub = Inf,
type = “real”, vars = c(“counts”,
“xmin”,
“xmax”))

stan_funs ← "
real paretocounts_lpdf(real Y, real mu, real xmin, real xmax, real counts){
if(mu != -1)
return(counts*(log((mu+1) / ( xmax^(mu+1) - xmin^(mu+1))) + mulog(Y)));
else
return(counts(log(log(xmin) - log(xmax)) + mu*log(Y)));
}
"
stanvars ← stanvar(scode = stan_funs, block = “functions”)

dw_short = dw %>% filter(sample <= 4)

fit2 ← brm(
dw | (counts + xmin + xmax) ~ mat_s + (1|sample),
data = dw_short,
family = paretocounts,
stanvars = stanvars
)

But that gives an error that:

6

y | vint() + vreal() ~ …

7

Thank you. If there are two vreals, how can we add them both?

This doesn’t seem to work:

y|vint(name1) + vreal(name2) + vreal(name3) ~ …

I get the error that addition terms can only occur once.

8

please read the docs more carefully. you are asking many questions that can be inferred from them. you do vreal(x1,x2)

9

Thanks. I’ll keep working on it through the docs. I appreciate the help. It is much closer than when I started.

10

Solved.

rparetocounts <- function(n = 300, mu = -1.2, vreal2 = 1, vreal3 = 1000) {
  samples <- numeric(n)
  {
    if(vreal2 <= 0 | vreal2 >= vreal3) stop("Parameters out of bounds in rPLB")
    u <- stats::runif(n)
    if(mu != -1)
    { y <- ( u*vreal3^(mu+1) +  (1-u) * vreal2^(mu+1) ) ^ (1/(mu+1))
    } else
    { y <- vreal3^u * vreal2^(1-u)
    }
    return(y)
  }
  return(samples)
}

dparetocounts <- function(x, mu, vreal2, vreal3) {
  if (vreal2 <= 0 || vreal2 >= vreal3)
    stop("Parameters out of bounds in dPLB")
  
  if (x < vreal2 || x > vreal3)
    return(0)
  
  if (mu != -1) {
    density <- (mu + 1) * (x^(mu+1)) / (vreal3^(mu+1) - vreal2^(mu+1))
  } else {
    density <- x^(-2) / (vreal2 * log(vreal3/vreal2))
  }
  
  density
}


log_lik_paretocounts <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  y <- prep$data$Y[i]
  dparetocounts(x, mu, vreal2, vreal3)
}


posterior_predict_paretocounts <- function(i, prep, ...) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  vreal2 = prep$data$vreal2[i]
  vreal3 = prep$data$vreal3[i]
  rparetocounts(prep$ndraws, mu, vreal2, vreal3)
}


posterior_epred_paretocounts <- function(prep) {
  mu <- prep$dpars$mu
  return(mu)
}

stan_funs <- "
real paretocounts_lpdf(real Y, real mu, real vreal1, real vreal2, real vreal3){
    if(mu != -1)
    return(vreal1*(log((mu+1) / ( vreal3^(mu+1) - vreal2^(mu+1))) + mu*log(Y)));
    else
    return(vreal1*(log(log(vreal2) - log(vreal3)) + mu*log(Y)));
}
"
stanvars <- brms::stanvar(scode = stan_funs, block = "functions")


paretocounts <- function(){brms::custom_family(
  "paretocounts",
  dpars = c("mu"),
  links = c("identity"),
  lb = -Inf,
  ub = Inf,
  type = "real",
  vars = c("vreal1[n]",
           "vreal2[n]",
           "vreal3[n]"),
  posterior_predict = posterior_predict_paretocounts,
  posterior_epred = posterior_epred_paretocounts,
  log_lik = log_lik_paretocounts)
}