Supplementary details supporting the analysis of OSL and C14 dates, lithics and shellfish from the 1989 excavations at Madjebebe (Malakanunja II)

Ben Marwick, benmarwick@gmail.com

This document contains the source code used to generate the plots and statistics presented in the manuscript. It is provided to give the reader additional insight into the data analysis workflow and allow anyone to reproduce our analysis and further explore the analyses and data for themselves.

OSL and C14 dates

chrono_data_tidied

library(mjb1989excavationpaper)
chrono_data <- get_chrono_data()

OK, let’s have a look at what get_chrono_data actually does. You can print the body of any function by calling it without a parenthesis.

get_chrono_data

function () 
{
    dates_1989 <- read.csv("data/Date_table_from_paper_on_1989_dig.csv")
    return(dates_1989)
}
<environment: namespace:mjb1989excavationpaper>

And let’s have a quick at the raw data too:

chrono_data_tidied

Let’s see the chrono_data_tidied data

interpolate_chrono_data

interpolate_chrono_data

function (dates_1989) 
{
    cal.date.lo <- loess(age2plot ~ Depth.m, dates_1989[dates_1989$Method %in% 
        c("C14", "TL", "OSL"), ], span = 0.4)
    cal.date.pr <- predict(cal.date.lo, data.frame(Depth.m = seq(0, 
        5, 0.01)))
    cal.date.pr <- data.frame(age = cal.date.pr, depth = names(cal.date.pr))
    oldest_depth = 287
    oldest <- cal.date.pr[cal.date.pr$depth == oldest_depth, 
        ]$age
    base_of_dense_depth = 250
    dense <- cal.date.pr[cal.date.pr$depth == base_of_dense_depth, 
        ]$age
    return(list(oldest_depth = oldest_depth, oldest = oldest, 
        base_of_dense_depth = base_of_dense_depth, dense = dense, 
        cal.date.lo = cal.date.lo))
}
<bytecode: 0x10d04fd10>
<environment: namespace:mjb1989excavationpaper>

plot_chrono_data

function (dates_1989, oldest_depth, oldest, base_of_dense_depth, 
    dense) 
{
    library(ggplot2)
    library(grid)
    limits <- aes(ymax = age2plot + Error/1000, ymin = age2plot - 
        Error/1000, colour = factor(Method))
    library(scales)
    p1 <- ggplot(dates_1989, aes(Depth.m, age2plot, label = Lab.Code)) + 
        geom_smooth(data = dates_1989[!(dates_1989$Lab.Code %in% 
            c("ANUA-9915")), ], method = "loess", se = FALSE, 
            span = 0.5, colour = "grey", alpha = 0.2, size = 1) + 
        geom_point(size = 2, aes(colour = factor(Method), shape = factor(Method))) + 
        geom_linerange(limits) + geom_text(angle = 0, hjust = -0.2, 
        vjust = 0.2, size = 2) + xlab("meters below surface") + 
        ylab("x 1000 years cal. BP") + scale_colour_discrete(name = "Dating\nmethod") + 
        scale_shape_discrete(name = "Dating\nmethod") + scale_y_continuous(labels = comma, 
        breaks = seq(0, 100, 20)) + scale_x_reverse() + annotate("segment", 
        x = oldest_depth/100, y = oldest, xend = 4.75, yend = oldest, 
        colour = "grey30") + annotate("segment", x = oldest_depth/100, 
        y = 0, xend = oldest_depth/100, yend = oldest, colour = "grey30") + 
        annotate("text", x = 3, y = 15, label = paste0("lowest artefacts (", 
            round(oldest, 0), " ka BP)"), size = 3) + annotate("segment", 
        x = base_of_dense_depth/100, y = dense, xend = 4.75, 
        yend = dense, colour = "grey30") + annotate("segment", 
        x = base_of_dense_depth/100, y = 0, xend = base_of_dense_depth/100, 
        yend = dense, , colour = "grey30") + annotate("text", 
        x = 2.5, y = 15, label = paste0("base of dense \noccupation layer (", 
            round(dense, 0), " ka BP)"), size = 3) + theme_bw() + 
        coord_flip() + theme(axis.title.x = element_text(size = 17), 
        axis.title.y = element_text(size = 17, angle = 90), axis.text.x = element_text(size = 15, 
            colour = "black"), axis.text.y = element_text(size = 15, 
            colour = "black"))
    ggsave(file = "figures/Fig_5_MKII_dates_from_1989_for_paper.svg", 
        width = 200, height = 200, units = "mm")
    return(p1)
}
<bytecode: 0x111866cd8>
<environment: namespace:mjb1989excavationpaper>

