# sample code for forest plot form Elise

if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  rio,           # import/export
  tidyverse,
  scales,        # to keep decimal points in numbers with str_glue
  janitor,
  ggplotify,
  finalfit,     # to keep trailing zeros
  forestplot    # forest plots
  )

####VLS
graph <- import("hr_table.csv") |> 
  mutate(model = case_when(model == "full" ~ "Full",
                           model == "lact 1" ~ "Lactation 1",
                           model == "lact>1" ~ "Lactation 2+",
                           model == "not lame" ~ "Not Lame at Enrollment",
                           model == "lame" ~ "Lame at Enrollment"),
         # keeps trailing 0
         risk_big = round_tidy(risk_big,1),
         hr_ci2 = str_glue("{number(hr, accuracy = 0.01)} ({number(lcl, accuracy = 0.01)},{number(ucl, accuracy = 0.01)})"),
  )

vls_forestplot <- graph |>
  forestplot(labeltext = c(model, n, risk_little, risk_big, 
                           hr_ci2),
             mean = hr,
             lower = lcl,
             upper = ucl,
             colgap = unit(3, "mm"),
             boxsize = 0.15,
             align = c("l", "r", "r", "r", "r"),
             clip = c(0.25, 2),
             # for some stupid reason this is in ln scale
             xticks = c(0.25, 1, 2.5),
             xlog = TRUE,
             # standard label
             #xlab = "HR and 95% CI",
             # explanatory label buut hard to do
             # 2 rows is better but doens't save well 
            # xlab = "                  Hazard Ratio          \n
            # Favors Big Model               Favors Little Model"
             
             xlab = "Favors Big Model          Hazard Ratio   Favors Little Model"
             ) |>
  fp_add_header(model = c(" ", "Model Type"),
                n = c(" " ,fp_align_center("N")),
                # fp_align centers the headings
                risk_little = c(fp_align_right("Lame"), 
                                fp_align_center("Little Model")),
                risk_big = c(fp_align_left("Risk %"),
                             fp_align_center("Big Model")),
                hr_ci2 = fp_align_center(c("Hazard Ratio" ,"(95% CI)"))
                )|>
  fp_set_style(box = "#7a0019",
               line = "black",
               txt_gp = fpTxtGp(label = gpar(cex = 0.9), # Adjust cex for size
                                ticks = gpar(cex = 0.9),
                                xlab  = gpar(cex = 0.8)
                                )
               ) |> 
  fp_append_row(model = fp_txt_italic("Effect Modifiers"),
                position = 4, is.summary = FALSE
  )
vls_forestplot
# save graph
# ggsave doesn't work due to error so saved using plot with 800 as width
# both these options aren't great as cut off bottom
png(file = "forestplot_vls.png", width = 600, height = 325)
print(vls_forestplot)
dev.off()

pdf(file = "forestplot_vls.pdf", width = 8, height = 3.25)
print(vls_forestplot)
dev.off()

pdf(file = "forestplot_vls.pdf")
print(vls_forestplot)
dev.off()


####Lesions
# clean up table
graph_les <- import("or_rd_table.csv") |> 
  clean_names() |> 
  mutate(model = case_when(model == "itt" ~ " Intent-to-Treat",
                           model == "ppp" ~ " Per-Protocol",
                           model == "random" ~ " Herd as Random Effect",
                           model == "lact 1" ~ " Lactation 1",
                           model == "lact 2" ~ " Lactation 2+",
                           model == "horn lact 1" ~ 
                             " Hoof Horn Lactation 1",
                           model == "horn lact 2" ~
                             " Hoof Horn Lactation 2+",
                           model == "infectious" ~ " Infectious",
                           model == "horn" ~ " Hoof Horn"
                           ),
         # compared to round round_tidy keeps trailing 0's
         across(contains("risk"), ~ round_tidy(.x * 100,1)
                         ),
         across(contains("rd"), ~ round(.x * 100,1)
                ),
         across(contains("or"), ~ round(.x, digits = 2)
                ),
         # to add columns together 
         # number maintains the 2 decimal places
         # needs to be in this long line so make sure it prints as 1 line
         or_ci = str_glue("{number(or, accuracy = 0.01)} ({number(or_lci, accuracy = 0.01)},{number(or_uci, accuracy = 0.01)})"),
         rd_ci = str_glue("{number(rd, accuracy = 0.1)} ({number(rd_lci, accuracy = 0.1)},{number(rd_uci, accuracy = 0.1)})")
         )

# then create graph with OR as the graph and the rest as text
les_forestplot <- graph_les |>
#  select(model, n, risk_little, risk_big, or, or_lci, or_uci) |> 
  forestplot(labeltext = c(model, n, risk_little, risk_big, 
                           rd_ci, 
                           or_ci),
             mean = or,
             lower = or_lci,
             upper = or_uci,
             colgap = unit(3, "mm"),
             boxsize = 0.15,
             align = c("l", "r", "r", "r", "r", "r"),
             clip = c(0.1, 2),
             # for some stupid reason this is in ln scale
             xticks = c(0.1 ,1, 2.5),
             xlog = TRUE,
             xlab = "          Favors Big Model             Odds Ratio                 Favors Little Model"
             )|> 
  fp_add_header(model = c(" ", "Outcome Model Type"),
                n = c(" " ,fp_align_center("N")),
                # fp_align centers the headings
                risk_little = c(fp_align_right("Lame"), 
                                fp_align_center("Little Model")),
                risk_big = c(fp_align_left("Risk %"),
                             fp_align_center("Big Model")),
                or_ci = fp_align_center(c("Odds Ratio" ,"(95% CI)")),
                rd_ci = fp_align_center(c("Risk Difference", "(95% CI)"))
  )|>
  fp_set_style(box = "#7a0019",
               line = "black",
               txt_gp = fpTxtGp(label = gpar(cex = 0.9), # Adjust cex for size
                                ticks = gpar(cex = 0.9),
                                xlab  = gpar(cex = 0.8)
               ) 
  ) |> 
  fp_append_row(model = fp_txt_italic("Any Lesion"),
                position = 3, is.summary = FALSE
                ) |> 
  fp_append_row(model = fp_txt_italic("Any Lesion Effect Modifiers"),
                position = 7, is.summary = FALSE
                ) |> 
  fp_append_row(model = fp_txt_italic("Lesion"),
                position = 10, is.summary = FALSE
                ) |> 
  fp_append_row(model = fp_txt_italic("Lesion Effect Modifiers"),
                position = 13, is.summary = FALSE
                )

les_forestplot