# R script for the publication
# Phillipps, R. et al. 2021. Quantification of stone artefacts assemblages in New Zealand. Submitted to Journal of Pacific Archaeology.

library(tidyverse)
library(ggpubr)
library(janitor)
library(fitdistrplus)


# Read in data ------------------------------------------------------------------------------------

tauroa_data <- read_csv("Tauroa.csv")
temataku_data <- read_csv("TeMataku.csv")



# Quantification comparison tables ----------------------------------------------------------------

# NAS, MNF, Average Weighted Method, Total Weight, Total Max Length

# Tauroa NAS
nas_tauroa <- tauroa_data %>% 
  tabyl(MAT.TYPE) %>% 
  adorn_totals("row") %>% 
  adorn_pct_formatting(digits = 1, affix_sign = FALSE)

# Tauroa MNF
mnf_tauroa <- tauroa_data %>% 
  group_by(MAT.TYPE) %>% 
  summarise(mnf1 = sum(c(FCLASS == "COMPFLAKE", FCLASS == "COMPTOOL", FCLASS == "DISTFLAKE", FCLASS == "DISTTOOL")),
            mnf2 = sum(c(FCLASS == "COMPSPLIT", FCLASS == "COMPSPLITTOOL", FCLASS == "PROXSPLIT", FCLASS == "PROXSPLITTOOL")) / 2) %>% 
  mutate(MNF = mnf1 + mnf2) %>% 
  dplyr::select(MAT.TYPE, MNF) %>% 
  mutate(MNF_percent = MNF / sum(MNF) * 100) %>% 
  adorn_totals("row")

# Tauroa average weighted method
av_wtd_method_tauroa <- tauroa_data %>% 
  mutate(COMP = case_when(FCLASS == "CORE" ~ "Core", 
                          FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "COMPBIPOLAR" ~ "Complete",
                          FCLASS != "COMPFLAKE" | FCLASS != "COMPTOOL" | FCLASS != "COMPBIPOLAR" | FCLASS != "CORE" ~ "Broken")) %>% 
  filter(COMP != "Core") %>% 
  group_by(MAT.TYPE, COMP) %>% 
  summarise(count = n(), total_wt = sum(WEIGHT, na.rm = TRUE), av_wt = mean(WEIGHT, na.rm = TRUE)) %>% 
  pivot_wider(names_from = COMP, values_from = c(count, total_wt, av_wt)) %>% 
  mutate(orgFE = round(count_Complete + (total_wt_Broken / av_wt_Complete))) %>% 
  dplyr::select(MAT.TYPE, count_Complete, total_wt_Broken, av_wt_Complete, orgFE) %>% 
  adorn_totals("row")

# Tauroa total weight
weight_tauroa <- tauroa_data %>% 
  group_by(MAT.TYPE) %>% 
  summarise(count = n(), total_wt = sum(WEIGHT, na.rm = TRUE)) %>% 
  adorn_totals("row")

# Tauroa total max length
max_l_tauroa <- tauroa_data %>% 
  group_by(MAT.TYPE) %>% 
  summarise(count = n(), total_max_l = sum(MAXLENGTH, na.rm = TRUE)) %>% 
  adorn_totals("row")

# Tauroa volume
volume_tauroa <- tauroa_data %>% 
  mutate(density = case_when(MAT.TYPE == "BASALT" ~ 2.9,
                             MAT.TYPE == "CHERT" ~ 2.49,
                             MAT.TYPE == "OBSIDIAN" ~ 2.38)) %>% 
  mutate(volume = WEIGHT / density) %>% 
  group_by(MAT.TYPE) %>% 
  summarise(count = n(), total_volume = sum(volume, na.rm = TRUE)) %>% 
  adorn_totals("row")
  

# Combine all measures into single table
comparison_tauroa <- cbind(dplyr::select(nas_tauroa, MAT.TYPE, n),
                             dplyr::select(mnf_tauroa, MNF),
                             dplyr::select(max_l_tauroa, total_max_l),
                             dplyr::select(weight_tauroa, total_wt),
                             dplyr::select(av_wtd_method_tauroa, orgFE),
                             dplyr::select(volume_tauroa, total_volume))
colnames(comparison_tauroa) <- c("Material", "NAS", "MNF", "Total Maximum Length (mm)", "Total Weight (g)", "Average Weighted Method", "Total volume (cm3)")



# Te Mataku NAS
nas_temataku <- temataku_data %>% 
  tabyl(MAT.TYPE) %>% 
  adorn_totals("row") %>% 
  adorn_pct_formatting(digits = 1, affix_sign = FALSE)

# Te Mataku MNF
mnf_temataku <- temataku_data %>% 
  group_by(MAT.TYPE) %>% 
  summarise(mnf1 = sum(c(FCLASS == "COMPFLAKE", FCLASS == "COMPTOOL", FCLASS == "DISTFLAKE", FCLASS == "DISTTOOL")),
            mnf2 = sum(c(FCLASS == "COMPSPLIT", FCLASS == "COMPSPLITTOOL", FCLASS == "PROXSPLIT", FCLASS == "PROXSPLITTOOL")) / 2) %>% 
  mutate(MNF = mnf1 + mnf2) %>% 
  dplyr::select(MAT.TYPE, MNF) %>% 
  mutate(MNF_percent = MNF / sum(MNF) * 100) %>% 
  adorn_totals("row")

