benmarwick/mjbnaturepaper
Code and data for analysis of the 2012 and 2015 MJB excavations

Package index
Search the benmarwick/mjbnaturepaper package
Vignettes
Functions
35
Source code
9
Man pages
24
Home
/
GitHub
/
benmarwick/mjbnaturepaper
/

In benmarwick/mjbnaturepaper: Code and data for analysis of the 2012 and 2015 MJB excavations 

knitr::opts_chunk$set(echo = FALSE,
                      warning = FALSE,
                      message = FALSE,
                      quiet = TRUE,
                      progress=FALSE)
# Submission preparation guidelines:
# http://www.nature.com/nature/authors/gta/
# Nature’s standard figure sizes are 89 mm (single column) and 183 mm (double column) and the full depth of the page is 247 mm.
library(tidyverse)
library(ggrepel)

library(knitr)
invisible(purl("supplementary_information.Rmd", output = "temp", quiet=TRUE))
read_chunk("temp")





unlink("temp")

# run osl_ages <- get_osl_ages(cleaned_rotated_points_in_main_excavation_area) from SI.Rmd
osl_ages$Sample <- ifelse(osl_ages$Sample == "NE1C", 
                                   "NE1B", osl_ages$Sample)
# Note that I seem to not to have the latest versions of the OSL ages, so let's get them in from the SI
library(readxl)
osl_ages_from_SI <-  read_excel("../../data/ages/OSL_Data_tables_all_samples_table_from_SI.xlsx")
osl_ages_from_SI_ages  <- 
  osl_ages_from_SI %>% 
  separate(`Age (kyr)‡§`, into = c("osl_age", "other1", "osl_error", "other2"), sep = " ")
osl_ages_from_SI_ages_depths <- left_join(osl_ages_from_SI_ages, osl_ages, by = "Sample")
osl_ages_from_SI_ages_depths$osl_age <- as.numeric(osl_ages_from_SI_ages_depths$osl_age.x)
osl_ages_from_SI_ages_depths$osl_error <- as.numeric(osl_ages_from_SI_ages_depths$osl_error.x)
osl_ages <- osl_ages_from_SI_ages_depths

# move depths for KTL ages
osl_ages$total_station_depth_below_surf <- ifelse(grepl("KTL", osl_ages$Sample), 
                                                  abs(osl_ages$`Depth below surface (m).x`),
                                                  osl_ages$total_station_depth_below_surf)

# make sure they are all -ve values
osl_ages$total_station_depth_below_surf <- -(abs(osl_ages$total_station_depth_below_surf))

OSL ages for lowest artefact band are consistent across the site

Kaifu and Shitaoka write:

"First, the OSL ages for the lowest dense artifact band (called as the ‘zone of initial human occupation’) are different as much as 10,000 years between the two main sections. The ages were ~53,000–55,000 years ago at square B4 (SW–A series), whereas ~57,000–65,000 years ago at a location one meter apart from there (SW–B/C series at square B5). The lower limit of this artifact band at square B4 has been lowered in Extended Data Figure 8c (the lower stippled line), but this must be an error because such a lower limit cannot be supported from the artifact distribution and the lithic artifact refitting analysis illustrated in Extended Data Figures 1, 2 and 6. The observed internal inconsistency indicates that either or both of the OSL ages have considerable errors."

So let's see what is the age of phase 2 in square B4 compared to B5.

Here is EDF 8C, showind OSL samples in the section of B4, the white stipple indicate the lowest dense layer of stone artefacts, the upper boundary of this layer passing between OSL samples SW9A and SW10A:

Here is EDF 1, the stipple lines here do not indicate artefact deposits, but instead indicate stratigraphic boundaries. The line passes between SW10A and SW11A, indicating that the stratigraphic boundary is below the top of the artefact layer. The artefact layer and stratigraphic boundaries do not align exactly. The 1989 pit is a stratigraphic feature, notated here with a stipple line, because of its distinctive texture and colour, in addition to having a high density of artefacts. This misinterpretation of the stipple lines in EDF8c and EDF 1 is the source of the first objection by K & S.

Here is EDF 2, below, showing the locations of the phases of artefacts. Phases 2 is the lowest dense layer, identified on EDF 1 as the 'Lower dense artefact band', and in EDF 8c by the stipple lines. Recalling that the stipple line in EDF 1 is a stratigraphic indicator rather than an indicator of an artefact layer, we see no inconsistency between EDF 2 and the other figures.

Here is EDF 6, showing the artefact refitting data. The depths in panel c are consistent with the depths of phase 2 in the above figures. 

stone_artefacts_only_BC <- 
  stone_artefacts_only %>% 
  filter(square %in% c("B4", "B5", "B6", "C4", "C5", "C6"))


# these are the vertical lines
row_c <- c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6)
row_mids <- row_c/2

library(viridis)
library(grid)
p <- ggplot(stone_artefacts_only,
       aes(Xnew_flipped,
           depth_below_ground_surface,
           colour = find)) +
  geom_point(size = 0.5) +
  scale_y_reverse(limits = c(3,0), 
                  breaks = seq(0, 3, 0.1)) +
  theme_minimal() +
  scale_x_continuous(breaks = row_c,
                     labels = NULL) +
  xlab("") +
  ylab("Depth below \nground surface (m)") +
  scale_color_viridis(discrete=TRUE,
                      "Artefact\ntype") +
  geom_segment(aes(x = -3.6,  
               y = 2.4,
               xend = 0.9,
               yend = 2.4),
               colour = "red") +
    geom_segment(aes(x = -3.6,    # B4
               y = 2.35,
               xend = -0.1,
               yend = 2.35),
               colour = "red") +
  coord_equal() +
  guides(colour = guide_legend(override.aes = list(size = 5))) +
  ggtitle("MJB SW Section with all lithics",
          "showing depths of base of phase 2 in B4 and B5")

