SI1.Rmd

---
title: 'Revisiting the concept of the “Neolithic Founder Crops” in southwest Asia'
subtitle: "Supplementary Information 1 – Quantitative Analysis"
author:
  - "Joe Roe <<joe@joeroe.io>>"
  - "Amaia Arranz-Otaegui"
date: "`r Sys.Date()`"
output: 
  bookdown::html_document2:
    toc: true
    toc_depth: 4
bibliography: references.bib
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE
)
```

```{r deps}
library("SWAsiaNeolithicFounderCrops") # Run devtools::build() first
library("cowplot")
library("dplyr")
library("forcats")
library("ggplot2")
library("ggrepel")
library("ggspatial")
library("glue")
library("gt")
library("here")
library("khroma")
library("patchwork")
library("purrr")
library("readr")
library("sf")
library("stringr")
library("tidyr")
```

# Introduction

This document describes the quantitative analysis supporting our paper 'Revisiting the concept of the “Neolithic Founder Crops” in southwest Asia' in *Vegetation History and Archaeobotany*.
It is generated from the compendium of R code and data available at <https://doi.org/10.5281/zenodo.5911219> or as a git repository at <https://github.com/joeroe/SWAsiaNeolithicFounderCrops>.

# Methods & Materials

Please refer to accompanying publication for an introduction and justification for the data and quantitative approaches selected.
Here, we explain in more detail the steps necessary to reproduce the figures, tables, and other results presented in the main text, and present some supplementary results that were not included due to length constraints.

## Study region and period

```{r data-spatial, fig.cap='The study region', results='hide'}
proj <- 22770
unproj <- 4326

sw_asia <- st_read(here("analysis/data/raw_data/regions/swasia.geojson"))

land <- st_read(here("analysis/data/raw_data/ne/ne_50m_land.shp"))
rivers <- st_read(here("analysis/data/raw_data/ne/ne_50m_rivers_lake_centerlines_scale_rank.shp"))
lakes <- st_read(here("analysis/data/raw_data/ne/ne_50m_lakes.shp"))
highland <- st_read(here("analysis/data/raw_data/srtm/srtm30_sw_asia_1000m.geojson"))
highland <- filter(highland, above_1000m == 1)

base_map <- ggplot() +
  shadow_spatial(st_bbox(sw_asia), fill = NA) +
  annotation_spatial(land, fill = "white") +
  annotation_spatial(highland, fill = "#f0f0f0", colour = NA) +
  annotation_spatial(rivers, colour = "#555555", size = 0.5) +
  annotation_spatial(lakes, fill = "lightblue") +
  coord_sf(expand = FALSE, crs = proj)

base_map
```

We aimed to collate archaeobotanical data from Southwest Asia (fig. \@ref(fig:data-spatial)) from the Neolithic period, c. 11,700 – 6500 BP (table \@ref(tab:periods)).
To mitigate against an "edge effect", we also included samples in the source databases dated to the preceding Late Epipalaeolithic (c. 15,000–11,700 BP) and succeeding Chalcolithic (c. 6500–5000 BP).

```{r periods}
periods <- tribble(
  ~period,             ~agriculture,                          ~start_bp, ~end_bp,
  "Late Epipal.",      "Foraging",                            15000,     11700,
  "PPNA",              "Pre-domestication cultivation",       11700,     10700,
  "EPPNB",             "Cultivation of domesticated species", 10700,     10200,
  "MPPNB",             "Cultivation of domesticated species", 10200,     9500,
  "LPPNB/C",           "Agriculture",                         9500,      8500,
  "Pottery Neolithic", "Agriculture",                         8500,      6500,
  "Chalcolithic",      "Agriculture",                         6500,      5000,
)

neolithic <- c("PPNA", "EPPNB", "MPPNB", "LPPNB/C", "Pottery Neolithic")
neolithic_start <- max(filter(periods, period %in% neolithic) %>% select(start_bp))
neolithic_end <- min(filter(periods, period %in% neolithic) %>% select(end_bp))

periods %>% 
  gt(caption = "Cultural periods used in the analysis") %>% 
  cols_label(period = "Period", agriculture = "Subsistence",
             start_bp = "Start (cal BP)", end_bp = "End (cal BP")
```

## Archaeobotanical data (Southwest Asia)

We collated data from three large-scale archaeobotanical databases: ORIGINS [@Wallace2018], and COMPAG [@Lucas2018, @Fuller2018; incorporating @Colledge2004, @Shennan2007], and ADEMNES [@Riehl2005].
To clean and make these datasets comparable, we:

* Standardised site names across databases (thesaurus available at <https://github.com/joeroe/swapdata>);
* Standardised taxonomic names across databases and applied various additional classifications described below;
* When sites were in multiple databases, preferred records from the most detailed/recent one (i.e. ORIGINS > COMPAG > ADEMNES);
* ORIGINS is recorded by sample, but COMPAG and ADEMNES by phase, so aggregated ORIGINS by phase for comparability – and because individual samples are more likely to be 'noisy';
* Calculated phase-level proportions for COMPAG, which only includes absolute frequencies;
* Excluded records from the source databases that:
  * Were not taxonomically determined to at least the level of genus;
  * Did not actually contain quantitative information (i.e. `n` is blank or missing);
  * Described wood remains.

Finally we filtered this data to include only samples from our region (Southwest Asia) and period (11.7–6.5 ka cal BP) of interest,
with a buffer of ±2,000 years to reduce "edge effects" in our time series analyses.

```{r data-flora}
start_bp <- first(periods$start_bp)
end_bp <- last(periods$end_bp)
flora <- collate_flora(
  origins_path = here("analysis/data/raw_data/origins/"),
  ademnes_path = here("analysis/data/raw_data/ademnes"),
  compag_path = here("analysis/data/raw_data/compag2018"),
  taxa_path = here("analysis/data/raw_data/founder_crops.tsv"),
  region = sw_asia,
  start_bp = start_bp,
  end_bp = end_bp
)

