library("clinicalsignificance")
library("tidyverse")
library("patchwork")

options(prompt = "R> ", continue = "+  ", width = 70,
        useFancyQuotes = FALSE)


# Custom theming for plots in the introduction. This uses the thematic package.
# If it is not installed, use the following line to install it:
# remotes::install_github("rstudio/thematic")
thematic::thematic_on(
  bg = "white",
  fg = "grey20",
  sequential = thematic::sequential_gradient(fg_weight = 0, bg_weight = 0.9, fg_low = FALSE),
  qualitative = c("#004358", "#BEDB39", "#FFE11A", "#FD7400", "#1F8A70"),
  font = thematic::font_spec("Roboto")
)


# Set ggplot2's light theme as the default theme
theme_set(
  theme_light() + 
    theme(
      legend.position = "bottom",
      panel.grid.minor = element_blank()
    )
  )




# Introduction ----
## Figure 1 ====
fig_group_anchor_categories <- tibble(
  effect = LETTERS[1:5],
  mean = c(1, 3, 5, 7, 9)
) |>
  mutate(
    ci_low = mean - 2,
    ci_high = mean + 2
  ) |>
  ggplot(aes(effect, mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2) +
  geom_hline(yintercept = 6, lty = "dotted") +
  geom_hline(yintercept = 0) +
  labs(x = "Effect Category", y = "Intervention Effect") +
  theme(
    panel.grid = element_blank()
  )
fig_group_anchor_categories
ggsave("../Figures/fig-anchor-effect.pdf", plot = fig_group_anchor_categories, width = 14, height = 8, units = "cm")


## Figure 2 ====
cutoff_a <- clinicalsignificance:::.calc_cutoff_jt(17, 7, 30, 7, type = "a")[["value"]]
cutoff_b <- clinicalsignificance:::.calc_cutoff_jt(17, 7, 30, 7, type = "b")[["value"]]
cutoff_c <- clinicalsignificance:::.calc_cutoff_jt(17, 7, 30, 7, type = "c")[["value"]]

fig_statistical_cutoffs <- ggplot() +
  geom_function(fun = dnorm, n = 400, args = list(mean = 17, sd = 9), aes(color = "Clinical")) +
  geom_function(fun = dnorm, n = 400, args = list(mean = 30, sd = 7), aes(color = "Functional")) +
  geom_vline(xintercept = cutoff_a, lty = "dotted") +
  geom_vline(xintercept = cutoff_b, lty = "dotted") +
  geom_vline(xintercept = cutoff_c, lty = "dotted") +
  geom_label(aes(label = "a"), x = cutoff_a, y = 0.005) +
  geom_label(aes(label = "b"), x = cutoff_b, y = 0.005) +
  geom_label(aes(label = "c"), x = cutoff_c, y = 0.005) +
  expand_limits(x = c(-10, 60)) +
  labs(x = "Instrument Score", y = "Density", color = "Population") +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank(),
    axis.ticks.y = element_blank()
  )
fig_statistical_cutoffs
ggsave("../Figures//fig-statistical-cutoffs.pdf", plot = fig_statistical_cutoffs, width = 14, height = 8, units = "cm")


## Example visualizations ====
set.seed(20230915)
uniform_data <- tibble(
  "1" = runif(1000, 0, 20),
  "2" = runif(1000, 0, 20)
) |> 
  mutate(id = row_number(), .before = 1) |> 
  pivot_longer("1":"2", names_to = "time") |> 
  mutate(time = as.numeric(time))


# Clinical significance plot
cs_plot <- uniform_data |> 
  cs_anchor(id, time, value, mid_improvement = 2) |> 
  plot(show = category, point_alpha = 0.4) +
  guides(color = guide_legend(nrow = 3))


# Clinical significance plot with cutoff (statistical approach)
cs_plot_cutoff <- uniform_data |> 
  cs_statistical(id, time, value, m_functional = 8, sd_functional = 3, cutoff_type = "c") |> 
  plot(show = category, point_alpha = 0.4) +
  guides(color = guide_legend(nrow = 3))


# Clinical significance plot for combined approach
cs_plot_combined <- uniform_data |> 
  cs_combined(id, time, value, m_functional = 8, sd_functional = 3, cutoff_type = "c", mid_improvement = 3) |> 
  plot(show = category, point_alpha = 0.4) +
  guides(color = guide_legend(nrow = 3))


# # Group level clinical significance plot
# antidepressants |> 
#   cs_anchor(patient, measurement, mom_di, group = condition, target = "group", mid_improvement = 9, effect = "between", post = "After") |> 
#   plot()


