Estimating prior sd for the parameter p of a beta binomial regression

Questions
1

I am trying to implement model checks for a Bayesian beta binomial model using the workflow suggested by Betancourt (https://betanalpha.github.io/assets/case_studies/principled_bayesian_workflow.html#51_first_iteration:_a_journey_of_a_thousand_models_begins_with_a_single_step).

I am currently struggling computing the prior_sd of the parameter p (which is required to compute the posterior shrinkage). I know how to compute the sd of a beta distribution based on parameters alpha and beta but I don’t know how to do it when alpha and beta are not numbers but distributions. Could anyone help or point me in the right direction? Below is the current implementation of the model (in pymc3).

with pm.Model() as beta_binomial:

# Priors
intercept = pm.Normal('intercept', mu=np.log(0.1), sd=0.05)
beta = pm.Normal('beta', mu=np.log(0.75), sd=0.05)

# Beta parameters
p_mu = pm.Deterministic('p_mu', pm.math.invlogit(intercept + beta * X_predictor))
p_kappa = pm.HalfNormal('p_kappa', 10)
p_alpha = p_mu * p_kappa
p_beta = (1-p_mu) * p_kappa

# Outcome definition
p = pm.Beta('p', alpha=p_alpha, beta=p_beta, shape=X_predictor.shape[0])
Y = pm.Binomial('num_successes', n=num_trials, p=p, observed=num_successes)
2

Hi,
I’m not sure I understand exactly what you want to do, but:

  • If you want to compute the standard deviation during sampling, you need to use theano operators, as the variables are tensors, not scalars.
  • If you want to compute it after sampling, you can just use the trace and do numpy operations on it, as the trace contains scalars.

Hope this helps :vulcan_salute:
PS: Any reason why you’re not using the built-in Beta-Binomial distribution?

3

My parameter of interest is p, I need to make sure that it gets estimated correctly. So I need to compute the posterior shrinkage and check if it is close to 1:

posterior shrinkage = 1 - post_sd**2/prior_sd**2

I am not sure how to compute the prior sd for p. For example, it is easy to compute it for the intercept or beta (it is 0.05 as per definition) but I don’t know how to propagate the prior variance to p.

Sometimes I use a betabinomial and then extract p from the posterior:

with beta_binomial_hierarchical:
    p = pm.Beta('p', alpha=p_alpha, beta=p_beta, shape=predictor.shape[0])
    posterior_trace = pm.sample_posterior_predictive(trace, samples=10000, var_names=['p'])

But even in this case I don’t know how to estimate prior_sd

4

You can also use

pm.sample_prior_predictive()
5

That is what I am doing at the moment, since I am new to Bayesian modelling I was not sure it was a valid option because all the examples I have seen used values taken directly from the prior.

6

I haven’t looked at the blog and don’t have your input data but making some up to show you the method

X_predictor = 1
num_trials = 8
with pm.Model() as model:
    # Priors
    intercept = pm.Normal('intercept', mu=np.log(0.1), sd=0.05)
    beta = pm.Normal('beta', mu=np.log(0.75), sd=0.05)

    # Beta parameters
    p_mu = pm.Deterministic('p_mu', pm.math.invlogit(intercept + beta * X_predictor))
    p_kappa = pm.HalfNormal('p_kappa', 10)
    p_alpha = p_mu * p_kappa
    p_beta = (1-p_mu) * p_kappa

    # Outcome definition
    p = pm.Beta('p', alpha=p_alpha, beta=p_beta)
    Y = pm.Binomial('num_successes', n=num_trials, p=p)
    
    prior = pm.sample_prior_predictive()
    
prior['p'].std() ### equals 0.1133

Is that what you want to achieve?
You can add an observed with proper predictors, then run a trace and posterior_predictive checks to get a posterior distribution for p to compare to prior.

1 Like
7

Yes, that was what I planned to do, just wanted to make sure it made sense. Thanks