# Final cleaning
# Drop records that are insufficiently taxonomically determined (taxon = NA) or
# which don't actually have quantitative information, to make sure they don't
# throw off our ubiquity measures. 
# We don't update prop, so they are still accounted for in abundance measures.
flora %>% 
  drop_na(taxon) %>% 
  drop_na(prop) %>% 
  filter(prop != 0)  %>% 
  filter(!wood | is.na(wood)) %>% 
  select(-wood) ->
  flora

# Summary totals
flora %>% 
  distinct(site_name) %>% 
  nrow() ->
  n_site

flora %>% 
  distinct(site_name, phase_code) %>% 
  nrow() ->
  n_assemb
```

The resulting collated dataset (`analysis/data/derived_data/swasia_neolithic_flora.tsv`) includes archaeobotanical assemblages representing `r n_assemb` distinct phases from `r n_site` sites.

## Archaeobotanical data (Europe)

```{r data-flora-europe}
flora_europe <- read_colledge2004(here("analysis/data/raw_data/colledge2004/"),
                                  europe_only = TRUE)

n_assemb_europe <- length(unique(flora_europe$phase))
n_site_europe <- length(unique(flora_europe$site))
```

We also used Colledge et al.'s [-@Colledge2004] compilation of European data [see also @Shennan2007].
Since this data did not need to be formally incorporated into the rest of the analyses, we simply imported Colledge et al.'s dataset as is, excluding non-European countries. 
It includes `r n_assemb_europe` assemblages from `r n_site_europe` Neolithic sites across Europe, but only presence data.
For completeness, Cyprus was included in both the European and Southwest Asian datasets.

We used this data to quantitatively assess early evidence for translocation of certain crops from Southwest Asia to Europe, reported in table 3 of the main text.

## Chronological data

All three databases include rough absolute date ranges (in the vast majority of cases summarised the available radiocarbon dates) for each sample.
Using this data, each phase was classified into one of the cultural periods defined above using its mid-point date (table \@ref(tab:sample-periods)).

```{r sample-periods}
flora <- mutate(
  flora,
  age_mid = age_end + ((age_start - age_end) / 2),
  period = cut(-age_mid, 
               breaks = -c(periods$start_bp, 0), 
               labels = periods$period,
               ordered_result = TRUE)
)

# MS Table 2
table2 <- flora %>% 
  distinct(period, site_name, phase) %>% 
  add_count(period, name ="n_phases") %>% 
  distinct(period, site_name, n_phases) %>% 
  add_count(period, name = "n_sites") %>% 
  distinct(period, n_phases, n_sites) %>% 
  select(period, n_sites, n_phases) %>% 
  arrange(period) %>%
  left_join(periods, by = "period") %>%
  transmute(
    period,
    ka_cal_bp = paste0(start_bp / 1000, "–", end_bp / 1000),
    agriculture,
    n_sites,
    n_phases
  ) %>% 
  gt(caption = "Sites and assemblages per period") %>% 
  cols_label(
    period = "Period", 
    ka_cal_bp = "ka cal BP",
    agriculture = "Subsistence",
    n_sites = "N sites", 
    n_phases = "N assemblages") %>% 
  tab_style(cell_text(weight = "bold"), 
            cells_body(columns = period, rows = period %in% neolithic))
table2
gtsave(table2, here("analysis/figures/table2.html"))

# Save key statistics for later
flora %>% 
  distinct(site_name, phase_code, period) %>% 
  filter(period %in% neolithic) %>% 
  nrow() ->
  n_assemb_neolithic

flora %>% 
  filter(period %in% neolithic) %>% 
  distinct(site_name) %>% 
  nrow() ->
  n_site_neolithic
```

For a finer-grained chronology, we also sliced the assemblages into century bins.
An assemblage was considered to belong to a bin (e.g. 5000–5099 BP) if any part of its absolute date range falls within that bin – this results in the duplication of data across bins, but better reflects the inherent imprecision of radiocarbon dating.

```{r sample-centuries}
centuries <- seq(from = first(periods$start_bp),
                 to = last(periods$end_bp),
                 by = -100)
names(centuries) <- centuries

flora %>% 
  mutate(century = map_dfc(centuries, ~. < age_start & . > age_end)) %>% 
  unpack(century) %>% 
  pivot_longer((ncol(flora)+1):ncol(.), names_to = "century",
               names_transform = list(century = as.integer)) %>% 
  filter(value) %>% 
  select(-value) %>% 
  mutate(period = cut(-century, 
               breaks = -c(periods$start_bp, last(periods$end_bp)), 
               labels = periods$period,
               ordered_result = TRUE)) ->
  flora_ts
```

## Sample coverage

We are now in a position to inspect the geographical (fig. \@ref(fig:sample-map)) and temporal \@ref(fig:sample-chronology) coverage of the data.
Unsurprisingly, the distribution is uneven. 
Notably, we have few Epipalaeolithic sites, but a very large number of Chalcolithic sites.
The regional coverage reflects broader trends in the research history of the region, with a large body of evidence from the Southern Levant, especially for earlier periods, and the rest of the region more patchily covered.
Nevertheless, we consider that we have a sufficient baseline number of samples for each time slice in our analysis—the lowest is 11 assemblages, for the EPPNB—especially for the key Neolithic periods.

```{r sample-summary}
flora %>% 
  group_by(site_name, phase_code) %>% 
  summarise(across(c(source, latitude, longitude, phase, age_start, age_end,
                     period), first),
            .groups = "drop") ->
  flora_phases
