med_iqr <- function(x){

  med_x <- round(median(x, na.rm = T), 1)
  iqr_x <- round(quantile(x, probs = c(0.25, 0.75), na.rm = T), 1)

  res <- "-"

  if(!is.na(med_x)){
    res <- paste0(med_x, " (", iqr_x[1], "-", iqr_x[2],")")
  }

  return(res)

}

number_percent <- function(df, gr1, gr2){

  df %>% group_by(!!sym(gr1), !!sym(gr2)) %>%
    summarise(n = n()) %>%
    mutate(perc = round(n/sum(n) * 100, 2)) %>%
    mutate(comb =  paste0(n, " (", perc,"%)")) %>%
    tidyr::pivot_wider(id_cols = c(gr1, gr2), names_from = gr1, values_from = "comb") %>%
    mutate(across(.cols = c(-!!sym(gr2)),.fns = function(x) replace(x, which(is.na(x)), "0 (0.00%)")))

}

number_percent_one_factor <- function(df, gr1){

  df %>% group_by(!!sym(gr1)) %>%
    summarise(n = n(), .groups = "drop") %>%
    mutate(perc = round(n/sum(n) * 100, 2)) %>%
    mutate(comb =  paste0(n, " (", perc,"%)")) %>%
    tidyr::pivot_wider(id_cols = c(gr1), names_from = gr1, values_from = "comb") #%>%
    # mutate(across(.cols = c(-!!sym(gr1)),.fns = function(x) replace(x, which(is.na(x)), "0 (0.00%)")))

}

mut_num_perc <- function(x, char_sel){
  
  abs <- sum(x == char_sel)
  rel <- round(abs/length(x)* 100, digits = 2)
  
  return(paste0(abs, " (", rel, "%)"))
}

fisher_test_table <- function(df, gr_row, gr_col){
  
  ft_table <- df %>% group_by(!!sym(gr_row), !!sym(gr_col)) %>% count(name = "Number of Subjects") %>% 
            tidyr::pivot_wider(id_cols = c(gr_col, "Number of Subjects"), names_from = gr_row, values_from = "Number of Subjects") 

  ft_res <- fisher.test(ft_table[,-1])

  kable_table <- df %>%  group_by(!!sym(gr_row), !!sym(gr_col)) %>% count(name = "n") %>% 
            group_by(!!sym(gr_col)) %>% mutate(val = paste0(n, " (", round(n/sum(n) * 100, 2), "%)" )) %>%
            tidyr::pivot_wider(id_cols = c(gr_row, "val", gr_col), names_from = gr_col, values_from = "val")

  kable_table <- cbind(kable_table, "Odds-Ratio" = round(ft_res$estimate, 2), "P-value" = ft_res$p.value)
  
  kable_table <- kable_table %>% mutate(`P-value` = dplyr::case_when(
                                                  `P-value` >= 0.01 ~ as.character(round(`P-value`, 2)),
                                                  `P-value` < 0.0001 ~ "< 0.0001",
                                                  `P-value` < 0.001 ~ "< 0.001",
                                                  `P-value` < 0.01 ~ "< 0.01")) 
  knitr::kable(kable_table, align = "c") %>%
    kableExtra::kable_styling(full_width = F, position = "left")%>%
    kableExtra::column_spec(1, bold = T) %>%
    kableExtra::collapse_rows(columns = 3:5, valign = "middle")
}

chisq_test_table <- function(df, gr_row, gr_col){
  
  chi_table <- df %>% group_by(!!sym(gr_row), !!sym(gr_col)) %>% count(name = "Number of Subjects") %>% 
            tidyr::pivot_wider(id_cols = c(gr_row, "Number of Subjects"), names_from = gr_col, values_from = "Number of Subjects")

  chisq_res <- chisq.test(chi_table[,-1])

  kable_table <- df %>% group_by(!!sym(gr_row), !!sym(gr_col)) %>% count(name = "n") %>% 
            group_by(!!sym(gr_col)) %>% mutate(val = paste0(n, " (", round(n/sum(n) * 100, 2), "%)" )) %>%
            tidyr::pivot_wider(id_cols = c(gr_row, "val", gr_col), names_from = gr_col, values_from = "val")

  knitr::kable(cbind(kable_table, "Odds-Ratio" =  "-", "P-value" = round(chisq_res$p.value, 2)), align = "c") %>%
    kableExtra::kable_styling(full_width = F, position = "left")%>%
    kableExtra::column_spec(1, bold = T) %>%
    kableExtra::collapse_rows(columns = dim(kable_table)[2]:(dim(kable_table)[2] + 2), valign = "middle")
}

model_specifics <- function(model, data){
  
  prds <- model %>% predict(data) 
  rmse <- caret::RMSE(prds, data$logIgG)
  r2 <- caret::R2(prds, data$logIgG)
  
  list(RMSE = rmse, R2 = r2)
  
}


fisher_tests_gr <- function(df, gr1, var_expr){

  df_subset <- df %>% select(!!sym(gr1), dplyr::matches(var_expr)) %>% na.omit()
  fisher_tests <- c()

  for(nw in colnames(df_subset %>% select(dplyr::matches(var_expr)))){

    conf_mat <- df_subset %>% dplyr::group_by(!!sym(gr1), !!sym(nw)) %>% 
                               count(name = "n") %>% ungroup() %>% 
                               tidyr::pivot_wider(id_cols = c(gr1, nw), values_from = c("n"), names_from = c(nw))
  
    fisher_res <-  broom::tidy(fisher.test(conf_mat[,-1], alternative = "two.sided"))
  
    colnames(conf_mat)[-1] <- c("gr2.1", "gr2.2")
    fisher_tests <- rbind(fisher_tests, cbind(group1 = conf_mat %>% pull(!!sym(gr1)) %>% as.character(), group2 = nw, conf_mat[,-1], fisher_res))
  }

  return(fisher_tests)
}

