Dietary intake, Nutritional status, and Health outcomes among Vegan, Vegetarian and Omnivore families: results from the observational study

Statistical report - mixed models summary

Authors and affiliations

Marina Heniková1,2, Anna Ouřadová1, Eliška Selinger1,3, Filip Tichanek4, Petra Polakovičová4, Dana Hrnčířová2, Pavel Dlouhý2, Martin Světnička5, Eva El-Lababidi5, Jana Potočková1, Tilman Kühn6, Monika Cahová4, Jan Gojda1


1 Department of Internal Medicine, Kralovske Vinohrady University Hospital and Third Faculty of Medicine, Charles University, Prague, Czech Republic.
2 Department of Hygiene, Third Faculty of Medicine, Charles University, Prague, Czech Republic.
3 National Health Institute, Prague, Czech Republic.
4 Institute for Clinical and Experimental Medicine, Prague, Czech Republic.
5 Department of Pediatrics, Kralovske Vinohrady University Hospital and Third Faculty of Medicine, Charles University, Prague, Czech Republic.
6 Department of Epidemiology, MedUni, Vienna, Austria.

This is a statistical report of the study currenlty under review in the Communications Medicine journal.

When using this code or data, cite the original publication:

TO BE ADDED

BibTex citation for the original publication:

TO BE ADDED

Original GitHub repository: https://github.com/filip-tichanek/kompas_clinical

Statistical reports can be found on the reports hub.

Data analysis is described in detail in the statistical methods report.

1 Introduction

This project is designed to evaluate and compare clinical outcomes across three distinct dietary strategy groups:

  • Vegans
  • Vegetarians
  • Omnivores

The dataset includes both adults and children, with data clustered within families.

1.1 Main Questions

The study addresses the following key questions:

Q1. Do clinical outcomes vary significantly across different diet strategies?

Q2. Beyond diet group, which factors (e.g., sex, age, breastfeeding status for children, or supplementation when applicable) most strongly influence clinical outcomes? How correlated (“clustered”) are these characteristics within the same family?

Q3. Could the clinical characteristics effectively discriminate between different diet groups?

1.2 Statistical Methods

For full methodological details, see this report. In brief:

  • Robust linear mixed-effects models (rLME) were used to estimate adjusted differences between diet groups (Q1) and assess the importance of other variables (Q2), including how much clinical characteristics tend to cluster within families. Covariates included age, sex, breastfeeding status for children, and relevant supplementation factors where applicable.

  • Elastic net logistic regression was employed to answer Q3, evaluating whether clinical characteristics provide a strong overall signal distinguishing between diet groups, incorporating a predictive perspective.

All analyses were conducted separately for adults and children.

2 Robust linear-mixed-effects models

Open code 
source('r/368_initiation.R')

2.1 Children

2.1.1 Diet effects

2.1.1.1 Unstandardized effects

Open code 
diet_child_all <- read.xlsx("gitignore/data/diet_child_all_mixeff.xlsx") %>%
  mutate(outcome = clean_outcome_names(outcome)) %>%
  mutate(outcome = case_when(
    outcome == "FOLAT" ~ "FOLATE",
    outcome == "VKFE" ~ "TIBC",
    outcome == "Cros" ~ "CTX",
    outcome == "log2_Cros" ~ "log2_CTX",
    outcome == "VIT_AKTB12" ~ "HOLOTC",
    outcome == "log2_VIT_AKTB12" ~ "log2_HOLOTC",
    outcome == "log2_Ur_Ca" ~ "log2_uCa",
    outcome == "log2_Ca_per_Krea" ~ "log2_uCCR",
    outcome == "log2_I_per_Krea" ~ "log2_uICR",
    outcome == "P_per_Krea" ~ "uPCR",
    outcome == "log2_P_per_Krea" ~ "log2_uPCR",
    TRUE ~ outcome
  ))

diet_child_all2 <- diet_child_all %>%
  mutate(across(.cols = c(VN_OM_diff, VG_OM_diff, VN_VG_diff), ~ round(.x, 2)))


kableExtra::kable(diet_child_all2,
  caption = "Differences in clinical outcomes (`outcome`) between diet groups in children, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`)."
) %>%
  kable_styling("striped", full_width = T) %>%
  add_header_above(c(
    " " = 2, "(A) VN vs OM" = 3,
    "(B) VG vs OM" = 3, "(C) VN vs VG" = 3
  )) %>%
  column_spec(1, width_min = "1.5in")
Differences in clinical outcomes (`outcome`) between diet groups in children, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
MASS_Perc mean -9.23 0.1368814 0.3220738 -3.60 0.6140181 0.9863118 -5.63 0.3798854 0.7436615
HEIGHT_Perc mean -7.97 0.2305660 0.3842767 0.86 0.9114004 0.9863118 -8.83 0.1967938 0.6361680
BMI_PERC mean -2.32 0.7047050 0.8290647 -2.79 0.6940089 0.9863118 0.47 0.9403351 0.9403351
M_per_H_PERC mean -3.17 0.5747725 0.7184656 -3.66 0.5716597 0.9863118 0.49 0.9331269 0.9403351
GLY mean 0.01 0.9195424 0.9566299 0.05 0.6481973 0.9863118 -0.04 0.6997053 0.8672515
TC mean -0.32 0.0370824 0.1236081 -0.18 0.3083592 0.9856067 -0.14 0.4090138 0.7436615
HDL mean -0.01 0.8826851 0.9542542 -0.06 0.4997270 0.9863118 0.05 0.5453018 0.8078545
LDL mean -0.25 0.0288423 0.1176745 0.08 0.5545440 0.9863118 -0.33 0.0089436 0.0715485
log2_TG mean -0.01 0.9566299 0.9566299 -0.31 0.0884329 0.9092226 0.30 0.0712944 0.3664534
Ca mean -0.02 0.4939300 0.6585733 0.01 0.7974114 0.9863118 -0.02 0.3579991 0.7436615
P mean -0.03 0.4692479 0.6472385 0.00 0.9409113 0.9863118 -0.03 0.4533815 0.7487340
Mg mean 0.03 0.0303454 0.1176745 0.03 0.0909223 0.9092226 0.00 0.8451016 0.9099736
log2_Se mean -0.17 0.1938190 0.3691790 -0.02 0.9174542 0.9863118 -0.15 0.2765187 0.6563408
log2_Zn mean -0.13 0.0481466 0.1374279 -0.08 0.2842636 0.9856067 -0.05 0.5226594 0.8040914
FE mean 1.81 0.1888387 0.3691790 -0.77 0.6250206 0.9863118 2.58 0.0874415 0.3664534
TIBC mean -4.02 0.0323605 0.1176745 -1.47 0.4931661 0.9863118 -2.54 0.2160265 0.6361680
log2_FERR mean -0.25 0.1295455 0.3220738 -0.19 0.3142560 0.9856067 -0.06 0.7459186 0.8775513
TRF mean -0.17 0.0242801 0.1176745 -0.07 0.4072626 0.9863118 -0.10 0.2390322 0.6374193
SATTRF mean 3.82 0.0515355 0.1374279 -0.03 0.9885952 0.9885952 3.85 0.0738694 0.3664534
log2_TRFINDEX mean 0.03 0.6950900 0.8290647 0.08 0.4471824 0.9863118 -0.04 0.6593668 0.8578668
log2_STRF mean -0.03 0.5724161 0.7184656 0.01 0.8399455 0.9863118 -0.04 0.4679587 0.7487340
HGB mean -0.44 0.8101738 0.9238143 1.14 0.5912508 0.9863118 -1.58 0.4337484 0.7487340
MCV mean 2.97 0.0001520 0.0020269 -0.13 0.8923264 0.9863118 3.10 0.0002773 0.0038755
PTH mean 0.38 0.2039043 0.3707351 -0.02 0.9616540 0.9863118 0.39 0.2226588 0.6361680
log2_CTX mean 0.06 0.4638317 0.6472385 0.08 0.3449623 0.9856067 -0.03 0.7154825 0.8672515
P1NP mean 41.88 0.2895625 0.4454808 63.17 0.1819924 0.9099621 -21.29 0.6218130 0.8576731
log2_UI mean -0.79 0.0005746 0.0057464 -0.37 0.1371521 0.9099621 -0.41 0.0725949 0.3664534
UREA mean -0.55 0.0104269 0.0595822 -0.17 0.4881839 0.9863118 -0.38 0.1075133 0.3909573
log2_CREA mean -0.11 0.0416995 0.1283063 -0.06 0.3350035 0.9856067 -0.05 0.4003741 0.7436615
UA mean -2.38 0.8314328 0.9238143 -4.94 0.7086404 0.9863118 2.56 0.8310370 0.9099736
log2_HOLOTC mean 0.30 0.1537853 0.3417452 -0.31 0.1449102 0.9099621 0.61 0.0002554 0.0038755
log2_HCY mean -0.14 0.2222796 0.3842767 0.02 0.9025532 0.9863118 -0.16 0.0916134 0.3664534
log2_MMA mean -0.58 0.0076679 0.0511195 -0.07 0.7609094 0.9863118 -0.51 0.0021082 0.0210824
log2_VIT_D mean 0.28 0.0064921 0.0511195 0.26 0.0328104 0.6562074 0.02 0.8644749 0.9099736
FOLATE mean 4.30 0.0000071 0.0001412 4.81 0.0000182 0.0007284 -0.51 0.6198255 0.8576731
log2_IGF1 mean -0.16 0.2883477 0.4454808 0.01 0.9512789 0.9863118 -0.17 0.3024966 0.6722147
log2_uCa mean -0.42 0.3236212 0.4794388 0.04 0.9403988 0.9863118 -0.46 0.2789448 0.6563408
log2_uCCR mean 0.02 0.9484718 0.9566299 0.19 0.6604621 0.9863118 -0.17 0.6648468 0.8578668
log2_uICR mean -0.38 0.1787353 0.3691790 -0.31 0.3233064 0.9856067 -0.07 0.7991710 0.9099736
uPCR mean -2.12 0.0000002 0.0000088 -0.63 0.1606090 0.9099621 -1.48 0.0002907 0.0038755

2.1.1.2 Standardized effects by 1SD

Open code 
meds <- diet_child_all %>% filter(estimand == "mean")

ids <- nrow(meds)
sds <- vector("double", ids)
sds 
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [39] 0 0

before_outcome_col <- which(colnames(dat_child_all) == 'aMASS_Perc')-1

sds <- vector("double", ids)

for(i in 1:ids){
  sds[i] <- ifelse(grepl("log2", diet_child_all[(i), 'outcome']),
                   sd(log2(dat_child_all[,(i + before_outcome_col)]), na.rm = TRUE),
                   sd(dat_child_all[,(i + before_outcome_col)], na.rm = TRUE))
}


diet_child_all_Z <- diet_child_all %>% 
  mutate(SD = sds) %>% 
  mutate(
    VN_OM_diff = VN_OM_diff/SD,
    VG_OM_diff = VG_OM_diff/SD,
    VN_VG_diff = VN_VG_diff/SD) %>% 
  select(-c(SD)) %>% 
  mutate(across(.cols = c(VG_OM_diff, VN_OM_diff, VN_VG_diff), ~ round(.x, 2)))


kbl(diet_child_all_Z, caption = 
      "Standardized differences in clinical outcomes (`outcome`) between diet groups in children, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).") %>% 
  kable_styling("striped",full_width = T) %>% 
  add_header_above(c(" " = 2, "(A) VN vs OM" = 3, 
                     "(B) VG vs OM" = 3, "(C) VN vs VG" =3)) %>% 
  column_spec(1, width_min = '1.5in')