row_c = c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6)
nums = paste0("B", 6:1)
row_mids <-  row_c[-length(row_c)] + diff(row_c)/2

for(i in 1:length(row_mids)){
p = p + annotation_custom(grob = textGrob(nums[i], gp=gpar(fontsize=10)),
                         xmin =  row_mids[i],
                         xmax =  row_mids[i],
                         ymin = -8.5,
                         ymax = 2)
}

# Code to override clipping
grid.newpage()
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$clip[gt$layout$name=="panel"] <- "off"
grid.draw(gt)

library(viridis)
library(grid)
p <- ggplot(stone_artefacts_only_BC,
       aes(Xnew_flipped,
           depth_below_ground_surface)) +
  geom_point(size = 0.5, colour = "black") +
  scale_y_reverse(limits = c(3,0), 
                  breaks = seq(0, 3, 0.1)) +
  theme_minimal() +
  scale_x_continuous(breaks = row_c,
                     labels = NULL) +
  xlab("") +
  ylab("Depth below \nground surface (m)") +
  scale_color_viridis(discrete=TRUE,
                      "Artefact\ntype") +
  geom_segment(aes(x = -3.6,  
               y = 2.4,
               xend = 0.9,
               yend = 2.4),
               colour = "red") +
    geom_segment(aes(x = -3.6,    # B4
               y = 2.35,
               xend = -0.1,
               yend = 2.35),
               colour = "red") +
  coord_equal() +
  guides(colour = guide_legend(override.aes = list(size = 5))) +
  ggtitle("MJB SW Section, lithics only from B-C squares",
          "showing depths of base of phase 2 in B4 and B5")

row_c = c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6)
nums = paste0("B", 6:1)
row_mids <-  row_c[-length(row_c)] + diff(row_c)/2

for(i in 1:length(row_mids)){
p = p + annotation_custom(grob = textGrob(nums[i], gp=gpar(fontsize=10)),
                         xmin =  row_mids[i],
                         xmax =  row_mids[i],
                         ymin = -8.5,
                         ymax = 2)
}

# Code to override clipping
grid.newpage()
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$clip[gt$layout$name=="panel"] <- "off"
grid.draw(gt)

# these are the vertical lines
row_c <- c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6)
row_mids <- row_c/2

library(ggrepel)
p1 <- 
ggplot() +
  coord_equal() +
    scale_y_reverse(limits = c(3,0), 
                  breaks = seq(0, 3, 0.1)) +
  scale_x_continuous(breaks = row_c,
                     labels = NULL) +
    geom_point(data = stone_artefacts_only_BC,
           aes(Xnew_flipped,
           depth_below_ground_surface),
           colour = "grey80",
           size = 0.5)  +
    geom_segment(aes(x = -0.5,          # C4
               y = 2.4,
               xend = 0.9,
               yend = 2.4),
               colour = "red") +
    geom_segment(aes(x = -0.5,        # B5
               y = 2.35,
               xend = -0.1,
               yend = 2.35),
               colour = "red") +
  geom_point(data = osl_ages[grepl("SW", osl_ages$Sample), ],
             aes(Xnew_flipped,
                 depth_below_ground_surface ),
             colour = "red") +
  geom_text_repel(data =  osl_ages[grepl("SW", osl_ages$Sample), ],
                  aes(Xnew_flipped,
                      depth_below_ground_surface,
                      label = paste0(osl_age, " (", Sample, ")" )),
                  size = 3) +
  ylab("Depth below \nground surface (m)") +

  xlab("") +
  theme_minimal() +
  ggtitle("MJB SW Section, lithics only from B-C squares",
          "showing OSL sample locations and depths of base of phase 2 in B4 and B5")

row_c = c(2.4, 1.4, 0.4, -0.6)
nums = paste0("B", 6:3)
row_mids <-  row_c[-length(row_c)] + diff(row_c)/2

for(i in 1:length(row_mids)){
p1 = p1 + annotation_custom(grob = textGrob(nums[i], gp=gpar(fontsize=10)),
                         xmin =  row_mids[i],
                         xmax =  row_mids[i],
                         ymin = -8.5,
                         ymax = 2)
}

# Code to override clipping
grid.newpage()
gt1 <- ggplot_gtable(ggplot_build(p1))
gt1$layout$clip[gt1$layout$name=="panel"] <- "off"
grid.draw(gt1)

Let's look at the artefact discard patterns in these squares:

library(readxl)
artefacts_per_litre <- read_excel("data/stone_artefact_data/Artefacts per litre discard BC4-6 by depth.xlsx")

C4 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 1:4]
names(C4) <- C4[1, ]
C4 <- C4[-1, ]
C4 <- map_df(C4, as.numeric)

B4 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 6:11]
names(B4) <- B4[1, ]
B4 <- B4[-1, ]
B4 <- map_df(B4, as.numeric)

B5 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 13:16]
names(B5) <- B5[1, ]
B5 <- B5[-1, ]
B5 <- map_df(B5, as.numeric)

B6 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 18:22]
names(B6) <- B6[1, ]
B6 <- B6[-1, ]
B6 <- map_df(B6, as.numeric)

C5 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 24:28]
names(C5) <- C5[1, ]
C5 <- C5[-1, ]
C5 <- map_df(C5, as.numeric)

C6 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 30:34]
names(C6) <- C6[1, ]
C6 <- C6[-1, ]
C6 <- map_df(C6, as.numeric)

# combine these into one df
the_list <- list("B4"=B4, "B5"=B5, "B6"=B6,
          "C4"=C4, "C5"=C5, "C6"=C6)
artefacts_per_litre_long <- 
bind_rows(the_list,
          .id = "id")

