# Code for SPDs in Solheim "Timing the emergence and development of arable farming in Southeastern Norway by using summed probability distribution of radiocarbon dates and a Bayesian age model"

#Summed probability distributions are modelled using Rcarbon package for R.
#Reference: Crema, E., Bevan,  A. 2020. Inference from large sets of radiocarbon dates: Software and methods. Radiocarbon. DOI:10.1017/RDC.2020.95

#Further resources and information:	
#https://cran.r-project.org/web/packages/rcarbon/vignettes/rcarbon.html
#https://github.com/ahb108/rcarbon


## Add libraries (you may need to install these first)

library(rcarbon)
library(rgdal)
library(drc)

# Load data
mydates <- read.csv("cereals_senorway.csv", header=TRUE, stringsAsFactors=FALSE, encoding="UTF-8",na.strings=c("NA",""),strip.white=TRUE)

############################# SPD radiocarbon dates ####################

# General SPD parameters
nsim <- 1000 # number of actual simulations 
ncores <- 4 # multi-core processing (set higher if available)
runm <- 100 #smoothing of SPDsnr
bins <- 100 #bin clustering
calendar <- 'BP'
#calendar <- 'BCAD'

### SPD ###
caldates=calibrate(x=mydates$nrca,errors=mydates$se,calCurves='intcal20',ncores=3)
spd = spd(caldates,timeRange=c(5700,700)) 
plot(spd, calendar=calendar)
plot(spd, runm=runm,calendar=calendar, add=TRUE,type="simple",col="indianred",lwd=2,lty=2) #using a rolling average of ## years for smoothing

### Binning ###
bins = binPrep(sites=mydates$site_id,ages=mydates$nrca,h=bins)
spd.bins = spd(caldates,bins=bins,timeRange=c(5700,700))
plot(spd.bins, calendar=calendar)

### Visualising bins ###
bins.med=binMed(x = caldates,bins=bins)
plot(spd.bins,runm=runm)
barCodes(bins.med,yrng = c(0,0.010))


### Testing against theorethical growth models ###
#Exponential model
expnull <- modelTest(caldates, errors=mydates$se, bins=bins, nsim=nsim, timeRange=c(5700,700), model="exponential",runm=runm, datenormalised=FALSE)
plot(expnull, calendar=calendar)
expnull$pval
summary(expnull)