# Longitudinal clinical significance plot (RCI method GLM)
cs_plot_longitudinal <- anxiety |> 
  cs_distribution(subject, measurement, anxiety, rci_method = "HLM") |> 
  plot()


cs_plot_layout <- "
11112222
11112222
##3333##
##3333##
44444444
44444444
"

fig_visualizations <- cs_plot + cs_plot_cutoff + cs_plot_combined + cs_plot_longitudinal + plot_layout(guides = "collect", design = cs_plot_layout) + plot_annotation(tag_levels = "A")
ggsave("../Figures//fig-visualizations.pdf", plot = fig_visualizations, width = 14, height = 22, units = "cm")


# Illustrations ----
## Anchor-based approach ####
### Individual ====
library("clinicalsignificance")
claus_2020
set.seed(20230920)

anchor_individual <- cs_anchor(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  mid_improvement = 7
)

anchor_individual
summary(anchor_individual)
plot(anchor_individual)
plot(anchor_individual, show = category)
cs_get_augmented_data(anchor_individual)

# Incorporating the grouping ion the data
anchor_individual_groups <- cs_anchor(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  mid_improvement = 7,
  group = treatment
)

anchor_individual_groups
plot(anchor_individual_groups)


fig_anchor_individual <- plot(anchor_individual)
fig_anchor_individual_groups <- plot(anchor_individual_groups)



### Group ====
# Whole group within
anchor_whole_group <- cs_anchor(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  mid_improvement = 7,
  target = "group"
)

anchor_whole_group
plot(anchor_whole_group)


# Considering the grouping in the data within
anchor_group <- cs_anchor(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  mid_improvement = 7,
  target = "group",
  group = treatment
)

anchor_group


set.seed(20230920)
anchor_group_between <- cs_anchor(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  post = 4,
  mid_improvement = 7,
  target = "group",
  group = treatment,
  effect = "between"
)

anchor_group_between


## Percentage-change approach ####
percentage <- cs_percentage(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  pct_improvement = 0.5
)

percentage
plot(percentage)
fig_percentage <- plot(percentage)




## Distribution-based approach ####
distribution_jt <- cs_distribution(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  reliability = 0.801
)

distribution_jt
plot(distribution_jt)
fig_distribution_jt <- plot(distribution_jt)


distribution_ha <- cs_distribution(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  reliability = 0.801,
  rci_method = "HA"
)

distribution_ha


distribution_hlm <- cs_distribution(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  rci_method = "HLM"
)

distribution_hlm
plot(distribution_hlm)
fig_distribution_hlm <- plot(distribution_hlm)




## Statistical approach ####
statistical_jt <- cs_statistical(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  m_functional = 7.69,
  sd_functional = 7.53,
  cutoff_type = "c"
)

statistical_jt
plot(statistical_jt)
fig_statistical <- plot(statistical_jt)



## Combined approaches ####
combined_cwb <- cs_combined(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  m_functional = 7.69,
  sd_functional = 7.53,
  cutoff_type = "c",
  mid_improvement = 7
)

combined_cwb
plot(combined_cwb)
fig_combined_cwb <- plot(combined_cwb)


combined_jt <- cs_combined(
  data = claus_2020,
  id = id,
  time = time,
  outcome = bdi,
  pre = 1,
  post = 4,
  m_functional = 7.69,
  sd_functional = 7.53,
  cutoff_type = "c",
  reliability = 0.801
)

combined_jt


## Example figures ####
all_layout <- "
123
456
777
"

fig_all_examples <- fig_anchor_individual + 
   fig_anchor_individual_groups + 
   fig_percentage + 
   fig_distribution_jt +  
   fig_statistical + 
   fig_combined_cwb + 
  fig_distribution_hlm +
  plot_layout(design = all_layout, guides = "collect") +
  plot_annotation(tag_levels = "A")

fig_all_examples
ggsave("../Figures//fig-all-examples.pdf", plot = fig_all_examples, width = 14, height = 22, units = "cm")




## Treating chronic pain ####
hechler_2014

hechler_results <- cs_combined(
  data = hechler_2014,
  id = patient,
  time = measurement,
  outcome = disability,
  m_functional = 26.7,
  sd_functional = 9.14,
  cutoff_type = "c",
  reliability = 0.865
)

summary(hechler_results)
fig_hechler <- plot(hechler_results)
ggsave("../Figures//fig-hechler.pdf", plot = fig_hechler, width = 14, height = 8, units = "cm")