Help with model reparametrization

1

Dear all,
I am trying to setup a Bayesian metanalysis of report rates of events. My target estimate is the probability that participants will report X in each study. Every participant provides ~20-40 response, where they can say X or non-X. Therefore, each study provides a mean rate X across its subjects. I am trying to setup a Bayesian model to estimate the mean rate.

So far I have reasoned the following.

  1. The distributions I want to work with have to be somehow constrained between 0-1, reflecting the probability to report X.
  2. It is better to work with logs instead of raw probabilities, given that my rates are typically close to 0 - .1.
  3. I will approximate the log standard error as SE / p (delta method).

Here are some simulated data and the model.

import numpy as np
import pandas as pd
import pymc as pm
import arviz as az

np.random.seed(123)

study_means = np.exp(np.random.normal(np.log(.1),.1,20))
study_se = np.random.uniform(0.02, 0.05, 20)

df = pd.DataFrame({
    "mean_rate":study_means,
    "se":study_se,
})

df["log_rate"] = np.log(df["mean_rate"])
df["log_se"] = df["se"] / df["mean_rate"]

# Naive random-effects model
coords = {"random_effects": df.index.values}
          
with pm.Model(coords=coords) as base_model:

    # Set model data
    log_rate = pm.Data("log_rate", df["log_rate"].values, dims="random_effects")
    log_se = pm.Data("log_se", df["log_se"].values, dims="random_effects")
    
    # Priors
    theta = pm.TruncatedNormal("theta", mu=-.7, sigma=1, upper=0)
    tau = pm.HalfCauchy("tau", beta=.5)

    # Transformations
    theta_k = pm.TruncatedNormal("theta_k", mu=theta, sigma=tau, upper=0, dims="random_effects") 

    # Likelihood
    obs = pm.TruncatedNormal("obs", mu = theta_k, sigma = log_se, upper=0, observed = log_rate, dims="random_effects")

    # Sample prior+posterior
    prior_ppc = pm.sample_prior_predictive(500)
    idata = pm.sample(5000,tune=5000,target_accept=0.999)
    posterior_ppc = pm.sample_posterior_predictive(trace=idata,var_names=["theta","tau"])

I have solved the range issue by forcing the normal distribution to have an upper bound of 0, which translates a probability of 1. I have centered the prior around -.5, which roughly centered the distribution at p=.5. The HalfCauchy(0,.5) is standard in modelling tau.

Depending on the simulated dataset (varying mean and sigma of study means), I have a few divergences.

Is there a way to reparametrize this model, or have a different model specification to make it easier to sample from the posterior? Increasing target_accept has lowered the number of divergences so far.

2

One nice transformation is the inverse logit. It will take numbers on the real line and bound them between 0 and 1. One advantage is you can parameterize the random effects in a more natural way and not have to worry about truncating everything just so. It also keeps the data simpler - no need to log transform the data

The other suggestion is to look into whether these random effects exhibit a funnel geometry. My favourite resource on this topic is just this video https://www.youtube.com/watch?v=n2aJYtuGu54&list=PLDcUM9US4XdMROZ57-OIRtIK0aOynbgZN&index=15 (roughly minute 38-45). If they do, you might benefit from a non-centered parameterization. Non-centered parameterizations tend to work better when the data is uninformative about the location of each random effect. Given that you have one observation per effect, you are probably in this case.

Last suggestion is to return the target accept back to the default value. It will cause you model to run much slower and can often hide the extent of a problems. Divergences might go down but your ess is probably suffering. It’s a good strategy to start by fixing model parameterization problems first and only adjust target accept as a last resort.

with pm.Model(coords=coords) as base_model:

    # Set model data
    rate = pm.Data("rate", df["rate"].values, dims="random_effects")
    se = pm.Data("se", df["se"].values, dims="random_effects")
    
    # Priors
    theta = pm.Normal("theta", mu=0, sigma=1)
    tau = pm.HalfCauchy("tau", beta=0.5)

    # Transformations
    theta_k_offset = pm.Normal("theta_k_offset", mu=0, sigma=1, dims="random_effects") 
    theta_k = theta_k_offset * tau + theta
    theta_k = pm.math.invlogit(theta_k)

    # Likelihood
    obs = pm.Normal("obs", mu = theta_k, sigma = se, observed = rate, dims="random_effects")
1 Like
3

Absolutely. You get funnels in hierarchical models whenever the prior plus data do not identify the parameter well. Adding more data usually removes the funnel and you want a centered parameterization.

The standard thing to do here is a hierarchical logistic regression, as @daniel-saunders-phil indicated. I wrote a case study on this model that was translated to PyMC here:

This is also the first serious hierarchical model example in Gelman et al.'s Bayesian Data Analysis at the start of Chapter 5, but Gelman throws in so many concepts from math stats in that one model to generate reparameterized priors that I found it very hard to follow as an introduction. There’s a free pdf on the book’s home page:

2 Likes