I have a dataset of $m$ individuals. For each individual $m$ I have $n_m$ (binomial ) observations with $s_m$ corresponding to the number of 'successes'.

I use this data to fit a beta-binomial Bayesian model (with a conjugate prior) and update it likewise. This way I can estimate $p_m$ (probability of success) for a (possible new) individual given the observations (if existing) on that (new) individual.

However, on many individuals in my dataset I have only very few observations. These are uncertain records in my dataset. However, there is a relation between the number of observations and $p_m$. So, I should somehow also take into account these 'uncertain individuals' in my dataset. In other words: individuals with a low number of observations tend to have lower $p_m$'s (which I have proven to be true for my specific dataset). If I cut those individuals out (in order to get a more certain dataset) I am getting biased outcomes.

I am thinking about using either Beta-Binomial regression, or fitting in a Poisson-Gamma model as well. However, $n_m$ (number of observations on individual $m$) can grow towards infinity, and, well, I don't know if this is the best way to go? (Overdispersion?)

I don't have much experience with these kinds of models, so I would like to ask what could be a good way forward (Beta-Binomial regression?), and if someone could explain this using my example in simple terms? Thank you!

  • 1
    $\begingroup$ I would go with the beta binomial $\endgroup$ – Demetri Pananos Jan 14 '20 at 16:13
  • $\begingroup$ Thank you Demetri! However, I am still a bit in the dark about how to actually implement the beta binomial regression for my problem (trying to account for the uncertainty in my dataset). Would you mind writing an answer for my question in which you could explain how to implement this model, beta-binomial with beta-binomial regression, for my problem? I know how to implement a simple Bayesian conjugate prior for a beta-binomial distribution (beta prior, binomial likelihood) $\endgroup$ – Damiaan Reijnaers Jan 14 '20 at 18:36
  • $\begingroup$ What software are you using to perform the beta binomial regression? $\endgroup$ – Demetri Pananos Jan 14 '20 at 19:08
  • 1
    $\begingroup$ Here is a tutorial I wrote for the beta binomial model in pymc3. docs.pymc.io/notebooks/GLM-hierarchical-binominal-model.html $\endgroup$ – Demetri Pananos Jan 14 '20 at 19:50
  • 1
    $\begingroup$ Ok. Let me code something up. $\endgroup$ – Demetri Pananos Jan 14 '20 at 20:25

OK, so this isn't the greatest implementation, but I think it serves our purpose.

The beta binomial regression makes the assumption that the data are generated from a binomial distribution

$$ y_i\vert p_i, \kappa, a, b \sim \operatorname{Binomial}(p_i, n) $$

Each person's probability of success is assumed to be iid Beta

$$ p_i \vert \kappa, a, b \sim \operatorname{Beta}(\mu\kappa, (1- \mu)\kappa)$$

Here, $\mu= \operatorname{inv-logit}(a + bn)$, so we are modelling the log odds of the mean of our Beta as a function of the number of times we observe each person. Now, we are free to put priors on $a$ an $b$.

Here is an implementation of this model in python.

import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt
from scipy.special import expit
import arviz as az
import matplotlib.pyplot as plt

#make data

N = 100
ids = np.arange(N)
n = np.random.randint(low=1, high=100, size = N)

true_betas = np.array([-2,0.05])
X = np.c_[np.ones_like(n),n]
true_mean = expit(X@true_betas)
true_kappa = 50

true_p = np.random.beta(true_mean*true_kappa, (1-true_mean)*true_kappa)
y = np.random.binomial(n=n, p = true_p)

with pm.Model() as model:

    intercept = pm.Normal('intercept',mu=0, sd=1)

    slope = pm.Normal('slope',mu=0, sd=1)

    mean = pm.math.invlogit(intercept + n*slope)

    ps = pm.Beta('ps', mean*50, (1-mean)*50, shape=N)

    Y = pm.Binomial('Y',n = n, p=ps, observed = y)

    trace = pm.sample(target_accept=0.9)

    data = az.from_pymc3(trace)

In the model, I haven't put a prior on $\kappa$ but you should. You also should use better priors, I just chose these for exposition.

Model recovers the predicted probabilities quite well (see below).

enter image description here

This is far from a complete bayesian analysis, but I think it should help you get started. Remember to choose priors justifiably and to check your model over and over.


What I do not understand: what do the second and third parameters in your pm.Beta() call mean?

This is just a different parameterization of the Beta distribution. $\mu$ is the mean of the beta and $\kappa$ (which I have set to 50 for illustrative purposes) is a measure of dispersion. When $\kappa$ is large, then the Beta is less variable.

  • $\begingroup$ Thank you for this thorough explanation. I am starting to understand it. What I do not understand: what do the second and third parameters in your pm.Beta() call mean? And why are they multiplied with 50? What exactly do the variables (especially kappa?) represent in your model? Sorry for the many questions! $\endgroup$ – Damiaan Reijnaers Jan 22 '20 at 15:12
  • 1
    $\begingroup$ @DamiaanReijnaers See my edit $\endgroup$ – Demetri Pananos Jan 22 '20 at 16:53
  • $\begingroup$ I'm starting to understand, but doing something wrong. What I do is fitting a Beta dist to the probabilities of individuals with n>=100 (by minimizing the neg-log-likelihood of the Beta's pdf), not taking n into account (my alpha=3.76;beta=10.31). I use my estimated value for Beta as the value for Kappa. I estimate the slope and intercept and find a new Beta distribution for each possible value of n. Is that correct? I think I should use my alpha as well? The results seem to be very wrong. Should I do something with 'ps'?Am I now correctly 'putting a prior on kappa' by using Beta from my fit? $\endgroup$ – Damiaan Reijnaers Jan 26 '20 at 0:44
  • 1
    $\begingroup$ @DamiaanReijnaers That's am Empirical Bayes technique. I'm familiar with David's work but not familiar enough with the technique to be of mite use beyond this point. $\endgroup$ – Demetri Pananos Jan 26 '20 at 21:22
  • 1
    $\begingroup$ @DamiaanReijnaers Kappa is a postitive parameter, so you would have to model kappa like $\log(\kappa) = \beta_0 + \beta_1x_1 + \beta_2x_2$ $\endgroup$ – Demetri Pananos Apr 14 '20 at 19:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.