# Please use the code below to install the required dependencies
# pkgs <- c("bayesnec", "bayesnec", "drc", "MASS", "corpcor", "dplyr", "purrr",
#           "ggdist", "ggplot2", "tidyr", "stringr", "ggpubr")
## install.packages("remotes")
# remotes::install_cran(pkgs)

library("tidyverse")
library("bayesnec")
options(prompt = 'R> ', continue = '+ ')
knitr::opts_chunk$set(
  purl = TRUE,
  warning = FALSE,
  message = FALSE,
  echo = TRUE,
  cache = TRUE,
  include = TRUE,
  eval = TRUE,
  fig.align = "center",
  fig.pos = "!ht"
)
fix_form <- function(x) {
  gsub("~", "\\textasciitilde{}", x, fixed = TRUE) |>
    gsub("^", "\\^{}", x = _, fixed = TRUE)
}


## response ~ crf(predictor, model = "a_model")


## response | trials(n_trials) ~ crf(log(predictor), model = "a_model")


library("bayesnec")
set.seed(17)
df <- data.frame(x = rgamma(100, 2, 0.1), y = rnorm(100))
form <- bnf(y ~ crf(log(x), model = "nec3param"))
head(model.frame(form, df))


## fit <- bnec(y ~ crf(x, model = "nec3param"), data = nec_data,
##             seed = 17)


timesimplefit <- system.time({
  fit <- bnec(y ~ crf(x, model = "nec3param"),
              data = nec_data,
              seed = 17)
})


summary(fit)


round_digits <- function(x)
  sprintf("%.2f", x)
autoplot(fit, xform = exp) +
  scale_x_continuous(trans = "log", labels = round_digits)


brms_fit <- pull_brmsfit(fit)
plot(brms_fit)


pull_brmsfit(fit) |> prior_summary() |> sample_priors()


check_priors(fit)


library("drc")
library("MASS")
library("corpcor")
library("dplyr")
library("purrr")
library("ggdist")
options(mc.cores = parallel::detectCores())

rounded <- function(value, precision = 1) {
  sprintf(paste0("%.", precision, "f"), round(value, precision))
}

nec3param <- function(beta, nec, top, x) {
  top * exp(-beta * (x - nec) * ifelse(x - nec < 0, 0, 1))
}

my_rmvnorm <- function(mu, Sigma) {
  r <- length(mu)
  L <- try(t(chol(Sigma)))
  if (inherits(L, "try-error")) {
    L <- t(chol(Sigma, pivot = TRUE))
  }
  Z <- rnorm(r)
  (L %*% Z + mu)[, 1]
}

ggdrc_data <- function(x) {
  nec_mean_drc <- summary(x)$coefficients["e:(Intercept)", "Estimate"]
  nec_cis_drc <- confint(x)["e:(Intercept)", ]
  drc_nec <- c(`50%` = nec_mean_drc,
               `2.5%` = nec_cis_drc[["2.5 %"]],
               `97.5%` = nec_cis_drc[["97.5 %"]])
  mod_df <- x$data
  df_x <- data.frame(log_x = seq(min(mod_df$log_x), max(mod_df$log_x), length =
                                 100))
  mu <- coef(x)
  Sigma <- drc:::vcov.drc(x)
  is_pos_def <- corpcor::is.positive.definite(Sigma)
  if (!is_pos_def) {
    Sigma <- corpcor::make.positive.definite(Sigma, tol = 1e-3)
  }
  Y <- matrix(0, 1000, nrow(df_x))
  pars_mv <- MASS::mvrnorm(n = nrow(Y), mu, Sigma)
  for (i in seq_len(nrow(Y))) {
    Y[i, ] <- nec3param(beta = pars_mv[i, 1],
                        nec = pars_mv[i, 3],
                        top = pars_mv[i, 2],
                        df_x$log_x)
  }
  preds <- apply(Y, 2, median_hdci, na.rm = TRUE) |>
    do.call("rbind.data.frame", args = _) |>
    dplyr::select(y, ymin, ymax)
  e_df <- data.frame(
    x_e = c(df_x$log_x, rev(df_x$log_x)),
    y_e = c(preds$y, rep(NA, nrow(preds))),
    y_ci = c(preds$ymin, rev(preds$ymax)),
    x_r = NA,
    y_r = NA,
    nec_vals = NA,
    nec_labs = NA,
    nec_labs_l = NA,
    nec_labs_u = NA
  )
  r_df <- data.frame(
    x_e = NA,
    y_e = NA,
    y_ci = NA,
    x_r = mod_df$log_x,
    y_r = mod_df$fvfm,
    nec_vals = NA,
    nec_labs = NA,
    nec_labs_l = NA,
    nec_labs_u = NA
  )
  n_df <- data.frame(
    x_e = NA,
    y_e = NA,
    y_ci = NA,
    x_r = NA,
    y_r = NA,
    nec_vals = drc_nec,
    nec_labs = c(drc_nec[[1]], NA, NA),
    nec_labs_l = c(drc_nec[[2]], NA, NA),
    nec_labs_u = c(drc_nec[[3]], NA, NA)
  )
  rbind(e_df, r_df, n_df)
}

