TODO:

knitr::opts_chunk$set(echo = TRUE,  
                      cache = TRUE,
                      warning = FALSE,
                      message = FALSE)

set.seed(2430024)

date <- lubridate::now()
format(date, "%a, %d. %B %Y")
## [1] "Wed, 12. February 2020"
# Load Packages
library(broman)
library(tidyverse)
library(knitr)
library(googlesheets4)
library(scales)
library(ggpubr)
library(magick)
library(gridExtra)
library(bayesplot)

# google sheets setup
googlesheet_id <- "1-0mfTM47nF2etqPqvcQzp9xspgTLNHjp-cZ0E1qGAjs"
sheets_auth(email = "jenshuesers@gmail.com")

# Import Local Package "PEDIS"
unloadNamespace("pedis")
devtools::load_all("~/r-projects/proj_wunde/pedis/")

# Data Import
pedis <- pedis::read_all_pedis("~/r-projects/proj_wunde/pedis-diabetic-care/datasets")

Descriptive Tables

desc_all <- descriptive_table_all(pedis)
sheets_write(data = desc_all,
             ss = googlesheet_id,
             sheet = "descr-table-all")

desc_any <- descriptive_table_outcome(pedis, outcome = "any_amputation")
sheets_write(data = desc_any,
             ss = googlesheet_id,
             sheet = "descr-table-any")

desc_maj <- descriptive_table_outcome(pedis, outcome = "major_amputation")
sheets_write(data = desc_maj,
             ss = googlesheet_id,
             sheet = "descr-table-maj")

Model Setup

# Setting model formula
f_any <- as.formula("any_amputation ~ p + e_ordinal_5 + d + i + s + alter_bei_aufnahme + gender")
f_major <- as.formula("major_amputation ~ p + e_ordinal_5 + d + i + s + alter_bei_aufnahme + gender")

# Setting priors

# Uninformed priors
uninformed_prior <- cauchy(location = c(rep(0, 5), 0, 0), scale = c(rep(1, 5), .5, .5))

source("./scripts/compute-informed-prior.R")

# Informed priors (from Pickwell et al.)
tibble(pooled_beta, pooled_sd) %>% 
  pivot_longer(cols = c("pooled_beta", "pooled_sd")) %>% 
  mutate(exp = exp(value)) %>% 
  set_names("statistic", "non-exponentiated", "exponentiated")
## # A tibble: 2 x 3
##   statistic   `non-exponentiated` exponentiated
##   <chr>                     <dbl>         <dbl>
## 1 pooled_beta               0.535          1.71
## 2 pooled_sd                 0.273          1.31
informed_prior <- cauchy(location = c(rep(pooled_beta, 5), 0, 0), scale = c(rep(pooled_sd, 5), .5, .5))

# Set number of MCMC interations
n_iter <-  4e4
seed <- 3412

Model Fitting

Outcome: Any Amputation

Non-Informed Prior

# Outcome: Any Amputation
# Approach: Non-Informed
noninformed_any <- fit_model(n_iter = n_iter, exponentiate = FALSE, seed = seed)

# Posterior Histogram with prior
errorbar_1 <- plot(noninformed_any[[1]]) + ggtitle("Non Informed Any")

# AUC Posterior
auc_noninformed_any <- posterior_auc(model = noninformed_any,
                                     data = pedis,
                                     outcome = "any_amputation")

Informed Prior

# Outcome: Any Amputation
# Approach: Informed Prior

informed_any <- fit_model(formula = f_any, prior = informed_prior, n_iter = n_iter, seed = seed)

# Posterior histograms with priors
errorbar_2 <- plot(informed_any[[1]]) + ggtitle("Informed Any")

# AUC Posterior
auc_informed_any <- posterior_auc(model = informed_any,
                                  data = pedis,
                                  outcome = "any_amputation")

Outcome: Major Amputation

Non-Informed Prior

# Outcome: Any Amputation
# Approach: Informed Prior

noninformed_major <- fit_model(formula = f_major, n_iter = n_iter, seed = seed)

# Posterior Histogram with priors
errorbar_3 <- plot(noninformed_major[[1]]) + ggtitle("Non-Informed Major")

# AUC Posterior
auc_noninformed_major <- posterior_auc(model = noninformed_major,
                                       data = pedis,
                                       outcome = "major_amputation")

Informed Prior

# Outcome: Major Amputation
# Approach: Informed Prior
informed_major <- fit_model(formula = f_major,
                            prior = informed_prior,
                            seed = seed,
                            n_iter = n_iter)

# Posterior histograms with priors
errorbar_4 <- plot(informed_major[[1]]) + ggtitle("Informed Major")