Standardized differences in clinical outcomes (`outcome`) between diet groups in children, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
MASS_Perc mean -0.32 0.1368814 0.3220738 -0.12 0.6140181 0.9863118 -0.20 0.3798854 0.7436615
HEIGHT_Perc mean -0.29 0.2305660 0.3842767 0.03 0.9114004 0.9863118 -0.32 0.1967938 0.6361680
BMI_PERC mean -0.09 0.7047050 0.8290647 -0.11 0.6940089 0.9863118 0.02 0.9403351 0.9403351
M_per_H_PERC mean -0.13 0.5747725 0.7184656 -0.15 0.5716597 0.9863118 0.02 0.9331269 0.9403351
GLY mean 0.02 0.9195424 0.9566299 0.09 0.6481973 0.9863118 -0.07 0.6997053 0.8672515
TC mean -0.49 0.0370824 0.1236081 -0.28 0.3083592 0.9856067 -0.21 0.4090138 0.7436615
HDL mean -0.04 0.8826851 0.9542542 -0.19 0.4997270 0.9863118 0.15 0.5453018 0.8078545
LDL mean -0.45 0.0288423 0.1176745 0.14 0.5545440 0.9863118 -0.58 0.0089436 0.0715485
log2_TG mean -0.01 0.9566299 0.9566299 -0.43 0.0884329 0.9092226 0.42 0.0712944 0.3664534
Ca mean -0.14 0.4939300 0.6585733 0.06 0.7974114 0.9863118 -0.21 0.3579991 0.7436615
P mean -0.16 0.4692479 0.6472385 0.02 0.9409113 0.9863118 -0.18 0.4533815 0.7487340
Mg mean 0.40 0.0303454 0.1176745 0.36 0.0909223 0.9092226 0.04 0.8451016 0.9099736
log2_Se mean -0.33 0.1938190 0.3691790 -0.03 0.9174542 0.9863118 -0.30 0.2765187 0.6563408
log2_Zn mean -0.43 0.0481466 0.1374279 -0.28 0.2842636 0.9856067 -0.16 0.5226594 0.8040914
FE mean 0.26 0.1888387 0.3691790 -0.11 0.6250206 0.9863118 0.37 0.0874415 0.3664534
TIBC mean -0.44 0.0323605 0.1176745 -0.16 0.4931661 0.9863118 -0.28 0.2160265 0.6361680
log2_FERR mean -0.28 0.1295455 0.3220738 -0.22 0.3142560 0.9856067 -0.07 0.7459186 0.8775513
TRF mean -0.47 0.0242801 0.1176745 -0.20 0.4072626 0.9863118 -0.27 0.2390322 0.6374193
SATTRF mean 0.39 0.0515355 0.1374279 0.00 0.9885952 0.9885952 0.39 0.0738694 0.3664534
log2_TRFINDEX mean 0.06 0.6950900 0.8290647 0.14 0.4471824 0.9863118 -0.08 0.6593668 0.8578668
log2_STRF mean -0.11 0.5724161 0.7184656 0.04 0.8399455 0.9863118 -0.15 0.4679587 0.7487340
HGB mean -0.05 0.8101738 0.9238143 0.12 0.5912508 0.9863118 -0.17 0.4337484 0.7487340
MCV mean 0.78 0.0001520 0.0020269 -0.03 0.8923264 0.9863118 0.81 0.0002773 0.0038755
PTH mean 0.30 0.2039043 0.3707351 -0.01 0.9616540 0.9863118 0.31 0.2226588 0.6361680
log2_CTX mean 0.15 0.4638317 0.6472385 0.24 0.3449623 0.9856067 -0.08 0.7154825 0.8672515
P1NP mean 0.14 0.2895625 0.4454808 0.21 0.1819924 0.9099621 -0.07 0.6218130 0.8576731
log2_UI mean -0.63 0.0005746 0.0057464 -0.30 0.1371521 0.9099621 -0.33 0.0725949 0.3664534
UREA mean -0.44 0.0104269 0.0595822 -0.14 0.4881839 0.9863118 -0.30 0.1075133 0.3909573
log2_CREA mean -0.27 0.0416995 0.1283063 -0.15 0.3350035 0.9856067 -0.12 0.4003741 0.7436615
UA mean -0.05 0.8314328 0.9238143 -0.10 0.7086404 0.9863118 0.05 0.8310370 0.9099736
log2_HOLOTC mean 0.37 0.1537853 0.3417452 -0.38 0.1449102 0.9099621 0.75 0.0002554 0.0038755
log2_HCY mean -0.32 0.2222796 0.3842767 0.03 0.9025532 0.9863118 -0.36 0.0916134 0.3664534
log2_MMA mean -0.64 0.0076679 0.0511195 -0.07 0.7609094 0.9863118 -0.56 0.0021082 0.0210824
log2_VIT_D mean 0.61 0.0064921 0.0511195 0.57 0.0328104 0.6562074 0.04 0.8644749 0.9099736
FOLATE mean 0.92 0.0000071 0.0001412 1.03 0.0000182 0.0007284 -0.11 0.6198255 0.8576731
log2_IGF1 mean -0.19 0.2883477 0.4454808 0.01 0.9512789 0.9863118 -0.20 0.3024966 0.6722147
log2_uCa mean -0.29 0.3236212 0.4794388 0.02 0.9403988 0.9863118 -0.31 0.2789448 0.6563408
log2_uCCR mean 0.02 0.9484718 0.9566299 0.13 0.6604621 0.9863118 -0.11 0.6648468 0.8578668
log2_uICR mean -0.29 0.1787353 0.3691790 -0.24 0.3233064 0.9856067 -0.06 0.7991710 0.9099736
uPCR mean -0.99 0.0000002 0.0000088 -0.30 0.1606090 0.9099621 -0.70 0.0002907 0.0038755

2.1.1.3 Vulcano plots

VN vs OM

Open code 
color_map <- c(
  "OM+" = rgb(213/255, 94/255, 0),
  "VN+" = rgb(0, 158/255, 75/255),
  "VG+" = rgb(0, 95/255, 195/255),
  "P>0.05" = "grey50"
)