df_list <- herbicide |>
  dplyr::filter(herbicide %in% c("hexazinone", "tebuthiuron")) |>
  dplyr::mutate(log_x = log(concentration)) |>
  dplyr::arrange(herbicide, concentration) |>
  split(f = ~ herbicide)

bnec_list <- purrr::map(df_list, ~ bnec(
  formula = fvfm ~ crf(log_x, model = "nec3param"),
  data = .x,
  seed = 17
))
drc_list <- purrr::map(df_list, ~ drm(
  formula = fvfm ~ log_x,
  fct = NEC.3(),
  data = .x
))

if (!all(names(bnec_list) == names(drc_list))) {
  stop("Incompatible name order.")
}

bnec_plot_data <- purrr::map_dfr(bnec_list, ggbnec_data, .id = "model") |>
  dplyr::mutate(source = "bayesnec")
drc_plot_data <- purrr::map_dfr(drc_list, ggdrc_data, .id = "model") |>
  dplyr::mutate(source = "drc")
plot_data <- rbind(bnec_plot_data, drc_plot_data) |>
  dplyr::mutate(
    nec_labs = as.character(signif(exp(
      as.numeric(nec_labs)
    ), 3)),
    nec_labs_l = as.character(signif(exp(
      as.numeric(nec_labs_l)
    ), 3)),
    nec_labs_u = as.character(signif(exp(
      as.numeric(nec_labs_u)
    ), 3))
  )


