In windyzn/urinaryDBP:

# Only set if the Rmd file is not in the parent directory (ie. 'projectname/')
knitr::opts_knit$set(root.dir = '../')

knitr::opts_chunk$set(collapse = TRUE, echo = FALSE, message = FALSE, warning = FALSE)
library(tidyverse)
devtools::load_all()
# load_data(update = TRUE)
set_options()
# extrafont::loadfonts(device="win")

# source('.Rprofile')
# run_setup()
# load_data('data/project_data.rda')
# dim(ds)
ds <- project_data

dsBase <- ds %>% 
  dplyr::filter(fVN == "Baseline")

# Subjects with measurements at all visit numbers

dsComplete <- ds %>%
  dplyr::group_by(SID) %>%
  dplyr::filter(n() == 3) %>%
  dplyr::ungroup()

## Include captions below using `captioner` package

fig_nums <- captioner::captioner(prefix = 'Figure')
supfig_nums <- captioner::captioner(prefix = 'Supplementary Figure')
tab_nums <- captioner::captioner(prefix = 'Table')
suptab_nums <- captioner::captioner(prefix = 'Supplementary Table')

PAPER 2: VITAMIN D RESULTS

Subject Characteristics

# Diet information is not available at baseline. Only VN 3 data is available at the present time.

tab_nums("subchar_vitd_baseline", "Subject characteristics across vitamin D status at baseline")

subchar_vitd_base <- tableone::CreateTableOne(
  vars = c("Age",
           "Sex",
           "Ethnicity",
           "BMI",
           "Waist",
           "eGFR",
           "ACR",
           "UrineCreatinine",
           "UrineMicroalbumin",
           "UrinaryCalcium",
           "UDBP",
           "udbpCr",
           "VitaminD",
           "PTH",
           "MET",
           "Systolic",
           "Diastolic",
           "MeanArtPressure",
           "diet_milk",
           "diet_cal",
           "diet_supp_cal",
           "diet_supp_vitd",
           "dmStatus"),
  strata = c("vitdStatus"),
  data = dsBase,
  factorVars = c("Sex", "Ethnicity", "dmStatus")
) %>% 
  print(nonnormal = c("UDBP",
                      "ACR",
                      "UrineMicroalbumin"),
        quote = FALSE,
        noSpaces = TRUE) %>% 
  knitr::kable()

# Diet information is not available at baseline. Only VN 3 data is available at the present time.

tab_nums("subchar_vitd", "Subject characteristics across visit numbers")

# for sensitivity analysis
m <- ds %>% 
  dplyr::filter(!is.na(VitaminD)) %>% 
  dplyr::filter(!VN == "6")

subchar_vitd_base <- tableone::CreateTableOne(
  vars = c("Age",
           "Sex",
           "Ethnicity",
           "BMI",
           "Waist",
           "MET",
           "VitaminD",
           "PTH",
           "ACR",
           "eGFR",
           "UDBP",
           "udbpCr",
           "Systolic",
           "Diastolic",
           "MeanArtPressure",
           "dmStatus"),
  strata = c("VN"),
  data = m,
  factorVars = c("Sex", "Ethnicity", "dmStatus")
) %>% 
  print(nonnormal = c("UDBP",
                      "udbpCr",
                      "ACR"),
        quote = FALSE,
        noSpaces = TRUE) %>% 
  knitr::kable()

ds %>%
  dplyr::select(VN, PTH) %>%
  dplyr::group_by(VN) %>%
  na.omit() %>%
  summarise(n = n())

Cross-Sectional

# Clean data
vitd <- dsBase %>% 
  dplyr::mutate(vitdLabels = ifelse(vitdStatus == "Sufficient", "Sufficient \n(≥75 nmol/L)",
                                   ifelse(vitdStatus == "Insufficient", "Insufficient \n(50 - 75 nmol/L)",
                                          "Deficient \n(<50 nmol/L)")),
                vitdLabels = factor(vitdLabels,
                                   levels = c("Sufficient \n(≥75 nmol/L)",
                                              "Insufficient \n(50 - 75 nmol/L)",
                                              "Deficient \n(<50 nmol/L)"),
                                   ordered = TRUE)) %>%
  select(vitdStatus, vitdLabels, udbpCr, UDBP) %>%
  na.omit()

# Box plot of uVDBP in different albuminuria categories
vitd %>% 
  box_plot_slides("vitdLabels", "log(udbpCr)", 
            "Vitamin D Status",
            "log uVDBP:cr (μg/mmol)")

# value
vitd %>%
  dplyr::group_by(vitdStatus) %>% 
  dplyr::summarise(n = n(),
                   Median = median(log(udbpCr)), 
                   lower = quantile(log(udbpCr), probs = 0.25),
                   higher = quantile(log(udbpCr), probs = 0.75),
                   IQR = IQR(log(udbpCr))
                   )

# n
vitd %>% 
  group_by(vitdStatus) %>% 
  summarise(n = n())

# ANOVA
anova <- aov(formula = log(udbpCr)~vitdStatus, data = vitd)
summary(anova)
TukeyHSD(anova)
rm(anova)