plotac <- "fig19"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_child_all_Z %>%
    filter(estimand != "log(OR)") %>%
    mutate(
      pv_adj = VN_OM_P_adj,
      pv_raw = VN_OM_P,
      b_coe = VN_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VN+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.3,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs. OM, all children") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VG vs OM

Open code 

plotac <- "fig20"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_child_all_Z %>%
    filter(estimand != "log(OR)") %>%
    mutate(VG_OM_P = if_else(VG_OM_P < 1e-5, 1e-5, VG_OM_P)) %>%
    mutate(
      pv_adj = VG_OM_P_adj,
      pv_raw = VG_OM_P,
      b_coe = VG_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VG+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VG+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.3,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VG vs. OM, all children") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VN vs VG

Open code 
plotac <- "fig21"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_child_all_Z %>%
    filter(estimand != "log(OR)") %>%
    
    mutate(
      pv_adj = VN_VG_P_adj,
      pv_raw = VN_VG_P,
      b_coe = VN_VG_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "VG+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c(
          "VG+",
          "VN+",
          "P>0.05"
        )
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.5,
      size = 3.1,
      seed = 16
    ) +
    labs(
      x = "Standardized difference",
      y = "-log10 (P-value)"
    ) +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs VG, all children") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

2.1.2 Covariates importance

2.1.2.1 Table

Open code 
AIC_child_all <- read.xlsx("gitignore/data/AIC_child_all_mixeff.xlsx") %>%
  mutate(across(where(is.numeric), ~ round(., 1)),
    outcome = clean_outcome_names(outcome)
  ) %>%
  rename(
    SEX = aSEX,
    BREASTFEED = Breast_Feed,
    `OTHER COV` = other_cov,
    FAMILY = family,
    DIET = diet,
    AGE = log2_age
  ) %>%
  mutate(outcome = case_when(
    outcome == "FOLAT" ~ "FOLATE",
    outcome == "VKFE" ~ "TIBC",
    outcome == "Cros" ~ "CTX",
    outcome == "log2_Cros" ~ "log2_CTX",
    outcome == "VIT_AKTB12" ~ "HOLOTC",
    outcome == "log2_VIT_AKTB12" ~ "log2_HOLOTC",
    
    outcome == 'log2_aUr_Ca' ~ 'log2_uCa',
    outcome == 'log2_Ca_per_Krea' ~ 'log2_uCCR',
    outcome == 'log2_I_per_Krea' ~ 'log2_uICR',
    outcome == 'P_per_Krea' ~ 'uPCR',
    outcome == 'log2_P_per_Krea' ~ 'log2_uPCR',
    
    TRUE ~ outcome
  ))

kbl(AIC_child_all,
  caption =
    "Change of the Akaike Information Criterion (AIC) after addition of a given variable to mixed effects models for children. The magnitude of the AIC change indicates the impact on out-of-sample predictive accuracy, with more negative values reflecting greater importance of the variable."
) %>%
  kable_styling("striped", full_width = T) %>%
  column_spec(1, width_min = "1.5in")
Change of the Akaike Information Criterion (AIC) after addition of a given variable to mixed effects models for children. The magnitude of the AIC change indicates the impact on out-of-sample predictive accuracy, with more negative values reflecting greater importance of the variable.
outcome estimand AGE SEX BREASTFEED OTHER COV DIET FAMILY
MASS_Perc mean NA 1.9 2.9 -25.4 1.8 -4.8
HEIGHT_Perc mean NA 2.3 5.1 -18.5 1.2 -9.2
BMI_PERC mean NA 0.8 4.4 -9.6 4.2 -2.6
M_per_H_PERC mean NA 1.5 4.6 -9.4 4.1 -1.1
GLY mean 2.3 2.0 5.7 NA 4.4 1.7
TC mean 1.8 2.3 3.2 NA -0.6 -1.1
HDL mean -5.9 -1.4 3.5 NA 3.8 -6.7
LDL mean 2.3 1.9 2.6 NA -3.1 -1.3
log2_TG mean -2.4 -0.8 3.6 NA 0.7 -5.4
Ca mean -25.0 -1.0 5.3 NA 3.8 -4.1
P mean -33.7 2.3 3.3 NA 3.7 -5.0
Mg mean -0.4 1.5 0.2 1.5 0.7 0.1
log2_Se mean 2.3 2.3 2.3 -2.4 2.8 -12.3
log2_Zn mean 0.8 2.3 6.3 2.2 2.2 -1.9
FE mean -6.8 2.0 6.3 2.3 2.6 -1.8
TIBC mean 0.4 2.1 -2.6 -2.5 0.5 -4.2
log2_FERR mean -13.3 1.7 -12.1 -1.3 2.5 -1.9
TRF mean 0.8 1.8 -3.0 -2.9 0.0 -3.4
SATTRF mean -7.5 1.9 5.0 1.7 1.3 -1.1
log2_TRFINDEX mean -19.8 1.2 -12.1 -2.5 4.2 -0.2
log2_STRF mean -18.3 1.8 -7.7 -3.4 3.8 -0.3
HGB mean -14.0 1.5 4.7 NA 4.2 -4.1
MCV mean -20.5 0.8 -0.6 NA -16.7 -4.8
PTH mean -0.6 1.6 2.8 NA 2.5 -8.3
log2_CTX mean 1.8 2.0 1.9 3.0 4.1 1.6
P1NP mean -77.2 -1.4 0.7 4.1 3.0 0.4
log2_UI mean 2.3 -0.2 5.0 NA -5.7 0.6
UREA mean -10.4 -0.3 -0.8 NA -4.3 -6.5
log2_CREA mean -84.9 1.0 2.1 NA 1.4 0.8
UA mean 2.2 2.1 3.3 NA 4.0 -6.1
log2_HOLOTC mean -6.0 1.7 1.4 -7.7 -5.5 -13.1
log2_HCY mean 2.2 -4.9 -13.0 -0.4 1.3 -7.0
log2_MMA mean 0.8 1.4 -14.6 0.0 -3.6 -3.5
log2_VIT_D mean -1.9 2.3 4.7 -0.1 -1.6 -19.7
FOLATE mean -11.4 2.3 5.4 1.6 -15.7 -6.0
log2_IGF1 mean -60.6 2.3 0.7 NA 2.9 -6.9
log2_Ur_Ca mean -0.6 1.3 -0.3 1.7 2.7 -6.7
log2_uCCR mean 2.3 2.4 -1.1 2.1 4.5 -3.8
log2_uICR mean -1.3 2.1 1.2 NA 1.8 2.2
uPCR mean -17.6 2.3 -18.2 NA -13.3 -2.2

2.1.2.2 Heatmap

Open code 
heat_map(AIC_child_all, mixef = TRUE)

Heatmap showing decrease in the Akaike Information Criterion (AIC) after the addition of a given variable to a mixed-effects model in children. The AIC change reflects the impact of the variable on out-of-sample predictive accuracy, with values bounded between -30 (high importance, shown in blue) and 0 (no or negligible importance, shown in red). Black indicates that the variable was not included in the model (e.g., inapplicable supplementation-related variables)

2.2 Adults

2.2.1 Diet effects

2.2.1.1 Unstandardized effects

Open code 
diet_adult <- read.xlsx("gitignore/data/diet_adult_mixeff.xlsx") %>%
  mutate(outcome = clean_outcome_names(outcome)) %>%
    mutate(outcome = case_when(
      outcome == "FOLAT" ~ "FOLATE",
      outcome == "log2_VKFE" ~ "log2_TIBC",
      outcome == "Cros" ~ "CTX",
      outcome == "log2_Cros" ~ "log2_CTX",
      outcome == "VIT_AKTB12" ~ "HOLOTC",
      outcome == "log2_VIT_AKTB12" ~ "log2_HOLOTC",
      
      outcome == "log2_Ur_Ca" ~ "log2_uCa",
      outcome == "log2_Ca_per_Krea" ~ "log2_uCCR",
      outcome == "log2_I_per_Krea" ~ "log2_uICR",
      outcome == "P_per_Krea" ~ "uPCR",
      outcome == "log2_P_per_Krea" ~ "log2_uPCR",
      TRUE ~ outcome
    ))

  diet_adult2 <- diet_adult %>%
    mutate(across(
      .cols = c(VN_OM_diff, VG_OM_diff, VN_VG_diff),
      ~ round(.x, 2)
    ))

  kbl(diet_adult2,
    caption =
      "Differences in clinical outcomes (`outcome`) between diet groups in adults, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`)."
  ) %>%
    kable_styling("striped", full_width = T) %>%
    add_header_above(c(
      " " = 2, "(A) VN vs OM" = 3,
      "(B) VG vs OM" = 3, "(C) VN vs VG" = 3
    )) %>%
    column_spec(1, width_min = "1.5in")
Differences in clinical outcomes (`outcome`) between diet groups in adults, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
BMI mean -1.32 0.02339 0.08455 -0.88 0.18962 0.49849 -0.44 0.46534 0.68346
WHR mean -0.01 0.14670 0.30981 0.00 0.98174 0.98786 -0.01 0.15240 0.41124
HEIGHT mean -0.01 0.57922 0.71114 0.00 0.77318 0.92348 -0.01 0.38939 0.60276
MASS mean -3.90 0.05222 0.15554 -1.34 0.56242 0.75241 -2.56 0.21763 0.44840
HG mean 0.52 0.67431 0.73703 1.06 0.45491 0.71270 -0.54 0.66808 0.87222
PFAT mean -1.80 0.13113 0.29349 -0.02 0.98786 0.98786 -1.78 0.14847 0.41124
log2_FAT mean -0.20 0.08307 0.20548 -0.01 0.92471 0.98786 -0.19 0.11576 0.41124
log2_FFM mean -0.06 0.03242 0.10884 -0.02 0.47136 0.71465 -0.03 0.20516 0.44840
SBP mean -4.13 0.11534 0.27104 -3.15 0.29789 0.60872 -0.98 0.71735 0.88481
DBP mean -2.32 0.15161 0.30981 1.96 0.29481 0.60872 -4.28 0.01050 0.09871
GLY mean -0.09 0.23697 0.40141 -0.07 0.37425 0.65147 -0.01 0.87870 0.95109
log2_TC mean -0.21 0.00003 0.00026 -0.15 0.01053 0.09898 -0.06 0.24578 0.44840
HDL mean -0.04 0.55224 0.70150 -0.07 0.31978 0.61085 0.03 0.59208 0.81847
LDL mean -0.50 0.00003 0.00026 -0.32 0.02052 0.12406 -0.18 0.14232 0.41124
log2_TG mean -0.12 0.27029 0.43806 -0.10 0.43165 0.71270 -0.02 0.85003 0.95109
Ca mean 0.01 0.63437 0.71256 0.01 0.60664 0.77060 0.00 0.90874 0.95109
P mean -0.03 0.20320 0.38201 0.00 0.96848 0.98786 -0.03 0.23478 0.44840
Mg mean 0.00 0.83778 0.87502 0.00 0.88936 0.98786 0.00 0.72386 0.88481
log2_Se mean -0.11 0.30324 0.47507 -0.16 0.19091 0.49849 0.05 0.63508 0.85282
log2_Zn mean -0.22 0.00001 0.00010 -0.13 0.02053 0.12406 -0.09 0.07562 0.35541
FE mean 0.18 0.89367 0.89367 -2.49 0.10966 0.37364 2.67 0.05666 0.33579
log2_TIBC mean 0.03 0.32440 0.49183 0.04 0.32492 0.61085 0.00 0.89083 0.95109
log2_FERR mean -0.74 0.00000 0.00010 -0.55 0.00340 0.04000 -0.20 0.23528 0.44840
log2_TRF mean 0.03 0.35077 0.50011 0.03 0.35959 0.65002 0.00 0.91062 0.95109
SATTRF mean -0.36 0.86237 0.88111 -4.20 0.07837 0.33487 3.84 0.07404 0.35541
log2_TRFINDEX mean 0.25 0.00312 0.01466 0.30 0.00189 0.02955 -0.05 0.55376 0.78868
log2_STRF mean 0.04 0.47448 0.63715 0.10 0.11130 0.37364 -0.06 0.28159 0.45637
HGB mean -1.09 0.52297 0.68277 -1.10 0.57631 0.75241 0.01 0.99550 0.99742
MCV mean 1.66 0.00706 0.03018 -0.82 0.25060 0.56088 2.48 0.00010 0.00467
PTH mean 0.37 0.05295 0.15554 -0.15 0.49051 0.71697 0.52 0.00818 0.09615
log2_CTX mean 0.05 0.71938 0.76843 -0.10 0.53386 0.73798 0.15 0.27323 0.45637
log2_P1NP mean 0.14 0.23914 0.40141 0.01 0.95007 0.98786 0.13 0.25759 0.44840
log2_UI mean -0.41 0.16036 0.31403 -0.51 0.13466 0.39556 0.10 0.73421 0.88481
log2_UREA mean -0.29 0.00002 0.00022 -0.10 0.20526 0.50774 -0.19 0.00627 0.09615
log2_CREA mean -0.13 0.00096 0.00562 -0.07 0.11007 0.37364 -0.06 0.15750 0.41124
UA mean -5.57 0.59209 0.71114 -8.03 0.50340 0.71697 2.46 0.81915 0.95109
log2_HOLOTC mean -0.09 0.60522 0.71114 -0.30 0.07083 0.33487 0.21 0.13576 0.41124
log2_HCY mean -0.04 0.63675 0.71256 0.11 0.23261 0.54663 -0.15 0.05716 0.33579
log2_MMA mean 0.14 0.22342 0.40141 0.41 0.00026 0.01229 -0.26 0.00521 0.09615
log2_VIT_D mean 0.24 0.00226 0.01178 0.14 0.13009 0.39556 0.10 0.18156 0.43844
FOLATE mean 3.27 0.00019 0.00125 3.27 0.00118 0.02761 0.00 0.99742 0.99742
log2_uCa mean -0.51 0.06027 0.16663 -0.15 0.63420 0.78441 -0.37 0.18657 0.43844
log2_uCCR mean -0.16 0.36178 0.50011 -0.01 0.97077 0.98786 -0.16 0.39756 0.60276
log2_uICR mean -0.29 0.35971 0.50011 -0.67 0.07375 0.33487 0.37 0.25137 0.44840
uPCR mean -0.62 0.00000 0.00010 -0.35 0.02112 0.12406 -0.28 0.04185 0.32786
AL_adult log(OR) -1.95 0.01704 0.06673 -0.61 0.44626 0.71270 -1.34 0.08705 0.37193
BREAKS log(OR) -0.69 0.07545 0.19700 -0.12 0.78594 0.92348 -0.57 0.15320 0.41124

2.2.1.2 Standardized effects by 1SD

Open code 
meds <- diet_adult %>% filter(estimand == "mean")

ids <- nrow(meds)
sds <- vector("double", ids)
sds 
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [39] 0 0 0 0 0 0 0

before_outcome_col <- which(colnames(dat_adult) == 'aBMI')-1

for (i in 1:ids) {
  sds[i] <- ifelse(grepl("log2", diet_adult[(i), 1]),
    sd(log2(dat_adult[, (i + before_outcome_col)]), na.rm = TRUE),
    sd(     dat_adult[, (i + before_outcome_col)] , na.rm = TRUE)
  )
}

sds[(ids+1) : (ids+2)] <- 1

diet_adult_Z <- diet_adult %>%
  mutate(SD = sds) %>%
  mutate(
    VN_OM_diff = VN_OM_diff / SD,
    VG_OM_diff = VG_OM_diff / SD,
    VN_VG_diff = VN_VG_diff / SD
  ) %>%
  select(-SD) %>%
  mutate(VG_OM_P = if_else(VG_OM_P < 1e-5, 1e-5, VG_OM_P)) %>%
  mutate(across(.cols = c(VN_OM_diff, VG_OM_diff, VN_VG_diff), ~ round(.x, 2)))


kbl(diet_adult_Z,
  caption =
    "Standardized differences in clinical outcomes (`outcome`) between diet groups in children, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`)."
) %>%
  kable_styling("striped", full_width = T) %>%
  add_header_above(c(
    " " = 2, "(A) VN vs OM" = 3,
    "(B) VG vs OM" = 3, "(C) VN vs VG" = 3
  )) %>%
  column_spec(1, width_min = "1.5in")
Standardized differences in clinical outcomes (`outcome`) between diet groups in children, estimated using robust linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
BMI mean -0.37 0.02339 0.08455 -0.25 0.18962 0.49849 -0.12 0.46534 0.68346
WHR mean -0.20 0.14670 0.30981 0.00 0.98174 0.98786 -0.21 0.15240 0.41124
HEIGHT mean -0.09 0.57922 0.71114 0.05 0.77318 0.92348 -0.14 0.38939 0.60276
MASS mean -0.28 0.05222 0.15554 -0.10 0.56242 0.75241 -0.18 0.21763 0.44840
HG mean 0.04 0.67431 0.73703 0.08 0.45491 0.71270 -0.04 0.66808 0.87222
PFAT mean -0.22 0.13113 0.29349 0.00 0.98786 0.98786 -0.22 0.14847 0.41124
log2_FAT mean -0.29 0.08307 0.20548 -0.02 0.92471 0.98786 -0.27 0.11576 0.41124
log2_FFM mean -0.19 0.03242 0.10884 -0.08 0.47136 0.71465 -0.12 0.20516 0.44840
SBP mean -0.24 0.11534 0.27104 -0.19 0.29789 0.60872 -0.06 0.71735 0.88481
DBP mean -0.25 0.15161 0.30981 0.21 0.29481 0.60872 -0.46 0.01050 0.09871
GLY mean -0.17 0.23697 0.40141 -0.15 0.37425 0.65147 -0.02 0.87870 0.95109
log2_TC mean -0.71 0.00003 0.00026 -0.51 0.01053 0.09898 -0.21 0.24578 0.44840
HDL mean -0.09 0.55224 0.70150 -0.18 0.31978 0.61085 0.08 0.59208 0.81847
LDL mean -0.70 0.00003 0.00026 -0.45 0.02052 0.12406 -0.25 0.14232 0.41124
log2_TG mean -0.17 0.27029 0.43806 -0.14 0.43165 0.71270 -0.03 0.85003 0.95109
Ca mean 0.08 0.63437 0.71256 0.10 0.60664 0.77060 -0.02 0.90874 0.95109
P mean -0.20 0.20320 0.38201 -0.01 0.96848 0.98786 -0.19 0.23478 0.44840
Mg mean 0.03 0.83778 0.87502 -0.03 0.88936 0.98786 0.06 0.72386 0.88481
log2_Se mean -0.22 0.30324 0.47507 -0.33 0.19091 0.49849 0.11 0.63508 0.85282
log2_Zn mean -0.81 0.00001 0.00010 -0.48 0.02053 0.12406 -0.33 0.07562 0.35541
FE mean 0.02 0.89367 0.89367 -0.31 0.10966 0.37364 0.33 0.05666 0.33579
log2_TIBC mean 0.16 0.32440 0.49183 0.18 0.32492 0.61085 -0.02 0.89083 0.95109
log2_FERR mean -0.60 0.00000 0.00010 -0.44 0.00340 0.04000 -0.16 0.23528 0.44840
log2_TRF mean 0.15 0.35077 0.50011 0.17 0.35959 0.65002 -0.02 0.91062 0.95109
SATTRF mean -0.03 0.86237 0.88111 -0.33 0.07837 0.33487 0.30 0.07404 0.35541
log2_TRFINDEX mean 0.47 0.00312 0.01466 0.56 0.00189 0.02955 -0.10 0.55376 0.78868
log2_STRF mean 0.09 0.47448 0.63715 0.24 0.11130 0.37364 -0.15 0.28159 0.45637
HGB mean -0.07 0.52297 0.68277 -0.07 0.57631 0.75241 0.00 0.99550 0.99742
MCV mean 0.39 0.00706 0.03018 -0.19 0.25060 0.56088 0.58 0.00010 0.00467
PTH mean 0.31 0.05295 0.15554 -0.13 0.49051 0.71697 0.44 0.00818 0.09615
log2_CTX mean 0.06 0.71938 0.76843 -0.13 0.53386 0.73798 0.19 0.27323 0.45637
log2_P1NP mean 0.23 0.23914 0.40141 0.01 0.95007 0.98786 0.22 0.25759 0.44840
log2_UI mean -0.26 0.16036 0.31403 -0.32 0.13466 0.39556 0.06 0.73421 0.88481
log2_UREA mean -0.69 0.00002 0.00022 -0.24 0.20526 0.50774 -0.46 0.00627 0.09615
log2_CREA mean -0.45 0.00096 0.00562 -0.25 0.11007 0.37364 -0.20 0.15750 0.41124
UA mean -0.08 0.59209 0.71114 -0.11 0.50340 0.71697 0.03 0.81915 0.95109
log2_HOLOTC mean -0.13 0.60522 0.71114 -0.43 0.07083 0.33487 0.31 0.13576 0.41124
log2_HCY mean -0.08 0.63675 0.71256 0.21 0.23261 0.54663 -0.29 0.05716 0.33579
log2_MMA mean 0.21 0.22342 0.40141 0.60 0.00026 0.01229 -0.39 0.00521 0.09615
log2_VIT_D mean 0.55 0.00226 0.01178 0.32 0.13009 0.39556 0.24 0.18156 0.43844
FOLATE mean 0.67 0.00019 0.00125 0.67 0.00118 0.02761 0.00 0.99742 0.99742
log2_uCa mean -0.39 0.06027 0.16663 -0.11 0.63420 0.78441 -0.28 0.18657 0.43844
log2_uCCR mean -0.16 0.36178 0.50011 -0.01 0.97077 0.98786 -0.15 0.39756 0.60276
log2_uICR mean -0.18 0.35971 0.50011 -0.41 0.07375 0.33487 0.23 0.25137 0.44840
uPCR mean -0.76 0.00000 0.00010 -0.42 0.02112 0.12406 -0.33 0.04185 0.32786
AL_adult log(OR) -1.95 0.01704 0.06673 -0.61 0.44626 0.71270 -1.34 0.08705 0.37193
BREAKS log(OR) -0.69 0.07545 0.19700 -0.12 0.78594 0.92348 -0.57 0.15320 0.41124

2.2.1.3 Vulcano plots

VN vs OM

Open code 

plotac <- "fig31"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_adult_Z %>%
    filter(estimand != "log(OR)") %>%
    
    mutate(
      pv_adj = VN_OM_P_adj,
      pv_raw = VN_OM_P,
      b_coe = VN_OM_diff
    ) %>%
    mutate(
      pv_adj = VN_OM_P_adj,
      pv_raw = VN_OM_P,
      b_coe = VN_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VN+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.7,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs. OM, adults") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VG vs OM

Open code 
plotac <- "fig32"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_adult_Z %>%
    filter(estimand != "log(OR)") %>%
    
    mutate(
      pv_adj = VG_OM_P_adj,
      pv_raw = VG_OM_P,
      b_coe = VG_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VG+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VG+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.6,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VG vs OM, adults") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VN vs VG

Open code 

plotac <- "fig33"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_adult_Z %>%
    filter(estimand != "log(OR)") %>%
    mutate(
      pv_adj = VN_VG_P_adj,
      pv_raw = VN_VG_P,
      b_coe = VN_VG_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "VG+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c(
          "VG+",
          "VN+",
          "P>0.05"
        )
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.6,
      size = 3.1,
      seed = 16
    ) +
    labs(
      x = "Standardized difference",
      y = "-log10 (P-value)"
    ) +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs VG, adults") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

2.2.2 Covariates importance

2.2.2.1 Table

Open code 

AIC_adult <- read.xlsx("gitignore/data/AIC_adult_mixeff.xlsx") %>%
  mutate(across(where(is.numeric), ~ round(., 1))) %>%
  mutate(outcome = clean_outcome_names(outcome)) %>%
  filter(!outcome %in% c("AL_adult", "BREAKS")) %>%
  rename(
    SEX = aSEX,
    `OTHER COV` = other_cov,
    FAMILY = family,
    DIET = diet,
    AGE = aAGE
  ) %>%
  mutate(outcome = case_when(
    outcome == "FOLAT" ~ "FOLATE",
    outcome == "VKFE" ~ "TIBC",
    outcome == "Cros" ~ "CTX",
    outcome == "log2_Cros" ~ "log2_CTX",
    outcome == "VIT_AKTB12" ~ "HOLOTC",
    outcome == "log2_VIT_AKTB12" ~ "log2_HOLOTC",
    
    outcome == 'log2_aUr_Ca' ~ 'log2_uCa',
    outcome == 'log2_Ca_per_Krea' ~ 'log2_uCCR',
    outcome == 'log2_I_per_Krea' ~ 'log2_uICR',
    outcome == 'P_per_Krea' ~ 'uPCR',
    outcome == 'log2_P_per_Krea' ~ 'log2_uPCR',
    TRUE ~ outcome
  ))

kbl(AIC_adult,
  caption =
    "Change of the Akaike Information Criterion (AIC) after addition of a given variable to mixed effects models for adults. The magnitude of the AIC change indicates the impact on out-of-sample predictive accuracy, with more negative values reflecting greater importance of the variable."
) %>%
  kable_styling("striped", full_width = T) %>%
  column_spec(1, width_min = "1.5in")
Change of the Akaike Information Criterion (AIC) after addition of a given variable to mixed effects models for adults. The magnitude of the AIC change indicates the impact on out-of-sample predictive accuracy, with more negative values reflecting greater importance of the variable.
outcome estimand AGE SEX OTHER COV DIET FAMILY
BMI mean 1.6 -11.9 NA 0.0 0.4
WHR mean -2.7 -79.2 NA 0.6 1.5
HEIGHT mean 1.6 -104.7 NA 3.4 -5.3
MASS mean 1.0 -71.7 NA -0.9 -1.3
HG mean 0.8 -228.4 NA 3.3 -2.0
PFAT mean 0.0 -90.2 NA 1.9 -0.6
log2_FAT mean 0.1 -18.7 NA 0.5 0.5
log2_FFM mean 2.1 -249.1 NA -1.2 2.0
SBP mean 1.5 -47.3 NA 0.8 1.7
DBP mean 1.2 -6.8 NA -4.4 1.6
GLY mean 1.9 -3.9 NA 2.6 1.5
log2_TC mean -1.3 -2.7 NA -11.5 2.0
HDL mean -2.4 -66.5 NA 3.1 -1.1
LDL mean 1.3 1.8 NA -9.8 1.6
log2_TG mean 1.8 -21.2 NA 2.8 -1.7
Ca mean 1.2 -19.7 NA 4.0 -0.3
P mean 1.6 -41.2 NA 2.4 1.9
Mg mean 2.0 1.2 2.0 3.4 2.2
log2_Se mean 2.1 0.4 0.0 2.1 -37.2
log2_Zn mean 1.1 -3.0 0.7 -8.6 -24.0
FE mean 1.9 -7.8 1.3 3.3 0.9
log2_VKFE mean 1.7 -13.0 2.1 2.6 -2.4
log2_FERR mean 1.8 -72.1 1.3 -14.0 -4.3
log2_TRF mean 1.6 -12.8 2.2 2.7 -2.7
SATTRF mean 1.8 -14.9 1.6 3.5 0.2
log2_TRFINDEX mean 2.1 -41.2 -0.7 -6.1 -1.4
log2_STRF mean 1.9 2.2 -1.1 2.3 2.1
HGB mean 0.5 -71.0 NA 3.2 0.5
MCV mean 1.1 -5.1 NA -2.8 2.1
PTH mean 0.9 1.0 NA -5.2 -0.4
log2_CTX mean -8.7 1.7 2.9 3.1 2.2
log2_P1NP mean -6.3 2.2 0.4 1.6 2.2
log2_UI mean -1.0 1.2 NA 1.4 -44.4
log2_UREA mean 0.4 -16.3 NA -12.3 0.3
log2_CREA mean -0.1 -74.7 NA -4.8 0.5
UA mean -2.8 -87.8 NA 3.7 1.3
log2_HOLOTC mean 2.1 -7.6 -0.4 -0.4 -13.3
log2_HCY mean 2.2 -55.6 2.2 2.2 -1.6
log2_MMA mean -0.6 -1.7 -9.9 -5.1 0.9
log2_VIT_D mean 1.4 2.1 0.7 -3.3 -6.4
FOLATE mean 1.2 -14.8 -2.7 -8.7 -12.5
log2_Ur_Ca mean 1.5 -14.4 1.4 0.3 -8.5
log2_uCCR mean 0.7 1.7 1.6 2.9 0.4
log2_uICR mean -3.5 -12.4 NA 1.6 -14.9
uPCR mean 2.0 -14.5 NA -13.5 -3.6

2.2.2.2 Heatmap

Open code 
heat_map(AIC_adult, mixef = TRUE)

Heatmap showing decrease in the Akaike Information Criterion (AIC) after the addition of a given variable to a mixed-effects model in adults. The AIC change reflects the impact of the variable on out-of-sample predictive accuracy, with values bounded between -30 (high importance, shown in blue) and 0 (no or negligible importance, shown in red). Black indicates that the variable was not included in the model (e.g., inapplicable supplementation-related variables)

2.3 Vulcano matrix

Prepare

Open code 
plotac <- "fig35"
path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  fig19 <- fig19 + theme(
    legend.position = "none",
    plot.title = element_blank()
  )
  fig20 <- fig20 + theme(
    legend.position = "none",
    plot.title = element_blank()
  )
  fig21 <- fig21 + theme(
    legend.position = "none",
    plot.title = element_blank()
  )

  fig31 <- fig31 + theme(
    legend.position = "bottom",
    plot.title = element_blank()
  )
  fig32 <- fig32 + theme(
    legend.position = "bottom",
    plot.title = element_blank()
  )
  fig33 <- fig33 + theme(
    legend.position = "bottom",
    plot.title = element_blank()
  )
}

Create subfigures

Open code 

if (file.exists(path) == FALSE) {
  children_volc <- ggarrange(
    fig19 +
      theme(
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank()
      ),
    fig20 +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank()
      ),
    fig21 +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank()
      ),
    nrow = 1,
    ncol = 3,
    widths = c(1.1, 1, 1),
    labels = c("A", "", "")
  )

  children_volc <- annotate_figure(
    children_volc,
    top = text_grob("Children", face = "bold", size = 10)
  )


  adult_volc <- ggarrange(
    fig31,
    fig32 +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank()
      ),
    fig33 +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank()
      ),
    nrow = 1,
    ncol = 3,
    widths = c(1.1, 1, 1),
    labels = c("B", "", "")
  )

  adult_volc <- annotate_figure(
    adult_volc,
    top = text_grob("Adults", face = "bold", size = 10)
  )
}