chrono_data_plotted <- plot_chrono_data(chrono_data_tidied, 
                        oldest_depth = chrono_data_interpolated$oldest_depth,
                        oldest = chrono_data_interpolated$oldest,
                        base_of_dense_depth = chrono_data_interpolated$base_of_dense_depth,
                        dense = chrono_data_interpolated$dense)
chrono_data_plotted

Based on the TL age estimates and the artefact distribution, Roberts et al. (1990a) suggested first occupation at the site began 55 ± 5 ka. Although the excavators were conservative in their interpretation s€”stressing the upper (i.e. 50 ka) limit of this age range and taking the lower limit of the high density band of artefacts (2.4 m bs) as the actual level of initial occupation €” we confirm here that the lowest artefacts occur in Spit 49, 2.76 €“ 2.8 m deep below surface and are bracketed by the original OSL age estimate of 55.5 ± 8.2 ka (KTL-162) and the TL estimate of 65 ± 14 ka (KTL-141). Subsequent redating of several of the lower samples at MJB using single grain and single aliquot OSL methods reduced the error ranges for the lower dates substantially, but did not alter the original results (Roberts et al. 1998).

At the time of Roberts et al. (1998) publication and commentary no calibration curve was available for radiocarbon dates greater than 11 ka, and hence the issue of underestimation of calendar years could not be resolved. A calibration curve is now available back to 50 ka (Reimer et al. 2013). The IntCal13 calibration curve results in a calibrated age for the 13.39 ka age (ANU-7006) as 15,001 €“ 17,429 cal BP, overlapping at 1sdƒ with the TL age of 15 ± 3 ka (KTL-165) at equivalent depth. Kamminga and Allen€™s (1973) 18.04 ka age (SUA-265) calibrates to c. 22.4 - 21.0 ka, overlapping at 1sdƒ with the TL age of 24 ± 5 ka (KTL-97) at equivalent depth. The 14.9 ka age (ANU-7007) provided by Roberts et al. (1998) calibrates to 18.5 - 17.8 ka, and remains within 1σ of the TL age of KTL-97. These new calibrations show that conventional radiocarbon years do substantially under-estimate the ages of the sediments and that calendar ages show a strong correlation with luminescence ages.

plotting_slopes_chrono_data

chrono_data_slopes_plot <- plotting_slopes_chrono_data(chrono_data_tidied)
# show plot
chrono_data_slopes_plot$the_plot

Chech the structure of the chrono_data_slopes_plot output.

str(chrono_data_slopes_plot, 1)

List of 2
 $ dates_1989:'data.frame': 24 obs. of  19 variables:
 $ the_plot  :List of 9
  ..- attr(*, "class")= chr [1:2] "gg" "ggplot"

Plotting the linear regression lines for OSL/TL ages against radiocarbon ages suggests that the slopes of the lines are not the same. We can explore this further by comparing the models directly.

anova_p_value <- testing_slopes_chrono_data(chrono_data_slopes_plot$dates_1989)

Computing a linear model with depth and age using the two dating methods as a cross to test for interaction between them gives a significant difference, with a p-value of 0.011.

ancova_p_values <- ancova_slopes_chrono_data(chrono_data_slopes_plot$dates_1989)
# show the plot
ancova_p_values$the_plot

If we restrict the sample to just the dates in the upper 2 m, which is the depth of the lowest in-sequence radiocarbon date (SUA-265), we can more robustly compare the ages from radiocabon and luminesensce methods.

In the ANCOVA model, age is modeled as the dependent variable with depth as the factor and dating method as the covariate. The summary of the results show a significant effect of depth (p = 0.00000), but no significant interaction between dating methods (p = 0.791). These results suggest that the slope of the regression between age and depth is similar for both dating methods.

anova_p_values <- anova_slopes_chrono_data(ancova_p_values$upper, ancova_p_values$m2)

