knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = 'C:\\Users\\camie\\Dropbox\\Dropbox\\Beg ar Loued\\Article AS_CP\\BAL R data and scripts\\clean_R_scripts\\data_for_figures')

##This R Markdown document provides the code for all figures in Prevost and Suryanarayan et al. to ensure transparency and reproducibility of data.

INSTRUCTIONS

#1) Please ensure you have downloaded the latest version of R Studio.

#2) Download data files from the data repository.

#3) Change file pathway to the name of the folder where you have saved the data files.

##publication order:

#FIG. 5

#FIG. 6

FIG. 5

#Load packages

library(readr)
library(readxl)
library(ggplot2)
library(patchwork)
library(ggrepel)

#Load data and verification

Pathway <- 'C:/Users/camie/Dropbox/Dropbox/Beg ar Loued/Article AS_CP/BAL R data and scripts/clean_R_scripts/data_for_figures'
setwd(Pathway)
tt <- read_excel ("data_fig4_iso_v2.xls", sheet = "archeo")
ref <- read_excel("data_fig4_iso_refs.xls", sheet = "ref")
dim (tt)
## [1] 24 12
head(ref)
head(tt)

#Convert variables to factors

tt$nomech <- as.factor(tt$nomech)
tt$gr <- as.factor(tt$gr)
tt$new_chrono <- as.factor(tt$US2)

#Define x, y and z

x=tt$d13C16MEcorr
y=tt$d13C18MEcorr
z=expression(list(Delta ^"13" ~"C"))
zz=expression("Sample identification clusters")
zzz=expression("Published references")
g=-3.3

#Preparing Figure 5

n3<-ggplot(tt, aes(x, y))  + 
    stat_density_2d(data=ref, 
                  aes(x=d13C16MEcorr,y=d13C18MEcorr,
                  alpha = (..level..) ^ 2, fill = cat), 
                  geom = "polygon", contour_var = "ndensity", bins=7)+
  scale_alpha_continuous(name=zz,range = c(0, 0.7))+
  guides(fill = guide_legend(title = zzz))+
  scale_fill_manual(values=c('turquoise2','#8cb9b2', 'royalblue','red','lightgreen','#ff89dc'))+
  labs(x=expression(list(delta^"13" ~"C"[16:0~MEcorr])), 
       y=expression(list(delta^"13" ~"C"[18:0~MEcorr])),
       size=1.5)+
  scale_y_continuous(limits=c(-40,-20))+
  theme_minimal() +
  facet_wrap(facets = vars(US2), strip.position = c("top"))+
  theme(legend.title = element_text(size = 10),
        strip.text = element_text(size=12, color = 'black'),
        strip.placement = "outside",
        strip.background = element_rect(fill = "#eaecee", colour = "black", size = 0.5)
        )
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#Preparing Figure 5: add reference data points

n4<-n3+ geom_point(aes(color=Delta13CcorrME, shape=gr), alpha=1, size=2.5)+
  scale_shape_manual(name=zz,values=c(15,17,19,18), labels=c("Milk","Ruminant", "Unspecific"))+
  geom_vline(xintercept = -30, linetype = "dotted",colour = "black") + 
  geom_vline(xintercept = -27, linetype = "dashed",colour = "grey", size=0.15) + 
  geom_vline(xintercept = -24, linetype = "dotted",colour = "black") + 
  scale_colour_gradient2(name=z,low = "#fc49ff",mid = "#05b400",high = "#111218",midpoint =-1.5)+
  geom_rug(aes(color=Delta13CcorrME))+ 
  scale_x_continuous(expand = expansion(0))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#Preparing Figure 5: add colours and annotations

n5<- n4+
  annotate(geom  ="rect", xmin = -40, xmax = -27, ymin = -40, ymax = -39,alpha = 0.1, fill = "green")
n6<- n5+
  annotate(geom  ="rect", xmin = -27, xmax = -20, ymin = -40, ymax = -39, alpha = .1, fill = "dodger blue")
n7<- n6+
  annotate("text", x = c(-23), y = (-39.5),label = c(">>Incr. Marine diet"),colour = c("royalblue"), size=3)
