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!