---
title: "Unmixing model evaluation"
author: "Thomas Chalaux Clergue"
date: "30/05/2023"
output: html_document
---

Setting up working direction

```{r setup}
general.path <- "C:/Users/tchalaux/ownCloud/Projets/Articles/En cours/Tracers_selection_comparison/Data/"

knitr::opts_knit$set(root.dir = general.path) #to modify according to computer
```



# Import data

## Tracer selection dataset

```{r}
library(readxl)

data <- readxl::read_xlsx("Database_Tracer_selection_TChalauxClergue.xlsx", sheet = 2)

#remove lab mixture as they lack OM analysis
data.prop <- as.data.frame(dplyr::filter(data, Class != "Lab Mixture"))

table(data.prop$Class)
```


# First sight

## Studied data

List all the studied properties and their error.

```{r}
# list of columns names with value and associated error

## ORGANIC MATTER
col.OM <- colnames(data.prop)[which(colnames(data.prop) == c("TOC_PrC")):which(colnames(data.prop) == c("TN_PrC"))]
col.OM.uncer <- colnames(data.prop)[which(colnames(data.prop) == c("TOC_SD")):which(colnames(data.prop) == c("TN_SD"))]

## Geochemistry
col.XRF <- colnames(data.prop)[which(colnames(data.prop) == c("Al_mg.kg")):which(colnames(data.prop) == c("Zr_mg.kg"))]
col.XRF.uncer <- colnames(data.prop)[which(colnames(data.prop) == c("Al_RMSE")):which(colnames(data.prop) == c("Zr_RMSE"))]

## VIS colour
col.colour <- colnames(data.prop)[c(which(colnames(data.prop) == c("L_star_D65")):which(colnames(data.prop) == c("h_D65")),
                                    which(colnames(data.prop) == c("A1")):which(colnames(data.prop) == c("Goethite_525nm_FstD")))]
col.colour.uncer <- colnames(data.prop)[c(which(colnames(data.prop) == c("L_star_D65_SD")):which(colnames(data.prop) == c("h_D65_SD")),
                                    which(colnames(data.prop) == c("A1_uncertainty")):which(colnames(data.prop) == c("Goethite_525nm_FstD_uncertainty")))]


prop.values <- c(col.OM, col.XRF, col.colour)
prop.uncertainties <- c(col.OM.uncer, col.XRF.uncer, col.colour.uncer)
```


## Descriptive statistics

```{r}
library(fingR)

fingR::data.watcher(data = data.prop, properties = prop.values, prop.uncer = prop.uncertainties)
```

According to the result we remove properties with some negative values and when uncertainty makes some sample values virtually impossible

```{r}
prop.values <- prop.values[!prop.values %in% c("Co_mg.kg", "Cr_mg.kg", "Cu_mg.kg", "Ni_mg.kg", "Rb_mg.kg", "A3", "Goethite_445nm_FstD")]
prop.uncertainties <- prop.uncertainties[!prop.uncertainties %in% c("Co_RMSE", "Cr_RMSE", "Cu_RMSE", "Ni_RMSE", "Rb_RMSE", "A3_uncertainty", "Goethite_445nm_FstD_uncertainty")]
```



# Tracer selection
## Three-steps method ***

```{r conv approach files}
save.dir.TSM <- "Three_steps_method/"
# create folder if it does not exist
if(!dir.exists(file.path(save.dir.TSM))){
  dir.create(file.path(dirname(save.dir.TSM), "Three_steps_method"))
}
```


### Conservative test

```{r conservative tests}
# Conventional approach
rt.results <- range.tests(data.prop, class = "Class", # variable where difference between source class/levels/groups and mixture is made
                          mixture = "Target", # how to identify mixtures in class
                          criteria = "all", # all range tests: min-max, hinge, whiskers, mean, mean, +/- sd, median,
                          sample.id = "Sample_name",
                          properties = prop.values, # tested properties
                          save.dir = save.dir.TSM, note = "TSM")

rt.results$results.df
```

```{r}
TSM.cons <- fingR::is.conservative(rt.results$results.df, "Property", "RT", 2, "_")
TSM.cons
```


### Discriminant power
```{r}
KW.results <- discriminant.test(data.prop,
                                class = "Class",
                                mixture = "Target",
                                test = "KW", # will perform Kruskal-Wallis (KW)
                                properties = prop.values,
                                # p.level = .05,
                                save.discrim.tests = T, # if save two-samples tests
                                save.dir = save.dir.TSM, note = "TSM")

KW.results
```


```{r}
TSM.discrim <- is.discriminant(KW.results, "Property", "Kruskal.Wallis_p.value", test.pos = 1, sep = "_", p.level = 0.05, note = "TSM")

TSM.discrim
```


### Stepwise selection and LDA
```{r}
# select properties that pass range test and Kruskal-Wallis
tracers.TSM.NSW <- list()
for(test in names(TSM.cons)){
  tracers.TSM.NSW[[test]] <- intersect(TSM.cons[[test]], TSM.discrim$TSM.Kruskal.Wallis)
}

equal <- matrix(nrow = length(tracers.TSM.NSW), ncol = length(tracers.TSM.NSW))
colnames(equal) <- names(tracers.TSM.NSW) ; rownames(equal) <- names(tracers.TSM.NSW)
for(i in 1:length(tracers.TSM.NSW)){
  for(j in 1:length(tracers.TSM.NSW)){
    equal[i,j] <- setequal(tracers.TSM.NSW[[i]], tracers.TSM.NSW[[j]])
  }
}
equal
```
Every tracer selection is unique

```{r}
### Stepwise selection
tracers.TSM.SW <- list()
for(test in names(tracers.TSM.NSW)){
  #do stepwise selection a reorder the selection following prop.values order
  tracers.TSM.SW[[test]] <- intersect(prop.values, fingR::stepwise.selection(data.prop, "Class", tracers.TSM.NSW[[test]], save.dir = save.dir.TSM, note = paste0("TSM_", test))) 
}

equal <- matrix(nrow = length(tracers.TSM.SW), ncol = length(tracers.TSM.SW))
colnames(equal) <- names(tracers.TSM.SW) ; rownames(equal) <- names(tracers.TSM.SW)
for(i in 1:length(tracers.TSM.SW)){
  for(j in 1:length(tracers.TSM.SW)){
    equal[i,j] <- setequal(tracers.TSM.SW[[i]], tracers.TSM.SW[[j]])
  }
}
equal
```
Every tracer selection after stepwise still unique

```{r}
equal <- matrix(nrow = length(tracers.TSM.SW), ncol = length(tracers.TSM.SW))
colnames(equal) <- names(tracers.TSM.SW) ; rownames(equal) <- names(tracers.TSM.SW)
for(i in 1:length(tracers.TSM.SW)){
  for(j in 1:length(tracers.TSM.SW)){
    equal[i,j] <- setequal(tracers.TSM.NSW[[i]], tracers.TSM.SW[[j]]) # is non-stepwise equal stepwise
  }
}
equal
```
mean and median Without and with stepwise selection are equal

```{r}
#we remove mean and median from stepwise list as they are the same in non stepwise version
tracers.TSM.SW <- tracers.TSM.SW[names(tracers.TSM.SW) %in% c("mean", "median") == FALSE]


tracers.TSM <- list(tracers.TSM.NSW, tracers.TSM.SW)
names(tracers.TSM) <- c("NSW", "SW")
tracers.TSM.all <- unlist(tracers.TSM, recursive = F)
#correct mean and median name
names(tracers.TSM.all)[which(names(tracers.TSM.all) %in% c("NSW.mean", "NSW.median"))] <- c("NSW.SW.mean", "NSW.SW.median")

tracers.TSM.all
```


Verify LDA assumptions
```{r eval = F}
library(QuantPsyc)
mult.norm(data.prop)
```


## Consensus Method ***

Ne tourne que sous R vers < 4.2.2

```{r install FingerPro, eval = F}
library(devtools)

install.packages("car", "Ternary", "pairs")

#devtools::install_github("eead-csic-eesa/fingerPro", ref = "master", force = F)

#Need to download the package the github path doesn't work
install.packages("C:/Users/tchalaux/ownCloud/Analyse_resultats_graphs/Eval_unmixing_models/Code/fingerPro-master/fingerPro_1.3.zip", repos = NULL)
```

```{r}
save.dir.CMA <- c("Consensus_method/")
# create folder if it does not exist
if(!dir.exists(file.path(save.dir.CMA))){
  dir.create(file.path(dirname(save.dir.CMA), "Consensus_method"))
}
```


### Formating data
```{r build data}
library(dplyr)

# selecting sample identifier, group class and studied properties
data.CMA <- dplyr::select(data.prop,
                          Sample_name, Class, all_of(c(prop.values)))

# split the main dataset into source and target subsets.
CMA.target <- dplyr::filter(data.CMA, Class == "Target")

CMA.sources <- dplyr::filter(data.CMA, Class != "Target")
CMA.sources <- CMA.sources[order(CMA.sources$Class),] # order source groups

data.CMA2 <- as.data.frame(rbind(CMA.sources, CMA.target))
```


### Conservativeness index (CI)

Note: CI_Method need a data.frame
```{r eval = FALSE}
library(fingerPro)

results_CI <- fingerPro::CI_Method(data.CMA2, points = 2000, Means = FALSE) #need to put means = F if using raw data (which means measured data)

CI.score <- results_CI$CI
CI.score$conservative <- ifelse(CI.score$V1 == 0, TRUE, FALSE)
colnames(CI.score) <- c("variable", "Conservative_index", "conservative")

write.csv(CI.score, paste(save.dir.CMA, "CMA_Conservativeness_index.csv"))

CI.score
```


```{r eval = FALSE}
# do one by one to have readable graphs
CMA.triangle <- "Data/Consensus_method/[graph] Conservativeness_index_plots/"

# create folder if it does not exist
if(!dir.exists(file.path(CMA.triangle))){
  dir.create(file.path(dirname(CMA.triangle), "[graph] Conservativeness_index_plots"))
}

for(i in 1:length(all.data)){
  Ternary_diagrVM(results_CI, tracers = c(i), n_row = 1, n_col = 1) # draw the ternary diagrVM for property i
  dev.print(device = png, filename = paste(CMA.triangle, "CMA", all.data[i], ".png", sep=""), width = 500, height = 500) # save the plot
}
```

### Consensus ranking (CR)
```{r eval = FALSE}
sgeo <- inputSource(data.CMA2) # mean for sources
mgeo <- inputSample(data.CMA2) # rename targets
#not a very versatile code, need to order data on a specific order to be used...

CR_resu <- data.frame("properties" = prop.values)

for(i in 1:nrow(mgeo)){
  crgeo <- cr_3s(source = sgeo, mixture = mgeo[i,], maxiter = 2000, seed = 1997) #CR values for 3 source groups. Long do not run
  
  CR_resu <- left_join(CR_resu, crgeo, by = c("properties" = "tracer"))
}

# overall score
crgeo <- cr_3s(source = sgeo, mixture = mgeo, maxiter = 2000, seed = 1997) #CR values for 3 source groups. Long do not run
CR_resu <- left_join(CR_resu, crgeo, by = c("properties" = "tracer"))
  
colnames(CR_resu) <- c("properties", paste("score", CMA.target$Sample_name, sep="."), "overall")

write.csv(CR_resu, file = paste(save.dir.CMA, "CMA_Consensus_Ranking_Scores.csv"), row.names = F)

CR_resu
```