```


```{r sample-map, fig.cap='Geographical distribution of sampled archaeobotanical assemblages'}
named_sites <- c(
  "Ahihud", # Not in flora data
  "Asikli Höyük", # Not in flora data
  "Cafer Höyük",
  "Chogha Golan",
  "Demirköy", # Not in flora data
  "Dhra'", # Not in flora data
  "Erbaba",
  "Gilgal I",
  "Gritille", # Not in flora data
  "Gusir Höyük", # Not in flora data
  "Hallan Çemi", # Not in flora data
  "Hacilar",
  "Höyücek",
  "Iraq ed-Dubb",
  "Jerf el Ahmar",
  "Kosak Shamali",
  "Körtik Tepe",
  "M'lefaat",
  "Nahal Oren",
  "Netiv Hagdud",
  "Nevali Çori", # Not in flora data
  "Sabi Abyad II",
  "Sheikh-e Abad",
  "Tell Aswad",
  "Tell Ain el-Kerkh",
  "Tell es-Sultan",
  "Tell Qarassa North", # Not in flora data
  "Tell Qarassa", # Not in flora data
  "Tell Qaramel", # Not in flora data
  "Tell Ramad",
  "Yiftah'el",
  "Zahrat adh-Dhra 2",
  "Çatalhöyük",
  "Çayönü"
)

flora_phases |> 
  distinct(site_name, period, latitude, longitude) |> 
  mutate(site_name = recode(site_name, `Hacılar` = "Hacilar")) |> # for EPS compatibility
  mutate(label = site_name %in% named_sites) ->
  p1_data

p1 <- base_map + 
  geom_spatial_point(data = p1_data, 
                     mapping = aes(longitude, latitude),
                     size = 0.5,
                     crs = unproj) +
  labs(x = NULL, y = NULL)
  
p1_top <- p1 + 
  geom_spatial_text_repel(
    data = p1_data |> 
      filter(label) |> 
      distinct(site_name, longitude, latitude),
    mapping = aes(longitude, latitude, label = site_name),
    size = 2,
    min.segment.length = 0
  ) +
  theme(text = element_text(size = 9))
p1_bottom <- p1 + 
  facet_wrap(vars(period), nrow = 1) + 
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        text = element_text(size = 9))

p1_top + inset_element(p1_bottom, 0, 0, 1, 0.25, align_to = "panel")

ggsave(here("analysis/figures/fig1.eps"), width = 174, height = 145,
       units = "mm")
ggsave(here("analysis/figures/fig1.png"), width = 174, height = 145,
       units = "mm", dpi = 300)
```


```{r sample-chronology, fig.cap='Distribution of sampled archaeobotanical assemblages by cultural period. Black lines indicate our period of interest (the Neolithic).'}
p2_top <- flora_phases %>% 
  count(period) %>% 
  left_join(periods) %>% 
  mutate(period = ordered(period, levels = periods$period)) %>% 
  ggplot(aes(xmin = start_bp, xmax = end_bp, ymin = 0, ymax = n,
             fill = period)) +
  geom_rect() +
  geom_vline(xintercept = c(neolithic_start, neolithic_end)) +
  scale_x_reverse(limits = c(15000, 5000),
                  breaks = seq(15000, 5000, -1000)) +
  scale_y_continuous(limits = c(0, 150)) +
  khroma::scale_fill_muted() +
  labs(x = NULL, y = "Samples by period", fill = "Period") +
  theme_minimal() +
  theme(legend.position = "right", 
        axis.text.x = element_blank(), axis.ticks.x = element_blank())

p2_bottom <- flora_ts %>% 
  group_by(site_name, phase_code, century) %>% 
  summarise(period = first(period), .groups = "drop") %>% 
  ggplot(aes(century, fill = period)) +
  geom_bar(width = 100, na.rm = TRUE) +
  geom_vline(xintercept = c(neolithic_start, neolithic_end)) +
  scale_x_reverse(limits = c(15000, 5000), 
                  breaks = seq(15000, 5000, -1000)) +
  scale_y_reverse(limits = c(150, 0)) +
  khroma::scale_fill_muted(guide = guide_none()) +
  labs(x = "cal BP", y = "Samples by century") +
  theme_minimal() +
  theme(legend.position = "bottom")

p2_top / p2_bottom
```
<!--
Table \@ref(tab:all-assemblages) lists all the assemblages included in our analysis.
-->

```{r all-assemblages}
# Breaks PDF generation :/
# flora_phases %>%
#   select(
#     site_name,
#     phase,
#     source,
#     age_start,
#     age_end,
#     period
#   ) %>%
#   mutate(
#     age_start = round(age_start / 1000, digits = 2),
#     age_end = round(age_end / 1000, digits = 2)
#   ) %>%
#   group_by(site_name) %>%
#   gt(caption = "Assemblages included in the analysis") %>%
#   cols_label(
#     phase = "Site/assemblage",
#     source = "Source",
#     age_start = "Reported age ka cal BP",
#     period = "Assigned period"
#   ) %>%
#   cols_merge(c(age_start, age_end), pattern = "{1}&ndash;{2}") %>%
#   cols_align("right", period) %>%
#   fmt_missing(everything(), missing_text = "—") %>%
#   tab_style(cell_text(weight = "bold"), cells_row_groups())
```

# Results

## Ubiquity of all taxa in the Neolithic

Table \@ref(tab:all-ubiquity) summarises the occurrence of each standardised taxon across Neolithic assemblages, broken down by plant category.

```{r all-ubiquity}
flora %>% 
  drop_na(category, taxon) %>% 
  filter(period %in% neolithic) %>% 
  group_by(category, taxon, founder_crop, phase_code) %>% 
  summarise(prop = sum(prop, na.rm = TRUE), .groups = "drop_last") %>% 
  summarise(
    n_present = sum(prop > 0),
    n_gthalf = sum(prop > 0.5),
    n_gtqrtr = sum(prop > 0.25),
    .groups = "drop"
  ) %>% 
  mutate(
    pct_present = n_present / n_assemb_neolithic,
    pct_gthalf = n_gthalf / n_assemb_neolithic,
    pct_gtqrtr = n_gtqrtr / n_assemb_neolithic
  ) %>% 
  group_by(category) %>% 
  arrange(-n_present, .by_group = TRUE) %>% 
  gt(caption = "Ubiquity of all taxa across Neolithic assemblages. Founder crops are highlighted in bold.") %>% 
  cols_hide(c(founder_crop)) %>% 
  fmt_percent(c(pct_present, pct_gthalf, pct_gtqrtr)) %>% 
  cols_merge_n_pct(n_present, pct_present) %>% 
  cols_merge_n_pct(n_gthalf, pct_gthalf) %>% 
  cols_merge_n_pct(n_gtqrtr, pct_gtqrtr) %>% 
  cols_label(taxon = "Taxa",
             n_present = "N (%) recorded present",
             n_gthalf = "N (%) making up more than half of assemblage",
             n_gtqrtr = "N (%) making up more than a quarter of assemblage") %>% 
  tab_style(cell_text(weight = "bold"), cells_body(rows = !is.na(founder_crop))) %>% 
  cols_width(taxon ~ pct(40))
