###Søndergaard and Villumsen. Radiocarbon-based models of earliest Iron Age village community. Supplementary materials SM1: Rcodes###

##ADD LIBRARIES
library(rcarbon)
library(spatstat)
library(sf)

##SET WORKING DIRECTORY AND READ DATA

setwd("...")#adjust according to directory
kbaydates <- read.csv('kbaydata.csv')

##General model parameters
runm <- 100 #100 years running mean smoothing of SPDs
binh <- 25 #bin clustering at site level
nsim <- 1000 #number of simulations
ncolours <- 100 #stkde colours

##Figure 3: SPDs all dates and shortlived dates only##
#Calibrate, bin, create SPDs
kbay_all.caldates <- calibrate(x=kbaydates$C14BP,errors=kbaydates$C14SD,calCurves='intcal20',normalised=FALSE)
kbay_all.bins <- binPrep(sites=kbaydates$SiteID,ages=kbaydates$C14BP,h=binh)
kbay_all.spd <- spd(kbay_all.caldates,timeRange=c(3000,900),bins=kbay_all.bins)
kbay_short <-kbaydates[grep(paste(c("grain","seed charred","plant charred","plant uncharred"),collapse="|"),kbaydates$SampleMaterial),] #subset shortlived samples
kbay_short.caldates <- calibrate(x=kbay_short$C14BP,errors=kbay_short$C14SD,normalised=FALSE)
kbay_short.bins <- binPrep(sites=kbay_short$SiteID,ages=kbay_short$C14BP,h=binh)
kbay_short.spd <- spd(kbay_short.caldates,timeRange=c(3000,900),bins=kbay_short.bins)

##Assess binning values for all dates dataset Supplement Figure S1
binsense(x=kbay_all.caldates,y=kbaydates$SiteID,h=seq(0,200,25),timeRange=c(3000,900),calendar="BCAD")

#Plot SPDs, also using rolling averages for smoothing)
plot(kbay_all.spd,fill.p="#696969",calendar="BCAD",xlim=c(-1000,1000))
plot(kbay_short.spd,fill.p="#A9A9A9",add=TRUE,calendar="BCAD")
plot(kbay_all.spd,runm=runm,type="simple",add=TRUE,col="#000080",lwd=2.5,lty=2,calendar="BCAD")
plot(kbay_short.spd,runm=runm,type="simple",add=TRUE,col="#1E90FF",lwd=2.5,lty=2,calendar="BCAD")
legend("topleft",bty="n",col=c("#696969","#A9A9A9","#000080","#1E90FF"),legend=c("all samples (n=5401)","shortlived samples subset (n=4096)","all samples 100 yrs running mean","shortlived samples subset 100 yrs running mean"),lwd=c(5,5,5,5),cex=0.85)

#Subset of kbay_all for timerange 2800BP-2000BP
kbaysub_all_2800_2000 <- subset(kbaydates,C14BP<=2800&C14BP>=2000)
#Calibrate, bin, create SPD of subset all 2800BP-2000BP
kbaysub_all_2800_2000.caldates <- calibrate(x=kbaysub_all_2800_2000$C14BP,errors=kbaysub_all_2800_2000$C14SD,normalised=FALSE)
kbaysub_all_2800_2000.bins <- binPrep(sites=kbaysub_all_2800_2000$SiteID,ages=kbaysub_all_2800_2000$C14BP,h=binh)
kbaysub_all_2800_2000.spd <- spd(kbaysub_all_2800_2000.caldates,timeRange=c(2800,2000),bins=kbaysub_all_2800_2000.bins)

#Subset of kbay_short for timerange 2800BP-2000BP
kbaysub_short_2800_2000 <- subset(kbay_short,C14BP<=2800&C14BP>=2000)
#Calibrate, bin, create SPD of subset shortlived 2800BP-2000BP
kbaysub_short_2800_2000.caldates <- calibrate(x=kbaysub_short_2800_2000$C14BP,errors=kbaysub_short_2800_2000$C14SD,normalised=FALSE)
kbaysub_short_2800_2000.bins <- binPrep(sites=kbaysub_short_2800_2000$SiteID,ages=kbaysub_short_2800_2000$C14BP,h=binh)
kbaysub_short_2800_2000.spd <- spd(kbaysub_short_2800_2000.caldates,timeRange=c(2800,2000),bins=kbaysub_short_2800_2000.bins)

##Calculate Pearson correlation values for Table 2 
#Pearson correlation coefficient test of kbay_short subset against kbay_all dataset
cor.test(kbay_all.spd$grid$PrDens,kbay_short.spd$grid$PrDens,method="pearson")
#Pearson correlation coefficient test of kbay_short subset against kbay_all dataset for timeframe 2800BP-2000BP
cor.test(kbaysub_all_2800_2000.spd$grid$PrDens,kbaysub_short_2800_2000.spd$grid$PrDens,method="pearson")

