Beta-Binomial regression or Poisson-Gamma model to account for uncertainty in (empricial Bayesian) prior? Explained in simple terms?

Question

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!

Answer 1

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).

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.

EDIT:

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.