Beta binomial regression rejecting initial values

Interfaces brms
,
1

I wanted to run a beta binomial regression, but sampling won’t start:

Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts.

Even with a minimal example it does not work. Providing initial values for Intercept and Intercept_phi did not resolve the issue, so I’m lost at what’s going wrong here. Did I specify something incorrectly?

library(brms)
N = 100
n.trial = sample(10:20, size=N, replace=T)
df = data.frame(n.success = rbinom(N, size=n.trial, prob=runif(N,0.4,0.6)),
                n.trial = n.trial)

fit.1 = brm(bf(n.success | trials(n.trial) ~ 1,
               phi ~ 1),
            family = beta_binomial(link="logit", link_phi="log"),
            data = df)
  • Operating System: Linux
  • brms Version: brms_2.22.0

Many thanks,
Benjamin

2

As a place to start, within brm() try setting init = 0, or something like init_r = 0.2.

3

That did not help.

I think found an error in the generated Stan code:

// generated with brms 2.22.0
functions {
}
data {
  int<lower=1> N;  // total number of observations
  array[N] int Y;  // response variable
  array[N] int trials;  // number of trials
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  real Intercept;  // temporary intercept for centered predictors
  real Intercept_phi;  // temporary intercept for centered predictors
}
transformed parameters {
  real lprior = 0;  // prior contributions to the log posterior
  lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
  lprior += student_t_lpdf(Intercept_phi | 3, 0, 2.5);
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = rep_vector(0.0, N);
    // initialize linear predictor term
    vector[N] phi = rep_vector(0.0, N);
    mu += Intercept;
    phi += Intercept_phi;
    mu = inv_logit(mu);
    phi = exp(phi);
    for (n in 1:N) {
      target += beta_binomial_lpmf(Y[n] | trials, mu[n] * phi[n], (1 - mu[n]) * phi[n]);
    }
  }
  // priors including constants
  target += lprior;
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
  // actual population-level intercept
  real b_phi_Intercept = Intercept_phi;
}

In the likelihood function, the variable trials is missing an index [n], it should read

    for (n in 1:N) {
      target += beta_binomial_lpmf(Y[n] | trials[n], mu[n] * phi[n], (1 - mu[n]) * phi[n]);
    }
1 Like
4

Hi Benjamin!

I think reinstalling from github should take care of that:

// generated with brms 2.22.8
functions {
}
data {
  int<lower=1> N;  // total number of observations
  array[N] int Y;  // response variable
  array[N] int trials;  // number of trials
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  real Intercept;  // temporary intercept for centered predictors
  real Intercept_phi;  // temporary intercept for centered predictors
}
transformed parameters {
  real lprior = 0;  // prior contributions to the log posterior
  lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
  lprior += student_t_lpdf(Intercept_phi | 3, 0, 2.5);
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = rep_vector(0.0, N);
    // initialize linear predictor term
    vector[N] phi = rep_vector(0.0, N);
    mu += Intercept;
    phi += Intercept_phi;
    mu = inv_logit(mu);
    phi = exp(phi);
    target += beta_binomial_lpmf(Y | trials, mu .* phi, (1 - mu) .* phi);
  }
  // priors including constants
  target += lprior;
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
  // actual population-level intercept
  real b_phi_Intercept = Intercept_phi;
}

I ran your example without problems.

4 Likes
5

Surprisingly, no-one seems to have used varying trials numbers before this. This bug happened to be fixed in November when switching from for loop to vectorized form as discussed at Use vectorized beta_binomial_lpmf · Issue #1703 · paul-buerkner/brms · GitHub

3 Likes
6

Thank you, will use the newest brms version

1 Like