library("ggplot2")
combine_fits <- function(x) {
  ltys <- rep(rep(c(1, 2, 2), length(unique(x$model))), 2)[-c(36:35)]
  lwds <- rep(rep(c(0.5, 0.2, 0.2), length(unique(x$model))), 2)[-c(36:35)]
  ggplot() +
    geom_polygon(
      data = x |> dplyr::filter(!is.na(y_ci)),
      mapping = aes(x = x_e, y = y_ci, fill = source),
      alpha = 0.5,
      show.legend = FALSE
    ) +
    geom_line(
      data = x |> dplyr::filter(!is.na(y_e)),
      mapping = aes(x = x_e, y = y_e, colour = source),
      linetype = 2,
      show.legend = FALSE
    ) +
    geom_point(
      data = x |> dplyr::filter(!is.na(y_r)),
      mapping = aes(x = x_r, y = y_r),
      fill = "grey30",
      shape = 21,
      show.legend = FALSE
    ) +
    geom_vline(
      data = x |> dplyr::filter(!is.na(nec_vals)),
      mapping = aes(xintercept = nec_vals, colour = source),
      linetype = ltys,
      lwd = lwds,
      show.legend = FALSE
    ) +
    geom_text(
      data = x |> dplyr::filter(!is.na(nec_labs), source == "bayesnec"),
      mapping = aes(
        label = paste0("bayesnec: ", nec_labs, " (", nec_labs_l, "-",
                       nec_labs_u, ")")
      ),
      x = Inf,
      y = Inf,
      hjust = 0.95,
      vjust = 1.5,
      size = 2.75,
      colour = "tomato"
    ) +
    geom_text(
      data = x |> dplyr::filter(!is.na(nec_labs), source == "drc"),
      mapping = aes(
        label = paste0("drc: ", nec_labs, " (", nec_labs_l, "-", nec_labs_u,
                       ")")
      ),
      x = Inf,
      y = Inf,
      hjust = 0.95,
      vjust = 4,
      size = 2.75,
      colour = "skyblue2"
    ) +
    coord_cartesian(clip = "off") +
    scale_colour_manual(values = c(bayesnec = "tomato", drc = "skyblue2")) +
    scale_fill_manual(values = c(bayesnec = "tomato", drc = "skyblue2")) +
    ylim(0, 0.8) +
    scale_x_continuous(breaks = log(c(0.08, 0.61, 4.5, 33, 240)),
                       labels = c("0.08", "0.61", "4.5", "33", "240")) +
    theme_classic() +
    facet_wrap( ~ model, scales = "free", ncol = 2) +
    theme(
      strip.text = element_text(hjust = 0),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      strip.background = element_blank(),
      strip.text.x = element_text(size = 12),
      panel.border = element_rect(colour = NA, fill = NA)
    ) +
    labs(x = expression("Concentration " ~ "(" * mu * "g/L)"),
         y = expression("" * Delta * "F" / "F m'"))
}
combine_fits(plot_data)


timeallametrynfits <- system.time({
  ametryn <- herbicide |>
    dplyr::mutate(concentration = log(concentration)) |>
    dplyr::filter(herbicide == "ametryn")
  manecfit_ametryn <- bayesnec::bnec(fvfm ~ crf(concentration, model =
                                                "decline"),
                                     data = ametryn,
                                     seed = 17)
})


## ametryn <- herbicide |>
##   dplyr::mutate(concentration = log(concentration)) |>
##   dplyr::filter(herbicide == "ametryn")
## manecfit_ametryn <- bayesnec::bnec(
##   fvfm ~ crf(concentration, model = "decline"), data = ametryn, seed = 17
## )


check_chains(manecfit_ametryn, filename = "ametryn_check_chains")
check_priors(manecfit_ametryn, filename = "ametryn_check_priors")


## rhat(manecfit_ametryn)


summary(manecfit_ametryn)


## ametryn_plot_all <- autoplot(manecfit_ametryn, all_models = TRUE)


## ametryn_plot <- autoplot(manecfit_ametryn)


ametryn_plot_all <- autoplot(manecfit_ametryn, all_models = TRUE, xform = exp)
ametryn_plot_all +
  ylab(expression("" * Delta * "F" / "F m'")) +
  xlab(expression("Concentration " ~ "(" * mu * "g/L)")) +
  scale_x_continuous(trans = "log", labels = round_digits)


ametryn_plot <- autoplot(manecfit_ametryn, xform = exp)
ametryn_plot +
  ylab(expression("" * Delta * "F" / "F m'")) +
  xlab(expression("Concentration " ~ "(" * mu * "g/L)")) +
  scale_x_continuous(trans = "log", labels = round_digits)


## manecfit <- herbicide |>
##   dplyr::mutate(concentration = log(concentration)) |>
##   split(f = ~ herbicide) |>
##   purrr::map(function(x) {
##     bayesnec::bnec(fvfm ~ crf(concentration, model = "decline"), data = x,
##                    seed = 17)
##   })
## save(manecfit, file = "manecfits.RData")


timeallherbfits <- system.time({
  manecfit <- herbicide |>
    dplyr::mutate(concentration = log(concentration)) |>
    split(f = ~ herbicide) |>
    purrr::map(function(x) {
      bayesnec::bnec(fvfm ~ crf(concentration, model = "decline"),
                     data = x,
                     seed = 17)
    })
})


