# Dirichlet-multinomial custom family in brms

**URL:** <https://discourse.mc-stan.org/t/dirichlet-multinomial-custom-family-in-brms/27590>\
**Category:** brms\
**Tags:** dirichlet-multinomial, brms\
**Created:** [May 22, 2022, 4:18am UTC](https://discourse.mc-stan.org/t/dirichlet-multinomial-custom-family-in-brms/27590 "2022-05-22T04:18:03Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![quantifish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/quantifish/32/14329_2.png) [@quantifish](https://discourse.mc-stan.org/u/quantifish)\
**Post date:** [May 22, 2022, 4:18am UTC](https://discourse.mc-stan.org/t/dirichlet-multinomial-custom-family-in-brms/27590/1 "2022-05-22T04:18:03Z")

</div>

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](https://discourse.mc-stan.org/t/transforming-a-multinomial-model-into-a-dirichlet-multinomial/26399/2)) 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.

```nohighlight
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

---

<div class="post-metadata">

**Author:** ![hrlai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hrlai/32/6419_2.png) [@hrlai](https://discourse.mc-stan.org/u/hrlai)\
**Post date:** [May 5, 2023, 5:12am UTC](https://discourse.mc-stan.org/t/dirichlet-multinomial-custom-family-in-brms/27590/2 "2023-05-05T05:12:32Z")

</div>

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](https://github.com/JHarrisonEcoEvo/CNVRG/blob/master/inst/stan/dm.stan)

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

---

<div class="post-metadata">

**Author:** ![quantifish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/quantifish/32/14329_2.png) [@quantifish](https://discourse.mc-stan.org/u/quantifish)\
**Post date:** [May 8, 2023, 11:57pm UTC](https://discourse.mc-stan.org/t/dirichlet-multinomial-custom-family-in-brms/27590/3 "2023-05-08T23:57:01Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![hrlai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hrlai/32/6419_2.png) [@hrlai](https://discourse.mc-stan.org/u/hrlai)\
**Post date:** [August 15, 2023, 9:36pm UTC](https://discourse.mc-stan.org/t/dirichlet-multinomial-custom-family-in-brms/27590/4 "2023-08-15T21:36:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![tom-peatman](https://avatars.discourse-cdn.com/v4/letter/t/87869e/32.png) [@tom-peatman](https://discourse.mc-stan.org/u/tom-peatman)\
**Post date:** [February 3, 2025, 4:42pm UTC](https://discourse.mc-stan.org/t/dirichlet-multinomial-custom-family-in-brms/27590/5 "2025-02-03T16:42:43Z")

</div>

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](https://github.com/paul-buerkner/brms/issues/1728)).

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