Model Comparison: How to constrain a model to positive mu-differences?

Questions
17

I tried adding a new level, but I’m stuck with the broadcasting issue. I really don’t know what I’m supposed to do with the dims or shape keyword arguments and the slicing of the parameters. Unfortunately, I didn’t find a comprehensible documentation on how this works anywhere. All the examples just have models with 1 or 2 levels, nothing nested deeper than that.

What I’d like to get is an estimate for intercept and slope of projection variance for Blocks 1,2, and 3, while each user can have their own intercepts and slopes in each block as different baselines (random effect).

In addition to the above I tried with blocks (0,1,2) and an ID for each single block.

coords['Block_ID'] = block_ids  # 0..n_users*n_blocks

with pm.Model(coords=coords) as by_block_model:
    user_idx = pm.Data('user_idx', user_indices, dims='obs_id')
    block_id_idx = pm.Data('block_id_idx', block_id_indices, dims='obs_id')
    block_idx = pm.Data('block_idx', block_indices, dims='obs_id')
    direction_idx = pm.Data('direction_idx', direction_indices, dims='obs_id')
    experience = pm.Data('experience', gaming_exp, dims='User')
    rating = pm.Data('rating', ratings, dims='Block_ID')  # Block has indices 0,1,2. But ratings are size of n_users*n_blocks. Try block_id?
    
    # Hyperpriors:
    global_intercept = pm.Gamma('global_intercept', mu=4.0, sigma=2.0)
    global_intercept_sd = pm.Exponential('global_intercept_sd', 1.0)
    global_slope = pm.Normal('global_slope', mu=1.5, sigma=1.0)  # We expect the slope to be positive most of the time.
    global_slope_sd = pm.Exponential('global_slope_sd', 0.5)

    # Users.
    # Varying intercepts per user:
    intercept_user = pm.Gamma('intercept_user', mu=global_intercept, sigma=global_intercept_sd, dims='User')
    intercept_user_sd = pm.Exponential('intercept_user_sd', 1.0, dims='User')
    # Varying slopes per user:
    slope_user = pm.Normal('slope_user', mu=global_slope, sigma=global_slope_sd, dims='User')
    slope_user_sd = pm.Exponential('slope_user_sd', 1.0, dims='User')

    # Blocks.
    # Varying intercepts per block:
    intercept_block = pm.Gamma('intercept_block', mu=intercept_user, sigma=intercept_user_sd, dims='Block')  # FixMe: ValueError: operands could not be broadcast together with shapes (3,) (30,) 
    # Varying slopes per block:
    slope_block = pm.Normal('slope_block', mu=slope_user, sigma=slope_user_sd, dims='Block')

    # Expected deviation per block:
    theta = intercept_block[block_idx] + slope_block[block_idx] * direction_idx

    projection = pm.Normal('projection', 0.0, sigma=theta, observed=df['projection'], dims='obs_id')

Once I know how to add more levels, I’d like to cluster users by condition, which may have an effect on standard deviation.