cleaned_fits <- purrr::map(manecfit, function(x) {
  bad_fits <- rhat(x)$failed
  out_fit <- x
  if (length(bad_fits) > 0) {
    out_fit <- bayesnec::amend(x, drop = bad_fits)
  }
  out_fit
})


library("tidyr")
library("stringr")
modtab <- purrr::map_dfr(cleaned_fits, function(x) {
  x$mod_stats |>
    dplyr::select(model, wi) |>
    dplyr::mutate(wi = round(wi, 3))
}, .id = "herbicide") |>
  tidyr::pivot_wider(names_from = herbicide, values_from = wi) |>
  data.frame()
colnames(modtab) <- stringr::str_to_title(colnames(modtab))


knitr::kable(modtab, caption = "Fitted valid models and their relative
             pseudo-BMA weights for CR curves for the effects of seven
             herbicides on maximum effective quantum yield ($\\Delta F / Fm'$)
                                                            of symbiotic
                                                            dinoflagellates of
                                                            the coral
                                                            \\textit{Seriatopora
                                                            hystrix}. Ame. =
                                                              Ametryn; Atr. =
                                                              Atrazine; Diu. =
                                                              Diuron; Hex. =
                                                              Hexazinone; Irg. =
                                                              Irgarol; Sim. =
                                                              Simazine; Teb. =
                                                              Tebuthiuron.")


library("ggpubr")
all_plots <- lapply(cleaned_fits, function(x) {
  autoplot(x, xform = exp) +
    scale_x_continuous(trans = "log", labels = round_digits) +
    theme(
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      strip.background = element_blank(),
      strip.text.x = element_blank()
    ) +
    ggtitle("")
})
figure <- ggpubr::ggarrange(
  plotlist = all_plots,
  nrow = 4,
  ncol = 2,
  labels = names(all_plots),
  align = "hv",
  font.label = list(
    color = "black",
    size = 12,
    face = "plain"
  )
)


post_comp <- compare_posterior(cleaned_fits, comparison = "n(s)ec")
prob_diff <- post_comp$prob_diff |>
  tidyr::separate(col = comparison, into = c("herbicide", "columns")) |>
  tidyr::pivot_wider(names_from = columns, values_from = prob) |>
  dplyr::mutate(across(where(is.numeric), ~round(.x, 3)))
colnames(prob_diff) <- stringr::str_to_sentence(colnames(prob_diff))
prob_diff$Herbicide <- stringr::str_to_sentence(prob_diff$Herbicide)


ggpubr::annotate_figure(
  figure,
  left = ggpubr::text_grob(expression("" * Delta * "F" / "F m'"), rot = 90),
  bottom = ggpubr::text_grob(expression("Concentration " ~ "(" * mu * "g/L)"))
)


knitr::kable(prob_diff, caption = "Probability of no-difference in no-effect
             toxicity for seven herbicides. Values are based on the proportional
             overlap in predicted posterior probability density of the N(S)EC
             estimates. Ame. = Ametryn; Atr. = Atrazine; Diu. = Diuron; Hex. =
               Hexazinone; Irg. = Irgarol; Sim. = Simazine; Teb. =
               Tebuthiuron.")


cb_pal <- c("#E69F00",
            "#56B4E9",
            "#009E73",
            "#F0E442",
            "#0072B2",
            "#D55E00",
            "#CC79A7")
post_comp$posterior_data |>
  dplyr::mutate(conc = exp(value)) |>
  ggplot(data = _, mapping = aes(x = conc)) +
  geom_density(mapping = aes(
    group = model,
    colour = model,
    fill = model
  ),
  alpha = 0.3) +
  labs(
    x = expression("Concentration " ~ "(" * mu * "g/L)"),
    y = "Density",
    colour = "Model",
    fill = "Model"
  ) +
  scale_fill_manual(values = cb_pal) +
  scale_colour_manual(values = cb_pal) +
  scale_x_continuous(
    trans = "log",
    breaks = c(0.3, 1, 3, 10, 30),
    labels = c(0.3, 1, 3, 10, 30)
  ) +
  theme_classic()


sessionInfo()