# Depth in two cols...
artefacts_per_litre_long$depth <- 
  with(artefacts_per_litre_long, ifelse(is.na(depth), 
                                        Depth, depth))

vlines <- data.frame(id = rep(names(the_list), 2), 
                     xintercept =c(rep(c(2.15, 
                                        2.15,
                                        2.15), 2),
                                  rep(c(2.6,
                                        2.6, 
                                        2.6), 2)))



ggplot(artefacts_per_litre_long,
       aes(depth, 
           `Artefacts/Litre`)) +
  geom_line() +
  facet_wrap(~id, scales = "free_y"
             ) +
  geom_vline(data = vlines, 
             aes(xintercept = xintercept),
             colour = "red") +
  geom_text(data = vlines, 
           aes(x = xintercept,
               y = -0.5,
               label = xintercept),
           size = 2) +
  theme_minimal() +
  xlab("Depth below surface (m)") +
  ggtitle("MJB artefact density by depth in squares B-C 4-5-6")

ggsave("figures/Artefact density by depth in squares B-C 4-5-6.png")

Show the dated samples from each square, by depth, so we can see the age constraints on each peak of artefacts

# We need to combine OSL and C14, with a col for square, a col for depth, a col for sample ID, and a col for approx age

c14_ages_rugplot_data <- 
c14_ages %>% 
  select(square,
         depth_below_ground_surface,
         Bchron_Median,
         Lab.ID) %>% 
  dplyr::rename(grid_square = square,
                age = Bchron_Median,
         depth = depth_below_ground_surface,
         id = Lab.ID) %>% 
  mutate(age = round(age / 1000, 1))

# For the OSL samples, what square was each sample taken from? We don
# have this in our data so far. We can look EDF 8 from the Nature article. 

# NE, NW, SWA, SWB, SWC,
# E2, C4, B4,  B5,  B5, 

osl_ages_rugplot_data <- 
osl_ages %>% 
  mutate(grid_square = case_when(
  grepl("SW8C", Sample) ~ "B6",
  grepl("NE1B", Sample) ~ "C6",
  grepl("NE", Sample) ~ "E2",
  grepl("NW8B", Sample) ~ "B5",
  grepl("NW9B", Sample) ~ "B5",
  grepl("NW", Sample) ~ "C5",
  grepl("SW.{1,2}A", Sample) ~ "B4",
  grepl("SW.{1,2}B", Sample) ~ "B5",
  grepl("SW.{1,2}C", Sample) ~ "B5",
  grepl("KTL158", Sample) ~ "B4",
  grepl("KTL162", Sample) ~ "B5",
  is.na(Sample) ~ "Unknown"
)) %>% 
  select(grid_square,
         total_station_depth_below_surf,
         osl_age,
         Sample) %>% 
 dplyr::rename(depth = total_station_depth_below_surf,
               age = osl_age,
               id = Sample) %>% 
  mutate(depth = abs(depth),
         age = as.numeric(age)) %>% 
  filter(grid_square != 'Unknown')



# combine the two sets of dates
date_for_plots <- 
bind_rows(c14_ages_rugplot_data,
          osl_ages_rugplot_data)

date_for_plots <- 
  date_for_plots %>% 
  dplyr::rename(id = grid_square,
                sample_id = id) %>% 
  filter(id %in% artefacts_per_litre_long$id )

# find max y-axis value for each sq to plot age points
date_for_plots_max_y <- 
artefacts_per_litre_long %>% 
  group_by(id) %>% 
  summarise(max_artefacts_litre = max(`Artefacts/Litre`, na.rm = TRUE))

# add max height for red lines
vlines <- 
  vlines %>% 
  left_join(date_for_plots_max_y)

date_for_plots <- 
  date_for_plots %>% 
  left_join(date_for_plots_max_y) %>% 
  filter(age > 50) %>% 
  filter(age < 90) %>% 
  mutate( label = paste0(sample_id, ": ", age, " ka"))

# plot with ages

p <- 
ggplot(artefacts_per_litre_long,
       aes(depth, 
           `Artefacts/Litre`)) +
  geom_line() +
  geom_segment(data = vlines, 
             aes(x = xintercept,
                 xend = xintercept,
                 y = 0,
                 yend = max_artefacts_litre + 2),
             colour = "red") +
  geom_text(data = vlines, 
           aes(x = xintercept,
               y = -0.5,
               label = xintercept),
           size = 2) +
  geom_point(data = date_for_plots,
             aes(depth,
                 max_artefacts_litre + 2)) +
  geom_text(data = date_for_plots,
             aes(depth,
                 max_artefacts_litre + 2,
                 label = label),
             size = 2,
            angle = 45, 
            hjust = -0.25,
            vjust = 0) +
  facet_wrap(~id, scales = "free_y"
             ) +
  theme_minimal() +
  xlab("Depth below surface (m)") +
  ggtitle("MJB artefact density by depth in squares B-C 4-5-6")

psvg("figures/Artefact density by depth in squares B-C 4-5-6.svg",
    width = 10, height = 7)
gt = ggplot_gtable(ggplot_build(p))
gt$layout$clip = "off"
grid.draw(gt)
dev.off()

# save, then tidy up in inkscape

OSL and C14 ages are in agreement with each other

Kaifu and Shitaoka write:

"Second, contrary to the authors’ statement, the reported radiocarbon and OSL chronologies are not in agreement with each other. Among the seven radiocarbon ages measured below the –110 cm level (> 20,000 years ago in the OSL time scale), six of them collectively shows its own time trajectory which is younger than the OSL time scale in Extended Data Figure 8g. Clarkson and colleagues suggested ~30–50 cm downward relocation for the two dated charcoals either by artificial or natural processes (Supplementary Information, p.18 in ref. 1). Even if we accept this possibility, the claimed consistency is unsound in the lower stratigraphic section because it is supported by single isolated charcoal (OZT592). Such a piece could be intruded from nearby older deposit in a sandy area like Madjedbebe."