```

This summary data underlies figure 5 in the main text.

```{r all-ubiquity-plot, fig.height=12, include=FALSE}
flora %>% 
  drop_na(category, taxon) %>% 
  filter(period %in% neolithic) %>% 
  group_by(category, taxon, founder_crop, phase_code) %>% 
  summarise(prop = sum(prop, na.rm = TRUE), .groups = "drop_last") %>% 
  summarise(
    n_present = sum(prop > 0),
    n_gthalf = sum(prop > 0.5),
    n_gtqrtr = sum(prop > 0.25),
    .groups = "drop"
  ) %>% 
  mutate(
    pct_present = n_present / n_assemb_neolithic,
    pct_gthalf = n_gthalf / n_assemb_neolithic,
    pct_gtqrtr = n_gtqrtr / n_assemb_neolithic
  ) %>% 
  group_by(category) %>% 
  slice_max(n_present, n = 15, with_ties = FALSE) %>% 
  mutate(
    taxon = str_trunc(taxon, 64),
    taxon = fct_reorder(taxon, n_present, .desc = FALSE)
    ) %>% 
  ggplot(aes(taxon, pct_present, fill = !is.na(founder_crop))) +
  facet_wrap(vars(category), ncol = 2, scales = "free_y",
             strip.position = "left") +
  geom_col(colour = "black", width = 0.4) +
  geom_text(aes(y = 0, label = taxon), 
            size = 7 / (14/5), hjust = 0, nudge_x = 0.5) + 
  coord_flip() +
  scale_x_discrete(labels = NULL, expand = expansion(c(0.05, 0), c(0, 0.7))) +
  scale_y_continuous(labels = scales::label_percent(), limits = c(0, 1),
                     expand = expansion(0.025, 0)) +
  scale_fill_manual(values = c("white", "black"), guide = "none") +
  labs(x = NULL, y = "Proportion of assemblages present") +
  theme_half_open(font_size = 9) + 
  theme(axis.ticks.y = element_blank(),
        axis.line.y = element_blank())

ggsave(here("analysis/figures/fig5.eps"), width = 174, height = 174,
       units = "mm")
ggsave(here("analysis/figures/fig5.png"), width = 174, height = 174,
       units = "mm", bg = "white")
```

## Founder crop ubiquity and abundance

### By period

Table \@ref(tab:founder-ubiquity) summarises the ubiquity of the founder crop species in assemblages broken down by period, i.e. for each period, it shows the number of assemblages where *N* distinct founder species are present, for varying values of *N*.
Given the inconsistent levels of identification between assemblages, emmer wheat and einkorn wheat are counted as one crop.
Hence, the maximum possible number of founder crops present is seven rather than eight.

```{r founder-ubiquity}
flora %>% 
  mutate(founder_crop = recode(founder_crop,
                               "emmer wheat" = "wheat",
                               "einkorn wheat" = "wheat",
                               .default = founder_crop)) %>% 
  group_by(period, site_name, phase_code) %>% 
  summarise(n_founder_crop = n_distinct(founder_crop, na.rm = TRUE),
            .groups = "drop") %>% 
  group_by(period) %>% 
  summarise(
    n_assemb = n(),
    n_atleast7 = sum(n_founder_crop >= 7),
    p_atleast7 = n_atleast7 / n_assemb,
    n_atleast6 = sum(n_founder_crop >= 6),
    p_atleast6 = n_atleast6 / n_assemb,
    n_atleast5 = sum(n_founder_crop >= 5),
    p_atleast5 = n_atleast5 / n_assemb,
    n_atleast4 = sum(n_founder_crop >= 4),
    p_atleast4 = n_atleast4 / n_assemb,
    n_atleast3 = sum(n_founder_crop >= 3),
    p_atleast3 = n_atleast3 / n_assemb,
    n_atleast2 = sum(n_founder_crop >= 2),
    p_atleast2 = n_atleast2 / n_assemb,
    n_atleast1 = sum(n_founder_crop >= 1),
    p_atleast1 = n_atleast1 / n_assemb
  ) %>% 
  gt(caption = "Ubiquity of the founder crops by period.") %>% 
  fmt_percent(starts_with("p_")) %>% 
  cols_merge_n_pct(n_atleast7, p_atleast7) %>% 
  cols_merge_n_pct(n_atleast6, p_atleast6) %>% 
  cols_merge_n_pct(n_atleast5, p_atleast5) %>% 
  cols_merge_n_pct(n_atleast4, p_atleast4) %>% 
  cols_merge_n_pct(n_atleast3, p_atleast3) %>% 
  cols_merge_n_pct(n_atleast2, p_atleast2) %>% 
  cols_merge_n_pct(n_atleast1, p_atleast1) %>% 
  tab_spanner("Number of assemblages with N founder crops present:",
              columns = starts_with("n_atleast"), id = "n_atleast") %>% 
  cols_label(
    period = "Period",
    n_assemb = "Total assemblages",
    n_atleast7 = "All",
    n_atleast6 = ">= 6",
    n_atleast5 = ">= 5",
    n_atleast4 = ">= 4",
    n_atleast3 = ">= 3",
    n_atleast2 = ">= 2",
    n_atleast1 = ">= 1"
  )