fisher_plot_odds <- function(fisher_test_df, ylab_add = "Coadministration", wrap = F, ylim = NULL){
  
  fisher_plt_df <- fisher_test_df %>%
                   select(group2, estimate, matches("conf|super"), p.value) %>%
                   unique() %>% mutate(p.adj = p.adjust(p.value, "BY"))
  
  odds_plot <- ggplot2::ggplot(data = fisher_plt_df, ggplot2::aes(x = group2, y = estimate))
  
  if(wrap){
    odds_plot <- odds_plot + ggplot2::facet_wrap(~super_group, nrow = 2)
  }
  odds_plot <- odds_plot + 
               ggplot2::geom_point(cex = 3) +
               ggplot2::geom_errorbar(ggplot2::aes(ymin = conf.low, ymax = conf.high), width = 0.5,
                                      position = ggplot2::position_dodge2(), show.legend = FALSE) +
               ggplot2::theme_bw() +
               ggplot2::scale_x_discrete(labels = c("Local reactions", "Headache", "Muscle Pain", "Fever and/or Chills", "Fatigue", "Other s.e.")) +
               ggsci::scale_color_nejm() +
               ggplot2::theme(axis.text = ggplot2::element_text(size = 12, colour = "black"),
                              axis.text.x = ggplot2::element_text(angle = 45, hjust = 1),
                              axis.title = ggplot2::element_text(face = "bold", size = 14),
                              legend.title = ggplot2::element_blank(),
                              legend.text = ggplot2::element_text(size = 12),
                              legend.position = "bottom",
                              strip.text = ggplot2::element_text(size = 14, face = "bold"), 
                              strip.background = ggplot2::element_rect(fill = "white"))+
               ggplot2::ylab(paste0("OR ", ylab_add)) +
               ggplot2::xlab("side effects")
  
  if(!is.null(ylim)){
    
    odds_plot <- odds_plot + ggplot2::ylim(ylim)
  }
  
  if(dim( fisher_plt_df %>% filter(p.adj < 0.05))[1] > 0){
     odds_plot <- odds_plot +  ggplot2::geom_text(data = fisher_plt_df %>% filter(p.adj < 0.05), ggplot2::aes(x = group2, y = conf.high, label = paste0("italic(p)== ",round(p.adj, 2))), vjust = -0.5, parse = T)
  }
  
  odds_plot
}


barplot_sideeffects_percentage <- function(p.adj_fisher_tests_df, wrap = F){
  
  plt_df <- c()
  
  if(wrap){
    plt_df <- p.adj_fisher_tests_df %>% group_by(super_group, group1, group2) %>% 
              mutate(n =  gr2.1 + gr2.2,
                     percentage = round(gr2.1/n * 100, 2)) %>% 
                     select(group1, group2, percentage, p.adjusted)
  }else{
    
    plt_df <- p.adj_fisher_tests_df %>% group_by(group1, group2) %>% 
              mutate(n =  gr2.1 + gr2.2,
                     percentage = round(gr2.1/n * 100, 2)) %>% 
                     select(group1, group2, percentage, p.adjusted)
  }
  
  p_val_data <- plt_df %>% mutate(group1 = plt_df$group1[1]) %>% 
                select(group2, p.adjusted, percentage)  %>% 
                mutate(percentage = max(percentage)) %>% 
                unique() %>% ungroup()
  
  bar_plt <- ggplot2::ggplot(data = plt_df, ggplot2::aes(x = group2, y = percentage, fill = group1))

  if(wrap){
    
    bar_plt <- bar_plt + ggplot2::facet_wrap(~super_group, nrow = 2)
  }

  bar_plt <- bar_plt + 
             ggplot2::geom_bar(stat = "identity", position = ggplot2::position_dodge2()) +
             ggplot2::xlab("Side effects") +
             ggplot2::ylab("Percentage of participants") +
             ggplot2::theme_bw() +
             ggsci::scale_fill_nejm() +
             ggplot2::scale_x_discrete(labels = c("Local reactions",
                                                  "Headache", "Muscle Pain", "Fever and/or Chills", "Fatigue", "Other")) +
             ggplot2::scale_y_continuous(expand = c(0,1), limits = c(0,100)) +
             ggplot2::guides(fill = ggplot2::guide_legend(title = "Effect")) +
             ggplot2::theme(axis.title = ggplot2::element_text(color = "black", size = 13, face = "bold"),
                            axis.text = ggplot2::element_text(color = "black", size = 12),
                            axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, vjust = 1),
                            strip.text = ggplot2::element_text(size = 14, face = "bold"), 
                            strip.background = ggplot2::element_rect(fill = "white"),
                            legend.title = ggplot2::element_blank(),
                            legend.text = ggplot2::element_text(size = 10),
                            axis.title.x = ggplot2::element_blank(),
                            legend.position = "bottom",
                            legend.justification = c("left", "top"),
                            legend.box.just = "left",
                            legend.margin = ggplot2::margin(6, 6, 6, 6),
                            legend.background = ggplot2::element_rect(fill = NA))+
              ggplot2::geom_text(data = p_val_data,
                            mapping = ggplot2::aes(x = group2, y = percentage + 1, label = paste0("italic(p)== ", round(p.adjusted, digits = 3))),
                            hjust   = 0.5, vjust = 0, size = 4, parse = T)
  
  bar_plt
}

library(dplyr)
library(ggplot2)
library(lubridate)
library(glmnet)
library(emmeans)
library(knitr)  # required for kable
library(kableExtra)  # required for kableExtra

# required libraries:
# ggstance
# broom.mixed
# jtools
# effects
# sjPlot
# rstatix

super_dir <- "./"
plot.dir <- file.path(super_dir, "plots")
data.dir <- file.path(super_dir, "data")

if(!file.exists(data.dir)){
  dir.create(data.dir, recursive = T)
}

for(plt_format in c("png", "pdf", "svg")){

  if(!file.exists(file.path(plot.dir, plt_format))){
    dir.create(file.path(plot.dir, plt_format), recursive = T)
  }
}

df_coad <- readRDS(file.path(data.dir, "covacser_data_coadministration.rds"))

fisher_test_table(df = df_coad, gr_col = "coadmin", gr_row = "vaccine")

fisher_test_table(df = df_coad, gr_col = "coadmin", gr_row = "gender")

fisher_test_table(df = df_coad, gr_col = "coadmin", gr_row = "smoking")

chisq_test_table(df = df_coad, gr_col = "coadmin", gr_row = "risk")

fisher_test_table(df = df_coad, gr_col = "coadmin", gr_row = "sideeffects")

df_coad %>% ggplot2::ggplot(ggplot2::aes(x = age)) +
  ggplot2::facet_wrap(~coadmin) +
  ggplot2::geom_histogram(bins = 40, col = "black") +
  ggplot2::xlab("Age of subjects") +
  ggplot2::ylab("Number of subjects") +
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text = ggplot2::element_text(size = 12),
                 axis.title = ggplot2::element_text(face = "bold"),
                 strip.text = ggplot2::element_text(face = "bold", size = 12)) +
  ggplot2::geom_text(data = df_coad %>% group_by(coadmin) %>% count(),
                     mapping = ggplot2::aes(x = 45, y = 50, label = paste0("Subjects: ",n)),
                     hjust   = 0.5,
                     vjust   = 0, size = 4) +
  ggplot2::geom_text(data = df_coad %>% group_by(coadmin) %>% summarise(age = median(age)),
                     mapping = ggplot2::aes(x = 45, y = 45, label = paste0("Median Age: ",round(age))),
                     hjust   = 0.5,
                     vjust   = 0, size = 4)

ggsave(filename = file.path(plot.dir,"Histogramm_Altersverteilung.pdf"), device = "pdf", width = 6, height = 5)
ggsave(filename = file.path(plot.dir,"svg/Histogramm_Altersverteilung.svg"), device = "svg", width = 6, height = 5)

