rrnmr.Rmd

---
title: "NMR metabolites and hypertension"
author: "Joonatan Palmu"
date: "`r format(Sys.time(), '%d.%m.%Y')`"
output: 
  html_document:
    number_sections: true
---

```{r options, include=FALSE}
knitr::opts_chunk$set(include = TRUE, echo = TRUE, message = FALSE, results='asis', cache = FALSE, warning=FALSE, error = FALSE)
options(knitr.kable.NA = "")
```

# Libraries

<details><summary>Libaries</summary>

```{r libraries, results = 'hide'}
# BASIC
library(magrittr)
library(knitr)
library(kableExtra)
library(tidyr)
library(broom)
# MODELING
library(dplyr)
library(parallel)
library(officer)
library(flextable)
library(lme4)
# PLOT
library(ggplot2)
library(ggrepel)
library(ggeffects)
library(isoband)
library(randomForest)
library(png)
library(RColorBrewer)
library(ggquiver)
library(ggcorrplot)
# OTHER
library(Boruta)
library(pander)
library(caret)
library(pROC)
library(broom.mixed)
library(haven)
library(readxl)
library(openxlsx)
library(ggforestplot)
library(patchwork)
library(ggfortify)
library(cluster)
library(tableone)
library(purrr)
library(nnet)
library(ordinal)
library(egg)
library(gridExtra)
library(sure)
library(stringr)
library(haven)
library(car)
library(xfun)
library(PredictABEL)
library(pROC)
library(stargazer)
library(nricens)
library(metafor)
library(psych)
```

</details>

## Functions

<details><summary>Functions</summary>

```{r Functions}
sourcefiles <- c("articles-functions.R",
                 "articles-models.R",
                 "articles-plots.R",
                 "articles-tables.R")
```

```{r import files, echo = FALSE}
for (f in sourcefiles) {
    source(f)
}
```

```{r embbed files, echo = FALSE}
xfun::embed_files(c("rrnmr.Rmd", sourcefiles))
```

</details>

#  Folder structure

```{r clean folders}
mkdir("cache", "report", "session")
```


# Definitions

```{r definitions}
dvars <- c("sys", "dias", "htn")
ivars <- c("age", "female", "bmi", "smoker", "diabetes", "leisure", "bptreat", "koltreat", "cohort")
```

# Data import

```{r import raw data}
dset <- readRDS("rds/dset-20210217-102531.rds")
metabolites <- getmetabolites(dset$crosssectional)
metabolites.subclass <- getsubclasses(metabolites)
```

## Characteristics

<details><summary>Cross-sectional sample</summary>

Median age is `r median(dset$crosssectional$age)`.

```{r tableone cross-sectional, echo = FALSE}
characteristicsTableFull(dset$crosssectional, test = TRUE) %>% mykable
```

</details>

<details><summary>Longitudinal sample</summary>

At baseline

```{r tableone longitudinal baseline, echo = FALSE}
dset$longitudinal %>%
    filter(cohort %in% c("F2007", "T2000")) %>%
    characteristicsTable %>%
    mykable
```

Both

```{r tableone longitudinal, echo = FALSE}
characteristicsTable(dset$longitudinal, "cohort") %>% mykable
```

</details>


## Correlation matrix

<details><summary>Correlation matrix</summary>

```{r heatmap for supplement, echo = FALSE}
lipoproteins_for_correlation_matrix <- tibble(id = metabolites,
                                              desc = bioproperty(id, length = 100),
                                              group = bioproperty(id, "group")) %>%
    filter(group != "Lipoprotein subclasses") %>%
    mutate(group_name = id)

spearman_correlation_matrix <- correlationmatrix(dset$crosssectional,
                                                 lipoproteins_for_correlation_matrix,
                                                 method = "spearman",
                                                 adjust = "fdr")

pheatmap::pheatmap(spearman_correlation_matrix[["r"]],
                   display_numbers = ifelse(spearman_correlation_matrix[["p"]] < 0.05, "", "×"),
                   fontsize = 6,
                   fontsize_row = 6,
                   fontsize_col = 5,
                   fontsize_number = 10,
                   cutree_rows = 6,
                   cutree_cols = 6,
                   labels_row = lipoproteins_for_correlation_matrix$desc,
                   labels_col = lipoproteins_for_correlation_matrix$group,
                   cellwidth = 6,
                   cellheight = 6,
                   angle_col = 45,
                   treeheight_col = 0,
                   silent = TRUE,
                   legend_breaks = c(-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1)) %>%
ggsave(file = "cache/crosssectional-correlationmatrix.png",
           plot = .,
           height = 5.6,
           width = 8,
           dpi = 300,
           unit = "in")
```

