library(knitr) opts_chunk$set(fig.path='figure/')
# make sure dplyr comes first select <- dplyr::select library(PiperetalGMHBone) # First load the data and code libraries. spit_depths <- read.csv("GMH_spit_depths.csv", stringsAsFactors = FALSE) dates_input <- read.csv("GMH_C14_dates.csv", stringsAsFactors = FALSE) my_libraries()
In this section we calibrate radiocarbon dates with OxCal 4.2 & IntCal13. We manaully send the raw radiocarbon ages to OxCal and merge the Oxcal output back in so we can work with the depths and ages.
# need to have: # Name = lab code, Date = radiocarbon age, Uncertainty = error dates_for_oxcal <- with(dates_input, data.frame(Name = DirectAMS.code, Date = BP, Uncertainty = X1s.error.1)) dates <- dates_for_oxcal oxcal_formatter(dates) # now follow instructions in oxcal_formatter.R to use OxCal website # now load output from OxCal website and put output in data folder dates <- read.csv("GMH_C14_dates_OxCal_output.csv", stringsAsFactors = FALSE) # combine with raw input table, prepare by making colnames useful colnames(dates) <- c("Name", as.character(dates[1,])[2:7]) dates$DirectAMS.code <- gsub("R_Date ", "", dates$Name) dates[,2:7] <- sapply(dates[,2:7], as.numeric) dates_calib <- merge(dates_input, dates, by = 'DirectAMS.code') # don't want to see the spit column as those are an early rough guess # not a computed result dates_calib <- select(dates_calib, -Excavation.unit..spit.) # tidy up a little, remove spaces dates_calib$material <- gsub("\\s", "", dates_calib$X)
Here's the OxCal output:
dates_calib
Compute separate linear models for the ages obtained from the shell and the organics, verify that the models show a non-random association of age and depth, then compute the difference in intercept.
# by material type # function to extract a, b, r2, and p values lm_eqn = function(df){ m = lm(median ~ Depth.below.surface..m., df); out <-data.frame(a = coef(m)[1], b = coef(m)[2], r2 = summary(m)$r.squared, p = anova(m)$'Pr(>F)'[1], n = nrow(df) ) out } # compute linear models models_x <- plyr::dlply(dates_calib, "material", function(i) lm_eqn(i)) # present in tidy data frame lm_coeffs <- do.call("rbind", models_x) # for organics only: age = 3030.2 * depth -342.5 # check for differences in intercept to get offset intercepts <- round(lm_coeffs$a,2) difference <- round(sum(abs(intercepts)), 1) # plot limits <- with(dates_calib, ggplot2::aes(ymax = from, ymin = to)) ggplot(dates_calib, aes(Depth.below.surface..m., median, colour = material)) + # all, by colour geom_point(size = 4, aes(shape = material)) + geom_pointrange(limits, width = 2) + geom_smooth(se = FALSE, method = "lm") + coord_flip() + theme_minimal() + scale_x_reverse() + xlab("meters below surface") + ylab("cal yrs BP") + theme(axis.title.x=element_text(colour="black", size = 15), # increase axis title size slightly axis.title.y=element_text(colour="black", angle=90, size = 15), # increase axis title size slightly and rotate axis.text.x=element_text(colour="black", size = 15), # increase size of numbers on x-axis axis.text.y=element_text(colour="black", size = 15)) # increase size of numbers on y-axis
Here's a table of the linear model coefficients:
lm_coeffs
The difference in intercepts in the linear models of the shell dates and organic dates is r difference
. This indicates an approximate difference of r difference
years between paired shell and organic dates.
If we want to use a linear model to predict age from depth, using only the ages from organic samples, the equation would be:
age = r round(lm_coeffs[1,1],2)
+ r round(lm_coeffs[1,2],2)
* depth
Convert the field measurements to be relative to the datum, then convert to to be relative to the ground surface (ie. the start depths of spit one).
# relate every spit depth back to the datum spit_depths[] <- sapply(spit_depths, as.numeric) relative_to_datum <- spit_depths$datum - spit_depths[,4:13] # take average of start levels of spit one for # ground surface relative to datum ground_surface <- mean(as.numeric(relative_to_datum[1,1:4])) # now compute relative to ground surface, this should give # values that we can read as actual depths below the surface relative_to_ground_surface <- relative_to_datum + abs(ground_surface) relative_to_ground_surface$spit <- spit_depths$Spit relative_to_ground_surface # plot spits_m <- melt(relative_to_ground_surface, id.var = 'spit') spits_m$start_or_end <- ifelse(grepl("end", spits_m$variable), "end", "start") remove <- c(".end", ".start") spits_m$spit_start_or_end <- with(spits_m, paste0(spit, "_", start_or_end)) spits_m$variable <- gsub(paste(remove, collapse = "|"), "", spits_m$variable) ggplot(spits_m, aes(variable, value, group = spit_start_or_end, linetype = start_or_end, colour = start_or_end)) + geom_line() + theme_minimal() + scale_y_continuous(breaks = seq(-3,0,0.25)) + xlab("Measurement location") + ylab("Depth below surface (m)") # how deep was spit nine? nine_start <- round(mean(filter(spits_m, spit == 9, start_or_end == "start")$value),2) nine_end <- round(mean(filter(spits_m, spit == 9, start_or_end == "end")$value),2)
The measured spit depths indicate a maximum depth below the surface of about 2.6 m. The rocky nature of the base of the excavation made it difficult to measure depth accurately and it's likely that the true depth is slightly greater. The section drawing indicates a depth of 2.7 m below the surface for the lowest points around the rocks, which would have been measured by tape measure, rather than the staff used for the spit depths. So 2.6-2.7m is a reliable maximum depth value.
The lowest pottery is in spit nine, which is between r -nine_start
and r -nine_end
meters below the surface.
The dated samples were taken from the south wall and their depths can be measured from the section drawing. With these depths we can find the spit depths from the SW and SE corners that will identify the spit that correlates with the dated sample.
dates_calib$section <- "SW." # plot spits with dated samples ggplot(spits_m, aes(variable, value)) + geom_line(aes(group = spit_start_or_end, linetype = start_or_end, colour = start_or_end)) + geom_text(data = dates_calib, aes(section, -Depth.below.surface..m., label = DirectAMS.code), size = 1, hjust = 1.2) + theme_minimal() + scale_y_continuous(breaks = seq(-3,0,0.25)) + xlab("Measurement location") + ylab("Depth below surface (m)")
Here we can see which spits the dated sample fall into.
# match sample depths up with a specific spit from SW-SE data relative_to_ground_surface$South_means <- with(relative_to_ground_surface, apply(cbind(SE..end, SW..end), 1, mean)) South_depth_means <- sort(as.vector(na.omit(apply(cbind(relative_to_ground_surface$SE..end, relative_to_ground_surface$SW..end), 1, mean))), decreasing = FALSE) intervals <- cbind(-dates_calib$Depth.below.surface..m., findInterval(-dates_calib$Depth.below.surface..m., South_depth_means)) means <- merge(data.frame(South_depth_means), relative_to_ground_surface, by.x = 'South_depth_means', by.y = 'South_means') means <- arrange(means, South_depth_means) # subset the table of spit depths with the intervals that the samples come from, # and these are the spits that correlate with the depths that the dated samples # were taken from in the South wall df <- (list(means[intervals[,2], ]$spit, dates_calib$Depth.below.surface..m., dates_calib$Name, dates_calib$material, dates_calib$BP, dates_calib$to, dates_calib$from)) df <- data.frame(sapply(df,'[',seq(max(sapply(df,length))))) names(df) <- c("spit", "depth below surface", "code", "material", "uncal age","to cal BP", "from cal BP") df
The implication for bone point phases is that it puts the start of phase C squarely in spit 17. Otherwise it's consistent with David's calculations.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.