2.3.1 Print

Open code 
if (file.exists(path) == FALSE) {
  assign(
    plotac,
    ggarrange(
      children_volc,
      adult_volc,
      nrow = 2,
      ncol = 1,
      heights = c(1, 1.2)
    )
  )

  saveRDS(get(plotac), file = path)
  get(plotac)
} else {
  assign(plotac, readRDS(path))
  get(plotac)
}

Vulcano plot showing standardized differences between vegan (VN) and omnivores (OM) (left), vegetarians (VG) vs OM (middle) and VN vs VG (right) across clinical outcomes in children (A, top) and adults (B, bottom). Estimated effects are based on robust linear mixed-effects models. P-values are not corrected for multiple comparisons to maximize sensitivity for a potential risks
Open code 

## save svg
plotac <- "fig35"
path <- paste0("gitignore/figures/", plotac, ".svg")

if(!file.exists(path)){
  ggsave(path, get(plotac), device = 'svg', width = 10, height = 7)
}

2.4 Heatmap matrix

Open code 
h1 <- heat_map(AIC_child_all, mixef = TRUE, main = 'Children', listt = TRUE)
h2 <- heat_map(AIC_adult, mixef = TRUE, main = 'Adults', listt = TRUE)

grid.arrange(grobs = list(h1$gtable, h2$gtable),
             nrow = 1, ncol = 2)

