#Title: "Smaller than a needle! Assessing technological behaviour and contextual integrity through lithic refittings at Cova Gran de Santa Linya (SE Pyrenees)"
#subtitle: Supplementary Material (SM)
#author: Javier Sánchez-Martínez*, Jorge Martínez-Moreno, Alfonso Benito-Calvo
#* Corresponding author: Javier Sánchez-Martínez (jsanchez@iphes.cat)



library(tidyverse)
library(readxl)
library(janitor)
library(RColorBrewer)
library(ggplot2)
library(circular)


# ------------------------ IMPORT and CLEANING --------------------

CG_7P_refits = read_excel("CG_7P_Ref.xlsx")

CG_7P_refits$artefact [CG_7P_refits$artefact == 'lit'] = 'lithic'
CG_7P_refits$artefact [CG_7P_refits$artefact == 'ose'] = 'bone'

CG_7P_refits = CG_7P_refits %>%
  filter(!artefact == "bone")

CG_7P_refits$category [CG_7P_refits$category == 'bpf'] = 'Frag'
CG_7P_refits$category [CG_7P_refits$category == 'bn1g'] = 'Cores'
CG_7P_refits$category [CG_7P_refits$category == 'bna'] = 'Cobbles'
CG_7P_refits$category [CG_7P_refits$category == 'bne'] = 'Hammerstone'
CG_7P_refits$category [CG_7P_refits$category == 'bp'] = 'Complete_blanks'
CG_7P_refits$category [CG_7P_refits$category == 'bpi'] = 'Indet_Frag.'
CG_7P_refits$category [CG_7P_refits$category == 'bp2g'] = 'Debris'
CG_7P_refits$category [CG_7P_refits$category == 'bn2g'] = 'Retouched'


CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'bpf'] = 'Flake_frag.'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'frg'] = 'Nodule'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'bp'] = 'Flake'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'bp(l)'] = 'Blade'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'canto'] = 'Cobble'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'fl'] = 'Blade frag.'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'laminar'] = 'Blade'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'lasca'] = 'Flake'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'ld'] = 'Distal_frag.'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'lm'] = 'Medial_frag.'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'lp'] = 'Proximal_frag.'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 's.def'] = 'Fragments'
CG_7P_refits$blank_type [CG_7P_refits$blank_type == 'reb'] = 'Flake'


CG_7P_refits$module [CG_7P_refits$module == 'F.Lamina'] = 'Blade_frag'
CG_7P_refits$module [CG_7P_refits$module == 'F.Laminita'] = 'Bladelet_frag'
CG_7P_refits$module [CG_7P_refits$module == 'Lamina'] = 'Blade'
CG_7P_refits$module [CG_7P_refits$module == 'Laminita'] = 'Bladelet'
CG_7P_refits$module [CG_7P_refits$module == 'Lasca'] = 'Flake'
CG_7P_refits$module [CG_7P_refits$module == 'Lasca laminar'] = 'Elongated_flake'
CG_7P_refits$module [CG_7P_refits$module == 'Sin clasificar'] = 'Flake_frag'
CG_7P_refits$module [CG_7P_refits$module == 'NA'] = 'Flake_frag'


CG_7P_refits$raw_material [CG_7P_refits$raw_material == 'sil'] = 'Flint'
CG_7P_refits$raw_material [CG_7P_refits$raw_material == 'cta'] = 'Quartzite'

CG_7P_refits = CG_7P_refits %>%
  mutate(RMG = raw_material_cat) %>%
  mutate(RMG = replace(RMG, raw_material_cat %in% c('a--002','a--003'), 'AA' )) %>%
  mutate(RMG = replace(RMG, raw_material_cat %in% c('m--201','m--203','m--204','m--203','m--205','m--207','m--208','m--220','m--2109','m--nc'), 'M' ))%>%
  mutate(RMG = replace(RMG, raw_material_cat %in% c('d--202','e--202','e--203','g--001','g--203','g--206','g--207','g-208','g--210','g--211','m--209','g--nc'), 'G' ))%>%
  mutate(RMG = replace(RMG, raw_material_cat %in% c('j--204','j--205'), 'J'))%>%
  mutate(RMG = replace(RMG, raw_material_cat %in% c('sll--1','sll--3','sll--4'), 'Lacustrine' )) %>%
  mutate(RMG = replace(RMG, raw_material_cat %in% c('a--001','a--201','a--206','a--207','a--216','g--201'), 'A' )) %>%
  mutate(RMG = replace(RMG, raw_material_cat %in% c('nc', 'nc--c', 's.def'), 'Other' )) 


CG_7P_refits$tech_type [CG_7P_refits$tech_type == 'cburil 1'] = 'burin_crest1'
CG_7P_refits$tech_type [CG_7P_refits$tech_type == 'cburil 2'] = 'burin_crest2'
CG_7P_refits$tech_type [CG_7P_refits$tech_type == 'gburil'] = 'burin_spall'
CG_7P_refits$tech_type [CG_7P_refits$tech_type == 'semicresta'] = 'semicrest'
CG_7P_refits$tech_type [CG_7P_refits$tech_type == 'semitableta'] = 'semitablet'
CG_7P_refits$tech_type [CG_7P_refits$tech_type == 'tb'] = 'burin_tablet'