```{r, eval = F}
tracers.CM <- c("TOC_PrC",	"TN_PrC", "b_star_D65",	"C_star_D65")
```


## List of all tracers selections

Create a list containing all selection of tracers.
Remove redundant lists from CA approach (should be 10)
```{r}
all.tracers.list <- tracers.TSM.all
all.tracers.list[["CM"]] <- tracers.CM

tracers.all <- unique(c(unlist(tracers.TSM), unlist(tracers.CM))) # all selected tracers
```



## Tracer property 
### Correlation
```{r, eval = F}
corr.tracers <- round(cor(dplyr::select(data.prop, all_of(tracers.all)), method = "pearson"), 2)

#corr.tracers

write.csv(corr.tracers, file = paste("Correlation_selected_tracers.csv"), fileEncoding = "latin1")

library(corrplot)
cor.pvalue <- cor.mtest(corr.tracers, method = "pearson")

write.csv(x = cor.pvalue$p, file = "Correlation_signif_selected_tracers.csv")
```

### d13C vs. C/N

```{r}
library(ggplot2) ; library(ggthemes)
library(dplyr)

graph.d13C <- data.prop %>%
  dplyr::filter(Material == "Target") %>%
  dplyr::select(Sample_name, Class, ratio_TOC_TN, d13C_PrM)

graph.d13C$depth_cm <- seq(6.5, 37.5, by=1)

graph.d13C.s <- data.prop %>%
  dplyr::filter(Material != "Target") %>%
  dplyr::select(Sample_name, Class, ratio_TOC_TN, d13C_PrM)

l.width <- .3
a <- .2
p.size <- 1.5
# THE GRAPH
graph.d13C <- ggplot() +
  # Organic matter inputs signatures (Lamb 2006)
  geom_rect(aes(xmin=12,xmax=30, ymin= -32, ymax= -20), fill = "#85cf00", alpha = a, colour = "#85cf00", linewidth = l.width) + # C3 OM SOILS
  geom_rect(aes(xmin=4.5,xmax=10, ymin= -33, ymax= -25), fill = "#2c5b92", alpha = a, colour = "#2c5b92", linewidth = l.width) + # Freshwater POC
  geom_rect(aes(xmin=4,xmax=6, ymin= -27, ymax= -20), fill = "purple", alpha = a, colour = "purple", linewidth = l.width) + # bacteria
  geom_rect(aes(xmin=7,xmax=27, ymin= -32.5, ymax= -26), fill = "white", linetype = "dashed", alpha = 0, colour = "#00501d") + # Freshwater DOC
  geom_rect(aes(xmin=5.5,xmax=9.5, ymin= -30, ymax= -24.5), fill = "white", linetype = "dotted", alpha = 0, colour = "#21538b") + # Freshwater algae
  
  # Sources samples
  geom_point(data = dplyr::filter(graph.d13C.s, Class == "Cropland"), aes(x = ratio_TOC_TN, y = d13C_PrM), size = p.size, shape = 1, color = "#ff6600") +
  geom_point(data = dplyr::filter(graph.d13C.s, Class == "Forest"), aes(x = ratio_TOC_TN, y = d13C_PrM), size = p.size, shape = 2, color = "#85cf00") +
  geom_point(data = dplyr::filter(graph.d13C.s, Class == "Subsoil"), aes(x = ratio_TOC_TN, y = d13C_PrM), size = p.size, shape = 0, color = "#00CCFF") +
  #sediments
  geom_point(data = graph.d13C,   aes(x = ratio_TOC_TN, y = d13C_PrM, colour = -depth_cm), pch = 16, size = 2.3) +
  # axis
  scale_y_continuous("", limits = c(-33, -20), breaks = seq(-34, -20, by = 2)) +
  scale_x_continuous("", limits = c(4, 30), breaks = seq(4, 30, by = 2)) +
  labs(color = "Depth (cm)") +
  
  # THEME
  theme_classic() + theme(axis.title = element_text(size = 9.5),
                          axis.text = element_text(size = 8.5), panel.border = element_blank())

ggsave(filename = "Core_d13C_CN.png", path = paste("Organic_matter/", sep=""), plot = graph.d13C, width = 15, height = 10, dpi = 300, units = c("cm"))
```



# Modelling - MixSIAR

```{r local install MixSIAR eval = F}
library(remotes)
remotes::install_local("C:/Users/tchalaux/ownCloud/Code_R/MixSIAR/MixSIAR-master.zip") # vers. 3.1.11
```

## csv file for MixSIAR
MixSIAR need unique file for each approach. It does not really deal well with single file where variables are not used. Therefore we need to generate a csv file for each approach.
```{r, eval = FALSE}
# WHERE TO SAVE
save.dir.properties <- "MixSIAR_datasets/"
# create folder if it does not exist
if(!dir.exists(file.path(save.dir.properties))){
  dir.create(file.path(dirname(save.dir.properties), "MixSIAR_datasets"))
}

# TARGETS - SEDIMENTS
target.sed <- data.prop %>%
  dplyr::filter(Class == "Target") %>%
  dplyr::select(Sample_name, all_of(tracers.all))

write.csv(target.sed, file = paste(save.dir.properties, "MixSIAR_core_targets.csv", sep=""), row.names = F)

# save targets for each approach
for(app in names(all.tracers.list)){
  write.csv(dplyr::select(target.sed, Sample_name, all_of(all.tracers.list[[app]])),
          file = paste0(save.dir.properties, "MixSIAR_core_targets_", app, ".csv"), row.names = F)
}



# SOURCES #
sources <- data.prop %>%
  dplyr::filter(Class != "Target") %>%
  dplyr::select(Class, all_of(tracers.all)) %>%
  dplyr::group_by(Class) %>%
  dplyr::summarise(across(everything(), list(Mean = mean, SD = sd))) %>% # calculate Mean and standard deviation
  as.data.frame()

# correct colnames
for(i in 1:ncol(sources)){
  if(grepl("Mean", colnames(sources)[i])){
    colnames(sources)[i] <- paste0("Mean",gsub("_Mean", "", x = colnames(sources)[i]))
  }
  if(grepl("SD", colnames(sources)[i])){
    colnames(sources)[i] <- paste0("SD",gsub("_SD", "", x = colnames(sources)[i]))
  }
}
  
sources <- dplyr::left_join(sources,
                            data.prop %>% dplyr::filter(Class != "Target") %>% dplyr::select(Class, all_of(tracers.all)) %>% dplyr::group_by(Class) %>% summarise(n = n()),
                            by = "Class") %>%
  dplyr::select(-Class)
  

rownames(sources) <- c("Cropland", "Forest", "Subsoil")

for(t in tracers.all){
  sources[which(grepl(t, colnames(sources)))] <- round(sources[which(grepl(t, colnames(sources)))], lvl <- lvl.signif(data.prop[[t]]))
}

write.csv(sources, file = paste(save.dir.properties, "MixSIAR_core_sources.csv", sep=""), row.names = TRUE)

# save targets for each approach
for(app in names(all.tracers.list)){
  write.csv(dplyr::select(sources, all_of(c(paste0("Mean", all.tracers.list[[app]]), paste0("SD", all.tracers.list[[app]]))), n),
          file = paste0(save.dir.properties, "MixSIAR_core_sources_", app, ".csv"), row.names = TRUE)
}



# Make discrimination file
# It reflects trophic relations, but there is nothing like that here so its an empty matrix.
s.discr <- matrix(data = 0, nrow = nrow(sources), ncol = (ncol(sources)-1))
colnames(s.discr) <- colnames(sources)[1:(ncol(sources)-1)]
row.names(s.discr) <- row.names(sources)

write.csv(s.discr, file = paste(save.dir.properties, "MixSIAR_core_discrimination.csv", sep=""), row.names = T)
```


## Generate Virtual Mixtures

```{r eval = F}
save.dir.VM <- "MixSIAR_datasets/Virtual_Mixtures/"
# create folder if it does not exist
if(!dir.exists(file.path(save.dir.VM))){
  dir.create(file.path(dirname(save.dir.VM), "Virtual_Mixtures"))
}

# generate contribution
VM.contrib <- VM.contrib.generator(n.sources = 3, groups = c("Cropland", "Forest", "Subsoil"), step = .05, save.dir = save.dir.VM)
# generate property values
VM.dt <- VM.builder(data.prop, material = "Material", source.name = "Source", class = "Class", tracers = tracers.all,
                    contributions = VM.contrib, VM.name = colnames(VM.contrib)[1], save.dir = save.dir.VM)

# save
for(app in names(all.tracers.list)){
  write.csv(dplyr::select(VM.dt$property, mix.names, all_of(all.tracers.list[[app]])),
            file = paste0(save.dir.VM, "VM_properties_", app, ".csv"), row.names = FALSE)
  write.csv(dplyr::select(VM.dt$uncertainty, mix.names, all_of(all.tracers.list[[app]])),
            file = paste0(save.dir.VM, "VM_properties_SD_", app, ".csv"), row.names = FALSE)
}
```

To load previously generated virtual mixtures not to generate new properties distribution
```{r, eval = FALSE}
save.dir.VM <- "MixSIAR_datasets/Virtual_Mixtures/"
# create folder if it does not exist
if(!dir.exists(file.path(save.dir.VM))){
  dir.create(file.path(dirname(save.dir.VM), "Virtual_Mixtures"))
}

VM.prop <- read.csv("MixSIAR_datasets/Virtual_Mixtures/VM_proportions_138.csv")

VM.property <- read.csv("MixSIAR_datasets/Virtual_Mixtures/VM_properties.csv", sep=";")

for(app in names(all.tracers.list)){
  write.csv(dplyr::select(VM.dt$property, Sample_name, all_of(all.tracers.list[[app]])),
            file = paste0(save.dir.VM, "VM_properties_", app, ".csv"), row.names = FALSE)
  write.csv(dplyr::select(VM.dt$uncertainty, Sample_name, all_of(all.tracers.list[[app]])),
            file = paste0(save.dir.VM, "VM_properties_SD_", app, ".csv"), row.names = FALSE)
}
```

##run MixSIAR

