Mixed-effect Bayesian network reveals personal effects of nutrition

Nutrition experts know by their experience that people can react very differently to the same nutrition. If we could systematically quantify these differences, it would enable more personal dietary understanding and guidance. This work proposes a mixed-effect Bayesian network as a method for modeling the multivariate system of nutrition effects. Estimation of this network reveals a system of both population-wide and personal correlations between nutrients and their biological responses. Fully Bayesian estimation in the method allows managing the uncertainty in parameters and incorporating the existing nutritional knowledge into the model. The method is evaluated by modeling data from a dietary intervention study, called Sysdimet, which contains personal observations from food records and the corresponding fasting concentrations of blood cholesterol, glucose, and insulin. The model’s usefulness in nutritional guidance is evaluated by predicting personally if a given diet increases or decreases future levels of concentrations. The proposed method is shown to be comparable with the well-performing Extreme Gradient Boosting (XGBoost) decision tree method in classifying the directions of concentration increases and decreases. In addition to classification, we can also predict the precise concentration level and use the biologically interpretable model parameters to understand what personal effects contribute to the concentration. We found considerable personal differences in the contributing nutrients, and while these nutritional effects are previously known at a population level, recognizing their personal differences would result in more accurate estimates and more effective nutritional guidance.

In this work we propose Bayesian network as an appealing way to model and predict the effects of nutrition. Bayesian networks are directed graphical models that describe a joint probability distribution where nodes of the graph are random variables and edges between the nodes form a conditional probability structure. In nutritional modeling we assign nutritional composition of subjects' diet and different blood concentrations as those random variables, and study the connections between them. The goal is then to find such connections for the graph that most probably describe the correct conditional relationships between nutrients and their effects.

