How to switch from `dBetaBinom_One()` to `dBetaBinom_v()` and `dBetaBinom_s()`?

[maxdrohde]

I'm trying to create a Nimble model with a beta-binomial distribution using the nimbleEcology package.

I started by trying the dBetaBinom_v() and dBetaBinom_s() functions but I couldn't get the MCMC to run properly because I think these functions expect a vector of observations, rather than a scalar. However, I saw that there used to be a function in this package, dBetaBinom_One(), which I used (by installing an old version of the this package) and the model worked correctly.

The relevant part of my Nimble code looks like this:

for (g in 1:G) {
    p[g]  <- expit(eta[g])
    s1[g] <- phi * p[g]
    s2[g] <- phi * (1 - p[g])

    x[g] ~ dBetaBinom_One(
        N      = x_size,
        shape1 = s1[g],
        shape2 = s2[g]
    )
}

Would you be able to let me know how to write this same model using the current version of the functions, either dBetaBinom_v() and dBetaBinom_s()? It may also be helpful to others who are having this issue and I would be happy to contribute a pull request to the documentation.

Thank you,
Max

[dochvam]

Hi @maxdrohde --

Thanks for the comment. You're right that the dBetaBinom_v and dBetaBinom_s distributions both expect vector-type x, and because of NIMBLE's strict typing that means you can't have scalar x. We removed dBetaBinom_One after doing some reorganization of the dNmixture code and finding that we no longer needed it internally, but it does cover a distinct case from the _s and _v distributions.

In this case, you can achieve what you're after with dBetaBinom_v by moving your call to dBetaBinom out of the for-loop. I believe the following model code is equivalent to what you've written:

  for (g in 1:G) {
    p[g]  <- expit(eta[g])
    s1[g] <- phi * p[g]
    s2[g] <- phi * (1 - p[g])
  }
  
  x[1:G] ~ dBetaBinom_v(N = x_size, shape1 = s1[1:G], shape2 = s2[1:G], len = G)

Here's a full test script you can use to confirm for yourself that these are doing the same thing:

library(nimbleEcology)
# Define dBetaBinom_One locally:
dBetaBinom_One <- nimbleFunction(
  run = function(x = double(0),
                 N = double(0),
                 shape1 = double(0),
                 shape2 = double(0),
                 log = integer(0, default = 0)) {
    
    logprob <- 0
    logprob <- logprob +
      nimBetaFun(a = x + shape1, b = N - x + shape2, log = TRUE) -
      nimBetaFun(a = shape1, b = shape2, log = TRUE) +
      lgamma(N+1) - (lgamma(x+1) + lgamma(N - x + 1))
    if (log) return(logprob)
    return(exp(logprob))
    returnType(double(0))
  }
)
rBetaBinom_One <- nimbleFunction(
  run = function(n = double(0),
                 N = double(0),
                 shape1 = double(0),
                 shape2 = double(0)) {
    p <- rbeta(1, shape1, shape2)
    x <- rbinom(1, N, p)
    return(x)
    returnType(double())
  })

# Your example:
test_code_One <- nimbleCode({
  for (g in 1:G) {
    p[g]  <- expit(eta[g])
    s1[g] <- phi * p[g]
    s2[g] <- phi * (1 - p[g])
    
    x[g] ~ dBetaBinom_One(
      N      = x_size,
      shape1 = s1[g],
      shape2 = s2[g]
    )
  }
})


mod_One <- nimbleModel(
  test_code_One,
  constants = list(G = 10, eta = 1:10/10, phi = 1, x_size = 10),
  data = list(x = rep(2, 10))
)
cmod_One <- compileNimble(mod_One)
lprob_One <- cmod_One$calculate()


# Vectorized version
test_code_v <- nimbleCode({
  for (g in 1:G) {
    p[g]  <- expit(eta[g])
    s1[g] <- phi * p[g]
    s2[g] <- phi * (1 - p[g])
  }
  
  x[1:G] ~ dBetaBinom_v(N = x_size, shape1 = s1[1:G], shape2 = s2[1:G], len = G)
})

mod_v <- nimbleModel(
  test_code_v,
  constants = list(G = 10, eta = 1:10/10, phi = 1, x_size = 10),
  data = list(x = rep(2, 10))
)
cmod_v <- compileNimble(mod_v)
lprob_v <- cmod_v$calculate()


# Confirm that they match
lprob_v
lprob_One

However, there might be similar cases where this solution isn't appropriate--for example, if each observation has its own x_size, or if there's some reason you don't want all elements of x to be calculated at the same time--so we should definitely consider bringing back dBetaBinom_One even though we don't need it internally for dNmixture anymore. For that reason I'm going to leave this issue open.
If you want to use dBetaBinom_One, I recommend defining dBetaBinom_One in your local session as a custom function as in the example above, rather than installing an older version of the package. NIMBLE will register it automatically when you attempt to use it in a model.

[maxdrohde]

Thanks so much for your help @dochvam!

I had tried the solution you proposed with

x[1:G] ~ dBetaBinom_v(N = x_size, shape1 = s1[1:G], shape2 = s2[1:G], len = G)

and it didn't work, I realize now, for the reason you mentioned

if there's some reason you don't want all elements of x to be calculated at the same time

because most of the values of x are missing since I'm using this model for imputation. I will use the local version of the dBetaBinom_One function that you showed.