```{r MixSIAR output options}
output_options <- list(summary_save = TRUE,
                       summary_name = "summary_statistics",
                       sup_post = FALSE,
                       plot_post_save_pdf = TRUE,
                       plot_post_name = "posterior_density",
                       sup_pairs = FALSE,
                       plot_pairs_save_pdf = TRUE,
                       plot_pairs_name = "pairs_plot",
                       sup_xy = TRUE,
                       plot_xy_save_pdf = FALSE,
                       plot_xy_name = "xy_plot",
                       gelman = TRUE,
                       heidel = FALSE,
                       geweke = TRUE,
                       diag_save = TRUE,
                       diag_name = "diagnostics",
                       indiv_effect = FALSE,
                       plot_post_save_png = FALSE,
                       plot_pairs_save_png = FALSE,
                       plot_xy_save_png = FALSE,
                       diag_save_ggmcmc = FALSE)
```


#### Load data

```{r}
path.mixsiar.VM <- "MixSIAR_datasets/Virtual_Mixtures/"
path.mixsiar.data <- "MixSIAR_datasets/"

MS.mix.VM <- list()
MS.source.VM <- list()
MS.discr.VM <- list()

library(MixSIAR)

for(app in names(all.tracers.list)){
  MS.mix.VM[[app]] <- load_mix_data(filename = paste0(path.mixsiar.VM, paste0("VM_properties_", app,".csv")), # target files
                           iso_names = all.tracers.list[[app]], # optimal tracers (conventionally isotopes)
                           factors = c("mix.names"), # how to differentiate samples
                           fac_random = FALSE, # if the factor is random effect
                           cont_effects = NULL) # continuous effect column
  
  MS.source.VM[[app]] <- load_source_data(filename = paste0(path.mixsiar.data, paste0("MixSIAR_core_sources_", app,".csv")),
                                          source_factors = NULL, conc_dep = FALSE, data_type = "means",  MS.mix.VM[[app]])
  
  MS.discr.VM[[app]] <- load_discr_data(filename = paste0(path.mixsiar.data, "MixSIAR_core_discrimination.csv"),
                               MS.mix.VM[[app]])
}
```

#### Writing JAGS model

```{r, eval = F}
path.model <- "MixSIAR_modelling/"
# create folder if it does not exist
if(!dir.exists(file.path(path.model))){
  dir.create(file.path(dirname(path.model), "MixSIAR_modelling"))
}

for(app in names(all.tracers.list)){
  
  setwd(general.path)
  if(!dir.exists(file.path( paste0(path.model, paste0("MixSIAR_VM_", app, "/")) ))){ #create the folder if it doesn't already exist
  dir.create(file.path(dirname(paste0(path.model, paste0("MixSIAR_VM_", app, "/"))),  paste0("MixSIAR_VM_", app, "/")))
  }
  setwd(paste0(general.path, path.model, paste0("MixSIAR_VM_", app, "/")))
  
  
  write_JAGS_model(paste0("MixSIAR_model_VM_", app, ".txt"),
                   resid_err = FALSE, process_err = TRUE,
                   mix = MS.mix.VM[[app]],
                   source = MS.source.VM[[app]])
  
  setwd(general.path)
}
```

#### Running models

```{r, eval = F}
path.model <- "MixSIAR_modelling/"

jags.VM <- list()

run.type <- "long"
#run.type <- "test"

library(MixSIAR)

for(app in names(all.tracers.list)){
  setwd(paste0(general.path, path.model, paste0("MixSIAR_VM_", app, "/")))
  print(app)
  jags.VM[[app]] <- run_model(run = run.type, mix = MS.mix.VM[[app]],
                              source =  MS.source.VM[[app]],
                              discr = MS.discr.VM[[app]],
                              model_filename = paste0("MixSIAR_model_VM_", app, ".txt"))
}
#  note if "Error: .onload ... 'rgags' -> it's because R version is too old need at least R.2.2
```


#### All statistics

Quite slow...
```{r}
library(fingR)

pred.VM <- list()

for(app in names(all.tracers.list)){
  print(app)
  pred.VM[[app]] <- JAGS.summary(jags.VM[[app]],
                           MS.mix.VM[[app]],
                           MS.source.VM[[app]],
                           path = paste0(general.path, path.model, paste0("MixSIAR_VM_", app, "/")),
                           note = paste0("VM_", app), save_pred = T)
}
```


```{r}
preds.VM <- list()
stats <- c("Median")

for(app in names(all.tracers.list)){
  print(app)
  preds.VM[[app]] <- JAGS.pred(path = paste0(path.model, paste0("MixSIAR_VM_", app, "/"), paste0("summary_statistics_VM_", app, ".csv")),
                               stats, save = T, note = paste0("VM_", app))
}
```

##### Eval model prediction quality

```{r}
obs <- as.data.frame(read.csv("MixSIAR_datasets/Virtual_Mixtures/VM_proportions_138.csv", sep=",")) #if not opened by Excel change sep to ","

error.VM <- list()

for(app in names(all.tracers.list)){
  print(app)
  error.VM[[app]] <- eval.groups(obs, preds.VM[[app]], by = c("mix.names" = "sample"), note = paste0("VM_", app),
                                 path = paste0(path.model, paste0("MixSIAR_VM_", app, "/")))
}

error.VM$NSW.SW.median
```

##### CRPS

Quite slow...
```{r calc CRPS}
crps.VM <- list()

for(app in names(all.tracers.list)){
  print(app)
  crps.VM[[app]] <- JAGS.crps(obs,
                              read.csv(paste0(path.model, paste0("MixSIAR_VM_", app, "/"), paste0("MixSIAR_prevision_VM_", app, ".csv"))),
                              sample.name = "mix.names", source.groups = c("Cropland", "Forest", "Subsoil"),
                              save.dir = paste0(path.model, paste0("MixSIAR_VM_", app, "/")),
                              note = paste0("VM_", app))
}

crps.vals <- data.frame()
for(sources in c("Cropland", "Forest", "Subsoil")){
  val <- c()
  for(app in names(all.tracers.list)){
    val <- append(val, round(mean(crps.VM[[which(names(crps.VM) == app)]][[sources]]), 2))
  }
  crps.vals <- rbind(crps.vals, val)
}
colnames(crps.vals) <- names(all.tracers.list)
rownames(crps.vals) <- c("Cropland", "Forest", "Subsoil")

crps.vals
```


##### Interval width

```{r W50 and W95}
pred.W <- list()

for(app in names(all.tracers.list)){
  print(app)
 pred.W[[app]] <- interval.width(path.to.prev = paste0(path.model, paste0("MixSIAR_VM_", app, "/"), paste0("MixSIAR_prevision_VM_", app, ".csv")),
                                      mean.cal = TRUE, save = T, note = paste0("VM_", app))
}
pred.W$CM$mean
```

#### [Graphs]
##### Statistics

Compile stats
```{r}
stats.VM <- list()
for(app in names(all.tracers.list)){
  stats.VM[[app]] <- as.data.frame(cbind(cbind(cbind(error.VM[[app]][1:2], pred.W[[app]]$mean$W50.mean), error.VM[[app]][3:ncol(error.VM[[app]])]), c(crps.vals[[app]])))
  colnames(stats.VM[[app]]) <- c("Type", "Group", "W50.mean", "ME", "RMSE", "r2", "NSE", "CRPS.mean")
  
  # Correct the unit
  stats.VM[[app]][, c("W50.mean", "ME", "RMSE", "CRPS.mean")] <- stats.VM[[app]][, c("W50.mean", "ME", "RMSE", "CRPS.mean")]*100
  
  # calculate the approach average values
  m <- colMeans(stats.VM[[app]][,c("W50.mean", "ME", "RMSE", "r2", "NSE", "CRPS.mean")])
  #m[c("W50.mean", "ME", "RMSE", "CRPS.mean")] <- round(m[c("W50.mean", "ME", "RMSE", "CRPS.mean")],0)
  #m[c("r2", "NSE")] <- round(m[c("r2", "NSE")],2)
  m <- unname(m)
  stats.VM[[app]][4,1:2] <- c("Median", "Average")
  stats.VM[[app]][4,c("W50.mean", "ME", "RMSE", "r2", "NSE", "CRPS.mean")] <- m
  
  stats.VM[[app]][c("W50.mean", "ME", "RMSE", "CRPS.mean")] <- round(stats.VM[[app]][c("W50.mean", "ME", "RMSE", "CRPS.mean")],0)
  stats.VM[[app]][c("r2", "NSE")] <- round(stats.VM[[app]][c("r2", "NSE")],2)
}
```

```{r}
library(ggplot2) ; library(ggthemes)

path.to.graph <- "MixSIAR_modelling/MixSIAR_VM_graphics/Stats_all/"

all.colours <- c("#ff8288", "#910b1c", "#fbbc70", "#963b12", "#85cf00", "#00501d", "#ffe300", "#ff6600", "#2c5b92", "#00CCFF", "#003e82", "#240691", "#BB00BB")

s <- "W50.mean"
zero.l.width <- 0.2
separator <- 0.4
bar.width <- 5.1
alpha.val <- .2
vline.linewidth <- 0.4
#min.app <- 1 ; max.app <- 8 ; sep.app <- 1
min.app <- 1 ; max.app <- 16 ; sep.app <- 2

for(s in c("W50.mean", "ME", "RMSE", "r2", "NSE", "CRPS.mean")){
  
  local.min <- min(c(sapply(stats.VM, "[[", 1, s), sapply(stats.VM, "[[", 2, s), sapply(stats.VM, "[[", 3, s)))
  if(local.min>0){ local.min <- 0 }
  if(!s %in% c("r2", "NSE")){
    local.max <- max(c(sapply(stats.VM, "[[", 1, s), sapply(stats.VM, "[[", 2, s), sapply(stats.VM, "[[", 3, s)))
  }else{ local.max <- 1}
  
  for(i in 1:4){
      graph.stats <- data.frame("app" = names(all.tracers.list), "value" = unname(sapply(stats.VM, "[[", i, s)))
      graph.stats <- graph.stats[c(1, 8, 2, 9,  4, 11, 3, 10, 5, 6, 12, 7, 13),]
      #graph.stats$pos.x <-  c(rep(1, 2), rep(2, 2), rep(3, 2), rep(4, 2), 5, rep(6,2), 7, 8) + c(rep(c(-separator, separator), 4), 0, c(-separator, separator), 0, 0) #position of bars
      graph.stats$pos.x <- c(rep(1, 2), rep(2, 2), rep(3, 2), rep(4, 2), 5, rep(6,2), 7, 8)*2 + c(rep(c(-separator, separator), 4), 0, c(-separator, separator), 0, 0)
      
      graph.back <- data.frame(y = seq(1,8,1), x = 0)
      
      # THE GRAPH
      plot.stats <- ggplot(data = graph.back) + aes(x = x, y = y) +
        geom_point(size = 0) +
        scale_y_continuous("", breaks = seq(min.app,max.app,sep.app), limits = c(max.app, min.app)) +
        theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),  axis.text.y = element_blank(),  axis.text.x = element_text(size = 12),
                         panel.grid.major.x = element_line(linewidth = 0.2), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(),
                         axis.ticks = element_line(linewidth = 0.3), axis.ticks.y = element_blank(),  plot.background = element_rect(fill='transparent', color=NA))
        
      for(j in rev(1:nrow(graph.stats))){
        plot.stats <- plot.stats +
          geom_segment(y = graph.stats$pos.x[j], yend = graph.stats$pos.x[j],
                       x = 0, xend = graph.stats$value[j],
                       colour = all.colours[j], linewidth = bar.width, alpha = alpha.val)
      }
      
      if(s %in% c("r2", "NSE")){ # ideal value
        plot.stats <- plot.stats +
          scale_x_continuous("", breaks = seq(0, 1,.2), limits = c(local.min, local.max), sec.axis = dup_axis()) +
          geom_vline(xintercept = 0, linetype = "dotted", linewidth = vline.linewidth) +
          geom_vline(xintercept = 1, linetype = "dotted", linewidth = vline.linewidth)
      }else{
        plot.stats <- plot.stats +
          scale_x_continuous("", limits = c(local.min, local.max), sec.axis = dup_axis()) +
          geom_vline(xintercept = 0, linetype = "dotted", linewidth = vline.linewidth)
      }
      
      #plot.stats
      
      #ggsave(filename = paste0(paste(s, c("Cropland", "Forest", "Subsoil")[i], "all", sep = "_"), ".pdf"), plot = plot.stats, path = path.to.graph, width = 6, height = 12, units = c("cm"))
      ggsave(filename = paste0(paste(s, c("Cropland", "Forest", "Subsoil", "Average")[i], "all", sep = "_"), ".png"), plot = plot.stats, path = path.to.graph, width = 6, height = 11, units = c("cm"), dpi = 300)
  }
}
```