"The data from the squares E and D series show more obvious inconsistency between the two dating methods, although this point was not mentioned by Clarkson and colleagues. Stratigraphic correlations could not be established between the opposite ends of the excavated grids, that is, between the area around square B4 and the area around square E2 (ref. 1). Because of this, the OSL ages plotted in Extended Data Figure 8g were restricted to the former (square B4, B5, and C4), although this chart includes radiocarbon ages from both areas (C3, C4, C5, E3, and E4) except for the data from square D2 (Wk43606). In order to look at tendency in each area, in Table 1, we compiled two regional data sets from Supplementary Tables 2, 5, and 6 (ref. 1), focusing on chronologies between the –90 and –200 cm levels from where the radiocarbon data are available. In the area around E2, two radiocarbon ages from the hearths at about –100 cm level (Wk43610 and Wk43606) are younger than the OSL age 10 cm above (NE3) by ~3,000–7,000 years. Overall, the radiocarbon data suggests a younger chronology than the OSL results, or more appropriately, the number of reported radiocarbon ages is too small to settle this question."

So, let's have a look at the correlation between the C14 and OSL ages at the site.

# prepare data
excluding <- c("Wk43605" , # pit feature
               #"Wk43606", 
               #"Wk43607", 
               #"Wk43610", 
               #"Wk43604",
               #"Wk43611",
               #"OZT591",
               "OZT593"
  )

# excluding <- NA

squares <- "B3|B4|B5|C3|C4|C5|D3|D4|D5|E3|E4|E5"

c14_ages_excludes <-  c14_ages[!grepl(paste0(excluding, collapse = "|"), c14_ages$Lab.ID), ]
c14_ages_excludes <- c14_ages_excludes[grepl(paste0(squares, collapse = "|"), c14_ages_excludes$square), ]

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(c14_ages_excludes$Bchron_Median/1000, 
                                        osl_ages$osl_age),
                               depths = c(c14_ages_excludes$depth_below_surface,
                                          osl_ages$total_station_depth_below_surf),
                               group = c(rep("c14", nrow(c14_ages_excludes)),
                                         rep("osl", nrow(osl_ages))))

ages_with_groups <- 
  ages_with_groups %>% 
  filter(!is.na(depths)) %>% 
  mutate(ages = as.numeric(as.character(ages)),
         depths = abs(depths))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
all_osl_c14_BCDE_3_5 <- anova(fit0, fit1)
pval <- round(all_osl_c14_BCDE_3_5$`Pr(>F)`[2], 4)
fval <- round(all_osl_c14_BCDE_3_5$F[2], 3)
df <- all_osl_c14_BCDE_3_5$Df[2]

library(ggrepel)

formula <- y ~ poly(x, 2)