# Te Mataku average weighted method
av_wtd_method_temataku <- temataku_data %>% 
  mutate(COMP = case_when(FCLASS == "CORE" ~ "Core", 
                          FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "COMPBIPOLAR" ~ "Complete",
                          FCLASS != "COMPFLAKE" | FCLASS != "COMPTOOL" | FCLASS != "COMPBIPOLAR" | FCLASS != "CORE" ~ "Broken")) %>% 
  filter(COMP != "Core") %>% 
  group_by(MAT.TYPE, COMP) %>% 
  summarise(count = n(), total_wt = sum(WEIGHT, na.rm = TRUE), av_wt = mean(WEIGHT, na.rm = TRUE)) %>% 
  pivot_wider(names_from = COMP, values_from = c(count, total_wt, av_wt)) %>% 
  mutate(orgFE = round(count_Complete + (total_wt_Broken / av_wt_Complete))) %>% 
  dplyr::select(MAT.TYPE, count_Complete, total_wt_Broken, av_wt_Complete, orgFE) %>% 
  adorn_totals("row")

# Te Mataku total weight
weight_temataku <- temataku_data %>% 
  group_by(MAT.TYPE) %>% 
  summarise(count = n(), total_wt = sum(WEIGHT, na.rm = TRUE)) %>% 
  adorn_totals("row")

# Te Mataku total max length
max_l_temataku <- temataku_data %>% 
  group_by(MAT.TYPE) %>% 
  summarise(count = n(), total_max_l = sum(MAXLENGTH, na.rm = TRUE)) %>% 
  adorn_totals("row")

# Te Mataku volume
volume_temataku <- temataku_data %>% 
  mutate(density = case_when(MAT.TYPE == "BASALT" ~ 2.9,
                             MAT.TYPE == "CHERT" ~ 2.49,
                             MAT.TYPE == "OBSIDIAN" ~ 2.38)) %>% 
  mutate(volume = WEIGHT / density) %>% 
  group_by(MAT.TYPE) %>% 
  summarise(count = n(), total_volume = sum(volume, na.rm = TRUE)) %>% 
  adorn_totals("row")


# Combine all measures into single table
comparison_temataku <- cbind(dplyr::select(nas_temataku, MAT.TYPE, n),
                             dplyr::select(mnf_temataku, MNF),
                             dplyr::select(max_l_temataku, total_max_l),
                             dplyr::select(weight_temataku, total_wt),
                             dplyr::select(av_wtd_method_temataku, orgFE),
                             dplyr::select(volume_temataku, total_volume))
colnames(comparison_temataku) <- c("Material", "NAS", "MNF", "Total Maximum Length (mm)", "Total Weight (g)", "Average Weighted Method", "Total volume (cm3)")

write_csv(comparison_tauroa, "quantification_comparison_tauroa.csv")
write_csv(comparison_temataku, "quantification_comparison_temataku.csv")




# Quantification methods comparison plot ----------------------------------------------------------

# Combine all measures (NAS, MNF, Total Max Length, Total Weight, Average Weighted Method) in one data frame
quant_methods <- comparison_tauroa %>% 
  filter(Material != "Total") %>% 
  mutate(`Total Weight` = `Total Weight (g)`) %>% 
  mutate(`Total Length` = `Total Maximum Length (mm)`) %>% 
  mutate(`Total Volume` = `Total volume (cm3)`) %>% 
  dplyr::select(Material, NAS, MNF, `Total Length`, `Total Weight`, `Average Weighted Method`, `Total Volume`) %>% 
  pivot_longer(-Material, names_to = "method", values_to = "value") %>% 
  cbind(Assemblage = "Tauroa", .) %>% 
  rbind(comparison_temataku %>% 
          filter(Material != "Total") %>% 
          mutate(`Total Weight` = `Total Weight (g)`) %>% 
          mutate(`Total Length` = `Total Maximum Length (mm)`) %>%
          mutate(`Total Volume` = `Total volume (cm3)`) %>% 
          dplyr::select(Material, NAS, MNF, `Total Length`, `Total Weight`, `Average Weighted Method`, `Total Volume`) %>% 
          pivot_longer(-Material, names_to = "method", values_to = "value") %>% 
          cbind(Assemblage = "Te Mataku", .))
quant_methods$method <- str_replace(quant_methods$method, "Average Weighted Method", "AWM")
quant_methods$method = factor(quant_methods$method, levels = c("NAS", "MNF", "Total Length", "Total Weight", "AWM", "Total Volume"))

# Plot measures as stacked percent chart split by assemblage
levels(quant_methods$method) <- gsub(" ", "\n", levels(quant_methods$method))
ggplot(quant_methods, aes(x = method, y = value, fill = Material)) +
  geom_bar(position = "fill", stat = "identity") +
  facet_wrap(vars(Assemblage), nrow = 1, ncol = 2) +
  labs(x = "Quantification Method", y = "Proportion", fill = "") +
  theme_bw() +
  scale_fill_grey() +
  scale_y_continuous(labels = scales::percent) +
  theme(legend.position="bottom")