CG_7P_refits = CG_7P_refits %>%
  mutate(retouch_group = retouch_type) %>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('a11','a13','a23'), 'Abrupt' )) %>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('b11','b12','b21','b22','b23','b32'), 'Burin' ))%>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('bc1','bc2'), 'Bec' ))%>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('g11','g12','g21','g23','g31'), 'End-scraper' ))%>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('ld11','ld12','ld21','ld22'), 'Armatures' ))%>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('r1','r2','r4'), 'Scrapers' ))%>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('d1','d2','d3'), 'Denticulates and notches'))%>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('e1'), 'Ecaillé')) %>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('t1','t1x','t2','t2x'), 'Truncatures')) %>%
  mutate(retouch_group = replace (retouch_group, retouch_type %in% c('NA'), 'NA'))

CG_7P_refits$ref_class = as.character(CG_7P_refits$ref_class)


#Creating observations 

x_inclination = CG_7P_refits %>%
  filter(!is.na(x_inclination))

y_inclination = CG_7P_refits %>%
  filter(!is.na(y_inclination)) 

incxmo = CG_7P_refits %>%
  filter(!is.na(incxmo)) 

incymo = CG_7P_refits %>%
  filter(!is.na(incymo)) 

horizontal_dist = CG_7P_refits %>%
  filter(!is.na(horizontal_dist)) 

azimuth = CG_7P_refits %>%
  filter(!is.na(azimuth)) 

axis = CG_7P_refits %>%
  filter(!is.na(axis)) 

inclination = CG_7P_refits %>%
  filter(!is.na(inclination)) 

vertical_dist = CG_7P_refits %>%
  filter(!is.na(vertical_dist)) 



ref_1 = CG_7P_refits %>%
  filter(ref_class == "1")

ref_2 = CG_7P_refits %>%
  filter(ref_class == "2")

ref_3 = CG_7P_refits %>%
  filter(ref_class == "3")

ref_4 = CG_7P_refits %>%
  filter(ref_class == "4")

ref_class_clean = CG_7P_refits %>%
  filter(ref_class %in% c("1", "2", "4"))



# ----------------------- STATISTICAL ANALYSIS ------------------------


summary(CG_7P_refits)

# NORMALITY TEST with METRIC VARIABLES 
# Shapiro–Wilk test 

shapiro.test(CG_7P_refits$vertical_dist)
shapiro.test(CG_7P_refits$horizontal_dist)
shapiro.test(CG_7P_refits$inclination)

# NORMALITY TEST with CIRCULAR VARIABLES 
# Rayleigh test and Kuiper test
CG_7P_refits_az <- circular(CG_7P_refits$azimuth,
                    units = "degrees",
                    template = "geographics",
                    modulo = "2pi")

rayleigh.test(CG_7P_refits_az)

kuiper.test(CG_7P_refits_az)

# QQ plots for metric variables

ggplot(CG_7P_refits, aes(sample = horizontal_dist)) +
  stat_qq() +
  stat_qq_line(col = "blue") +
  theme_bw()

ggplot(CG_7P_refits, aes(sample = vertical_dist)) +
  stat_qq() +
  stat_qq_line(col = "orange") +
  theme_bw()

ggplot(CG_7P_refits, aes(sample = inclination)) +
  stat_qq() +
  stat_qq_line(col = "green3") +
  theme_bw()

#All together 

p = CG_7P_refits %>%
  select(horizontal_dist, vertical_dist, inclination) %>%
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "value") %>%
  ggplot(aes(sample = value)) +
  stat_qq(size = 1, alpha = 0.6) +
  stat_qq_line(colour = "red3", linewidth = 0.8) +
  facet_wrap(~ variable, scales = "free", ncol = 3) +
  theme_bw() +
  theme(strip.text = element_text(size = 12, face = "bold"),
        axis.title = element_blank())
print(p)

ggsave("QQ plots.png", plot = p, bg = "white", dpi = 300)
ggsave("QQ plots.tiff", plot = p, bg = "white", dpi = 300)



# Non-parametric tests 
# Kruskal-Wallis

kruskal.test(horizontal_dist ~ ref_class, data = ref_class_clean)
kruskal.test(vertical_dist ~ ref_class, data = ref_class_clean)
kruskal.test(inclination ~ ref_class, data = ref_class_clean)



# Watson two-sample

az_ref1 = circular(ref_1$azimuth,
                    units = "degrees",
                    template = "geographics",
                    modulo = "2pi")

az_ref2 = circular(ref_2$azimuth,
                    units = "degrees",
                    template = "geographics",
                    modulo = "2pi")

az_ref4 = circular(ref_3$azimuth,
                   units = "degrees",
                   template = "geographics",
                   modulo ="2pi")