df_coad %>% group_by(coadmin, vaccine) %>% count() %>% group_by(coadmin) %>% mutate(n = n/sum(n)) %>%
  ggplot2::ggplot(ggplot2::aes(x = vaccine, y = n, fill = coadmin)) +
  ggplot2::geom_bar(stat = "identity", position = ggplot2::position_dodge()) +
  ggplot2::xlab("Vaccine") +
  ggplot2::ylab("relative frequency of subjects") +
  ggplot2::theme_bw() +
  ggsci::scale_fill_nejm()+
  ggplot2::theme(axis.text = ggplot2::element_text(size = 12, colour = "black"),
                 axis.title = ggplot2::element_text(face = "bold", size = 14),
                 legend.title = ggplot2::element_blank(),
                 legend.text = ggplot2::element_text(size = 12),
                 legend.position = "bottom") +
  ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(shape = 22, size = 8, stroke = 1), ncol = 2))

ggsave(filename = file.path(plot.dir,"png/Barplot_Vaccine.png"), device = "png", width = 6, height = 6, dpi = 300)
ggsave(filename = file.path(plot.dir,"svg/Barplot_Vaccine.svg"), device = "svg", width = 6, height = 6, dpi = 300)

df_coad %>%
  ggplot2::ggplot(ggplot2::aes(y = coadmin, x = date_vacc3, col = vaccine)) +
  ggplot2::geom_point(position=ggplot2::position_jitterdodge(jitter.height = 0.2)) +
  ggplot2::theme_bw() +
  ggplot2::ylab("") +
  ggplot2::xlab("Date of vaccination") +
  ggsci::scale_color_nejm() +
  ggplot2::theme(axis.text = ggplot2::element_text(size = 12, colour = "black"),
                 axis.text.x = ggplot2::element_text(angle = 45, vjust = 1, hjust = 1),
                 axis.title = ggplot2::element_text(face = "bold", size = 14),
                 legend.title = ggplot2::element_blank(),
                 legend.text = ggplot2::element_text(size = 12),
                 legend.position = "bottom") +
  ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(shape = 20, size = 8, stroke = 1), ncol = 2))

ggsave(filename = file.path(plot.dir,"png/Impfdaten_Vakzin_scatter.png"), device = "png", width = 6, height = 6, dpi = 300)
ggsave(filename = file.path(plot.dir,"svg/Impfdaten_Vakzin_scatter.svg"), device = "svg", width = 6, height = 6, dpi = 300)

gr_stats <- rbind(df_coad %>% group_by(vaccine) %>% summarise("Med ± IQR" = med_iqr(time_serum_coadmin)) %>% rename(gr = "vaccine"),
                  df_coad %>% group_by(coadmin) %>% summarise("Med ± IQR" = med_iqr(time_serum_coadmin)) %>% rename(gr = "coadmin"),
                  df_coad %>% group_by(vaccine, coadmin) %>% summarise("Med ± IQR" = med_iqr(time_serum_coadmin)) %>%
                              mutate(gr = paste0(vaccine, ":", coadmin)) %>% ungroup() %>% select(-vaccine, -coadmin))

time_test_table <- rbind(df_coad %>% rstatix::wilcox_test(time_serum_coadmin ~ vaccine) %>% select(-.y.),
      df_coad %>% rstatix::wilcox_test(time_serum_coadmin ~ coadmin) %>% select(-.y.),

df_coad %>% group_by(vaccine) %>% rstatix::wilcox_test(time_serum_coadmin ~ coadmin) %>% 
            mutate(group1 = paste0(vaccine, ":", group1),
                   group2 = paste0(vaccine, ":", group2)) %>% 
            select(-vaccine, -.y.),
df_coad %>% group_by(coadmin) %>% rstatix::wilcox_test(time_serum_coadmin ~ vaccine) %>% 
            mutate(group1 = paste0(group1, ":", coadmin),
                   group2 = paste0(group2, ":", coadmin)) %>%
            select(-coadmin, -.y.)) %>% inner_join(gr_stats, by = c("group1" = "gr")) %>%
            inner_join(gr_stats, by = c("group2" = "gr")) %>% 
  rename("Median(gr1) (IQR)" = `Med ± IQR.x`,
         "Median(gr2) (IQR)" = `Med ± IQR.y`) %>%
  relocate(group1, group2, matches("Med"))

time_test_table <- cbind(time_test_table, p.adjusted = p.adjust(time_test_table$p, method = "holm"))
time_test_table %>% knitr::kable() %>%
  kableExtra::kable_paper() %>%
  kableExtra::pack_rows("Subgroups: ", 3, 6) 

df_coad %>% group_by(coadmin) %>% mutate(IgG = log10(IgG)) %>%
           summarise(geo_mean = 10^(mean(IgG)),
                     minIgG = min(10^IgG),
                     maxIgG = max(10^IgG),
                     IQR_IgG = IQR(10^IgG),
                     `N(IgG > 31.5)` = sum(10^IgG > 31.5),
                     N = n()) %>%
          mutate(sumN = sum(N),
                 `% N/sumN` = N/sumN * 100,
                 `% N(IgG > 31.5)/N` = `N(IgG > 31.5)`/N * 100) %>% 
          mutate(across(.fns = ~ format(x = .x,  digits = 2, nsmall = 2))) %>% 
  dplyr::relocate(coadmin, N, sumN, geo_mean,   minIgG, maxIgG, IQR_IgG, `% N/sumN`,  `N(IgG > 31.5)`, `% N(IgG > 31.5)/N`) %>% knitr::kable() %>%
  kableExtra::kable_paper()

plot_density <- function(){
  
  igg_df <- df_coad %>% select(IgG, logIgG) %>% tidyr::pivot_longer(cols = c(IgG, logIgG), names_to = "transform", values_to = "values")
  theo_params <- igg_df %>% group_by(transform) %>% summarise(mu_est = mean(values),
                                                             sd_est = sd(values),
                                                             xmin = min(values),
                                                             xmax = max(values))
  dens_plts <- list()
  igg_names <- unique(igg_df$transform)

  mains <- c("Density of IgG titres", "Density of logarithmized IgG titres")
  xlabs <- c("IgG", "log10(IgG)")

  m <- rbind(c(0, 0.5, 0, 1), 
             c(0.5, 1, 0, 1))

  split.screen(m)
  
  for(i in seq_along(igg_names)){
  
    screen(i)
  
    if(i == 1){
      par(mar = c(4.1, 4.1, 2.1, .5))
    }
    
    if(i == 2){
      par(mar = c(4.1, 4.1, 2.1, .5))
    }
  
    values <- df_coad %>% pull(!!sym(igg_names[i]))
    x_vals <- seq(min(values), max(values), length.out = 1000)

    hist(values, xlab = xlabs[i], prob = T, col = "white", main = mains[i], breaks = 30)
    lines(x = x_vals, y = dnorm(x = x_vals, mean = mean(values), sd = sd(values)), lty = 1, lwd = 2)
    box()
  }

  close.screen(all.screens = TRUE)

}

