# Figure 2, panel a: carbon_nitrogen plot, color by sex, comvex hulls by sex

library(dplyr)
library(ggplot2)
library(ggrepel)

iso_dat <- read.csv(file.choose())

iso_clean <- iso_dat %>%
  mutate(
    Sex = as.character(Sex),           # convert factor → character
    Sex = na_if(Sex, ""),             # remove empty strings
    Sex = ifelse(Sex %in% c("M","F"), Sex, NA),   # keep only M/F
    Sex = factor(Sex, levels = c("M","F")),       # convert back to factor
    Sample = factor(Sample, levels = c("Human","Herbivore","S. domesticus"))
  ) %>%
  filter(!is.na(d13C), !is.na(d15N))


# carbon_nitrogen plot, color by sex, comvex hulls by sex

hulls_sex <- iso_clean %>%
  filter(Sample == "Human") %>%        
  group_by(Sex) %>%
  slice(chull(d13C, d15N)) %>%
  ungroup()

ggplot() +
  geom_polygon(
    data = hulls_sex,
    aes(
      x = d13C,
      y = d15N,
      group = Sex,
      fill  = Sex
    ),
    alpha  = 0.2,     
    colour = NA
  ) +
  geom_point(
    data = iso_clean,
    aes(
      x = d13C,
      y = d15N,
      shape = Sample,
      fill  = Sex
    ),
    size   = 3,
    alpha  = 0.6,
    colour = "black",
    stroke = 0.3
  ) +
  scale_shape_manual(
    name   = "Sample",
    values = c(
      "Human"         = 21,
      "Herbivore"     = 24,
      "S. domesticus" = 23
    )
  ) +
  scale_fill_manual(
    name   = "Sex",
    values = c(
      "M" = "#1f78b4",
      "F" = "#e31a1c"
    )
  ) +
  labs(
    x = expression(delta^13*C~"(‰)"),
    y = expression(delta^15*N~"(‰)")
  ) +
  theme_classic()


















# Figure 2, panel b: Novilara vs comparative dataset


library(dplyr)
library(ggrepel)
library(ggplot2)

## Your data frame: iso_dat
## Columns: Context, Chronology, d13C, d15N, Reference, SampleType

# 1. Summarise: mean & SD for each Context × Chronology × Reference × SampleType


iso_dat <- read.csv("file.choose()")

iso_summary <- iso_dat %>%
  group_by(Context, Chronology, Reference, SampleType) %>%
  summarise(
    mean_d13C = mean(d13C, na.rm = TRUE),
    sd_d13C   = sd(d13C,   na.rm = TRUE),
    mean_d15N = mean(d15N, na.rm = TRUE),
    sd_d15N   = sd(d15N,   na.rm = TRUE),
    n         = n(),
    .groups   = "drop"
  )

# 2. Label for legend: Context + Chronology + Reference
iso_summary <- iso_summary %>%
  mutate(
    ContextLabel = paste0(Context, " (", Chronology, "; ", Reference, ")")
  )

# 3. Single plot: d13C (x) vs d15N (y) with SD error bars



iso_summary <- iso_dat %>%
  group_by(Context, Chronology, Reference, SampleType) %>%
  summarise(
    mean_d13C = mean(d13C, na.rm = TRUE),
    sd_d13C   = sd(d13C,   na.rm = TRUE),
    mean_d15N = mean(d15N, na.rm = TRUE),
    sd_d15N   = sd(d15N,   na.rm = TRUE),
    n         = n(),
    .groups   = "drop"
  ) %>%
  mutate(
    LegendLabel = paste0(Context, " (", Chronology, "; ", Reference, ")"),
    is_novilara = ifelse(Context == "Novilara", "Novilara", "Other")
  )



ggplot(
  iso_summary,
  aes(
    x = mean_d13C,
    y = mean_d15N,
    color = LegendLabel,   # now exists
    shape = SampleType,
    alpha = is_novilara
  )
) +
  geom_errorbarh(
    aes(xmin = mean_d13C - sd_d13C,
        xmax = mean_d13C + sd_d13C),
    height = 0,
    linetype = "dotted"
  ) +
  geom_errorbar(
    aes(ymin = mean_d15N - sd_d15N,
        ymax = mean_d15N + sd_d15N),
    width = 0,
    linetype = "dotted"
  ) +
  geom_point(size = 3) +
  geom_label_repel(
    data = subset(iso_summary, SampleType == "Human"),
    aes(label = Context),
    size = 3,
    label.size = 0.1,
    max.overlaps = 20,
    show.legend = FALSE
  ) +
  scale_alpha_manual(values = c("Other" = 0.35, "Novilara" = 1)) +
  labs(
    x = expression(delta^13*C~"(‰, mean ± SD)"),
    y = expression(delta^15*N~"(‰, mean ± SD)"),
    color = "Context (chronology; reference)",
    shape = "Sample type",
    alpha = NULL
  ) +
  theme_classic()