#### [HAND MADE Graphs]
##### Obs vs. Pred
###### Clear
```{r}
library(dplyr)
path.model <- "MixSIAR_modelling/"

# Loading observed and predicted data for all approaches
app <- names(all.tracers.list)[1]
OP <- dplyr::select(read.csv(paste(path.model, paste0("MixSIAR_VM_", app), paste0("ObsPred_VM_", app, ".csv"), sep="/")), c("mix.names", "Obs_Cropland", "Obs_Forest", "Obs_Subsoil"))
OP.names <- c("mix.names", "Obs_Cropland", "Obs_Forest", "Obs_Subsoil")
for(app in names(all.tracers.list)){
  OP <- left_join(OP,
                  read.csv(paste(path.model, paste0("MixSIAR_VM_", app), paste0("ObsPred_VM_", app, ".csv"), sep="/")),
                  by = c("mix.names", "Obs_Cropland", "Obs_Forest", "Obs_Subsoil"))
  OP.names <- c(OP.names, paste0("Pred_Cropland_", app), paste0("Pred_Forest_", app), paste0("Pred_Subsoil_", app))
}
colnames(OP) <- OP.names
```

Colour ladder
```{r eval = F}
# _NSW.MM$ + _NSW.MM$
c("#ff8288" + 21 , "#910b1c" + 1,
# _NSW.MMe + _SW.MMe
"#fbbc70" + 21 + 4, "#963b12" + 1 + 4,
# _NSW.hinge + _SW.hinge
"#85cf00" + 15, "#00501d" + 0,
# _NSW.whiskers + _SW.whiskers
"#ffe300" + 23, "#ff6600" + 5,
# _NSW.SW.mean
"#2c5b92" + 24,
# _NSW.mean.sd + _SW.mean.sd
"#00CCFF" + 25 , "#003e82" + 6,
# _NSW.SW.median
"#240691" + 0 + 4,
# _CM
"#BB00BB" + 10)
```


```{r eval = F}
library(ggplot2)
library(ggthemes)

path.to.graph <- "MixSIAR_modelling/MixSIAR_VM_graphics/ObsPred_all/"
pt.size <- 2
lwidth <- .8
acceptance <- 5 #corresponding to 5% error

for(grp in c("Cropland", "Forest", "Subsoil")){
  ###### Minimum-Maximum ######
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MM$", colnames(OP))]]]*100), shape = 21, size = pt.size, colour = "#ff8288", fill = "#ff8288") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.MM$", colnames(OP))]]]*100), shape = 1 , size = pt.size, colour = "#910b1c") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MM$", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#ff8288", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.MM$", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#910b1c", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "MinMax", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "MinMax", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Minimum-Maximum ######
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MMe", colnames(OP))]]]*100), shape = 21, size = pt.size, colour = "#fbbc70", fill = "#fbbc70") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MMe", colnames(OP))]]]*100), shape = 4, size = pt.size, colour = "#fbbc70") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.MMe", colnames(OP))]]]*100), shape = 1 , size = pt.size, colour = "#963b12") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.MMe", colnames(OP))]]]*100), shape = 4 , size = pt.size, colour = "#963b12") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MMe", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#fbbc70", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.MMe", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#963b12", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "MinMaxError", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "MinMaxError", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Hinge ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.hinge", colnames(OP))]]]*100), shape = 15, size = pt.size, colour = "#85cf00", fill = "#85cf00") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.hinge", colnames(OP))]]]*100), shape = 0 , size = pt.size, colour = "#00501d") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.hinge", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#85cf00", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.hinge", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#00501d", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "Hinge", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "Hinge", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
   ###### whiskers ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.whiskers", colnames(OP))]]]*100), shape = 23, size = pt.size, colour = "#ffe300", fill = "#ffe300") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.whiskers", colnames(OP))]]]*100), shape = 5 , size = pt.size, colour = "#ff6600") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.whiskers", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#ffe300", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.whiskers", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#ff6600", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "Whiskers", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "Whiskers", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
   ###### Mean ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_NSW.SW.mean", colnames(OP))]]]*100), shape = 24 , size = pt.size, colour = "#2c5b92") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.SW.mean", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#2c5b92", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "Mean", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "Mean", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Mean +/- SD ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.mean.sd", colnames(OP))]]]*100), shape = 5, size = pt.size, colour = "#00CCFF", fill = "#00CCFF") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.mean.sd", colnames(OP))]]]*100), shape = 6 , size = pt.size, colour = "#003e82") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.mean.sd", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#00CCFF", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.mean.sd", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#003e82", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "MeanSD", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "MeanSD", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)

  
  ###### Median ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_NSW.SW.median", colnames(OP))]]]*100), shape = 0 , size = pt.size, colour = "#240691") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_NSW.SW.median", colnames(OP))]]]*100), shape = 4 , size = pt.size, colour = "#240691") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.SW.median", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#240691", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "Median", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "Median", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  
  ###### Consensus Method ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_CM$", colnames(OP))]]]*100), shape = 10, size = pt.size, colour = "#BB00BB") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_CM", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#BB00BB", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_median", "CM", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_median", "CM", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
}
```

###### Core
Here we need to have the sediment core contributions predictions !
```{r eval = F}
library(ggplot2) ; library(ggthemes) ; library(ggExtra)
path.to.graph <- "MixSIAR_modelling/MixSIAR_VM_graphics/ObsPred_all/"

pt.size <- 2
lwidth <- .8
acceptance <- 5 #corresponding to 5% error
decal <- 2

for(grp in c("Cropland", "Forest", "Subsoil")){
  
  ###### Minimum-Maximum ######
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "NSW.MM")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#ff8288") +
    geom_point(aes(x = 100-decal, y = preds.core[[which(names(preds.core) == "SW.MM")]][[paste0("Median_", grp)]]*100), shape = 4 , color = "#910b1c") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MM$", colnames(OP))]]]*100), shape = 21, size = pt.size, colour = "#ff8288", fill = "#ff8288") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.MM$", colnames(OP))]]]*100), shape = 1 , size = pt.size, colour = "#910b1c") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MM$", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#ff8288", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.MM$", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#910b1c", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "MinMax", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "MinMax", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Minimum-Maximum ######
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "NSW.MMe")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#fbbc70") +
    geom_point(aes(x = 100-decal, y = preds.core[[which(names(preds.core) == "SW.MMe")]][[paste0("Median_", grp)]]*100), shape = 4 , color = "#963b12") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MMe", colnames(OP))]]]*100), shape = 21, size = pt.size, colour = "#fbbc70", fill = "#fbbc70") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MMe", colnames(OP))]]]*100), shape = 4, size = pt.size, colour = "#fbbc70") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.MMe", colnames(OP))]]]*100), shape = 1 , size = pt.size, colour = "#963b12") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.MMe", colnames(OP))]]]*100), shape = 4 , size = pt.size, colour = "#963b12") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.MMe", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#fbbc70", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.MMe", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#963b12", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "MinMaxError", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "MinMaxError", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Hinge ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "NSW.hinge")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#85cf00") +
    geom_point(aes(x = 100-decal, y = preds.core[[which(names(preds.core) == "SW.hinge")]][[paste0("Median_", grp)]]*100), shape = 4 , color = "#00501d") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.hinge", colnames(OP))]]]*100), shape = 15, size = pt.size, colour = "#85cf00", fill = "#85cf00") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.hinge", colnames(OP))]]]*100), shape = 0 , size = pt.size, colour = "#00501d") + 
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.hinge", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#85cf00", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.hinge", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#00501d", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "Hinge", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "Hinge", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
   ###### whiskers ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "NSW.whiskers")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#ffe300") +
    geom_point(aes(x = 100-decal, y = preds.core[[which(names(preds.core) == "SW.whiskers")]][[paste0("Median_", grp)]]*100), shape = 4 , color = "#ff6600") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.whiskers", colnames(OP))]]]*100), shape = 23, size = pt.size, colour = "#ffe300", fill = "#ffe300") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.whiskers", colnames(OP))]]]*100), shape = 5 , size = pt.size, colour = "#ff6600") + 
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.whiskers", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#ffe300", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.whiskers", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#ff6600", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "Whiskers", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "Whiskers", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
   ###### Mean ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "NSW.SW.mean")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#2c5b92") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_NSW.SW.mean$", colnames(OP))]]]*100), shape = 24 , size = pt.size, colour = "#2c5b92") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.SW.mean$", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#2c5b92", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "Mean", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "Mean", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Mean +/- SD ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "NSW.mean.sd")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#00CCFF") +
    geom_point(aes(x = 100-decal, y = preds.core[[which(names(preds.core) == "SW.mean.sd")]][[paste0("Median_", grp)]]*100), shape = 4 , color = "#003e82") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.mean.sd", colnames(OP))]]]*100), shape = 5, size = pt.size, colour = "#00CCFF", fill = "#00CCFF") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_SW.mean.sd", colnames(OP))]]]*100), shape = 6 , size = pt.size, colour = "#003e82") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.mean.sd", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#00CCFF", alpha = 0.2, linewidth = lwidth) +
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_SW.mean.sd", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#003e82", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "MeanSD", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "MeanSD", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)

  
  ###### Median ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "NSW.SW.median")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#240691") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_NSW.SW.median", colnames(OP))]]]*100), shape = 0 , size = pt.size, colour = "#240691") +
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_NSW.SW.median", colnames(OP))]]]*100), shape = 4 , size = pt.size, colour = "#240691") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_NSW.SW.median", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#240691", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "Median", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "Median", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  
  ###### Consensus Method ###### 
  plt.OP <- ggplot() + aes(x = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("Obs", colnames(OP))]]]*100) +
    # 1:1 line
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    # Prediction interval
    geom_point(aes(x = 100, y = preds.core[[which(names(preds.core) == "CM")]][[paste0("Median_", grp)]]*100), shape = 4, color = "#BB00BB") +
    # points
    geom_point(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl( "_CM$", colnames(OP))]]]*100), shape = 10, size = pt.size, colour = "#BB00BB") +
    scale_size_manual(values = c(rep(pt.size, 4))) + # manage points size
    ## regression line
    geom_smooth(aes(y = OP[[colnames(OP)[grepl(grp, colnames(OP)) & grepl("_CM", colnames(OP))]]]*100), method = "lm", formula = y~x, se = FALSE, colour = "#BB00BB", alpha = 0.2, linewidth = lwidth) +
    # scale graph
    scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
    scale_y_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Predicted
    # acceptance
    geom_abline(slope = 1, intercept = acceptance, linetype = "dotted") +
    geom_abline(slope = 1, intercept = -acceptance, linetype = "dotted") +
    # theme
    theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                       axis.text.y = element_text(size = 12), 
                       axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                       panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("obspred_core_median_", "CM", grp, "all", sep = "_"), ".pdf"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("obspred_core_median_", "CM", grp, "all", sep = "_"), ".png"), plot = plt.OP, path = path.to.graph, width = 12, height = 12, units = c("cm"), dpi = 200)
}
```