png(filename = file.path(plot.dir, paste0("png/Distribution_IgG_titres.png")),units="px", width=2500, height=1000, res=300)
plot_density()
close.screen(all.screens = TRUE)
dev.off()

svg(filename = file.path(plot.dir, paste0("svg/Distribution_IgG_titres.svg")), width=10, height=4)
plot_density()
close.screen(all.screens = TRUE)
dev.off()

plot_density()

df_test <- df_coad %>% select(IgG, logIgG, AgeCategory, gender, coadmin, time_serum_coadmin, smoking, risk, vaccine, sideeffects)

gls_mod <- nlme::gls(logIgG ~ AgeCategory + gender +  smoking + coadmin  + vaccine + time_serum_coadmin + risk + sideeffects,
                      data = df_test,
                      weights = nlme::varPower())

mod_means <- modelbased::get_emmeans(gls_mod, at = c("vaccine", "coadmin")) %>% as_tibble()

int_plt_mm <- mod_means %>% ggplot2::ggplot(ggplot2::aes(x = coadmin, y = 10^emmean, col = vaccine)) + 
                  ggplot2::geom_errorbar(ggplot2::aes(ymax = 10^upper.CL, ymin = 10^lower.CL), width = 0.4) + 
                  ggplot2::geom_point(size = 2) + 
                  ggplot2::geom_line(ggplot2::aes(group = vaccine)) + 
                  ggplot2::scale_y_log10(limits = c(1000, 4000), breaks = c(1000, 2000, 3000, 4000), labels = c(1000, 2000, 3000, 4000)) +
                  ggplot2::annotation_logticks(sides = "l") +
                  ggplot2::xlab("") + 
                  ggplot2::ylab("Anti-SARS-CoV-2-Spike IgG [BAU/ml]") +
                  ggplot2::theme_bw() +
                  ggplot2::theme(axis.text = ggplot2::element_text(size = 14, colour = "black"),
                                 axis.title = ggplot2::element_text(face = "bold", size = 12),
                                 legend.title = ggplot2::element_blank(),
                                 legend.text = ggplot2::element_text(size = 16),
                                 legend.position = "bottom")

int_plt_mm

ggsave(plot = int_plt_mm, filename = file.path(plot.dir,"png/Interaction_marginal_means_1.png"), device = "png", width = 6, height = 6, dpi = 300)
ggsave(plot = int_plt_mm, filename = file.path(plot.dir,"svg/Interaction_marginal_means_1.svg"), device = "svg", width = 6, height = 6, dpi = 300)

model_residuals <- as.numeric(residuals(gls_mod))
dw_stat <- car::durbinWatsonTest(model_residuals)

# Calculate P-value based on permutation
n_obs <- length(model_residuals)
n_samplings <- 1000
mu <- fitted.values(gls_mod)
X_model <- model.matrix(logIgG ~  AgeCategory + gender +  smoking + coadmin  + vaccine + time_serum_coadmin + risk + sideeffects, 
                        data = df_test)

sampled_data <-  matrix(sample(model_residuals, n_obs * n_samplings, replace = TRUE), n_obs, n_samplings) + matrix(mu, n_obs, n_samplings)

E <- residuals(nlme::gls(sampled_data ~ X_model - 1, weights = nlme::varPower()))
dw_sampled <- apply(E, 2, car::durbinWatsonTest)

p_value_dw <- (sum(dw_sampled > dw_stat))/n_samplings
p_value_dw <- 2 * (min(p_value_dw, 1 - p_value_dw))

hist(dw_sampled, breaks = 50, freq = T, main = "Histogram sampled DW statistics", xlab = "Sampled DW values", col = "white")
abline(v = dw_stat, lwd = 3, col = "red")

x_vals <- seq(min(model_residuals), max(model_residuals), length.out = 1000)

hist(model_residuals, xlab = "Residuals", prob = T, col = "white", main = expression('Density of residuals'), breaks = 30)
lines(x = x_vals, y = dnorm(x = x_vals, mean = mean(model_residuals), sd = sd(model_residuals)), lty = 1, lwd = 2)

png(filename = file.path(plot.dir, "png/histogram_residuals.png"), units="px", width=1600, height=1600, res=300)
hist(model_residuals, xlab = "Residuals", prob = T, col = "white", main = expression('Density of residuals'), breaks = 30)
lines(x = x_vals, y = dnorm(x = x_vals, mean = mean(model_residuals), sd = sd(model_residuals)), lty = 1, lwd = 2)
dev.off()

svg(filename = file.path(plot.dir, "svg/histogram_residuals.svg"), width=6, height=6)
hist(model_residuals, xlab = "Residuals", prob = T, col = "white", main = expression('Density of residuals'), breaks = 30)
lines(x = x_vals, y = dnorm(x = x_vals, mean = mean(model_residuals), sd = sd(model_residuals)), lty = 1, lwd = 2)
dev.off()

data.frame(fit = mu, residuals = model_residuals) %>% 
  ggplot2::ggplot(aes(x = fit, y = residuals)) + 
  ggplot2::geom_point(alpha = 0.75, pch = 20, size = 2) + 
  ggplot2::xlab("Fitted values based on linear model") + 
  ggplot2::ylab("Residuals from linear model") + 
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text = ggplot2::element_text(size = 14, colour = "black"),
                 axis.title = ggplot2::element_text(face = "bold", size = 12),
                 legend.title = ggplot2::element_blank(),
                 legend.text = ggplot2::element_text(size = 16),
                 legend.position = "bottom")

ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "png", "Regression_model_residual_variances.png"), device = "png",  width = 10, height = 8, dpi = 300)
ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "svg", "Regression_model_residual_variances.svg"), device = "svg",  width = 10, height = 8, dpi = 300)

sum_gls <- summary(gls_mod)
glm_p_values <- cbind(sum_gls$tTable, p.adj = p.adjust(sum_gls$tTable[,4], method = "BY"))

anova_gls <- as.data.frame(anova(gls_mod))
anova_gls <- cbind(anova_gls, p.adj = p.adjust(anova_gls$`p-value`, method = "BY"))
anova_gls <- anova_gls %>% mutate(p_plot = dplyr::case_when(p.adj < 1e-4 ~ "< 0.0001",
                                             p.adj > 1e-4 &  p.adj < 1e-3  ~"< 0.001",
                                             p.adj > 1e-3 &  p.adj < 1e-2  ~ "< 0.01",
                                             p.adj > 1e-2 &  p.adj < 0.05  ~ "< 0.05",
                                             TRUE ~ format(round(p.adj, 2))
                                             ))

y_hat_igg <- predict(gls_mod, interval = "prediction") 
rsquare_cv <- cor(df_test$logIgG, y_hat_igg)^2

anova_gls %>% knitr::kable() %>%
  kableExtra::kable_paper()

