Anova - MathAch dataset of 160 schools with math achievements of students

Questions v5
,
3

Maybe the sigma parameter of your mu is causing the issue? A 0.01 SD parameter for a Gaussian is quite narrow, so maybe a standard Gaussian provides better coverage (always better to conduct prior predictive checks). Just by changing that parameter you can recover the value of the central tendency of simulated data:

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

np.random.seed(27)

M = np.random.normal(12, 3, 160*32) #mean 12 and SD 3

school_idx = np.array([np.repeat(i, 32) for i in range(160)]).flatten()

schools = np.arange(160).astype("str")

coords = {"school": schools,"obs_id":np.arange(len(M))}

with pm.Model(coords=coords) as model:
    school = pm.ConstantData("shool_idx", school_idx, dims="obs_id")
    sigmas = pm.Uniform('sigmas', 0, 10, shape=2)
    tau = pm.Deterministic('tau', sigmas**(-2))
    mu = pm.Normal('mu', mu=0, sigma=1)
    alphas = pm.Normal('alphas', mu=mu, sigma=tau[1]**0.5, dims='school')
    math = pm.Normal('math', mu=alphas[school], sigma=tau[0]**0.5, observed=M, dims='obs_id')
    
with model:
    idata = pm.sample(nuts_sampler="numpyro", random_seed=27)
    
az.summary(idata, var_names=['mu','sigmas','tau'])

Running chain 0: 100%|██████████| 2000/2000 [00:06<00:00, 325.59it/s]
Running chain 1: 100%|██████████| 2000/2000 [00:06<00:00, 325.70it/s]
Running chain 2: 100%|██████████| 2000/2000 [00:06<00:00, 325.86it/s]
Running chain 3: 100%|██████████| 2000/2000 [00:06<00:00, 326.02it/s]
Sampling time =  0:00:06.318887
Transforming variables...
Transformation time =  0:00:00.116303
Out[31]: 
             mean     sd  hdi_3%  hdi_97%  ...  mcse_sd  ess_bulk  ess_tail  r_hat
mu         11.997  0.045  11.916   12.079  ...    0.002     389.0     623.0   1.01
sigmas[0]   0.335  0.003   0.328    0.341  ...    0.000    5458.0    3010.0   1.00
sigmas[1]   7.631  1.544   5.048    9.997  ...    0.082     184.0     318.0   1.03
tau[0]      8.897  0.177   8.556    9.227  ...    0.002    5458.0    3010.0   1.00
tau[1]      0.020  0.010   0.010    0.039  ...    0.001     184.0     318.0   1.03

[5 rows x 9 columns]

May be also possible to remove the sigmas and assign tau directly as a HalfNormal? That may help to better recover the data total SD (if that’s desired):

with pm.Model(coords=coords) as model:
    school = pm.ConstantData("shool_idx", school_idx, dims="obs_id")
    tau = pm.HalfNormal('tau', 1, shape=2) 
    mu = pm.Normal('mu', mu=0, sigma=1)
    alphas = pm.Normal('alphas', mu=mu, sigma=tau[1], dims='school')
    math = pm.Normal('math', mu=alphas[school], sigma=tau[0], observed=M, dims='obs_id')
    
with model:
    idata = pm.sample(nuts_sampler="numpyro", random_seed=27)
    
az.summary(idata, var_names=['mu','tau'])

Running chain 0: 100%|██████████| 2000/2000 [00:07<00:00, 277.07it/s]
Running chain 1: 100%|██████████| 2000/2000 [00:07<00:00, 277.15it/s]
Running chain 2: 100%|██████████| 2000/2000 [00:07<00:00, 277.26it/s]
Running chain 3: 100%|██████████| 2000/2000 [00:07<00:00, 277.42it/s]
Sampling time =  0:00:07.397331
Transforming variables...
Transformation time =  0:00:00.078383
Out[34]: 
          mean     sd  hdi_3%  hdi_97%  ...  mcse_sd  ess_bulk  ess_tail  r_hat
mu      11.999  0.047  11.908   12.085  ...    0.003     134.0     401.0   1.04
tau[0]   2.982  0.030   2.922    3.035  ...    0.001    1507.0    1982.0   1.00
tau[1]   0.115  0.059   0.025    0.216  ...    0.012      12.0      22.0   1.30

[3 rows x 9 columns]

Hope this makes sense and helps a bit.

3 Likes