# my.data <- power.data.Pcytokine.p1
# group.factor <- "TYPE"
# facet.factor <- "TYPE"
# time.points <- unique(power.data.Pcytokine.p1$TIMELINE)
# time.points.var <- "TIMELINE"
# x.val <- "VALUE"
# id.var <- "ID"
# sym("TIMELINE")

#############################################
ggBoxPlotPaired <- function(my.data,
                      group.factor,
                      facet.factor,
                      time.points.var,
                      time.points,
                      x.val,
                      id.var,
                      box.plot.title) {
  require(ggplot2); require(dplyr); require(ggpubr)
  # conduct paired t test
  paired.t.test.results <- my.data %>%
    # omit missing values for getting the p-values
    # na.omit() %>%
    # Group by both group.factor and id.var
    group_by(.data[[group.factor]], .data[[id.var]]) %>%
    # Keep only subjects with BOTH time points (baseline and six_month)
    filter(n() == 2) %>%
    # Regroup by the original grouping factor
    ungroup() %>%
    group_by(.data[[group.factor]]) %>%
    summarise(
      p_value = t.test(
        .data[[x.val]][.data[[time.points.var]] == time.points[1]], 
        .data[[x.val]][.data[[time.points.var]] == time.points[2]], 
        paired = TRUE
      )$p.value
    ) %>%
    mutate(
      ##### Add if statement to print p = or p<0.001
      label = paste("p =", format.pval(p_value, digits = 2))
    )
  
  paired.wilcox.test.results <- my.data %>%
    # omit missing values for getting the p-values
    # na.omit() %>%
    # Group by both group.factor and id.var
    group_by(.data[[group.factor]], .data[[id.var]]) %>%
    # Keep only subjects with BOTH time points (baseline and six_month)
    filter(n() == 2) %>%
    # Regroup by the original grouping factor
    ungroup() %>%
    group_by(.data[[group.factor]]) %>%
    summarise(
      p_value = wilcox.test(
        .data[[x.val]][.data[[time.points.var]] == time.points[1]], 
        .data[[x.val]][.data[[time.points.var]] == time.points[2]], 
        paired = TRUE
      )$p.value
    ) %>%
    mutate(
      label = paste("p =", format.pval(p_value, digits = 2))
    )
  
  boxplot.1 <- ggplot(my.data, 
                      aes(x = .data[[time.points.var]], 
                          y = .data[[x.val]])) +
    # Paired lines (gray)
    geom_line(aes(group = .data[[id.var]])
              , color = "gray", alpha = 1,linetype="dashed",
              linewidth=1.5) +
    # add box plot
    geom_violin(aes(fill=.data[[time.points.var]]),
                 alpha=0.5) +
    # Points colored by time point
    geom_point(aes(color = .data[[time.points.var]]), size = 2.5) +
    # Split by test type
    facet_wrap(~.data[[facet.factor]], scales = "free_y") +
    # Customize
    labs(
      title = paste0(box.plot.title," with t (Paired) test"),
      x = "Time Point",
      y = "Value"
    ) +
    # Add p-value annotations
    geom_text(
      data = paired.t.test.results,
      aes(x = 1.5, 
          y = Inf, 
          label = label),  # Center between time points
      vjust = 5, size = 5, color = "black"
    ) +
    theme_bw() +
    theme(legend.position = "top",
          strip.text = element_text(size=16)) 
  
  boxplot.2 <- ggplot(my.data, 
                      aes(x = .data[[time.points.var]], 
                          y = .data[[x.val]])) +
    # Paired lines (gray)
    geom_line(aes(group = .data[[id.var]])
              , color = "gray", alpha = 1,linetype="dashed",
              linewidth=1.5) +
    # add box plot
    geom_violin(aes(fill=.data[[time.points.var]]),
                alpha=0.5) +
    # Points colored by time point
    geom_point(aes(color = .data[[time.points.var]]), size = 2.5) +
    # Split by test type
    facet_wrap(~.data[[facet.factor]], scales = "free_y") +
    # Customize
    labs(
      title = paste0(box.plot.title," with Wilcoxon (Paired) test"),
      x = "Time Point",
      y = "Value"
    ) +
    # Add p-value annotations
    geom_text(
      data = paired.wilcox.test.results,
      aes(x = 1.5, 
          y = Inf, 
          label = label),  # Center between time points
      vjust = 5, size = 5, color = "black"
    ) +
    theme_bw() +
    theme(legend.position = "top",
          strip.text = element_text(size=16)) 
  
  # return a plot
  res.list <- list(boxplot.t = boxplot.1,
                   boxplot.w = boxplot.2)
  return(res.list)
}