p1 <- 
ggplot() +
    geom_point(data = osl_ages,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              method = "lm", 
              formula = formula, 
              #fullrange = TRUE,
             colour = viridis(3)[1]) +
    geom_text_repel(data = osl_ages,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 1.5) +
    geom_point(data =c14_ages_excludes,
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[3]) +
    geom_smooth(data = c14_ages_excludes,
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
             # fullrange = TRUE,
              colour = viridis(3)[3]) +
      geom_text_repel(data = c14_ages_excludes,
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 1.5) +
   annotate("text", x = 2.7, y = 10, 
            label = paste0("Excluding:\n", paste0(excluding, collapse = "\n")), 
            size = 2)  +
   annotate("text",  x = 0.7, y = 75, 
            label = paste0("Age-depth model comparison ANOVA:\n", 
                           "F = ", fval, ", df = ", df, ", p = ", pval, " (quadratic fits)",
                           "\nModels are ", ifelse(pval <= 0.05, "", "not"), " significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("All OSL ages with all C14 ages from B-C-D-E 3-5")

ggsave("../../figures/All OSL ages with all C14 ages from B-C-D-E 3-5.png")

For the BCA: 

library(cowplot)
plot_grid(p4, p1, 
          labels = c("A", "B"), 
          align = "v", 
          ncol = 1)

Let's check for differences between the two age-depth models. As a general approach to detecting a significant difference between age-depth models, we first make linear models of the age-depth relationship for sub-sets of the samples (i.e. one model for all OSL samples, and one model for all C14 samples). We found that a quadratic model (i.e. second-order polynomial) is a simple model that fits all the subsets of our data very well. Then we compute a model that combines the data from two sub-sets, and adds a third variable, the grouping variable, that identifies the two models. We then compute an ANOVA to investigate the effect of the grouping variable on the fit of the models. If the ANOVA returns a significant difference, we conclude that the grouping variable has a significant effect, and thus the two sub-sets of the our data have significantly different age-depth models.

As we can see, the two models are not significantly different, so the effect of the grouping variable is not significant. Since the grouping variable introduced when combining th C14 and OSL ages into a single data set represents the respective datasets, the polynomial fits to the two datasets can be considered not significantly different. We only need to exclude one anomously deep C14 sample to find a non-significant difference between the age-depth curves of these two sets.

all_osl_c14_BCDE_3_5

OSL SWABC and all C14:

# prepare data
excluding <- c("Wk43605" , # pit feature
               #"Wk43606", 
               #"Wk43607", 
               #"Wk43610", 
               #"Wk43604",
               #"Wk43611",
               #"OZT591",
               "OZT593"
  )

# excluding <- NA

squares <- "B3|B4|B5|C3|C4|C5|D3|D4|D5|E3|E4|E5"

c14_ages_excludes <-  c14_ages[!grepl(paste0(excluding, collapse = "|"), c14_ages$Lab.ID), ]
c14_ages_excludes <- c14_ages_excludes[grepl(paste0(squares, collapse = "|"), c14_ages_excludes$square), ]

osl_ages_SWABC <-  osl_ages[grepl("SW.*A|SW.*B|SW.*C", osl_ages$Sample), ]

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(c14_ages_excludes$Bchron_Median/1000, 
                                        osl_ages_SWABC$osl_age),
                               depths = c(c14_ages_excludes$depth_below_surface,
                                          osl_ages_SWABC$total_station_depth_below_surf),
                               group = c(rep("c14", nrow(c14_ages_excludes)),
                                         rep("osl SWABC", nrow(osl_ages_SWABC))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
all_osl_SWABC_c14_BCDE_3_5 <- anova(fit0, fit1)
pval <- round(all_osl_SWABC_c14_BCDE_3_5$`Pr(>F)`[2], 4)
fval <- round(all_osl_SWABC_c14_BCDE_3_5$F[2], 3)
df <- all_osl_SWABC_c14_BCDE_3_5$Df[2]

library(ggrepel)

formula <- y ~ poly(x, 2)

p2 <- 
ggplot() +
    geom_point(data = osl_ages_SWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages_SWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              method = "lm", 
              formula = formula, 
              #fullrange = TRUE,
             colour = viridis(3)[1]) +
    geom_text_repel(data = osl_ages_SWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 1.5) +
    geom_point(data =c14_ages_excludes,
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[3]) +
    geom_smooth(data = c14_ages_excludes,
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
             # fullrange = TRUE,
              colour = viridis(3)[3]) +
      geom_text_repel(data = c14_ages_excludes,
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 1.5) +
   annotate("text", x = 2.7, y = 10, 
            label = paste0("Excluding:\n", paste0(excluding, collapse = "\n")), 
            size = 2)  +
   annotate("text",  x = 0.7, y = 75, 
            label = paste0("Age-depth model comparison ANOVA:\n", 
                            "F = ", fval, ", df = ", df, ", p = ", pval, " (quadratic fits)",
                           "\nModels are ", ifelse(pval <= 0.05, "", "not "), "significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages from SW A-B-C and C14 ages from B-C-D-E 3-5")

ggsave("../../figures/OSL ages from SW A-B-C and C14 ages from B-C-D-E 3-5.png")

OSL NWSWABC and all C14:

# prepare data
excluding <- c("Wk43605" , # pit feature
               #"Wk43606", 
               #"Wk43607", 
               #"Wk43610", 
               #"Wk43604",
               #"Wk43611",
               #"OZT591",
               "OZT593"
  )

# excluding <- NA

squares <- "B3|B4|B5|C3|C4|C5|D3|D4|D5|E3|E4|E5"

c14_ages_excludes <-  c14_ages[!grepl(paste0(excluding, collapse = "|"), c14_ages$Lab.ID), ]
c14_ages_excludes <- c14_ages_excludes[grepl(paste0(squares, collapse = "|"), c14_ages_excludes$square), ]

osl_ages_NWSWABC <-  osl_ages[grepl("NW.|SW.*A|SW.*B|SW.*C", osl_ages$Sample), ]
osl_ages_NE <- osl_ages[grepl("NE.", osl_ages$Sample), ]

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(c14_ages_excludes$Bchron_Median/1000, 
                                        osl_ages_NWSWABC$osl_age),
                               depths = c(c14_ages_excludes$depth_below_surface,
                                          osl_ages_NWSWABC$total_station_depth_below_surf),
                               group = c(rep("c14", nrow(c14_ages_excludes)),
                                         rep("osl NWSWABC", nrow(osl_ages_NWSWABC))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
all_osl_NWSWABC_c14_BCDE_3_5 <- anova(fit0, fit1)
pval <- round(all_osl_NWSWABC_c14_BCDE_3_5$`Pr(>F)`[2], 4)
fval <- round(all_osl_NWSWABC_c14_BCDE_3_5$F[2], 3)
df <- all_osl_NWSWABC_c14_BCDE_3_5$Df[2]

library(ggrepel)

formula <- y ~ poly(x, 2)

p3 <- 
ggplot() +
        geom_point(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[2]) +
      geom_text_repel(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 1.5) +
    geom_point(data = osl_ages_NWSWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages_NWSWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              method = "lm", 
              formula = formula, 
              #fullrange = TRUE,
             colour = viridis(3)[1]) +
    geom_text_repel(data = osl_ages_NWSWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 1.5) +
    geom_point(data =c14_ages_excludes,
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[3]) +
    geom_smooth(data = c14_ages_excludes,
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
             # fullrange = TRUE,
              colour = viridis(3)[3]) +
      geom_text_repel(data = c14_ages_excludes,
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 1.5) +
   annotate("text", x = 2.7, y = 10, 
            label = paste0("Excluding:\n", paste0(excluding, collapse = "\n")), 
            size = 2)  +
   annotate("text",  x = 0.5, y = 75, hjust = 0,
            label = paste0("Age-depth model comparison ANOVA\n",
                           "Comparing OSL ages from NW & SW A-B-C to\n",
                            "C14 ages from B-C-D-E 3-5:\n",
                            "F = ", fval, ", df = ", df, ", p = ", pval, " (quadratic fits)",
                           "\nModels are ", ifelse(pval <= 0.05, "", "not "), "significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("Age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages and\nC14 ages from B-C-D-E 3-5")

ggsave("../../figures/OSL ages from NW, SW A-B-C and C14 ages from B-C-D-E 3-5.png")

Comparing OSL for B4 (SW A) and B5 (SW B & C):

# prepare data
osl_ages_SWA <-  osl_ages[grepl("SW.*A", osl_ages$Sample), ]
osl_ages_SWBC <-  osl_ages[grepl("SW.*B|SW.*C", osl_ages$Sample), ]
# osl_ages_SWBC <- osl_ages_SWBC %>% filter(!Sample %in% c("SW4C",  X
#                                                          "SW7B", #
#                                                          "SW6B", X
#                                                          "SW7C", X
#                                                          "SW8C")) X

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(osl_ages_SWA$osl_age,
                                        osl_ages_SWBC$osl_age),
                               depths = c(osl_ages_SWA$total_station_depth_below_surf,
                                          osl_ages_SWBC$total_station_depth_below_surf),
                               group = c(rep("B4 (SW A)", nrow(osl_ages_SWA)),
                                         rep("B5 (SW B & C)", nrow(osl_ages_SWBC))))

# subset only phase 2, 2.6 and 2.15 m
ages_with_groups <- 
  ages_with_groups %>% 
  filter(depths <= -2.15) %>% 
  filter(depths >= -2.6)

#model without grouping variable
n <- 2
fit0 <- lm(ages ~ poly(depths, n), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, n) * group, data = ages_with_groups)

#compare models 
all_osl_B4_B5 <- anova(fit0, fit1)
pval <- round(all_osl_B4_B5$`Pr(>F)`[2], 3)
fval <- round(all_osl_B4_B5$F[2], 3)
df <- all_osl_B4_B5$Df[2]

library(ggrepel)
formula <- y ~ poly(x, n)

p4 <- 
ggplot() +
  geom_rect(data = data_frame(xmin = 2.15, 
                            xmax = 2.6, 
                            ymin = 0, 
                            ymax = 100),
          aes(xmin=xmin, 
              xmax=xmax, 
              ymin=ymin, 
              ymax=ymax), 
          fill = NA,
          colour = "black") +
  annotate("text", 
           x=2.375, 
           y=25, 
           hjust = 0.5,
           label = "Phase 2") +
    geom_point(data = osl_ages_SWA,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages_SWA,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_SWA,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
    geom_point(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[3]) +
    geom_smooth(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[3],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
   annotate("text",  x = 0.5, y = 75, hjust = 0,
            label = paste0("Age-depth model comparison ANOVA for phase 2 only\n",
                           "Comparing OSL ages from B4 (SW A) and B5 (SW B & C)\n",

                            "F = ", 
                           fval, 
                           ", df = ", 
                           df, 
                           ", p = ", 
                           pval, " (quadratic fits)",
                           "\nModels are ", 
                           ifelse(pval <= 0.05, "", "not "), 
                           "significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("Age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("Phase 2 OSL ages from B4 (SW A) and B5 (SW B & C)") 


ggsave("../../figures/OSL ages from B4 (SW A) and B5 (SW B & C).png")

Comparing OSL SWA B-C and NW:

# prepare data
osl_ages_SWBC <-  osl_ages[grepl("SW.*B|SW.*C", osl_ages$Sample), ]
osl_ages_NW <-  osl_ages[grepl("NW", osl_ages$Sample), ]

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(osl_ages_SWBC$osl_age,
                                        osl_ages_NW$osl_age),
                               depths = c(osl_ages_SWBC$total_station_depth_below_surf,
                                          osl_ages_NW$total_station_depth_below_surf),
                               group = c(rep("SW B & C", nrow(osl_ages_SWBC)),
                                         rep("NW", nrow(osl_ages_NW))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
all_osl_SWBC_NW <- anova(fit0, fit1)
pval <- round(all_osl_SWBC_NW$`Pr(>F)`[2], 3)

library(ggrepel)
formula <- y ~ poly(x, 2)

p5 <- 
ggplot() +
    geom_point(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
    geom_point(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[3]) +
    geom_smooth(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[3],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
     annotate("text",  x = 0.7, y = 75, 
            label = paste0("Age-depth model comparison ANOVA:\n", 
                           "p = ", pval, " (quadratic fits)",
                           "\nModels are ", ifelse(pval <= 0.05, "", "not"), " significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages from C4 (NW) and B5 (SW B & C)")

ggsave("../../figures/OSL ages from C4 (NW) and B5 (SW B & C).png")

OSL front and back, NE and SWBC:

# prepare data
osl_ages_SWBC <-  osl_ages[grepl("SW.*B|SW.*C", osl_ages$Sample), ]
osl_ages_NE <-  osl_ages[grepl("NE.", osl_ages$Sample), ]

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(osl_ages_SWBC$osl_age,
                                        osl_ages_NE$osl_age),
                               depths = c(osl_ages_SWBC$total_station_depth_below_surf,
                                          osl_ages_NE$total_station_depth_below_surf),
                               group = c(rep("SW B & C", nrow(osl_ages_SWBC)),
                                         rep("NE", nrow(osl_ages_NE))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
all_osl_SWBC_NE <- anova(fit0, fit1)
pval <- round(all_osl_SWBC_NE$`Pr(>F)`[2], 3)

library(ggrepel)
formula <- y ~ poly(x, 2)

p6 <- 
ggplot() +
    geom_point(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
    geom_point(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[3]) +
    geom_smooth(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[3],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_SWBC,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
     annotate("text",  x = 0.7, y = 75, 
            label = paste0("Age-depth model comparison ANOVA:\n", 
                           "p = ", pval, " (quadratic fits)",
                           "\nModels are ", ifelse(pval <= 0.05, "", "not"), " significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages from front (NE) and back (SW B & C)")

ggsave("../../figures/OSL ages from front (NE) and back (SW B & C).png")

OSL SWABC and NW:

# prepare data
osl_ages_SWABC <-  osl_ages[grepl("SW.*A|SW.*B|SW.*C", osl_ages$Sample), ]
osl_ages_NW <-  osl_ages[grepl("NW.", osl_ages$Sample), ]

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(osl_ages_SWABC$osl_age,
                                        osl_ages_NW$osl_age),
                               depths = c(osl_ages_SWABC$total_station_depth_below_surf,
                                          osl_ages_NW$total_station_depth_below_surf),
                               group = c(rep("SW A, B & C", nrow(osl_ages_SWABC)),
                                         rep("NW", nrow(osl_ages_NW))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
all_osl_SWABC_NW <- anova(fit0, fit1)
pval <- round(all_osl_SWABC_NW$`Pr(>F)`[2], 3)

library(ggrepel)
formula <- y ~ poly(x, 2)

p7 <- 
ggplot() +
    geom_point(data = osl_ages_SWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages_SWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_SWABC,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
    geom_point(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[3]) +
    geom_smooth(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[3],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
     annotate("text",  x = 0.7, y = 75, 
            label = paste0("Age-depth model comparison ANOVA:\n", 
                           "p = ", pval, " (quadratic fits)",
                           "\nModels are ", ifelse(pval <= 0.05, "", "not"), " significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages from NW and SW A, B & C")

ggsave("../../figures/OSL ages from NW and SW A, B & C.png")

OSL NE and NW:

# prepare data
osl_ages_NE <-  osl_ages[grepl("NE.", osl_ages$Sample), ]
osl_ages_NW <-  osl_ages[grepl("NW.", osl_ages$Sample), ]

# test for difference in age-depth models

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(osl_ages_NE$osl_age,
                                        osl_ages_NW$osl_age),
                               depths = c(osl_ages_NE$total_station_depth_below_surf,
                                          osl_ages_NW$total_station_depth_below_surf),
                               group = c(rep("NE", nrow(osl_ages_NE)),
                                         rep("NW", nrow(osl_ages_NW))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
all_osl_NE_NW <- anova(fit0, fit1)
pval <- round(all_osl_NE_NW$`Pr(>F)`[2], 3)

library(ggrepel)
formula <- y ~ poly(x, 2)

p8 <- 
ggplot() +
    geom_point(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_NE,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
    geom_point(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[3]) +
    geom_smooth(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[3],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages_NW,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
     annotate("text",  x = 0.7, y = 75, 
            label = paste0("Age-depth model comparison ANOVA:\n", 
                           "p = ", pval, " (quadratic fits)",
                           "\nModels are ", ifelse(pval <= 0.05, "", "not "), "significantly different" ), 
            size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages from NW and NE")

ggsave("../../figures/OSL ages from NW and NE.png")

library(ggrepel)
formula <- y ~ poly(x, 2)

p9 <- 
ggplot() +
    geom_point(data = osl_ages,
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages,
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages,
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
    geom_point(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[3]) +
    geom_smooth(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
              fullrange = TRUE,
              colour = viridis(3)[3]) +
      geom_text_repel(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("All OSL with C14 ages only from squares D and E")

Unlike the above models, these two age-depth models of all OSL ages and C14 ages only from squares D and E are significantly different. A significant interaction is reported by the ANOVA comparing the null model with the model that includes the grouping variable (i.e. the dating method). This indicates that the dating method has a significant effect on the age-depth relationship in this case. The radiocarbon ages are signficantly younger than the OSL ages.

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(c14_ages[grepl("D|E", c14_ages$square_spit), ]$Bchron_Median/1000, 
                                        osl_ages$osl_age),
                               depths = c(c14_ages[grepl("D|E", c14_ages$square_spit), ]$depth_below_surface,
                                          osl_ages$total_station_depth_below_surf),
                               group = c(rep("c14", nrow(c14_ages[grepl("D|E", c14_ages$square_spit), ])),
                                         rep("osl", nrow(osl_ages))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
anova(fit0, fit1)

ggplot() +
  geom_point(data = c14_ages[grepl("B|C", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[1]) +
  geom_smooth(data = c14_ages[grepl("B|C", c14_ages$square_spit), ],
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
              fullrange = TRUE,
              colour = viridis(3)[1]) +
  geom_text_repel(data = c14_ages[grepl("B|C", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 2) +
  geom_point(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[3]) +
  geom_smooth(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
              fullrange = TRUE,
              colour = viridis(3)[3]) +
  geom_text_repel(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("C14 ages from squares D & E vs B & C")

In comparing C14 ages from the two sides of the excavation, we see no significant difference in the age-depth models for samples taken from squares D & E compared to squares B & C: 

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(c14_ages[grepl("B|C", c14_ages$square_spit), ]$Bchron_Median/1000, 
                                        c14_ages[grepl("D|E", c14_ages$square_spit), ]$Bchron_Median/1000),
                               depths = c(c14_ages[grepl("B|C", c14_ages$square_spit), ]$depth_below_surface,
                                          c14_ages[grepl("D|E", c14_ages$square_spit), ]$depth_below_surface),
                               group = c(rep("c14 BC", nrow(c14_ages[grepl("B|C", c14_ages$square_spit), ])),
                                         rep("osl DE", nrow(c14_ages[grepl("D|E", c14_ages$square_spit), ]))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
anova(fit0, fit1)

ggplot() +
      geom_point(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
  geom_point(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[3]) +
  geom_smooth(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
              fullrange = TRUE,
              colour = viridis(3)[3]) +
  geom_text_repel(data = c14_ages[grepl("D|E", c14_ages$square_spit), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages from NE column & C14 ages from squares D & E")

OSL ages from the NE column present a significantly different age-depth curve to the C14 ages from squares D and E. They do seem to have a point about the difficulty of correlating ages from the two different methods across the site. 

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(c14_ages[grepl("D|E", c14_ages$square_spit), ]$Bchron_Median/1000, 
                                        osl_ages[grepl("NE", osl_ages$Sample), ]$osl_age),
                               depths = c(c14_ages[grepl("D|E", c14_ages$square_spit), ]$depth_below_surface,
                                          osl_ages[grepl("NE", osl_ages$Sample), ]$total_station_depth_below_surf),
                               group = c(rep("c14", nrow(c14_ages[grepl("D|E", c14_ages$square_spit), ])),
                                         rep("osl", nrow(osl_ages[grepl("NE", osl_ages$Sample), ]))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
anova(fit0, fit1)

back <- "1|2"
front <- "3|4|5"

ggplot() +
  # back
    geom_point(data = c14_ages[grepl(back, c14_ages$square), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[1]) +
  geom_smooth(data = c14_ages[grepl(back, c14_ages$square), ],
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
              fullrange = TRUE,
              colour = viridis(3)[1]) +
  geom_text_repel(data = c14_ages[grepl(back, c14_ages$square), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 2) +
  # front
      geom_point(data = c14_ages[grepl(front, c14_ages$square), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000),
             colour = viridis(3)[3]) +
  geom_smooth(data = c14_ages[grepl(front, c14_ages$square), ],
              aes(-depth_below_surface, 
                 Bchron_Median/1000),
              method = "lm", 
              formula = formula, 
              fullrange = TRUE,
              colour = viridis(3)[3]) +
  geom_text_repel(data = c14_ages[grepl(front, c14_ages$square), ],
             aes(-depth_below_surface, 
                 Bchron_Median/1000,
                 label = Lab.ID),
             size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("C14 ages from back (rows 1, 2) and front (rows 3, 4, 5)")

No significant different in age-depth curves of C14 ages from the front of the site compared to C14 ages from the back of the site. 

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(c14_ages[grepl(front, c14_ages$square), ]$Bchron_Median/1000, 
                                        c14_ages[grepl(back, c14_ages$square), ]$Bchron_Median/1000),
                               depths = c(c14_ages[grepl(front, c14_ages$square), ]$depth_below_surface,
                                          c14_ages[grepl(back, c14_ages$square), ]$depth_below_surface),
                               group = c(rep("c14 front", nrow(c14_ages[grepl(front, c14_ages$square), ])),
                                         rep("c14 back", nrow(c14_ages[grepl(back, c14_ages$square), ]))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
anova(fit0, fit1)

nudge <- 0
ggplot() +
      geom_point(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf + nudge, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf + nudge, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf + nudge, 
                 osl_age,
                 label = Sample),
             size = 2) +
  geom_point(data = osl_ages[grepl("SW", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[3]) +
    geom_smooth(data = osl_ages[grepl("SW", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[3],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages[grepl("SW", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle("OSL ages from NE column & OSL ages from SW columns")

There is a very significant difference in the age-depth curves of the NE OSL ages compared to the SW OSL ages. We do not seem to have a good understanding of how the NE side of the deposit relates to the SW side.

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(osl_ages[grepl("NE", osl_ages$Sample), ]$osl_age, 
                                         osl_ages[grepl("SW", osl_ages$Sample), ]$osl_age),
                               depths = c(osl_ages[grepl("NE", osl_ages$Sample), ]$total_station_depth_below_surf + nudge,
                                           osl_ages[grepl("SW", osl_ages$Sample), ]$total_station_depth_below_surf),
                               group = c(rep("osl NE", nrow(osl_ages[grepl("NE", osl_ages$Sample), ])),
                                         rep("osl SW", nrow( osl_ages[grepl("SW", osl_ages$Sample), ]))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
anova(fit0, fit1)

If we nudge the depths of the samples from the NE side up about 0.45 m, then they line up very well. However the age-depth models are still not equivalent, according to the ANOVA. 

nudge <- 0.45
ggplot() +
      geom_point(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf + nudge, 
                 osl_age),
             colour = viridis(3)[1]) +
    geom_smooth(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf + nudge, 
                 osl_age),
              colour = viridis(3)[1],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages[grepl("NE", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf + nudge, 
                 osl_age,
                 label = Sample),
             size = 2) +
  geom_point(data = osl_ages[grepl("SW", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age),
             colour = viridis(3)[3]) +
    geom_smooth(data = osl_ages[grepl("SW", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age),
              colour = viridis(3)[3],
              method = "lm", 
              formula = formula, 
              fullrange = TRUE) +
    geom_text_repel(data = osl_ages[grepl("SW", osl_ages$Sample), ],
             aes(-total_station_depth_below_surf, 
                 osl_age,
                 label = Sample),
             size = 2) +
  theme_minimal() +
  ylab("age ka") +
  xlab("Depth below surface (m)") +
  ggtitle(paste0("OSL ages from NE column (nudged up ", nudge, "m \nfor demonstration purposes) & OSL ages from SW columns"))

# combine data with grouping variable
ages_with_groups <- data.frame(ages = c(osl_ages[grepl("NE", osl_ages$Sample), ]$osl_age, 
                                         osl_ages[grepl("SW", osl_ages$Sample), ]$osl_age),
                               depths = c(osl_ages[grepl("NE", osl_ages$Sample), ]$total_station_depth_below_surf + nudge,
                                           osl_ages[grepl("SW", osl_ages$Sample), ]$total_station_depth_below_surf),
                               group = c(rep("osl NE", nrow(osl_ages[grepl("NE", osl_ages$Sample), ])),
                                         rep("osl SW", nrow( osl_ages[grepl("SW", osl_ages$Sample), ]))))

#model without grouping variable
fit0 <- lm(ages ~ poly(depths, 2), data = ages_with_groups)

#model with grouping variable
fit1 <- lm(ages ~ poly(depths, 2) * group, data = ages_with_groups)

#compare models 
anova(fit0, fit1)


benmarwick/mjbnaturepaper documentation built on May 9, 2022, 6:25 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close