entities <- c("AgeCategory", "gender", "coadmin", "vaccine", "smoking", "sideeffects", "risk", "vaccine + coadmin")
marg_list <- list()
for(i in seq_along(entities)){
  marginal <- emmeans(gls_mod, eval(parse( text = paste0("~ ", entities[i]))), rg.limit = 160000, weights = "proportional")
  marg_list[[i]] <- as.data.frame(pairs(marginal))
}

names(marg_list) <- entities

post_hoc_comp <- do.call("rbind", marg_list)

post_hoc_comp$contrast <- gsub(pattern = "\\(|\\)", replacement = "", x = post_hoc_comp$contrast)
post_hoc_comp <- cbind(post_hoc_comp, p.adj = p.adjust(post_hoc_comp$p.value, method = "BY"))

risk_labels <- paste0("Be at Risk: ", apply(combn(x = c("False", "Rather not true ", "Partly", "Rather true", "True"), m = 2),  2, function(col_el){
  paste0(col_el, collapse = " vs. ")
}))

coadmin_labels <- paste0(apply(combn(x = c("BNT162b2mRNA:Coadministration", "mRNA-1273:Coadministration", "BNT162b2mRNA:No Coadministration", "mRNA-1273:No Coadministration"), m = 2),  2, function(col_el){
  paste0(col_el, collapse = " vs. ")
}))

contrast_labels <- c("Age (18-40 vs. ≥41)", "Gender (female vs. male)", "Coadministration (yes vs. no)", "Vaccine (BNT162b2mRNA vs. mRNA-1273)", "Smoking (yes vs. no)", "Side effects (yes vs. no)", risk_labels, coadmin_labels)
post_hoc_comp <- cbind(comparison = contrast_labels, post_hoc_comp)

post_hoc_comp %>% select(-contrast) %>% knitr::kable(row.names = F) %>%
  kableExtra::kable_paper() %>%
  kableExtra::scroll_box(width = "1000px")

var_factors <- list("AgeCategory", "gender", "coadmin", "vaccine", "smoking", "risk", "sideeffects", c("vaccine","coadmin"))

gr_stat <- do.call("rbind", lapply(var_factors, function(var_fac){
  
  df_sum <- c()
  if(length(var_fac) == 1){

    df_sum <- df_test %>% group_by(!!sym(var_fac)) %>% summarise(value = med_iqr(IgG),
                                                                 Min = min(IgG),
                                                                 Max = max(IgG),
                                                                 N = n()) %>% 
                          rename(gr = !!sym(var_fac))
    
  }else{
    
    df_sum <- df_test %>% group_by(vaccine, coadmin) %>% summarise(value = med_iqr(IgG),
                                                                   Min = min(IgG),
                                                                   Max = max(IgG),
                                                                   N = n()) %>% 
                          mutate(gr = paste0(vaccine, " ", coadmin)) %>% ungroup() %>% dplyr::select(gr, value, Min, Max, N)
  }
  
  df_sum
}))

 
gr_mat <- do.call("rbind", strsplit(x = as.character(post_hoc_comp$contrast), split =  " - "))
colnames(gr_mat) <- c("gr1", "gr2")

df_table_igg <- cbind(post_hoc_comp, gr_mat) %>% dplyr::inner_join(gr_stat, by = c("gr1" = "gr")) %>% dplyr::inner_join(gr_stat, by = c("gr2" = "gr")) %>% 
                mutate(p.short = dplyr::case_when(p.value >= 0.01 ~ as.character(round(p.value, 2)),
                                                  p.value < 0.0001 ~ "< 0.0001",
                                                  p.value < 0.001 ~ "< 0.001",
                                                  p.value < 0.01 ~ "< 0.01"),
                       p.ad.short = dplyr::case_when(p.adj >= 0.01 ~ as.character(round(p.adj, 2)),
                                                  p.adj < 0.0001 ~ "< 0.0001",
                                                  p.adj < 0.001 ~ "< 0.001",
                                                  p.adj < 0.01 ~ "< 0.01"),
                       across(.cols = where(is.numeric), .fns = function(x) format(round(x, digits = 2), 
                                                                                    digits = 2, nsmall = 2, 
                                                                                    format = "f", scientific = F)),
                       gr1 = if_else(grepl(x = comparison, pattern = "Side"), true = "side effects", false = gr1),
                       gr2 = if_else(grepl(x = comparison, pattern = "Side"), true = "no side effects", false = gr2))

df_table_igg <- df_table_igg %>% rename(p.adjusted = p.adj,
                                        "Median(Gr1) ± IQR" = value.x,
                                        "Median(Gr2) ± IQR" = value.y,
                                        "Min(Gr1)" = Min.x,
                                        "Min(Gr2)" = Min.y,
                                        "Max(Gr1)" = Max.x,
                                        "Max(Gr2)" = Max.y,
                                        "N(Gr1)" = N.x,
                                        "N(Gr2)" = N.y,
                                        "Group 1" = gr1,
                                        "Group 2" = gr2) %>% 
                relocate("Group 1", "Group 2", "N(Gr1)", "Min(Gr1)", "Max(Gr1)", "Median(Gr1) ± IQR", "N(Gr2)", "Min(Gr2)", "Max(Gr2)", "Median(Gr2) ± IQR", p.value, p.adjusted, p.short, p.ad.short) %>% 
                dplyr::select(-comparison, -contrast, -estimate, -SE, -df, -t.ratio)

df_table_igg %>% knitr::kable() %>%
  kableExtra::kable_paper() %>% 
  kableExtra::pack_rows("Be at risk: ", 7, 16) %>%
  kableExtra::pack_rows("Vaccine third doses and Coadministration: ", 17, 22) %>%
  kableExtra::add_header_above(c(" " = 2, "Group 1 (IgG statistics)" = 4, "Group 2 (IgG statistics)" = 4, "Group 1 vs. Group 2" = 4)) %>%
  kableExtra::scroll_box(width = "1400px")

post_hoc_comp <- post_hoc_comp[rev(seq_len(nrow(post_hoc_comp))),]
post_hoc_comp$comparison <- factor(post_hoc_comp$comparison, levels = post_hoc_comp$comparison) 

post_hoc_comp %>% ggplot2::ggplot() +
  ggplot2::geom_point(ggplot2::aes(y = comparison, x = estimate)) +
  ggplot2::geom_errorbar(aes(y = comparison, xmin=estimate-SE, xmax=estimate+SE), width=.2,
                         position=position_dodge(.9)) +
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text = ggplot2::element_text(size = 12, colour = "black"),
                 axis.title = ggplot2::element_text(face = "bold", size = 12),
                 legend.title = ggplot2::element_blank(),
                 legend.text = ggplot2::element_text(size = 12),
                 legend.position = "bottom") +
  ggplot2::xlab("Effect on IgG (estimated coefficient difference)") + 
  ggplot2::ylab("Comparison")

ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "png", "Regression_model_mean_se_estimates.png"), device = "png",  width = 10, height = 12, dpi = 300)
ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "svg", "Regression_model_mean_se_estimates.svg"), device = "svg",  width = 10, height = 12, dpi = 300)