```

### By century

For a more fine-grained look at the time series, we can use the century bins calculated above.
Figure \@ref(fig:ts-founder-crops) shows the importance of the individual founder crops measured by both cross-assemblage ubiquity and by relative abundance.
We measure ubiquity by calculated the proportion of assemblages dated to each century where one or more, two or more, etc., of the founder crops were reported present;
and abundance by the taking the mean proportion of the assemblage made up by each founder species for each century (hence NB., the y axis of the lower plot will not sum to 100%).
These different measures reveal different trends, as discussed in the main text.

```{r ts-founder-crops, fig.cap='Ubiquity and abundance of the founder crops by century', fig.height=6}
flora_ts %>%
  # Aggregate wheat
  mutate(founder_crop = recode(founder_crop,
                               "einkorn wheat" = "wheat",
                               "emmer wheat" = "wheat",
                               .default = founder_crop)) %>%
  # Count number of assemblages with N founder crops per century
  group_by(century) %>%
  mutate(n_assemb = length(unique(phase_code))) %>% 
  drop_na(founder_crop) %>% 
  group_by(century, n_assemb, phase_code) %>%
  summarise(n_founder_crop = length(unique(founder_crop)), .groups = "drop_last") %>% 
  summarise(
    atleast1 = sum(n_founder_crop >= 1) / n_assemb,
    atleast2 = sum(n_founder_crop >= 2) / n_assemb,
    atleast3 = sum(n_founder_crop >= 3) / n_assemb,
    atleast4 = sum(n_founder_crop >= 4) / n_assemb,
    atleast5 = sum(n_founder_crop >= 5) / n_assemb,
    atleast6 = sum(n_founder_crop >= 6) / n_assemb,
    atleast7 = sum(n_founder_crop >= 7) / n_assemb,
  ) %>% 
  pivot_longer(
    atleast1:atleast7,
    names_to = "founder_crops",
    values_to = "n"
  ) %>% 
  mutate(founder_crops = str_replace(founder_crops, "atleast", ">=")) %>% 
  # Plot
  ggplot() +
  geom_area(aes(century / 1000, n, fill = founder_crops), position = "identity") +
  geom_vline(data = periods, aes(xintercept = start_bp / 1000),
             linetype = "dashed") +
  geom_label(data = periods,
             aes(x = start_bp / 1000, y = Inf, label = period),
             hjust = 0, vjust = 1, label.size = NA, fill = NA, size = 7 * 0.36) +
  scale_x_reverse(breaks = scales::breaks_width(-1), limits = c(11.7, 6.5),
                  expand = expansion(0)) +
  scale_y_continuous(labels = scales::label_percent(accuracy = 1),
                     expand = expansion(mult = c(0, 0.05))) +
  khroma::scale_fill_YlOrBr(discrete = TRUE) +
  labs(x = NULL, y = "Proportion of assemblages", 
       fill = "Founder crops\npresent") +
  theme_half_open(font_size = 9) ->
  p1

# By average proportion
flora_ts %>%
  # Aggregate wheat and drop non-founder crops
  mutate(founder_crop = recode(founder_crop,
                               "einkorn wheat" = "wheat",
                               "emmer wheat" = "wheat",
                               .default = founder_crop)) %>% 
  # Aggregate by founder crops
  group_by(century, site_name, phase_code, founder_crop) %>% 
  summarise(prop = sum(prop, na.rm = TRUE), .groups = "drop") %>% 
  # Add number of assemblages per century
  group_by(century) %>% 
  mutate(n_assemb = length(unique(phase_code))) %>% 
  # Calculate average proportion
  group_by(century, founder_crop) %>%
  summarise(avg_prop = sum(prop) / first(n_assemb)) %>%
  # Refactor categories for plot
  mutate(founder_crop = fct_reorder(founder_crop, avg_prop, .fun = sum, .desc = FALSE)) %>%
  drop_na(founder_crop) %>% 
  # Plot
  ggplot(aes(century / 1000, avg_prop, fill = founder_crop)) +
  geom_area() +
  geom_vline(data = periods, aes(xintercept = start_bp / 1000),
             linetype = "dashed") +
  geom_label(data = periods, 
             aes(x = start_bp / 1000, y = Inf, label = period),
             hjust = 0, vjust = 1, label.size = NA, fill = NA, size = 7 * 0.36) +
  scale_x_reverse(breaks = scales::breaks_width(-1), limits = c(11.7, 6.5),
                  expand = expansion(0)) +
  scale_y_continuous(limits = c(0, 0.5), breaks = scales::breaks_width(0.1),
                     labels = scales::label_percent(accuracy = 1),
                     expand = expansion(mult = c(0, 0.05))) +
  khroma::scale_fill_muted() +
  labs(x = "ka cal BP", y = "Mean proportion of assemblages", fill = NULL) +
  theme_half_open(font_size = 9) ->
  p2

p1 / p2

ggsave(here("analysis/figures/fig4.eps"), width = 174, height = 117,
       units = "mm")
ggsave(here("analysis/figures/fig4.png"), width = 174, height = 117,
       units = "mm")
