Bmm: R package for easy and flexible Bayesian Measurement Modeling (v1.0.0)

Publicity
, , ,
3

here’s the brms code that adds a custom link function to a binomial model:

library(brms)

# alldata is a dataframe of per-category binomial outcomes: number of tests, number of positive tests, categorical covariates sex, eth, age, time

alldata$spec = 0.95
alldata$sens = 0.7

hyper_data <- list(
  intercept_prior_mean = 0,
  intercept_prior_scale = 5,
  coef_prior_scale = 2.5
)

# define a *vectorized* custom family (no loop over observations)
binomial_sens_spec_vec <- custom_family(
"binomial_sens_spec", dpars = c("mu"),
links = c("logit"), lb = c(NA),
type="int", vars= c("trials", "vreal1", "vreal2"), loop = FALSE
)

# define the corresponding Stan density function
stan_density_vec <- "
real binomial_sens_spec_lpmf(int[] y, vector mu, int[] N, real[] vreal1, real[] vreal2) {
  return binomial_lpmf(y | N, mu * vreal1[1] + (1 - mu) * (1 - vreal2[1]));
}
"

stanvars_vec <- stanvar(scode = stan_density_vec, block = "functions")

priors <- c(
  set_prior(paste("normal(", hyper_data$intercept_prior_mean, ",", hyper_data$intercept_prior_scale, ")"), class = "Intercept"),
  set_prior(paste("normal(0,", hyper_data$coef_prior_scale, ")"), class = "b", coef = "sex"),
  set_prior(paste("normal(0,", hyper_data$coef_prior_scale, ")"), class = "sd", group = "eth"),
  set_prior(paste("normal(0,", hyper_data$coef_prior_scale, ")"), class = "sd", group = "age"),
  set_prior(paste("normal(0,", hyper_data$coef_prior_scale, ")"), class = "sd", group = "time")
)


fit = brm(cases | trials(tests) + vreal(sens, spec) ~ sex + (1 | eth) + (1 | age) + (1 | time),
          data = alldata,
          family = binomial_sens_spec_vec,
          prior = priors,
          stanvars = stanvars_vec)

More code was required to do proj_pred checks. The project team included @jonah, and it took a lot of digging on both our parts to get things running.

3 Likes