# Coadmin
gr_matrix <- as.data.frame(t(combn(levels(df_test$coadmin), 2)))
colnames(gr_matrix) <- c("group1", "group2")
                           
gr_count <- df_test %>% group_by(coadmin) %>% count() 

stat.comp <- gr_matrix %>% full_join(gr_count, by = c("group1" = "coadmin"), suffix = c("1", "2")) %>% 
             inner_join(gr_count, by = c("group2" = "coadmin"), suffix = c("1", "2")) %>% 
             mutate(contrast = paste0(group1, " - ", group2)) %>% 
             inner_join(post_hoc_comp %>% select(contrast, p.value, p.adj)) %>% mutate(".y." = "IgG")

stat.comp <- stat.comp %>% rstatix::add_significance() %>% filter(p.adj < 0.05) %>% rstatix::add_x_position()
box_upper <- df_test %>% group_by(coadmin) %>% summarise(igg_max = quantile(x = 10^logIgG, probs = 0.995))
stat.comp <- stat.comp %>% inner_join(box_upper, by = c("group1" = "coadmin")) %>% rename(ygr1 = igg_max) %>% inner_join(box_upper, by = c("group2" = "coadmin"))  %>% rename(ygr2 = igg_max)

stat.comp <- stat.comp %>% group_by(group1, group2) %>%
  mutate(diff_comp = xmax - xmin,
         ymax = round(log10(max(ygr1, ygr2)), 2)) %>%
  arrange(diff_comp, ymax) %>% group_by(diff_comp, ymax) %>%
  mutate(id = cur_group_id()) %>% as.data.frame()


# Vaccine + Coadmin
gr_matrix <- as.data.frame(t(combn(x = c("BNT162b2mRNA Coadministration", "mRNA-1273 Coadministration", "BNT162b2mRNA No Coadministration", "mRNA-1273 No Coadministration"), m = 2)))
colnames(gr_matrix) <- c("group1", "group2")

gr_count <- df_test %>% mutate(`vaccine:coadmin` = paste0(vaccine," ", coadmin)) %>% group_by(`vaccine:coadmin`) %>% count() 

stat.comp <- gr_matrix %>% full_join(gr_count, by = c("group1" = "vaccine:coadmin"), suffix = c("1", "2")) %>% 
             inner_join(gr_count, by = c("group2" = "vaccine:coadmin"), suffix = c("1", "2")) %>% 
             mutate(contrast = paste0(group1, " - ", group2)) %>% 
             inner_join(post_hoc_comp %>% select(contrast, p.value, p.adj)) %>% mutate(".y." = "IgG")

stat.comp <- stat.comp %>% rstatix::add_significance() %>% filter(p.adj < 0.05) %>% rstatix::add_x_position()
box_upper <- df_test %>% mutate(`vaccine:coadmin` = paste0(vaccine," ", coadmin)) %>% group_by(`vaccine:coadmin`) %>% summarise(igg_max = quantile(x = 10^logIgG, probs = 0.995))
stat.comp <- stat.comp %>% inner_join(box_upper, by = c("group1" = "vaccine:coadmin")) %>% rename(ygr1 = igg_max) %>% inner_join(box_upper, by = c("group2" = "vaccine:coadmin"))  %>% rename(ygr2 = igg_max)

stat.comp <- stat.comp %>% group_by(group1, group2) %>%
  mutate(diff_comp = xmax - xmin,
         ymax = round(log10(max(ygr1, ygr2)), 2)) %>%
  arrange(diff_comp, ymax) %>% group_by(diff_comp, ymax) %>%
  mutate(id = cur_group_id()) %>% as.data.frame()

stat.comp <- stat.comp %>% mutate(xmin = case_when(xmin <= 2 ~  0.8 + (xmin-1) * 0.4,
                                     xmin > 2 ~   1 + (xmin-1) * 0.4),
                    xmax = case_when(xmax <= 2 ~  0.8 + (xmax-1) * 0.4,
                                     xmax > 2 ~   1 + (xmax-1) * 0.4))

bx_plt <- df_test %>% mutate(`vaccine:coadmin` = paste0(vaccine," ", coadmin)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = vaccine, y = 10^logIgG, fill = `vaccine:coadmin`)) +
  ggplot2::stat_boxplot(geom = "errorbar", lwd = 1, show.legend = FALSE, width = 0.4, position = position_dodge(width = 0.8)) +
  ggplot2::geom_boxplot(aes(x = vaccine, y = 10^logIgG, fill = `vaccine:coadmin`), width = 0.8, color = "black", lwd = 1, show.legend = FALSE, outlier.size = 0) +
  ggbeeswarm::geom_quasirandom(dodge.width=.75, alpha = 0.8,  pch = 20, size = 0.5) +
  ggplot2::theme_bw() +
  ggplot2::scale_fill_manual(labels = rep(c("Coadministration", "No Coadministration"), times = 2), values = c("#42B540FF", "darkgreen", "#ED0000FF", "darkred")) + 
  ggplot2::theme(axis.text = ggplot2::element_text(size = 17, colour = "black"),
                 axis.title = ggplot2::element_text(face = "bold", size = 16),
                 legend.title = ggplot2::element_blank(),
                 legend.text = ggplot2::element_text(size = 16),
                 legend.position = "bottom") +
  ggplot2::ylab("Anti-SARS-CoV-2-Spike IgG [BAU/ml]") +
  ggplot2::scale_y_log10(limits = c(10, 70000), breaks = 10^(0:4), labels = format(10^(0:4), scientific = F)) +
  ggplot2::annotation_logticks(sides = "l") +
  ggplot2::xlab("") +
  # ggplot2::geom_hline(yintercept = 31.5, lty = 2) +
  ggsignif::geom_signif(xmin = stat.comp$xmin,
                        xmax = stat.comp$xmax,
                        y_position = 0.1 * 1.4^(seq_along(stat.comp$ymax)) + (stat.comp$ymax),
                        annotation = stat.comp$p.adj.signif,
                        tip_length = 0.01, vjust = 0.75, lwd = 0.75, textsize = 7) + 
  ggplot2::guides(fill = ggplot2::guide_legend(override.aes = list(shape = 22, size = 8, stroke = 1), ncol = 2))

bx_plt

ggplot2::ggsave(plot = bx_plt, filename = file.path(plot.dir, "png", "Boxplot_coadministration_comparison.png"), device = "png",  width = 8, height = 6, dpi = 300)

ggplot2::ggsave(plot = bx_plt, filename = file.path(plot.dir, "svg", "Boxplot_coadministration_comparison.svg"), device = "svg",  width = 8, height = 6, dpi = 300)

plot_vars <- c("AgeCategory", "gender", "smoking", "sideeffects", "risk")