We can fit a more parsimonious model without the interaction to test for significant differences in slope. We can compare the two models for each dating method with ANOVA to assess if removing the interaction significantly affects the fit of the model. The result is that removing the interaction does not significantly effect the fit of the model. There is still no significant effect of dating method on the slope of the age-depth lines (p = 0.791)

bayes_values <- bayes_slopes_chrono_data(ancova_p_values$upper)

Loading required package: coda
Loading required package: MASS
##
## Markov Chain Monte Carlo Package (MCMCpack)
## Copyright (C) 2003-2017 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
##
## Support provided by the U.S. National Science Foundation
## (Grants SES-0350646 and SES-0350613)
##

# Warning - very time consuming, >30 min
# bayes_test_values <- bayes_test_slopes_chrono_data(bayes_values$x1_post, bayes_values$x2_post)

# since the MCMC is quite time consuming,I've included the output
# from a run so we can use that rather than redo the analysis. To
# redo the analysis, simply run the code in the previous chunk.
load("data/slopes_test.rda")
# This requires the summary.BEST from BEST
library("BEST")

Loading required package: HDInterval

Attaching package: ‘BEST’

The following object is masked from ‘package:mjb1989excavationpaper’:

    plotPost

summary(slopes_test)

             mean  median    mode HDI%   HDIlo   HDIup compVal %>compVal
mu1       10.8219 10.8220 10.8231   95 10.7351 10.9060                  
mu2       10.7623 10.7623 10.7618   95 10.7358 10.7880                  
muDiff     0.0596  0.0598  0.0620   95 -0.0309  0.1482       0      90.4
sigma1     3.6205  3.6204  3.6230   95  3.5240  3.7115                  
sigma2     1.1386  1.1385  1.1386   95  1.1148  1.1625                  
sigmaDiff  2.4819  2.4818  2.4817   95  2.3938  2.5696       0     100.0
nu         4.6537  4.6485  4.6495   95  4.3116  5.0095                  
log10nu    0.6675  0.6673  0.6676   95  0.6352  0.7003                  
effSz      0.0222  0.0223  0.0233   95 -0.0108  0.0559       0      90.4

best1 <- round(BEST:::summary.BEST(slopes_test)[3], 2)
best2 <- round(unname(BEST:::summary.BEST(slopes_test)[3, 5:6]),2)

Smith and Jones obtained several more radiocarbon dates that have not been previously published (Table 1), increasing the number of available radiocarbon ages for the site to 13. With the exception of one anomalous age from the base of the sequence already discussed above (ANUA-9915), the additional 14C ages provide a picture of consistent depth-age relationships between radiocarbon and luminescence techniques down to two meters below the surface (the depth of the lowest in-sequence radiocarbon date). A test of the difference between the correlation coefficients for the linear regression slopes for 14C and luminescence ages indicates no significant difference between slopes, indicating that both provide effectively identical age-depth relationships for the uppermost two meters of deposit (Bayesian estimation of difference in means = 0.06, 95% HDI = -0.03, 0.15, the interval includes zero, indicating no credible difference). Concerns over the degree of fit between 14C and luminescence chronologies for the upper half of the deposit can no longer be sustained (cf. Bowdler 1990; Hiscock 1990).

#### bayesian change point analysis on sedimentation rates
bayes_cp_result <- bayes_cp_test(chrono_data_slopes_plot$dates_1989)

Bayesian change point analysis indicates that sedimentation rates slow substantially below 2.0 m to the base of occupation (posterior probability of change is 0.948 at 22 ka cal BP) and then accelerate again below the lowest occupation (posterior probability of change is 0.91 at 65 ka).

# check 15 €“ 20 ka sedimentation rates
# examine slope from KTL97 through KTL165 compare C14 to OSL
# now limit for depth range KTL97 through KTL165
sed_rates_15_to_20_result <- sed_rates_15_to_20(chrono_data_slopes_plot$dates_1989)
# plot
sed_rates_15_to_20_result$the_plot

bayes_slope_posteriors <- bayes_slope_difference(sed_rates_15_to_20_result$sub_OSL,
                                                 sed_rates_15_to_20_result$sub_C14)

# rather time consuming, >30 min
# bayes_slope_difference_result <- bayes_slope_test(bayes_slope_difference_result$x1_post,
#                                                bayes_slope_difference_result$x2_post)