##### CRPS
```{r}
library(ggthemes)
library(ggplot2)

obs <- as.data.frame(read.csv("MixSIAR_datasets/Virtual_Mixtures/VM_proportions_138.csv", sep=","))
library(dplyr)

# Loading observed and predicted data for all approaches
OCRPS <- obs
colnames(OCRPS)[2:4] <- paste0(colnames(OCRPS), "_Obs")[2:4]
OCRPS.names <- c(colnames(OCRPS))
for(app in names(all.tracers.list)){
  OCRPS <- left_join(OCRPS,
                     crps.VM[[which(names(crps.VM) == app)]],
                     by = c("mix.names" = "sample_name"))
  OCRPS.names <- c(OCRPS.names, paste0("Cropland_", app), paste0("Forest_", app), paste0("Subsoil_", app))
}
colnames(OCRPS) <- OCRPS.names

summary(OCRPS)
```

Max : 0.56500

```{r}
path.to.graph.CRPS <- "MixSIAR_modelling/MixSIAR_VM_graphics/CRPS_all/"

pt.size <- 2
lwidth <- .8
y.max <- 57
y.breaks.steps <- 10

for(grp in c("Cropland", "Forest", "Subsoil")){
  ###### Minimum-Maximum ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.MM$", colnames(OCRPS))]]]*100), shape = 21, size = pt.size, colour = "#ff8288", fill = "#ff8288") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.MM$", colnames(OCRPS))]]]*100), shape = 1, size = pt.size, colour = "#910b1c") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.MM$", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#ff8288", linewidth = lwidth) +
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.MM$", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#910b1c", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "MinMax", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "MinMax", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Minimum-Maximum error ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.MMe", colnames(OCRPS))]]]*100), shape = 21, size = pt.size, colour = "#fbbc70", fill = "#fbbc70") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.MMe", colnames(OCRPS))]]]*100), shape = 4, size = pt.size, colour = "#fbbc70") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.MMe", colnames(OCRPS))]]]*100), shape = 1, size = pt.size, colour = "#963b12") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.MMe", colnames(OCRPS))]]]*100), shape = 4, size = pt.size, colour = "#963b12") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.MMe", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#fbbc70", linewidth = lwidth) +
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.MMe", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#963b12", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "MinMaxError", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "MinMaxError", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Hinge ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.hinge", colnames(OCRPS))]]]*100), shape = 15, size = pt.size, colour = "#85cf00", fill = "#85cf00") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.hinge", colnames(OCRPS))]]]*100), shape = 0, size = pt.size, colour = "#00501d") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.hinge", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#85cf00", linewidth = lwidth) +
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.hinge", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#00501d", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "Hinge", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "Hinge", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Whiskers ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.whiskers", colnames(OCRPS))]]]*100), shape = 23, size = pt.size, colour = "#ffe300", fill = "#ffe300") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.whiskers", colnames(OCRPS))]]]*100), shape = 5, size = pt.size, colour = "#ff6600") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.whiskers", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#ffe300", linewidth = lwidth) +
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.whiskers", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#ff6600", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "Whiskers", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "Whiskers", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Mean ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.SW.mean", colnames(OCRPS))]]]*100), shape = 24, size = pt.size, colour = "#2c5b92") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.SW.mean", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#2c5b92", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "Mean", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "Mean", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Mean SD ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.mean.sd", colnames(OCRPS))]]]*100), shape = 25, size = pt.size, colour = "#00CCFF", fill = "#00CCFF") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.mean.sd", colnames(OCRPS))]]]*100), shape = 6, size = pt.size, colour = "#003e82") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.mean.sd", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#00CCFF", linewidth = lwidth) +
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_SW.mean.sd", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#003e82", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "MeanSD", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "MeanSD", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Median ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.SW.median", colnames(OCRPS))]]]*100), shape = 0, size = pt.size, colour = "#240691") +
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.SW.median", colnames(OCRPS))]]]*100), shape = 4, size = pt.size, colour = "#240691") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_NSW.SW.median", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#240691", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "Median", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "Median", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### CM ######
  plt.CRPS <- ggplot() + aes(x = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("Obs", colnames(OCRPS))]]]*100) +
  # points
  geom_point(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_CM", colnames(OCRPS))]]]*100), shape = 10, size = pt.size, colour = "#BB00BB") +
  # lines
  geom_smooth(aes(y = OCRPS[[colnames(OCRPS)[grepl(grp, colnames(OCRPS)) & grepl("_CM", colnames(OCRPS))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#BB00BB", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("CRPS_median", "CM", grp, "all", sep = "_"), ".pdf"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("CRPS_median", "CM", grp, "all", sep = "_"), ".png"), plot = plt.CRPS, path = path.to.graph.CRPS, width = 12, height = 12, units = c("cm"), dpi = 200)
}
```





##### Prediction inteval width
```{r}
path.to.graph.W50 <- c("MixSIAR_modelling/MixSIAR_VM_graphics")

obs <- as.data.frame(read.csv("MixSIAR_datasets/Virtual_Mixtures/VM_proportions_138.csv", sep=","))
library(dplyr)

# Loading observed and predicted data for all approaches
OW50 <- obs
colnames(OW50)[2:4] <- paste0(colnames(OW50), "_Obs")[2:4]
OW50.names <- c(colnames(OW50))
for(app in names(all.tracers.list)){
  OW50 <- left_join(OW50,
                    tidyr::pivot_wider(dplyr::select(pred.W[[which(names(crps.VM) == app)]]$samples, sample, source, W50), names_from = "source", values_from = c("W50")),
                    by = c("mix.names" = "sample"))
  OW50.names <- c(OW50.names, paste0("Cropland_", app), paste0("Forest_", app), paste0("Subsoil_", app))
}
colnames(OW50) <- OW50.names

summary(OW50)
```
Maximum <- 0.4460

```{r}
library(ggplot2)
path.to.graph.W50 <- "MixSIAR_modelling/MixSIAR_VM_graphics/W50_all/"
pt.size <- 1.5
lwidth <- .8
y.max <- 45
y.breaks.steps <- 10

for(grp in c("Cropland", "Forest", "Subsoil")){
  ###### Minimum-Maximum ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.MM$", colnames(OW50))]]]*100), shape = 21, size = pt.size, colour = "#ff8288", fill = "#ff8288") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.MM$", colnames(OW50))]]]*100), shape = 1, size = pt.size, colour = "#910b1c") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.MM$", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#ff8288", linewidth = lwidth) +
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.MM$", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#910b1c", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "MinMax", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "MinMax", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Minimum-Maximum Error ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.MMe", colnames(OW50))]]]*100), shape = 21, size = pt.size, colour = "#fbbc70", fill = "#fbbc70") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.MMe", colnames(OW50))]]]*100), shape = 4, size = pt.size, colour = "#fbbc70") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.MMe", colnames(OW50))]]]*100), shape = 1, size = pt.size, colour = "#963b12") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.MMe", colnames(OW50))]]]*100), shape = 4, size = pt.size, colour = "#963b12") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.MMe", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#fbbc70", linewidth = lwidth) +
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.MMe", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#963b12", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "MinMaxError", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "MinMaxError", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Hinge ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.hinge", colnames(OW50))]]]*100), shape = 15, size = pt.size, colour = "#85cf00", fill = "#85cf00") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.hinge", colnames(OW50))]]]*100), shape = 0, size = pt.size, colour = "#00501d") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.hinge", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#85cf00", linewidth = lwidth) +
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.hinge", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#00501d", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "Hinge", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "Hinge", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Whiskers ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.whiskers", colnames(OW50))]]]*100), shape = 23, size = pt.size, colour = "#ffe300", fill = "#ffe300") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.whiskers", colnames(OW50))]]]*100), shape = 5, size = pt.size, colour = "#ff6600") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.whiskers", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#ffe300", linewidth = lwidth) +
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.whiskers", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#ff6600", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "Whiskers", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "Whiskers", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Mean ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.SW.mean", colnames(OW50))]]]*100), shape = 2, size = pt.size, colour = "#2c5b92") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.SW.mean", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#2c5b92", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "Mean", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "Mean", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Mean SD ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.mean.sd", colnames(OW50))]]]*100), shape = 25, size = pt.size, colour = "#00CCFF", fill = "#00CCFF") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.mean.sd", colnames(OW50))]]]*100), shape = 6, size = pt.size, colour = "#003e82") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.mean.sd", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#00CCFF", linewidth = lwidth) +
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_SW.mean.sd", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#003e82", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "MeanSD", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "MeanSD", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### Median ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.SW.median", colnames(OW50))]]]*100), shape = 0, size = pt.size, colour = "#240691") +
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.SW.median", colnames(OW50))]]]*100), shape = 4, size = pt.size, colour = "#240691") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_NSW.SW.median", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#240691", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "Median", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "Median", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
  
  ###### CM ######
  plt.W50 <- ggplot() + aes(x = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("Obs", colnames(OW50))]]]*100) +
  # points
  geom_point(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_CM", colnames(OW50))]]]*100), shape = 10, size = pt.size, colour = "#BB00BB") +
  # lines
  geom_smooth(aes(y = OW50[[colnames(OW50)[grepl(grp, colnames(OW50)) & grepl("_CM", colnames(OW50))]]]*100), method = "loess", formula = y~poly(x,2), se = FALSE, colour = "#BB00BB", linewidth = lwidth) +
  #geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  # scale graph
  scale_x_continuous("", limits = c(0, 100), breaks = seq(0, 100, 20), sec.axis = dup_axis()) + #Observed
  scale_y_continuous("", limits = c(0, y.max), breaks = seq(0, 100, y.breaks.steps), sec.axis = dup_axis()) + #Predicted
  # theme
  theme_bw() + theme(legend.position = "none", axis.title = element_text(size = 13),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12), #, angle = 45, hjust = 1
                     panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("W50_median", "CM", grp, "all", sep = "_"), ".pdf"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"))
  ggsave(filename = paste0(paste("W50_median", "CM", grp, "all", sep = "_"), ".png"), plot = plt.W50, path = path.to.graph.W50, width = 12, height = 12, units = c("cm"), dpi = 200)
  
}
```



