Specifying a simple heirarchical model for nutritional exposures

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:

R(X, W) = expit(\alpha + X\boldsymbol{\beta} + W\boldsymbol{\gamma})

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:

\beta_j = \pi_1 z_{1j} + \pi_2 z_{2j}...\pi_p z_{pj} + \delta_j

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:

BP_i ~ Normal(\mu_i, sigma)
\mu_i = \alpha + \beta_1 * milk
\beta_1 = \pi_1 * fat + \pi_2 * sugar + \pi_3 * calcium + \pi_3 * protein + \delta

And after reading a bit about the brms multivariate syntax, the model code could potentially be:

  brm(bf(BP ~ 1 + milk) +
        bf(milk ~ fat + sugar + calcium + protein) +
        set_rescor(rescor = FALSE)) 

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

X_1\beta_1 + X_2 \beta_2 = X_1\delta_1 + X_2\delta_2 + (X_1z_{1j} + X_2z_{2j})\pi_1 + (X_1z_{1j} + X_2z_{2j})\pi_2

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.

Best of luck with your model!

1 Like