```

## Abundance of plant categories

To move beyond the concept of "founder crops" vs. "wild plants", we classified each taxonomic record in our dataset on a number of other axes:

* Broad plant category, i.e. grasses, legumes, "wild plants", or fruits/nuts
* Its edibility
* For grasses and legumes, whether it is large/medium- or small-seeded

```{r category-summary}
flora %>% 
  filter(period %in% neolithic) %>% 
  mutate(category = replace_na(category, "Unclassified")) %>% 
  group_by(site_name, phase_code, category) %>%
  summarise(
    n = sum(n, na.rm = TRUE),
    prop = sum(prop, na.rm = TRUE), 
    .groups = "drop"
  ) %>% 
  group_by(category) %>% 
  summarise(
    prop_avg = mean(prop),
    n_present = n(),
    n_gthalf = sum(prop > 0.5),
    n_gtqrtr = sum(prop > 0.25),
    pct_present = n_present / n_assemb_neolithic,
    pct_gthalf = n_gthalf / n_assemb_neolithic,
    pct_gtqrtr = n_gtqrtr / n_assemb_neolithic,
    .groups = "drop"
  ) %>% 
  arrange(desc(prop_avg)) %>% 
  gt(caption = "Abundance and ubiquity of broad plant categories across Neolithic assemblages.") %>% 
  cols_label(
    category = "Category",
    prop_avg = "Average proportion",
    n_present = "Present",
    n_gthalf = "More than ½",
    n_gtqrtr = "More than ¼"
  ) %>% 
  cols_merge_n_pct(n_present, pct_present) %>% 
  cols_merge_n_pct(n_gthalf, pct_gthalf) %>% 
  cols_merge_n_pct(n_gtqrtr, pct_gtqrtr) %>% 
  fmt_percent(starts_with("pct")) %>% 
  fmt_percent(prop_avg) %>% 
  cols_align("left", category)
```

```{r category-summary-byperiod}
flora %>% 
  mutate(category = replace_na(category, "Unclassified")) %>% 
  group_by(period, site_name, phase_code, category) %>%
  summarise(
    n = sum(n, na.rm = TRUE),
    prop = sum(prop, na.rm = TRUE), 
    .groups = "drop"
  ) %>% 
  group_by(period) %>% 
  mutate(n_assemb = length(unique(phase_code))) %>% 
  group_by(period, category) %>% 
  summarise(
    prop_avg = mean(prop),
    n_present = n(),
    n_gthalf = sum(prop > 0.5),
    n_gtqrtr = sum(prop > 0.25),
    pct_present = n_present / first(n_assemb),
    pct_gthalf = n_gthalf / first(n_assemb),
    pct_gtqrtr = n_gtqrtr / first(n_assemb),
    .groups = "drop_last"
  ) %>% 
  arrange(desc(prop_avg), .by_group = TRUE) %>% 
  gt(caption = "Abundance and ubiquity of broad plant categories, by period.") %>% 
  cols_label(
    category = "Category",
    prop_avg = "Average proportion",
    n_present = "Present",
    n_gthalf = "More than ½",
    n_gtqrtr = "More than ¼"
  ) %>% 
  cols_merge_n_pct(n_present, pct_present) %>% 
  cols_merge_n_pct(n_gthalf, pct_gthalf) %>% 
  cols_merge_n_pct(n_gtqrtr, pct_gtqrtr) %>% 
  fmt_percent(starts_with("pct")) %>% 
  fmt_percent(prop_avg) %>% 
  cols_align("left", category)
```

Table \@ref(tab:category-summary) summarises the abundance of the broad plant categories across Neolithic assemblages, \@ref(tab:category-summary-byperiod) the same data broken down by period, and figure \@ref(fig:category-ts) how this changed through time.

```{r category-ts, fig.cap='Abundance of broad plant categories through time'}
flora_ts %>% 
  # Aggregate by category
  group_by(century, site_name, phase_code, category) %>% 
  summarise(prop = sum(prop, na.rm = TRUE), .groups = "drop") %>% 
  # Add total number of assemblages per century
  group_by(century) %>% 
  mutate(n_assemb = length(unique(phase_code))) %>% 
  # Calculate average proportion
  group_by(century, category) %>%
  summarise(avg_prop = sum(prop) / first(n_assemb)) %>%
  # Refactor categories for plot
  mutate(
    category = fct_reorder(category, avg_prop, .fun = sum, .desc = FALSE),
    category = fct_explicit_na(category, "Unclassified")
  ) %>% 
  # Plot
  ggplot(aes(century / 1000, avg_prop, fill = category)) +
  geom_area() +
  geom_vline(data = periods, aes(xintercept = start_bp / 1000)) +
  geom_label(data = periods, 
             aes(x = start_bp / 1000, y = Inf, label = period),
             hjust = 0, vjust = 1, label.size = NA, fill = NA, size = 9 * 0.36) +
  scale_x_reverse(breaks = scales::breaks_width(-1), limits = c(11.7, 6.5),
                  expand = expansion()) +
  scale_y_continuous(breaks = scales::breaks_pretty(5), 
                     labels = scales::percent,
                     limits = c(0, 1),
                     expand = expansion()) +
  khroma::scale_fill_muted() +
  labs(x = "ka cal BP", y = "mean proportion of assemblage", fill = NULL) +
  theme_half_open(font_size = 9) +
  theme(legend.position = "bottom")

ggsave(here("analysis/figures/fig2.eps"), width = 174, height = 87,
       units = "mm")
ggsave(here("analysis/figures/fig2.png"), width = 174, height = 87,
       units = "mm", dpi = 600, bg = "#ffffff")
```

In the following sections, we will analyse trends within these categories in more detail.

### Grasses

```{r common-grasses}
flora %>% 
  filter(category == "Grasses", period %in% neolithic) %>% 
  group_by(taxon) %>% 
  summarise(n_present = n()) %>% 
  arrange(desc(n_present)) %>% 
  slice(1:8) %>% 
  pull(taxon) ->
  top_grasses