## Core sediments

### run MixSIAR
```{r MixSIAR output options}
output_options <- list(summary_save = TRUE,
                       summary_name = "summary_statistics",
                       sup_post = FALSE,
                       plot_post_save_pdf = TRUE,
                       plot_post_name = "posterior_density",
                       sup_pairs = FALSE,
                       plot_pairs_save_pdf = TRUE,
                       plot_pairs_name = "pairs_plot",
                       sup_xy = TRUE,
                       plot_xy_save_pdf = FALSE,
                       plot_xy_name = "xy_plot",
                       gelman = TRUE,
                       heidel = FALSE,
                       geweke = TRUE,
                       diag_save = TRUE,
                       diag_name = "diagnostics",
                       indiv_effect = FALSE,
                       plot_post_save_png = FALSE,
                       plot_pairs_save_png = FALSE,
                       plot_xy_save_png = FALSE,
                       diag_save_ggmcmc = FALSE)
```


#### Load data

```{r CA}
path.mixsiar <- "MixSIAR_datasets/"

MS.mix.core <- list()
MS.source.core <- list()
MS.discr.core <- list()

library(MixSIAR)

for(app in names(all.tracers.list)){
  MS.mix.core[[app]] <- load_mix_data(filename = paste0(path.mixsiar, paste0("MixSIAR_core_targets_", app, ".csv")), # target files
                                           iso_names = all.tracers.list[[app]], # optimal tracers (conventionally isotopes)
                                           factors = c("Sample_name"), # how to differentiate samples
                                           fac_random = FALSE, # if the factor is random effect
                                           cont_effects = NULL) # continuous effect column
  
  MS.source.core[[app]] <- load_source_data(filename = paste0(path.mixsiar,paste0("MixSIAR_core_sources_", app,".csv")),
                                             source_factors = NULL, conc_dep = FALSE, data_type = "means",  MS.mix.core[[app]])
  
  MS.discr.core[[app]] <- load_discr_data(filename = paste0(path.mixsiar, "MixSIAR_core_discrimination.csv"),
                               MS.mix.core[[app]])
}
```

#### Writing JAGS model
```{r CAPA, eval = F}
path.model <- "MixSIAR_modelling/"
# create folder if it does not exist
if(!dir.exists(file.path(path.model))){
  dir.create(file.path(dirname(path.model), "MixSIAR_modelling"))
}

for(app in names(all.tracers.list)){
  
  # change working location in order to save at the right place
  setwd(general.path)
  if(!dir.exists(file.path( paste0(path.model, paste0("MixSIAR_core_", app, "/")) ))){ #create the folder if it doesn't already exist
  dir.create(file.path(dirname(paste0(path.model, paste0("MixSIAR_core_", app, "/")) ),  paste0("MixSIAR_core_", app, "/")))
  }
  setwd(paste0(general.path, path.model, paste0("MixSIAR_core_", app, "/")))
  
  write_JAGS_model(paste0("MixSIAR_model_core", app, ".txt"), resid_err = FALSE, process_err = TRUE,
                   mix = MS.mix.core[[app]],
                   source = MS.source.core[[app]])
}
```

#### Running models
```{r CAPA, eval = F}
path.model <- "MixSIAR_modelling/"

jags.core <- list()

run.type <- "long"
#run.type <- "test"

for(app in names(all.tracers.list)){
  setwd(paste0(general.path, path.model, paste0("MixSIAR_core_", app, "/")))
  print(app)
  jags.core[[app]] <- MixSIAR::run_model(run = run.type,
                                         mix = MS.mix.core[[app]],
                                         source =  MS.source.core[[app]],
                                         discr = MS.discr.core[[app]],
                                         model_filename = paste0("MixSIAR_model_core", app, ".txt"))
}
```


#### All statistics
```{r}
pred.core <- list()

for(app in names(all.tracers.list)){
  print(app)
  pred.core[[app]] <- fingR::JAGS.summary(jags.core[[app]], # jag model
                                          mix = MS.mix.core[[app]], # mix
                                          sources = MS.source.core[[app]],
                                          path = paste0(general.path, path.model, paste0("MixSIAR_core_", app, "/")), #where to save
                                          note = paste0("core_", app),
                                          save_pred = T)
}
```

```{r}
preds.core <- list()
stats <- c("Median")

for(app in names(all.tracers.list)){
  print(app)
  preds.core[[app]] <- fingR::JAGS.pred(path = paste0(path.model, paste0("MixSIAR_core_", app, "/"), paste0("summary_statistics_core_", app, ".csv") ),
                               stats, save = T, note = paste0("core_", app))
}
```

##### Interval width

```{r W50 and W95}
pred.core.W <- list()

for(app in names(all.tracers.list)){
  print(app)
 pred.core.W[[app]] <- fingR::interval.width(path.to.prev = paste0(path.model, paste0("MixSIAR_core_", app, "/"), paste0("MixSIAR_prevision_core_", app, ".csv")),
                                      mean.cal = TRUE, save = T, note = paste0("core_", app))
}
```

##### Encompassed sample prediction

```{r}
library(dplyr)

CovPer <- list()
for(app in names(all.tracers.list)){
  app.cp <- list()
  for(grp in colnames(preds.core[[app]])[2:4]){
    app.cp[[grp]] <- round(fingR::ESP(preds.core[[app]][[grp]], preds.VM[[app]][[grp]]), 0)
  }
  CovPer[[app]] <- app.cp
}
data.frame(bind_rows(CovPer), row.names = names(all.tracers.list))
```

##### Min, max, mean, SD
```{r}
general.stats <- list()
for(app in names(all.tracers.list)){
  app.stats <- list()
  app.stats[["min"]] <- c(round(min(round(preds.core[[app]][[2]]*100,0)),0),round(min(round(preds.core[[app]][[3]]*100,0)),0), round(min(round(preds.core[[app]][[4]]*100,0)),0))
  app.stats[["max"]] <- c(round(max(round(preds.core[[app]][[2]]*100,0)),0),round(max(round(preds.core[[app]][[3]]*100,0)),0), round(max(round(preds.core[[app]][[4]]*100,0)),0))
  app.stats[["mean"]] <- c(round(mean(round(preds.core[[app]][[2]]*100,0)),0),round(mean(round(preds.core[[app]][[3]]*100,0)),0), round(mean(round(preds.core[[app]][[4]]*100,0)),0))
  app.stats[["sd"]] <- c(round(sd(round(preds.core[[app]][[2]]*100,0)),0),round(sd(round(preds.core[[app]][[3]]*100,0)),0), round(sd(round(preds.core[[app]][[4]]*100,0)),0))
  general.stats[[app]] <- app.stats
}
general.stats
```

### [Graphs]

#### Predictions

```{r adding depth}
graph.core <- preds.core

## adding RMSE value as prediction error
for(app in names(all.tracers.list)){
  for(grp in c("Cropland", "Forest", "Subsoil")){
    graph.core[[app]] <- cbind(graph.core[[app]], error.VM[[app]][grepl("Median", error.VM[[app]]$Type) & grepl(grp, error.VM[[app]]$Group), "RMSE"])
  }
  colnames(graph.core[[app]])[(ncol(graph.core[[app]])-2):ncol(graph.core[[app]])] <- paste0("RMSE_", c("Cropland", "Forest", "Subsoil"))
  graph.core[[app]]$depth <- seq(6.5, 37.5, 1)
}
```


```{r}
for(app in names(all.tracers.list)){
  if(!dir.exists(file.path( paste0(path.model, paste0("MixSIAR_core_", app, "/"), "graphs/") ))){ #create the folder if it doesn't already exist
  dir.create(file.path(dirname(paste0(path.model, paste0("MixSIAR_core_", app, "/"), "graphs/")),  "graphs"))
  }
}
```

##### Vertical contributions