tick_labels <- list(Age = c("18-40", "41+"),
                    Gender = c("female", "male"),
                    Smoking = c("smoker", "non-smoker"),
                    "Side effects" =  c("yes", "no"),
                    "Be at Risk" = c("False", "Rather not true ", "Partly", "Rather true", "True")
                    )
plot_list <- list()

for(i in seq_along(plot_vars)){
  
    gr_matrix <- as.data.frame(t(combn(levels(df_test[[plot_vars[i]]]), 2)))
    colnames(gr_matrix) <- c("group1", "group2")
                           
    gr_count <- df_test %>% group_by(!!sym(plot_vars[i])) %>% count() 

    stat.comp <- gr_matrix %>% full_join(gr_count, by = c("group1" = plot_vars[i]), suffix = c("1", "2")) %>% 
                 inner_join(gr_count, by = c("group2" = plot_vars[i]), suffix = c("1", "2")) %>% 
                 mutate(contrast = paste0(group1, " - ", group2)) %>% 
                 inner_join(post_hoc_comp %>% select(comparison, contrast, p.value, p.adj) %>% filter(grepl(pattern = names(tick_labels)[i], x = comparison))) %>% mutate(".y." = "IgG")
    
    stat.comp <- stat.comp %>% rstatix::add_significance() %>% filter(p.adj < 0.05) %>% rstatix::add_x_position()
    box_upper <- df_test %>% group_by(!!sym(plot_vars[i])) %>% summarise(igg_max = quantile(x = IgG, probs = 0.995))

    stat.comp <- stat.comp %>% inner_join(box_upper, by = c("group1" = plot_vars[i])) %>% rename(ygr1 = igg_max) %>% 
      inner_join(box_upper, by = c("group2" = plot_vars[i]))  %>% rename(ygr2 = igg_max)

    stat.comp <- stat.comp %>% group_by(group1, group2) %>%
                 mutate(diff_comp = xmax - xmin,
                        ymax = round(log10(max(ygr1, ygr2)), 2)) %>%
                 arrange(diff_comp, ymax) %>% group_by(diff_comp, ymax) %>%
                 mutate(id = cur_group_id()) %>% as.data.frame()

    cat('### Comparison: ', unique(as.character(stat.comp$comparison)), '\n')
    
    bx_plt <- df_test %>% ggplot2::ggplot(ggplot2::aes_string(x = plot_vars[i], y = "IgG")) +
      ggplot2::stat_boxplot(geom = "errorbar", lwd = 1, width = 0.5, show.legend = FALSE) +
      ggplot2::geom_boxplot(aes_string(fill = plot_vars[i]), color = "black", lwd = 1, show.legend = FALSE, outlier.size = 0) +
      ggbeeswarm::geom_quasirandom(alpha = 0.8, width = 0.4, pch = 20, size = 0.5) +
      ggplot2::scale_y_log10() +
      ggplot2::theme_bw() +
      ggsci::scale_fill_lancet() +
      ggplot2::scale_x_discrete(labels = tick_labels[[i]])+
      ggplot2::theme(axis.text = ggplot2::element_text(size = 14, colour = "black"),
                     axis.title = ggplot2::element_text(face = "bold", size = 12),
                     legend.title = ggplot2::element_blank(),
                     legend.text = ggplot2::element_text(size = 12),
                     legend.position = "none") +
      ggplot2::xlab(names(tick_labels)[i]) +
      ggplot2::ylab("Anti-SARS-CoV-2-Spike IgG [BAU/ml]") +
      ggplot2::annotation_logticks(sides = "l") 
      # ggplot2::geom_hline(yintercept = 31.5, lty = 2)

    if(dim(stat.comp)[1] > 0 ){

      box_upper <- df_test %>% group_by(!!sym(plot_vars[i])) %>% summarise(igg_max = quantile(x = IgG, probs = 0.995))

      
      if(df_test %>% pull(!!sym(plot_vars[i])) %>% unique() %>% length() > 2){

        stat.comp <- stat.comp %>% group_by(group1, group2) %>%
          mutate(diff_comp = xmax - xmin,
                 ymax = round(log10(max(ygr1, ygr2)), 2),
                 p_val = p.adj.signif) %>%
          arrange(diff_comp, ymax) %>% group_by(diff_comp, ymax) %>%
          mutate(id = cur_group_id()) %>% as.data.frame()
      }else{

        stat.comp <- stat.comp %>% group_by(group1, group2) %>%
          mutate(diff_comp = xmax - xmin,
                 ymax = round(log10(max(ygr1, ygr2)), 2),
                 p_val = p.adj.signif) %>%
          arrange(diff_comp, ymax) %>% group_by(diff_comp, ymax) %>%
          mutate(id = cur_group_id()) %>% as.data.frame()

      }

      bx_plt <-  bx_plt + ggsignif::geom_signif(xmin = stat.comp$group1,
                                                xmax = stat.comp$group2,
                                                y_position = 0.175 * seq_along(stat.comp$ymax) + (stat.comp$ymax),
                                                annotation = stat.comp$p.adj.signif,
                                                tip_length = 0.01, vjust = 0.75, lwd = 0.75, textsize = 7)
    }

    plot_list[[i]] <- bx_plt
}

## Generate combined plot for confounding factors associated to IgG titres

plot_list[[1]] <- plot_list[[1]] + 
                  ggplot2::xlab("") + 
                  ggsci::scale_fill_jco() + 
                  ggplot2::facet_wrap(~., labeller = ggplot2::as_labeller(c(`(all)` = "Age [Years]"))) + 
                  ggplot2::theme(axis.title.x = ggplot2::element_blank(),
                                 strip.text = ggplot2::element_text(size = 14, face = "bold"), 
                                 strip.background = ggplot2::element_rect(fill = "white"))
plot_list[[2]] <- plot_list[[2]] + 
                  ggplot2::xlab("") + 
                  ggplot2::facet_wrap(~., labeller = ggplot2::as_labeller(c(`(all)` = "Gender"))) + 
                  ggplot2::theme(axis.title = ggplot2::element_blank(),
                                 strip.text = ggplot2::element_text(size = 14, face = "bold"), 
                                 strip.background = ggplot2::element_rect(fill = "white"))
plot_list[[3]] <- plot_list[[3]] + 
                  ggplot2::xlab("") + 
                  ggsci::scale_fill_jama() + 
                  ggplot2::facet_wrap(~., labeller = ggplot2::as_labeller(c(`(all)` = "Smoking"))) + 
                  ggplot2::theme(axis.title = ggplot2::element_blank(),
                                 strip.text = ggplot2::element_text(size = 14, face = "bold"), 
                                 strip.background = ggplot2::element_rect(fill = "white"))
plot_list[[4]] <- plot_list[[4]] + 
                  ggplot2::xlab("") + 
                  ggsci::scale_fill_lancet() + 
                  ggplot2::facet_wrap(~., labeller = ggplot2::as_labeller(c(`(all)` = "Side effects"))) + 
                  ggplot2::theme(axis.title = ggplot2::element_blank(),
                                 strip.text = ggplot2::element_text(size = 14, face = "bold"), 
                                 strip.background = ggplot2::element_rect(fill = "white"))
