Good afternoon! I have lesion proportion data in which I have the number of lesions in different regions/levels within the same person and the total number identified. I am interested in knowing whether lesion proportions differ between regions/levels. Therefore, I am attempting to use the custom beta binomial family (Define Custom Response Distributions with brms) with the following code:

lesion_proportions_total_level_complete<-read.csv(“Lesion Proportions Total Level Complete.csv”,header=T)

as.factor(lesion_proportions_total_level_complete$ID)->lesion_proportions_total_level_complete$ID

as.factor(lesion_proportions_total_level_complete$Level)->lesion_proportions_total_level_complete$Level

beta_binomial2 ← custom_family(

“beta_binomial2”, dpars = c(“mu”, “phi”),

links = c(“logit”, “log”), lb = c(NA, 0),

type = “int”, vars = “vint1[n]”

)

stan_funs ← "

real beta_binomial2_lpmf(int y, real mu, real phi, int T) {

return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);

}

int beta_binomial2_rng(real mu, real phi, int T) {

return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);

}

"

stanvars ← stanvar(scode = stan_funs, block = “functions”)

model.lesion.proportions.total.level.complete<-brm(formula=LesionNumber | vint(Total) ~ Level + (1 | ID), data=lesion_proportions_total_level_complete, family=beta_binomial2, stanvars = stanvars, warmup=2000,iter=4000,seed=34)

The model runs and I would then like to use emmeans and bayestestR to extract marginal means and perform contrasts of pairwise comparisons with the following code (as I’ve done with built-in brms families in the past):

emmeans.model.lesion.proportions.total.level.complete<- emmeans(model.lesion.proportions.total.level.complete, ~ Level, epred=TRUE, re_formula=NULL, nesting=NULL)

emmeans.model.lesion.proportions.total.level.complete.summary<- summary(emmeans.model.lesion.proportions.total.level.complete, point.est = mean)

contrasts.model.lesion.proportions.total.level.complete<- contrast(emmeans.model.lesion.proportions.total.level.complete, “pairwise”)

tests.model.lesion.proportions.total.level.complete<-describe_posterior(contrasts.model.lesion.proportions.total.level.complete, centrality=“all”, test=c(“pd”,“ps”,“rope”,“p_map”,“equivalence_test”))

However, I receive the following error after my first emmeans call:

Error: ‘vint’ requires whole numbers as input.

Would you be able to provide any insight into how I can make the model from the custom family compatible with emmeans and bayestestR?

Thanks in advance for your help!