ggsave("quantification_comparison_figure.tiff", units="mm", width=200, height=140, dpi=600, compression = 'lzw')




# Artefact type breakdown by material -------------------------------------------------------------
artefact_types_tauroa <- tauroa_data %>%
  tabyl(FCLASS, MAT.TYPE) %>% 
  adorn_totals(c("row", "col")) %>% 
  adorn_percentages("col") %>% 
  adorn_pct_formatting(digits = 1, affix_sign = FALSE) %>% 
  adorn_ns(position = "front")
colnames(artefact_types_tauroa) <- c("Artefact Type", "Basalt", "Chert", "Obsidian", "Total")


artefact_types_temataku <- temataku_data %>%
  tabyl(FCLASS, MAT.TYPE) %>% 
  adorn_totals(c("row", "col")) %>% 
  adorn_percentages("col") %>% 
  adorn_pct_formatting(digits = 1, affix_sign = FALSE) %>% 
  adorn_ns(position = "front")
colnames(artefact_types_temataku) <- c("Artefact Type", "Basalt", "Chert", "Obsidian", "Total")


write_csv(artefact_types_tauroa, "artefact_types_tauroa.csv")
write_csv(artefact_types_temataku, "artefact_types_temataku.csv")




# Fragmentation ratio (complete and proximal artefacts) -------------------------------------------

# Complete = complete flakes + complete tools
# Proximal = proximal flakes + proximal tools
frag_ratio_tauroa <- tauroa_data %>% 
  filter(FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "PROXFLAKE" | FCLASS == "PROXTOOL") %>% 
  tabyl(MAT.TYPE, FCLASS) %>% 
  mutate(Complete = COMPFLAKE + COMPTOOL, Proximal = PROXFLAKE + PROXTOOL, Ratio = Complete / Proximal) %>% 
  mutate(Material = MAT.TYPE) %>% 
  dplyr::select(Material, Complete, Proximal, Ratio)


frag_ratio_temataku <- temataku_data %>% 
  filter(FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "PROXFLAKE" | FCLASS == "PROXTOOL") %>% 
  tabyl(MAT.TYPE, FCLASS) %>% 
  mutate(Complete = COMPFLAKE + COMPTOOL, Proximal = PROXFLAKE + PROXTOOL, Ratio = Complete / Proximal) %>% 
  mutate(Material = MAT.TYPE) %>% 
  dplyr::select(Material, Complete, Proximal, Ratio)


write_csv(frag_ratio_tauroa, "frag_ratio_tauroa.csv")
write_csv(frag_ratio_temataku, "frag_ratio_temataku.csv")




# Weibull analysis --------------------------------------------------------------------------------

# Lin et al. (2016) analyse max length of flakes with platform: COMPFLAKE, PROXFLAKE, COMPSPLIT, PROXSPLIT
# Therefore Tauroa and Te Mataku data are limited to these classes for comparability (plus tools)

# ~ Plot separate Weibull curves for Tauroa and Te Mataku data ------------------------------------

# Obtain complete and proximal flakes, tools and splits then set up for fitting (subtract 19.99 from max dimension)
tauroa_comp <- tauroa_data %>%
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>% 
  filter(MAXDIM >= 20) %>% 
  filter(FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "PROXFLAKE" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLIT" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLIT" | FCLASS == "PROXSPLITTOOL") %>% 
  dplyr::select(MAT.TYPE, MAXDIM) %>% 
  cbind(site = "Tauroa") %>% 
  unite("assem", c(site, MAT.TYPE), sep = " ") %>% 
  mutate(MAXDIM2 = MAXDIM - 19.99)


temataku_comp <- temataku_data %>%
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>% 
  filter(MAXDIM >= 20) %>% 
  filter(FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "PROXFLAKE" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLIT" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLIT" | FCLASS == "PROXSPLITTOOL") %>% 
  dplyr::select(MAT.TYPE, MAXDIM) %>% 
  cbind(site = "Te Mataku") %>% 
  unite("assem", c(site, MAT.TYPE), sep = " ") %>% 
  mutate(MAXDIM2 = MAXDIM - 19.99)


# Fit Weibull distribution and obtain shape and scale estimates
tauroa_basalt_weibull <- fitdist(tauroa_comp$MAXDIM2[which(tauroa_comp$assem == "Tauroa BASALT")], "weibull") 
tauroa_chert_weibull <- fitdist(tauroa_comp$MAXDIM2[which(tauroa_comp$assem == "Tauroa CHERT")], "weibull") 
tauroa_obsidian_weibull <- fitdist(tauroa_comp$MAXDIM2[which(tauroa_comp$assem == "Tauroa OBSIDIAN")], "weibull")

temataku_basalt_weibull <- fitdist(temataku_comp$MAXDIM2[which(temataku_comp$assem == "Te Mataku BASALT")], "weibull") 
temataku_chert_weibull <- fitdist(temataku_comp$MAXDIM2[which(temataku_comp$assem == "Te Mataku CHERT")], "weibull") 
temataku_obsidian_weibull <- fitdist(temataku_comp$MAXDIM2[which(temataku_comp$assem == "Te Mataku OBSIDIAN")], "weibull") 