##Figure 6: cKDE and SPD with model test##
#Sample dates and create cKDE for timerange 2800-900 BP
kbay.randates <- sampleDates(kbay_all.caldates,bins=kbay_all.bins,nsim=nsim,boot=TRUE)
kbay.ckde <- ckde(kbay.randates,timeRange=c(2800,900),bw=50,normalised=FALSE)

#Generate smoothed SPD for timerange 2800-900 BP
kbay_all.spd.smoothed=spd(kbay_all.caldates,timeRange=c(2800,900),bins=kbay_all.bins,runm=runm)

#Start values adjusted to observed SPD
logFit <- nls(PrDens~SSlogis(calBP,Asym,xmid,scale),data=kbay_all.spd.smoothed$grid,control=nls.control(maxiter=200),start=list(Asym=0.2,xmid=median(seq(2800,900,-1)),scale=-100))

#Generate data frame containing fitted values
logFitDens=data.frame(calBP=kbay_all.spd.smoothed$grid$calBP,PrDens=SSlogis(input=kbay_all.spd.smoothed$grid$calBP,Asym=coefficients(logFit)[1],xmid=coefficients(logFit)[2],scal=coefficients(logFit)[3]))

#Use modelTest function (returning raw simulation output)
LogNull <-modelTest(kbay_all.caldates,errors=kbaydates$C14SD,bins=kbay_all.bins,nsim=nsim,timeRange=c(2800,900),model="custom",predgrid=logFitDens,runm=runm,raw=TRUE)

#Plot results
plot(kbay.ckde,type='envelope',fill.col="#87CEFA",calendar="BCAD",xlim=c(-800,1000))
title("a",adj=0.05,line=-2)

plot(LogNull,lwd.obs=1.9,calendar="BCAD",xlim=c(-800,1000))
title("b",adj=0.05,line=-2)
legend("bottomleft",bty="y",bg="white",col=c("black","lightgrey",rgb(0.7,0,0,0.2),rgb(0,0,0.7,0.2)),legend=c("SPD","MC 95% envelope","positive deviation","negative deviation"),lwd=c(1.9,5,5,5),cex=0.85)

#Plot rate of change
plot(LogNull,type='roc',lwd.obs=1.9,calendar="BCAD",xlim=c(-800,1000))
title("c",adj=0.05,line=-2)
legend("bottomleft",bty="y",bg="white",col=c("black","lightgrey",rgb(0.7,0,0,0.2),rgb(0,0,0.7,0.2)),legend=c("ROC","MC 95% envelope","positive deviation","negative deviation"),lwd=c(1.9,5,5,5),cex=0.85)
#plot xlim=c(2749,950) for BP axis

#Global p-value and summary
LogNull$pval
#0.000999001
summary(LogNull)
#'modelTest()' function summary:

#Number of radiocarbon dates: 5401
#Number of bins: 1143

#Statistical Significance computed using 1000 simulations. 
#Global p-value: 0.001.

##Figure 7: stkde##

#Create clas owin object from study area shp
kb <- st_read("kb.shp")
kb <- as.owin(kb)

focalBC <- seq(550,0,-50) # BC years
focalBP <- focalBC + 1950 # convert to cal BP

#Subset of all samples 2550BP-1950BP
kbaysub_all_2550_1950 <- subset(kbaydates,C14BP<=2550&C14BP>=1950)

#Calibrate timerange subset, bin and create 50 yrs timeslices 
x <- calibrate(x=kbaysub_all_2550_1950$C14BP,errors=kbaysub_all_2550_1950$C14SD,normalised=FALSE)
binskbsub_all_2550_1950 <- binPrep(sites=kbaysub_all_2550_1950$SiteID,ages=kbaysub_all_2550_1950$C14BP,h=binh)

stkdekbsub <- stkde(x=x,coords=kbaysub_all_2550_1950[,c("Eastings","Northings")],win=kb,sbw=500,cellres=10,focalyears=focalBP,tbw=30,bins=binskbsub_all_2550_1950,backsight=50,outdir=tempdir(),amount=1)

#Loop over BC labels while plotting the corresponding BP slices:

par(mfrow=c(6,4),mar=c(0.5,0.5,2.5,0.5))

for(bc in focalBC) {bp <- bc + 1950

plot(stkdekbsub,bp,type="focal",ramptype="std",ncolours=ncolours, main="")
title(paste0(bc," BC 
Focal Intensity"), cex.main = 0.8)

plot(stkdekbsub,bp,type="change",ramptype="std",ncolours=ncolours, main="")
title(paste0(bc," BC 
Focal Change (50 yrs backsight)"), cex.main = 0.8)}

#Plot supplementary Figure S3.2 in SM3
par(mfrow=c(6,4),mar=c(0.5,0.5,2.5,0.5))

for(bc in focalBC) {bp <- bc + 1950

plot(stkdekbsub,bp,type="focal",ramptype="std",ncolours=ncolours, main="")
title(paste0(bc," BC 
Focal Intensity"), cex.main = 0.8)

plot(stkdekbsub,bp,type="proportion",ramptype="std",ncolours=ncolours, main="")
title(paste0(bc," BC 
Focal Proportion"), cex.main = 0.8)}