# Supplementary File 3. R code for publication forest plot
# Predictors of Higher Patient Satisfaction (Adjusted Odds Ratios)

library(ggplot2)
library(dplyr)
library(tibble)
library(stringr)

df <- tribble(
  ~group, ~predictor, ~comparison, ~aor, ~lower, ~upper, ~p_value,
  "Modifiable service factors", "Staff attitude", "per 1-point increase", 2.76, 2.05, 3.73, "<0.001",
  "Modifiable service factors", "Waiting-time satisfaction", "per 1-point increase", 1.87, 1.39, 2.51, "<0.001",
  "Modifiable service factors", "Medicine availability", "per 1-point increase", 1.43, 1.07, 1.91, "0.016",
  "Respondent and access characteristics", "Female sex", "vs male", 0.71, 0.44, 1.13, "0.143",
  "Respondent and access characteristics", "Visit frequency", "per additional visit", 1.14, 1.03, 1.26, "0.013",
  "Respondent and access characteristics", "Travel-time category", "longer = higher", 0.91, 0.78, 1.06, "0.231",
  "Respondent and access characteristics", "Age", "per year", 0.98, 0.97, 1.00, "0.028",
  "Respondent and access characteristics", "Hospital or other non-CHC facility", "vs CHC", 1.21, 0.84, 1.72, "0.301",
  "District effects (reference = Kambia)", "Kailahun vs Kambia", "", 0.27, 0.17, 0.43, "<0.001",
  "District effects (reference = Kambia)", "Pujehun vs Kambia", "", 0.33, 0.22, 0.50, "<0.001"
) |>
  mutate(
    predictor_wrapped = str_wrap(predictor, width = 30),
    comparison_wrapped = str_wrap(comparison, width = 20),
    or_ci = sprintf("%.2f (%.2f to %.2f)", aor, lower, upper)
  ) |>
  slice(n():1)

group_cols <- c(
  "Modifiable service factors" = "#2C6BAA",
  "Respondent and access characteristics" = "#4C9A8A",
  "District effects (reference = Kambia)" = "#B3563F"
)

df$y <- seq(0, by = 1.35, length.out = nrow(df))
group_pos <- df |>
  group_by(group) |>
  summarise(y = max(y), .groups = "drop")

p <- ggplot(df, aes(x = aor, y = y, colour = group)) +
  geom_hline(aes(yintercept = y), linewidth = 0.3, colour = "#E5E7EB", inherit.aes = FALSE) +
  geom_vline(xintercept = 1, linetype = "dashed", linewidth = 0.4, colour = "#6B7280") +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.15, linewidth = 0.6) +
  geom_point(size = 2.8) +
  geom_text(data = group_pos, aes(x = 0.065, y = y + 0.55, label = group, colour = group),
            inherit.aes = FALSE, hjust = 0, size = 4.1, fontface = "bold") +
  geom_text(aes(x = 0.065, label = predictor_wrapped), hjust = 0, size = 4.0, lineheight = 1.0, colour = "#1F2937") +
  geom_text(aes(x = 0.50, label = comparison_wrapped), hjust = 0, size = 4.0, lineheight = 1.0, colour = "#1F2937") +
  geom_text(aes(x = 8.7, label = or_ci), hjust = 1, size = 4.0, colour = "#1F2937") +
  geom_text(aes(x = 12.1, label = p_value), hjust = 1, size = 4.0, colour = "#1F2937") +
  annotate("text", x = 0.065, y = max(df$y) + 1.95,
           label = "Predictors of Higher Patient Satisfaction (Adjusted Odds Ratios)",
           hjust = 0, size = 4.9, fontface = "bold", colour = "#111827") +
  annotate("text", x = 0.065, y = max(df$y) + 1.10, label = "Predictor", hjust = 0, size = 4.3, fontface = "bold", colour = "#111827") +
  annotate("text", x = 0.50, y = max(df$y) + 1.10, label = "Comparison", hjust = 0, size = 4.3, fontface = "bold", colour = "#111827") +
  annotate("text", x = 8.7, y = max(df$y) + 1.10, label = "Adjusted odds ratio (95% CI)", hjust = 1, size = 4.3, fontface = "bold", colour = "#111827") +
  annotate("text", x = 12.1, y = max(df$y) + 1.10, label = "P value", hjust = 1, size = 4.3, fontface = "bold", colour = "#111827") +
  scale_x_log10(
    limits = c(0.06, 12.5),
    breaks = c(0.25, 0.50, 1.00, 2.00, 4.00),
    labels = c("0.25", "0.50", "1.00", "2.00", "4.00")
  ) +
  scale_y_continuous(limits = c(-0.8, max(df$y) + 2.4), breaks = NULL) +
  scale_colour_manual(values = group_cols) +
  labs(x = "Adjusted odds ratio, log scale", y = NULL) +
  theme_classic(base_size = 13) +
  theme(
    axis.text.x = element_text(size = 13, colour = "#374151"),
    axis.title.x = element_text(size = 14, colour = "#111827"),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.position = "none",
    plot.margin = margin(12, 18, 12, 12)
  )

ggsave("forest_plot_bmc_publication_4k.png", p, width = 16, height = 9, dpi = 300)
ggsave("forest_plot_bmc_publication_vector.pdf", p, width = 16, height = 9, device = cairo_pdf)
ggsave("forest_plot_bmc_publication_vector.svg", p, width = 16, height = 9)