<img src = "cache/crosssectional-correlationmatrix.png" />

</details>

<details><summary>Correlation matrix lipoprotein subclasses</summary>

```{r lipoprotein subclasses heatmap for supplement, echo = FALSE}
lipoproteins_for_correlation_matrix_subclass <- tibble(id = metabolites,
                                                       desc = bioproperty(id, length = 100),
                                                       name = bioproperty(id, property = "abbreviation", length = 100),
                                                       group = bioproperty(id, "group")) %>%
    filter(group == "Lipoprotein subclasses") %>%
    mutate(group_name = id)

spearman_correlation_matrix_subclass <- correlationmatrix(dset$crosssectional,
                                                          lipoproteins_for_correlation_matrix_subclass,
                                                          method = "spearman",
                                                          adjust = "fdr")

pheatmap::pheatmap(spearman_correlation_matrix_subclass[["r"]],
                                     display_numbers = ifelse(spearman_correlation_matrix_subclass[["p"]] < 0.05, "", "×"),
                                     fontsize = 6,
                                     fontsize_row = 6,
                                     fontsize_col = 5,
                                     fontsize_number = 10,
                                     cutree_rows = 4,
                                     cutree_cols = 4,
                                     labels_row = lipoproteins_for_correlation_matrix_subclass$name,
                                     labels_col = lipoproteins_for_correlation_matrix_subclass$desc,
                                     cellwidth = 6,
                                     cellheight = 6,
                                     angle_col = 45,
                                     treeheight_col = 0,
                                     silent = TRUE,
                                     legend_breaks = c(-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1)) %>%
                   ggsave(file = "cache/crosssectional-correlationmatrix-subclass.png",
                          plot = .,
                          height = 8.7,
                          width = 8.9,
                          dpi = 300,
                          unit = "in")
```

<img src = "cache/crosssectional-correlationmatrix-subclass.png" />

</details>

## Cross-sectional associations

### Associations between metabolites and htn

```{r plot cross-sectional models}
glm.sys <- loop.lm(dset$crosssectional, "sys", metabolites, ivars)
glm.dias <- loop.lm(dset$crosssectional, "dias", metabolites, ivars)
glm.htn <- loop.binomial(dset$crosssectional, "htn", metabolites, ivars)
results.crosssectional <- loop.results(glm.sys, glm.dias, glm.htn)
```

q<details><summary>Logistic model</summary>

```{r ggforestplot for logistic models, echo = FALSE}
myforestplot(results.crosssectional,
             xlim = c(0.7, 1.5),
             dims = c(18, 23),
             logodds = TRUE,
             responses = c("htn"),
             file = "cache/glmforest-crosssectional-htn.png")
```

<img src="cache/glmforest-crosssectional-htn.png" />

</details>

<details><summary>Linear models</summary>

```{r ggforestplot for linear models, echo = FALSE}
myforestplot(results.crosssectional,
             xlim = c(-2, 3),
             dims = c(18, 23),
             responses = c("sys", "dias"),
             file = "cache/glmforest-crosssectional-bp.png")
```

<img src="cache/glmforest-crosssectional-bp.png" />

</details>

<details><summary>Lipoprotein subclasses histogram</summary>