n8<- n7+
  annotate("text", x = c(-33), y = (-39.5),label = c("Incr.C3 diet <<"),colour = c("palegreen4"), size=3)
n9<- n8+
  annotate(geom  ="rect", xmin = -24, xmax = -20, ymin = -39, ymax = -38,alpha = 0.2, fill = "blue")
n10<- n9+
  annotate(geom  ="rect", xmin = -40, xmax = -27, ymin = -39, ymax = -38, alpha = .2, fill = "turquoise2")
n11<- n10+
  annotate(geom  ="rect", xmin = -30, xmax = -26, ymin = -39, ymax = -38, alpha = .2, fill = "rosybrown4")
n12<- n11+
  annotate("text", x = c(-34), y = (-38.5),label = c("Freshwater"),colour = c("turquoise4"),size=3)
n13<- n12+
  annotate("text", x = c(-22), y = (-38.5),label = c("Marine"),colour = c("royalblue3"),size=3)
n14<- n13+
  annotate("text", x = c(-27), y = (-38.5),label = c("Terrestrial"),colour = c("rosybrown4"),size=3)
n14
## Warning: The dot-dot notation (`..level..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(level)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 9 rows containing non-finite values (`stat_density2d()`).

FIG. 6

#Load packages

library(dplyr) 
## 
## Attachement du package : 'dplyr'
## Les objets suivants sont masqués depuis 'package:stats':
## 
##     filter, lag
## Les objets suivants sont masqués depuis 'package:base':
## 
##     intersect, setdiff, setequal, union
library(plyr)
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attachement du package : 'plyr'
## Les objets suivants sont masqués depuis 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ plyr::arrange()   masks dplyr::arrange()
## ✖ purrr::compact()  masks plyr::compact()
## ✖ plyr::count()     masks dplyr::count()
## ✖ plyr::desc()      masks dplyr::desc()
## ✖ plyr::failwith()  masks dplyr::failwith()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ plyr::id()        masks dplyr::id()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ plyr::mutate()    masks dplyr::mutate()
## ✖ plyr::rename()    masks dplyr::rename()
## ✖ plyr::summarise() masks dplyr::summarise()
## ✖ plyr::summarize() masks dplyr::summarize()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggthemes)

#Load data and verification

Pathway <- "C:/Users/camie/Dropbox/Dropbox/Beg ar Loued/Article AS_CP/BAL R data and scripts/clean_R_scripts/data_for_figures"
setwd(Pathway)
df <- read_excel("data_fig_supp_2.xls", sheet = "slip")
df$type <- factor(df$type)
df_sorted <- arrange(df, id) 
str(df)
## tibble [14 × 6] (S3: tbl_df/tbl/data.frame)
##  $ id     : chr [1:14] "MR6507" "MR6507" "MR6507" "MR6510" ...
##  $ nom_ech: chr [1:14] "1 (SE)" "2 (T)" "3 (SI)" "1 (SE)" ...
##  $ qt     : num [1:14] 4368 452 678 303 109 ...
##  $ type   : Factor w/ 2 levels "coating","surface": 1 1 1 2 2 1 1 1 1 1 ...
##  $ qtprc  : num [1:14] 100 10.3 15.5 100 36.1 ...
##  $ type2  : chr [1:14] "SE" "T" "SI" "SE" ...
plot(df)

#Plotting Figure 6: Geom_bar chart

n5<-ggplot(df_sorted, aes (x=nom_ech, y=qtprc, fill=type2))+
  scale_fill_manual(values=c("#9d91a6", "#4d3560", "#6a5c76"))+
  geom_bar(stat= "identity", show.legend = FALSE, width = 0.80)+
  scale_y_continuous(limits=c(0,100))+
  facet_wrap(~id, ncol=6, scales = "free") +
  labs(title = "TLE of samples from sherds with surface/coating layers", y = "Percentage")+
  xlab("")+
  theme_minimal()
leg<-element_text(color = "black", size = 6, angle=90, vjust = 0.4)
n6<-n5+theme(axis.text.x =leg, legend.background = element_rect(colour = NA))
n7<-n6+theme(panel.grid.minor = element_line(colour="white"))
n7