# Partial pooling: Understanding sampling efficiency and "scale dependence"

I haven’t been able to fully wrap my head around the “simpson paradox” example yet. For the data generating process in the example it seems to be quite challenging to set up a straightforward (hierarchical) model that has no divergences/low ESS values for certain variables, cf. discussion in Simpsons paradox and mixed models - large number of groups - #16 by aseyboldt

In order to understand this I have tried to simplify the data generation process as follows. I am simply drawing (univariate) random normal variates having a group-specific mean \mu_g and a variance \sigma^2 (the same for all groups for the sake of simplicity). The group-specific mean is random normally distributed with location (mean) \mu and variance \omega^2.

In other words I have an intra-group variability \sigma and an inter-group variability \omega in my model.

To start with I set \mu=2, \sigma=0.75, and \omega=0.1, i.e. the inter-group variability (of the mean) is significantly smaller than the intra-group variability of observations (around the group mean).

For later experiments I can change num_groups, num_samples_per_group, and omega in my data generating process.

def generate_data(num_groups: int = 5, num_samples_per_group: int = 20, omega: float = 0.1) -> pd.DataFrame:
rng = np.random.default_rng(seed=41)
mu = 2.0
# intra-group varability of observations, same for all groups
sigma = 0.75
# inter-group variability in location mu_group
# omega = 0.1
result = {"group": [], "observation": []}
for i in range(num_groups):
mu_group = rng.normal(loc=mu, scale=omega)
result["group"] += [i] * num_samples_per_group
result["observation"] += rng.normal(loc=mu_group, scale=sigma, size=num_samples_per_group).tolist()
return pd.DataFrame(result)


The hierarchical model has population level variables \mu, \sigma, and \omega. It does not model a group specific variability of observations, i.e. \sigma_g. The group specific mean is modelled in the “standard” way (i.e. essentially copy & paste from the examples above).

def make_model_v1(frame: pd.DataFrame) -> pm.Model:
"""No group-specific volatility in the observation generating process, i.e. 'sigma_group = sigma_population'."""
coords = {"group": frame["group"].unique()}
with pm.Model(coords=coords) as model:
# Data
y = pm.MutableData("y", frame["observation"])
g = pm.MutableData("g", frame["group"])

# Population level parameters
mu = pm.Normal("mu")
sigma = pm.HalfNormal("sigma")
omega = pm.HalfNormal("omega", 10.0)

# Group level parameters
epsilon_mu_g = pm.ZeroSumNormal("epsilon_mu_g", dims="group")
mu_g = pm.Deterministic("mu_g", mu + omega * epsilon_mu_g)

# Likelihood
pm.Normal("y_obs", mu=mu_g[g], sigma=sigma, observed=y)
return model


For omega=0.1 the model seems to behave perfectly reasonable, e.g. increasing num_groups “narrows” the posterior on the true values for \mu, \sigma, and \omega, cf. table below

omega (true value) num_groups num_samples_per_group parameter mean hdi_3% hdi_97% ess_bulk ess_tail r_hat
0 0.1 5 20 mu 1.902 1.762 2.024 2382 1913 1
1 0.1 5 20 omega 0.238 0 0.587 883 980 1.01
2 0.1 5 20 sigma 0.706 0.607 0.801 2842 2585 1
3 0.1 50 20 mu 1.985 1.94 2.031 4798 2500 1
4 0.1 50 20 omega 0.071 0 0.137 1266 2026 1
5 0.1 50 20 sigma 0.772 0.739 0.806 4429 2337 1
6 0.1 500 20 mu 1.995 1.98 2.009 8185 2322 1
7 0.1 500 20 omega 0.096 0.07 0.119 1482 1985 1
8 0.1 500 20 sigma 0.751 0.74 0.761 5512 2852 1

One can also observe that the ESS increases with increasing number of groups. My interpretation is that the only evidence for omega is the between group variability and for a small number of groups there is just not much evidence that can be gathered. This in turn leads to a small ESS, especially for omega.

This is somewhat counterintuitive to the “there is no minimum sample size for Bayesian data analysis” statement since the presence of somewhat low ESS leaves me (as someone getting just started) doubtful that I will be able to trust my hierarchical model (due to problems with posterior sampling) in situations where I don’t control the data generating process. I guess I am overlooking something here though, but any input would be very welcome.

Now if I increase omega by a factor of ten and run the same experiment, then I observe that the posterior “hdi” does not contain the true value for 5 and 50 groups and the ESS values for omega drop significantly compared to the ones for the experiment with omega=0.1,

omega (true value) num_groups num_samples_per_group parameter mean hdi_3% hdi_97% ess_bulk ess_tail r_hat
0 1 5 20 mu 1.746 1.606 1.872 3210 1771 1
1 1 5 20 omega 1.112 0.413 2.077 660 841 1.01
2 1 5 20 sigma 0.706 0.61 0.802 2523 2146 1
3 1 50 20 mu 1.861 1.817 1.909 6323 2901 1
4 1 50 20 omega 0.922 0.737 1.117 501 670 1.01
5 1 50 20 sigma 0.773 0.741 0.807 4813 2549 1
6 1 500 20 mu 1.988 1.975 2.002 9388 2736 1
7 1 500 20 omega 0.991 0.93 1.051 264 665 1
8 1 500 20 sigma 0.75 0.741 0.761 7152 2465 1

For omega=10, the results look as follows

omega (true value) num_groups num_samples_per_group parameter mean hdi_3% hdi_97% ess_bulk ess_tail r_hat
0 10 5 20 mu 0.149 0.018 0.288 1905 1915 1
1 10 5 20 omega 8.656 3.897 13.736 594 879 1.01
2 10 5 20 sigma 0.705 0.615 0.802 1541 1921 1
3 10 50 20 mu 0.623 0.578 0.672 3588 2597 1
4 10 50 20 omega 8.89 7.173 10.57 132 127 1.03
5 10 50 20 sigma 0.772 0.741 0.804 3405 2963 1
6 10 500 20 mu 1.922 1.909 1.936 939 997 1.01
7 10 500 20 omega 9.888 9.433 10.931 5 12 2.11
8 10 500 20 sigma 0.75 0.741 0.761 414 853 1

I am not entirely sure what I should make of this. I could say that the standard error of the estimate for mu is essentially omega / sqrt(num_groups) and the deviation between the posterior mean and the true value seems to be compatible with that number. Also if omega is larger than sigma` does it make sense to talk about a (group-level) random effect, there is not much that those groups are sharing (in terms of the average of the observations within each group).

I’d be glad for any suggestions how to reconcile this intuition with the observed posterior summaries.