Heatmap showing decrease in the Akaike Information Criterion (AIC) after the addition of a given variable to a mixed-effects model, separately for children (left) and adults (right). The AIC change reflects the impact of the variable on out-of-sample predictive accuracy, with values bounded between -30 (high importance, shown in blue) and 0 (no or negligible importance, shown in red). Black indicates that the variable was not included in the model (e.g., inapplicable supplementation-related variables)
Open code 

plotac <- 'heatmaps_mixef_combo'
path = paste0('gitignore/figures/', plotac, '.svg')
if(file.exists(path) == FALSE){
  svg(path, width = 7, height = 6)  
  grid.arrange(grobs = list(h1$gtable, h2$gtable),
             nrow = 1, ncol = 2)
  dev.off()
}

path <- paste0('gitignore/figures/', plotac, '.pdf')
if (!file.exists(path)) {
  pdf(path, width = 7, height = 6)  
  grid.arrange(grobs = list(h1$gtable, h2$gtable),
             nrow = 1, ncol = 2)
  dev.off()
}

3 Linear-mixed-effects models (non-robust variant)

3.1 Children

3.1.1 Diet effects

3.1.1.1 Unstandardized effects

Open code 
diet_child_all <- read.xlsx("gitignore/data/diet_child_all_mixeff_non_robust.xlsx") %>% 
  mutate(outcome = clean_outcome_names(outcome)) %>%
  mutate(outcome = case_when(
    outcome == 'FOLAT' ~ 'FOLATE',
    outcome == 'VKFE' ~ 'TIBC',
    outcome == 'Cros' ~ 'CTX',
    outcome == 'log2_Cros' ~ 'log2_CTX',
    outcome == 'VIT_AKTB12' ~ 'HOLOTC',
    outcome == 'log2_VIT_AKTB12' ~ 'log2_HOLOTC',
    
    outcome == 'log2_aUr_Ca' ~ 'log2_uCa',
    outcome == 'log2_Ca_per_Krea' ~ 'log2_uCCR',
    outcome == 'log2_I_per_Krea' ~ 'log2_uICR',
    outcome == 'P_per_Krea' ~ 'uPCR',
    outcome == 'log2_P_per_Krea' ~ 'log2_uPCR',
    
    TRUE ~ outcome
  ))

diet_child_all2 <- diet_child_all %>%   
  mutate(across(.cols = c(VN_OM_diff, VG_OM_diff, VN_VG_diff), ~ round(.x, 2)))


kableExtra::kable(diet_child_all2, 
                caption = "Differences in clinical outcomes (`outcome`) between diet groups in children, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`)."
                ) %>% 
  kable_styling("striped",full_width = T) %>% 
  add_header_above(c(" " = 2, "(A) VN vs OM" = 3, 
                     "(B) VG vs OM" = 3, "(C) VN vs VG" =3)) %>% 
  column_spec(1, width_min = '1.5in')
Differences in clinical outcomes (`outcome`) between diet groups in children, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
MASS_Perc mean -9.03 0.1224106 0.3139879 -3.48 0.6045777 0.9971723 -5.55 0.3552860 0.7105720
HEIGHT_Perc mean -8.70 0.1412946 0.3139879 0.02 0.9971723 0.9971723 -8.72 0.1520505 0.5320605
BMI_PERC mean -2.51 0.6471659 0.7396182 -2.93 0.6428129 0.9971723 0.42 0.9402838 0.9754998
M_per_H_PERC mean -3.05 0.5724270 0.7155338 -3.39 0.5854662 0.9971723 0.34 0.9511123 0.9754998
GLY mean -0.03 0.8000192 0.8505275 0.02 0.9175935 0.9971723 -0.05 0.7323108 0.8876495
TC mean -0.33 0.0298531 0.1492653 -0.20 0.2682746 0.9971723 -0.14 0.4077395 0.7766466
HDL mean -0.02 0.8080011 0.8505275 -0.07 0.4158826 0.9971723 0.06 0.5063240 0.8623993
LDL mean -0.30 0.0206126 0.1177862 -0.01 0.9715679 0.9971723 -0.30 0.0350391 0.2803128
log2_TG mean 0.00 0.9937607 0.9937607 -0.31 0.0985864 0.9971723 0.31 0.0737402 0.3910889
Ca mean -0.01 0.6338424 0.7396182 0.01 0.7469291 0.9971723 -0.02 0.4273929 0.7770781
P mean -0.03 0.4078436 0.5826337 -0.01 0.8957263 0.9971723 -0.02 0.5344134 0.8623993
Mg mean 0.03 0.0789224 0.2630747 0.03 0.1282706 0.9971723 0.00 0.9760101 0.9760101
log2_Se mean -0.15 0.2469283 0.4489605 -0.01 0.9547127 0.9971723 -0.14 0.3180938 0.7068751
log2_Zn mean -0.10 0.1342529 0.3139879 -0.07 0.4099470 0.9971723 -0.04 0.6354120 0.8876495
FE mean 1.47 0.3367963 0.4989575 -0.61 0.7349709 0.9971723 2.08 0.2126189 0.5669838
TIBC mean -3.77 0.0728617 0.2630747 -0.69 0.7784954 0.9971723 -3.09 0.1729197 0.5320605
log2_FERR mean -0.23 0.2216702 0.4489605 -0.26 0.2368430 0.9971723 0.03 0.8756598 0.9729554
TRF mean -0.16 0.0512531 0.2050123 -0.04 0.6850538 0.9971723 -0.12 0.1712928 0.5320605
SATTRF mean 3.24 0.1341627 0.3139879 -0.01 0.9957130 0.9971723 3.25 0.1668047 0.5320605
log2_TRFINDEX mean 0.03 0.7547654 0.8386283 0.08 0.5273110 0.9971723 -0.05 0.6925916 0.8876495
log2_STRF mean -0.04 0.4314124 0.5950516 -0.01 0.9199755 0.9971723 -0.04 0.5389996 0.8623993
HGB mean 0.96 0.6462453 0.7396182 1.36 0.5857707 0.9971723 -0.40 0.8603510 0.9729554
MCV mean 3.18 0.0000784 0.0013022 0.24 0.7920987 0.9971723 2.94 0.0006757 0.0270265
PTH mean 0.32 0.3000010 0.4653449 -0.10 0.7844067 0.9971723 0.42 0.2112858 0.5669838
log2_CTX mean 0.06 0.4783561 0.6255021 0.05 0.6274840 0.9971723 0.01 0.9083010 0.9754998
P1NP mean 41.07 0.3021109 0.4653449 55.74 0.2460170 0.9971723 -14.67 0.7315742 0.8876495
log2_UI mean -0.96 0.0031127 0.0311274 -0.40 0.2663855 0.9971723 -0.56 0.0782178 0.3910889
UREA mean -0.57 0.0073351 0.0586807 -0.10 0.6769059 0.9971723 -0.47 0.0425533 0.2836889
log2_CREA mean -0.10 0.0887183 0.2729793 -0.04 0.5546942 0.9971723 -0.06 0.3507577 0.7105720
UA mean -7.85 0.4847641 0.6255021 -5.10 0.7005236 0.9971723 -2.75 0.8209653 0.9658415
log2_HOLOTC mean 0.27 0.2420910 0.4489605 -0.30 0.1955345 0.9971723 0.57 0.0021550 0.0431006
log2_HCY mean -0.17 0.2379736 0.4489605 0.02 0.8671099 0.9971723 -0.19 0.0969613 0.4309390
log2_MMA mean -0.51 0.0464757 0.2050123 -0.02 0.9508927 0.9971723 -0.49 0.0124297 0.1242974
log2_VIT_D mean 0.26 0.0171848 0.1145657 0.22 0.0893232 0.9971723 0.04 0.7185425 0.8876495
FOLATE mean 3.94 0.0000928 0.0013022 4.49 0.0001577 0.0063064 -0.55 0.6002629 0.8876495
log2_IGF1 mean -0.14 0.3024742 0.4653449 0.02 0.9130673 0.9971723 -0.16 0.2940446 0.6918698
log2_Ur_Ca mean -0.39 0.2951999 0.4653449 0.04 0.9254037 0.9971723 -0.43 0.2455164 0.6137910
log2_uCCR mean 0.02 0.9638138 0.9885270 0.15 0.7138363 0.9971723 -0.13 0.7103345 0.8876495
log2_uICR mean -0.47 0.1131462 0.3139879 -0.33 0.3135735 0.9971723 -0.14 0.6416148 0.8876495
uPCR mean -1.88 0.0000977 0.0013022 -0.53 0.3183705 0.9971723 -1.35 0.0041575 0.0554332