# since the MCMC is quite time consuming,I've included the output
# from a run so we can use that rather than redo the analysis. To
# redo the analysis, simply run the code in the previous chunk.
load("data/slopes_test_sub.rda")
BEST:::summary.BEST(slopes_test_sub)

              mean   median      mode HDI%     HDIlo    HDIup compVal %>compVal
mu1       25.94944 25.94812 25.942351   95  2.55e+01 26.42176                  
mu2       25.73194 25.72959 25.731983   95  2.49e+01 26.54876                  
muDiff     0.21750  0.22006  0.252338   95 -7.46e-01  1.14991       0      67.6
sigma1     4.90146  4.88684  4.839097   95  4.20e+00  5.62274                  
sigma2     8.56242  8.53835  8.504374   95  7.36e+00  9.84656                  
sigmaDiff -3.66096 -3.64934 -3.705678   95 -5.14e+00 -2.26751       0       0.0
nu         1.00374  1.00257  1.000514   95  1.00e+00  1.01130                  
log10nu    0.00162  0.00112  0.000247   95  9.72e-09  0.00488                  
effSz      0.03137  0.03164  0.035904   95 -1.05e-01  0.16664       0      67.6

# print(slopes_test_sub)
# plot(slopes_test_sub)
# windows()
# plotAll(slopes_test_sub, credMass=0.8, ROPEm=c(-0.1,0.1), 
#         ROPEeff=c(-0.2,0.2), compValm=0.5) 
# pairs(slopes_test_sub)

Sedimentation rates indicated by 14C and OSL ages during the period 15 €“ 20 ka, just before the first major change in sedimentation rates, are not credibly different (Bayesian estimation of difference in means = 0.22, 95% HDI = -0.75, 1.15, the interval includes zero, indicating no credible difference). Although there are no reliable 14C ages older than 20 ka, the lack of difference in rates between the two methods before this time suggests that the change in sedimentation rates was a real event, rather than an effect of changes in the age-depth relationship between 14C and luminescence methods. The depth-age curve also allows us to estimate the likely age of the lowest artefact as approximately 64 ka, and the base of the dense artefact layer as approximately 55 ka (though it must be acknowledged, that the error ranges on these lowest luminescence ages are large). Any possible refinement of temporal resolution cannot be addressed prior to analysis of new samples collected in 2012. Non-local stone in the form of silcrete is present from Spit 47 (2.9m bs), suggesting human activity at this level, but artefact numbers do not increase markedly until Spit 45.

Lithics

# Now onto the lithic data...
lithics_data_cleaned <- clean_lithics_data(lithics_data)

lithics_data_plotted <- plots_lithics_data(lithics_data_cleaned,  bayes_cp_result$cal.date.lo)

Saving 7 x 7 in image

lithics_data_plotted$p1

# figures/Fig_6_lithics_over_time_from_1989_for_paper.svg

Stone artefacts are present in every spit of the 1989 excavation at MJB (Figure 6, Table 2). These show distinct pulses of accumulation, centred around 5, 7, 12.5, 18.4, 36, and 45 €“ 52 ka (Figure 6). Between Spits 37 €“ 39 (2.2 €“ 2.4 m bs) there are 1900 flaked stone artefacts from the 7 mm sieve (including 41 of exotic chert), numerous pieces of high-grade haematite (totalling 4.92 kg, including seven ground pieces), 143 g of red or yellow ochre, pieces of dolerite (presumed to be fragments of edge-ground axes) and fragments of grindstones; a further 568 flaked stone artefacts and 344 pieces of haematite were recovered from the pit fill at this level. All of these are associated with, or are slightly beneath, a TL age of 45 ± 9 ka. Beneath this, in levels dating 52 ± 11 to 61 ± 13 ka there were an additional 82 flaked stone artefacts and 52 pieces of haematite recovered from the 7 mm sieve. This is a comparatively large assemblage from a trench measuring only 1 x 1.5 m. It is unlikely the quartzite stone artefacts derive from downward displacement of artefacts from the predominately quartz assemblage that overlies it.

lithics_rawmaterials_plotted <- plots_lithics_rawmaterials(lithics_data, 
                                                           lithics_data_plotted$spit, 
                                                           bayes_cp_result$cal.date.pr)                                                     

Saving 7 x 7 in image

