```{r m-prediction-domain-1}
### assess the prediction domain of the original models

## Holocellulose

# orginal model
d_prediction_domain_hol_or <-
  tibble::tibble(
    variable_y = "Holocellulose",
    model = "Original",
    data_type = 
      c(
        rep("Training", sum(d$sample_type != "old magazines")),
        rep("Peat", nrow(d_peat_res)),
        rep("Vegetation", nrow(d_veg_res))
      ),
    predictor_value = 
      c(
        d$carb[d$sample_type != "old magazines"],
        d_peat_res$carb,
        d_veg_res$carb
      )
  )

d_prediction_domain_hol_or <-
  dplyr::bind_cols(
    d_prediction_domain_hol_or,
    m1.2 %>%
      rstanarm::posterior_predict(newdata = d_prediction_domain_hol_or %>%
                                    dplyr::mutate(carb = predictor_value)) %>%
      as.data.frame() %>%
      purrr::map_df(function(x) {
        tibble::tibble(
          fit = median(x),
          lwr = quantile(x, prob = 0.05),
          upr = quantile(x, prob = 0.95)
        )
      })
    )

# binned spectra, best model
d_prediction_domain_hol_bin <-
  dplyr::bind_rows(
    # training data
    d_underfitting[[2]] %>% 
      dplyr::select(-kl, -hol) %>%
      tidyr::pivot_longer(cols = dplyr::everything(),
                          names_to = "variable_x", values_to = "value_x") %>%
      dplyr::mutate(
        wn = 
          d_underfitting_flat[[2]]$x %>%
          rep(nrow(d_underfitting[[2]])),
        sample_id = 
          seq_len(nrow(d_underfitting[[2]])) %>% 
          rep(each = nrow(d_underfitting_flat[[2]])) %>% 
          as.character(),
        data_type = "Training"
      ),
    # peat data
     d_peat_bin20 %>% 
      tidyr::pivot_longer(cols = dplyr::everything(),
                          names_to = "variable_x", values_to = "value_x") %>%
      dplyr::mutate(
        wn = 
          d_peat_bin20_flat$x %>%
          rep(nrow(d_peat_bin20)),
        sample_id = 
          seq_len(nrow(d_peat_bin20)) %>% 
          rep(each = nrow(d_peat_bin20_flat)) %>% 
          magrittr::add(nrow(d_underfitting[[2]])) %>%
          as.character(),
        data_type = "Peat"
      ),
    # vegetation data
     d_veg_bin20 %>% 
      tidyr::pivot_longer(cols = dplyr::everything(),
                          names_to = "variable_x", values_to = "value_x") %>%
      dplyr::mutate(
        wn = 
          d_veg_bin20_flat$x %>%
          rep(nrow(d_veg_bin20)),
        sample_id = 
          seq_len(nrow(d_veg_bin20)) %>% 
          rep(each = nrow(d_veg_bin20_flat)) %>% 
          magrittr::add(nrow(d_underfitting[[2]]) + nrow(d_peat_bin20)) %>%
          as.character(),
        data_type = "Vegetation"
      )
  ) %>%
  dplyr::group_by(data_type, wn) %>%
  dplyr::summarise(value_x_min = min(value_x),
                   value_x_max = max(value_x),
                   .groups = "keep") %>%
  dplyr::mutate(y_variable = "Holocellulose",
                model = "Binned")

## Klason lignin

# original model
d_prediction_domain_kl_or <-
  tibble::tibble(
    variable_y = "Klason lignin",
    model = "Original",
    data_type = 
      c(
        rep("Training", sum(d$sample_type != "office paper")),
        rep("Peat", nrow(d_peat_res)),
        rep("Vegetation", nrow(d_veg_res))
      ),
    predictor_value = 
      c(
        d$arom15arom16[d$sample_type != "office paper"],
        d_peat_res$arom15arom16,
        d_veg_res$arom15arom16
      )
  )

d_prediction_domain_kl_or <-
  dplyr::bind_cols(
    d_prediction_domain_kl_or,
    m2.2 %>%
      rstanarm::posterior_predict(newdata = d_prediction_domain_kl_or %>%
                                    dplyr::mutate(arom15arom16 = predictor_value)) %>%
      as.data.frame() %>%
      purrr::map_df(function(x) {
        tibble::tibble(
          fit = median(x),
          lwr = quantile(x, prob = 0.05),
          upr = quantile(x, prob = 0.95)
        )
      })
    )

# binned spectra, best model
d_prediction_domain_kl_bin <-
  dplyr::bind_rows(
    # training data
    d_underfitting[[2]] %>% 
      dplyr::select(-kl, -hol) %>%
      tidyr::pivot_longer(cols = dplyr::everything(),
                          names_to = "variable_x", values_to = "value_x") %>%
      dplyr::mutate(
        wn = 
          d_underfitting_flat[[2]]$x %>%
          rep(nrow(d_underfitting[[2]])),
        sample_id = 
          seq_len(nrow(d_underfitting[[2]])) %>% 
          rep(each = nrow(d_underfitting_flat[[2]])) %>% 
          as.character(),
        data_type = "Training"
      ),
    # peat data
     d_peat_bin20 %>% 
      tidyr::pivot_longer(cols = dplyr::everything(),
                          names_to = "variable_x", values_to = "value_x") %>%
      dplyr::mutate(
        wn = 
          d_peat_bin20_flat$x %>%
          rep(nrow(d_peat_bin20)),
        sample_id = 
          seq_len(nrow(d_peat_bin20)) %>% 
          rep(each = nrow(d_peat_bin20_flat)) %>% 
          magrittr::add(nrow(d_underfitting[[2]])) %>%
          as.character(),
        data_type = "Peat"
      ),
    # vegetation data
     d_veg_bin20 %>% 
      tidyr::pivot_longer(cols = dplyr::everything(),
                          names_to = "variable_x", values_to = "value_x") %>%
      dplyr::mutate(
        wn = 
          d_veg_bin20_flat$x %>%
          rep(nrow(d_veg_bin20)),
        sample_id = 
          seq_len(nrow(d_veg_bin20)) %>% 
          rep(each = nrow(d_veg_bin20_flat)) %>% 
          magrittr::add(nrow(d_underfitting[[2]]) + nrow(d_peat_bin20)) %>%
          as.character(),
        data_type = "Vegetation"
      )
  ) %>%
  dplyr::group_by(data_type, wn) %>%
  dplyr::summarise(value_x_min = min(value_x),
                   value_x_max = max(value_x),
                   .groups = "keep") %>%
  dplyr::mutate(y_variable = "Klason lignin",
                model = "Binned")

## plots

# colors_prediction_domain <- scales::hue_pal()(3)
# colors_prediction_domain <- scales::grey_pal(0.1, 0.9)(3)[c(3, 2, 1)]
# colors_prediction_domain <- scales::viridis_pal(option = "plasma")(3)
# colors_prediction_domain <- c("grey", "orange1", "lightcoral") %>% rev()
colors_prediction_domain <- scales::brewer_pal(palette = "Paired")(3)
colors_prediction_domain[[3]] <- "salmon"
data_type_levels <- c("Training", "Vegetation", "Peat") %>% rev()

# make colors darker: source: https://gist.github.com/Jfortin1/72ef064469d1703c6b30
colors_prediction_domain_darker <- 
  colors_prediction_domain %>%
  col2rgb() %>%
  magrittr::divide_by(2) %>%
  t() %>%
  grDevices::rgb(maxColorValue = 255)

p1 <-
  d_prediction_domain_hol_or %>%
  dplyr::mutate(
    data_type = factor(data_type, levels = data_type_levels)
  ) %>%
  ggplot(aes(x = predictor_value, fill = data_type)) +
  geom_histogram(colour = NA, bins = 30) +
  labs(x = "carb", y = "Count", title = "Holocellulose") +
  guides(fill = guide_legend(title = "Dataset")) +
  scale_fill_manual(
    breaks = data_type_levels, 
    values = colors_prediction_domain
  )

p2 <-
  d_prediction_domain_kl_or %>%
  dplyr::mutate(
    data_type = factor(data_type, levels = data_type_levels)
  ) %>%
  ggplot(aes(x = predictor_value, fill = data_type)) +
  geom_histogram(colour = NA, bins = 30) +
  labs(x = "arom15arom16", y = "Count", title = "Klason lignin") +
  guides(fill = guide_legend(title = "Dataset")) +
  scale_fill_manual(
    breaks = data_type_levels, 
    values = colors_prediction_domain
  )

p3 <-
  d_prediction_domain_hol_bin %>%
  dplyr::mutate(
    data_type = factor(data_type, levels = data_type_levels)
  ) %>%
  ggplot(aes(x = wn, ymin = value_x_min, ymax = value_x_max, fill = data_type)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_ribbon(colour = NA, alpha = 1) +
  geom_vline(xintercept = 
               m_pars %>%
               dplyr::filter(stringr::str_detect(variable_y, "Holocellulose") & abs(median) >= 0.02) %>%
               dplyr::pull(wn),
             colour = "grey", linetype = 2) + 
  geom_vline(xintercept =  
               m_pars %>%
               dplyr::filter(stringr::str_detect(variable_y, "Klason lignin") & abs(median) >= 0.02) %>%
               dplyr::pull(wn),
             colour = "grey", linetype = 2) + 
  labs(x = expression("Bin wavenumber ["*cm^{-1}*"]"),
       y = "Predictor value",
       title = "Bin width = 20") +
  guides(fill = guide_legend(title = "Dataset")) +
  scale_fill_manual(
    breaks = data_type_levels, 
    values = colors_prediction_domain
  )

p4 <-
  d_prediction_domain_kl_bin %>%
  dplyr::mutate(
    data_type = factor(data_type, levels = data_type_levels)
  ) %>%
  ggplot(aes(x = wn, ymin = value_x_min, ymax = value_x_max, fill = data_type)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_ribbon(colour = NA, alpha = 1) +
  geom_vline(xintercept =  
               m_pars %>%
               dplyr::filter(stringr::str_detect(variable_y, "Klason lignin") & abs(median) >= 0.02) %>%
               dplyr::pull(wn),
             colour = "grey", linetype = 2) + 
  labs(x = expression("Bin wavenumber ["*cm^{-1}*"]"),
       y = "Predictor value",
       title = "Bin width = 10") +
  guides(fill = guide_legend(title = "Dataset")) +
  scale_fill_manual(
    breaks = data_type_levels, 
    values = colors_prediction_domain
  )

# combine plots
p_prediction_domain <-
  (p1 | p2 | p3) + 
  plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") +
  plot_layout(widths = c(1, 1, 3), guides = "collect") &
  theme(legend.position = "bottom",
        legend.box = "vertical")

## old version wit two different plots for holocellulose and Klason lignin
## combine plots
#pp1 <-
#  (p1 | p2) +
#  plot_annotation(title = "Original models")
#
#pp2 <-
#  p3 / p4 +
#  plot_annotation(title = "Models using binned spectra")
#
#p_prediction_domain <- 
#  pp1 / pp2 +
#  plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") +
#  plot_layout(heights = c(0.4, 2), guides = "collect") &
#  theme(legend.position = "bottom",
#        legend.box = "vertical")
```