```{r density plot subcalsses crosssectional, echo = FALSE}
{ dset$crosssectional %>%
      gather(Metabolite, Value, colnames(dset$crosssectional) %>% mygrep(word = "NMR_.*_.*DL_.*")) %>%
      ggplot(aes(x = Value, col = htn)) +
      facet_wrap(~Metabolite, ncol=5) +
      geom_histogram(binwidth = 1) +
      scale_y_continuous(trans='log10') } %>%
    ggsave(file = "cache/crosssectional-subclasses-densityplot.png",
           plot = .,
           height = 16,
           width = 12,
           dpi = 300,
           unit = "in")
```

<img src = "cache/crosssectional-subclasses-densityplot.png" />

</details>

<details><summary>Lipoprotein subclasses</summary>

```{r ggforestplot for linear model subclasses, echo = FALSE}
myforestplot(results.crosssectional,
             xlim = c(-2, 6),
             dims = c(18, 30),
             legend.position = "bottom",
             responses = c("sys", "dias"),
             panelone = c("Very low-density lipoprotein",
                          "Intermediate-density lipoprotein"),
             paneltwo = c("Low-density lipoprotein",
                          "High-density lipoprotein",
                          "General lipoprotein subclasses"),
             preprocess = subclasspreprocess,
             file = "cache/glmforest-crosssectional-liposubclasses-bp.png")
```

<img src = "cache/glmforest-crosssectional-liposubclasses-bp.png" />

```{r ggforestplot for logistic model subclasses, echo = FALSE}
myforestplot(results.crosssectional,
             xlim = c(0.7, 1.5),
             dims = c(18, 28),
             logodds = TRUE,
             responses = c("htn"),
             panelone = c("Very low-density lipoprotein",
                          "Intermediate-density lipoprotein"),
             paneltwo = c("Low-density lipoprotein",
                          "High-density lipoprotein",
                          "General lipoprotein subclasses"),
             preprocess = subclasspreprocess,
             file = "cache/glmforest-crosssectional-liposubclasses-htn.png")
```

<img src = "cache/glmforest-crosssectional-liposubclasses-htn.png" />

</details>

<details><summary>Residuals</summary>

```{r residuals, echo = FALSE}
ggsave(file="cache/glmforest-crosssectional-residuals-sys.png",
       arrangeGrob(grobs = loop.residuals(glm.sys, ignore = metabolites.subclass), ncol = 6),
       height = 16,
       width = 12,
       dpi = 300,
       unit = "in")
```

<img src="cache/glmforest-crosssectional-residuals-sys.png" />

</details>

<details><summary>Numerical values</summary>

```{r ggforestplot numeric values, echo = FALSE}
results.crosssectional %>% 
    results.table %>%
    mykable
```

```{r size of association, include = FALSE, eval = FALSE}
results.crosssectional %>%
    filter(response == "sys") %>%
    filter(qval < 0.05) %>%
    mutate(order = sign(estimate)) %>%
    group_by(-order) %>%
    summarise(estimate_median = median(estimate),
              estimate_mean = mean(estimate),
              estimate_sd = sd(estimate),
              estimate_min = min(estimate),
              estimate_max = max(estimate))
```

</details>

### Compare young and old

```{r crosssectional young and old}
glm.age.young <- dset$crosssectional %>%
    filter(agegroup == "young") %>%
    rename(htn.young = htn) %>%
    loop.binomial(., "htn.young", metabolites, ivars)

glm.age.old <- dset$crosssectional %>%
    filter(agegroup == "old") %>%
    rename(htn.old = htn) %>%
    loop.binomial(., "htn.old", metabolites, ivars)

glm.age <- loop.results(glm.age.young, glm.age.old)
```

<details><summary>Logistic models</summary>

```{r ggforestplot crosssectional age model, echo = FALSE}
myforestplot(glm.age,
             xlim = c(0.7, 1.5),
             dims = c(18, 23),
             logodds = TRUE,
             scale = agescale,
             responses = c("htn.old", "htn.young"),
             file = "cache/glmforest-crosssectional-age.png")
```

<img src = "cache/glmforest-crosssectional-age.png" />

</details>

<details><summary>Numerical values</summary>

```{r ggforestplot age diff numeric values, echo = FALSE}
glm.age %>% 
    results.table %>%
    mykable
```