3.1.1.2 Standardized effects by 1SD

Open code 
meds <- diet_child_all %>% filter(estimand == "mean")

ids <- nrow(meds)
sds <- vector("double", ids)
sds 
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [39] 0 0

before_outcome_col <- which(colnames(dat_child_all) == 'aMASS_Perc')-1

sds <- vector("double", ids)

for(i in 1:ids){
  sds[i] <- ifelse(grepl("log2", diet_child_all[(i), 'outcome']),
                   sd(log2(dat_child_all[,(i + before_outcome_col)]), na.rm = TRUE),
                   sd(dat_child_all[,(i + before_outcome_col)], na.rm = TRUE))
}


diet_child_all_Z <- diet_child_all %>% 
  mutate(SD = sds) %>% 
  mutate(
    VN_OM_diff = VN_OM_diff/SD,
    VG_OM_diff = VG_OM_diff/SD,
    VN_VG_diff = VN_VG_diff/SD) %>% 
  select(-c(SD)) %>% 
  mutate(across(.cols = c(VG_OM_diff, VN_OM_diff, VN_VG_diff), ~ round(.x, 2)))


kbl(diet_child_all_Z, caption = 
      "Standardized differences in clinical outcomes (`outcome`) between diet groups in children, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).") %>% 
  kable_styling("striped",full_width = T) %>% 
  add_header_above(c(" " = 2, "(A) VN vs OM" = 3, 
                     "(B) VG vs OM" = 3, "(C) VN vs VG" =3)) %>% 
  column_spec(1, width_min = '1.5in')
Standardized differences in clinical outcomes (`outcome`) between diet groups in children, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
MASS_Perc mean -0.31 0.1224106 0.3139879 -0.12 0.6045777 0.9971723 -0.19 0.3552860 0.7105720
HEIGHT_Perc mean -0.32 0.1412946 0.3139879 0.00 0.9971723 0.9971723 -0.32 0.1520505 0.5320605
BMI_PERC mean -0.10 0.6471659 0.7396182 -0.12 0.6428129 0.9971723 0.02 0.9402838 0.9754998
M_per_H_PERC mean -0.12 0.5724270 0.7155338 -0.14 0.5854662 0.9971723 0.01 0.9511123 0.9754998
GLY mean -0.06 0.8000192 0.8505275 0.03 0.9175935 0.9971723 -0.08 0.7323108 0.8876495
TC mean -0.50 0.0298531 0.1492653 -0.30 0.2682746 0.9971723 -0.21 0.4077395 0.7766466
HDL mean -0.06 0.8080011 0.8505275 -0.22 0.4158826 0.9971723 0.17 0.5063240 0.8623993
LDL mean -0.54 0.0206126 0.1177862 -0.01 0.9715679 0.9971723 -0.53 0.0350391 0.2803128
log2_TG mean 0.00 0.9937607 0.9937607 -0.44 0.0985864 0.9971723 0.43 0.0737402 0.3910889
Ca mean -0.09 0.6338424 0.7396182 0.07 0.7469291 0.9971723 -0.17 0.4273929 0.7770781
P mean -0.16 0.4078436 0.5826337 -0.03 0.8957263 0.9971723 -0.13 0.5344134 0.8623993
Mg mean 0.38 0.0789224 0.2630747 0.39 0.1282706 0.9971723 -0.01 0.9760101 0.9760101
log2_Se mean -0.29 0.2469283 0.4489605 -0.02 0.9547127 0.9971723 -0.27 0.3180938 0.7068751
log2_Zn mean -0.36 0.1342529 0.3139879 -0.23 0.4099470 0.9971723 -0.13 0.6354120 0.8876495
FE mean 0.21 0.3367963 0.4989575 -0.09 0.7349709 0.9971723 0.30 0.2126189 0.5669838
TIBC mean -0.41 0.0728617 0.2630747 -0.08 0.7784954 0.9971723 -0.34 0.1729197 0.5320605
log2_FERR mean -0.26 0.2216702 0.4489605 -0.29 0.2368430 0.9971723 0.04 0.8756598 0.9729554
TRF mean -0.45 0.0512531 0.2050123 -0.11 0.6850538 0.9971723 -0.34 0.1712928 0.5320605
SATTRF mean 0.33 0.1341627 0.3139879 0.00 0.9957130 0.9971723 0.33 0.1668047 0.5320605
log2_TRFINDEX mean 0.06 0.7547654 0.8386283 0.15 0.5273110 0.9971723 -0.09 0.6925916 0.8876495
log2_STRF mean -0.16 0.4314124 0.5950516 -0.02 0.9199755 0.9971723 -0.14 0.5389996 0.8623993
HGB mean 0.10 0.6462453 0.7396182 0.14 0.5857707 0.9971723 -0.04 0.8603510 0.9729554
MCV mean 0.83 0.0000784 0.0013022 0.06 0.7920987 0.9971723 0.77 0.0006757 0.0270265
PTH mean 0.26 0.3000010 0.4653449 -0.08 0.7844067 0.9971723 0.34 0.2112858 0.5669838
log2_CTX mean 0.16 0.4783561 0.6255021 0.13 0.6274840 0.9971723 0.03 0.9083010 0.9754998
P1NP mean 0.14 0.3021109 0.4653449 0.19 0.2460170 0.9971723 -0.05 0.7315742 0.8876495
log2_UI mean -0.77 0.0031127 0.0311274 -0.32 0.2663855 0.9971723 -0.45 0.0782178 0.3910889
UREA mean -0.46 0.0073351 0.0586807 -0.08 0.6769059 0.9971723 -0.38 0.0425533 0.2836889
log2_CREA mean -0.24 0.0887183 0.2729793 -0.10 0.5546942 0.9971723 -0.14 0.3507577 0.7105720
UA mean -0.16 0.4847641 0.6255021 -0.11 0.7005236 0.9971723 -0.06 0.8209653 0.9658415
log2_HOLOTC mean 0.33 0.2420910 0.4489605 -0.37 0.1955345 0.9971723 0.70 0.0021550 0.0431006
log2_HCY mean -0.37 0.2379736 0.4489605 0.05 0.8671099 0.9971723 -0.42 0.0969613 0.4309390
log2_MMA mean -0.56 0.0464757 0.2050123 -0.02 0.9508927 0.9971723 -0.54 0.0124297 0.1242974
log2_VIT_D mean 0.55 0.0171848 0.1145657 0.47 0.0893232 0.9971723 0.08 0.7185425 0.8876495
FOLATE mean 0.85 0.0000928 0.0013022 0.96 0.0001577 0.0063064 -0.12 0.6002629 0.8876495
log2_IGF1 mean -0.17 0.3024742 0.4653449 0.02 0.9130673 0.9971723 -0.20 0.2940446 0.6918698
log2_Ur_Ca mean -0.27 0.2951999 0.4653449 0.03 0.9254037 0.9971723 -0.30 0.2455164 0.6137910
log2_uCCR mean 0.01 0.9638138 0.9885270 0.10 0.7138363 0.9971723 -0.09 0.7103345 0.8876495
log2_uICR mean -0.36 0.1131462 0.3139879 -0.26 0.3135735 0.9971723 -0.11 0.6416148 0.8876495
uPCR mean -0.88 0.0000977 0.0013022 -0.25 0.3183705 0.9971723 -0.63 0.0041575 0.0554332

3.1.1.3 Vulcano plots

VN vs OM

Open code 