# Plot Weibull curves
tiff("weibull_curves.tiff", width = 200, height = 140, units = 'mm', res = 600)

par(mfrow = c(2, 3))

denscomp(tauroa_basalt_weibull, xlab = "Maximum Dimension (mm)", main = "Tauroa Basalt", addlegend = FALSE)
denscomp(tauroa_chert_weibull, xlab = "Maximum Dimension (mm)", main = "Tauroa Chert", addlegend = FALSE)
denscomp(tauroa_obsidian_weibull, xlab = "Maximum Dimension (mm)", main = "Tauroa Obsidian", addlegend = FALSE)
denscomp(temataku_basalt_weibull, xlab = "Maximum Dimension (mm)", main = "Te Mataku Basalt", addlegend = FALSE)
denscomp(temataku_chert_weibull, xlab = "Maximum Dimension (mm)", main = "Te Mataku Chert", addlegend = FALSE)
denscomp(temataku_obsidian_weibull, xlab = "Maximum Dimension (mm)", main = "Te Mataku Obsidian", addlegend = FALSE)

dev.off()


# ~ Shape and scale graph with Lin et al. (2016) data ---------------------------------------------

# Read in Lin et al. (2016) data
lin_chert_data <- read_csv("experimental_chert_platform_flakes_only.csv")
lin_chert_data$Event_ID = paste(lin_chert_data$Event, "-", lin_chert_data$Core.ID, sep="")

lin_obsidian_data <- read_csv("experimental_obsidian_platform_flakes_only.csv")
lin_obsidian_data <- filter(lin_obsidian_data, Event_ID != 53)


# Create functions for obtaining Weibull estimates

# Function for obtaining shape estimates
weibull_shape <- function (data) {
  a = fitdist(data, "weibull")
  shape = a$estimate[1]
  return(shape)
}

# Function for obtaining scale estimates
weibull_scale <- function (data) {
  a = fitdist(data, "weibull")
  scale = a$estimate[2]
  return(scale)
}



# Get Weibull shape and scale estimates for each reduction event in Lin et al. (2016) data

# ~ Chert data set

# Split data by reduction event
lin_chert_data_split <- split(lin_chert_data, lin_chert_data$Event_ID)

# Apply Weibull functions over each reduction event
chert_shape <- sapply (lin_chert_data_split, function (lin_chert_data) weibull_shape(lin_chert_data$MAXLENGTH - 24.99))
chert_scale <- sapply (lin_chert_data_split, function (lin_chert_data) weibull_scale(lin_chert_data$MAXLENGTH - 24.99))

# Combine into data frame
experimental_chert_weibull <- as_tibble(cbind(assem = "Experimental CHERT", shape = chert_shape, scale = chert_scale))



# ~ Obsidian data set

# Split data by reduction event
lin_obsidian_data_split <- split(lin_obsidian_data, lin_obsidian_data$Event_ID)

# Apply Weibull functions over each reduction event
obsidian_shape <- sapply (lin_obsidian_data_split, function (lin_obsidian_data) weibull_shape(lin_obsidian_data$MAXLENGTH - 24.99))
obsidian_scale <- sapply (lin_obsidian_data_split, function (lin_obsidian_data) weibull_scale(lin_obsidian_data$MAXLENGTH - 24.99))

# Combine into data frame
experimental_obsidian_weibull <- as_tibble(cbind(assem = "Experimental OBSIDIAN", shape = obsidian_shape, scale = obsidian_scale))



# Get Weibull shape and scale estimates for Tauroa and Te Mataku raw materials

# ~ Tauroa data set

# Filter data to only contain COMPFLAKE, PROXFLAKE, COMPSPLIT, PROXSPLIT plus tools
tauroa_comp_prox <- tauroa_data %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>% 
  filter(FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "PROXFLAKE" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLIT" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLIT" | FCLASS == "PROXSPLITTOOL") %>% 
  dplyr::select(MAT.TYPE, MAXDIM) %>% 
  cbind(site = "Tauroa") %>% 
  unite("assem", c(site, MAT.TYPE), sep = " ")

# Split data by raw material type
tauroa_comp_prox_split <- split(tauroa_comp_prox, tauroa_comp_prox$assem)

# Apply Weibull functions over each raw material type
tauroa_shape <- sapply (tauroa_comp_prox_split, function (tauroa_comp_prox) weibull_shape(tauroa_comp_prox$MAXDIM - 19.99))
tauroa_scale <- sapply (tauroa_comp_prox_split, function (tauroa_comp_prox) weibull_scale(tauroa_comp_prox$MAXDIM - 19.99))

# Combine into data frame
tauroa_weibull <- tibble(assem = substr(names(tauroa_shape), 1, nchar(names(tauroa_shape))-6), shape = tauroa_shape, scale = tauroa_scale)



# ~ Te Mataku data set