</details>

### Compare sex

```{r crosssectional sex}
glm.sex.female <- dset$crosssectional %>%
    filter(female == 1) %>%
    rename(htn.female = htn) %>%
    loop.binomial(., "htn.female", metabolites, ivars %difference% "female")

glm.sex.male <- dset$crosssectional %>%
    filter(female == 0) %>%
    mutate(htn.male = htn) %>%
    loop.binomial(., "htn.male", metabolites, ivars %difference% "female")

glm.sex <- loop.results(glm.sex.male, glm.sex.female)
```

<details><summary>Logistic models</summary>

```{r ggforestplot crosssectional sex model, echo = FALSE}
myforestplot(glm.sex,
             xlim = c(0.7, 1.5),
             dims = c(18, 23),
             logodds = TRUE,
             scale = sexscale,
             responses = c("htn.female", "htn.male"),
             file = "cache/glmforest-crosssectional-sex.png")
```

<img src = "cache/glmforest-crosssectional-sex.png" />

</details>

<details><summary>Numerical values</summary>

Numerical values

```{r ggforestplot sex diff numeric values, echo = FALSE}
glm.sex %>% 
    results.table %>%
    mykable
```

</details>

### Comparing cross-sectional samples

```{r Comparing cross-sectional samples}
glm.compare <- dset$crosssectional %>%
    group_by(cohort) %>% 
    dplyr::do(loop.binomial(., "htn", metabolites, ivars %difference% "cohort") %>%
              loop.results) %>%
    ungroup %>%
    mutate(response = cohort)
```

<details><summary>Comparing cross-sectional samples</summary>

```{r comparing crosssectional samples plot, echo = FALSE}
myforestplot(glm.compare,
             xlim = c(0.4, 2.5),
             dims = c(18, 23),
             logodds = TRUE,
             responses = c("1997", "2002", "2007", "2012", "T2000"),
             file = "cache/glmforest-crosssectional-cohort.png",
             scale = cohortscale)
```

<img src = "cache/glmforest-crosssectional-cohort.png" />

</details>

## Longitudinal set

```{r longitudinal model}
glm.longitudinal.sys <- loop.lm(dset = getincidental(dset$longitudinal),
                                response = "sys.followup",
                                loops = metabolites,
                                covariates = ivars %union% "sys")
```

<details><summary>Numerical values for systole</summary>

```{r longitudinal model numerical results sys, echo = FALSE}
( results.longitudinal <- loop.results(glm.longitudinal.sys) ) %>%
    results.table %>%
    mykable    
```

</details>

<details><summary>Plot systole</summary>

```{r ggforestplot for incidental model all }
myforestplot(results.longitudinal,
             responses = "sys.followup",
             logodds = FALSE,
             xlim = c(-2, 3),
             dims = c(18, 23),
             file = "cache/glmforest-longitudinal-sys.png")
```

<img src = "cache/glmforest-longitudinal-sys.png" />

</details>

<details><summary>Lipoprotein subclasses histogram</summary>

```{r density plot subcalsses longitudinal, echo = FALSE}
{ getincidental(dset$longitudinal) %>%
    gather(Metabolite, Value, colnames(dset$longitudinal) %>% mygrep(word = "NMR_.*_.*DL_.*")) %>%
      ggplot(aes(x = Value, col = htn)) +
      facet_wrap(~Metabolite, ncol=5) +
      geom_histogram(binwidth = 1) +
      scale_y_continuous(trans='log10') } %>%
    ggsave(file = "cache/longitudinal-subclasses-densityplot.png",
       plot = .,
       height = 16,
       width = 12,
       dpi = 300,
       unit = "in")
```

<img src = "cache/longitudinal-subclasses-densityplot.png" />

</details>

<details><summary>Lipoprotein subclasses</summary>