plot_list[[5]] <- plot_list[[5]] + 
                  ggplot2::xlab("") + 
                  ggsci::scale_fill_jco() + 
                  ggplot2::facet_wrap(~., labeller = ggplot2::as_labeller(c(`(all)` = "Be at Risk"))) + 
                  ggplot2::theme(axis.title.x= ggplot2::element_blank(),
                                 strip.text = ggplot2::element_text(size = 14, face = "bold"), 
                                 strip.background = ggplot2::element_rect(fill = "white"))

comb_plt_top <- cowplot::plot_grid(plotlist = plot_list[1:3], ncol = 3, rel_widths = c(0.31, 0.3, 0.3), labels = LETTERS[1:3], label_size = 16)
comb_plt_bottom <- cowplot::plot_grid(plotlist = plot_list[5:4], ncol = 2, rel_widths = c(0.61, 0.3), labels = LETTERS[4:5], label_size = 16)


comb_plt <- cowplot::plot_grid(comb_plt_top, comb_plt_bottom, nrow = 2, ncol = 1)

comb_plt

ggplot2::ggsave(plot = comb_plt, filename = file.path(plot.dir, "png", "Boxplot_other_categorical_factors.png"), device = "png",  width = 12, height = 8, dpi = 300)
ggplot2::ggsave(plot = comb_plt, filename = file.path(plot.dir, "png/Combined_boxplotsl.png"), 
                device = "png",  width = 12, height = 8, dpi = 300)

ggplot2::ggsave(plot = comb_plt, filename = file.path(plot.dir, "svg", "Boxplot_other_categorical_factors.svg"), device = "svg",  width = 12, height = 8, dpi = 300)
ggplot2::ggsave(plot = comb_plt, filename = file.path(plot.dir, "svg/Combined_boxplotsl.svg"), 
                device = "svg",  width = 12, height = 8, dpi = 300)

pred_int <- predict(gls_mod, interval = "prediction") 
pred_time_igg <- cbind(df_test %>% select(logIgG, time_serum_coadmin), Prediction = pred_int)

y_hat_igg <- predict(gls_mod, interval = "prediction") 
rsquare_cv <- cor(df_test$logIgG, y_hat_igg)^2

cor_time_igg <- round(cor(df_test$time_serum_coadmin, df_test$logIgG, method = "spearman"), 2)

line_plt_time <- df_test %>%
  ggplot2::ggplot(ggplot2::aes(x = time_serum_coadmin, y = IgG)) +
  ggplot2::geom_point(size = 1, alpha = 0.75) +
  ggplot2::scale_y_log10(limits = c(10, 30000), expand = c(0,0)) +
  ggplot2::geom_smooth(data = pred_time_igg, ggplot2::aes(x = time_serum_coadmin, y = 10^Prediction), method = "lm") +
  ggsci::scale_color_lancet() +
  ggplot2::scale_x_continuous(breaks = seq(0,90,15), labels = seq(0,90,15)) +
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text = ggplot2::element_text(size = 14, colour = "black"),
                 axis.title = ggplot2::element_text(face = "bold", size = 12),
                 legend.title = ggplot2::element_blank(),
                 legend.text = ggplot2::element_text(size = 16),
                 legend.position = "bottom") +
  ggplot2::xlab("Days since vaccination") +
  ggplot2::ylab("Anti-SARS-CoV-2-Spike IgG [BAU/ml]") +
  ggplot2::annotation_logticks(sides = "l") +
  ggplot2::annotate(geom = "text", label = paste0("italic(c)== ", cor_time_igg), x = 50, y = 10000, cex = 5, parse = TRUE, hjust = 0) + 
  ggplot2::annotate(geom = "text", label = paste0("italic(p) ", anova_gls$p_plot[7]), x = 50, y = 13000, cex = 5, parse = TRUE, hjust = 0) 

line_plt_time

ggplot2::ggsave(plot = line_plt_time, filename = file.path(plot.dir, "png/Scatterplot_time_serum.png"), 
                device = "png",  width = 8, height = 6, dpi = 300)
ggplot2::ggsave(plot = line_plt_time, filename = file.path(plot.dir, "svg/Scatterplot_time_serum.svg"), 
                device = "svg",  width = 8, height = 6, dpi = 300)

fisher_tests_coad <- fisher_tests_gr(df = df_coad, gr1 = "coadmin", var_expr = "se[1-6]")

p.adj_fisher_tests_coad <- cbind(fisher_tests_coad[,-c(9,10)], "p.adjusted" = round(p.adjust(fisher_tests_coad$p.value, "BY"), 2))

knitr::kable(p.adj_fisher_tests_coad, align = "c") %>%
    kableExtra::kable_styling(full_width = F, position = "left")%>%
    kableExtra::column_spec(1, bold = T) %>%
    kableExtra::collapse_rows(columns = 4:9, valign = "middle")

fisher_plot_odds(fisher_test_df = fisher_tests_coad, ylab_add = "coadministration")

ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "png/Side_effects_OR_only_coadmin.png"), 
                device = "png",  width = 8, height = 6, dpi = 300)

ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "svg/Side_effects_OR_only_coadmin.svg"), 
                device = "svg",  width = 8, height = 6, dpi = 300)

barplot_sideeffects_percentage(p.adj_fisher_tests_df = p.adj_fisher_tests_coad)

ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "png/Side_effects_Barplot_percentage_only_coadmin.png"), 
                device = "png",  width = 8, height = 6, dpi = 300)

ggplot2::ggsave(plot = last_plot(), filename = file.path(plot.dir, "svg/Side_effects_Barplot_percentage_only_coadmin.svg"), 
                device = "svg",  width = 8, height = 6, dpi = 300)

fisher_tests_sep <- sapply(levels(df_coad$vaccine), function(vacc_name){
                          fisher_tests_gr(df = df_coad %>% filter(vaccine %in% vacc_name), gr1 = "coadmin", var_expr = "se[1-6]")
                    }, simplify = F)

df_vacc_sep <- data.frame(super_group = rep(levels(df_coad$vaccine), each = dim(fisher_tests_sep[[1]])[1]), do.call("rbind", fisher_tests_sep))
df_vacc_sep_p_adj <- cbind(df_vacc_sep, "p.adjusted" = round(p.adjust(df_vacc_sep$p.value, "BY"), 2))

knitr::kable(df_vacc_sep_p_adj[,-c(1,10,11)], align = "l", row.names = F) %>%
    kableExtra::collapse_rows(columns = 6:9) %>%
    kableExtra::pack_rows("BNT162b2mRNA", 1, 12) %>%
    kableExtra::pack_rows("mRNA-1273", 13, 24) %>%
    kableExtra::kable_styling(full_width = F, position = "left")%>%
    kableExtra::column_spec(1, bold = T)