# Filter data to only contain COMPFLAKE, PROXFLAKE, COMPSPLIT, PROXSPLIT plus tools
temataku_comp_prox <- temataku_data %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>% 
  filter(FCLASS == "COMPFLAKE" | FCLASS == "COMPTOOL" | FCLASS == "PROXFLAKE" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLIT" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLIT" | FCLASS == "PROXSPLITTOOL") %>% 
  dplyr::select(MAT.TYPE, MAXDIM) %>% 
  cbind(site = "Te Mataku") %>% 
  unite("assem", c(site, MAT.TYPE), sep = " ")

# Split data by raw material type
temataku_comp_prox_split <- split(temataku_comp_prox, temataku_comp_prox$assem)

# Apply Weibull functions over each raw material type
temataku_shape <- sapply (temataku_comp_prox_split, function (temataku_comp_prox) weibull_shape(temataku_comp_prox$MAXDIM - 19.99))
temataku_scale <- sapply (temataku_comp_prox_split, function (temataku_comp_prox) weibull_scale(temataku_comp_prox$MAXDIM - 19.99))

# Combine into data frame
temataku_weibull <- tibble(assem = substr(names(temataku_shape), 1, nchar(names(temataku_shape))-6), shape = temataku_shape, scale = temataku_scale)


# Plot Weibull shape and scale estimates

# ~ First, combine all data into single data frame
all_weibull_estimates <- rbind(experimental_chert_weibull, experimental_obsidian_weibull, tauroa_weibull, temataku_weibull)
all_weibull_estimates$shape = as.numeric(all_weibull_estimates$shape)
all_weibull_estimates$scale = as.numeric(all_weibull_estimates$scale)


# ~ Plot the data
ggplot(all_weibull_estimates, aes(x = scale, y = shape)) +
  geom_point(aes(shape = assem, color = assem), size = 3) +
  labs(x = "Weibull scale (distribution spread)", y = "Weibull shape (distribution slope)", shape = "Assemblage", color = "Assemblage") +
  theme_bw() +
  scale_shape_manual(values = c(1, 2, 15, 16, 17, 15, 16, 17)) +
  scale_color_manual(values = c("grey39", "grey39", "blue", "blue", "blue", "orange", "orange", "orange"))

ggsave("weibull_shape_scale_plot.tiff", units="mm", width=200, height=140, dpi=600, compression = 'lzw')




# Histograms for Tauroa and Te Mataku COMPFLAKE + PROXFLAKE + COMPSPLIT + PROXSPLIT + TOOLS -------
tauroa_temataku_comp_prox <- rbind(tauroa_comp_prox, temataku_comp_prox)
tauroa_temataku_comp_prox$assem = factor(tauroa_temataku_comp_prox$assem, levels = c("Tauroa BASALT", "Tauroa CHERT", "Tauroa OBSIDIAN", "Te Mataku BASALT", "Te Mataku CHERT", "Te Mataku OBSIDIAN"))

ggplot(tauroa_temataku_comp_prox, aes(x = MAXDIM)) +
  geom_histogram(binwidth = 3) +
  facet_wrap(vars(assem), nrow = 2, ncol = 3, scales = "free_y") +
  labs(x = "Maximum length (mm)", y = "Frequency") +
  theme_bw()

ggsave("comp_prox_flake_tool_splits_histograms.tiff", units="mm", width=200, height=140, dpi=600, compression = 'lzw')


# Fractal distribution analysis -------------------------------------------------------------------

# Following Brown (2001)

# Create table with fractal plot variables: 
# 1) size class intervals, 
# 2) lower bound of the interval (r)
# 3) frequency of cases
# 4) cumulative frequency starting with largest size class (N(>r))
# 5) ln(r)
# 6) ln(N(>r))

# Then create fractal plot by plotting ln(r) against ln(N(>r))


# ~ Tauroa Basalt  