```{r ggforestplot for longitudinal model subclasses,  results = 'hide', echo = FALSE}
myforestplot(results.longitudinal,
             xlim = c(-2, 2),
             dims = c(18, 28),
             responses = "sys.followup",
             panelone = c("Very low-density lipoprotein",
                          "Intermediate-density lipoprotein"),
             paneltwo = c("Low-density lipoprotein",
                          "High-density lipoprotein",
                          "General lipoprotein subclasses"),
             preprocess = subclasspreprocess,
             file = "cache/glmforest-longitudinal-liposubclasses-bp.png")
```

<img src = "cache/glmforest-longitudinal-liposubclasses-bp.png" />

</details>

<details><summary>Residuals systole</summary>

```{r longitudinal residuals sys, echo = FALSE}
ggsave(file = "cache/residuals-longitudinal-sys.png",
       plot = grid.arrange(grobs = loop.residuals(glm.longitudinal.sys, ignore = metabolites.subclass),
                           ncol = 6),
       height = 16,
       width = 12,
       dpi = 300,
       unit = "in")
```

<img src = "cache/residuals-longitudinal-sys.png" />

```{r longitudinal qq sys, echo = FALSE}
ggsave(file = "cache/qq-longitudinal-sys.png",
       plot = grid.arrange(grobs = loop.qq(glm.longitudinal.sys, ignore = metabolites.subclass),
                           ncol = 6),
       height = 16,
       width = 12,
       dpi = 300,
       unit = "in")
```

<img src = "cache/qq-longitudinal-sys.png" />

</details>

## Meta-analyses for main results

### Cross-sectional


```{r meta-analysis for cross-sectional}
glm.meta.sys <- loop.rma(dset = dset$crosssectional,
                         response = "sys",
                         loops = metabolites,
                         covariates = ivars %difference% "cohort") %>%
    loop.results(idtoterm = TRUE) %>%
    mutate(response = "sys")

glm.meta.dias <- loop.rma(dset = dset$crosssectional,
                         response = "dias",
                         loops = metabolites,
                         covariates = ivars %difference% "cohort") %>%
    loop.results(idtoterm = TRUE) %>%
    mutate(response = "dias")

glm.meta <- rbind(glm.meta.sys, glm.meta.dias)
```

<details><summary>Linear models</summary>

```{r ggforestplot for linear meta-analysis models, echo = FALSE}
myforestplot(glm.meta,
             xlim = c(-2, 3),
             dims = c(18, 23),
             responses = c("sys", "dias"),
             file = "cache/glmforest-crosssectional-bp-meta.png")
```

<img src="cache/glmforest-crosssectional-bp-meta.png" />

</details>


### Longitudinal

```{r meta-analysis for longitudinal}
glm.longitudinal.meta.sys <- loop.rma(dset = getincidental(dset$longitudinal),
                                      response = "sys.followup",
                                      loops = metabolites,
                                      covariates = ivars %union% "sys" %difference% "cohort") %>%
    loop.results(idtoterm = TRUE) %>%
    mutate(response = "sys.followup")
```

<details><summary>Linear models</summary>

```{r ggforestplot for longitudinal models, echo = FALSE}
myforestplot(glm.longitudinal.meta.sys,
             xlim = c(-2, 3),
             dims = c(18, 23),
             responses = "sys.followup",
             file = "cache/glmforest-longitudinal-sys-meta.png")
```

<img src="cache/glmforest-longitudinal-sys-meta.png" />

</details>


# Supplement

<details><summary>Metabolite names</summary>

```{r metabolite names, echo = FALSE}
table.metabolites <- lapply(c2l(metabolites), function(x)
    data_frame(abbreviation = bioproperty(x),
               description = bioproperty(x, "description"),
               group = bioproperty(x, "group"),
               N = sum(!is.na(dset$crosssectional[[x]])),
               mean = attr(dset$crosssectional[[x]], "scaled:center"),
               sd = attr(dset$crosssectional[[x]], "scaled:scale"),
               unit = bioproperty(x, "unit"))) %>%
    purrr::map_df(., ~as.data.frame(.x)) %>%
    mutate_at(vars(mean, sd), round, 3) %>%
    arrange(abbreviation) 
```

```{r Metabolites names }
table.metabolites %>%
    filter(group != "Lipoprotein subclasses") %>%
    mykable
```

</details>