Hierarchical models with several of coefficients

I’m trying to fit a hierarchical model using data from 30 features and 5 groups (code below). Since I do not want to define each of the 30 features’ distributions, I’m using the shape argument to speed up things. This will require to create a matrix of size (n_groups, n_features) to store the coefficients of each feature per group. I’ve two questions:

  1. When I run this, I get an error: Exception: ("Compilation failed (return status=1). I guess I’m doing something wrong when defining the model, but I cannot figure out what’s the mistake.
  2. My understanding is that the code below defines a model that assumes that the coefficients are independent of each other. If I wanted to consider dependence among coefficients, what would need to be modified?

Any reference / guidance for these questions is highly appreciated.

# X is a matrix of size (N, 30)
# y is a vector of length N
# group_id is a vector of length N with n_groups unique values
n_features = 30
n_groups = 5

with pm.Model() as hierarchical_model:
    # Hyperpriors for group nodes
    mu_a = pm.Normal('mu_a', mu=0., sigma=10)
    sigma_a = pm.HalfNormal('sigma_a', sigma=20)
    mu_b = pm.Normal('mu_b', mu=0., sigma=10, shape=n_features)
    sigma_b = pm.HalfNormal('sigma_b', sigma=20, shape=n_features)
    
    a_offset = pm.Normal('a_offset', mu=0, sd=1, shape=n_groups)
    a = pm.Deterministic("a", mu_a + a_offset * sigma_a)
    
    b_offset = pm.Normal('b_offset', mu=0, sd=1, shape=(n_groups, n_features))
    b = pm.Deterministic("b", mu_b + b_offset * sigma_b)
    
    # Model Error
    noise = pm.Gamma('noise', alpha=2, beta=1)
    
    estimate = a[group_id] + X @ b[group_id]
    
    y_observed = pm.Normal('y_observed',
                           mu=estimate,
                           sigma=noise,
                           observed=y)

Hi Roberto,
I don’t see any obvious issue in your model definition – although I’d probably use more informative hyperpriors.
Is the error appearing when you’re trying to compile the model, or when you’re trying to sample? Also, what are the PyMC3 and Theano versions on your virtual env?

Hi Alex,
The version of PyMC3 is 3.8 and Theano is 1.0.4. The error appears when I compile the model.

Exception: ("Compilation failed (return status=1): /Users/user/.theano/compiledir_macOS-10.14.6-x86_64-i386-64bit-i386-3.8.5-64/tmpbygutbi4/mod.cpp:510:27: 
error: non-constant-expression cannot be narrowed from type 'npy_intp' (aka 'long') to 'int' in initializer list [-Wc++11-narrowing].     int init_totals[2] = {V3_n0, V3_n1};.                          
^~~~~. /Users/user/.theano/compiledir_macOS-10.14.6-x86_64-i386-64bit-i386-3.8.5-64/tmpbygutbi4/mod.cpp:510:27: 
note: insert an explicit cast to silence this issue.     int init_totals[2] = {V3_n0, V3_n1};.  

Btw, can you expand on the more informative hyperpriors comment?

Thanks Roberto. Can you update them to the latest versions and retry please?
conda update pymc3 theano -c conda-forge
It should be 3.9.3 and 1.0.5 after upgrading.

Regarding priors, it’s very important to be careful about them when using a hierarchical model, for various reasons that would be too long to detail here. You can take a look these recommendations, and at Richard McElreath’s Statistical Rethinking book (PyMC3 port).
Spoiler alert: I’ll talk about choice of priors in my PyMCon talk :wink:

Hi Alex. Just updated both PyMC3 and Theano to the latest versions, nevertheless, a ‘shorter’ error appears:

Exception: ('Compilation failed (return status=254): clang-10: error: 
unable to execute command: Segmentation fault: 11. clang-10: error: 
dsymutil command failed due to signal (use -v to see invocation). ', 
'[Elemwise{sub,no_inplace}(b_offset, <TensorType(float64, (True, True))>)]')

The error points out this line b_offset = pm.Normal('b_offset', mu=0, sd=1, shape=(n_stores, n_features)), but I cannot see what may be wrong.

Are you sure n_stores and n_features are int and not float?
Also, try casting mu and sigma to floats: pm.Normal('b_offset', mu=0.0, sigma=1.0, shape=(n_stores, n_features))

Yes, both n_stores and n_features are int. I also modified the parameters of the distributions to be floats. Still, the same error as above is shown during model compilation.

Ok, thanks. It looks like it could be a local compilation error, so not a Theano or PyMC3 issue – pinging @brandonwillard, as he knows much more about this than me!

Sorry to ask you that, but could you create a brand new conda virtual env using only the conda-forge channel? Meaning that when you do conda env create ... and conda install ... you make sure to add -c conda-forge. That way, we’ll be sure python, the C compiler, PyMC3 and Theano all come from the same source, because this could be the issue.

If this still doesn’t solve it, could you share some data so that we can run your model and try to reproduce the issue?

Looks like this is the issue you’re facing: https://github.com/pymc-devs/Theano-PyMC/issues/127