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

Questions
15

I agree with your opinion about a less restrictive hierarchical model having the potential to reveal more about the underlying data. I’ve tried to work my way up to the full model after reading this informative Primer on Bayesian Methods for Multilevel Modeling, but I’m stuck early on after getting to the user-level and trying to add the block-level.
Using the dims keyword argument and coords seems more elegant than the indexing scheme from this post.

Example data with multilevel predictors

    user	experience condition  block rating	direction	projection
0	9	    0.0        A      	  1	    2.0     parallel	-1.058687
1	9	    0.0        A      	  1	    2.0     parallel	0.201144
2	9	    0.0        A      	  1	    2.0     parallel	-0.806420
...

experience is a user-level predictor and rating a block-level “predictor” (actually, rating is measured afterwards, so it’s correlated with block).

The indexing/design matrix still poses a problem to me. I don’t understand how to use the predictors when doing outcome, predictors = dmatrices("projection ~ condition * user * block * direction", df). The example didn’t utilize the design matrix, stating the glm() function does not play nice with hierarchical models yet.
This is an excerpt from what I get with aforementioned formula.

Terms:
‘Intercept’ (column 0)
‘condition’ (columns 1:3)
‘user’ (columns 3:32)
‘block’ (columns 32:34)
‘direction’ (column 34)
‘condition:user’ (columns 35:93)
‘condition:block’ (columns 93:97)
‘user:block’ (columns 97:155)
‘condition:direction’ (columns 155:157)
‘user:direction’ (columns 157:186)
‘block:direction’ (columns 186:188)
‘condition:user:block’ (columns 188:304)
‘condition:user:direction’ (columns 304:362)
‘condition:block:direction’ (columns 362:366)
‘user:block:direction’ (columns 366:424)
‘condition:user:block:direction’ (columns 424:540)

From what I gathered from a tutorial in R about rat growth that @OriolAbril provided, in R this formula would be kind of this \sigma_{proj} \sim (condition * block * direction)/(1|trial) + (1|user), meaning trials are viewed as nested random effects, and users as random effects. It seems patsy doesn’t support nesting and random effects.

Following the primer I’ve come so far as to include users, as they are what counties are in the radon example. The difference is that I model the standard deviation sigma, not the mu. When trying to add rating as a pendent to uranium in the radon example, I ran into broadcasting issues.

Model with only direction x user:

conditions = df['condition'].cat.categories.tolist()
n_conditions = df['condition'].nunique()
condition_indices = df['condition'].cat.codes.values # Condition index for every single value.
condition_lookup = pd.Series(*pd.factorize(df['condition'].cat.categories), name='condition_idx_lookup')

# We need to be able to index by user. But the user IDs of the sample don't start at 0. We create an index and a mapping between them.
users = df['user'].cat.categories.tolist()
n_users = df['user'].nunique()
user_indices = df['user'].cat.codes.values  # User index for every single value.
user_lookup = pd.Series(*pd.factorize(df['user'].cat.categories), name='user_idx_lookup')
# User level predictors. More options: age group, gender.
# experience is not on an interval scale (no equidistance between items), but ordinal.
experience = df.drop_duplicates(subset="user")['experience']  # Shape = user_indices

blocks = df['block'].cat.categories.tolist()
n_blocks = df['block'].nunique()  # Same as using cat.categories.size
block_indices = df['block'].cat.codes.values # Block index for every single value.
block_lookup = pd.Series(*pd.factorize(df['block'].cat.categories), name='block_idx_lookup')
# Block level predictors.
# Note: Ratings are correlated with block index.
ratings = df.drop_duplicates(subset=["user", 'block'])['rating']  # Shape = block_indices

# Coded direction of projection: 0 = orthogonal, 1 = parallel.
directions = df['direction'].cat.categories.tolist()
n_directions = df['direction'].nunique()
direction_indices = df['direction'].cat.codes.values
direction_lookup = pd.Series(*pd.factorize(df['direction'].cat.categories), name='direction_idx_lookup')

coords = {'Direction': directions, 'User': users, 'Block': blocks, 'Condition': conditions, 'obs_id': np.arange(direction_indices
.size)}

with pm.Model(coords=coords) as varying_intercept_slope_by_user:
    direction_idx = pm.Data('direction_idx', direction_indices, dims='obs_id')
    user_idx = pm.Data('user_idx', user_indices, dims='obs_id')

    # Hyperpriors:
    intercept = pm.Gamma('intercept', mu=4.0, sigma=2.0)
    intercept_sd = pm.Exponential('sigma_intercept', 1.0)
    slope = pm.Normal('slope', mu=0.0, sigma=1.0)
    slope_sd = pm.Exponential('sigma_slope', 0.5)

    # Varying intercepts:
    intercept_user = pm.Gamma('intercept_user', mu=intercept, sigma=intercept_sd, dims='User')
    # Varying slopes:
    slope_user = pm.Normal('slope_user', mu=slope, sigma=slope_sd, dims='User')

    # Expected deviation per user:
    theta = intercept_user[user_idx] + slope_user[user_idx] * direction_idx

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

How do I add the other levels? Condition would be like state in the radon example. I try, but I keep hitting walls. In the meantime I try again and share my progress. I guess I have to formulate linear expressions for each level and use those as inputs to the levels below.