```

The eight most ubiquitous grass taxonomic groups, measured by presence across Neolithic assemblages, are: `r paste(top_grasses, collapse = ", ")`.
Figure \@ref(fig:grasses-ts) shows how the abundance of these taxa changed through time.

```{r grasses-ts, fig.cap='Abundance of grass taxa through time'}
p_ts_grasses <- flora_ts %>% 
  # Include only top 8 grasses
  filter(category == "Grasses") %>% 
  mutate(taxon = if_else(taxon %in% top_grasses, taxon, NA_character_)) %>% 
  # Count total number of assemblages per century
  group_by(century) %>% 
  mutate(n_assemb = length(unique(phase_code))) %>% 
  # Calculate average proportion per century
  group_by(century, taxon) %>%
  summarise(avg_prop = sum(prop) / first(n_assemb), .groups = "drop_last") %>%
  # Make gaps explicit
  full_join(expand.grid(century = unique(.$century), taxon = unique(.$taxon)),
            by = c("taxon", "century")) %>%
  mutate(avg_prop = replace_na(avg_prop, 0)) %>% 
  # Refactor & order categories for plotting
  mutate(
    taxon = factor(taxon, levels = rev(top_grasses)),
    taxon = fct_explicit_na(taxon, "Other"),
    taxon = fct_relabel(taxon, str_wrap, width = 40)
  ) %>% 
  # Plot
  ggplot(aes(century / 1000, avg_prop, fill = taxon)) +
  geom_area(na.rm = TRUE) +
  geom_vline(data = periods, aes(xintercept = start_bp / 1000),
             linetype = "dashed") +
  geom_label(data = periods, 
             aes(x = start_bp / 1000, y = Inf, label = period),
             hjust = 0, vjust = 1, label.size = NA, fill = NA, size = 7 * 0.36) +
  scale_x_reverse(breaks = scales::breaks_width(-1), limits = c(11.7, 6.5),
                  expand = expansion()) +
  scale_y_continuous(breaks = scales::breaks_width(0.1), 
                     labels = scales::label_percent(accuracy = 1),
                     expand = expansion(mult = c(0, 0.05))) +
  khroma::scale_fill_muted(guide = guide_legend(ncol = 1)) +
  labs(x = "ka cal BP", y = NULL, fill = "Grasses") +
  theme_half_open(font_size = 9) +
  theme(legend.position = "right", legend.title = element_text(face = "bold"))

p_ts_grasses
```
 
### Pulses

```{r common-pulses}
flora %>% 
  filter(category == "Pulses", period %in% neolithic) %>% 
  group_by(taxon) %>% 
  summarise(n_present = n()) %>% 
  arrange(desc(n_present)) %>% 
  slice(1:8) %>% 
  pull(taxon) ->
  top_pulses
```

The eight most ubiquitous pulse taxonomic groups, measured by presence across Neolithic assemblages, are: `r paste(top_pulses, collapse = ", ")`.
Figure \@ref(fig:pulses-ts) shows how the abundance of these taxa changed through time.

```{r pulses-ts, fig.cap='Abundance of pulse taxa through time'}
p_ts_pulses <- flora_ts %>% 
  # Include only top 8 grasses
  filter(category == "Pulses") %>% 
  mutate(taxon = if_else(taxon %in% top_pulses, taxon, NA_character_)) %>% 
  # Count total number of assemblages per century
  group_by(century) %>% 
  mutate(n_assemb = length(unique(phase_code))) %>% 
  # Calculate average proportion per century
  group_by(century, taxon) %>%
  summarise(avg_prop = sum(prop) / first(n_assemb), .groups = "drop_last") %>%
  # Make gaps explicit
  full_join(expand.grid(century = unique(.$century), taxon = unique(.$taxon)),
            by = c("taxon", "century")) %>%
  mutate(avg_prop = replace_na(avg_prop, 0)) %>% 
  # Refactor & order categories for plotting
  mutate(
    taxon = factor(taxon, levels = rev(top_pulses)),
    taxon = fct_explicit_na(taxon, "Other"),
    taxon = fct_relabel(taxon, str_wrap, width = 40)
  ) %>% 
  # Plot
  ggplot(aes(century / 1000, avg_prop, fill = taxon)) +
  geom_area() +
  geom_vline(data = periods, aes(xintercept = start_bp / 1000),
             linetype = "dashed") +
  geom_label(data = periods, 
             aes(x = start_bp / 1000, y = Inf, label = period),
             hjust = 0, vjust = 1, label.size = NA, fill = NA, size = 7 * 0.36) +
  scale_x_reverse(breaks = scales::breaks_width(-1), limits = c(11.7, 6.5),
                  expand = expansion()) +
  scale_y_continuous(breaks = scales::breaks_width(0.1), 
                     labels = scales::percent,
                     expand = expansion(mult = c(0, 0.05))) +
  khroma::scale_fill_muted(guide = guide_legend(ncol = 1)) +
  labs(x = "ka cal BP", y = NULL, fill = "Pulses") +
  theme_half_open(font_size = 9) +
  theme(legend.position = "right", legend.title = element_text(face = "bold"))

p_ts_pulses
```

### Wild plants

```{r common-wild-plants}
flora %>% 
  filter(category == "Wild plants", period %in% neolithic) %>% 
  group_by(taxon) %>% 
  summarise(n_present = n()) %>% 
  arrange(desc(n_present)) %>% 
  slice(1:8) %>% 
  pull(taxon) ->
  top_wild_plants