Personal nutrition data
We analyze a dataset from the Sysdimet study (Lankinen et al. 2011) that contains four repeated measurements of 17 nutrients, some basic information about the subject (gender, medication) and their corresponding blood concentrations, from 106 individuals.
Following spaghetti plots show the progress of blood concentrations during the 12-week study. Measurements of only ten from 106 subjects are shown at the plots for clarity. Our goal is to predict next personal blood concentration when subject's diet from past weeks is known.
When the future direction of the blood concentrations can be predicted, the typical and personal importance of nutrients contributing this change be estimated from the model parameters. Study timeline with laboratory measurements in weeks 0,4,8, and 12 Glucose (mmol/l) Supplementary Figure S 1: Progress of concentration levels during the 12-week study. Measurements were taken at weeks 0, 4, 8, and 12, and the food records were kept during the week before the measurement. The figure is plotted with ggplot2 package for R language (v 3.3.2, https://ggplot2.tidyverse.org).

Mixed-effect Bayesian network
We focus only on two-level, bipartite, graphs where nutrients and personal details are assumed to affect blood concentrations, and only the magnitude of the effect is left to be estimated. We assume a connection from every nutrient and personal information to every blood concentration value at the dataset. The model search is then focused on to estimating optimal coefficients describing the connections.

Model development
Our starting assumption is that both nutrients and the blood concentrations are normally distributed. This assumption may be naive as it allows the blood concentrations to have negative values. We normalize all the input values to same scale, so that the estimated regression coefficients can be used as a indicators of connection strenght.

Developing the model beyond normal distributions
More realistic probability distribution for blood concentrations would be Log-Normal or Gamma distributions (Limpert 2011). They both allow only positive values and model better the right tail of individual larger values. For further development, we choose Gamma distribution with identity link function. This is important as it keeps the regression coefficients of the nutrients on the same absolute scale than with Normal models. The regression coefficients are used as weights of the edges at the Bayesian network and they all should be on the same scale regardless of the random variables' distribution.

Modeling correlated observations with an autocorrelation process
In the dataset the observations from subjects' food records and blood concentrations have one week response time. This might not be optimal time for all the responses and the successive observations might be correlated. For modeling the correlated observations, we add an autocorrelation process to the model by adding an estimated portion of previous measurement to the linear predictor This coefficient 'ar1' is estimated separately for each response The autocorrelation structure seems to decrease the variance of the model and correcting the systematical estimation errors from the previous model candidate It is interesting also to compare the AR(1) coefficients of different local distributions at the graph. Blood insulin (fsins) has largest AR(1) coefficient indicating that successive measurements are most correlated and the values are not changing as rapidly as with other measurements. Credible interval for cholesterol measurements includes zero indicating that autocorrelation is very low or non-existent between two measurements. Supplementary

Skrinkage prior model
In our analysis, we use two models: 1) typical effects are better shown with previous model that has vague prior on beta coefficients, and 2) personal effects became clearer in a model that uses shrinkage prior on beta coefficients. This shrinkaged model is also faster to evaluate in k-fold setting and is possibly a better predictor.
The predictive ability of the model might be reduced if the model overfits to small nuances of the dataset.
To overcome the possible overfitting, we apply a shrinkage prior "Regularized Horseshoe Shrinkage (RHS)" (Piironen and Vehtari 2017), also called Finnish horseshoe, on regression coefficients. It allows specifying a prior knowledge, or at least an educated guess, about the number of significant predictors for a target. Here we assume that one third of the nutrients might be relevant for any given blood concenration, and apply following parameters for the shrinkage shrinkage_parameters <-within(list(), { scale_icept <-1 # prior std for the intercept scale_global <-0.02266 # scale for the half-t prior for tau: # ((p0=7) / (D=22-7)) * (sigma = 1 / sqrt(n=106*4)) nu_global <-1 # degrees of freedom for the half-t priors for tau nu_local <-1 # degrees of freedom for the half-t priors for lambdas slab_scale <-1 # slab scale for the regularized horseshoe slab_df <-1 # slab degrees of freedom for the regularized horseshoe }) Next we fit the previous model with the Finnish horseshoe added. Besides the shrinkage parameters, we provide now the data in two sets: input data for estimation and target data to hold out for prediction.

Evaluation of the models' fit and predictive performance
Different model versions are compared with a normalized root mean squared error (NRMSE) and also with a classification test that allows comparison with other classification methods. 8

Comparison with in-sample data
Following compares different model candidates by using normalized root mean squared error (NRMSE) metrics. Normalization allows comparing the prediction accuracy of responses with different scales.

Predictive accuracy cross-validated MEBN
Execution of the k-fold cross-validation is done in a separated notebook (cross-validation.Rmd) as it takes substantial amount of time to run. The results are stored and load for analysis here. For the cross-validation we partitioned the observations in data to 12 partitions for weeks 4,8 and 12. Then each of the data partitions in turn is left out from model estimation. This allows us to make predictions for unseen data.

Evaluation with classification metrics
Normalized root mean squared error (NRMSE) is good for comparing model candidates, but it does not say much about the actual usefulness of the model. For personal nutrition guiding, it is important to predict how given diet will affect the future blood concentration values. In the standard NRMSE comparison, the error between true and predicted values were compared, but a custom classification test can be also used. It tells if the given diet raises or lowers the future blood concentrations compared to past measurement.

Summary of the model comparison
This table summarizes how accurately we can predict if blood concenration increase or decrease when the current level and a diet are known.

Typical effects of nutrition and personal variance
We have now reached a reasonably good fit for the local distributions of the Bayesian network. Our code has extracted a graphical model from these posterior distributions of random variables and parameters. The graph consists mean values of the posteriors and also their credible intervals. This is a lighter data structure for further analysis and allows also graph operations besides the regression modeling.
Besides the random variables, the graph also includes the regression coefficients and that denote the typical and personal magnitudes of the effects for different nutrients. For better overview, let us plot the graph of typical nutritional effects.
We can examine closer the effects in the previous graph. In following, we show 90% credible interval for each effect and sort them by variance between persons. These effects with high variance between persons are interesting to us.
These are the effect with most difference between subjects. Let us the examine closer those of the effects that have most probable variance over 0.30 and plot them as a graph

Sensitivity analysis
In Table S4, we provide a sensitivity analysis on how the graph estimates with either vague or shrinkage priors on beta parameters. It can be seen that standard deviations of personal effects are not sensitive to prior choice. In the actual personal effects, the shrinkage prior puts more weight on personal variations and less weight on general part of the effect.
In addition, the effect parameters with most personal variation are provided in Table S5 HDL

Finding the clusters of personal reaction types
Although there are personal differences, it is likely that everyone is not behaving uniquely but there might exist similar groups of behavior. We can analyze personal differences in two ways. We can look the absolute magnitudes of effects and we can also look how persons differ from typical behavior or mean of the effect.
We can now take the personal estimations of the effects and see, if they form clusters The elbow plot of k-means clustering shows that difference between clusters decreases when more than four clusters are considered. These clusters include all the predictors, nutrients and the cholesterol lowering medication. The figure is plotted with ggplot2 package for R language (v 3.3.2, https://ggplot2.tidyverse.org).
By looking the previous elbow-plot diagram we see that there are four clearly identifiable groups between subjects.

The effect of cholesterol medication
At this plot, zero of X-axis denotes a typical behavior of that particular reaction, and bars denote a deviation from this typical mean.
Cholesterol medication seems to dominate the clusters. There are subjects for who cholesterol medication raises the blood insulin levels more than on average and for other group the raise is less than on average. By looking at the absolute effect of the cholesterol medication it seems that there separates two groups where one it is a significant factor in explaining blood insulin level, and other group where it does not play virtually any role. But does these subjects even take cholesterol medications? Let's create a cross tabulation of observations to explore this effect more closely.
This indicates that subjects in clusters 2 and 4 are indeed taking cholesterol medication.
The average insulin concentrations in these clusters are 18 cluster 3 (5/106, 4.7%) cluster 4 (7/106, 6.6%) cluster 1 (83/106, 78.3%) cluster 2 (11/106, 10. The effect of cholesterol-lowering medication to the insulin concentration dominates the personally varying effects. Besides decreasing the cholesterol concentration, the medication has a side effect of increasing insulin concentration for most of the subjects. Yet, there exists a large variation. The increasing effect is dramatic for 6,6% of the subjects in cluster 4, but for 10,4% of subjects in cluster 2 the effect is non-existing. The figure is plotted with ggplot2 package for R language (v 3.3.2, https://ggplot2.tidyverse.org). There are 9 (=36/4) and 7 (=28/4) subjects in clusters 2 and 4 who are taking cholesterol medication. For subjects in the cluster 4 the average insulin level is the medication seems to increase the insulin concentration quite much, but for subjects in the cluster 2 not much at all.

Clustering with nutritional effects only
Both gender and cholesterol medication are unchanging factors at the dataset. Let us next remove those from clustering features and see how the clusters form based on nutritional effects only. In clusters 1 and 2 there are subjects whose insulin concentrations decrease and increase more than on average.

Personally predicted reaction types
During the cross-validation we estimated personal models for each subject. Next we illustrate the differences within personal reaction types by picking few most distant ones.
We can then compare the visualizations of the personal reaction graphs. The effect of cholesterol medication is most significant difference, but other differences exist as well.