watson.two.test(az_ref1, az_ref2)
watson.two.test(az_ref1, az_ref4)
watson.two.test(az_ref2, az_ref4)



# --------------------------- DATA ANALYSIS-----------------------------

# LINEAR METRIC VARIABLES
# Horizontal dist

p = ggplot(CG_7P_refits, aes(x = horizontal_dist)) +
  geom_histogram(aes(y = after_stat(count / sum(count) * 100),
                     fill = after_stat(count)),
                 bins = 20,
                 colour = "black") +
  scale_fill_distiller(palette = "Blues", direction = 1) +
  labs(title = "Horizontal distance between refitted pieces",
       x = "Horizontal distance (mm)",
       y = "Percentage (%)",
       fill = "Frequency") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))
print(p)

ggsave("7P_horizontal_dist.png", plot = p, bg = "white", dpi = 300)
ggsave("7P_horizontal_dist.tiff", plot = p, bg = "white", dpi = 300)


# Vertical dist

p = ggplot(CG_7P_refits, aes(x = vertical_dist)) +
  geom_histogram(aes(y = after_stat(count / sum(count) * 100),
                     fill = after_stat(count)),
                 bins = 20,
                 colour = "black") +
  scale_fill_distiller(palette = "Oranges", direction = 1) +
  labs(title = "Vertical distance between refitted pieces",
       x = "Vertical distance (mm)",
       y = "Percentage (%)",
       fill = "Frequency") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))
print(p)

ggsave("7P_vertical_dist.png", plot = p, bg = "white", dpi = 300)
ggsave("7P_vertical_dist.tiff", plot = p, bg = "white", dpi = 300)


# Inclination

p = ggplot(CG_7P_refits, aes(x = inclination)) +
  geom_histogram(aes(y = after_stat(count / sum(count) * 100),
                     fill = after_stat(count)),
                 bins = 30,
                 colour = "black") +
  scale_fill_distiller(palette = "Greens", direction = 1) +
  labs(title = "Inclination between refitted pieces",
       x = "Inclination (º)",
       y = "Percentage (%)",
       fill = "Frequency") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))
print(p)

ggsave("7P_inclination.png", plot = p, bg = "white", dpi = 300)
ggsave("7P_inclination.tiff", plot = p, bg = "white", dpi = 300)



#CIRCULAR METRIC VARIABLES 
#REF 1

p = ggplot(ref_1, aes(x = azimuth, y = after_stat(count / sum(count)))) +
  geom_histogram(breaks = seq(0, 360, by = 10), boundary = 0, closed = "left",
                 fill = "red4",colour = "black") +
  coord_polar(start = 0, direction = +1) +  
  scale_x_continuous(limits = c(0, 360),breaks = seq(0, 330, by = 30)) +
  labs(title = "Rose diagram – REF 1", x = NULL, y = "Frequency") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
print(p)

ggsave("RD_ref1.png", plot = p, bg = "white", dpi = 300)
ggsave("RD_ref1.tiff", plot = p, bg = "white", dpi = 300)


#REF 2

p = ggplot(ref_2, aes(x = azimuth, y = after_stat(count / sum(count)))) +
  geom_histogram(breaks = seq(0, 360, by = 10), boundary = 0, closed = "left",
                 fill = "green3",colour = "black") +
  coord_polar(start = 0, direction = +1) +  
  scale_x_continuous(limits = c(0, 360),breaks = seq(0, 330, by = 30)) +
  labs(title = "Rose diagram – REF 2", x = NULL, y = "Frequency") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
print(p)

ggsave("RD_ref2.png", plot = p, bg = "white", dpi = 300)
ggsave("RD_ref2.tiff", plot = p, bg = "white", dpi = 300)


#REF 3

p = ggplot(ref_3, aes(x = azimuth, y = after_stat(count / sum(count)))) +
  geom_histogram(breaks = seq(0, 360, by = 10), boundary = 0, closed = "left",
                 fill = "grey",colour = "black") +
  coord_polar(start = 0, direction = +1) +  
  scale_x_continuous(limits = c(0, 360),breaks = seq(0, 330, by = 30)) +
  labs(title = "Rose diagram – REF 3", x = NULL, y = "Frequency") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
print(p)

ggsave("RD_ref3.png", plot = p, bg = "white", dpi = 300)
ggsave("RD_ref3.tiff", plot = p, bg = "white", dpi = 300)



#REF 4

p = ggplot(ref_4, aes(x = azimuth, y = after_stat(count / sum(count)))) +
  geom_histogram(breaks = seq(0, 360, by = 10), boundary = 0, closed = "left",
                 fill = "orange2",colour = "black") +
  coord_polar(start = 0, direction = +1) +  
  scale_x_continuous(limits = c(0, 360),breaks = seq(0, 330, by = 30)) +
  labs(title = "Rose diagram – REF 4", x = NULL, y = "Frequency") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
print(p)

ggsave("RD_ref4.png", plot = p, bg = "white", dpi = 300)
ggsave("RD_ref4.tiff", plot = p, bg = "white", dpi = 300)