HI all, I posted this question recently on cross-validated but have realised this is a much better place for this question. There’s a hierarchical model for the effect of food intake on disease risk described in the textbook ‘Modern Epidemiology’ (3rd ed, pages 435-437) that I would love to learn how to implement.

Here’s the typical (non-hierarchical) model for studying the effect of intake of a given food on disease risk:

Where \alpha represents the intercept, X is a matrix of food items, and W is a matrix of confounders.

We can include more detail in this model and improve our effect estimates by including a model for some or all of the coefficients in \mathbf{\beta}. Specifically for each food item X_j we have a row vector \mathbf{Z_j} = (z_{1j}, z_{2j}... z_{pj}) containing the nutrient content of that food (for instance, fats, carbohydrates, protein, etc). We can then model the food effects \beta_j as a function of the nutrient content of the foods:

In this hierarchical model, the terms \pi_k z_{kj} represent the contribution of nutrient k to the effect \beta_j of food item j on risk of disease, while \delta_j represents the residual effect of food j not captured by the sum of its constitutive nutrients.

I really like this model and would love to try it out. The problem is though, I have no idea how to implement it in Stan or brms. If anyone can help me clarify this, that would be greatly appreciated!

For a concrete example, I’ll use a continuous outcome for simplicity. Say, blood pressure. We can consider the food item to be serves of milk. Ignoring priors, I believe the model looks like:

But I am very much not sure if this is correct. I have had very little exposure to multivariate models and need to do some reading. Let me know if i’m way off the mark!

EDIT: Now that I think about it, that won’t work because the nutrients from food tables will be perfectly correlated with the food item. I’m really lost with this model…

If I understand you correctly, this is actually exactly the same as adding the sum of nutritional contents of all the foods (weighted by X) as additional predictors, because

so you would have something like BP ~ 1 + milk + butter + ... + total_sugar + total_calcium + .... Now the brms coefficents b_milk and b_butter correspond to the \delta_j and b_total_sugar and b_total_calcium to the \pi_k.

Does that make sense?

I am not sure if you are having additional problems or if this was the main issue - feel free to ask for clarifications.