```{r, eval = F}
path.to.graph.core <- "MixSIAR_modelling/MixSIAR_core_graphics/"
marge <- 19
zero.l.width <- 0.55

for(grp in c("Cropland", "Forest", "Subsoil")){
  
  grp.val <- paste0("Median_", grp)
  grp.RMSE <- paste0("RMSE_", grp)
  
  ###### Minimum-maximum ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["NSW.MM"]]$depth, ymin = round(graph.core[["NSW.MM"]][[grp.val]]*100)-round(graph.core[["NSW.MM"]][[grp.RMSE]]*100), ymax = round(graph.core[["NSW.MM"]][[grp.val]]*100)+round(graph.core[["NSW.MM"]][[grp.RMSE]]*100) ), fill = "#ff8288", alpha = 0.1) +
    geom_ribbon(aes(x = graph.core[["SW.MM"]]$depth, ymin = round(graph.core[["SW.MM"]][[grp.val]]*100)-round(graph.core[["SW.MM"]][[grp.RMSE]]*100), ymax = round(graph.core[["SW.MM"]][[grp.val]]*100)+round(graph.core[["SW.MM"]][[grp.RMSE]]*100) ), fill = "#910b1c", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["NSW.MM"]]$depth, y = round(graph.core[["NSW.MM"]][[grp.val]]*100)), colour = "#ff8288", fill = "#ff8288", shape = 21, size = 2) +
    geom_line(aes(x = graph.core[["NSW.MM"]]$depth, y = round(graph.core[["NSW.MM"]][[grp.val]]*100)),  colour = "#ff8288", linewidth = .7) +
    geom_point(aes(x = graph.core[["SW.MM"]]$depth, y = round(graph.core[["SW.MM"]][[grp.val]]*100)), colour = "#910b1c", shape = 1, size = 2) +
    geom_line(aes(x = graph.core[["SW.MM"]]$depth, y = round(graph.core[["SW.MM"]][[grp.val]]*100)),  colour = "#910b1c", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["NSW.MM"]]$depth, y = round(graph.core[["NSW.MM"]][[grp.val]]*100)-round(graph.core[["NSW.MM"]][[grp.RMSE]]*100)),  colour = "#ff8288", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["NSW.MM"]]$depth, y = round(graph.core[["NSW.MM"]][[grp.val]]*100)+round(graph.core[["NSW.MM"]][[grp.RMSE]]*100)),  colour = "#ff8288", linewidth = .4, alpha = .25) +
    geom_line(aes(x = graph.core[["SW.MM"]]$depth, y = round(graph.core[["SW.MM"]][[grp.val]]*100)-round(graph.core[["SW.MM"]][[grp.RMSE]]*100)),  colour = "#910b1c", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["SW.MM"]]$depth, y = round(graph.core[["SW.MM"]][[grp.val]]*100)+round(graph.core[["SW.MM"]][[grp.RMSE]]*100)),  colour = "#910b1c", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "MinMax", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "MinMax", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
  
  ###### Minimum-maximum Error ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["NSW.MMe"]]$depth, ymin = round(graph.core[["NSW.MMe"]][[grp.val]]*100)-round(graph.core[["NSW.MMe"]][[grp.RMSE]]*100), ymax = round(graph.core[["NSW.MMe"]][[grp.val]]*100)+round(graph.core[["NSW.MMe"]][[grp.RMSE]]*100) ), fill = "#fbbc70", alpha = 0.1) +
    geom_ribbon(aes(x = graph.core[["SW.MMe"]]$depth, ymin = round(graph.core[["SW.MMe"]][[grp.val]]*100)-round(graph.core[["SW.MMe"]][[grp.RMSE]]*100), ymax = round(graph.core[["SW.MMe"]][[grp.val]]*100)+round(graph.core[["SW.MMe"]][[grp.RMSE]]*100) ), fill = "#963b12", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["NSW.MMe"]]$depth, y = round(graph.core[["NSW.MMe"]][[grp.val]]*100)), colour = "#fbbc70", fill = "#fbbc70", shape = 21, size = 2) +
    geom_point(aes(x = graph.core[["NSW.MMe"]]$depth, y = round(graph.core[["NSW.MMe"]][[grp.val]]*100)), colour = "#fbbc70", shape = 4, size = 2) +
    geom_line(aes(x = graph.core[["NSW.MMe"]]$depth, y = round(graph.core[["NSW.MMe"]][[grp.val]]*100)),  colour = "#fbbc70", linewidth = .7) +
    geom_point(aes(x = graph.core[["SW.MMe"]]$depth, y = round(graph.core[["SW.MMe"]][[grp.val]]*100)), colour = "#963b12", shape = 1, size = 2) +
    geom_point(aes(x = graph.core[["SW.MMe"]]$depth, y = round(graph.core[["SW.MMe"]][[grp.val]]*100)), colour = "#963b12", shape = 4, size = 2) +
    geom_line(aes(x = graph.core[["SW.MMe"]]$depth, y = round(graph.core[["SW.MMe"]][[grp.val]]*100)),  colour = "#963b12", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["NSW.MMe"]]$depth, y = round(graph.core[["NSW.MMe"]][[grp.val]]*100)-round(graph.core[["NSW.MMe"]][[grp.RMSE]]*100)),  colour = "#fbbc70", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["NSW.MMe"]]$depth, y = round(graph.core[["NSW.MMe"]][[grp.val]]*100)+round(graph.core[["NSW.MMe"]][[grp.RMSE]]*100)),  colour = "#fbbc70", linewidth = .4, alpha = .25) +
    geom_line(aes(x = graph.core[["SW.MMe"]]$depth, y = round(graph.core[["SW.MMe"]][[grp.val]]*100)-round(graph.core[["SW.MMe"]][[grp.RMSE]]*100)),  colour = "#963b12", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["SW.MMe"]]$depth, y = round(graph.core[["SW.MMe"]][[grp.val]]*100)+round(graph.core[["SW.MMe"]][[grp.RMSE]]*100)),  colour = "#963b12", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "MinMaxError", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "MinMaxError", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
  
  ###### Hinge ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["NSW.hinge"]]$depth, ymin = round(graph.core[["NSW.hinge"]][[grp.val]]*100)-round(graph.core[["NSW.hinge"]][[grp.RMSE]]*100), ymax = round(graph.core[["NSW.hinge"]][[grp.val]]*100)+round(graph.core[["NSW.hinge"]][[grp.RMSE]]*100) ), fill = "#85cf00", alpha = 0.1) +
    geom_ribbon(aes(x = graph.core[["SW.hinge"]]$depth, ymin = round(graph.core[["SW.hinge"]][[grp.val]]*100)-round(graph.core[["SW.hinge"]][[grp.RMSE]]*100), ymax = round(graph.core[["SW.hinge"]][[grp.val]]*100)+round(graph.core[["SW.hinge"]][[grp.RMSE]]*100) ), fill = "#00501d", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["NSW.hinge"]]$depth, y = round(graph.core[["NSW.hinge"]][[grp.val]]*100)), colour = "#85cf00", fill = "#85cf00", shape = 15, size = 2) +
    geom_line(aes(x = graph.core[["NSW.hinge"]]$depth, y = round(graph.core[["NSW.hinge"]][[grp.val]]*100)),  colour = "#85cf00", linewidth = .7) +
    geom_point(aes(x = graph.core[["SW.hinge"]]$depth, y = round(graph.core[["SW.hinge"]][[grp.val]]*100)), colour = "#00501d", shape = 0, size = 2) +
    geom_line(aes(x = graph.core[["SW.hinge"]]$depth, y = round(graph.core[["SW.hinge"]][[grp.val]]*100)),  colour = "#00501d", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["NSW.hinge"]]$depth, y = round(graph.core[["NSW.hinge"]][[grp.val]]*100)-round(graph.core[["NSW.hinge"]][[grp.RMSE]]*100)),  colour = "#85cf00", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["NSW.hinge"]]$depth, y = round(graph.core[["NSW.hinge"]][[grp.val]]*100)+round(graph.core[["NSW.hinge"]][[grp.RMSE]]*100)),  colour = "#85cf00", linewidth = .4, alpha = .25) +
    geom_line(aes(x = graph.core[["SW.hinge"]]$depth, y = round(graph.core[["SW.hinge"]][[grp.val]]*100)-round(graph.core[["SW.hinge"]][[grp.RMSE]]*100)),  colour = "#00501d", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["SW.hinge"]]$depth, y = round(graph.core[["SW.hinge"]][[grp.val]]*100)+round(graph.core[["SW.hinge"]][[grp.RMSE]]*100)),  colour = "#00501d", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "Hinge", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "Hinge", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
  
  
  ###### Whiskers ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["NSW.whiskers"]]$depth, ymin = round(graph.core[["NSW.whiskers"]][[grp.val]]*100)-round(graph.core[["NSW.whiskers"]][[grp.RMSE]]*100), ymax = round(graph.core[["NSW.whiskers"]][[grp.val]]*100)+round(graph.core[["NSW.whiskers"]][[grp.RMSE]]*100) ), fill = "#ffe300", alpha = 0.1) +
    geom_ribbon(aes(x = graph.core[["SW.whiskers"]]$depth, ymin = round(graph.core[["SW.whiskers"]][[grp.val]]*100)-round(graph.core[["SW.whiskers"]][[grp.RMSE]]*100), ymax = round(graph.core[["SW.whiskers"]][[grp.val]]*100)+round(graph.core[["SW.whiskers"]][[grp.RMSE]]*100) ), fill = "#ff6600", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["NSW.whiskers"]]$depth, y = round(graph.core[["NSW.whiskers"]][[grp.val]]*100)), colour = "#ffe300", fill = "#ffe300", shape = 23, size = 2) +
    geom_line(aes(x = graph.core[["NSW.whiskers"]]$depth, y = round(graph.core[["NSW.whiskers"]][[grp.val]]*100)),  colour = "#ffe300", linewidth = .7) +
    geom_point(aes(x = graph.core[["SW.whiskers"]]$depth, y = round(graph.core[["SW.whiskers"]][[grp.val]]*100)), colour = "#ff6600", shape = 5, size = 2) +
    geom_line(aes(x = graph.core[["SW.whiskers"]]$depth, y = round(graph.core[["SW.whiskers"]][[grp.val]]*100)),  colour = "#ff6600", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["NSW.whiskers"]]$depth, y = round(graph.core[["NSW.whiskers"]][[grp.val]]*100)-round(graph.core[["NSW.whiskers"]][[grp.RMSE]]*100)),  colour = "#ffe300", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["NSW.whiskers"]]$depth, y = round(graph.core[["NSW.whiskers"]][[grp.val]]*100)+round(graph.core[["NSW.whiskers"]][[grp.RMSE]]*100)),  colour = "#ffe300", linewidth = .4, alpha = .25) +
    geom_line(aes(x = graph.core[["SW.whiskers"]]$depth, y = round(graph.core[["SW.whiskers"]][[grp.val]]*100)-round(graph.core[["SW.whiskers"]][[grp.RMSE]]*100)),  colour = "#ff6600", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["SW.whiskers"]]$depth, y = round(graph.core[["SW.whiskers"]][[grp.val]]*100)+round(graph.core[["SW.whiskers"]][[grp.RMSE]]*100)),  colour = "#ff6600", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "Whiskers", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "Whiskers", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
  
  
  ###### Mean ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["NSW.SW.mean"]]$depth, ymin = round(graph.core[["NSW.SW.mean"]][[grp.val]]*100)-round(graph.core[["NSW.SW.mean"]][[grp.RMSE]]*100), ymax = round(graph.core[["NSW.SW.mean"]][[grp.val]]*100)+round(graph.core[["NSW.SW.mean"]][[grp.RMSE]]*100) ), fill = "#2c5b92", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["NSW.SW.mean"]]$depth, y = round(graph.core[["NSW.SW.mean"]][[grp.val]]*100)), colour = "#2c5b92", shape = 24, size = 2) +
    geom_line(aes(x = graph.core[["NSW.SW.mean"]]$depth, y = round(graph.core[["NSW.SW.mean"]][[grp.val]]*100)),  colour = "#2c5b92", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["NSW.SW.mean"]]$depth, y = round(graph.core[["NSW.SW.mean"]][[grp.val]]*100)-round(graph.core[["NSW.SW.mean"]][[grp.RMSE]]*100)),  colour = "#2c5b92", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["NSW.SW.mean"]]$depth, y = round(graph.core[["NSW.SW.mean"]][[grp.val]]*100)+round(graph.core[["NSW.SW.mean"]][[grp.RMSE]]*100)),  colour = "#2c5b92", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "Mean", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "Mean", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
  
  
  ###### Mean SD ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["NSW.mean.sd"]]$depth, ymin = round(graph.core[["NSW.mean.sd"]][[grp.val]]*100)-round(graph.core[["NSW.mean.sd"]][[grp.RMSE]]*100), ymax = round(graph.core[["NSW.mean.sd"]][[grp.val]]*100)+round(graph.core[["NSW.mean.sd"]][[grp.RMSE]]*100) ), fill = "#00CCFF", alpha = 0.1) +
    geom_ribbon(aes(x = graph.core[["SW.mean.sd"]]$depth, ymin = round(graph.core[["SW.mean.sd"]][[grp.val]]*100)-round(graph.core[["SW.mean.sd"]][[grp.RMSE]]*100), ymax = round(graph.core[["SW.mean.sd"]][[grp.val]]*100)+round(graph.core[["SW.mean.sd"]][[grp.RMSE]]*100) ), fill = "#003e82", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["NSW.mean.sd"]]$depth, y = round(graph.core[["NSW.mean.sd"]][[grp.val]]*100)), colour = "#00CCFF", fill = "#00CCFF", shape = 25, size = 2) +
    geom_line(aes(x = graph.core[["NSW.mean.sd"]]$depth, y = round(graph.core[["NSW.mean.sd"]][[grp.val]]*100)),  colour = "#00CCFF", linewidth = .7) +
    geom_point(aes(x = graph.core[["SW.mean.sd"]]$depth, y = round(graph.core[["SW.mean.sd"]][[grp.val]]*100)), colour = "#003e82", shape = 6, size = 2) +
    geom_line(aes(x = graph.core[["SW.mean.sd"]]$depth, y = round(graph.core[["SW.mean.sd"]][[grp.val]]*100)),  colour = "#003e82", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["NSW.mean.sd"]]$depth, y = round(graph.core[["NSW.mean.sd"]][[grp.val]]*100)-round(graph.core[["NSW.mean.sd"]][[grp.RMSE]]*100)),  colour = "#00CCFF", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["NSW.mean.sd"]]$depth, y = round(graph.core[["NSW.mean.sd"]][[grp.val]]*100)+round(graph.core[["NSW.mean.sd"]][[grp.RMSE]]*100)),  colour = "#00CCFF", linewidth = .4, alpha = .25) +
    geom_line(aes(x = graph.core[["SW.mean.sd"]]$depth, y = round(graph.core[["SW.mean.sd"]][[grp.val]]*100)-round(graph.core[["SW.mean.sd"]][[grp.RMSE]]*100)),  colour = "#003e82", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["SW.mean.sd"]]$depth, y = round(graph.core[["SW.mean.sd"]][[grp.val]]*100)+round(graph.core[["SW.mean.sd"]][[grp.RMSE]]*100)),  colour = "#003e82", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "MeanSD", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "MeanSD", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
  
  ###### Median ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["NSW.SW.median"]]$depth, ymin = round(graph.core[["NSW.SW.median"]][[grp.val]]*100)-round(graph.core[["NSW.SW.median"]][[grp.RMSE]]*100), ymax = round(graph.core[["NSW.SW.median"]][[grp.val]]*100)+round(graph.core[["NSW.SW.median"]][[grp.RMSE]]*100) ), fill = "#240691", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["NSW.SW.median"]]$depth, y = round(graph.core[["NSW.SW.median"]][[grp.val]]*100)), colour = "#240691", shape = 0, size = 2) +
    geom_point(aes(x = graph.core[["NSW.SW.median"]]$depth, y = round(graph.core[["NSW.SW.median"]][[grp.val]]*100)), colour = "#240691", shape = 4, size = 2) +
    geom_line(aes(x = graph.core[["NSW.SW.median"]]$depth, y = round(graph.core[["NSW.SW.median"]][[grp.val]]*100)),  colour = "#240691", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["NSW.SW.median"]]$depth, y = round(graph.core[["NSW.SW.median"]][[grp.val]]*100)-round(graph.core[["NSW.SW.median"]][[grp.RMSE]]*100)),  colour = "#240691", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["NSW.SW.median"]]$depth, y = round(graph.core[["NSW.SW.median"]][[grp.val]]*100)+round(graph.core[["NSW.SW.median"]][[grp.RMSE]]*100)),  colour = "#240691", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "Median", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "Median", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
  
  ###### CM ###### 
  plot.core <- ggplot() +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    geom_hline(yintercept = 100, linetype = "dashed", linewidth = zero.l.width, color = "#b5b8b7") +
    # error
    geom_ribbon(aes(x = graph.core[["CM"]]$depth, ymin = round(graph.core[["CM"]][[grp.val]]*100)-round(graph.core[["CM"]][[grp.RMSE]]*100), ymax = round(graph.core[["CM"]][[grp.val]]*100)+round(graph.core[["CM"]][[grp.RMSE]]*100) ), fill = "#BB00BB", alpha = 0.1) +
    # points
    geom_point(aes(x = graph.core[["CM"]]$depth, y = round(graph.core[["CM"]][[grp.val]]*100)), colour = "#BB00BB", shape = 10, size = 2) +
    geom_line(aes(x = graph.core[["CM"]]$depth, y = round(graph.core[["CM"]][[grp.val]]*100)),  colour = "#BB00BB", linewidth = .7) +
    # error line
    geom_line(aes(x = graph.core[["CM"]]$depth, y = round(graph.core[["CM"]][[grp.val]]*100)-round(graph.core[["CM"]][[grp.RMSE]]*100)),  colour = "#BB00BB", linewidth = .4, alpha = .25) + 
    geom_line(aes(x = graph.core[["CM"]]$depth, y = round(graph.core[["CM"]][[grp.val]]*100)+round(graph.core[["CM"]][[grp.RMSE]]*100)),  colour = "#BB00BB", linewidth = .4, alpha = .25) +
    # scale graph
    scale_x_reverse("", breaks = seq(7, 38, 2), sec.axis = dup_axis()) +
    scale_y_continuous("", breaks = seq(0, 100, 20), limits = c(0-marge, 100+marge), sec.axis = dup_axis()) +
    coord_flip() +
    # theme
    theme_clean() +
    theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
          panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))
  
  #ggsave(filename = paste0(paste("MixSIAR_core", "CM", grp, sep = "_"), ".pdf"), plot = plot1.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"))
  ggsave(filename = paste0(paste("MixSIAR_core", "CM", grp, sep = "_"), ".png"), plot = plot.core, path = path.to.graph.core, height = 17, width = 9, units = c("cm"), dpi = 200)
}
```