plotac <- "fig19_nr"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_child_all_Z %>%
    filter(estimand != "log(OR)") %>%
    mutate(
      pv_adj = VN_OM_P_adj,
      pv_raw = VN_OM_P,
      b_coe = VN_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VN+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.3,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs. OM, all children") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VG vs OM

Open code 
plotac <- "fig20_nr"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_child_all_Z %>%
    filter(estimand != "log(OR)") %>%
    mutate(VG_OM_P = if_else(VG_OM_P < 1e-5, 1e-5, VG_OM_P)) %>%
    mutate(
      pv_adj = VG_OM_P_adj,
      pv_raw = VG_OM_P,
      b_coe = VG_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VG+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VG+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.3,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VG vs. OM, all children") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VN vs VG

Open code 
plotac <- "fig21_nr"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_child_all_Z %>%
    filter(estimand != "log(OR)") %>%
    
    mutate(
      pv_adj = VN_VG_P_adj,
      pv_raw = VN_VG_P,
      b_coe = VN_VG_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "VG+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c(
          "VG+",
          "VN+",
          "P>0.05"
        )
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.5,
      size = 3.1,
      seed = 16
    ) +
    labs(
      x = "Standardized difference",
      y = "-log10 (P-value)"
    ) +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs VG, all children") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

3.2 Adults

3.2.1 Diet effects

3.2.1.1 Unstandardized effects

Open code 
diet_adult <- read.xlsx("gitignore/data/diet_adult_non_robust_mixeff_non_robust.xlsx") %>%
  mutate(outcome = clean_outcome_names(outcome)) %>%
  mutate(outcome = case_when(
    outcome == "FOLAT" ~ "FOLATE",
    outcome == "log2_VKFE" ~ "log2_TIBC",
    outcome == "Cros" ~ "CTX",
    outcome == "log2_Cros" ~ "log2_CTX",
    outcome == "VIT_AKTB12" ~ "HOLOTC",
    outcome == "log2_VIT_AKTB12" ~ "log2_HOLOTC",
    
    outcome == 'log2_aUr_Ca' ~ 'log2_uCa',
    outcome == 'log2_Ca_per_Krea' ~ 'log2_uCCR',
    outcome == 'log2_I_per_Krea' ~ 'log2_uICR',
    outcome == 'P_per_Krea' ~ 'uPCR',
    outcome == 'log2_P_per_Krea' ~ 'log2_uPCR',
    TRUE ~ outcome
  ))

diet_adult2 <- diet_adult %>%
  mutate(across(.cols = c(VN_OM_diff, VG_OM_diff, VN_VG_diff), ~ round(.x, 2)))

kbl(diet_adult2,
  caption =
    "Differences in clinical outcomes (`outcome`) between diet groups in adults, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`)."
) %>%
  kable_styling("striped", full_width = T) %>%
  add_header_above(c(
    " " = 2, "(A) VN vs OM" = 3,
    "(B) VG vs OM" = 3, "(C) VN vs VG" = 3
  )) %>%
  column_spec(1, width_min = "1.5in")
Differences in clinical outcomes (`outcome`) between diet groups in adults, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
BMI mean -1.27 0.05644 0.17544 -0.38 0.62265 0.81153 -0.90 0.19057 0.42652
WHR mean -0.02 0.12747 0.27232 0.00 0.94638 0.98844 -0.02 0.12097 0.42652
HEIGHT mean -0.01 0.71233 0.77859 0.01 0.62498 0.81153 -0.01 0.36602 0.57343
MASS mean -4.45 0.04604 0.15455 -0.65 0.79770 0.89062 -3.80 0.09800 0.42652
HG mean 0.83 0.54439 0.62405 1.52 0.33940 0.79760 -0.68 0.62794 0.77666
PFAT mean -1.60 0.20782 0.34506 0.02 0.98770 0.99349 -1.63 0.21572 0.44534
log2_FAT mean -0.22 0.08157 0.19170 -0.05 0.74859 0.87960 -0.17 0.18161 0.42652
log2_FFM mean -0.06 0.03110 0.11687 -0.02 0.58978 0.81153 -0.04 0.13328 0.42652
SBP mean -5.01 0.06837 0.18663 -3.55 0.26104 0.72113 -1.46 0.60348 0.77666
DBP mean -2.64 0.11151 0.24958 2.24 0.24169 0.72113 -4.88 0.00498 0.08969
GLY mean -0.08 0.36725 0.49317 0.02 0.81483 0.89062 -0.11 0.25703 0.46464
log2_TC mean -0.20 0.00010 0.00151 -0.13 0.02523 0.23715 -0.07 0.16557 0.42652
HDL mean -0.04 0.56114 0.62794 -0.08 0.30203 0.74713 0.04 0.55223 0.77666
LDL mean -0.48 0.00027 0.00250 -0.24 0.10186 0.51319 -0.24 0.07216 0.42652
log2_TG mean -0.15 0.24651 0.37498 -0.10 0.53313 0.81153 -0.06 0.66848 0.78546
Ca mean 0.00 0.76142 0.81334 0.01 0.59798 0.81153 0.00 0.76706 0.85837
P mean -0.03 0.26930 0.38813 0.00 0.99349 0.99349 -0.03 0.28863 0.50243
Mg mean -0.01 0.51192 0.60897 -0.01 0.38038 0.81153 0.00 0.73404 0.84146
log2_Se mean -0.14 0.20119 0.34506 -0.16 0.20654 0.72113 0.02 0.85731 0.87595
log2_Zn mean -0.20 0.00052 0.00404 -0.10 0.12011 0.51319 -0.10 0.08854 0.42652
FE mean 0.19 0.89833 0.93826 -1.30 0.44987 0.81153 1.49 0.33540 0.54357
log2_TIBC mean 0.05 0.20427 0.34506 0.03 0.50558 0.81153 0.02 0.62114 0.77666
log2_FERR mean -0.86 0.00003 0.00107 -0.58 0.01227 0.14415 -0.28 0.17286 0.42652
log2_TRF mean 0.05 0.21291 0.34506 0.03 0.52717 0.81153 0.02 0.61144 0.77666
SATTRF mean 0.00 0.99979 0.99979 -2.03 0.45260 0.81153 2.03 0.40363 0.61195
log2_TRFINDEX mean 0.27 0.00385 0.02011 0.28 0.00898 0.14061 -0.01 0.91768 0.91768
log2_STRF mean 0.09 0.24733 0.37498 0.11 0.22009 0.72113 -0.02 0.80476 0.85963
HGB mean -2.40 0.31410 0.43420 -1.31 0.63448 0.81153 -1.09 0.65602 0.78546
MCV mean 1.05 0.16119 0.31567 -0.94 0.27485 0.72113 1.98 0.01093 0.12840
PTH mean 0.48 0.03233 0.11687 -0.17 0.51452 0.81153 0.65 0.00553 0.08969
log2_CTX mean 0.00 0.99695 0.99979 -0.15 0.38453 0.81153 0.15 0.31117 0.52232
log2_P1NP mean 0.15 0.16018 0.31567 0.02 0.84126 0.89862 0.13 0.22741 0.44534
log2_UI mean -0.49 0.19890 0.34506 -0.72 0.11270 0.51319 0.22 0.56837 0.77666
log2_UREA mean -0.28 0.00022 0.00250 -0.07 0.42385 0.81153 -0.21 0.00573 0.08969
log2_CREA mean -0.13 0.00349 0.02011 -0.06 0.20547 0.72113 -0.07 0.14378 0.42652
UA mean -7.83 0.45945 0.56826 -5.85 0.63171 0.81153 -1.98 0.85613 0.87595
log2_HOLOTC mean -0.15 0.41789 0.54558 -0.37 0.04226 0.33106 0.22 0.15627 0.42652
log2_HCY mean -0.09 0.43623 0.55413 0.05 0.63887 0.81153 -0.15 0.15720 0.42652
log2_MMA mean 0.30 0.05973 0.17544 0.46 0.00308 0.14061 -0.16 0.22122 0.44534
log2_VIT_D mean 0.24 0.00857 0.04027 0.11 0.27618 0.72113 0.13 0.15364 0.42652
FOLATE mean 3.29 0.00070 0.00468 3.03 0.00622 0.14061 0.26 0.79120 0.85963
log2_Ur_Ca mean -0.48 0.07264 0.18663 -0.11 0.70675 0.85960 -0.37 0.17532 0.42652
log2_uCCR mean -0.22 0.27251 0.38813 -0.08 0.71328 0.85960 -0.14 0.49930 0.73335
log2_uICR mean -0.23 0.51827 0.60897 -0.67 0.11644 0.51319 0.43 0.24063 0.45239
uPCR mean -0.65 0.00005 0.00107 -0.32 0.06937 0.46575 -0.33 0.03866 0.36338
AL_adult log(OR) -1.95 0.01704 0.07280 -0.61 0.44626 0.81153 -1.34 0.08705 0.42652
BREAKS log(OR) -0.69 0.07545 0.18663 -0.12 0.78594 0.89062 -0.57 0.15320 0.42652

3.2.1.2 Standardized effects by 1SD

Open code 
meds <- diet_adult %>% filter(estimand == "mean")

ids <- nrow(meds)
sds <- vector("double", ids)
sds 
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [39] 0 0 0 0 0 0 0

before_outcome_col <- which(colnames(dat_adult) == 'aBMI')-1

for (i in 1:ids) {
  sds[i] <- ifelse(grepl("log2", diet_adult[(i), 1]),
    sd(log2(dat_adult[, (i + before_outcome_col)]), na.rm = TRUE),
    sd(     dat_adult[, (i + before_outcome_col)] , na.rm = TRUE)
  )
}

sds[(ids+1) : (ids+2)] <- 1

diet_adult_Z <- diet_adult %>%
  mutate(SD = sds) %>%
  mutate(
    VN_OM_diff = VN_OM_diff / SD,
    VG_OM_diff = VG_OM_diff / SD,
    VN_VG_diff = VN_VG_diff / SD
  ) %>%
  select(-SD) %>%
  mutate(VG_OM_P = if_else(VG_OM_P < 1e-5, 1e-5, VG_OM_P)) %>%
  mutate(across(.cols = c(VN_OM_diff, VG_OM_diff, VN_VG_diff), ~ round(.x, 2)))


kbl(diet_adult_Z,
  caption =
    "Standardized differences in clinical outcomes (`outcome`) between diet groups in adults, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`)."
) %>%
  kable_styling("striped", full_width = T) %>%
  add_header_above(c(
    " " = 2, "(A) VN vs OM" = 3,
    "(B) VG vs OM" = 3, "(C) VN vs VG" = 3
  )) %>%
  column_spec(1, width_min = "1.5in")
Standardized differences in clinical outcomes (`outcome`) between diet groups in adults, estimated using linear mixed-effects models. For each outcome, the table reports the estimated difference (`[...]_diff`), raw p-value (`[...]_P`), and FDR-corrected p-value (`[...]_P_adj`) for multiple comparisons. Comparisons include (A) vegans vs. omnivores (`VN_OM`), (B) vegetarians vs. omnivores (`VG_OM`), and (C) vegans vs. vegetarians (`VN_VG`).
(A) VN vs OM
(B) VG vs OM
(C) VN vs VG
outcome estimand VN_OM_diff VN_OM_P VN_OM_P_adj VG_OM_diff VG_OM_P VG_OM_P_adj VN_VG_diff VN_VG_P VN_VG_P_adj
BMI mean -0.35 0.05644 0.17544 -0.10 0.62265 0.81153 -0.25 0.19057 0.42652
WHR mean -0.22 0.12747 0.27232 0.01 0.94638 0.98844 -0.23 0.12097 0.42652
HEIGHT mean -0.06 0.71233 0.77859 0.09 0.62498 0.81153 -0.14 0.36602 0.57343
MASS mean -0.32 0.04604 0.15455 -0.05 0.79770 0.89062 -0.27 0.09800 0.42652
HG mean 0.06 0.54439 0.62405 0.12 0.33940 0.79760 -0.05 0.62794 0.77666
PFAT mean -0.19 0.20782 0.34506 0.00 0.98770 0.99349 -0.20 0.21572 0.44534
log2_FAT mean -0.32 0.08157 0.19170 -0.07 0.74859 0.87960 -0.25 0.18161 0.42652
log2_FFM mean -0.20 0.03110 0.11687 -0.06 0.58978 0.81153 -0.14 0.13328 0.42652
SBP mean -0.30 0.06837 0.18663 -0.21 0.26104 0.72113 -0.09 0.60348 0.77666
DBP mean -0.28 0.11151 0.24958 0.24 0.24169 0.72113 -0.53 0.00498 0.08969
GLY mean -0.17 0.36725 0.49317 0.05 0.81483 0.89062 -0.22 0.25703 0.46464
log2_TC mean -0.70 0.00010 0.00151 -0.45 0.02523 0.23715 -0.25 0.16557 0.42652
HDL mean -0.10 0.56114 0.62794 -0.20 0.30203 0.74713 0.10 0.55223 0.77666
LDL mean -0.68 0.00027 0.00250 -0.34 0.10186 0.51319 -0.34 0.07216 0.42652
log2_TG mean -0.22 0.24651 0.37498 -0.13 0.53313 0.81153 -0.08 0.66848 0.78546
Ca mean 0.06 0.76142 0.81334 0.11 0.59798 0.81153 -0.06 0.76706 0.85837
P mean -0.18 0.26930 0.38813 0.00 0.99349 0.99349 -0.18 0.28863 0.50243
Mg mean -0.12 0.51192 0.60897 -0.18 0.38038 0.81153 0.06 0.73404 0.84146
log2_Se mean -0.29 0.20119 0.34506 -0.33 0.20654 0.72113 0.04 0.85731 0.87595
log2_Zn mean -0.74 0.00052 0.00404 -0.38 0.12011 0.51319 -0.37 0.08854 0.42652
FE mean 0.02 0.89833 0.93826 -0.16 0.44987 0.81153 0.19 0.33540 0.54357
log2_TIBC mean 0.25 0.20427 0.34506 0.15 0.50558 0.81153 0.10 0.62114 0.77666
log2_FERR mean -0.69 0.00003 0.00107 -0.47 0.01227 0.14415 -0.23 0.17286 0.42652
log2_TRF mean 0.24 0.21291 0.34506 0.14 0.52717 0.81153 0.10 0.61144 0.77666
SATTRF mean 0.00 0.99979 0.99979 -0.16 0.45260 0.81153 0.16 0.40363 0.61195
log2_TRFINDEX mean 0.51 0.00385 0.02011 0.53 0.00898 0.14061 -0.02 0.91768 0.91768
log2_STRF mean 0.21 0.24733 0.37498 0.26 0.22009 0.72113 -0.05 0.80476 0.85963
HGB mean -0.16 0.31410 0.43420 -0.09 0.63448 0.81153 -0.07 0.65602 0.78546
MCV mean 0.25 0.16119 0.31567 -0.22 0.27485 0.72113 0.47 0.01093 0.12840
PTH mean 0.41 0.03233 0.11687 -0.14 0.51452 0.81153 0.55 0.00553 0.08969
log2_CTX mean 0.00 0.99695 0.99979 -0.19 0.38453 0.81153 0.19 0.31117 0.52232
log2_P1NP mean 0.26 0.16018 0.31567 0.04 0.84126 0.89862 0.22 0.22741 0.44534
log2_UI mean -0.31 0.19890 0.34506 -0.45 0.11270 0.51319 0.14 0.56837 0.77666
log2_UREA mean -0.66 0.00022 0.00250 -0.16 0.42385 0.81153 -0.50 0.00573 0.08969
log2_CREA mean -0.45 0.00349 0.02011 -0.22 0.20547 0.72113 -0.23 0.14378 0.42652
UA mean -0.11 0.45945 0.56826 -0.08 0.63171 0.81153 -0.03 0.85613 0.87595
log2_HOLOTC mean -0.21 0.41789 0.54558 -0.53 0.04226 0.33106 0.31 0.15627 0.42652
log2_HCY mean -0.18 0.43623 0.55413 0.11 0.63887 0.81153 -0.29 0.15720 0.42652
log2_MMA mean 0.45 0.05973 0.17544 0.69 0.00308 0.14061 -0.24 0.22122 0.44534
log2_VIT_D mean 0.55 0.00857 0.04027 0.26 0.27618 0.72113 0.29 0.15364 0.42652
FOLATE mean 0.68 0.00070 0.00468 0.62 0.00622 0.14061 0.05 0.79120 0.85963
log2_Ur_Ca mean -0.37 0.07264 0.18663 -0.09 0.70675 0.85960 -0.28 0.17532 0.42652
log2_uCCR mean -0.22 0.27251 0.38813 -0.08 0.71328 0.85960 -0.13 0.49930 0.73335
log2_uICR mean -0.14 0.51827 0.60897 -0.41 0.11644 0.51319 0.26 0.24063 0.45239
uPCR mean -0.79 0.00005 0.00107 -0.39 0.06937 0.46575 -0.40 0.03866 0.36338
AL_adult log(OR) -1.95 0.01704 0.07280 -0.61 0.44626 0.81153 -1.34 0.08705 0.42652
BREAKS log(OR) -0.69 0.07545 0.18663 -0.12 0.78594 0.89062 -0.57 0.15320 0.42652

3.2.1.3 Vulcano plots

VN vs OM

Open code 

plotac <- "fig31_nr"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_adult_Z %>%
    filter(estimand != "log(OR)") %>%
    
    mutate(
      pv_adj = VN_OM_P_adj,
      pv_raw = VN_OM_P,
      b_coe = VN_OM_diff
    ) %>%
    mutate(
      pv_adj = VN_OM_P_adj,
      pv_raw = VN_OM_P,
      b_coe = VN_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VN+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.7,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs. OM, adults") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VG vs OM

Open code 
plotac <- "fig32_nr"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_adult_Z %>%
    filter(estimand != "log(OR)") %>%
    
    mutate(
      pv_adj = VG_OM_P_adj,
      pv_raw = VG_OM_P,
      b_coe = VG_OM_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VG+",
          "OM+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c("OM+", "VG+", "P>0.05")
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.6,
      size = 3.1,
      seed = 17
    ) +
    labs(x = "Standardized difference", y = "-log10 (P-value)") +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VG vs OM, adults") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

VN vs VG

Open code 
plotac <- "fig33_nr"

path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  assign(plotac, diet_adult_Z %>%
    filter(estimand != "log(OR)") %>%
    
    mutate(
      pv_adj = VN_VG_P_adj,
      pv_raw = VN_VG_P,
      b_coe = VN_VG_diff
    ) %>%
    mutate(pv_raw = if_else(pv_raw < 1e-5, 1e-5, pv_raw)) %>%
    mutate(
      m_ylab = -log10(pv_raw),
      m_xlab = b_coe,
      rescat = if_else(
        pv_raw >= 0.05,
        "P>0.05",
        if_else(
          b_coe > 0,
          "VN+",
          "VG+"
        )
      )
    ) %>%
    mutate(
      rescat = factor(rescat,
        levels = c(
          "VG+",
          "VN+",
          "P>0.05"
        )
      ),
      meta_labels = if_else(pv_raw < 0.05, outcome, "")
    ) %>%
    ggplot(aes(
      x = b_coe,
      y = m_ylab,
      color = rescat,
      alpha = rescat,
      label = meta_labels
    )) +
    geom_point(size = 2, shape = 16, alpha = 0.6) +
    geom_text_repel(
      aes(
        x = m_xlab,
        y = m_ylab
      ),
      show.legend = FALSE,
      alpha = 1,
      box.padding = 0.6,
      size = 3.1,
      seed = 16
    ) +
    labs(
      x = "Standardized difference",
      y = "-log10 (P-value)"
    ) +
    scale_color_manual(values = color_map, name = "") +
    ggtitle("VN vs VG, adults") +
    theme(
      legend.text = element_text(size = 11),
      axis.title = element_text(size = 12)
    ) +
    coord_cartesian(
      ylim = c(0, 5),
      xlim = c(-1.7, 1.7)
    ))

  saveRDS(get(plotac), file = path)
} else {
  assign(plotac, readRDS(path))
}

3.3 Vulcano matrix

Prepare

Open code 
plotac <- "fig35_nr"
path <- paste0("gitignore/figures/", plotac, ".rds")

if (file.exists(path) == FALSE) {
  fig19_nr <- fig19_nr + theme(
    legend.position = "none",
    plot.title = element_blank()
  )
  fig20_nr <- fig20_nr + theme(
    legend.position = "none",
    plot.title = element_blank()
  )
  fig21_nr <- fig21_nr + theme(
    legend.position = "none",
    plot.title = element_blank()
  )

  fig31_nr <- fig31_nr + theme(
    legend.position = "bottom",
    plot.title = element_blank()
  )
  fig32_nr <- fig32_nr + theme(
    legend.position = "bottom",
    plot.title = element_blank()
  )
  fig33_nr <- fig33_nr + theme(
    legend.position = "bottom",
    plot.title = element_blank()
  )
}

Create subfigures

Open code 

if (file.exists(path) == FALSE) {
  children_volc <- ggarrange(
    fig19_nr +
      theme(
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank()
      ),
    fig20_nr +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank()
      ),
    fig21_nr +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank()
      ),
    nrow = 1,
    ncol = 3,
    widths = c(1.1, 1, 1),
    labels = c("A", "", "")
  )

  children_volc <- annotate_figure(
    children_volc,
    top = text_grob("Children", face = "bold", size = 10)
  )


  adult_volc <- ggarrange(
    fig31_nr,
    fig32_nr +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank()
      ),
    fig33_nr +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank()
      ),
    nrow = 1,
    ncol = 3,
    widths = c(1.1, 1, 1),
    labels = c("B", "", "")
  )

  adult_volc <- annotate_figure(
    adult_volc,
    top = text_grob("Adults", face = "bold", size = 10)
  )
}