```

The eight most ubiquitous 'wild plant' taxonomic groups, measured by presence across Neolithic assemblages, are: `r paste(top_wild_plants, collapse = ", ")`.
Figure \@ref(fig:wild-plants-ts) shows how the abundance of these taxa changed through time.

```{r wild-plants-ts, fig.cap='Abundance of wild plant taxa through time'}
p_ts_wild_plants <- flora_ts %>% 
  # Include only top 8 grasses
  filter(category == "Wild plants") %>% 
  mutate(taxon = if_else(taxon %in% top_wild_plants, taxon, NA_character_)) %>% 
  # Count total number of assemblages per century
  group_by(century) %>% 
  mutate(n_assemb = length(unique(phase_code))) %>% 
  # Calculate average proportion per century
  group_by(century, taxon) %>%
  summarise(avg_prop = sum(prop) / first(n_assemb), .groups = "drop_last") %>%
  # Make gaps explicit
  full_join(expand.grid(century = unique(.$century), taxon = unique(.$taxon)),
            by = c("taxon", "century")) %>%
  mutate(avg_prop = replace_na(avg_prop, 0)) %>% 
  # Refactor & order categories for plotting
  mutate(
    taxon = factor(taxon, levels = rev(top_wild_plants)),
    taxon = fct_explicit_na(taxon, "Other"),
    taxon = fct_relabel(taxon, str_wrap, width = 40)
  ) %>% 
  # Plot
  ggplot(aes(century / 1000, avg_prop, fill = taxon)) +
  geom_area() +
  geom_vline(data = periods, aes(xintercept = start_bp / 1000),
             linetype = "dashed") +
  geom_label(data = periods, 
             aes(x = start_bp / 1000, y = Inf, label = period),
             hjust = 0, vjust = 1, label.size = NA, fill = NA, size = 7 * 0.36) +
  scale_x_reverse(breaks = scales::breaks_width(-1), limits = c(11.7, 6.5),
                  expand = expansion()) +
  scale_y_continuous(breaks = scales::breaks_width(0.1), 
                     labels = scales::percent,
                     expand = expansion(mult = c(0, 0.05))) +
  khroma::scale_fill_muted(guide = guide_legend(ncol = 1)) +
  labs(x = "ka cal BP", y = NULL, fill = "Wild plants") +
  theme_half_open(font_size = 9) +
  theme(legend.position = "right", legend.title = element_text(face = "bold"))

p_ts_wild_plants
```

### Fruits and nuts

```{r common-fruits-nuts}
flora %>% 
  filter(category == "Fruits/nuts", period %in% neolithic) %>% 
  group_by(taxon) %>% 
  summarise(n_present = n()) %>% 
  arrange(desc(n_present)) %>% 
  slice(1:8) %>% 
  pull(taxon) ->
  top_fruits_nuts
```

The eight most ubiquitous fruit/nut taxonomic groups, measured by presence across Neolithic assemblages, are: `r paste(top_fruits_nuts, collapse = ", ")`.
Figure \@ref(fig:fruits-nuts-ts) shows how the abundance of these taxa changed through time.

```{r fruits-nuts-ts, fig.cap='Abundance of fruit/nut taxa through time'}
p_ts_fruits_nuts <- flora_ts %>% 
  # Include only top 8 grasses
  filter(category == "Fruits/nuts") %>% 
  mutate(taxon = if_else(taxon %in% top_fruits_nuts, taxon, NA_character_)) %>% 
  # Count total number of assemblages per century
  group_by(century) %>% 
  mutate(n_assemb = length(unique(phase_code))) %>% 
  # Calculate average proportion per century
  group_by(century, taxon) %>%
  summarise(avg_prop = sum(prop) / first(n_assemb), .groups = "drop_last") %>%
  # Make gaps explicit
  full_join(expand.grid(century = unique(.$century), taxon = unique(.$taxon)),
            by = c("taxon", "century")) %>%
  mutate(avg_prop = replace_na(avg_prop, 0)) %>% 
  # Refactor & order categories for plotting
  mutate(
    taxon = factor(taxon, levels = rev(top_fruits_nuts)),
    taxon = fct_explicit_na(taxon, "Other"),
    taxon = fct_relabel(taxon, str_wrap, width = 40)
  ) %>% 
  # Plot
  ggplot(aes(century / 1000, avg_prop, fill = taxon)) +
  geom_area() +
  geom_vline(data = periods, aes(xintercept = start_bp / 1000),
             linetype = "dashed") +
  geom_label(data = periods, 
             aes(x = start_bp / 1000, y = Inf, label = period),
             hjust = 0, vjust = 1, label.size = NA, fill = NA, size = 7 * 0.36) +
  scale_x_reverse(breaks = scales::breaks_width(-1), limits = c(11.7, 6.5),
                  expand = expansion()) +
  scale_y_continuous(breaks = scales::breaks_width(0.1), 
                     labels = scales::percent,
                     expand = expansion(mult = c(0, 0.05))) +
  khroma::scale_fill_muted(guide = guide_legend(ncol = 1)) +
  labs(x = "ka cal BP", y = NULL, fill = "Fruits & nuts") +
  theme_half_open(font_size = 9) +
  theme(legend.position = "right", legend.title = element_text(face = "bold"))

p_ts_fruits_nuts
```

```{r combined-category-ts, include=FALSE}
p_ts_grasses /
p_ts_pulses /
p_ts_wild_plants /
p_ts_fruits_nuts

ggsave(here("analysis/figures/fig3.eps"), width = 174, height = 234,
       units = "mm")
ggsave(here("analysis/figures/fig3.png"), width = 174, height = 234,
       units = "mm")
```

# References

<!-- auto-generated -->

```{r print-pdf, include=FALSE}
# Unfortunately gt's LaTeX output is still too limited for us to knit to PDF
# directly. Specifically, it lacks support for:
# cross-references to tables
# wrapping text in wide tables
#
# So for the PDF version of this file, we print the HTML output using Chrome:
#
# pagedown::chrome_print(here("analysis/SI1.html"), here("analysis/SI1.pdf"),
#                        options = list(paperWidth = 8.27,
#                                       paperHeight = 11.69))
```

```{css}
@media print {
  a:link:after, a:visited:after {
      content: "";
  }
}
```