# Table with fractal variables
tauroa_basalt_fractal <- tauroa_data %>%
  # Filter material and size > 20 mm
  filter(MAT.TYPE == "BASALT") %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>%
  # Filter artefacts of interest
  filter(FCLASS == "COMPFLAKE" | FCLASS == "PROXFLAKE" | FCLASS == "COMPSPLIT" | FCLASS == "PROXSPLIT" | 
         FCLASS == "COMPTOOL" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLITTOOL") %>% 
  # Create size class intervals
  mutate(SIZE_INT = case_when(MAXDIM >= 20 & MAXDIM < 30 ~ "20-30",
                              MAXDIM >= 30 & MAXDIM < 40 ~ "30-40",
                              MAXDIM >= 40 & MAXDIM < 50 ~ "40-50",
                              MAXDIM >= 50 & MAXDIM < 60 ~ "50-60",
                              MAXDIM >= 60 ~ ">60")) %>% 
  mutate(SIZE_INT = factor(SIZE_INT, levels = c("20-30", "30-40", "40-50", "50-60", ">60"))) %>% 
  # Calculate frequencies for each size class
  dplyr::select(SIZE_INT) %>% 
  group_by(SIZE_INT) %>% 
  summarise(count = n()) %>% 
  # Add lower bound column
  cbind(lower_bound_r = c(20, 30, 40, 50, 60)) %>% 
  # Calculate cumulative frequency
  mutate("cumul_freq_N>r" = rev(cumsum(rev(count)))) %>% 
  # Calculate ln(r)
  mutate(ln_r = log(lower_bound_r)) %>% 
  # Calculate ln(N(>r))
  mutate("ln_N>r" = log(`cumul_freq_N>r`))
  
# Fractal plot
ggplot(tauroa_basalt_fractal, aes(x = ln_r, y = `ln_N>r`)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  stat_regline_equation(label.y = 3, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 2.5, aes(label = ..rr.label..)) +
  labs(x = "ln(r) (mm)", y = "ln(N(>r))") +
  theme_bw()


# ~ Tauroa Chert

# Table with fractal variables
tauroa_chert_fractal <- tauroa_data %>%
  # Filter material and size > 20 mm
  filter(MAT.TYPE == "CHERT") %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>%
  # Filter artefacts of interest
  filter(FCLASS == "COMPFLAKE" | FCLASS == "PROXFLAKE" | FCLASS == "COMPSPLIT" | FCLASS == "PROXSPLIT" |
         FCLASS == "COMPTOOL" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLITTOOL") %>% 
  # Create size class intervals
  mutate(SIZE_INT = case_when(MAXDIM >= 20 & MAXDIM < 30 ~ "20-30",
                              MAXDIM >= 30 & MAXDIM < 40 ~ "30-40",
                              MAXDIM >= 40 & MAXDIM < 50 ~ "40-50",
                              MAXDIM >= 50 & MAXDIM < 60 ~ "50-60",
                              MAXDIM >= 60 & MAXDIM < 70 ~ "60-70",
                              MAXDIM >= 70 ~ ">70")) %>% 
  mutate(SIZE_INT = factor(SIZE_INT, levels = c("20-30", "30-40", "40-50", "50-60", "60-70", ">70"))) %>% 
  # Calculate frequencies for each size class
  dplyr::select(SIZE_INT) %>% 
  group_by(SIZE_INT) %>% 
  summarise(count = n()) %>% 
  # Add lower bound column
  cbind(lower_bound_r = c(20, 30, 40, 50, 60, 70)) %>% 
  # Calculate cumulative frequency
  mutate("cumul_freq_N>r" = rev(cumsum(rev(count)))) %>% 
  # Calculate ln(r)
  mutate(ln_r = log(lower_bound_r)) %>% 
  # Calculate ln(N(>r))
  mutate("ln_N>r" = log(`cumul_freq_N>r`))

# Fractal plot
ggplot(tauroa_chert_fractal, aes(x = ln_r, y = `ln_N>r`)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  stat_regline_equation(label.y = 3, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 2.5, aes(label = ..rr.label..)) +
  labs(x = "ln(r) (mm)", y = "ln(N(>r))") +
  theme_bw()


# ~ Tauroa Obsidian

# Table with fractal variables
tauroa_obsidian_fractal <- tauroa_data %>%
  # Filter material and size > 20 mm
  filter(MAT.TYPE == "OBSIDIAN") %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>%
  # Filter artefacts of interest
  filter(FCLASS == "COMPFLAKE" | FCLASS == "PROXFLAKE" | FCLASS == "COMPSPLIT" | FCLASS == "PROXSPLIT" |
         FCLASS == "COMPTOOL" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLITTOOL") %>% 
  # Create size class intervals
  mutate(SIZE_INT = case_when(MAXDIM >= 20 & MAXDIM < 25 ~ "20-25",
                              MAXDIM >= 25 & MAXDIM < 30 ~ "25-30",
                              MAXDIM >= 30 & MAXDIM < 35 ~ "30-35",
                              MAXDIM >= 35 & MAXDIM < 40 ~ "35-40",
                              MAXDIM >= 40 & MAXDIM < 50 ~ "40-50",
                              MAXDIM >= 50 ~ ">50")) %>% 
  mutate(SIZE_INT = factor(SIZE_INT, levels = c("20-25", "25-30", "30-35", "35-40", "40-50", ">50"))) %>% 
  # Calculate frequencies for each size class
  dplyr::select(SIZE_INT) %>% 
  group_by(SIZE_INT) %>% 
  summarise(count = n()) %>% 
  # Add lower bound column
  cbind(lower_bound_r = c(20, 25, 30, 35, 40, 50)) %>% 
  # Calculate cumulative frequency
  mutate("cumul_freq_N>r" = rev(cumsum(rev(count)))) %>% 
  # Calculate ln(r)
  mutate(ln_r = log(lower_bound_r)) %>% 
  # Calculate ln(N(>r))
  mutate("ln_N>r" = log(`cumul_freq_N>r`))

# Fractal plot
ggplot(tauroa_obsidian_fractal, aes(x = ln_r, y = `ln_N>r`)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  stat_regline_equation(label.y = 3, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 2.5, aes(label = ..rr.label..)) +
  labs(x = "ln(r) (mm)", y = "ln(N(>r))") +
  theme_bw()


# ~ Te Mataku Basalt

# Table with fractal variables
temataku_basalt_fractal <- temataku_data %>%
  # Filter material and size > 20 mm
  filter(MAT.TYPE == "BASALT") %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>%
  # Filter artefacts of interest
  filter(FCLASS == "COMPFLAKE" | FCLASS == "PROXFLAKE" | FCLASS == "COMPSPLIT" | FCLASS == "PROXSPLIT" | 
         FCLASS == "COMPTOOL" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLITTOOL") %>% 
  # Create size class intervals
  mutate(SIZE_INT = case_when(MAXDIM >= 20 & MAXDIM < 30 ~ "20-30",
                              MAXDIM >= 30 & MAXDIM < 40 ~ "30-40",
                              MAXDIM >= 40 & MAXDIM < 50 ~ "40-50",
                              MAXDIM >= 50 & MAXDIM < 60 ~ "50-60",
                              MAXDIM >= 60 & MAXDIM < 70 ~ "60-70",
                              MAXDIM >= 70 & MAXDIM < 80 ~ "70-80",
                              MAXDIM >= 80 ~ ">80")) %>% 
  mutate(SIZE_INT = factor(SIZE_INT, levels = c("20-30", "30-40", "40-50", "50-60", "60-70", "70-80", ">80"))) %>% 
  # Calculate frequencies for each size class
  dplyr::select(SIZE_INT) %>% 
  group_by(SIZE_INT) %>% 
  summarise(count = n()) %>% 
  # Add lower bound column
  cbind(lower_bound_r = c(20, 30, 40, 50, 60, 70, 80)) %>% 
  # Calculate cumulative frequency
  mutate("cumul_freq_N>r" = rev(cumsum(rev(count)))) %>% 
  # Calculate ln(r)
  mutate(ln_r = log(lower_bound_r)) %>% 
  # Calculate ln(N(>r))
  mutate("ln_N>r" = log(`cumul_freq_N>r`))

# Fractal plot
ggplot(temataku_basalt_fractal, aes(x = ln_r, y = `ln_N>r`)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  stat_regline_equation(label.y = 3, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 2.5, aes(label = ..rr.label..)) +
  labs(x = "ln(r) (mm)", y = "ln(N(>r))") +
  theme_bw()


# ~ Te Mataku Chert 

# Table with fractal variables
temataku_chert_fractal <- temataku_data %>%
  # Filter material and size > 20 mm
  filter(MAT.TYPE == "CHERT") %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>%
  # Filter artefacts of interest
  filter(FCLASS == "COMPFLAKE" | FCLASS == "PROXFLAKE" | FCLASS == "COMPSPLIT" | FCLASS == "PROXSPLIT" | 
         FCLASS == "COMPTOOL" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLITTOOL") %>% 
  # Create size class intervals
  mutate(SIZE_INT = case_when(MAXDIM >= 20 & MAXDIM < 30 ~ "20-30",
                              MAXDIM >= 30 & MAXDIM < 40 ~ "30-40",
                              MAXDIM >= 40 & MAXDIM < 50 ~ "40-50",
                              MAXDIM >= 50 & MAXDIM < 60 ~ "50-60",
                              MAXDIM >= 60 & MAXDIM < 70 ~ "60-70",
                              MAXDIM >= 70 & MAXDIM < 80 ~ "70-80",
                              MAXDIM >= 80 ~ ">80")) %>% 
  mutate(SIZE_INT = factor(SIZE_INT, levels = c("20-30", "30-40", "40-50", "50-60", "60-70", "70-80", ">80"))) %>% 
  # Calculate frequencies for each size class
  dplyr::select(SIZE_INT) %>% 
  group_by(SIZE_INT) %>% 
  summarise(count = n()) %>% 
  # Add lower bound column
  cbind(lower_bound_r = c(20, 30, 40, 50, 60, 70, 80)) %>% 
  # Calculate cumulative frequency
  mutate("cumul_freq_N>r" = rev(cumsum(rev(count)))) %>% 
  # Calculate ln(r)
  mutate(ln_r = log(lower_bound_r)) %>% 
  # Calculate ln(N(>r))
  mutate("ln_N>r" = log(`cumul_freq_N>r`))

# Fractal plot
ggplot(temataku_chert_fractal, aes(x = ln_r, y = `ln_N>r`)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  stat_regline_equation(label.y = 3, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 2.5, aes(label = ..rr.label..)) +
  labs(x = "ln(r) (mm)", y = "ln(N(>r))") +
  theme_bw()


# ~ Te Mataku Obsidian 

# Table with fractal variables
temataku_obsidian_fractal <- temataku_data %>%
  # Filter material and size > 20 mm
  filter(MAT.TYPE == "OBSIDIAN") %>% 
  mutate(MAXDIM = apply(.[,38:40], MARGIN = 1, FUN = max)) %>%
  filter(MAXDIM >= 20) %>%
  # Filter artefacts of interest
  filter(FCLASS == "COMPFLAKE" | FCLASS == "PROXFLAKE" | FCLASS == "COMPSPLIT" | FCLASS == "PROXSPLIT" | 
           FCLASS == "COMPTOOL" | FCLASS == "PROXTOOL" | FCLASS == "COMPSPLITTOOL" | FCLASS == "PROXSPLITTOOL") %>% 
  # Create size class intervals
  mutate(SIZE_INT = case_when(MAXDIM >= 20 & MAXDIM < 25 ~ "20-25",
                              MAXDIM >= 25 & MAXDIM < 30 ~ "25-30",
                              MAXDIM >= 30 & MAXDIM < 35 ~ "30-35",
                              MAXDIM >= 35 & MAXDIM < 40 ~ "35-40",
                              MAXDIM >= 40 & MAXDIM < 45 ~ "40-45",
                              MAXDIM >= 45 & MAXDIM < 55 ~ "45-55",
                              MAXDIM >= 55 ~ ">55")) %>% 
  mutate(SIZE_INT = factor(SIZE_INT, levels = c("20-25", "25-30", "30-35", "35-40", "40-45", "45-55", ">55"))) %>% 
  # Calculate frequencies for each size class
  dplyr::select(SIZE_INT) %>% 
  group_by(SIZE_INT) %>% 
  summarise(count = n()) %>% 
  # Add lower bound column
  cbind(lower_bound_r = c(20, 25, 30, 35, 40, 45, 55)) %>% 
  # Calculate cumulative frequency
  mutate("cumul_freq_N>r" = rev(cumsum(rev(count)))) %>% 
  # Calculate ln(r)
  mutate(ln_r = log(lower_bound_r)) %>% 
  # Calculate ln(N(>r))
  mutate("ln_N>r" = log(`cumul_freq_N>r`))

# Fractal plot
ggplot(temataku_obsidian_fractal, aes(x = ln_r, y = `ln_N>r`)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  stat_regline_equation(label.y = 3, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 2.5, aes(label = ..rr.label..)) +
  labs(x = "ln(r) (mm)", y = "ln(N(>r))") +
  theme_bw()


# ~ Save fractal tables as csv files
colnames(tauroa_basalt_fractal) <- c("Size interval (mm)", "Frequency", "Lower bound (r)", "Cumulative frequency (N>r)", "ln(r)", "ln(N>r)")
write_csv(tauroa_basalt_fractal, "tauroa_basalt_fractal_table.csv")

colnames(tauroa_chert_fractal) <- c("Size interval (mm)", "Frequency", "Lower bound (r)", "Cumulative frequency (N>r)", "ln(r)", "ln(N>r)")
write_csv(tauroa_chert_fractal, "tauroa_chert_fractal_table.csv")

colnames(tauroa_obsidian_fractal) <- c("Size interval (mm)", "Frequency", "Lower bound (r)", "Cumulative frequency (N>r)", "ln(r)", "ln(N>r)")
write_csv(tauroa_obsidian_fractal, "tauroa_obsidian_fractal_table.csv")

colnames(temataku_basalt_fractal) <- c("Size interval (mm)", "Frequency", "Lower bound (r)", "Cumulative frequency (N>r)", "ln(r)", "ln(N>r)")
write_csv(temataku_basalt_fractal, "temataku_basalt_fractal_table.csv")

colnames(temataku_chert_fractal) <- c("Size interval (mm)", "Frequency", "Lower bound (r)", "Cumulative frequency (N>r)", "ln(r)", "ln(N>r)")
write_csv(temataku_chert_fractal, "temataku_chert_fractal_table.csv")

colnames(temataku_obsidian_fractal) <- c("Size interval (mm)", "Frequency", "Lower bound (r)", "Cumulative frequency (N>r)", "ln(r)", "ln(N>r)")
write_csv(temataku_obsidian_fractal, "temataku_obsidian_fractal_table.csv")


# ~ Plot fractal distributions for all materials together

# Combine all data into a single data frame for plotting
all_fractals <- rbind(cbind(assem = "Tauroa BASALT", tauroa_basalt_fractal),
                      cbind(assem = "Tauroa CHERT", tauroa_chert_fractal),
                      cbind(assem = "Tauroa OBSIDIAN", tauroa_obsidian_fractal),
                      cbind(assem = "Te Mataku BASALT", temataku_basalt_fractal),
                      cbind(assem = "Te Mataku CHERT", temataku_chert_fractal),
                      cbind(assem = "Te Mataku OBSIDIAN", temataku_obsidian_fractal))

# Fractal plot
ggplot(all_fractals, aes(x = `ln(r)`, y = `ln(N>r)`, shape = assem, color = assem)) +
  geom_point(size = 2) +
  geom_smooth(method = lm, se = FALSE) +
  stat_regline_equation(label.x = 3, label.y = 2, aes(label = ..eq.label..)) +
  stat_regline_equation(label.x = 3, label.y = 1.5, aes(label = ..rr.label..)) +
  facet_wrap(vars(assem), nrow = 2, ncol = 3) +
  labs(x = "ln(r) (mm)", y = "ln(N(>r))") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_shape_manual(name = "Assemblage", values = c(15, 16, 17, 15, 16, 17)) +
  scale_color_manual(name = "Assemblage", values = c("blue", "blue", "blue", "orange", "orange", "orange"))

ggsave("all_fractals_plot.tiff", units="mm", width=200, height=140, dpi=600, compression = 'lzw')