lithics_rawmaterials_stratplotted <- stratiplot_lithics(lithics_rawmaterials_plotted$plot2)

Loading required package: analogue
Loading required package: vegan
Loading required package: permute
Loading required package: lattice
This is vegan 2.4-4

Attaching package: ‘vegan’

The following object is masked from ‘package:MCMCpack’:

    procrustes

analogue version 0.17-0
Loading required package: rioja
This is rioja 0.9-15

Attaching package: ‘rioja’

The following objects are masked from ‘package:analogue’:

    crossval, performance

# figures/Fig_9_stratplot.png

lens_differences_raw_plot <- lens_differences_raw()

Attaching package: ‘dplyr’

The following object is masked from ‘package:MASS’:

    select

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

NAs introduced by coercionNAs introduced by coercionUsing Ages.kya, Axe.Flake.Frag, Ochre.Stained.Grindstone, part as id variables
Saving 11.3 x 4.36 in image

raw_fish_pval <- formatC(round(lens_differences_raw_plot$raw_chi_fisher$p.value, 4), big.mark=",",format="f", drop0trailing = TRUE)

lens_differences_tech_plot <- lens_differences_tech()

invalid factor level, NA generatedUsing Spit, part as id variables
Using Spit, part as id variables
Using Spit, part as id variables
Using percentage as value column: use value.var to override.
Saving 11.3 x 4.36 in image

tech_fish_pval <- round(lens_differences_tech_plot$tech_chi_fisher$p.value, 3)

lens_tech_and_raw_plot <- combine_lens_plots(lens_differences_tech_plot$p1,lens_differences_raw_plot$p1)
lens_tech_and_raw_plot

null device 1

cat("<p><img src='figures/Fig_10_combine_lens_plots.png', width='100%'></p>")

Chi square tests indicate that the differences in raw material proportions between all three assemblages are highly significant (p = 0.0005), although the differences in technological composition are not significant (p = 0.605), most likely due to small sample size (n = 10 for types in the assemblage below the lens).

Shellfish

shell_data_plotted <- plot_shell_data()

Using Spit, Taxon as id variables
Saving 7 x 7 in image

shell_data_plotted

R Session Information

sessionInfo()

R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2                      dplyr_0.7.4                      
 [3] rioja_0.9-15                      analogue_0.17-0                  
 [5] vegan_2.4-4                       lattice_0.20-35                  
 [7] permute_0.9-4                     RColorBrewer_1.1-2               
 [9] bcp_4.0.0                         BEST_0.5.0                       
[11] HDInterval_0.1.3                  MCMCpack_1.4-0                   
[13] MASS_7.3-47                       coda_0.19-1                      
[15] gdtools_0.1.6                     reshape2_1.4.2                   
[17] scales_0.5.0                      ggplot2_2.2.1                    
[19] mjb1989excavationpaper_0.0.0.9000

loaded via a namespace (and not attached):
 [1] colorspace_1.3-2   htmltools_0.3.6    yaml_2.1.14        mgcv_1.8-20       
 [5] base64enc_0.1-3    rlang_0.1.4        glue_1.2.0         withr_1.0.2       
 [9] bindr_0.1          plyr_1.8.4         stringr_1.2.0      MatrixModels_0.4-1
[13] munsell_0.4.3      gtable_0.2.0       devtools_1.13.4    memoise_1.1.0     
[17] evaluate_0.10.1    labeling_0.3       knitr_1.17         SparseM_1.77      
[21] quantreg_5.34      parallel_3.4.2     curl_3.0           highr_0.6         
[25] Rcpp_0.12.13       backports_1.1.0    jsonlite_1.5       princurve_1.1-12  
[29] gridExtra_2.3      mcmc_0.9-5         profileModel_0.5-9 rjags_4-6         
[33] digest_0.6.12      stringi_1.1.5      brglm_0.6.1        rprojroot_1.2     
[37] tools_3.4.2        magrittr_1.5       lazyeval_0.2.0     tibble_1.3.4      
[41] cluster_2.0.6      pkgconfig_2.0.1    Matrix_1.2-11      rsconnect_0.8.5   
[45] assertthat_0.2.0   rmarkdown_1.8      svglite_1.2.1      httr_1.2.1        
[49] R6_2.2.2           nlme_3.1-131       git2r_0.19.0       compiler_3.4.2