# Scatterplot of ACR and uVDBP ----------------------------------

dsBase %>% 
  # dplyr::filter(acrStatus != "Normoalbuminuria") %>% 
  scatter_plot_transparent("log(udbpCr)", "log(VitaminD)", 
                "log uVDBP:cr (μg/mmol)",
               "log serum 25(OH)D (nmol/L)")

# Spearman Correlation ------------------------------------------

dsBase %>% 
  cor.test(formula = ~ log(udbpCr) + log(VitaminD), data = ., method = "spearman")

# Linear Regression ---------------------------------------------

dsBase %>% 
  dplyr::mutate(lVitD = log(VitaminD)) %>% 
  prep_mason_data() %>% 
  mason_glm(y = "logudbpCr",
            x = "lVitD",
            covars = c("ageBase", "Sex", "Ethnicity", "MET", "Season", "dmStatus")
            ) %>% 
  gee_results_table() %>% 
  pander::pander()

Possible Effect Modifiers

# Scatterplot

dsBase %>% 
  scatter_plot("PTH", "VitaminD", 
               "Parathyroid Hormone (pmol/L)",
               "Serum 25(OH)D (nmol/L)")

# Linear Regression ---------------------------------------------

dsBase %>% 
  prep_mason_data() %>% 
  mason_glm(y = "VitaminD",
            x = "PTH"
            ) %>% 
  dplyr::filter(!term == "(Intercept)") %>%
  dplyr::mutate(p = round(p.value, 2),
    p = ifelse(p == "0", "<0.001", p),
    estCI = paste0(round(estimate, 2), " (",
                               round(conf.low, 2), ", ",
                               round(conf.high, 2), ")")) %>% 
  dplyr::select(Yterms, Xterms, term, estCI, p) %>% 
  # tidyr::spread(Yterms, estCI) %>%
  pander::pander()

# Box plot of 25(OH)D in different seasons
dsBase %>% 
  box_plot("Season", "VitaminD", 
            "Seasons",
            "Serum 25(OH)D (nmol/L)")

# n
dsBase %>% 
  group_by(Season) %>% 
  summarise(n = n())

# ANOVA
anova <- aov(formula = VitaminD~Season, data = dsBase)
summary(anova)
TukeyHSD(anova)
rm(anova)

Progression

ds %>% 
  plot_progress(yvar = "VitaminD",
                ylab = "Serum 25(OH)D")

# Complete data

dsComplete %>% 
  plot_progress(yvar = "VitaminD",
                ylab = "Serum 25(OH)D")

# ANOVA ------------------------------------------------------------

anova <- aov(formula = VitaminD~fVN, data = dsComplete)
summary(anova)
TukeyHSD(anova)
rm(anova)

Generalized Estimating Equations

# Predictor is baseline UDBP ----------------------------------------

gee_vitd_baseline <- ds %>% 
  prep_mason_data_vitd() %>% 
  mason_geeplot(y = "lVitD",
            x = "udbpCrBase",
            covars = c("YearsFromBaseline", "ageBase", "Sex", 
                       "Ethnicity", "BMI", "MET", "dmStatus", "Season"),
            intvar = "BMI") %>%
  mason::polish_renaming(rename_gee_vitd)

# GEE table ----------------------------------------------------------

gee_vitd_baseline %>% 
  gee_results_table() %>% 
  # tidyr::spread(Yterms, estCI) %>%
  pander::pander()

# Plot ---------------------------------------------------------------

plot_gee_results_vitd_base(gee_vitd_baseline,
                 yvars = "Serum 25(OH)D (nmol/L)")

# Plot interaction ---------------------------------------------------

# dsBase %>% 
#   interac_plot(yvar = "PTH",
#                intvar = "Season",
#                ylab = "Parathyroid Hormone (pmol/L)") #Serum 25(OH)D (nmol/L) #Parathyroid Hormone (pmol/L)
# 
# rms::ols(VitaminD ~ udbpCr * Season, data = dsBase)

# Predictor is UDBP ----------------------------------------

gee_vitd <- ds %>% 
  prep_mason_data_vitd() %>% 
  mason_geeplot(y = "lVitD",
            x = "udbpCr",
            covars = c("YearsFromBaseline", "ageBase", "Sex", "Ethnicity", "BMI",
                       "MET", "dmStatus", "Season"), 
            intvar = "BMI") %>%
  mason::polish_renaming(rename_gee_vitd)

# GEE table ----------------------------------------------------------

gee_vitd %>% 
  gee_results_table() %>% 
  # tidyr::spread(Yterms, estCI) %>%
  pander::pander()

# Plot ---------------------------------------------------------------

plot_gee_results_vitd(gee_vitd,
                 yvars = "Serum 25(OH)D (nmol/L)")

# Plot interaction ---------------------------------------------------

# ds %>% 
#   interac_plot(yvar = "VitaminD",
#                intvar = "Season",
#                ylab = "Serum 25(OH)D (nmol/L)")
# 
# rms::ols(VitaminD ~ udbpCr * Season, data = ds)