Dirichlet-multinomial custom family in brms

Interfaces brms
,
1

Hi all

Below I have produced a small simulated data set and used brms to fit a multinomial without any problems. But I am doing something wrong trying to get a custom_family working for the Dirichlet-multinomial (DM). I found the a closed-form PMF for the DM here (Transforming a multinomial model into a dirichlet-multinomial - #2 by nhuurre) and I have been fiddling around with the variables in custom_family and read up as much as I can but just can’t find what I need. Any pointers would be much appreciated.

library(tidyverse)
library(brms)

# Simulate some data

n <- 100
rdata <- t(rmultinom(n = n, size = 5, prob = 1:10))
data <- data.frame(x = factor(sample.int(n = 5, size = n, replace = TRUE)), 
                   size = rowSums(rdata))
data$Y <- rdata

# Multinomial model - this works OK

fit1 <- brm(formula = Y | trials(size) ~ x, data = data, 
            iter = 10, refresh = 1, family = multinomial())

# Dirichlet-multinomial model - cannot get to run

dirichlet_multinomial <- custom_family(
  name = "dirichlet_multinomial", 
  dpars = c("mu", "phi"),
  # multi_dpars = "mu",
  links = "logit", 
  type = "int", 
  lb = c(NA, 0),
  # vars = "vint1",
  # vars = "trials[n]",
  loop = TRUE
)

stan_density <- "
  real dirichlet_multinomial_lpmf(int[] y, vector mu, real phi) {
    real sum_mu = sum(mu);
    real output = lgamma(sum_mu) - lgamma(sum(y) + sum_mu) +
           // + lgamma(sum(y) + 1) - sum(lgamma(to_vector(y) + 1) // constant, may omit
           sum(lgamma(to_vector(y) + mu)) - sum(lgamma(mu));
    return phi * output;
  }
"

stanvars <- stanvar(scode = stan_density, block = "functions")

fit2 <- brm(Y | trials(size) ~ x, data = data, iter = 10, refresh = 1, 
            family = dirichlet_multinomial, stanvars = stanvars)

Thanks in advance - and many thanks for such an amazing R package, brms has become indispensable to me over the past few years!

  • Operating System: Pop!_OS 22.04 LTS (linux)
  • brms Version: 2.17.0
1 Like
2

Hi @quantifish not sure if you’re still pursuing this. I also couldn’t get your code to run. Have you tried to ditto from here? CNVRG/dm.stan at master · JHarrisonEcoEvo/CNVRG · GitHub

I might take a closer look and let’s see if we can get brms to fit this together…

3

Hi @hrlai - I ended up just using the multinomial in the end. But I would be keen to get this up and running for future use. I took a look at the GitHub repo that you sent me, it all makes sense but I really would like to get this working from within brms and I am sure that it is just an issue with the way I am trying to implement a custom family within the brms framework…

4

Hi @quantifish I still couldn’t get it work, but some thoughts:

  • in dirichlet_multinomial the links argument may need to be c("logit", "log") to include the link for phi too
  • should int[] Y be int y because loop = TRUE…?
  • do we need to unmute vars = "trials[n]" to define trial size, and if so…
  • …I don’t know whether we need something like int N in the dirichlet_multinomial_lpmf. If we add N I don’t know how this goes into the PMF
5

Hi @quantifish @hrlai,

I stumbled across this thread whilst trying to create a Dirichlet-multinomial custom family in brms.

As I understand it, custom_family does not currently support families with categorical-like responses (creating a custom_family with a categorical-like response (e.g., Dirichlet-multinomial) · Issue #1728 · paul-buerkner/brms · GitHub).

A native Dirichlet-multinomial family is now available on the master branch of brms on GitHub though, in case this is still of interest

2 Likes