4 Reproducibility

Open code 
sessionInfo()
## R version 4.4.3 (2025-02-28)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=cs_CZ.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=cs_CZ.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=cs_CZ.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=cs_CZ.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Prague
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] mice_3.17.0        patchwork_1.2.0    ggrepel_0.9.5      robustlmm_3.3-1   
##  [5] gridExtra_2.3      pheatmap_1.0.12    performance_0.12.2 quantreg_5.98     
##  [9] SparseM_1.81       bayesplot_1.8.1    ggdist_3.3.2       kableExtra_1.4.0  
## [13] lubridate_1.8.0    corrplot_0.92      arm_1.12-2         MASS_7.3-65       
## [17] projpred_2.0.2     glmnet_4.1-8       boot_1.3-31        cowplot_1.1.1     
## [21] pROC_1.18.0        mgcv_1.9-1         nlme_3.1-168       openxlsx_4.2.8    
## [25] flextable_0.9.6    sjPlot_2.8.16      car_3.1-2          carData_3.0-5     
## [29] gtsummary_2.0.2    emmeans_1.10.4     ggpubr_0.4.0       lme4_1.1-35.5     
## [33] Matrix_1.7-0       forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
## [37] purrr_1.0.2        readr_2.1.2        tidyr_1.3.1        tibble_3.2.1      
## [41] ggplot2_3.5.1      tidyverse_1.3.1   
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.4.3           later_1.3.0             gamm4_0.2-6            
##   [4] cellranger_1.1.0        datawizard_0.12.2       rpart_4.1.24           
##   [7] reprex_2.0.1            lifecycle_1.0.4         rstatix_0.7.0          
##  [10] lattice_0.22-5          insight_0.20.2          backports_1.5.0        
##  [13] magrittr_2.0.3          rmarkdown_2.27          yaml_2.3.5             
##  [16] httpuv_1.6.5            zip_2.2.0               askpass_1.1            
##  [19] DBI_1.1.2               minqa_1.2.4             RColorBrewer_1.1-2     
##  [22] multcomp_1.4-18         abind_1.4-5             rvest_1.0.2            
##  [25] nnet_7.3-20             TH.data_1.1-0           sandwich_3.0-1         
##  [28] gdtools_0.3.7           crul_1.5.0              MatrixModels_0.5-3     
##  [31] svglite_2.1.3           codetools_0.2-19        xml2_1.3.3             
##  [34] tidyselect_1.2.1        shape_1.4.6             farver_2.1.0           
##  [37] ggeffects_1.7.0         httpcode_0.3.0          matrixStats_1.3.0      
##  [40] jsonlite_1.8.8          mitml_0.4-3             ellipsis_0.3.2         
##  [43] ggridges_0.5.3          survival_3.7-0          iterators_1.0.14       
##  [46] systemfonts_1.0.4       foreach_1.5.2           tools_4.4.3            
##  [49] ragg_1.2.1              Rcpp_1.0.13             glue_1.7.0             
##  [52] pan_1.6                 xfun_0.46               distributional_0.4.0   
##  [55] loo_2.4.1               withr_3.0.1             fastmap_1.2.0          
##  [58] fansi_1.0.6             openssl_1.4.6           digest_0.6.37          
##  [61] R6_2.5.1                mime_0.12               estimability_1.5.1     
##  [64] textshaping_0.3.6       colorspace_2.0-2        utf8_1.2.4             
##  [67] generics_0.1.3          fontLiberation_0.1.0    data.table_1.15.4      
##  [70] robustbase_0.93-9       httr_1.4.2              htmlwidgets_1.6.4      
##  [73] pkgconfig_2.0.3         gtable_0.3.0            htmltools_0.5.8.1      
##  [76] fontBitstreamVera_0.1.1 scales_1.3.0            knitr_1.48             
##  [79] rstudioapi_0.16.0       tzdb_0.2.0              uuid_1.0-3             
##  [82] coda_0.19-4             curl_4.3.2              nloptr_2.0.0           
##  [85] zoo_1.8-9               sjlabelled_1.2.0        parallel_4.4.3         
##  [88] pillar_1.9.0            vctrs_0.6.5             promises_1.2.0.1       
##  [91] jomo_2.7-3              dbplyr_2.1.1            xtable_1.8-4           
##  [94] evaluate_1.0.0          fastGHQuad_1.0.1        mvtnorm_1.1-3          
##  [97] cli_3.6.3               compiler_4.4.3          rlang_1.1.4            
## [100] crayon_1.5.0            rstantools_2.1.1        ggsignif_0.6.3         
## [103] modelr_0.1.8            plyr_1.8.6              sjmisc_2.8.10          
## [106] fs_1.6.4                stringi_1.7.6           viridisLite_0.4.0      
## [109] assertthat_0.2.1        munsell_0.5.0           fontquiver_0.2.1       
## [112] sjstats_0.19.0          hms_1.1.1               gfonts_0.2.0           
## [115] shiny_1.9.1             highr_0.11              haven_2.4.3            
## [118] broom_1.0.6             DEoptimR_1.0-10         readxl_1.3.1           
## [121] officer_0.6.6