#### Tracers graphics

```{r}
library(ggplot2) ; library(ggthemes) ; library(dplyr)

dir.core <- "Core_tracer_distribution"
if(!dir.exists(file.path( dir.core ))){ #create the folder if it doesn't already exist
  dir.create(file.path(dirname(dir.core),  dir.core))
}

dt.graph <- dplyr::filter(data.prop, Class == "Target")
dt.graph$depth_cm <- seq(6.5, 37.5, 1)

dt.sources  <- dplyr::filter(data.prop, Class != "Target")

dt.srcs <- dt.sources %>%
  dplyr::select(Class, all_of(tracers.all)) %>%
  group_by(Class) %>%
  summarise(across(everything(), list(Mean=mean, sd=sd)))

for(prop in tracers.all){
    dt.srcs[[paste0(prop, "_Mean")]] <- round(dt.srcs[[paste0(prop, "_Mean")]] , fingR::lvl.signif(dt.sources[[prop]]))
    dt.srcs[[paste0(prop, "_sd")]] <- round(dt.srcs[[paste0(prop, "_sd")]] , fingR::lvl.signif(dt.sources[[prop]]))
}

dt.srcs.2 <- dt.sources %>%
  dplyr::select(Class, all_of(tracers.all))
dt.srcs.lg <- dt.srcs.2
dt.srcs.lg[2:ncol(dt.srcs.2)] <- log(dt.srcs.2[2:ncol(dt.srcs.2)])

dt.srcs.lg <- dt.srcs.lg %>%
  group_by(Class) %>%
  summarise(across(everything(), list(Mean=mean, sd=sd)))
dt.srcs.lg[2:ncol(dt.srcs.lg)] <- round(dt.srcs.lg[2:ncol(dt.srcs.lg)], 4)


dt.srcs
dt.srcs.lg


rect.alpha <- 0.15
sources.col <- c("#ff6600", "#85cf00", "#003e82")
```

faire même code pour non log
```{r}
tracers.all.err <- c("TOC_SD", "TN_SD", "Al_RMSE", "Ca_RMSE", "K_RMSE", "Si_RMSE", "Sr_RMSE", "Zn_RMSE", "Zr_RMSE", "L_star_D65_SD", "a_star_D65_SD", "b_star_D65_SD", "C_star_D65_SD", "A2_uncertainty", "Q7.4_uncertainty", "Goethite_525nm_FstD_uncertainty", "h_D65_SD", "A1_uncertainty", "Gt_uncertainty")

for(i in 1:length(tracers.all)){
  trcr <- tracers.all[i]
  trcr.err <- tracers.all.err[i]
  
  plot.graph <- ggplot(dt.graph) + aes(x = depth_cm) +
  # source properties
  geom_segment(data = dt.srcs.lg, aes(x = 6, xend = 38, y = .data[[paste0(trcr, "_Mean")]], yend = .data[[paste0(trcr, "_Mean")]]), colour = sources.col, linewidth = 0.5) +
  geom_rect(data = dt.srcs.lg, aes(xmin = 6, xmax = 38, ymin = .data[[paste0(trcr, "_Mean")]] - .data[[paste0(trcr, "_sd")]],  ymax = .data[[paste0(trcr, "_Mean")]] + .data[[paste0(trcr, "_sd")]]), fill = sources.col, alpha = rect.alpha, inherit.aes = FALSE) +
  #geom_text(data = dt.srcs.lg, aes(x = 5.2, y = .data[[paste0(trcr, "_Mean")]], label = Class), colour = sources.col) +
  # error
  geom_ribbon(aes(ymin = log(.data[[trcr]] - .data[[trcr.err]]), ymax = log(.data[[trcr]] + .data[[trcr.err]])), alpha = .7, fill = "#f4f8fb") +
  # Measure point and line
  geom_point(aes(y = log(.data[[trcr]])), pch = 16) +
  geom_line(aes(y = log(.data[[trcr]])), linewidth = 0.45) +
  # Flip the graph
  scale_x_reverse("", limits = c(38, 6), breaks = seq(7, 38, 2), sec.axis = dup_axis()) + 
  scale_y_continuous("", position = "right", sec.axis = dup_axis()) +
  #scale_y_continuous("", limits = c(min(dt.srcs.lg[[paste0(trcr, "_Mean")]] - dt.srcs.lg[[paste0(trcr, "_sd")]]), max(dt.srcs.lg[[paste0(trcr, "_Mean")]] + dt.srcs.lg[[paste0(trcr, "_sd")]])), position = "right", sec.axis = dup_axis()) +
  coord_flip() +
  # THEME
  theme_clean() +
  theme(legend.position = "none", axis.title = element_text(size = 15), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 11), #, angle = 45, hjust = 1
        panel.grid = element_line(linewidth = 0.5), axis.ticks = element_line(linewidth = 0.3), plot.background = element_rect(fill='transparent', color=NA))

#ggsave(filename = "Core_TOC.pdf", plot = plot.graph, path = dir.core, height = 17, width = 9, units = c("cm"))
ggsave(filename = paste0("Core_signature_", trcr, ".png"), plot = plot.graph, path = dir.core, height = 17, width = 9, units = c("cm"), dpi = 200)
}
```