# AUC Posterior
auc_informed_major <- posterior_auc(model = informed_major,
                                    data = pedis,
                                    outcome = "major_amputation")

AUCs

AUC Descriptive Tables

auc <- list("Any Amputation - Non-Informed" = auc_noninformed_any, 
            "Any Amputation - Informed" = auc_informed_any, 
            "Major Amputation - Non-Informed" = auc_noninformed_major,
            "Major Amputation - Informed" = auc_informed_any)

# Descriptive Summary
auc_summarised <- map(auc, auc_summary) %>% 
  bind_rows(.id = "Model") %>% 
  mutate(HDI = paste0("[", myround(lower, 3), "-", myround(upper, 3), "]")) %>% 
  mutate_if(is.numeric, myround, 3)

auc_summarised %>%
  select(Model, med, HDI) %>% 
  kable(digits = 3)
Model med HDI
Any Amputation - Non-Informed 0.793 [0.778-0.801]
Any Amputation - Informed 0.790 [0.774-0.802]
Major Amputation - Non-Informed 0.765 [0.725-0.779]
Major Amputation - Informed 0.790 [0.774-0.802]
sheets_write(data = auc_summarised,
             ss = googlesheet_id,
             sheet = "auc")

AUC Posterior Distributions

# Plot Posterior AUC Distributions

posterior_auc_distributions <- auc %>% map(hist, bw = .005) 

get_y_axis <- function(plt_obj) ggplot_build(plt_obj)$layout$panel_scales_y[[1]]$range$range

y_max <- map_df(posterior_auc_distributions, get_y_axis) %>% 
  pivot_longer(cols = names(.)) %>% 
  pull(value) %>% 
  max %>% 
  `+`(2)

x_min <- min(map_dbl(auc, quantile, .02)) 
x_max <- max(map_dbl(auc, max)) + .01

x_axis <- scale_x_continuous(name = "AUC value",
                             limits = c(x_min, x_max),
                             breaks = seq(0, 1, by = .02))

y_axis <- scale_y_continuous(name = "Density", limits = c(0, y_max))

g <- ggpubr::ggarrange(
  posterior_auc_distributions[[1]] + 
    x_axis + 
    y_axis +
    rremove("x.title"),
  posterior_auc_distributions[[2]] + 
    x_axis + 
    y_axis +
    rremove("x.title") +
    rremove("y.title"),
  posterior_auc_distributions[[3]] +
    x_axis + 
    y_axis,
  posterior_auc_distributions[[4]] +
    x_axis + 
    y_axis +
    rremove("y.title"),
  ncol = 2, 
  nrow = 2,
  labels = c("A", "B", "C", "D"))

g

ggsave(plot = g,
       filename = "img/posterior-auc-values.png",
       unit = "cm", width = 30, height = 30,
       dpi = 320,
       scale = .5)

Posterior Parameter Distribution

Errorbars

posterior_errorbars <- gridExtra::grid.arrange(
  errorbar_1 + scale_y_continuous(breaks = 0:10) + ggtitle("Any Amputation - Non-Informed."),
  errorbar_2 + scale_y_continuous(breaks = 0:10) + ggtitle("Any Amputation - Informed"),
  errorbar_3 + scale_y_continuous(breaks = 0:10) + ggtitle("Major Amputation - Non-Informed"),
  errorbar_4 + scale_y_continuous(breaks = 0:10) + ggtitle("Major Amputation - Informed"))

ggsave(posterior_errorbars,
       filename = "img/posterior-parameter-errorbars.png",
       unit = "cm",
       width = 20,
       height = 20)

Descriptive Statistics

models <- list("Any Amputation - Non-Informed" = noninformed_any, 
               "Any Amputation - Informed" = informed_any, 
               "Major Amputation - Non-Informed" = noninformed_major,
               "Major Amputation - Informed" = informed_major) 

summary_models <- map(models, summary_coefs) %>% bind_rows(.id = "model")

sheets_write(data = summary_models,
             ss = googlesheet_id,
             sheet = "coefficients")

Parameter Posterior Distributions

xlimit <- c(-1.2, 3)
ylimit <- c(0, 2.4)

hist_posterior <- models %>% 
  map(pluck, 3) %>% 
  map(function(x) x +
        # scale_x_continuous("Value of Coefficients", limits = c(-1.2, 3)) +
        theme(strip.text = element_text(size = 8),
              axis.text = element_text(size = 8)) +
        rremove("x.title") +
        rremove("y.title") +
        coord_cartesian(ylim = ylimit)
  )
  
hist_plot <- ggarrange(hist_posterior[[1]],
                       hist_posterior[[2]],
                       hist_posterior[[3]],
                       hist_posterior[[4]],
                       labels = LETTERS[1:4],
                       ncol = 4)

hist_plot