library(rcarbon)
alldates=read.csv("c14_burials.csv")
weighted=read.csv("c14_burials_weighted.csv")

#Marine curve construction
Curve1.2=mixCurves('intcal20', 'marine20', p=0.988, resOffsets=-47, resErrors=52)
Curve2.4=mixCurves('intcal20','marine20',p=0.976,resOffsets=-47,resErrors=52)
Curve3.5=mixCurves('intcal20','marine20',p=0.965,resOffsets=-47,resErrors=52)
Curve4.7=mixCurves('intcal20','marine20',p=0.953,resOffsets=-47,resErrors=52)
Curve5.9=mixCurves('intcal20','marine20',p=0.941,resOffsets=-47,resErrors=52)
Curve6.4=mixCurves('intcal20','marine20',p=0.936,resOffsets=-47,resErrors=52)
Curve6.7=mixCurves('intcal20','marine20',p=0.933,resOffsets=-47,resErrors=52)
Curve7.1=mixCurves('intcal20','marine20',p=0.929,resOffsets=-47,resErrors=52)
Curve8=mixCurves('intcal20','marine20',p=0.92,resOffsets=-47,resErrors=52)
Curve8.2=mixCurves('intcal20','marine20',p=0.918,resOffsets=-47,resErrors=52)
Curve8.3=mixCurves('intcal20','marine20',p=0.917,resOffsets=-47,resErrors=52)
Curve8.7=mixCurves('intcal20','marine20',p=0.913,resOffsets=-47,resErrors=52)
Curve9.4=mixCurves('intcal20','marine20',p=0.906,resOffsets=-47,resErrors=52)
Curve9.5=mixCurves('intcal20','marine20',p=0.905,resOffsets=-47,resErrors=52)
Curve9.6=mixCurves('intcal20','marine20',p=0.904,resOffsets=-47,resErrors=52)
Curve10=mixCurves('intcal20','marine20',p=0.9,resOffsets=-47,resErrors=52)
Curve10.1=mixCurves('intcal20','marine20',p=0.899,resOffsets=-47,resErrors=52)
Curve10.3=mixCurves('intcal20','marine20',p=0.897,resOffsets=-47,resErrors=52)
Curve10.5=mixCurves('intcal20','marine20',p=0.895,resOffsets=-47,resErrors=52)
Curve10.6=mixCurves('intcal20','marine20',p=0.894,resOffsets=-47,resErrors=52)
Curve11.5=mixCurves('intcal20','marine20',p=0.885,resOffsets=-47,resErrors=52)
Curve11.8=mixCurves('intcal20','marine20',p=0.882,resOffsets=-47,resErrors=52)
Curve12.2=mixCurves('intcal20','marine20',p=0.878,resOffsets=-47,resErrors=52)
Curve12.5=mixCurves('intcal20','marine20',p=0.875,resOffsets=-47,resErrors=52)
Curve12.9=mixCurves('intcal20','marine20',p=0.871,resOffsets=-47,resErrors=52)
Curve13.9=mixCurves('intcal20','marine20',p=0.861,resOffsets=-47,resErrors=52)
Curve14.1=mixCurves('intcal20','marine20',p=0.859,resOffsets=-47,resErrors=52)
Curve14.4=mixCurves('intcal20','marine20',p=0.856,resOffsets=-47,resErrors=52)
Curve14.6=mixCurves('intcal20','marine20',p=0.854,resOffsets=-47,resErrors=52)
Curve15.3=mixCurves('intcal20','marine20',p=0.847,resOffsets=-47,resErrors=52)
Curve16.5=mixCurves('intcal20','marine20',p=0.835,resOffsets=-47,resErrors=52)
Curve17.6=mixCurves('intcal20','marine20',p=0.824,resOffsets=-47,resErrors=0)
Curve17.8=mixCurves('intcal20','marine20',p=0.822,resOffsets=-47,resErrors=0)
Curve18.6=mixCurves('intcal20','marine20',p=0.814,resOffsets=-47,resErrors=52)
Curve18.8=mixCurves('intcal20','marine20',p=0.812,resOffsets=-47,resErrors=52)
Curve18.9=mixCurves('intcal20','marine20',p=0.811,resOffsets=-47,resErrors=52)
Curve20=mixCurves('intcal20','marine20',p=0.8,resOffsets=-47,resErrors=52)
Curve20.4=mixCurves('intcal20','marine20',p=0.796,resOffsets=-47,resErrors=52)
Curve21.1=mixCurves('intcal20','marine20',p=0.789,resOffsets=-47,resErrors=52)
Curve21.2=mixCurves('intcal20','marine20',p=0.788,resOffsets=-47,resErrors=52)
Curve21.9=mixCurves('intcal20','marine20',p=0.781,resOffsets=-47,resErrors=52)
Curve22.4=mixCurves('intcal20','marine20',p=0.776,resOffsets=-47,resErrors=52)
Curve23.5=mixCurves('intcal20','marine20',p=0.765,resOffsets=-47,resErrors=52)
Curve24.7=mixCurves('intcal20','marine20',p=0.753,resOffsets=-47,resErrors=52)
Curve25.9=mixCurves('intcal20','marine20',p=0.741,resOffsets=-47,resErrors=52)
Curve27.1=mixCurves('intcal20','marine20',p=0.729,resOffsets=-47,resErrors=52)
Curve28.2=mixCurves('intcal20','marine20',p=0.718,resOffsets=-47,resErrors=52)
Curve28.4=mixCurves('intcal20','marine20',p=0.716,resOffsets=-47,resErrors=52)
Curve29.4=mixCurves('intcal20','marine20',p=0.706,resOffsets=-47,resErrors=52)
Curve30.6=mixCurves('intcal20','marine20',p=0.694,resOffsets=-47,resErrors=52)
Curve31.8=mixCurves('intcal20','marine20',p=0.682,resOffsets=-47,resErrors=52)
Curve32=mixCurves('intcal20','marine20',p=0.68,resOffsets=-47,resErrors=52)
Curve32.7=mixCurves('intcal20','marine20',p=0.673,resOffsets=-47,resErrors=52)
Curve32.9=mixCurves('intcal20','marine20',p=0.671,resOffsets=-47,resErrors=52)
Curve34.1=mixCurves('intcal20','marine20',p=0.659,resOffsets=-47,resErrors=52)
Curve37.6=mixCurves('intcal20','marine20',p=0.624,resOffsets=-47,resErrors=52)
Curve38.8=mixCurves('intcal20','marine20',p=0.612,resOffsets=-4752,resErrors=52)

#All dated graves
dates1.2=subset(alldates, Marine.protein==1.2)
cal1.2=calibrate(x=dates1.2$CRA, errors=dates1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
dates2.4=subset(alldates, Marine.protein==2.4)
cal2.4=calibrate(x=dates2.4$CRA, errors=dates2.4$Error, normalised=FALSE, calCurves = Curve2.4)
dates3.5=subset(alldates, Marine.protein==3.5)
cal3.5=calibrate(x=dates3.5$CRA, errors=dates3.5$Error, normalised=FALSE, calCurves = Curve3.5)
dates4.7=subset(alldates, Marine.protein==4.7)
cal4.7=calibrate(x=dates4.7$CRA, errors=dates4.7$Error, normalised=FALSE, calCurves = Curve4.7)
dates5.9=subset(alldates, Marine.protein==5.9)
cal5.9=calibrate(x=dates5.9$CRA, errors=dates5.9$Error, normalised=FALSE, calCurves = Curve5.9)
dates6.4=subset(alldates, Marine.protein==6.4)
cal6.4=calibrate(x=dates6.4$CRA, errors=dates6.4$Error, normalised=FALSE, calCurves = Curve6.4)
dates6.7=subset(alldates, Marine.protein==6.7)
cal6.7=calibrate(x=dates6.7$CRA, errors=dates6.7$Error, normalised=FALSE, calCurves = Curve6.7)
dates7.1=subset(alldates, Marine.protein==7.1)
cal7.1=calibrate(x=dates7.1$CRA, errors=dates7.1$Error, normalised=FALSE, calCurves = Curve7.1)
dates8=subset(alldates, Marine.protein==8)
cal8=calibrate(x=dates8$CRA, errors=dates8$Error, normalised=FALSE, calCurves = Curve8)
dates8.2=subset(alldates, Marine.protein==8.2)
cal8.2=calibrate(x=dates8.2$CRA, errors=dates8.2$Error, normalised=FALSE, calCurves = Curve8.2)
dates8.3=subset(alldates, Marine.protein==8.3)
cal8.3=calibrate(x=dates8.3$CRA, errors=dates8.3$Error, normalised=FALSE, calCurves = Curve8.3)
dates8.7=subset(alldates, Marine.protein==8.7)
cal8.7=calibrate(x=dates8.7$CRA, errors=dates8.7$Error, normalised=FALSE, calCurves = Curve8.7)
dates9.4=subset(alldates, Marine.protein==9.4)
cal9.4=calibrate(x=dates9.4$CRA, errors=dates9.4$Error, normalised=FALSE, calCurves = Curve9.4)
dates9.5=subset(alldates, Marine.protein==9.5)
cal9.5=calibrate(x=dates9.5$CRA, errors=dates9.5$Error, normalised=FALSE, calCurves = Curve9.5)
dates9.6=subset(alldates, Marine.protein==9.6)
cal9.6=calibrate(x=dates9.6$CRA, errors=dates9.6$Error, normalised=FALSE, calCurves = Curve9.6)
dates10=subset(alldates, Marine.protein==10)
cal10=calibrate(x=dates10$CRA, errors=dates10$Error, normalised=FALSE, calCurves = Curve10)
dates10.1=subset(alldates, Marine.protein==10.1)
cal10.1=calibrate(x=dates10.1$CRA, errors=dates10.1$Error, normalised=FALSE, calCurves = Curve10.1)
dates10.3=subset(alldates, Marine.protein==10.3)
cal10.3=calibrate(x=dates10.3$CRA, errors=dates10.3$Error, normalised=FALSE, calCurves = Curve10.3)
dates10.5=subset(alldates, Marine.protein==10.5)
cal10.5=calibrate(x=dates10.5$CRA, errors=dates10.5$Error, normalised=FALSE, calCurves = Curve10.5)
dates10.6=subset(alldates, Marine.protein==10.6)
cal10.6=calibrate(x=dates10.6$CRA, errors=dates10.6$Error, normalised=FALSE, calCurves = Curve10.6)
dates11.5=subset(alldates, Marine.protein==11.5)
cal11.5=calibrate(x=dates11.5$CRA, errors=dates11.5$Error, normalised=FALSE, calCurves = Curve11.5)
dates11.8=subset(alldates, Marine.protein==11.8)
cal11.8=calibrate(x=dates11.8$CRA, errors=dates11.8$Error, normalised=FALSE, calCurves = Curve11.8)
dates12.2=subset(alldates, Marine.protein==12.2)
cal12.2=calibrate(x=dates12.2$CRA, errors=dates12.2$Error, normalised=FALSE, calCurves = Curve12.2)
dates12.5=subset(alldates, Marine.protein==12.5)
cal12.5=calibrate(x=dates12.5$CRA, errors=dates12.5$Error, normalised=FALSE, calCurves = Curve12.5)
dates12.9=subset(alldates, Marine.protein==12.9)
cal12.9=calibrate(x=dates12.9$CRA, errors=dates12.9$Error, normalised=FALSE, calCurves = Curve12.9)
dates13.9=subset(alldates, Marine.protein==13.9)
cal13.9=calibrate(x=dates13.9$CRA, errors=dates13.9$Error, normalised=FALSE, calCurves = Curve13.9)
dates14.1=subset(alldates, Marine.protein==14.1)
cal14.1=calibrate(x=dates14.1$CRA, errors=dates14.1$Error, normalised=FALSE, calCurves = Curve14.1)
dates14.6=subset(alldates, Marine.protein==14.6)
cal14.6=calibrate(x=dates14.6$CRA, errors=dates14.6$Error, normalised=FALSE, calCurves = Curve14.6)
dates15.3=subset(alldates, Marine.protein==15.3)
cal15.3=calibrate(x=dates15.3$CRA, errors=dates15.3$Error, normalised=FALSE, calCurves = Curve15.3)
dates16.5=subset(alldates, Marine.protein==16.5)
cal16.5=calibrate(x=dates16.5$CRA, errors=dates16.5$Error, normalised=FALSE, calCurves = Curve16.5)
dates17.6=subset(alldates, Marine.protein==17.6)
cal17.6=calibrate(x=dates17.6$CRA, errors=dates17.6$Error, normalised=FALSE, calCurves = Curve17.6)
dates17.8=subset(alldates, Marine.protein==17.8)
cal17.8=calibrate(x=dates17.8$CRA, errors=dates17.8$Error, normalised=FALSE, calCurves = Curve17.8)
dates18.6=subset(alldates, Marine.protein==18.6)
cal18.6=calibrate(x=dates18.6$CRA, errors=dates18.6$Error, normalised=FALSE, calCurves = Curve18.6)
dates18.8=subset(alldates, Marine.protein==18.8)
cal18.8=calibrate(x=dates18.8$CRA, errors=dates18.8$Error, normalised=FALSE, calCurves = Curve18.8)
dates18.9=subset(alldates, Marine.protein==18.9)
cal18.9=calibrate(x=dates18.9$CRA, errors=dates18.9$Error, normalised=FALSE, calCurves = Curve18.9)
dates20=subset(alldates, Marine.protein==20)
cal20=calibrate(x=dates20$CRA, errors=dates20$Error, normalised=FALSE, calCurves = Curve20)
dates20.4=subset(alldates, Marine.protein==20.4)
cal20.4=calibrate(x=dates20.4$CRA, errors=dates20.4$Error, normalised=FALSE, calCurves = Curve20.4)
dates21.1=subset(alldates, Marine.protein==21.1)
cal21.1=calibrate(x=dates21.1$CRA, errors=dates21.1$Error, normalised=FALSE, calCurves = Curve21.1)
dates21.2=subset(alldates, Marine.protein==21.2)
cal21.2=calibrate(x=dates21.2$CRA, errors=dates21.2$Error, normalised=FALSE, calCurves = Curve21.2)
dates21.9=subset(alldates, Marine.protein==21.9)
cal21.9=calibrate(x=dates21.9$CRA, errors=dates21.9$Error, normalised=FALSE, calCurves = Curve21.9)
dates22.4=subset(alldates, Marine.protein==22.4)
cal22.4=calibrate(x=dates22.4$CRA, errors=dates22.4$Error, normalised=FALSE, calCurves = Curve22.4)
dates23.5=subset(alldates, Marine.protein==23.5)
cal23.5=calibrate(x=dates23.5$CRA, errors=dates23.5$Error, normalised=FALSE, calCurves = Curve23.5)
dates24.7=subset(alldates, Marine.protein==24.7)
cal24.7=calibrate(x=dates24.7$CRA, errors=dates24.7$Error, normalised=FALSE, calCurves = Curve24.7)
dates25.9=subset(alldates, Marine.protein==25.9)
cal25.9=calibrate(x=dates25.9$CRA, errors=dates25.9$Error, normalised=FALSE, calCurves = Curve25.9)
dates27.1=subset(alldates, Marine.protein==27.1)
cal27.1=calibrate(x=dates27.1$CRA, errors=dates27.1$Error, normalised=FALSE, calCurves = Curve27.1)
dates28.2=subset(alldates, Marine.protein==28.2)
cal28.2=calibrate(x=dates28.2$CRA, errors=dates28.2$Error, normalised=FALSE, calCurves = Curve28.2)
dates28.4=subset(alldates, Marine.protein==28.4)
cal28.4=calibrate(x=dates28.4$CRA, errors=dates28.4$Error, normalised=FALSE, calCurves = Curve28.4)
dates29.4=subset(alldates, Marine.protein==29.4)
cal29.4=calibrate(x=dates29.4$CRA, errors=dates29.4$Error, normalised=FALSE, calCurves = Curve29.4)
dates30.6=subset(alldates, Marine.protein==30.6)
cal30.6=calibrate(x=dates30.6$CRA, errors=dates30.6$Error, normalised=FALSE, calCurves = Curve30.6)
dates31.8=subset(alldates, Marine.protein==31.8)
cal31.8=calibrate(x=dates31.8$CRA, errors=dates31.8$Error, normalised=FALSE, calCurves = Curve31.8)
dates32=subset(alldates, Marine.protein==32)
cal32=calibrate(x=dates32$CRA, errors=dates32$Error, normalised=FALSE, calCurves = Curve32)
dates32.7=subset(alldates, Marine.protein==32.7)
cal32.7=calibrate(x=dates32.7$CRA, errors=dates32.7$Error, normalised=FALSE, calCurves = Curve32.7)
dates32.9=subset(alldates, Marine.protein==32.9)
cal32.9=calibrate(x=dates32.9$CRA, errors=dates32.9$Error, normalised=FALSE, calCurves = Curve32.9)
dates34.1=subset(alldates, Marine.protein==34.1)
cal34.1=calibrate(x=dates34.1$CRA, errors=dates34.1$Error, normalised=FALSE, calCurves = Curve34.1)
dates37.6=subset(alldates, Marine.protein==37.6)
cal37.6=calibrate(x=dates37.6$CRA, errors=dates37.6$Error, normalised=FALSE, calCurves = Curve37.6)
dates38.8=subset(alldates, Marine.protein==38.8)
cal38.8=calibrate(x=dates38.8$CRA, errors=dates38.8$Error, normalised=FALSE, calCurves = Curve38.8)
terrestrialall=subset(alldates, is.na(alldates$Marine.protein)| Marine.protein=="0")
calterrestrial=calibrate(x=terrestrialall$CRA, errors=terrestrialall$Error, normalised=FALSE)
combdate=combine(cal1.2, cal2.4, cal3.5, cal4.7, cal5.9, cal6.4, cal6.7, cal7.1, cal8, cal8.2, cal8.3, cal8.7, cal9.4, cal9.5, cal9.6, cal10, cal10.1, cal10.3, cal10.5, cal10.6, cal11.5, cal11.8, cal12.2, cal12.5, cal12.9, cal13.9, cal14.1, cal14.6, cal15.3, cal16.5, cal17.6, cal17.8, cal18.6, cal18.8, cal18.9, cal20, cal20.4, cal21.1, cal21.2, cal21.9, cal22.4, cal23.5, cal24.7, cal25.9, cal27.1, cal28.2, cal28.4, cal29.4, cal30.6, cal31.8,cal32, cal32.7, cal32.9, cal34.1, cal37.6, cal38.8, calterrestrial, fixIDs = TRUE)
bins=binPrep(sites=alldates$SiteName, ages=alldates$CRA, h=50)
samplemarineall=sampleDates(combdate, bins=bins, nsim=1000, boot=TRUE)
kdemarineall30=ckde(samplemarineall, timeRange = c(1650,1050), bw=30)
plot(kdemarineall30, calendar="BCAD") +title("All radiocarbon-dated graves, n=1085")

#Furnished graves
furnished=subset(alldates, Furnished.=="Furnished")
furnished1.2=subset(furnished, Marine.protein==1.2)
furnishedcal1.2=calibrate(x=furnished1.2$CRA, errors=furnished1.2$Error, normalised=FALSE, calCurves=Curve1.2)
furnished2.4=subset(furnished, Marine.protein==2.4)
furnishedcal2.4=calibrate(x=furnished2.4$CRA, errors=furnished2.4$Error, normalised=FALSE, calCurves = Curve2.4)
furnished3.5=subset(furnished, Marine.protein==3.5)
furnishedcal3.5=calibrate(x=furnished3.5$CRA, errors=furnished3.5$Error, normalised=FALSE, calCurves = Curve3.5)
furnished4.7=subset(furnished, Marine.protein==4.7)
furnishedcal4.7=calibrate(x=furnished4.7$CRA, errors=furnished4.7$Error, normalised=FALSE, calCurves = Curve4.7)
furnished5.9=subset(furnished, Marine.protein==5.9)
furnishedcal5.9=calibrate(x=furnished5.9$CRA, errors=furnished5.9$Error, normalised=FALSE, calCurves = Curve5.9)
furnished6.7=subset(furnished, Marine.protein==6.7)
furnishedcal6.7=calibrate(x=furnished6.7$CRA, errors=furnished6.7$Error, normalised=FALSE, calCurves = Curve6.7)
furnished7.1=subset(furnished, Marine.protein==7.1)
furnishedcal7.1=calibrate(x=furnished7.1$CRA, errors=furnished7.1$Error, normalised=FALSE, calCurves = Curve7.1)
furnished8=subset(furnished, Marine.protein==8)
furnishedcal8=calibrate(x=furnished8$CRA, errors=furnished8$Error, normalised=FALSE, calCurves = Curve8)
furnished8.2=subset(furnished, Marine.protein==8.2)
furnishedcal8.2=calibrate(x=furnished8.2$CRA, errors=furnished8.2$Error, normalised=FALSE, calCurves = Curve8.2)
furnished8.3=subset(furnished, Marine.protein==8.3)
furnishedcal8.3=calibrate(x=furnished8.3$CRA, errors=furnished8.3$Error, normalised=FALSE, calCurves = Curve8.3)
furnished8.7=subset(furnished, Marine.protein==8.7)
furnishedcal8.7=calibrate(x=furnished8.7$CRA, errors=furnished8.7$Error, normalised=FALSE, calCurves = Curve8.7)
furnished9.4=subset(furnished, Marine.protein==9.4)
furnishedcal9.4=calibrate(x=furnished9.4$CRA, errors=furnished9.4$Error, normalised=FALSE, calCurves = Curve9.4)
furnished9.5=subset(furnished, Marine.protein==9.5)
furnishedcal9.5=calibrate(x=furnished9.5$CRA, errors=furnished9.5$Error, normalised=FALSE, calCurves = Curve9.5)
furnished9.6=subset(furnished, Marine.protein==9.6)
furnishedcal9.6=calibrate(x=furnished9.6$CRA, errors=furnished9.6$Error, normalised=FALSE, calCurves = Curve9.6)
furnished10.1=subset(furnished, Marine.protein==10.1)
furnishedcal10.1=calibrate(x=furnished10.1$CRA, errors=furnished10.1$Error, normalised=FALSE, calCurves = Curve10.1)
furnished10.6=subset(furnished, Marine.protein==10.6)
furnishedcal10.6=calibrate(x=furnished10.6$CRA, errors=furnished10.6$Error, normalised=FALSE, calCurves = Curve10.6)
furnished11.8=subset(furnished, Marine.protein==11.8)
furnishedcal11.8=calibrate(x=furnished11.8$CRA, errors=furnished11.8$Error, normalised=FALSE, calCurves = Curve11.8)
furnished12.9=subset(furnished, Marine.protein==12.9)
furnishedcal12.9=calibrate(x=furnished12.9$CRA, errors=furnished12.9$Error, normalised=FALSE, calCurves = Curve12.9)
furnished13.9=subset(furnished, Marine.protein==13.9)
furnishedcal13.9=calibrate(x=furnished13.9$CRA, errors=furnished13.9$Error, normalised=FALSE, calCurves = Curve13.9)
furnished14.1=subset(furnished, Marine.protein==14.1)
furnishedcal14.1=calibrate(x=furnished14.1$CRA, errors=furnished14.1$Error, normalised=FALSE, calCurves = Curve14.1)
furnished14.6=subset(furnished, Marine.protein==14.6)
furnishedcal14.6=calibrate(x=furnished14.6$CRA, errors=furnished14.6$Error, normalised=FALSE, calCurves = Curve14.6)
furnished15.3=subset(furnished, Marine.protein==15.3)
furnishedcal15.3=calibrate(x=furnished15.3$CRA, errors=furnished15.3$Error, normalised=FALSE, calCurves = Curve15.3)
furnished16.5=subset(furnished, Marine.protein==16.5)
furnishedcal16.5=calibrate(x=furnished16.5$CRA, errors=furnished16.5$Error, normalised=FALSE, calCurves = Curve16.5)
furnished17.6=subset(furnished, Marine.protein==17.6)
furnishedcal17.6=calibrate(x=furnished17.6$CRA, errors=furnished17.6$Error, normalised=FALSE, calCurves = Curve17.6)
furnished18.8=subset(furnished, Marine.protein==18.8)
furnishedcal18.8=calibrate(x=furnished18.8$CRA, errors=furnished18.8$Error, normalised=FALSE, calCurves = Curve18.8)
furnished20=subset(furnished, Marine.protein==20)
furnishedcal20=calibrate(x=furnished20$CRA, errors=furnished20$Error, normalised=FALSE, calCurves = Curve20)
furnished21.2=subset(furnished, Marine.protein==21.2)
furnishedcal21.2=calibrate(x=furnished21.2$CRA, errors=furnished21.2$Error, normalised=FALSE, calCurves = Curve21.2)
furnished22.4=subset(furnished, Marine.protein==22.4)
furnishedcal22.4=calibrate(x=furnished22.4$CRA, errors=furnished22.4$Error, normalised=FALSE, calCurves = Curve22.4)
furnished23.5=subset(furnished, Marine.protein==23.5)
furnishedcal23.5=calibrate(x=furnished23.5$CRA, errors=furnished23.5$Error, normalised=FALSE, calCurves = Curve23.5)
furnished25.9=subset(furnished, Marine.protein==25.9)
furnishedcal25.9=calibrate(x=furnished25.9$CRA, errors=furnished25.9$Error, normalised=FALSE, calCurves = Curve25.9)
furnished27.1=subset(furnished, Marine.protein==27.1)
furnishedcal27.1=calibrate(x=furnished27.1$CRA, errors=furnished27.1$Error, normalised=FALSE, calCurves = Curve27.1)
furnished28.2=subset(furnished, Marine.protein==28.2)
furnishedcal28.2=calibrate(x=furnished28.2$CRA, errors=furnished28.2$Error, normalised=FALSE, calCurves = Curve28.2)
furnished28.4=subset(furnished, Marine.protein==28.4)
furnishedcal28.4=calibrate(x=furnished28.4$CRA, errors=furnished28.4$Error, normalised=FALSE, calCurves = Curve28.4)
furnished32.7=subset(furnished, Marine.protein==32.7)
furnishedcal32.7=calibrate(x=furnished32.7$CRA, errors=furnished32.7$Error, normalised=FALSE, calCurves = Curve32.7)
furnished32.9=subset(furnished, Marine.protein==32.9)
furnishedcal32.9=calibrate(x=furnished32.9$CRA, errors=furnished32.9$Error, normalised=FALSE, calCurves = Curve32.9)
furnished37.6=subset(furnished, Marine.protein==37.6)
furnishedcal37.6=calibrate(x=furnished37.6$CRA, errors=furnished37.6$Error, normalised=FALSE, calCurves = Curve37.6)
furnterrestrial=subset(furnished, is.na(furnished$Marine.protein)| Marine.protein=="0")
furnishedcalterrestrial=calibrate(x=furnterrestrial$CRA, errors=furnterrestrial$Error, normalised=FALSE)
combfurnisheddate=combine(furnishedcal1.2, furnishedcal2.4, furnishedcal3.5, furnishedcal4.7, furnishedcal5.9, furnishedcal6.7, furnishedcal7.1, furnishedcal8, furnishedcal8.2, furnishedcal8.3, furnishedcal8.7, furnishedcal9.4, furnishedcal9.5, furnishedcal9.6, furnishedcal10.1, furnishedcal10.6, furnishedcal11.8, furnishedcal12.9, furnishedcal13.9, furnishedcal14.1, furnishedcal14.6, furnishedcal15.3, furnishedcal16.5, furnishedcal17.6, furnishedcal18.8, furnishedcal20, furnishedcal21.2, furnishedcal22.4, furnishedcal23.5, furnishedcal25.9, furnishedcal27.1, furnishedcal28.2, furnishedcal28.4, furnishedcal32.7, furnishedcal32.9, furnishedcal37.6, furnishedcalterrestrial, fixIDs = TRUE)
binsf=binPrep(sites = furnished$SiteName, ages=furnished$CRA, h=50)
samplefurnmarine=sampleDates(combfurnisheddate, bins=binsf, nsim=1000, boot=TRUE)
kdemarinefurnished30=ckde(samplefurnmarine, timeRange=c(1650,1050), bw=30)

#Furnished graves weighted by number of objects
weighted1.2=subset(weighted, Marine.protein==1.2)
weightedcal1.2=calibrate(x=weighted1.2$CRA, errors=weighted1.2$Error, normalised=FALSE, calCurves=Curve1.2)
weighted2.4=subset(weighted, Marine.protein==2.4)
weightedcal2.4=calibrate(x=weighted2.4$CRA, errors=weighted2.4$Error, normalised=FALSE, calCurves = Curve2.4)
weighted3.5=subset(weighted, Marine.protein==3.5)
weightedcal3.5=calibrate(x=weighted3.5$CRA, errors=weighted3.5$Error, normalised=FALSE, calCurves = Curve3.5)
weighted4.7=subset(weighted, Marine.protein==4.7)
weightedcal4.7=calibrate(x=weighted4.7$CRA, errors=weighted4.7$Error, normalised=FALSE, calCurves = Curve4.7)
weighted5.9=subset(weighted, Marine.protein==5.9)
weightedcal5.9=calibrate(x=weighted5.9$CRA, errors=weighted5.9$Error, normalised=FALSE, calCurves = Curve5.9)
weighted6.7=subset(weighted, Marine.protein==6.7)
weightedcal6.7=calibrate(x=weighted6.7$CRA, errors=weighted6.7$Error, normalised=FALSE, calCurves = Curve6.7)
weighted7.1=subset(weighted, Marine.protein==7.1)
weightedcal7.1=calibrate(x=weighted7.1$CRA, errors=weighted7.1$Error, normalised=FALSE, calCurves = Curve7.1)
weighted8=subset(weighted, Marine.protein==8)
weightedcal8=calibrate(x=weighted8$CRA, errors=weighted8$Error, normalised=FALSE, calCurves = Curve8)
weighted8.2=subset(weighted, Marine.protein==8.2)
weightedcal8.2=calibrate(x=weighted8.2$CRA, errors=weighted8.2$Error, normalised=FALSE, calCurves = Curve8.2)
weighted8.3=subset(weighted, Marine.protein==8.3)
weightedcal8.3=calibrate(x=weighted8.3$CRA, errors=weighted8.3$Error, normalised=FALSE, calCurves = Curve8.3)
weighted8.7=subset(weighted, Marine.protein==8.7)
weightedcal8.7=calibrate(x=weighted8.7$CRA, errors=weighted8.7$Error, normalised=FALSE, calCurves = Curve8.7)
weighted9.4=subset(weighted, Marine.protein==9.4)
weightedcal9.4=calibrate(x=weighted9.4$CRA, errors=weighted9.4$Error, normalised=FALSE, calCurves = Curve9.4)
weighted9.5=subset(weighted, Marine.protein==9.5)
weightedcal9.5=calibrate(x=weighted9.5$CRA, errors=weighted9.5$Error, normalised=FALSE, calCurves = Curve9.5)
weighted9.6=subset(weighted, Marine.protein==9.6)
weightedcal9.6=calibrate(x=weighted9.6$CRA, errors=weighted9.6$Error, normalised=FALSE, calCurves = Curve9.6)
weighted10.1=subset(weighted, Marine.protein==10.1)
weightedcal10.1=calibrate(x=weighted10.1$CRA, errors=weighted10.1$Error, normalised=FALSE, calCurves = Curve10.1)
weighted10.6=subset(weighted, Marine.protein==10.6)
weightedcal10.6=calibrate(x=weighted10.6$CRA, errors=weighted10.6$Error, normalised=FALSE, calCurves = Curve10.6)
weighted11.8=subset(weighted, Marine.protein==11.8)
weightedcal11.8=calibrate(x=weighted11.8$CRA, errors=weighted11.8$Error, normalised=FALSE, calCurves = Curve11.8)
weighted12.9=subset(weighted, Marine.protein==12.9)
weightedcal12.9=calibrate(x=weighted12.9$CRA, errors=weighted12.9$Error, normalised=FALSE, calCurves = Curve12.9)
weighted13.9=subset(weighted, Marine.protein==13.9)
weightedcal13.9=calibrate(x=weighted13.9$CRA, errors=weighted13.9$Error, normalised=FALSE, calCurves = Curve13.9)
weighted14.1=subset(weighted, Marine.protein==14.1)
weightedcal14.1=calibrate(x=weighted14.1$CRA, errors=weighted14.1$Error, normalised=FALSE, calCurves = Curve14.1)
weighted14.6=subset(weighted, Marine.protein==14.6)
weightedcal14.6=calibrate(x=weighted14.6$CRA, errors=weighted14.6$Error, normalised=FALSE, calCurves = Curve14.6)
weighted15.3=subset(weighted, Marine.protein==15.3)
weightedcal15.3=calibrate(x=weighted15.3$CRA, errors=weighted15.3$Error, normalised=FALSE, calCurves = Curve15.3)
weighted16.5=subset(weighted, Marine.protein==16.5)
weightedcal16.5=calibrate(x=weighted16.5$CRA, errors=weighted16.5$Error, normalised=FALSE, calCurves = Curve16.5)
weighted17.6=subset(weighted, Marine.protein==17.6)
weightedcal17.6=calibrate(x=weighted17.6$CRA, errors=weighted17.6$Error, normalised=FALSE, calCurves = Curve17.6)
weighted18.8=subset(weighted, Marine.protein==18.8)
weightedcal18.8=calibrate(x=weighted18.8$CRA, errors=weighted18.8$Error, normalised=FALSE, calCurves = Curve18.8)
weighted20=subset(weighted, Marine.protein==20)
weightedcal20=calibrate(x=weighted20$CRA, errors=weighted20$Error, normalised=FALSE, calCurves = Curve20)
weighted21.2=subset(weighted, Marine.protein==21.2)
weightedcal21.2=calibrate(x=weighted21.2$CRA, errors=weighted21.2$Error, normalised=FALSE, calCurves = Curve21.2)
weighted22.4=subset(weighted, Marine.protein==22.4)
weightedcal22.4=calibrate(x=weighted22.4$CRA, errors=weighted22.4$Error, normalised=FALSE, calCurves = Curve22.4)
weighted23.5=subset(weighted, Marine.protein==23.5)
weightedcal23.5=calibrate(x=weighted23.5$CRA, errors=weighted23.5$Error, normalised=FALSE, calCurves = Curve23.5)
weighted25.9=subset(weighted, Marine.protein==25.9)
weightedcal25.9=calibrate(x=weighted25.9$CRA, errors=weighted25.9$Error, normalised=FALSE, calCurves = Curve25.9)
weighted27.1=subset(weighted, Marine.protein==27.1)
weightedcal27.1=calibrate(x=weighted27.1$CRA, errors=weighted27.1$Error, normalised=FALSE, calCurves = Curve27.1)
weighted28.2=subset(weighted, Marine.protein==28.2)
weightedcal28.2=calibrate(x=weighted28.2$CRA, errors=weighted28.2$Error, normalised=FALSE, calCurves = Curve28.2)
weighted28.4=subset(weighted, Marine.protein==28.4)
weightedcal28.4=calibrate(x=weighted28.4$CRA, errors=weighted28.4$Error, normalised=FALSE, calCurves = Curve28.4)
weighted32.7=subset(weighted, Marine.protein==32.7)
weightedcal32.7=calibrate(x=weighted32.7$CRA, errors=weighted32.7$Error, normalised=FALSE, calCurves = Curve32.7)
weighted32.9=subset(weighted, Marine.protein==32.9)
weightedcal32.9=calibrate(x=weighted32.9$CRA, errors=weighted32.9$Error, normalised=FALSE, calCurves = Curve32.9)
weighted37.6=subset(weighted, Marine.protein==37.6)
weightedcal37.6=calibrate(x=weighted37.6$CRA, errors=weighted37.6$Error, normalised=FALSE, calCurves = Curve37.6)
weightedterrestrial=subset(weighted, is.na(weighted$Marine.protein)| Marine.protein=="0")
weightedcalterrestrial=calibrate(x=weightedterrestrial$CRA, errors=weightedterrestrial$Error, normalised=FALSE)
combweighteddate=combine(weightedcal1.2, weightedcal2.4, weightedcal3.5, weightedcal4.7, weightedcal5.9, weightedcal6.7, weightedcal7.1, weightedcal8, weightedcal8.2, weightedcal8.3, weightedcal8.7, weightedcal9.4, weightedcal9.5, weightedcal9.6, weightedcal10.1, weightedcal10.6, weightedcal11.8, weightedcal12.9, weightedcal13.9, weightedcal14.1, weightedcal14.6, weightedcal15.3, weightedcal16.5, weightedcal17.6, weightedcal18.8, weightedcal20, weightedcal21.2, weightedcal22.4, weightedcal23.5, weightedcal25.9, weightedcal27.1, weightedcal28.2, weightedcal28.4, weightedcal32.7, weightedcal32.9, weightedcal37.6, weightedcalterrestrial, fixIDs = TRUE)
sampleweightedmarine=sampleDates(combweighteddate, nsim=1000, boot=TRUE)
kdemarineweighted30=ckde(sampleweightedmarine, timeRange=c(1650,1050), bw=30)

#Plot figure 4a
plot(kdemarinefurnished30, calendar="BCAD", ylim=c(0,0.004)) + title("Furnished graves, n=353")
par(new=TRUE)
plot(kdemarineweighted30, calendar="BCAD", ylim=c(0,0.004), fill.col="#FFFF007D")
legend("topright", legend=c("Furnished graves","Weighted by number of grave goods"), fill=c("lightgrey", "#FFFF007D"))

#Unfurnished graves
unfurnished=subset(alldates, Furnished.=="Unfurnished")
unfurnished1.2=subset(unfurnished, Marine.protein==1.2)
unfurnishedcal1.2=calibrate(x=unfurnished1.2$CRA, errors=unfurnished1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
unfurnished2.4=subset(unfurnished, Marine.protein==2.4)
unfurnishedcal2.4=calibrate(x=unfurnished2.4$CRA, errors=unfurnished2.4$Error, normalised=FALSE, calCurves = Curve2.4)
unfurnished3.5=subset(unfurnished, Marine.protein==3.5)
unfurnishedcal3.5=calibrate(x=unfurnished3.5$CRA, errors=unfurnished3.5$Error, normalised=FALSE, calCurves = Curve3.5)
unfurnished4.7=subset(unfurnished, Marine.protein==4.7)
unfurnishedcal4.7=calibrate(x=unfurnished4.7$CRA, errors=unfurnished4.7$Error, normalised=FALSE, calCurves = Curve4.7)
unfurnished5.9=subset(unfurnished, Marine.protein==5.9)
unfurnishedcal5.9=calibrate(x=unfurnished5.9$CRA, errors=unfurnished5.9$Error, normalised=FALSE, calCurves = Curve5.9)
unfurnished6.4=subset(unfurnished, Marine.protein==6.4)
unfurnishedcal6.4=calibrate(x=unfurnished6.4$CRA, errors=unfurnished6.4$Error, normalised=FALSE, calCurves = Curve6.4)
unfurnished7.1=subset(unfurnished, Marine.protein==7.1)
unfurnishedcal7.1=calibrate(x=unfurnished7.1$CRA, errors=unfurnished7.1$Error, normalised=FALSE, calCurves = Curve7.1)
unfurnished8.2=subset(unfurnished, Marine.protein==8.2)
unfurnishedcal8.2=calibrate(x=unfurnished8.2$CRA, errors=unfurnished8.2$Error, normalised=FALSE, calCurves = Curve8.2)
unfurnished9.4=subset(unfurnished, Marine.protein==9.4)
unfurnishedcal9.4=calibrate(x=unfurnished9.4$CRA, errors=unfurnished9.4$Error, normalised=FALSE, calCurves = Curve9.4)
unfurnished10=subset(unfurnished, Marine.protein==10)
unfurnishedcal10=calibrate(x=unfurnished10$CRA, errors=unfurnished10$Error, normalised=FALSE, calCurves = Curve10)
unfurnished10.1=subset(unfurnished, Marine.protein==10.1)
unfurnishedcal10.1=calibrate(x=unfurnished10.1$CRA, errors=unfurnished10.1$Error, normalised=FALSE, calCurves = Curve10.1)
unfurnished10.5=subset(unfurnished, Marine.protein==10.5)
unfurnishedcal10.5=calibrate(x=unfurnished10.5$CRA, errors=unfurnished10.5$Error, normalised=FALSE, calCurves = Curve10.5)
unfurnished10.6=subset(unfurnished, Marine.protein==10.6)
unfurnishedcal10.6=calibrate(x=unfurnished10.6$CRA, errors=unfurnished10.6$Error, normalised=FALSE, calCurves = Curve10.6)
unfurnished11.5=subset(unfurnished, Marine.protein==11.5)
unfurnishedcal11.5=calibrate(x=unfurnished11.5$CRA, errors=unfurnished11.5$Error, normalised=FALSE, calCurves=Curve11.5)
unfurnished11.8=subset(unfurnished, Marine.protein==11.8)
unfurnishedcal11.8=calibrate(x=unfurnished11.8$CRA, errors=unfurnished11.8$Error, normalised=FALSE, calCurves = Curve11.8)
unfurnished12.2=subset(unfurnished, Marine.protein==12.2)
unfurnishedcal12.2=calibrate(x=unfurnished12.2$CRA, errors=unfurnished12.2$Error, normalised=FALSE, calCurves = Curve12.2)
unfurnished12.5=subset(unfurnished, Marine.protein==12.5)
unfurnishedcal12.5=calibrate(x=unfurnished12.5$CRA, errors=unfurnished12.5$Error, normalised=FALSE, calCurves = Curve12.5)
unfurnished12.9=subset(unfurnished, Marine.protein==12.9)
unfurnishedcal12.9=calibrate(x=unfurnished12.9$CRA, errors=unfurnished12.9$Error, normalised=FALSE, calCurves = Curve12.9)
unfurnished14.1=subset(unfurnished, Marine.protein==14.1)
unfurnishedcal14.1=calibrate(x=unfurnished14.1$CRA, errors=unfurnished14.1$Error, normalised=FALSE, calCurves = Curve14.1)
unfurnished15.3=subset(unfurnished, Marine.protein==15.3)
unfurnishedcal15.3=calibrate(x=unfurnished15.3$CRA, errors=unfurnished15.3$Error, normalised=FALSE, calCurves = Curve15.3)
unfurnished16.5=subset(unfurnished, Marine.protein==16.5)
unfurnishedcal16.5=calibrate(x=unfurnished16.5$CRA, errors=unfurnished16.5$Error, normalised=FALSE, calCurves = Curve16.5)
unfurnished17.6=subset(unfurnished, Marine.protein==17.6)
unfurnishedcal17.6=calibrate(x=unfurnished17.6$CRA, errors=unfurnished17.6$Error, normalised=FALSE, calCurves = Curve17.6)
unfurnished17.8=subset(unfurnished, Marine.protein==17.8)
unfurnishedcal17.8=calibrate(x=unfurnished17.8$CRA, errors=unfurnished17.8$Error, normalised=FALSE, calCurves = Curve17.8)
unfurnished18.6=subset(unfurnished, Marine.protein==18.6)
unfurnishedcal18.6=calibrate(x=unfurnished18.6$CRA, errors=unfurnished18.6$Error, normalised=FALSE, calCurves = Curve18.6)
unfurnished18.8=subset(unfurnished, Marine.protein==18.8)
unfurnishedcal18.8=calibrate(x=unfurnished18.8$CRA, errors=unfurnished18.8$Error, normalised=FALSE, calCurves = Curve18.8)
unfurnished18.9=subset(unfurnished, Marine.protein==18.9)
unfurnishedcal18.9=calibrate(x=unfurnished18.9$CRA, errors=unfurnished18.9$Error, normalised=FALSE, calCurves = Curve18.9)
unfurnished20=subset(unfurnished, Marine.protein==20)
unfurnishedcal20=calibrate(x=unfurnished20$CRA, errors=unfurnished20$Error, normalised=FALSE, calCurves = Curve20)
unfurnished20.4=subset(unfurnished, Marine.protein==20.4)
unfurnishedcal20.4=calibrate(x=unfurnished20.4$CRA, errors=unfurnished20.4$Error, normalised=FALSE, calCurves = Curve20.4)
unfurnished21.1=subset(unfurnished, Marine.protein==21.1)
unfurnishedcal21.1=calibrate(x=unfurnished21.1$CRA, errors=unfurnished21.1$Error, normalised=FALSE, calCurves = Curve21.1)
unfurnished21.2=subset(unfurnished, Marine.protein==21.2)
unfurnishedcal21.2=calibrate(x=unfurnished21.2$CRA, errors=unfurnished21.2$Error, normalised=FALSE, calCurves = Curve21.2)
unfurnished21.9=subset(unfurnished, Marine.protein==21.9)
unfurnishedcal21.9=calibrate(x=unfurnished21.9$CRA, errors=unfurnished21.9$Error, normalised=FALSE, calCurves = Curve21.9)
unfurnished22.4=subset(unfurnished, Marine.protein==22.4)
unfurnishedcal22.4=calibrate(x=unfurnished22.4$CRA, errors=unfurnished22.4$Error, normalised=FALSE, calCurves = Curve22.4)
unfurnished23.5=subset(unfurnished, Marine.protein==23.5)
unfurnishedcal23.5=calibrate(x=unfurnished23.5$CRA, errors=unfurnished23.5$Error, normalised=FALSE, calCurves = Curve23.5)
unfurnished24.7=subset(unfurnished, Marine.protein==24.7)
unfurnishedcal24.7=calibrate(x=unfurnished24.7$CRA, errors=unfurnished24.7$Error, normalised=FALSE, calCurves = Curve24.7)
unfurnished25.9=subset(unfurnished, Marine.protein==25.9)
unfurnishedcal25.9=calibrate(x=unfurnished25.9$CRA, errors=unfurnished25.9$Error, normalised=FALSE, calCurves = Curve25.9)
unfurnished27.1=subset(unfurnished, Marine.protein==27.1)
unfurnishedcal27.1=calibrate(x=unfurnished27.1$CRA, errors=unfurnished27.1$Error, normalised=FALSE, calCurves = Curve27.1)
unfurnished28.2=subset(unfurnished, Marine.protein==28.2)
unfurnishedcal28.2=calibrate(x=unfurnished28.2$CRA, errors=unfurnished28.2$Error, normalised=FALSE, calCurves = Curve28.2)
unfurnished29.4=subset(unfurnished, Marine.protein==29.4)
unfurnishedcal29.4=calibrate(x=unfurnished29.4$CRA, errors=unfurnished29.4$Error, normalised=FALSE, calCurves = Curve29.4)
unfurnished30.6=subset(unfurnished, Marine.protein==30.6)
unfurnishedcal30.6=calibrate(x=unfurnished30.6$CRA, errors=unfurnished30.6$Error, normalised=FALSE, calCurves = Curve30.6)
unfurnished31.8=subset(unfurnished, Marine.protein==31.8)
unfurnishedcal31.8=calibrate(x=unfurnished31.8$CRA, errors=unfurnished31.8$Error, normalised=FALSE, calCurves = Curve31.8)
unfurnished32=subset(unfurnished, Marine.protein==32)
unfurnishedcal32=calibrate(x=unfurnished32$CRA, errors=unfurnished32$Error, normalised=FALSE, calCurves = Curve32)
unfurnished32.9=subset(unfurnished, Marine.protein==32.9)
unfurnishedcal32.9=calibrate(x=unfurnished32.9$CRA, errors=unfurnished32.9$Error, normalised=FALSE, calCurves = Curve32.9)
unfurnished34.1=subset(unfurnished, Marine.protein==34.1)
unfurnishedcal34.1=calibrate(x=unfurnished34.1$CRA, errors=unfurnished34.1$Error, normalised=FALSE, calCurves = Curve34.1)
unfurnished37.6=subset(unfurnished, Marine.protein==37.6)
unfurnishedcal37.6=calibrate(x=unfurnished37.6$CRA, errors=unfurnished37.6$Error, normalised=FALSE, calCurves = Curve37.6)
unfurnished38.8=subset(unfurnished, Marine.protein==38.8)
unfurnishedcal38.8=calibrate(x=unfurnished38.8$CRA, errors=unfurnished38.8$Error, normalised=FALSE, calCurves = Curve38.8)
unfterrestrial=subset(unfurnished, is.na(unfurnished$Marine.protein)| Marine.protein=="0")
unfurnishedcalterrestrial=calibrate(x=unfterrestrial$CRA, errors=unfterrestrial$Error, normalised=FALSE)
combdateunf=combine(unfurnishedcal1.2, unfurnishedcal2.4, unfurnishedcal3.5, unfurnishedcal4.7, unfurnishedcal5.9, unfurnishedcal6.4, unfurnishedcal7.1, unfurnishedcal8.2, unfurnishedcal9.4, unfurnishedcal10, unfurnishedcal10.1, unfurnishedcal10.5, unfurnishedcal10.6, unfurnishedcal11.5, unfurnishedcal11.8, unfurnishedcal12.2, unfurnishedcal12.5, unfurnishedcal12.9, unfurnishedcal14.1, unfurnishedcal15.3, unfurnishedcal16.5, unfurnishedcal17.6, unfurnishedcal17.8, unfurnishedcal18.6, unfurnishedcal18.8, unfurnishedcal18.9, unfurnishedcal20, unfurnishedcal20.4, unfurnishedcal21.1, unfurnishedcal21.2, unfurnishedcal21.9, unfurnishedcal22.4, unfurnishedcal23.5, unfurnishedcal24.7, unfurnishedcal25.9, unfurnishedcal27.1, unfurnishedcal28.2, unfurnishedcal29.4, unfurnishedcal30.6, unfurnishedcal31.8,unfurnishedcal32, unfurnishedcal32.9, unfurnishedcal34.1, unfurnishedcal37.6, unfurnishedcal38.8, unfurnishedcalterrestrial, fixIDs = TRUE)
binsu=binPrep(sites=unfurnished$SiteName, ages=unfurnished$CRA, h=50)
sampleunfurnishedmarine=sampleDates(combdateunf, bins=binsu, nsim=1000, boot=TRUE)
kdeunfurnished30marine=ckde(sampleunfurnishedmarine, timeRange=c(1650,1050), bw=30)
plot(kdeunfurnished30marine, calendar="BCAD", ylim=c(0,0.002)) + title("Unfurnished graves, n=715")

#Female furnished graves
femfurnished=subset(furnished, Sex=="F" | Sex=="F?")
femfurdates1.2=subset(femfurnished, Marine.protein==1.2)
calfemfurnished1.2=calibrate(x=femfurdates1.2$CRA, errors=femfurdates1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
femfurdates2.4=subset(femfurnished, Marine.protein==2.4)
calfemfurnished2.4=calibrate(x=femfurdates2.4$CRA, errors=femfurdates2.4$Error, normalised=FALSE, calCurves = Curve2.4)
femfurdates3.5=subset(femfurnished, Marine.protein==3.5)
calfemfurnished3.5=calibrate(x=femfurdates3.5$CRA, errors=femfurdates3.5$Error, normalised=FALSE, calCurves = Curve3.5)
femfurdates4.7=subset(femfurnished, Marine.protein==4.7)
calfemfurnished4.7=calibrate(x=femfurdates4.7$CRA, errors=femfurdates4.7$Error, normalised=FALSE, calCurves = Curve4.7)
femfurdates5.9=subset(femfurnished, Marine.protein==5.9)
calfemfurnished5.9=calibrate(x=femfurdates5.9$CRA, errors=femfurdates5.9$Error, normalised=FALSE, calCurves = Curve5.9)
femfurdates7.1=subset(femfurnished, Marine.protein==7.1)
calfemfurnished7.1=calibrate(x=femfurdates7.1$CRA, errors=femfurdates7.1$Error, normalised=FALSE, calCurves = Curve7.1)
femfurdates8=subset(femfurnished, Marine.protein==8)
calfemfurnished8=calibrate(x=femfurdates8$CRA, errors=femfurdates8$Error, normalised=FALSE, calCurves = Curve8)
femfurdates8.2=subset(femfurnished, Marine.protein==8.2)
calfemfurnished8.2=calibrate(x=femfurdates8.2$CRA, errors=femfurdates8.2$Error, normalised=FALSE, calCurves = Curve8.2)
femfurdates8.3=subset(femfurnished, Marine.protein==8.3)
calfemfurnished8.3=calibrate(x=femfurdates8.3$CRA, errors=femfurdates8.3$Error, normalised=FALSE, calCurves = Curve8.3)
femfurdates9.4=subset(femfurnished, Marine.protein==9.4)
calfemfurnished9.4=calibrate(x=femfurdates9.4$CRA, errors=femfurdates9.4$Error, normalised=FALSE, calCurves = Curve9.4)
femfurdates9.6=subset(femfurnished, Marine.protein==9.6)
calfemfurnished9.6=calibrate(x=femfurdates9.6$CRA, errors=femfurdates9.6$Error, normalised=FALSE, calCurves = Curve9.6)
femfurdates10.1=subset(femfurnished, Marine.protein==10.1)
calfemfurnished10.1=calibrate(x=femfurdates10.1$CRA, errors=femfurdates10.1$Error, normalised=FALSE, calCurves = Curve10.1)
femfurdates10.6=subset(femfurnished, Marine.protein==10.6)
calfemfurnished10.6=calibrate(x=femfurdates10.6$CRA, errors=femfurdates10.6$Error, normalised=FALSE, calCurves = Curve10.6)
femfurdates11.8=subset(femfurnished, Marine.protein==11.8)
calfemfurnished11.8=calibrate(x=femfurdates11.8$CRA, errors=femfurdates11.8$Error, normalised=FALSE, calCurves = Curve11.8)
femfurdates12.9=subset(femfurnished, Marine.protein==12.9)
calfemfurnished12.9=calibrate(x=femfurdates12.9$CRA, errors=femfurdates12.9$Error, normalised=FALSE, calCurves = Curve12.9)
femfurdates13.9=subset(femfurnished, Marine.protein==13.9)
calfemfurnished13.9=calibrate(x=femfurdates13.9$CRA, errors=femfurdates13.9$Error, normalised=FALSE, calCurves = Curve13.9)
femfurdates14.1=subset(femfurnished, Marine.protein==14.1)
calfemfurnished14.1=calibrate(x=femfurdates14.1$CRA, errors=femfurdates14.1$Error, normalised=FALSE, calCurves = Curve14.1)
femfurdates15.3=subset(femfurnished, Marine.protein==15.3)
calfemfurnished15.3=calibrate(x=femfurdates15.3$CRA, errors=femfurdates15.3$Error, normalised=FALSE, calCurves = Curve15.3)
femfurdates16.5=subset(femfurnished, Marine.protein==16.5)
calfemfurnished16.5=calibrate(x=femfurdates16.5$CRA, errors=femfurdates16.5$Error, normalised=FALSE, calCurves = Curve16.5)
femfurdates17.6=subset(femfurnished, Marine.protein==17.6)
calfemfurnished17.6=calibrate(x=femfurdates17.6$CRA, errors=femfurdates17.6$Error, normalised=FALSE, calCurves = Curve17.6)
femfurdates21.2=subset(femfurnished, Marine.protein==21.2)
calfemfurnished21.2=calibrate(x=femfurdates21.2$CRA, errors=femfurdates21.2$Error, normalised=FALSE, calCurves = Curve21.2)
femfurdates22.4=subset(femfurnished, Marine.protein==22.4)
calfemfurnished22.4=calibrate(x=femfurdates22.4$CRA, errors=femfurdates22.4$Error, normalised=FALSE, calCurves = Curve22.4)
femfurdates25.9=subset(femfurnished, Marine.protein==25.9)
calfemfurnished25.9=calibrate(x=femfurdates25.9$CRA, errors=femfurdates25.9$Error, normalised=FALSE, calCurves = Curve25.9)
femfurdates27.1=subset(femfurnished, Marine.protein==27.1)
calfemfurnished27.1=calibrate(x=femfurdates27.1$CRA, errors=femfurdates27.1$Error, normalised=FALSE, calCurves = Curve27.1)
femfurdates32.9=subset(femfurnished, Marine.protein==32.9)
calfemfurnished32.9=calibrate(x=femfurdates32.9$CRA, errors=femfurdates32.9$Error, normalised=FALSE, calCurves = Curve32.9)
femfurdates37.6=subset(femfurnished, Marine.protein==37.6)
calfemfurnished37.6=calibrate(x=femfurdates37.6$CRA, errors=femfurdates37.6$Error, normalised=FALSE, calCurves = Curve37.6)
femfurterrestrial=subset(femfurnished, is.na(femfurnished$Marine.protein)| Marine.protein=="0")
calfemfurnishedterrestrial=calibrate(x=femfurterrestrial$CRA, errors=femfurterrestrial$Error, normalised=FALSE)
combdatefemfurnished=combine(calfemfurnished1.2, calfemfurnished2.4, calfemfurnished3.5, calfemfurnished4.7, calfemfurnished5.9, calfemfurnished7.1, calfemfurnished8, calfemfurnished8.2, calfemfurnished8.3, calfemfurnished9.4, calfemfurnished9.6, calfemfurnished10.1, calfemfurnished10.6, calfemfurnished11.8, calfemfurnished12.9, calfemfurnished13.9, calfemfurnished14.1, calfemfurnished15.3, calfemfurnished16.5, calfemfurnished17.6, calfemfurnished21.2, calfemfurnished22.4, calfemfurnished25.9, calfemfurnished27.1, calfemfurnished32.9, calfemfurnished37.6, calfemfurnishedterrestrial, fixIDs = TRUE)
binsff=binPrep(sites = femfurnished$SiteName, ages = femfurnished$CRA, h=50)
samplefemfurnmarine=sampleDates(combdatefemfurnished, bins=binsff, nsim=1000, boot=TRUE)
kdefemfurnished30marine=ckde(samplefemfurnmarine, timeRange=c(1650,1050), bw=30)

#Female furnished graves, weighted by number of objects
femweighted=subset(weighted, Sex=="F" | Sex=="F?")
femweightdates1.2=subset(femweighted, Marine.protein==1.2)
calfemweighted1.2=calibrate(x=femweightdates1.2$CRA, errors=femweightdates1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
femweightdates2.4=subset(femweighted, Marine.protein==2.4)
calfemweighted2.4=calibrate(x=femweightdates2.4$CRA, errors=femweightdates2.4$Error, normalised=FALSE, calCurves = Curve2.4)
femweightdates3.5=subset(femweighted, Marine.protein==3.5)
calfemweighted3.5=calibrate(x=femweightdates3.5$CRA, errors=femweightdates3.5$Error, normalised=FALSE, calCurves = Curve3.5)
femweightdates4.7=subset(femweighted, Marine.protein==4.7)
calfemweighted4.7=calibrate(x=femweightdates4.7$CRA, errors=femweightdates4.7$Error, normalised=FALSE, calCurves = Curve4.7)
femweightdates5.9=subset(femweighted, Marine.protein==5.9)
calfemweighted5.9=calibrate(x=femweightdates5.9$CRA, errors=femweightdates5.9$Error, normalised=FALSE, calCurves = Curve5.9)
femweightdates7.1=subset(femweighted, Marine.protein==7.1)
calfemweighted7.1=calibrate(x=femweightdates7.1$CRA, errors=femweightdates7.1$Error, normalised=FALSE, calCurves = Curve7.1)
femweightdates8=subset(femweighted, Marine.protein==8)
calfemweighted8=calibrate(x=femweightdates8$CRA, errors=femweightdates8$Error, normalised=FALSE, calCurves = Curve8)
femweightdates8.2=subset(femweighted, Marine.protein==8.2)
calfemweighted8.2=calibrate(x=femweightdates8.2$CRA, errors=femweightdates8.2$Error, normalised=FALSE, calCurves = Curve8.2)
femweightdates8.3=subset(femweighted, Marine.protein==8.3)
calfemweighted8.3=calibrate(x=femweightdates8.3$CRA, errors=femweightdates8.3$Error, normalised=FALSE, calCurves = Curve8.3)
femweightdates9.4=subset(femweighted, Marine.protein==9.4)
calfemweighted9.4=calibrate(x=femweightdates9.4$CRA, errors=femweightdates9.4$Error, normalised=FALSE, calCurves = Curve9.4)
femweightdates9.6=subset(femweighted, Marine.protein==9.6)
calfemweighted9.6=calibrate(x=femweightdates9.6$CRA, errors=femweightdates9.6$Error, normalised=FALSE, calCurves = Curve9.6)
femweightdates10.1=subset(femweighted, Marine.protein==10.1)
calfemweighted10.1=calibrate(x=femweightdates10.1$CRA, errors=femweightdates10.1$Error, normalised=FALSE, calCurves = Curve10.1)
femweightdates10.6=subset(femweighted, Marine.protein==10.6)
calfemweighted10.6=calibrate(x=femweightdates10.6$CRA, errors=femweightdates10.6$Error, normalised=FALSE, calCurves = Curve10.6)
femweightdates11.8=subset(femweighted, Marine.protein==11.8)
calfemweighted11.8=calibrate(x=femweightdates11.8$CRA, errors=femweightdates11.8$Error, normalised=FALSE, calCurves = Curve11.8)
femweightdates12.9=subset(femweighted, Marine.protein==12.9)
calfemweighted12.9=calibrate(x=femweightdates12.9$CRA, errors=femweightdates12.9$Error, normalised=FALSE, calCurves = Curve12.9)
femweightdates13.9=subset(femweighted, Marine.protein==13.9)
calfemweighted13.9=calibrate(x=femweightdates13.9$CRA, errors=femweightdates13.9$Error, normalised=FALSE, calCurves = Curve13.9)
femweightdates14.1=subset(femweighted, Marine.protein==14.1)
calfemweighted14.1=calibrate(x=femweightdates14.1$CRA, errors=femweightdates14.1$Error, normalised=FALSE, calCurves = Curve14.1)
femweightdates15.3=subset(femweighted, Marine.protein==15.3)
calfemweighted15.3=calibrate(x=femweightdates15.3$CRA, errors=femweightdates15.3$Error, normalised=FALSE, calCurves = Curve15.3)
femweightdates16.5=subset(femweighted, Marine.protein==16.5)
calfemweighted16.5=calibrate(x=femweightdates16.5$CRA, errors=femweightdates16.5$Error, normalised=FALSE, calCurves = Curve16.5)
femweightdates17.6=subset(femweighted, Marine.protein==17.6)
calfemweighted17.6=calibrate(x=femweightdates17.6$CRA, errors=femweightdates17.6$Error, normalised=FALSE, calCurves = Curve17.6)
femweightdates21.2=subset(femweighted, Marine.protein==21.2)
calfemweighted21.2=calibrate(x=femweightdates21.2$CRA, errors=femweightdates21.2$Error, normalised=FALSE, calCurves = Curve21.2)
femweightdates22.4=subset(femweighted, Marine.protein==22.4)
calfemweighted22.4=calibrate(x=femweightdates22.4$CRA, errors=femweightdates22.4$Error, normalised=FALSE, calCurves = Curve22.4)
femweightdates25.9=subset(femweighted, Marine.protein==25.9)
calfemweighted25.9=calibrate(x=femweightdates25.9$CRA, errors=femweightdates25.9$Error, normalised=FALSE, calCurves = Curve25.9)
femweightdates27.1=subset(femweighted, Marine.protein==27.1)
calfemweighted27.1=calibrate(x=femweightdates27.1$CRA, errors=femweightdates27.1$Error, normalised=FALSE, calCurves = Curve27.1)
femweightdates32.9=subset(femweighted, Marine.protein==32.9)
calfemweighted32.9=calibrate(x=femweightdates32.9$CRA, errors=femweightdates32.9$Error, normalised=FALSE, calCurves = Curve32.9)
femweightdates37.6=subset(femweighted, Marine.protein==37.6)
calfemweighted37.6=calibrate(x=femweightdates37.6$CRA, errors=femweightdates37.6$Error, normalised=FALSE, calCurves = Curve37.6)
femweightterrestrial=subset(femweighted, is.na(femweighted$Marine.protein)| Marine.protein=="0")
calfemweightedterrestrial=calibrate(x=femweightterrestrial$CRA, errors=femweightterrestrial$Error, normalised=FALSE)
combdatefemweighted=combine(calfemweighted1.2, calfemweighted2.4, calfemweighted3.5, calfemweighted4.7, calfemweighted5.9, calfemweighted7.1, calfemweighted8, calfemweighted8.2, calfemweighted8.3, calfemweighted9.4, calfemweighted9.6, calfemweighted10.1, calfemweighted10.6, calfemweighted11.8, calfemweighted12.9, calfemweighted13.9, calfemweighted14.1, calfemweighted15.3, calfemweighted16.5, calfemweighted17.6, calfemweighted21.2, calfemweighted22.4, calfemweighted25.9, calfemweighted27.1, calfemweighted32.9, calfemweighted37.6, calfemweightedterrestrial, fixIDs = TRUE)
samplefemweightnmarine=sampleDates(combdatefemweighted, nsim=1000, boot=TRUE)
kdefemweighted30marine=ckde(samplefemweightnmarine, timeRange=c(1650,1050), bw=30)

#Plot figure 5a
plot(kdefemfurnished30marine, calendar="BCAD", ylim=c(0,0.004)) + title("Female furnished graves, n=162")
par(new=TRUE)
plot(kdefemweighted30marine, calendar="BCAD", ylim=c(0,0.004), fill.col="#FFFF007D")
legend("topright", legend=c("Female furnished graves","Weighted by number of grave goods"), fill=c("lightgrey", "#FFFF007D"))

#Female unfurnished graves
femunfurnished=subset(unfurnished, Sex=="F" | Sex=="F?")
femunfurdates2.4=subset(femunfurnished, Marine.protein==2.4)
calfemunfurnished2.4=calibrate(x=femunfurdates2.4$CRA, errors=femunfurdates2.4$Error, normalised=FALSE, calCurves = Curve2.4)
femunfurdates3.5=subset(femunfurnished, Marine.protein==3.5)
calfemunfurnished3.5=calibrate(x=femunfurdates3.5$CRA, errors=femunfurdates3.5$Error, normalised=FALSE, calCurves = Curve3.5)
femunfurdates4.7=subset(femunfurnished, Marine.protein==4.7)
calfemunfurnished4.7=calibrate(x=femunfurdates4.7$CRA, errors=femunfurdates4.7$Error, normalised=FALSE, calCurves = Curve4.7)
femunfurdates5.9=subset(femunfurnished, Marine.protein==5.9)
calfemunfurnished5.9=calibrate(x=femunfurdates5.9$CRA, errors=femunfurdates5.9$Error, normalised=FALSE, calCurves = Curve5.9)
femunfurdates6.4=subset(femunfurnished, Marine.protein==6.4)
calfemunfurnished6.4=calibrate(x=femunfurdates6.4$CRA, errors=femunfurdates6.4$Error, normalised=FALSE, calCurves = Curve6.4)
femunfurdates7.1=subset(femunfurnished, Marine.protein==7.1)
calfemunfurnished7.1=calibrate(x=femunfurdates7.1$CRA, errors=femunfurdates7.1$Error, normalised=FALSE, calCurves = Curve7.1)
femunfurdates8.2=subset(femunfurnished, Marine.protein==8.2)
calfemunfurnished8.2=calibrate(x=femunfurdates8.2$CRA, errors=femunfurdates8.2$Error, normalised=FALSE, calCurves = Curve8.2)
femunfurdates9.4=subset(femunfurnished, Marine.protein==9.4)
calfemunfurnished9.4=calibrate(x=femunfurdates9.4$CRA, errors=femunfurdates9.4$Error, normalised=FALSE, calCurves = Curve9.4)
femunfurdates10.5=subset(femunfurnished, Marine.protein==10.5)
calfemunfurnished10.5=calibrate(x=femunfurdates10.5$CRA, errors=femunfurdates10.5$Error, normalised=FALSE, calCurves = Curve10.5)
femunfurdates10.6=subset(femunfurnished, Marine.protein==10.6)
calfemunfurnished10.6=calibrate(x=femunfurdates10.6$CRA, errors=femunfurdates10.6$Error, normalised=FALSE, calCurves = Curve10.6)
femunfurdates11.8=subset(femunfurnished, Marine.protein==11.8)
calfemunfurnished11.8=calibrate(x=femunfurdates11.8$CRA, errors=femunfurdates11.8$Error, normalised=FALSE, calCurves = Curve11.8)
femunfurdates12.9=subset(femunfurnished, Marine.protein==12.9)
calfemunfurnished12.9=calibrate(x=femunfurdates12.9$CRA, errors=femunfurdates12.9$Error, normalised=FALSE, calCurves = Curve12.9)
femunfurdates14.1=subset(femunfurnished, Marine.protein==14.1)
calfemunfurnished14.1=calibrate(x=femunfurdates14.1$CRA, errors=femunfurdates14.1$Error, normalised=FALSE, calCurves = Curve14.1)
femunfurdates15.3=subset(femunfurnished, Marine.protein==15.3)
calfemunfurnished15.3=calibrate(x=femunfurdates15.3$CRA, errors=femunfurdates15.3$Error, normalised=FALSE, calCurves = Curve15.3)
femunfurdates16.5=subset(femunfurnished, Marine.protein==16.5)
calfemunfurnished16.5=calibrate(x=femunfurdates16.5$CRA, errors=femunfurdates16.5$Error, normalised=FALSE, calCurves = Curve16.5)
femunfurdates17.6=subset(femunfurnished, Marine.protein==17.6)
calfemunfurnished17.6=calibrate(x=femunfurdates17.6$CRA, errors=femunfurdates17.6$Error, normalised=FALSE, calCurves = Curve17.6)
femunfurdates18.6=subset(femunfurnished, Marine.protein==18.6)
calfemunfurnished18.6=calibrate(x=femunfurdates18.6$CRA, errors=femunfurdates18.6$Error, normalised=FALSE, calCurves = Curve18.6)
femunfurdates18.8=subset(femunfurnished, Marine.protein==18.8)
calfemunfurnished18.8=calibrate(x=femunfurdates18.8$CRA, errors=femunfurdates18.8$Error, normalised=FALSE, calCurves = Curve18.8)
femunfurdates20=subset(femunfurnished, Marine.protein==20)
calfemunfurnished20=calibrate(x=femunfurdates20$CRA, errors=femunfurdates20$Error, normalised=FALSE, calCurves = Curve20)
femunfurdates21.2=subset(femunfurnished, Marine.protein==21.2)
calfemunfurnished21.2=calibrate(x=femunfurdates21.2$CRA, errors=femunfurdates21.2$Error, normalised=FALSE, calCurves = Curve21.2)
femunfurdates22.4=subset(femunfurnished, Marine.protein==22.4)
calfemunfurnished22.4=calibrate(x=femunfurdates22.4$CRA, errors=femunfurdates22.4$Error, normalised=FALSE, calCurves = Curve22.4)
femunfurdates23.5=subset(femunfurnished, Marine.protein==23.5)
calfemunfurnished23.5=calibrate(x=femunfurdates23.5$CRA, errors=femunfurdates23.5$Error, normalised=FALSE, calCurves = Curve23.5)
femunfurdates24.7=subset(femunfurnished, Marine.protein==24.7)
calfemunfurnished24.7=calibrate(x=femunfurdates24.7$CRA, errors=femunfurdates24.7$Error, normalised=FALSE, calCurves = Curve24.7)
femunfurdates25.9=subset(femunfurnished, Marine.protein==25.9)
calfemunfurnished25.9=calibrate(x=femunfurdates25.9$CRA, errors=femunfurdates25.9$Error, normalised=FALSE, calCurves = Curve25.9)
femunfurdates27.1=subset(femunfurnished, Marine.protein==27.1)
calfemunfurnished27.1=calibrate(x=femunfurdates27.1$CRA, errors=femunfurdates27.1$Error, normalised=FALSE, calCurves = Curve27.1)
femunfurdates29.4=subset(femunfurnished, Marine.protein==29.4)
calfemunfurnished29.4=calibrate(x=femunfurdates29.4$CRA, errors=femunfurdates29.4$Error, normalised=FALSE, calCurves = Curve29.4)
femunfurdates30.6=subset(femunfurnished, Marine.protein==30.6)
calfemunfurnished30.6=calibrate(x=femunfurdates30.6$CRA, errors=femunfurdates30.6$Error, normalised=FALSE, calCurves = Curve30.6)
femunfurdates37.6=subset(femunfurnished, Marine.protein==37.6)
calfemunfurnished37.6=calibrate(x=femunfurdates37.6$CRA, errors=femunfurdates37.6$Error, normalised=FALSE, calCurves = Curve37.6)
femunfterrestrial=subset(femunfurnished, is.na(femunfurnished$Marine.protein)| Marine.protein=="0")
calfemunfurnishedterrestrial=calibrate(x=femunfterrestrial$CRA, errors=femunfterrestrial$Error, normalised=FALSE)
combdatefemunfurnished=combine(calfemunfurnished2.4, calfemunfurnished3.5, calfemunfurnished4.7, calfemunfurnished5.9, calfemunfurnished6.4, calfemunfurnished7.1, calfemunfurnished8.2, calfemunfurnished9.4, calfemunfurnished10.5, calfemunfurnished10.6, calfemunfurnished11.8, calfemunfurnished12.9, calfemunfurnished14.1, calfemunfurnished15.3, calfemunfurnished16.5, calfemunfurnished17.6, calfemunfurnished18.6, calfemunfurnished18.8, calfemunfurnished20, calfemunfurnished21.2, calfemunfurnished22.4, calfemunfurnished23.5, calfemunfurnished24.7, calfemunfurnished25.9, calfemunfurnished27.1, calfemunfurnished29.4, calfemunfurnished30.6, calfemunfurnished37.6, calfemunfurnishedterrestrial, fixIDs = TRUE)
binsfu=binPrep(sites = femunfurnished$SiteName, ages = femunfurnished$CRA, h=50)
samplefemunfurnmarine=sampleDates(combdatefemunfurnished, bins=binsfu, nsim=1000, boot=TRUE)
kdefemunfurnished30marine=ckde(samplefemunfurnmarine, timeRange=c(1650,1050), bw=30)
plot(kdefemunfurnished30marine, calendar="BCAD", ylim=c(0,0.0025)) + title("Female unfurnished graves, n=202") 

#Male furnished graves
malfurnished=subset(furnished, Sex=="M" | Sex=="M?")
malfurnished1.2=subset(malfurnished, Marine.protein==1.2)
malfurnishedcal1.2=calibrate(x=malfurnished1.2$CRA, errors=malfurnished1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
malfurnished3.5=subset(malfurnished, Marine.protein==3.5)
malfurnishedcal3.5=calibrate(x=malfurnished3.5$CRA, errors=malfurnished3.5$Error, normalised=FALSE, calCurves = Curve3.5)
malfurnished4.7=subset(malfurnished, Marine.protein==4.7)
malfurnishedcal4.7=calibrate(x=malfurnished4.7$CRA, errors=malfurnished4.7$Error, normalised=FALSE, calCurves = Curve4.7)
malfurnished5.9=subset(malfurnished, Marine.protein==5.9)
malfurnishedcal5.9=calibrate(x=malfurnished5.9$CRA, errors=malfurnished5.9$Error, normalised=FALSE, calCurves = Curve5.9)
malfurnished6.7=subset(malfurnished, Marine.protein==6.7)
malfurnishedcal6.7=calibrate(x=malfurnished6.7$CRA, errors=malfurnished6.7$Error, normalised=FALSE, calCurves = Curve6.7)
malfurnished7.1=subset(malfurnished, Marine.protein==7.1)
malfurnishedcal7.1=calibrate(x=malfurnished7.1$CRA, errors=malfurnished7.1$Error, normalised=FALSE, calCurves = Curve7.1)
malfurnished8.2=subset(malfurnished, Marine.protein==8.2)
malfurnishedcal8.2=calibrate(x=malfurnished8.2$CRA, errors=malfurnished8.2$Error, normalised=FALSE, calCurves = Curve8.2)
malfurnished8.7=subset(malfurnished, Marine.protein==8.7)
malfurnishedcal8.7=calibrate(x=malfurnished8.7$CRA, errors=malfurnished8.7$Error, normalised=FALSE, calCurves = Curve8.7)
malfurnished9.4=subset(malfurnished, Marine.protein==9.4)
malfurnishedcal9.4=calibrate(x=malfurnished9.4$CRA, errors=malfurnished9.4$Error, normalised=FALSE, calCurves = Curve9.4)
malfurnished10.6=subset(malfurnished, Marine.protein==10.6)
malfurnishedcal10.6=calibrate(x=malfurnished10.6$CRA, errors=malfurnished10.6$Error, normalised=FALSE, calCurves = Curve10.6)
malfurnished11.8=subset(malfurnished, Marine.protein==11.8)
malfurnishedcal11.8=calibrate(x=malfurnished11.8$CRA, errors=malfurnished11.8$Error, normalised=FALSE, calCurves = Curve11.8)
malfurnished12.9=subset(malfurnished, Marine.protein==12.9)
malfurnishedcal12.9=calibrate(x=malfurnished12.9$CRA, errors=malfurnished12.9$Error, normalised=FALSE, calCurves = Curve12.9)
malfurnished14.1=subset(malfurnished, Marine.protein==14.1)
malfurnishedcal14.1=calibrate(x=malfurnished14.1$CRA, errors=malfurnished14.1$Error, normalised=FALSE, calCurves = Curve14.1)
malfurnished14.6=subset(malfurnished, Marine.protein==14.6)
malfurnishedcal14.6=calibrate(x=malfurnished14.6$CRA, errors=malfurnished14.6$Error, normalised=FALSE, calCurves = Curve14.6)
malfurnished15.3=subset(malfurnished, Marine.protein==15.3)
malfurnishedcal15.3=calibrate(x=malfurnished15.3$CRA, errors=malfurnished15.3$Error, normalised=FALSE, calCurves = Curve15.3)
malfurnished16.5=subset(malfurnished, Marine.protein==16.5)
malfurnishedcal16.5=calibrate(x=malfurnished16.5$CRA, errors=malfurnished16.5$Error, normalised=FALSE, calCurves = Curve16.5)
malfurnished17.6=subset(malfurnished, Marine.protein==17.6)
malfurnishedcal17.6=calibrate(x=malfurnished17.6$CRA, errors=malfurnished17.6$Error, normalised=FALSE, calCurves = Curve17.6)
malfurnished18.8=subset(malfurnished, Marine.protein==18.8)
malfurnishedcal18.8=calibrate(x=malfurnished18.8$CRA, errors=malfurnished18.8$Error, normalised=FALSE, calCurves = Curve18.8)
malfurnished20=subset(malfurnished, Marine.protein==20)
malfurnishedcal20=calibrate(x=malfurnished20$CRA, errors=malfurnished20$Error, normalised=FALSE, calCurves = Curve20)
malfurnished21.2=subset(malfurnished, Marine.protein==21.2)
malfurnishedcal21.2=calibrate(x=malfurnished21.2$CRA, errors=malfurnished21.2$Error, normalised=FALSE, calCurves = Curve21.2)
malfurnished22.4=subset(malfurnished, Marine.protein==22.4)
malfurnishedcal22.4=calibrate(x=malfurnished22.4$CRA, errors=malfurnished22.4$Error, normalised=FALSE, calCurves = Curve22.4)
malfurnished23.5=subset(malfurnished, Marine.protein==23.5)
malfurnishedcal23.5=calibrate(x=malfurnished23.5$CRA, errors=malfurnished23.5$Error, normalised=FALSE, calCurves = Curve23.5)
malfurnished25.9=subset(malfurnished, Marine.protein==25.9)
malfurnishedcal25.9=calibrate(x=malfurnished25.9$CRA, errors=malfurnished25.9$Error, normalised=FALSE, calCurves = Curve25.9)
malfurnished27.1=subset(malfurnished, Marine.protein==27.1)
malfurnishedcal27.1=calibrate(x=malfurnished27.1$CRA, errors=malfurnished27.1$Error, normalised=FALSE, calCurves = Curve27.1)
malfurnished28.2=subset(malfurnished, Marine.protein==28.2)
malfurnishedcal28.2=calibrate(x=malfurnished28.2$CRA, errors=malfurnished28.2$Error, normalised=FALSE, calCurves = Curve28.2)
malfurterrestrial=subset(malfurnished, is.na(malfurnished$Marine.protein)| Marine.protein=="0")
malfurnishedcalterrestrial=calibrate(x=malfurterrestrial$CRA, errors=malfurterrestrial$Error, normalised=FALSE)
combdatemalfurnished=combine(malfurnishedcal1.2, malfurnishedcal3.5, malfurnishedcal4.7, malfurnishedcal5.9, malfurnishedcal6.7, malfurnishedcal7.1, malfurnishedcal8.2, malfurnishedcal8.7, malfurnishedcal9.4, malfurnishedcal10.6, malfurnishedcal11.8, malfurnishedcal12.9, malfurnishedcal14.1, malfurnishedcal14.6, malfurnishedcal15.3, malfurnishedcal16.5, malfurnishedcal17.6, malfurnishedcal18.8, malfurnishedcal20, malfurnishedcal21.2, malfurnishedcal22.4, malfurnishedcal23.5, malfurnishedcal25.9, malfurnishedcal27.1, malfurnishedcal28.2, malfurnishedcalterrestrial, fixIDs = TRUE)
binsmf=binPrep(sites = malfurnished$SiteName, ages = malfurnished$CRA, h=50)
samplemalfurnmarine=sampleDates(combdatemalfurnished, bins=binsmf, nsim=1000, boot=TRUE)
kdemalfurnished30marine=ckde(samplemalfurnmarine, timeRange=c(1650,1050), bw=30)

#Male furnished graves, weighted by number of objects
malweighted=subset(weighted, Sex=="M" | Sex=="M?")
malweighted1.2=subset(malweighted, Marine.protein==1.2)
malweightedcal1.2=calibrate(x=malweighted1.2$CRA, errors=malweighted1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
malweighted3.5=subset(malweighted, Marine.protein==3.5)
malweightedcal3.5=calibrate(x=malweighted3.5$CRA, errors=malweighted3.5$Error, normalised=FALSE, calCurves = Curve3.5)
malweighted4.7=subset(malweighted, Marine.protein==4.7)
malweightedcal4.7=calibrate(x=malweighted4.7$CRA, errors=malweighted4.7$Error, normalised=FALSE, calCurves = Curve4.7)
malweighted5.9=subset(malweighted, Marine.protein==5.9)
malweightedcal5.9=calibrate(x=malweighted5.9$CRA, errors=malweighted5.9$Error, normalised=FALSE, calCurves = Curve5.9)
malweighted6.7=subset(malweighted, Marine.protein==6.7)
malweightedcal6.7=calibrate(x=malweighted6.7$CRA, errors=malweighted6.7$Error, normalised=FALSE, calCurves = Curve6.7)
malweighted7.1=subset(malweighted, Marine.protein==7.1)
malweightedcal7.1=calibrate(x=malweighted7.1$CRA, errors=malweighted7.1$Error, normalised=FALSE, calCurves = Curve7.1)
malweighted8.2=subset(malweighted, Marine.protein==8.2)
malweightedcal8.2=calibrate(x=malweighted8.2$CRA, errors=malweighted8.2$Error, normalised=FALSE, calCurves = Curve8.2)
malweighted8.7=subset(malweighted, Marine.protein==8.7)
malweightedcal8.7=calibrate(x=malweighted8.7$CRA, errors=malweighted8.7$Error, normalised=FALSE, calCurves = Curve8.7)
malweighted9.4=subset(malweighted, Marine.protein==9.4)
malweightedcal9.4=calibrate(x=malweighted9.4$CRA, errors=malweighted9.4$Error, normalised=FALSE, calCurves = Curve9.4)
malweighted10.6=subset(malweighted, Marine.protein==10.6)
malweightedcal10.6=calibrate(x=malweighted10.6$CRA, errors=malweighted10.6$Error, normalised=FALSE, calCurves = Curve10.6)
malweighted11.8=subset(malweighted, Marine.protein==11.8)
malweightedcal11.8=calibrate(x=malweighted11.8$CRA, errors=malweighted11.8$Error, normalised=FALSE, calCurves = Curve11.8)
malweighted12.9=subset(malweighted, Marine.protein==12.9)
malweightedcal12.9=calibrate(x=malweighted12.9$CRA, errors=malweighted12.9$Error, normalised=FALSE, calCurves = Curve12.9)
malweighted14.1=subset(malweighted, Marine.protein==14.1)
malweightedcal14.1=calibrate(x=malweighted14.1$CRA, errors=malweighted14.1$Error, normalised=FALSE, calCurves = Curve14.1)
malweighted14.6=subset(malweighted, Marine.protein==14.6)
malweightedcal14.6=calibrate(x=malweighted14.6$CRA, errors=malweighted14.6$Error, normalised=FALSE, calCurves = Curve14.6)
malweighted15.3=subset(malweighted, Marine.protein==15.3)
malweightedcal15.3=calibrate(x=malweighted15.3$CRA, errors=malweighted15.3$Error, normalised=FALSE, calCurves = Curve15.3)
malweighted16.5=subset(malweighted, Marine.protein==16.5)
malweightedcal16.5=calibrate(x=malweighted16.5$CRA, errors=malweighted16.5$Error, normalised=FALSE, calCurves = Curve16.5)
malweighted17.6=subset(malweighted, Marine.protein==17.6)
malweightedcal17.6=calibrate(x=malweighted17.6$CRA, errors=malweighted17.6$Error, normalised=FALSE, calCurves = Curve17.6)
malweighted18.8=subset(malweighted, Marine.protein==18.8)
malweightedcal18.8=calibrate(x=malweighted18.8$CRA, errors=malweighted18.8$Error, normalised=FALSE, calCurves = Curve18.8)
malweighted20=subset(malweighted, Marine.protein==20)
malweightedcal20=calibrate(x=malweighted20$CRA, errors=malweighted20$Error, normalised=FALSE, calCurves = Curve20)
malweighted21.2=subset(malweighted, Marine.protein==21.2)
malweightedcal21.2=calibrate(x=malweighted21.2$CRA, errors=malweighted21.2$Error, normalised=FALSE, calCurves = Curve21.2)
malweighted22.4=subset(malweighted, Marine.protein==22.4)
malweightedcal22.4=calibrate(x=malweighted22.4$CRA, errors=malweighted22.4$Error, normalised=FALSE, calCurves = Curve22.4)
malweighted23.5=subset(malweighted, Marine.protein==23.5)
malweightedcal23.5=calibrate(x=malweighted23.5$CRA, errors=malweighted23.5$Error, normalised=FALSE, calCurves = Curve23.5)
malweighted25.9=subset(malweighted, Marine.protein==25.9)
malweightedcal25.9=calibrate(x=malweighted25.9$CRA, errors=malweighted25.9$Error, normalised=FALSE, calCurves = Curve25.9)
malweighted27.1=subset(malweighted, Marine.protein==27.1)
malweightedcal27.1=calibrate(x=malweighted27.1$CRA, errors=malweighted27.1$Error, normalised=FALSE, calCurves = Curve27.1)
malweighted28.2=subset(malweighted, Marine.protein==28.2)
malweightedcal28.2=calibrate(x=malweighted28.2$CRA, errors=malweighted28.2$Error, normalised=FALSE, calCurves = Curve28.2)
malweightedterrestrial=subset(malweighted, is.na(malweighted$Marine.protein)| Marine.protein=="0")
malweightedcalterrestrial=calibrate(x=malweightedterrestrial$CRA, errors=malweightedterrestrial$Error, normalised=FALSE)
combdatemalweighted=combine(malweightedcal1.2, malweightedcal3.5, malweightedcal4.7, malweightedcal5.9, malweightedcal6.7, malweightedcal7.1, malweightedcal8.2, malweightedcal8.7, malweightedcal9.4, malweightedcal10.6, malweightedcal11.8, malweightedcal12.9, malweightedcal14.1, malweightedcal14.6, malweightedcal15.3, malweightedcal16.5, malweightedcal17.6, malweightedcal18.8, malweightedcal20, malweightedcal21.2, malweightedcal22.4, malweightedcal23.5, malweightedcal25.9, malweightedcal27.1, malweightedcal28.2, malweightedcalterrestrial, fixIDs = TRUE)
samplemalfurnmarine=sampleDates(combdatemalweighted, nsim=1000, boot=TRUE)
kdemalweighted30marine=ckde(samplemalfurnmarine, timeRange=c(1650,1050), bw=30)

#plot figure 6a
plot(kdemalfurnished30marine, calendar="BCAD", ylim=c(0,0.005)) + title("Male furnished graves, n=156")
par(new=TRUE)
plot(kdemalweighted30marine, calendar="BCAD", ylim=c(0,0.005), fill.col="#FFFF007D")
legend("topright", legend=c("Male furnished graves","Weighted by number of grave goods"), fill=c("lightgrey", "#FFFF007D"))

#Male unfurnished graves
malunfurnished=subset(unfurnished, Sex=="M" | Sex=="M?")
malunfurnished1.2=subset(malunfurnished, Marine.protein==1.2)
malunfurnishedcal1.2=calibrate(x=malunfurnished1.2$CRA, errors=malunfurnished1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
malunfurnished2.4=subset(malunfurnished, Marine.protein==2.4)
malunfurnishedcal2.4=calibrate(x=malunfurnished2.4$CRA, errors=malunfurnished2.4$CRA, normalised=FALSE, calCurves=Curve2.4)
malunfurnished3.5=subset(malunfurnished, Marine.protein==3.5)
malunfurnishedcal3.5=calibrate(x=malunfurnished3.5$CRA, errors=malunfurnished3.5$Error, normalised=FALSE, calCurves = Curve3.5)
malunfurnished4.7=subset(malunfurnished, Marine.protein==4.7)
malunfurnishedcal4.7=calibrate(x=malunfurnished4.7$CRA, errors=malunfurnished4.7$Error, normalised=FALSE, calCurves = Curve4.7)
malunfurnished5.9=subset(malunfurnished, Marine.protein==5.9)
malunfurnishedcal5.9=calibrate(x=malunfurnished5.9$CRA, errors=malunfurnished5.9$Error, normalised=FALSE, calCurves = Curve5.9)
malunfurnished7.1=subset(malunfurnished, Marine.protein==7.1)
malunfurnishedcal7.1=calibrate(x=malunfurnished7.1$CRA, errors=malunfurnished7.1$Error, normalised=FALSE, calCurves = Curve7.1)
malunfurnished8.2=subset(malunfurnished, Marine.protein==8.2)
malunfurnishedcal8.2=calibrate(x=malunfurnished8.2$CRA, errors=malunfurnished8.2$Error, normalised=FALSE, calCurves = Curve8.2)
malunfurnished9.4=subset(malunfurnished, Marine.protein==9.4)
malunfurnishedcal9.4=calibrate(x=malunfurnished9.4$CRA, errors=malunfurnished9.4$Error, normalised=FALSE, calCurves = Curve9.4)
malunfurnished10=subset(malunfurnished, Marine.protein==10)
malunfurnishedcal10=calibrate(x=malunfurnished10$CRA, errors=malunfurnished10$Error, normalised=FALSE, calCurves = Curve10)
malunfurnished10.6=subset(malunfurnished, Marine.protein==10.6)
malunfurnishedcal10.6=calibrate(x=malunfurnished10.6$CRA, errors=malunfurnished10.6$Error, normalised=FALSE, calCurves = Curve10.6)
malunfurnished11.5=subset(malunfurnished, Marine.protein==11.5)
malunfurnishedcal11.5=calibrate(x=malunfurnished11.5$CRA, errors=malunfurnished11.5$Error, normalised=FALSE, calCurves = Curve11.5)
malunfurnished11.8=subset(malunfurnished, Marine.protein==11.8)
malunfurnishedcal11.8=calibrate(x=malunfurnished11.8$CRA, errors=malunfurnished11.8$Error, normalised=FALSE, calCurves = Curve11.8)
malunfurnished12.2=subset(malunfurnished, Marine.protein==12.2)
malunfurnishedcal12.2=calibrate(x=malunfurnished12.2$CRA, errors=malunfurnished12.2$Error, normalised=FALSE, calCurves = Curve12.2)
malunfurnished12.9=subset(malunfurnished, Marine.protein==12.9)
malunfurnishedcal12.9=calibrate(x=malunfurnished12.9$CRA, errors=malunfurnished12.9$Error, normalised=FALSE, calCurves = Curve12.9)
malunfurnished14.1=subset(malunfurnished, Marine.protein==14.1)
malunfurnishedcal14.1=calibrate(x=malunfurnished14.1$CRA, errors=malunfurnished14.1$Error, normalised=FALSE, calCurves = Curve14.1)
malunfurnished15.3=subset(malunfurnished, Marine.protein==15.3)
malunfurnishedcal15.3=calibrate(x=malunfurnished15.3$CRA, errors=malunfurnished15.3$Error, normalised=FALSE, calCurves = Curve15.3)
malunfurnished16.5=subset(malunfurnished, Marine.protein==16.5)
malunfurnishedcal16.5=calibrate(x=malunfurnished16.5$CRA, errors=malunfurnished16.5$Error, normalised=FALSE, calCurves = Curve16.5)
malunfurnished17.6=subset(malunfurnished, Marine.protein==17.6)
malunfurnishedcal17.6=calibrate(x=malunfurnished17.6$CRA, errors=malunfurnished17.6$Error, normalised=FALSE, calCurves = Curve17.6)
malunfurnished17.8=subset(malunfurnished, Marine.protein==17.8)
malunfurnishedcal17.8=calibrate(x=malunfurnished17.8$CRA, errors=malunfurnished17.8$Error, normalised=FALSE, calCurves = Curve17.8)
malunfurnished18.8=subset(malunfurnished, Marine.protein==18.8)
malunfurnishedcal18.8=calibrate(x=malunfurnished18.8$CRA, errors=malunfurnished18.8$Error, normalised=FALSE, calCurves = Curve18.8)
malunfurnished20=subset(malunfurnished, Marine.protein==20)
malunfurnishedcal20=calibrate(x=malunfurnished20$CRA, errors=malunfurnished20$Error, normalised=FALSE, calCurves = Curve20)
malunfurnished20.4=subset(malunfurnished, Marine.protein==20.4)
malunfurnishedcal20.4=calibrate(x=malunfurnished20.4$CRA, errors=malunfurnished20.4$Error, normalised=FALSE, calCurves = Curve20.4)
malunfurnished21.1=subset(malunfurnished, Marine.protein==21.1)
malunfurnishedcal21.1=calibrate(x=malunfurnished21.1$CRA, errors=malunfurnished21.1$Error, normalised=FALSE, calCurves = Curve21.1)
malunfurnished21.2=subset(malunfurnished, Marine.protein==21.2)
malunfurnishedcal21.2=calibrate(x=malunfurnished21.2$CRA, errors=malunfurnished21.2$Error, normalised=FALSE, calCurves = Curve21.2)
malunfurnished22.4=subset(malunfurnished, Marine.protein==22.4)
malunfurnishedcal22.4=calibrate(x=malunfurnished22.4$CRA, errors=malunfurnished22.4$Error, normalised=FALSE, calCurves = Curve22.4)
malunfurnished23.5=subset(malunfurnished, Marine.protein==23.5)
malunfurnishedcal23.5=calibrate(x=malunfurnished23.5$CRA, errors=malunfurnished23.5$Error, normalised=FALSE, calCurves = Curve23.5)
malunfurnished24.7=subset(malunfurnished, Marine.protein==24.7)
malunfurnishedcal24.7=calibrate(x=malunfurnished24.7$CRA, errors=malunfurnished24.7$Error, normalised=FALSE, calCurves = Curve24.7)
malunfurnished25.9=subset(malunfurnished, Marine.protein==25.9)
malunfurnishedcal25.9=calibrate(x=malunfurnished25.9$CRA, errors=malunfurnished25.9$Error, normalised=FALSE, calCurves = Curve25.9)
malunfurnished27.1=subset(malunfurnished, Marine.protein==27.1)
malunfurnishedcal27.1=calibrate(x=malunfurnished27.1$CRA, errors=malunfurnished27.1$Error, normalised=FALSE, calCurves = Curve27.1)
malunfurnished28.2=subset(malunfurnished, Marine.protein==28.2)
malunfurnishedcal28.2=calibrate(x=malunfurnished28.2$CRA, errors=malunfurnished28.2$Error, normalised=FALSE, calCurves = Curve28.2)
malunfurnished29.4=subset(malunfurnished, Marine.protein==29.4)
malunfurnishedcal29.4=calibrate(x=malunfurnished29.4$CRA, errors=malunfurnished29.4$Error, normalised=FALSE, calCurves = Curve29.4)
malunfurnished31.8=subset(malunfurnished, Marine.protein==31.8)
malunfurnishedcal31.8=calibrate(x=malunfurnished31.8$CRA, errors=malunfurnished31.8$Error, normalised=FALSE, calCurves = Curve31.8)
malunfurnished32=subset(malunfurnished, Marine.protein==32)
malunfurnishedcal32=calibrate(x=malunfurnished32$CRA, errors=malunfurnished32$Error, normalised=FALSE, calCurves = Curve32)
malunfurnished32.9=subset(malunfurnished, Marine.protein==32.9)
malunfurnishedcal32.9=calibrate(x=malunfurnished32.9$CRA, errors=malunfurnished32.9$Error, normalised=FALSE, calCurves = Curve32.9)
malunfurnished38.8=subset(malunfurnished, Marine.protein==38.8)
malunfurnishedcal38.8=calibrate(x=malunfurnished38.8$CRA, errors=malunfurnished38.8$Error, normalised=FALSE, calCurves = Curve38.8)
terrestrialmalunfur=subset(malunfurnished, is.na(malunfurnished$Marine.protein)| Marine.protein=="0")
malunfurnishedcalterrestrial=calibrate(x=terrestrialmalunfur$CRA, errors=terrestrialmalunfur$Error, normalised=FALSE)
combdatemalunfurnished=combine(malunfurnishedcal1.2, malunfurnishedcal2.4, malunfurnishedcal3.5, malunfurnishedcal4.7, malunfurnishedcal5.9, malunfurnishedcal7.1, malunfurnishedcal8.2, malunfurnishedcal9.4, malunfurnishedcal10, malunfurnishedcal10.6, malunfurnishedcal11.5, malunfurnishedcal11.8, malunfurnishedcal12.2, malunfurnishedcal12.9, malunfurnishedcal14.1, malunfurnishedcal15.3, malunfurnishedcal16.5, malunfurnishedcal17.6, malunfurnishedcal17.8, malunfurnishedcal18.8, malunfurnishedcal20, malunfurnishedcal20.4, malunfurnishedcal21.1, malunfurnishedcal21.2, malunfurnishedcal22.4, malunfurnishedcal23.5, malunfurnishedcal24.7, malunfurnishedcal25.9, malunfurnishedcal27.1, malunfurnishedcal28.2, malunfurnishedcal29.4, malunfurnishedcal31.8, malunfurnishedcal32, malunfurnishedcal32.9, malunfurnishedcal38.8, malunfurnishedcalterrestrial, fixIDs = TRUE)
binsmu=binPrep(sites = malunfurnished$SiteName, ages = malunfurnished$CRA, h=50)
samplemalunfurnmarine=sampleDates(combdatemalunfurnished, bins=binsmu, nsim=1000, boot=TRUE)
kdemalunfurnished30marine=ckde(samplemalunfurnmarine, timeRange=c(1650,1050), bw=30)
plot(kdemalunfurnished30marine, calendar="BCAD", ylim=c(0,0.002)) + title("Male unfurnished graves, n=269")

#Furnished burials in the north and west
furnishednw=subset(furnished, Region=="NW")
furnishednw1.2=subset(furnishednw, Marine.protein==1.2)
furnishednwcal1.2=calibrate(x=furnishednw1.2$CRA, errors=furnishednw1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
furnishednw2.4=subset(furnishednw, Marine.protein==2.4)
furnishednwcal2.4=calibrate(x=furnishednw2.4$CRA, errors=furnishednw2.4$Error, normalised=FALSE, calCurves = Curve2.4)
furnishednw4.7=subset(furnishednw, Marine.protein==4.7)
furnishednwcal4.7=calibrate(x=furnishednw4.7$CRA, errors=furnishednw4.7$Error, normalised=FALSE, calCurves = Curve4.7)
furnishednw5.9=subset(furnishednw, Marine.protein==5.9)
furnishednwcal5.9=calibrate(x=furnishednw5.9$CRA, errors=furnishednw5.9$Error, normalised=FALSE, calCurves = Curve5.9)
furnishednw8.2=subset(furnishednw, Marine.protein==8.2)
furnishednwcal8.2=calibrate(x=furnishednw8.2$CRA, errors=furnishednw8.2$Error, normalised=FALSE, calCurves = Curve8.2)
furnishednw9.4=subset(furnishednw, Marine.protein==9.4)
furnishednwcal9.4=calibrate(x=furnishednw9.4$CRA, errors=furnishednw9.4$Error, normalised=FALSE, calCurves = Curve9.4)
furnishednw9.5=subset(furnishednw, Marine.protein==9.5)
furnishednwcal9.5=calibrate(x=furnishednw9.5$CRA, errors=furnishednw9.5$Error, normalised=FALSE, calCurves = Curve9.5)
furnishednw10.6=subset(furnishednw, Marine.protein==10.6)
furnishednwcal10.6=calibrate(x=furnishednw10.6$CRA, errors=furnishednw10.6$Error, normalised=FALSE, calCurves = Curve10.6)
furnishednw11.8=subset(furnishednw, Marine.protein==11.8)
furnishednwcal11.8=calibrate(x=furnishednw11.8$CRA, errors=furnishednw11.8$Error, normalised=FALSE, calCurves = Curve11.8)
furnishednw12.9=subset(furnishednw, Marine.protein==12.9)
furnishednwcal12.9=calibrate(x=furnishednw12.9$CRA, errors=furnishednw12.9$Error, normalised=FALSE, calCurves = Curve12.9)
furnishednw14.1=subset(furnishednw, Marine.protein==14.1)
furnishednwcal14.1=calibrate(x=furnishednw14.1$CRA, errors=furnishednw14.1$Error, normalised=FALSE, calCurves = Curve14.1)
furnishednw15.3=subset(furnishednw, Marine.protein==15.3)
furnishednwcal15.3=calibrate(x=furnishednw15.3$CRA, errors=furnishednw15.3$Error, normalised=FALSE, calCurves = Curve15.3)
furnishednw16.5=subset(furnishednw, Marine.protein==16.5)
furnishednwcal16.5=calibrate(x=furnishednw16.5$CRA, errors=furnishednw16.5$Error, normalised=FALSE, calCurves = Curve16.5)
furnishednw17.6=subset(furnishednw, Marine.protein==17.6)
furnishednwcal17.6=calibrate(x=furnishednw17.6$CRA, errors=furnishednw17.6$Error, normalised=FALSE, calCurves = Curve17.6)
furnishednw22.4=subset(furnishednw, Marine.protein==22.4)
furnishednwcal22.4=calibrate(x=furnishednw22.4$CRA, errors=furnishednw22.4$Error, normalised=FALSE, calCurves = Curve22.4)
furnishednw23.5=subset(furnishednw, Marine.protein==23.5)
furnishednwcal23.5=calibrate(x=furnishednw23.5$CRA, errors=furnishednw23.5$Error, normalised=FALSE, calCurves = Curve23.5)
furnishednw27.1=subset(furnishednw, Marine.protein==27.1)
furnishednwcal27.1=calibrate(x=furnishednw27.1$CRA, errors=furnishednw27.1$Error, normalised=FALSE, calCurves = Curve27.1)
furnishednw28.2=subset(furnishednw, Marine.protein==28.2)
furnishednwcal28.2=calibrate(x=furnishednw28.2$CRA, errors=furnishednw28.2$Error, normalised=FALSE, calCurves = Curve28.2)
furnishednw32.9=subset(furnishednw, Marine.protein==32.9)
furnishednwcal32.9=calibrate(x=furnishednw32.9$CRA, errors=furnishednw32.9$Error, normalised=FALSE, calCurves = Curve32.9)
terrestrialfurnw=subset(furnishednw, is.na(furnishednw$Marine.protein)| Marine.protein=="0")
furnishednwcalterrestrial=calibrate(x=terrestrialfurnw$CRA, errors=terrestrialfurnw$Error, normalised=FALSE)
combfurnishednwdate=combine(furnishednwcal1.2, furnishednwcal2.4, furnishednwcal4.7, furnishednwcal5.9, furnishednwcal8.2, furnishednwcal9.4, furnishednwcal9.5, furnishednwcal10.6, furnishednwcal11.8, furnishednwcal12.9, furnishednwcal14.1, furnishednwcal15.3, furnishednwcal16.5, furnishednwcal17.6, furnishednwcal22.4, furnishednwcal23.5, furnishednwcal27.1, furnishednwcal28.2, furnishednwcal32.9, furnishednwcalterrestrial, fixIDs = TRUE)
binsfnw=binPrep(sites = furnishednw$SiteName, ages=furnishednw$CRA, h=50)
samplefurnnwmarine=sampleDates(combfurnishednwdate, bins=binsfnw, nsim=1000, boot=TRUE)
kdemarinefurnishednw30=ckde(samplefurnnwmarine, timeRange=c(1650,1050), bw=30)
plot(kdemarinefurnishednw30, calendar="BCAD", ylim=c(0,0.004)) +title("Furnished graves in the north and west, n=72")

#Furnished burials in the north and west, weighted by number of grave goods
weightednw=subset(weighted, Region=="NW")
weightednw1.2=subset(weightednw, Marine.protein==1.2)
weightednwcal1.2=calibrate(x=weightednw1.2$CRA, errors=weightednw1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
weightednw2.4=subset(weightednw, Marine.protein==2.4)
weightednwcal2.4=calibrate(x=weightednw2.4$CRA, errors=weightednw2.4$Error, normalised=FALSE, calCurves = Curve2.4)
weightednw4.7=subset(weightednw, Marine.protein==4.7)
weightednwcal4.7=calibrate(x=weightednw4.7$CRA, errors=weightednw4.7$Error, normalised=FALSE, calCurves = Curve4.7)
weightednw5.9=subset(weightednw, Marine.protein==5.9)
weightednwcal5.9=calibrate(x=weightednw5.9$CRA, errors=weightednw5.9$Error, normalised=FALSE, calCurves = Curve5.9)
weightednw8.2=subset(weightednw, Marine.protein==8.2)
weightednwcal8.2=calibrate(x=weightednw8.2$CRA, errors=weightednw8.2$Error, normalised=FALSE, calCurves = Curve8.2)
weightednw9.4=subset(weightednw, Marine.protein==9.4)
weightednwcal9.4=calibrate(x=weightednw9.4$CRA, errors=weightednw9.4$Error, normalised=FALSE, calCurves = Curve9.4)
weightednw9.5=subset(weightednw, Marine.protein==9.5)
weightednwcal9.5=calibrate(x=weightednw9.5$CRA, errors=weightednw9.5$Error, normalised=FALSE, calCurves = Curve9.5)
weightednw10.6=subset(weightednw, Marine.protein==10.6)
weightednwcal10.6=calibrate(x=weightednw10.6$CRA, errors=weightednw10.6$Error, normalised=FALSE, calCurves = Curve10.6)
weightednw11.8=subset(weightednw, Marine.protein==11.8)
weightednwcal11.8=calibrate(x=weightednw11.8$CRA, errors=weightednw11.8$Error, normalised=FALSE, calCurves = Curve11.8)
weightednw12.9=subset(weightednw, Marine.protein==12.9)
weightednwcal12.9=calibrate(x=weightednw12.9$CRA, errors=weightednw12.9$Error, normalised=FALSE, calCurves = Curve12.9)
weightednw14.1=subset(weightednw, Marine.protein==14.1)
weightednwcal14.1=calibrate(x=weightednw14.1$CRA, errors=weightednw14.1$Error, normalised=FALSE, calCurves = Curve14.1)
weightednw15.3=subset(weightednw, Marine.protein==15.3)
weightednwcal15.3=calibrate(x=weightednw15.3$CRA, errors=weightednw15.3$Error, normalised=FALSE, calCurves = Curve15.3)
weightednw16.5=subset(weightednw, Marine.protein==16.5)
weightednwcal16.5=calibrate(x=weightednw16.5$CRA, errors=weightednw16.5$Error, normalised=FALSE, calCurves = Curve16.5)
weightednw17.6=subset(weightednw, Marine.protein==17.6)
weightednwcal17.6=calibrate(x=weightednw17.6$CRA, errors=weightednw17.6$Error, normalised=FALSE, calCurves = Curve17.6)
weightednw22.4=subset(weightednw, Marine.protein==22.4)
weightednwcal22.4=calibrate(x=weightednw22.4$CRA, errors=weightednw22.4$Error, normalised=FALSE, calCurves = Curve22.4)
weightednw23.5=subset(weightednw, Marine.protein==23.5)
weightednwcal23.5=calibrate(x=weightednw23.5$CRA, errors=weightednw23.5$Error, normalised=FALSE, calCurves = Curve23.5)
weightednw27.1=subset(weightednw, Marine.protein==27.1)
weightednwcal27.1=calibrate(x=weightednw27.1$CRA, errors=weightednw27.1$Error, normalised=FALSE, calCurves = Curve27.1)
weightednw28.2=subset(weightednw, Marine.protein==28.2)
weightednwcal28.2=calibrate(x=weightednw28.2$CRA, errors=weightednw28.2$Error, normalised=FALSE, calCurves = Curve28.2)
weightednw32.9=subset(weightednw, Marine.protein==32.9)
weightednwcal32.9=calibrate(x=weightednw32.9$CRA, errors=weightednw32.9$Error, normalised=FALSE, calCurves = Curve32.9)
terrestrialnwweighted=subset(weightednw, is.na(weightednw$Marine.protein)| Marine.protein=="0")
weightednwcalterrestrial=calibrate(x=terrestrialnwweighted$CRA, errors=terrestrialnwweighted$Error, normalised=FALSE)
combweightednwdate=combine(weightednwcal1.2, weightednwcal2.4, weightednwcal4.7, weightednwcal5.9, weightednwcal8.2, weightednwcal9.4, weightednwcal9.5, weightednwcal10.6, weightednwcal11.8, weightednwcal12.9, weightednwcal14.1, weightednwcal15.3, weightednwcal16.5, weightednwcal17.6, weightednwcal22.4, weightednwcal23.5, weightednwcal27.1, weightednwcal28.2, weightednwcal32.9, weightednwcalterrestrial, fixIDs = TRUE)
sampleweightednwmarine=sampleDates(combweightednwdate, nsim=1000, boot=TRUE)
kdemarineweightednw30=ckde(sampleweightednwmarine, timeRange=c(1650,1050), bw=30)
plot(kdemarineweightednw30, calendar="BCAD",ylim=c(0,0.004)) +title("Furnished graves in the north and west, weighted by number of objects per grave")

#plot figure 8a
plot(kdemarinefurnishednw30, calendar="BCAD", ylim=c(0,0.004)) +title("Furnished graves in the north and west, n=72")
par(new=TRUE)
plot(kdemarineweightednw30, calendar="BCAD",ylim=c(0,0.004), fill.col="#FFFF007D")
legend("topright", legend=c("Furnished graves in the north and west","Weighted by number of grave goods"), fill=c("lightgrey", "#FFFF007D"))

#Unfurnished burials in the north and west
unfurnishednw=subset(unfurnished, Region=="NW")
unfurnishednw1.2=subset(unfurnishednw, Marine.protein==1.2)
unfurnishednwcal1.2=calibrate(x=unfurnishednw1.2$CRA, errors=unfurnishednw1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
unfurnishednw2.4=subset(unfurnishednw, Marine.protein==2.4)
unfurnishednwcal2.4=calibrate(x=unfurnishednw2.4$CRA, errors=unfurnishednw2.4$Error, normalised=FALSE, calCurves = Curve2.4)
unfurnishednw3.5=subset(unfurnishednw, Marine.protein==3.5)
unfurnishednwcal3.5=calibrate(x=unfurnishednw3.5$CRA, errors=unfurnishednw3.5$Error, normalised=FALSE, calCurves = Curve3.5)
unfurnishednw4.7=subset(unfurnishednw, Marine.protein==4.7)
unfurnishednwcal4.7=calibrate(x=unfurnishednw4.7$CRA, errors=unfurnishednw4.7$Error, normalised=FALSE, calCurves = Curve4.7)
unfurnishednw5.9=subset(unfurnishednw, Marine.protein==5.9)
unfurnishednwcal5.9=calibrate(x=unfurnishednw5.9$CRA, errors=unfurnishednw5.9$Error, normalised=FALSE, calCurves = Curve5.9)
unfurnishednw7.1=subset(unfurnishednw, Marine.protein==7.1)
unfurnishednwcal7.1=calibrate(x=unfurnishednw7.1$CRA, errors=unfurnishednw7.1$Error, normalised=FALSE, calCurves = Curve7.1)
unfurnishednw8.2=subset(unfurnishednw, Marine.protein==8.2)
unfurnishednwcal8.2=calibrate(x=unfurnishednw8.2$CRA, errors=unfurnishednw8.2$Error, normalised=FALSE, calCurves = Curve8.2)
unfurnishednw9.4=subset(unfurnishednw, Marine.protein==9.4)
unfurnishednwcal9.4=calibrate(x=unfurnishednw9.4$CRA, errors=unfurnishednw9.4$Error, normalised=FALSE, calCurves = Curve9.4)
unfurnishednw10.6=subset(unfurnishednw, Marine.protein==10.6)
unfurnishednwcal10.6=calibrate(x=unfurnishednw10.6$CRA, errors=unfurnishednw10.6$Error, normalised=FALSE, calCurves = Curve10.6)
unfurnishednw11.8=subset(unfurnishednw, Marine.protein==11.8)
unfurnishednwcal11.8=calibrate(x=unfurnishednw11.8$CRA, errors=unfurnishednw11.8$Error, normalised=FALSE, calCurves = Curve11.8)
unfurnishednw12.9=subset(unfurnishednw, Marine.protein==12.9)
unfurnishednwcal12.9=calibrate(x=unfurnishednw12.9$CRA, errors=unfurnishednw12.9$Error, normalised=FALSE, calCurves = Curve12.9)
unfurnishednw14.1=subset(unfurnishednw, Marine.protein==14.1)
unfurnishednwcal14.1=calibrate(x=unfurnishednw14.1$CRA, errors=unfurnishednw14.1$Error, normalised=FALSE, calCurves = Curve14.1)
unfurnishednw15.3=subset(unfurnishednw, Marine.protein==15.3)
unfurnishednwcal15.3=calibrate(x=unfurnishednw15.3$CRA, errors=unfurnishednw15.3$Error, normalised=FALSE, calCurves = Curve15.3)
unfurnishednw16.5=subset(unfurnishednw, Marine.protein==16.5)
unfurnishednwcal16.5=calibrate(x=unfurnishednw16.5$CRA, errors=unfurnishednw16.5$Error, normalised=FALSE, calCurves = Curve16.5)
unfurnishednw17.6=subset(unfurnishednw, Marine.protein==17.6)
unfurnishednwcal17.6=calibrate(x=unfurnishednw17.6$CRA, errors=unfurnishednw17.6$Error, normalised=FALSE, calCurves = Curve17.6)
unfurnishednw18.8=subset(unfurnishednw, Marine.protein==18.8)
unfurnishednwcal18.8=calibrate(x=unfurnishednw18.8$CRA, errors=unfurnishednw18.8$Error, normalised=FALSE, calCurves = Curve18.8)
unfurnishednw21.1=subset(unfurnishednw, Marine.protein==21.1)
unfurnishednwcal21.1=calibrate(x=unfurnishednw21.1$CRA, errors=unfurnishednw21.1$Error, normalised=FALSE, calCurves = Curve21.1)
unfurnishednw21.2=subset(unfurnishednw, Marine.protein==21.2)
unfurnishednwcal21.2=calibrate(x=unfurnishednw21.2$CRA, errors=unfurnishednw21.2$Error, normalised=FALSE, calCurves = Curve21.2)
unfurnishednw22.4=subset(unfurnishednw, Marine.protein==22.4)
unfurnishednwcal22.4=calibrate(x=unfurnishednw22.4$CRA, errors=unfurnishednw22.4$Error, normalised=FALSE, calCurves = Curve22.4)
unfurnishednw23.5=subset(unfurnishednw, Marine.protein==23.5)
unfurnishednwcal23.5=calibrate(x=unfurnishednw23.5$CRA, errors=unfurnishednw23.5$Error, normalised=FALSE, calCurves = Curve23.5)
unfurnishednw24.7=subset(unfurnishednw, Marine.protein==24.7)
unfurnishednwcal24.7=calibrate(x=unfurnishednw24.7$CRA, errors=unfurnishednw24.7$Error, normalised=FALSE, calCurves = Curve24.7)
unfurnishednw25.9=subset(unfurnishednw, Marine.protein==25.9)
unfurnishednwcal25.9=calibrate(x=unfurnishednw25.9$CRA, errors=unfurnishednw25.9$Error, normalised=FALSE, calCurves = Curve25.9)
unfurnishednw28.2=subset(unfurnishednw, Marine.protein==28.2)
unfurnishednwcal28.2=calibrate(x=unfurnishednw28.2$CRA, errors=unfurnishednw28.2$Error, normalised=FALSE, calCurves = Curve28.2)
unfurnishednw29.4=subset(unfurnishednw, Marine.protein==29.4)
unfurnishednwcal29.4=calibrate(x=unfurnishednw29.4$CRA, errors=unfurnishednw29.4$Error, normalised=FALSE, calCurves = Curve29.4)
unfurnishednw31.8=subset(unfurnishednw, Marine.protein==31.8)
unfurnishednwcal31.8=calibrate(x=unfurnishednw31.8$CRA, errors=unfurnishednw31.8$Error, normalised=FALSE, calCurves = Curve31.8)
unfurnishednw34.1=subset(unfurnishednw, Marine.protein==34.1)
unfurnishednwcal34.1=calibrate(x=unfurnishednw34.1$CRA, errors=unfurnishednw34.1$Error, normalised=FALSE, calCurves = Curve34.1)
unfurnishednw38.8=subset(unfurnishednw, Marine.protein==38.8)
unfurnishednwcal38.8=calibrate(x=unfurnishednw38.8$CRA, errors=unfurnishednw38.8$Error, normalised=FALSE, calCurves = Curve38.8)
terrestrialunfnw=subset(unfurnishednw, is.na(unfurnishednw$Marine.protein)| Marine.protein=="0")
unfurnishednwcalterrestrial=calibrate(x=terrestrialunfnw$CRA, errors=terrestrialunfnw$Error, normalised=FALSE)
combdateunfnw=combine(unfurnishednwcal1.2, unfurnishednwcal2.4, unfurnishednwcal3.5, unfurnishednwcal4.7, unfurnishednwcal5.9, unfurnishednwcal7.1, unfurnishednwcal8.2, unfurnishednwcal9.4, unfurnishednwcal10.6, unfurnishednwcal11.8, unfurnishednwcal12.9, unfurnishednwcal14.1, unfurnishednwcal15.3, unfurnishednwcal16.5, unfurnishednwcal17.6, unfurnishednwcal18.8, unfurnishednwcal21.1, unfurnishednwcal21.2, unfurnishednwcal22.4, unfurnishednwcal23.5, unfurnishednwcal24.7, unfurnishednwcal25.9, unfurnishednwcal28.2, unfurnishednwcal29.4, unfurnishednwcal31.8, unfurnishednwcal34.1, unfurnishednwcal38.8, unfurnishednwcalterrestrial, fixIDs = TRUE)
binsunw=binPrep(sites = unfurnishednw$SiteName, ages=unfurnishednw$CRA, h=50)
sampleunfurnishednwmarine=sampleDates(combdateunfnw, bins=binsunw, nsim=1000, boot=TRUE)
kdeunfurnishednw30marine=ckde(sampleunfurnishednwmarine, timeRange=c(1650,1050), bw=30)
plot(kdeunfurnishednw30marine, calendar="BCAD",ylim=c(0,0.002)) +title("Unfurnished graves in the north and west, n=348")

#furnished burial in the south and east
furnishedse=subset(furnished, Region=="SE")
furnishedse1.2=subset(furnishedse, Marine.protein==1.2)
furnishedsecal1.2=calibrate(x=furnishedse1.2$CRA, errors=furnishedse1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
furnishedse2.4=subset(furnishedse, Marine.protein==2.4)
furnishedsecal2.4=calibrate(x=furnishedse2.4$CRA, errors=furnishedse2.4$Error, normalised=FALSE, calCurves = Curve2.4)
furnishedse3.5=subset(furnishedse, Marine.protein==3.5)
furnishedsecal3.5=calibrate(x=furnishedse3.5$CRA, errors=furnishedse3.5$Error, normalised=FALSE, calCurves = Curve3.5)
furnishedse4.7=subset(furnishedse, Marine.protein==4.7)
furnishedsecal4.7=calibrate(x=furnishedse4.7$CRA, errors=furnishedse4.7$Error, normalised=FALSE, calCurves = Curve4.7)
furnishedse5.9=subset(furnishedse, Marine.protein==5.9)
furnishedsecal5.9=calibrate(x=furnishedse5.9$CRA, errors=furnishedse5.9$Error, normalised=FALSE, calCurves = Curve5.9)
furnishedse6.7=subset(furnishedse, Marine.protein==6.7)
furnishedsecal6.7=calibrate(x=furnishedse6.7$CRA, errors=furnishedse6.7$Error, normalised=FALSE, calCurves = Curve6.7)
furnishedse7.1=subset(furnishedse, Marine.protein==7.1)
furnishedsecal7.1=calibrate(x=furnishedse7.1$CRA, errors=furnishedse7.1$Error, normalised=FALSE, calCurves = Curve7.1)
furnishedse8=subset(furnishedse, Marine.protein==8)
furnishedsecal8=calibrate(x=furnishedse8$CRA, errors=furnishedse8$Error, normalised=FALSE, calCurves = Curve8)
furnishedse8.2=subset(furnishedse, Marine.protein==8.2)
furnishedsecal8.2=calibrate(x=furnishedse8.2$CRA, errors=furnishedse8.2$Error, normalised=FALSE, calCurves = Curve8.2)
furnishedse8.3=subset(furnishedse, Marine.protein==8.3)
furnishedsecal8.3=calibrate(x=furnishedse8.3$CRA, errors=furnishedse8.3$Error, normalised=FALSE, calCurves = Curve8.3)
furnishedse8.7=subset(furnishedse, Marine.protein==8.7)
furnishedsecal8.7=calibrate(x=furnishedse8.7$CRA, errors=furnishedse8.7$Error, normalised=FALSE, calCurves = Curve8.7)
furnishedse9.4=subset(furnishedse, Marine.protein==9.4)
furnishedsecal9.4=calibrate(x=furnishedse9.4$CRA, errors=furnishedse9.4$Error, normalised=FALSE, calCurves = Curve9.4)
furnishedse9.6=subset(furnishedse, Marine.protein==9.6)
furnishedsecal9.6=calibrate(x=furnishedse9.6$CRA, errors=furnishedse9.6$Error, normalised=FALSE, calCurves = Curve9.6)
furnishedse10.1=subset(furnishedse, Marine.protein==10.1)
furnishedsecal10.1=calibrate(x=furnishedse10.1$CRA, errors=furnishedse10.1$Error, normalised=FALSE, calCurves = Curve10.1)
furnishedse10.6=subset(furnishedse, Marine.protein==10.6)
furnishedsecal10.6=calibrate(x=furnishedse10.6$CRA, errors=furnishedse10.6$Error, normalised=FALSE, calCurves = Curve10.6)
furnishedse11.8=subset(furnishedse, Marine.protein==11.8)
furnishedsecal11.8=calibrate(x=furnishedse11.8$CRA, errors=furnishedse11.8$Error, normalised=FALSE, calCurves = Curve11.8)
furnishedse12.9=subset(furnishedse, Marine.protein==12.9)
furnishedsecal12.9=calibrate(x=furnishedse12.9$CRA, errors=furnishedse12.9$Error, normalised=FALSE, calCurves = Curve12.9)
furnishedse13.9=subset(furnishedse, Marine.protein==13.9)
furnishedsecal13.9=calibrate(x=furnishedse13.9$CRA, errors=furnishedse13.9$Error, normalised=FALSE, calCurves = Curve13.9)
furnishedse14.1=subset(furnishedse, Marine.protein==14.1)
furnishedsecal14.1=calibrate(x=furnishedse14.1$CRA, errors=furnishedse14.1$Error, normalised=FALSE, calCurves = Curve14.1)
furnishedse14.6=subset(furnishedse, Marine.protein==14.6)
furnishedsecal14.6=calibrate(x=furnishedse14.6$CRA, errors=furnishedse14.6$Error, normalised=FALSE, calCurves = Curve14.6)
furnishedse15.3=subset(furnishedse, Marine.protein==15.3)
furnishedsecal15.3=calibrate(x=furnishedse15.3$CRA, errors=furnishedse15.3$Error, normalised=FALSE, calCurves = Curve15.3)
furnishedse16.5=subset(furnishedse, Marine.protein==16.5)
furnishedsecal16.5=calibrate(x=furnishedse16.5$CRA, errors=furnishedse16.5$Error, normalised=FALSE, calCurves = Curve16.5)
furnishedse17.6=subset(furnishedse, Marine.protein==17.6)
furnishedsecal17.6=calibrate(x=furnishedse17.6$CRA, errors=furnishedse17.6$Error, normalised=FALSE, calCurves = Curve17.6)
furnishedse18.8=subset(furnishedse, Marine.protein==18.8)
furnishedsecal18.8=calibrate(x=furnishedse18.8$CRA, errors=furnishedse18.8$Error, normalised=FALSE, calCurves = Curve18.8)
furnishedse20=subset(furnishedse, Marine.protein==20)
furnishedsecal20=calibrate(x=furnishedse20$CRA, errors=furnishedse20$Error, normalised=FALSE, calCurves = Curve20)
furnishedse21.2=subset(furnishedse, Marine.protein==21.2)
furnishedsecal21.2=calibrate(x=furnishedse21.2$CRA, errors=furnishedse21.2$Error, normalised=FALSE, calCurves = Curve21.2)
furnishedse22.4=subset(furnishedse, Marine.protein==22.4)
furnishedsecal22.4=calibrate(x=furnishedse22.4$CRA, errors=furnishedse22.4$Error, normalised=FALSE, calCurves = Curve22.4)
furnishedse23.5=subset(furnishedse, Marine.protein==23.5)
furnishedsecal23.5=calibrate(x=furnishedse23.5$CRA, errors=furnishedse23.5$Error, normalised=FALSE, calCurves = Curve23.5)
furnishedse25.9=subset(furnishedse, Marine.protein==25.9)
furnishedsecal25.9=calibrate(x=furnishedse25.9$CRA, errors=furnishedse25.9$Error, normalised=FALSE, calCurves = Curve25.9)
furnishedse27.1=subset(furnishedse, Marine.protein==27.1)
furnishedsecal27.1=calibrate(x=furnishedse27.1$CRA, errors=furnishedse27.1$Error, normalised=FALSE, calCurves = Curve27.1)
furnishedse28.4=subset(furnishedse, Marine.protein==28.4)
furnishedsecal28.4=calibrate(x=furnishedse28.4$CRA, errors=furnishedse28.4$Error, normalised=FALSE, calCurves = Curve28.4)
furnishedse32.7=subset(furnishedse, Marine.protein==32.7)
furnishedsecal32.7=calibrate(x=furnishedse32.7$CRA, errors=furnishedse32.7$Error, normalised=FALSE, calCurves = Curve32.7)
furnishedse37.6=subset(furnishedse, Marine.protein==37.6)
furnishedsecal37.6=calibrate(x=furnishedse37.6$CRA, errors=furnishedse37.6$Error, normalised=FALSE, calCurves = Curve37.6)
terrestrialfurnse=subset(furnishedse, is.na(furnishedse$Marine.protein)| Marine.protein=="0")
furnishedsecalterrestrial=calibrate(x=terrestrialfurnse$CRA, errors=terrestrialfurnse$Error, normalised=FALSE)
combfurnishedsedate=combine(furnishedsecal1.2, furnishedsecal2.4, furnishedsecal3.5, furnishedsecal4.7, furnishedsecal5.9, furnishedsecal6.7, furnishedsecal7.1, furnishedsecal8, furnishedsecal8.2, furnishedsecal8.3, furnishedsecal8.7, furnishedsecal9.4, furnishedsecal9.6, furnishedsecal10.1, furnishedsecal10.6, furnishedsecal11.8, furnishedsecal12.9, furnishedsecal13.9, furnishedsecal14.1, furnishedsecal14.6, furnishedsecal15.3, furnishedsecal16.5, furnishedsecal17.6, furnishedsecal18.8, furnishedsecal20, furnishedsecal21.2, furnishedsecal22.4, furnishedsecal23.5, furnishedsecal25.9, furnishedsecal27.1, furnishedsecal28.4, furnishedsecal32.7, furnishedsecal37.6, furnishedsecalterrestrial, fixIDs = TRUE)
binsfse=binPrep(sites = furnishedse$SiteName, ages=furnishedse$CRA, h=50)
samplefurnsemarine=sampleDates(combfurnishedsedate, bins=binsfse, nsim=1000, boot=TRUE)
kdemarinefurnishedse30=ckde(samplefurnsemarine, timeRange=c(1650,1050), bw=30)

#Furnished burials in the south and east, weighted by number of grave goods
weightedse=subset(weighted, Region=="SE")
weightedse1.2=subset(weightedse, Marine.protein==1.2)
weightedsecal1.2=calibrate(x=weightedse1.2$CRA, errors=weightedse1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
weightedse2.4=subset(weightedse, Marine.protein==2.4)
weightedsecal2.4=calibrate(x=weightedse2.4$CRA, errors=weightedse2.4$Error, normalised=FALSE, calCurves = Curve2.4)
weightedse3.5=subset(weightedse, Marine.protein==3.5)
weightedsecal3.5=calibrate(x=weightedse3.5$CRA, errors=weightedse3.5$Error, normalised=FALSE, calCurves = Curve3.5)
weightedse4.7=subset(weightedse, Marine.protein==4.7)
weightedsecal4.7=calibrate(x=weightedse4.7$CRA, errors=weightedse4.7$Error, normalised=FALSE, calCurves = Curve4.7)
weightedse5.9=subset(weightedse, Marine.protein==5.9)
weightedsecal5.9=calibrate(x=weightedse5.9$CRA, errors=weightedse5.9$Error, normalised=FALSE, calCurves = Curve5.9)
weightedse6.7=subset(weightedse, Marine.protein==6.7)
weightedsecal6.7=calibrate(x=weightedse6.7$CRA, errors=weightedse6.7$Error, normalised=FALSE, calCurves = Curve6.7)
weightedse7.1=subset(weightedse, Marine.protein==7.1)
weightedsecal7.1=calibrate(x=weightedse7.1$CRA, errors=weightedse7.1$Error, normalised=FALSE, calCurves = Curve7.1)
weightedse8=subset(weightedse, Marine.protein==8)
weightedsecal8=calibrate(x=weightedse8$CRA, errors=weightedse8$Error, normalised=FALSE, calCurves = Curve8)
weightedse8.2=subset(weightedse, Marine.protein==8.2)
weightedsecal8.2=calibrate(x=weightedse8.2$CRA, errors=weightedse8.2$Error, normalised=FALSE, calCurves = Curve8.2)
weightedse8.3=subset(weightedse, Marine.protein==8.3)
weightedsecal8.3=calibrate(x=weightedse8.3$CRA, errors=weightedse8.3$Error, normalised=FALSE, calCurves = Curve8.3)
weightedse8.7=subset(weightedse, Marine.protein==8.7)
weightedsecal8.7=calibrate(x=weightedse8.7$CRA, errors=weightedse8.7$Error, normalised=FALSE, calCurves = Curve8.7)
weightedse9.4=subset(weightedse, Marine.protein==9.4)
weightedsecal9.4=calibrate(x=weightedse9.4$CRA, errors=weightedse9.4$Error, normalised=FALSE, calCurves = Curve9.4)
weightedse9.6=subset(weightedse, Marine.protein==9.6)
weightedsecal9.6=calibrate(x=weightedse9.6$CRA, errors=weightedse9.6$Error, normalised=FALSE, calCurves = Curve9.6)
weightedse10.1=subset(weightedse, Marine.protein==10.1)
weightedsecal10.1=calibrate(x=weightedse10.1$CRA, errors=weightedse10.1$Error, normalised=FALSE, calCurves = Curve10.1)
weightedse10.6=subset(weightedse, Marine.protein==10.6)
weightedsecal10.6=calibrate(x=weightedse10.6$CRA, errors=weightedse10.6$Error, normalised=FALSE, calCurves = Curve10.6)
weightedse11.8=subset(weightedse, Marine.protein==11.8)
weightedsecal11.8=calibrate(x=weightedse11.8$CRA, errors=weightedse11.8$Error, normalised=FALSE, calCurves = Curve11.8)
weightedse12.9=subset(weightedse, Marine.protein==12.9)
weightedsecal12.9=calibrate(x=weightedse12.9$CRA, errors=weightedse12.9$Error, normalised=FALSE, calCurves = Curve12.9)
weightedse13.9=subset(weightedse, Marine.protein==13.9)
weightedsecal13.9=calibrate(x=weightedse13.9$CRA, errors=weightedse13.9$Error, normalised=FALSE, calCurves = Curve13.9)
weightedse14.1=subset(weightedse, Marine.protein==14.1)
weightedsecal14.1=calibrate(x=weightedse14.1$CRA, errors=weightedse14.1$Error, normalised=FALSE, calCurves = Curve14.1)
weightedse14.6=subset(weightedse, Marine.protein==14.6)
weightedsecal14.6=calibrate(x=weightedse14.6$CRA, errors=weightedse14.6$Error, normalised=FALSE, calCurves = Curve14.6)
weightedse15.3=subset(weightedse, Marine.protein==15.3)
weightedsecal15.3=calibrate(x=weightedse15.3$CRA, errors=weightedse15.3$Error, normalised=FALSE, calCurves = Curve15.3)
weightedse16.5=subset(weightedse, Marine.protein==16.5)
weightedsecal16.5=calibrate(x=weightedse16.5$CRA, errors=weightedse16.5$Error, normalised=FALSE, calCurves = Curve16.5)
weightedse17.6=subset(weightedse, Marine.protein==17.6)
weightedsecal17.6=calibrate(x=weightedse17.6$CRA, errors=weightedse17.6$Error, normalised=FALSE, calCurves = Curve17.6)
weightedse18.8=subset(weightedse, Marine.protein==18.8)
weightedsecal18.8=calibrate(x=weightedse18.8$CRA, errors=weightedse18.8$Error, normalised=FALSE, calCurves = Curve18.8)
weightedse20=subset(weightedse, Marine.protein==20)
weightedsecal20=calibrate(x=weightedse20$CRA, errors=weightedse20$Error, normalised=FALSE, calCurves = Curve20)
weightedse21.2=subset(weightedse, Marine.protein==21.2)
weightedsecal21.2=calibrate(x=weightedse21.2$CRA, errors=weightedse21.2$Error, normalised=FALSE, calCurves = Curve21.2)
weightedse22.4=subset(weightedse, Marine.protein==22.4)
weightedsecal22.4=calibrate(x=weightedse22.4$CRA, errors=weightedse22.4$Error, normalised=FALSE, calCurves = Curve22.4)
weightedse23.5=subset(weightedse, Marine.protein==23.5)
weightedsecal23.5=calibrate(x=weightedse23.5$CRA, errors=weightedse23.5$Error, normalised=FALSE, calCurves = Curve23.5)
weightedse25.9=subset(weightedse, Marine.protein==25.9)
weightedsecal25.9=calibrate(x=weightedse25.9$CRA, errors=weightedse25.9$Error, normalised=FALSE, calCurves = Curve25.9)
weightedse27.1=subset(weightedse, Marine.protein==27.1)
weightedsecal27.1=calibrate(x=weightedse27.1$CRA, errors=weightedse27.1$Error, normalised=FALSE, calCurves = Curve27.1)
weightedse28.4=subset(weightedse, Marine.protein==28.4)
weightedsecal28.4=calibrate(x=weightedse28.4$CRA, errors=weightedse28.4$Error, normalised=FALSE, calCurves = Curve28.4)
weightedse32.7=subset(weightedse, Marine.protein==32.7)
weightedsecal32.7=calibrate(x=weightedse32.7$CRA, errors=weightedse32.7$Error, normalised=FALSE, calCurves = Curve32.7)
weightedse37.6=subset(weightedse, Marine.protein==37.6)
weightedsecal37.6=calibrate(x=weightedse37.6$CRA, errors=weightedse37.6$Error, normalised=FALSE, calCurves = Curve37.6)
terrestrialweightedse=subset(weightedse, is.na(weightedse$Marine.protein)| Marine.protein=="0")
weightedsecalterrestrial=calibrate(x=terrestrialweightedse$CRA, errors=terrestrialweightedse$Error, normalised=FALSE)
combweightedsedate=combine(weightedsecal1.2, weightedsecal2.4, weightedsecal3.5, weightedsecal4.7, weightedsecal5.9, weightedsecal6.7, weightedsecal7.1, weightedsecal8, weightedsecal8.2, weightedsecal8.3, weightedsecal8.7, weightedsecal9.4, weightedsecal9.6, weightedsecal10.1, weightedsecal10.6, weightedsecal11.8, weightedsecal12.9, weightedsecal13.9, weightedsecal14.1, weightedsecal14.6, weightedsecal15.3, weightedsecal16.5, weightedsecal17.6, weightedsecal18.8, weightedsecal20, weightedsecal21.2, weightedsecal22.4, weightedsecal23.5, weightedsecal25.9, weightedsecal27.1, weightedsecal28.4, weightedsecal32.7, weightedsecal37.6, weightedsecalterrestrial, fixIDs = TRUE)
sampleweightedsemarine=sampleDates(combweightedsedate, nsim=1000, boot=TRUE)
kdemarineweightedse30=ckde(sampleweightedsemarine, timeRange=c(1650,1050), bw=30)

#plot figure 9a
plot(kdemarinefurnishedse30, calendar="BCAD",ylim=c(0,0.004)) + title("Furnished graves in the south and east, n=281")
par(new=TRUE)
plot(kdemarineweightedse30, calendar="BCAD",ylim=c(0,0.004), fill.col="#FFFF007D")
legend("topright", legend=c("Furnished graves in the south and east","Weighted by number of grave goods"), fill=c("lightgrey", "#FFFF007D"))

#unfurnished burials south and east
seunfurnished=subset(unfurnished, Region=="SE")
seunfurnished1.2=subset(seunfurnished, Marine.protein==1.2)
seunfurnishedcal1.2=calibrate(x=seunfurnished1.2$CRA, errors=seunfurnished1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
seunfurnished2.4=subset(seunfurnished, Marine.protein==2.4)
seunfurnishedcal2.4=calibrate(x=seunfurnished2.4$CRA, errors=seunfurnished2.4$Error, normalised=FALSE, calCurves = Curve2.4)
seunfurnished3.5=subset(seunfurnished, Marine.protein==3.5)
seunfurnishedcal3.5=calibrate(x=seunfurnished3.5$CRA, errors=seunfurnished3.5$Error, normalised=FALSE, calCurves = Curve3.5)
seunfurnished4.7=subset(seunfurnished, Marine.protein==4.7)
seunfurnishedcal4.7=calibrate(x=seunfurnished4.7$CRA, errors=seunfurnished4.7$Error, normalised=FALSE, calCurves = Curve4.7)
seunfurnished5.9=subset(seunfurnished, Marine.protein==5.9)
seunfurnishedcal5.9=calibrate(x=seunfurnished5.9$CRA, errors=seunfurnished5.9$Error, normalised=FALSE, calCurves = Curve5.9)
seunfurnished6.4=subset(seunfurnished, Marine.protein==6.4)
seunfurnishedcal6.4=calibrate(x=seunfurnished6.4$CRA, errors=seunfurnished6.4$Error, normalised=FALSE, calCurves = Curve6.4)
seunfurnished7.1=subset(seunfurnished, Marine.protein==7.1)
seunfurnishedcal7.1=calibrate(x=seunfurnished7.1$CRA, errors=seunfurnished7.1$Error, normalised=FALSE, calCurves = Curve7.1)
seunfurnished8.2=subset(seunfurnished, Marine.protein==8.2)
seunfurnishedcal8.2=calibrate(x=seunfurnished8.2$CRA, errors=seunfurnished8.2$Error, normalised=FALSE, calCurves = Curve8.2)
seunfurnished9.4=subset(seunfurnished, Marine.protein==9.4)
seunfurnishedcal9.4=calibrate(x=seunfurnished9.4$CRA, errors=seunfurnished9.4$Error, normalised=FALSE, calCurves = Curve9.4)
seunfurnished10=subset(seunfurnished, Marine.protein==10)
seunfurnishedcal10=calibrate(x=seunfurnished10$CRA, errors=seunfurnished10$Error, normalised=FALSE, calCurves = Curve10)
seunfurnished10.1=subset(seunfurnished, Marine.protein==10.1)
seunfurnishedcal10.1=calibrate(x=seunfurnished10.1$CRA, errors=seunfurnished10.1$Error, normalised=FALSE, calCurves = Curve10.1)
seunfurnished10.5=subset(seunfurnished, Marine.protein==10.5)
seunfurnishedcal10.5=calibrate(x=seunfurnished10.5$CRA, errors=seunfurnished10.5$Error, normalised=FALSE, calCurves = Curve10.5)
seunfurnished10.6=subset(seunfurnished, Marine.protein==10.6)
seunfurnishedcal10.6=calibrate(x=seunfurnished10.6$CRA, errors=seunfurnished10.6$Error, normalised=FALSE, calCurves = Curve10.6)
seunfurnished11.5=subset(seunfurnished, Marine.protein==11.5)
seunfurnishedcal11.5=calibrate(x=seunfurnished11.5$CRA, errors=seunfurnished11.5$Error, normalised=FALSE, calCurves=Curve11.5)
seunfurnished11.8=subset(seunfurnished, Marine.protein==11.8)
seunfurnishedcal11.8=calibrate(x=seunfurnished11.8$CRA, errors=seunfurnished11.8$Error, normalised=FALSE, calCurves = Curve11.8)
seunfurnished12.2=subset(seunfurnished, Marine.protein==12.2)
seunfurnishedcal12.2=calibrate(x=seunfurnished12.2$CRA, errors=seunfurnished12.2$Error, normalised=FALSE, calCurves = Curve12.2)
seunfurnished12.5=subset(seunfurnished, Marine.protein==12.5)
seunfurnishedcal12.5=calibrate(x=seunfurnished12.5$CRA, errors=seunfurnished12.5$Error, normalised=FALSE, calCurves = Curve12.5)
seunfurnished12.9=subset(seunfurnished, Marine.protein==12.9)
seunfurnishedcal12.9=calibrate(x=seunfurnished12.9$CRA, errors=seunfurnished12.9$Error, normalised=FALSE, calCurves = Curve12.9)
seunfurnished14.1=subset(seunfurnished, Marine.protein==14.1)
seunfurnishedcal14.1=calibrate(x=seunfurnished14.1$CRA, errors=seunfurnished14.1$Error, normalised=FALSE, calCurves = Curve14.1)
seunfurnished15.3=subset(seunfurnished, Marine.protein==15.3)
seunfurnishedcal15.3=calibrate(x=seunfurnished15.3$CRA, errors=seunfurnished15.3$Error, normalised=FALSE, calCurves = Curve15.3)
seunfurnished16.5=subset(seunfurnished, Marine.protein==16.5)
seunfurnishedcal16.5=calibrate(x=seunfurnished16.5$CRA, errors=seunfurnished16.5$Error, normalised=FALSE, calCurves = Curve16.5)
seunfurnished17.6=subset(seunfurnished, Marine.protein==17.6)
seunfurnishedcal17.6=calibrate(x=seunfurnished17.6$CRA, errors=seunfurnished17.6$Error, normalised=FALSE, calCurves = Curve17.6)
seunfurnished17.8=subset(seunfurnished, Marine.protein==17.8)
seunfurnishedcal17.8=calibrate(x=seunfurnished17.8$CRA, errors=seunfurnished17.8$Error, normalised=FALSE, calCurves = Curve17.8)
seunfurnished18.6=subset(seunfurnished, Marine.protein==18.6)
seunfurnishedcal18.6=calibrate(x=seunfurnished18.6$CRA, errors=seunfurnished18.6$Error, normalised=FALSE, calCurves = Curve18.6)
seunfurnished18.8=subset(seunfurnished, Marine.protein==18.8)
seunfurnishedcal18.8=calibrate(x=seunfurnished18.8$CRA, errors=seunfurnished18.8$Error, normalised=FALSE, calCurves = Curve18.8)
seunfurnished18.9=subset(seunfurnished, Marine.protein==18.9)
seunfurnishedcal18.9=calibrate(x=seunfurnished18.9$CRA, errors=seunfurnished18.9$Error, normalised=FALSE, calCurves = Curve18.9)
seunfurnished20=subset(seunfurnished, Marine.protein==20)
seunfurnishedcal20=calibrate(x=seunfurnished20$CRA, errors=seunfurnished20$Error, normalised=FALSE, calCurves = Curve20)
seunfurnished20.4=subset(seunfurnished, Marine.protein==20.4)
seunfurnishedcal20.4=calibrate(x=seunfurnished20.4$CRA, errors=seunfurnished20.4$Error, normalised=FALSE, calCurves = Curve20.4)
seunfurnished21.2=subset(seunfurnished, Marine.protein==21.2)
seunfurnishedcal21.2=calibrate(x=seunfurnished21.2$CRA, errors=seunfurnished21.2$Error, normalised=FALSE, calCurves = Curve21.2)
seunfurnished21.9=subset(seunfurnished, Marine.protein==21.9)
seunfurnishedcal21.9=calibrate(x=seunfurnished21.9$CRA, errors=seunfurnished21.9$Error, normalised=FALSE, calCurves = Curve21.9)
seunfurnished22.4=subset(seunfurnished, Marine.protein==22.4)
seunfurnishedcal22.4=calibrate(x=seunfurnished22.4$CRA, errors=seunfurnished22.4$Error, normalised=FALSE, calCurves = Curve22.4)
seunfurnished23.5=subset(seunfurnished, Marine.protein==23.5)
seunfurnishedcal23.5=calibrate(x=seunfurnished23.5$CRA, errors=seunfurnished23.5$Error, normalised=FALSE, calCurves = Curve23.5)
seunfurnished24.7=subset(seunfurnished, Marine.protein==24.7)
seunfurnishedcal24.7=calibrate(x=seunfurnished24.7$CRA, errors=seunfurnished24.7$Error, normalised=FALSE, calCurves = Curve24.7)
seunfurnished25.9=subset(seunfurnished, Marine.protein==25.9)
seunfurnishedcal25.9=calibrate(x=seunfurnished25.9$CRA, errors=seunfurnished25.9$Error, normalised=FALSE, calCurves = Curve25.9)
seunfurnished27.1=subset(seunfurnished, Marine.protein==27.1)
seunfurnishedcal27.1=calibrate(x=seunfurnished27.1$CRA, errors=seunfurnished27.1$Error, normalised=FALSE, calCurves = Curve27.1)
seunfurnished28.2=subset(seunfurnished, Marine.protein==28.2)
seunfurnishedcal28.2=calibrate(x=seunfurnished28.2$CRA, errors=seunfurnished28.2$Error, normalised=FALSE, calCurves = Curve28.2)
seunfurnished29.4=subset(seunfurnished, Marine.protein==29.4)
seunfurnishedcal29.4=calibrate(x=seunfurnished29.4$CRA, errors=seunfurnished29.4$Error, normalised=FALSE, calCurves = Curve29.4)
seunfurnished30.6=subset(seunfurnished, Marine.protein==30.6)
seunfurnishedcal30.6=calibrate(x=seunfurnished30.6$CRA, errors=seunfurnished30.6$Error, normalised=FALSE, calCurves = Curve30.6)
seunfurnished31.8=subset(seunfurnished, Marine.protein==31.8)
seunfurnishedcal31.8=calibrate(x=seunfurnished31.8$CRA, errors=seunfurnished31.8$Error, normalised=FALSE, calCurves = Curve31.8)
seunfurnished32=subset(seunfurnished, Marine.protein==32)
seunfurnishedcal32=calibrate(x=seunfurnished32$CRA, errors=seunfurnished32$Error, normalised=FALSE, calCurves = Curve32)
seunfurnished32.9=subset(seunfurnished, Marine.protein==32.9)
seunfurnishedcal32.9=calibrate(x=seunfurnished32.9$CRA, errors=seunfurnished32.9$Error, normalised=FALSE, calCurves = Curve32.9)
seunfurnished37.6=subset(seunfurnished, Marine.protein==37.6)
seunfurnishedcal37.6=calibrate(x=seunfurnished37.6$CRA, errors=seunfurnished37.6$Error, normalised=FALSE, calCurves = Curve37.6)
terrestrialunfse=subset(seunfurnished, is.na(seunfurnished$Marine.protein)| Marine.protein=="0")
seunfurnishedcalterrestrial=calibrate(x=terrestrialunfse$CRA, errors=terrestrialunfse$Error, normalised=FALSE)
combdateunfse=combine(seunfurnishedcal1.2, seunfurnishedcal2.4, seunfurnishedcal3.5, seunfurnishedcal4.7, seunfurnishedcal5.9, seunfurnishedcal6.4, seunfurnishedcal7.1, seunfurnishedcal8.2, seunfurnishedcal9.4, seunfurnishedcal10, seunfurnishedcal10.1, seunfurnishedcal10.5, seunfurnishedcal10.6, seunfurnishedcal11.5, seunfurnishedcal11.8, seunfurnishedcal12.2, seunfurnishedcal12.5, seunfurnishedcal12.9, seunfurnishedcal14.1, seunfurnishedcal15.3, seunfurnishedcal16.5, seunfurnishedcal17.6, seunfurnishedcal17.8, seunfurnishedcal18.6, seunfurnishedcal18.8, seunfurnishedcal18.9, seunfurnishedcal20, seunfurnishedcal20.4, seunfurnishedcal21.2, seunfurnishedcal21.9, seunfurnishedcal22.4, seunfurnishedcal23.5, seunfurnishedcal24.7, seunfurnishedcal25.9, seunfurnishedcal27.1, seunfurnishedcal28.2, seunfurnishedcal29.4, seunfurnishedcal30.6, seunfurnishedcal31.8,seunfurnishedcal32, seunfurnishedcal32.9, seunfurnishedcal37.6, seunfurnishedcalterrestrial, fixIDs = TRUE)
binsuse=binPrep(sites = seunfurnished$SiteName, ages=seunfurnished$CRA, h=50)
sampleseunfurnishedmarine=sampleDates(combdateunfse, bins=binsuse, nsim=1000, boot=TRUE)
kdeseunfurnished30marine=ckde(sampleseunfurnishedmarine, timeRange=c(1650,1050), bw=30)
plot(kdeseunfurnished30marine, calendar="BCAD") +title("Unfurnished graves in the south and east, n=367")

#Richly furnished graves
rich=subset(alldates, Rich.materials=="Y")
rich2.4=subset(rich, Marine.protein==2.4)
richcal2.4=calibrate(x=rich2.4$CRA, errors=rich2.4$Error, normalised=FALSE, calCurves = Curve2.4)
rich4.7=subset(rich, Marine.protein==4.7)
richcal4.7=calibrate(x=rich4.7$CRA, errors=rich4.7$Error, normalised=FALSE, calCurves = Curve4.7)
rich5.9=subset(rich, Marine.protein==5.9)
richcal5.9=calibrate(x=rich5.9$CRA, errors=rich5.9$Error, normalised=FALSE, calCurves = Curve5.9)
rich7.1=subset(rich, Marine.protein==7.1)
richcal7.1=calibrate(x=rich7.1$CRA, errors=rich7.1$Error, normalised=FALSE, calCurves = Curve7.1)
rich8.2=subset(rich, Marine.protein==8.2)
richcal8.2=calibrate(x=rich8.2$CRA, errors=rich8.2$Error, normalised=FALSE, calCurves = Curve8.2)
rich8.3=subset(rich, Marine.protein==8.3)
richcal8.3=calibrate(x=rich8.3$CRA, errors=rich8.3$Error, normalised=FALSE, calCurves = Curve8.3)
rich9.4=subset(rich, Marine.protein==9.4)
richcal9.4=calibrate(x=rich9.4$CRA, errors=rich9.4$Error, normalised=FALSE, calCurves = Curve9.4)
rich10.6=subset(rich, Marine.protein==10.6)
richcal10.6=calibrate(x=rich10.6$CRA, errors=rich10.6$Error, normalised=FALSE, calCurves = Curve10.6)
rich11.8=subset(rich, Marine.protein==11.8)
richcal11.8=calibrate(x=rich11.8$CRA, errors=rich11.8$Error, normalised=FALSE, calCurves = Curve11.8)
rich12.9=subset(rich, Marine.protein==12.9)
richcal12.9=calibrate(x=rich12.9$CRA, errors=rich12.9$Error, normalised=FALSE, calCurves = Curve12.9)
rich14.1=subset(rich, Marine.protein==14.1)
richcal14.1=calibrate(x=rich14.1$CRA, errors=rich14.1$Error, normalised=FALSE, calCurves = Curve14.1)
rich15.3=subset(rich, Marine.protein==15.3)
richcal15.3=calibrate(x=rich15.3$CRA, errors=rich15.3$Error, normalised=FALSE, calCurves = Curve15.3)
rich16.5=subset(rich, Marine.protein==16.5)
richcal16.5=calibrate(x=rich16.5$CRA, errors=rich16.5$Error, normalised=FALSE, calCurves = Curve16.5)
rich17.6=subset(rich, Marine.protein==17.6)
richcal17.6=calibrate(x=rich17.6$CRA, errors=rich17.6$Error, normalised=FALSE, calCurves = Curve17.6)
rich21.2=subset(rich, Marine.protein==21.2)
richcal21.2=calibrate(x=rich21.2$CRA, errors=rich21.2$Error, normalised=FALSE, calCurves = Curve21.2)
richterrestrial=subset(rich, is.na(rich$Marine.protein)| Marine.protein=="0")
richcalterrestrial=calibrate(x=richterrestrial$CRA, errors=richterrestrial$Error, normalised=FALSE)
combrichdate=combine(richcal2.4, richcal4.7, richcal5.9, richcal7.1, richcal8.2, richcal8.3, richcal9.4, richcal10.6, richcal11.8, richcal12.9, richcal14.1, richcal15.3, richcal16.5, richcal17.6, richcal21.2, richcalterrestrial, fixIDs = TRUE)
binsr=binPrep(sites= rich$SiteName, ages= rich$CRA, h=50)
samplerichmarine=sampleDates(combrichdate, bins=binsr, nsim=1000, boot=TRUE)
kdemarinerich30=ckde(samplerichmarine, timeRange=c(1650,1050), bw=30)
plot(kdemarinerich30, calendar="BCAD", ylim=c(0,0.005)) +title("Richly furnished graves, n=49")

#Isolated burials
isoburial=subset(alldates, Burial.Type=="Isolated burial")
isoburial2.4=subset(isoburial, Marine.protein==2.4)
isoburialcal2.4=calibrate(x=isoburial2.4$CRA, errors=isoburial2.4$Error, normalised=FALSE, calCurves = Curve2.4)
isoburial3.5=subset(isoburial, Marine.protein==3.5)
isoburialcal3.5=calibrate(x=isoburial3.5$CRA, errors=isoburial3.5$Error, normalised=FALSE, calCurves = Curve3.5)
isoburial4.7=subset(isoburial, Marine.protein==4.7)
isoburialcal4.7=calibrate(x=isoburial4.7$CRA, errors=isoburial4.7$Error, normalised=FALSE, calCurves = Curve4.7)
isoburial5.9=subset(isoburial, Marine.protein==5.9)
isoburialcal5.9=calibrate(x=isoburial5.9$CRA, errors=isoburial5.9$Error, normalised=FALSE, calCurves = Curve5.9)
isoburial7.1=subset(isoburial, Marine.protein==7.1)
isoburialcal7.1=calibrate(x=isoburial7.1$CRA, errors=isoburial7.1$Error, normalised=FALSE, calCurves = Curve7.1)
isoburial8=subset(isoburial, Marine.protein==8)
isoburialcal8=calibrate(x=isoburial8$CRA, errors=isoburial8$Error, normalised=FALSE, calCurves = Curve8)
isoburial8.2=subset(isoburial, Marine.protein==8.2)
isoburialcal8.2=calibrate(x=isoburial8.2$CRA, errors=isoburial8.2$Error, normalised=FALSE, calCurves = Curve8.2)
isoburial9.4=subset(isoburial, Marine.protein==9.4)
isoburialcal9.4=calibrate(x=isoburial9.4$CRA, errors=isoburial9.4$Error, normalised=FALSE, calCurves = Curve9.4)
isoburial10=subset(isoburial, Marine.protein==10)
isoburialcal10=calibrate(x=isoburial10$CRA, errors=isoburial10$Error, normalised=FALSE, calCurves = Curve10)
isoburial10.6=subset(isoburial, Marine.protein==10.6)
isoburialcal10.6=calibrate(x=isoburial10.6$CRA, errors=isoburial10.6$Error, normalised=FALSE, calCurves = Curve10.6)
isoburial11.5=subset(isoburial, Marine.protein==11.5)
isoburialcal11.5=calibrate(x=isoburial11.5$CRA, errors=isoburial11.5$Error, normalised=FALSE, calCurves=Curve11.5)
isoburial11.8=subset(isoburial, Marine.protein==11.8)
isoburialcal11.8=calibrate(x=isoburial11.8$CRA, errors=isoburial11.8$Error, normalised=FALSE, calCurves = Curve11.8)
isoburial12.9=subset(isoburial, Marine.protein==12.9)
isoburialcal12.9=calibrate(x=isoburial12.9$CRA, errors=isoburial12.9$Error, normalised=FALSE, calCurves = Curve12.9)
isoburial13.9=subset(isoburial, Marine.protein==13.9)
isoburialcal13.9=calibrate(x=isoburial13.9$CRA, errors=isoburial13.9$Error, normalised=FALSE, calCurves = Curve13.9)
isoburial14.1=subset(isoburial, Marine.protein==14.1)
isoburialcal14.1=calibrate(x=isoburial14.1$CRA, errors=isoburial14.1$Error, normalised=FALSE, calCurves = Curve14.1)
isoburial16.5=subset(isoburial, Marine.protein==16.5)
isoburialcal16.5=calibrate(x=isoburial16.5$CRA, errors=isoburial16.5$Error, normalised=FALSE, calCurves = Curve16.5)
isoburial20=subset(isoburial, Marine.protein==20)
isoburialcal20=calibrate(x=isoburial20$CRA, errors=isoburial20$Error, normalised=FALSE, calCurves = Curve20)
isoburial20.4=subset(isoburial, Marine.protein==20.4)
isoburialcal20.4=calibrate(x=isoburial20.4$CRA, errors=isoburial20.4$Error, normalised=FALSE, calCurves = Curve20.4)
isoburial21.2=subset(isoburial, Marine.protein==21.2)
isoburialcal21.2=calibrate(x=isoburial21.2$CRA, errors=isoburial21.2$Error, normalised=FALSE, calCurves = Curve21.2)
isoburial23.5=subset(isoburial, Marine.protein==23.5)
isoburialcal23.5=calibrate(x=isoburial23.5$CRA, errors=isoburial23.5$Error, normalised=FALSE, calCurves = Curve23.5)
isoburial28.2=subset(isoburial, Marine.protein==28.2)
isoburialcal28.2=calibrate(x=isoburial28.2$CRA, errors=isoburial28.2$Error, normalised=FALSE, calCurves = Curve28.2)
isoburial29.4=subset(isoburial, Marine.protein==29.4)
isoburialcal29.4=calibrate(x=isoburial29.4$CRA, errors=isoburial29.4$Error, normalised=FALSE, calCurves = Curve29.4)
isoburial32.9=subset(isoburial, Marine.protein==32.9)
isoburialcal32.9=calibrate(x=isoburial32.9$CRA, errors=isoburial32.9$Error, normalised=FALSE, calCurves = Curve32.9)
isoterrestrial=subset(isoburial, is.na(isoburial$Marine.protein)| Marine.protein=="0")
isoburialcalterrestrial=calibrate(x=isoterrestrial$CRA, errors=isoterrestrial$Error, normalised=FALSE)
combdateiso=combine(isoburialcal2.4, isoburialcal3.5, isoburialcal4.7, isoburialcal5.9, isoburialcal7.1, isoburialcal8, isoburialcal8.2, isoburialcal9.4, isoburialcal10, isoburialcal10.6, isoburialcal11.5, isoburialcal11.8, isoburialcal12.9, isoburialcal13.9, isoburialcal14.1, isoburialcal16.5, isoburialcal20, isoburialcal20.4, isoburialcal21.2, isoburialcal23.5, isoburialcal28.2, isoburialcal29.4, isoburialcal32.9, isoburialcalterrestrial, fixIDs = TRUE)
sampleisoburialmarine=sampleDates(combdateiso, nsim=1000, boot=TRUE)
kdeisoburialmarine=ckde(sampleisoburialmarine, timeRange = c(1650, 1050), bw=30)
plot(kdeisoburialmarine, calendar = "BCAD", ylim=c(0,0.0025)) + title("Isolated burials, n=72")

#Small groups of burials
smallgroup=subset(alldates, Burial.Type=="Small group of burials")
smallgroup3.5=subset(smallgroup, Marine.protein==3.5)
smallgroupcal3.5=calibrate(x=smallgroup3.5$CRA, errors=smallgroup3.5$Error, normalised=FALSE, calCurves = Curve3.5)
smallgroup4.7=subset(smallgroup, Marine.protein==4.7)
smallgroupcal4.7=calibrate(x=smallgroup4.7$CRA, errors=smallgroup4.7$Error, normalised=FALSE, calCurves = Curve4.7)
smallgroup5.9=subset(smallgroup, Marine.protein==5.9)
smallgroupcal5.9=calibrate(x=smallgroup5.9$CRA, errors=smallgroup5.9$Error, normalised=FALSE, calCurves = Curve5.9)
smallgroup7.1=subset(smallgroup, Marine.protein==7.1)
smallgroupcal7.1=calibrate(x=smallgroup7.1$CRA, errors=smallgroup7.1$Error, normalised=FALSE, calCurves = Curve7.1)
smallgroup8.2=subset(smallgroup, Marine.protein==8.2)
smallgroupcal8.2=calibrate(x=smallgroup8.2$CRA, errors=smallgroup8.2$Error, normalised=FALSE, calCurves = Curve8.2)
smallgroup9.4=subset(smallgroup, Marine.protein==9.4)
smallgroupcal9.4=calibrate(x=smallgroup9.4$CRA, errors=smallgroup9.4$Error, normalised=FALSE, calCurves = Curve9.4)
smallgroup10.1=subset(smallgroup, Marine.protein==10.1)
smallgroupcal10.1=calibrate(x=smallgroup10.1$CRA, errors=smallgroup10.1$Error, normalised=FALSE, calCurves = Curve10.1)
smallgroup10.6=subset(smallgroup, Marine.protein==10.6)
smallgroupcal10.6=calibrate(x=smallgroup10.6$CRA, errors=smallgroup10.6$Error, normalised=FALSE, calCurves = Curve10.6)
smallgroup11.8=subset(smallgroup, Marine.protein==11.8)
smallgroupcal11.8=calibrate(x=smallgroup11.8$CRA, errors=smallgroup11.8$Error, normalised=FALSE, calCurves = Curve11.8)
smallgroup12.9=subset(smallgroup, Marine.protein==12.9)
smallgroupcal12.9=calibrate(x=smallgroup12.9$CRA, errors=smallgroup12.9$Error, normalised=FALSE, calCurves = Curve12.9)
smallgroup16.5=subset(smallgroup, Marine.protein==16.5)
smallgroupcal16.5=calibrate(x=smallgroup16.5$CRA, errors=smallgroup16.5$Error, normalised=FALSE, calCurves = Curve16.5)
smallgroup17.6=subset(smallgroup, Marine.protein==17.6)
smallgroupcal17.6=calibrate(x=smallgroup17.6$CRA, errors=smallgroup17.6$Error, normalised=FALSE, calCurves = Curve17.6)
smallgroup17.8=subset(smallgroup, Marine.protein==17.8)
smallgroupcal17.8=calibrate(x=smallgroup17.8$CRA, errors=smallgroup17.8$Error, normalised=FALSE, calCurves = Curve17.8)
smallgroup18.6=subset(smallgroup, Marine.protein==18.6)
smallgroupcal18.6=calibrate(x=smallgroup18.6$CRA, errors=smallgroup18.6$Error, normalised=FALSE, calCurves = Curve18.6)
smallgroup20=subset(smallgroup, Marine.protein==20)
smallgroupcal20=calibrate(x=smallgroup20$CRA, errors=smallgroup20$Error, normalised=FALSE, calCurves = Curve20)
smallgroup21.9=subset(smallgroup, Marine.protein==21.9)
smallgroupcal21.9=calibrate(x=smallgroup21.9$CRA, errors=smallgroup21.9$Error, normalised=FALSE, calCurves = Curve21.9)
smallgroup22.4=subset(smallgroup, Marine.protein==22.4)
smallgroupcal22.4=calibrate(x=smallgroup22.4$CRA, errors=smallgroup22.4$Error, normalised=FALSE, calCurves = Curve22.4)
smallgroup23.5=subset(smallgroup, Marine.protein==23.5)
smallgroupcal23.5=calibrate(x=smallgroup23.5$CRA, errors=smallgroup23.5$Error, normalised=FALSE, calCurves = Curve23.5)
smallgroup24.7=subset(smallgroup, Marine.protein==24.7)
smallgroupcal24.7=calibrate(x=smallgroup24.7$CRA, errors=smallgroup24.7$Error, normalised=FALSE, calCurves = Curve24.7)
smallgroup25.9=subset(smallgroup, Marine.protein==25.9)
smallgroupcal25.9=calibrate(x=smallgroup25.9$CRA, errors=smallgroup25.9$Error, normalised=FALSE, calCurves = Curve25.9)
smallgroup28.2=subset(smallgroup, Marine.protein==28.2)
smallgroupcal28.2=calibrate(x=smallgroup28.2$CRA, errors=smallgroup28.2$Error, normalised=FALSE, calCurves = Curve28.2)
smallgroup38.8=subset(smallgroup, Marine.protein==38.8)
smallgroupcal38.8=calibrate(x=smallgroup38.8$CRA, errors=smallgroup38.8$Error, normalised=FALSE, calCurves = Curve38.8)
smallgroupterrestrial=subset(smallgroup, is.na(smallgroup$Marine.protein)| Marine.protein=="0")
smallgroupcalterrestrial=calibrate(x=smallgroupterrestrial$CRA, errors=smallgroupterrestrial$Error, normalised=FALSE)
combdatesmallgroup=combine(smallgroupcal3.5, smallgroupcal4.7, smallgroupcal5.9, smallgroupcal7.1, smallgroupcal8.2, smallgroupcal9.4, smallgroupcal10.1, smallgroupcal10.6, smallgroupcal11.8, smallgroupcal12.9, smallgroupcal16.5, smallgroupcal17.6, smallgroupcal17.8, smallgroupcal18.6, smallgroupcal20, smallgroupcal21.9, smallgroupcal22.4, smallgroupcal23.5, smallgroupcal24.7, smallgroupcal25.9, smallgroupcal28.2, smallgroupcal38.8, smallgroupcalterrestrial, fixIDs = TRUE)
binssgc=binPrep(sites = smallgroup$SiteName, ages = smallgroup$CRA, h=50)
samplesmallgroupmarine=sampleDates(combdatesmallgroup, bins=binssgc, nsim=1000, boot=TRUE)
kdesmallgroupmarine=ckde(samplesmallgroupmarine, timeRange = c(1650, 1050), bw=30)
plot(kdesmallgroupmarine, calendar="BCAD", ylim=c(0,0.003)) + title("Small groups of burials, n=81")

#Large cemeteries
cemeteries=subset(alldates, Burial.Type== "Unfurnished cemetery" | Burial.Type=="Furnished cemetery")
cemeteries1.2=subset(cemeteries, Marine.protein==1.2)
cemeteriescal1.2=calibrate(x=cemeteries1.2$CRA, errors=cemeteries1.2$CRA, normalised=FALSE, calCurves=Curve1.2)
cemeteries2.4=subset(cemeteries, Marine.protein==2.4)
cemeteriescal2.4=calibrate(x=cemeteries2.4$CRA, errors=cemeteries2.4$Error, normalised=FALSE, calCurves = Curve2.4)
cemeteries3.5=subset(cemeteries, Marine.protein==3.5)
cemeteriescal3.5=calibrate(x=cemeteries3.5$CRA, errors=cemeteries3.5$Error, normalised=FALSE, calCurves = Curve3.5)
cemeteries4.7=subset(cemeteries, Marine.protein==4.7)
cemeteriescal4.7=calibrate(x=cemeteries4.7$CRA, errors=cemeteries4.7$Error, normalised=FALSE, calCurves = Curve4.7)
cemeteries5.9=subset(cemeteries, Marine.protein==5.9)
cemeteriescal5.9=calibrate(x=cemeteries5.9$CRA, errors=cemeteries5.9$Error, normalised=FALSE, calCurves = Curve5.9)
cemeteries6.4=subset(cemeteries, Marine.protein==6.4)
cemeteriescal6.4=calibrate(x=cemeteries6.4$CRA, errors=cemeteries6.4$Error, normalised=FALSE, calCurves = Curve6.4)
cemeteries6.7=subset(cemeteries, Marine.protein==6.7)
cemeteriescal6.7=calibrate(x=cemeteries6.7$CRA, errors=cemeteries6.7$Error, normalised=FALSE, calCurves = Curve6.7)
cemeteries7.1=subset(cemeteries, Marine.protein==7.1)
cemeteriescal7.1=calibrate(x=cemeteries7.1$CRA, errors=cemeteries7.1$Error, normalised=FALSE, calCurves = Curve7.1)
cemeteries8.2=subset(cemeteries, Marine.protein==8.2)
cemeteriescal8.2=calibrate(x=cemeteries8.2$CRA, errors=cemeteries8.2$Error, normalised=FALSE, calCurves = Curve8.2)
cemeteries8.3=subset(cemeteries, Marine.protein==8.3)
cemeteriescal8.3=calibrate(x=cemeteries8.3$CRA, errors=cemeteries8.3$Error, normalised=FALSE, calCurves = Curve8.3)
cemeteries8.7=subset(cemeteries, Marine.protein==8.7)
cemeteriescal8.7=calibrate(x=cemeteries8.7$CRA, errors=cemeteries8.7$Error, normalised=FALSE, calCurves = Curve8.7)
cemeteries9.4=subset(cemeteries, Marine.protein==9.4)
cemeteriescal9.4=calibrate(x=cemeteries9.4$CRA, errors=cemeteries9.4$Error, normalised=FALSE, calCurves = Curve9.4)
cemeteries9.5=subset(cemeteries, Marine.protein==9.5)
cemeteriescal9.5=calibrate(x=cemeteries9.5$CRA, errors=cemeteries9.5$Error, normalised=FALSE, calCurves = Curve9.5)
cemeteries9.6=subset(cemeteries, Marine.protein==9.6)
cemeteriescal9.6=calibrate(x=cemeteries9.6$CRA, errors=cemeteries9.6$Error, normalised=FALSE, calCurves = Curve9.6)
cemeteries10.1=subset(cemeteries, Marine.protein==10.1)
cemeteriescal10.1=calibrate(x=cemeteries10.1$CRA, errors=cemeteries10.1$Error, normalised=FALSE, calCurves = Curve10.1)
cemeteries10.3=subset(cemeteries, Marine.protein==10.3)
cemeteriescal10.3=calibrate(x=cemeteries10.3$CRA, errors=cemeteries10.3$Error, normalised=FALSE, calCurves = Curve10.3)
cemeteries10.5=subset(cemeteries, Marine.protein==10.5)
cemeteriescal10.5=calibrate(x=cemeteries10.5$CRA, errors=cemeteries10.5$Error, normalised=FALSE, calCurves = Curve10.5)
cemeteries10.6=subset(cemeteries, Marine.protein==10.6)
cemeteriescal10.6=calibrate(x=cemeteries10.6$CRA, errors=cemeteries10.6$Error, normalised=FALSE, calCurves = Curve10.6)
cemeteries11.8=subset(cemeteries, Marine.protein==11.8)
cemeteriescal11.8=calibrate(x=cemeteries11.8$CRA, errors=cemeteries11.8$Error, normalised=FALSE, calCurves = Curve11.8)
cemeteries12.2=subset(cemeteries, Marine.protein==12.2)
cemeteriescal12.2=calibrate(x=cemeteries12.2$CRA, errors=cemeteries12.2$Error, normalised=FALSE, calCurves = Curve12.2)
cemeteries12.5=subset(cemeteries, Marine.protein==12.5)
cemeteriescal12.5=calibrate(x=cemeteries12.5$CRA, errors=cemeteries12.5$Error, normalised=FALSE, calCurves = Curve12.5)
cemeteries12.9=subset(cemeteries, Marine.protein==12.9)
cemeteriescal12.9=calibrate(x=cemeteries12.9$CRA, errors=cemeteries12.9$Error, normalised=FALSE, calCurves = Curve12.9)
cemeteries14.1=subset(cemeteries, Marine.protein==14.1)
cemeteriescal14.1=calibrate(x=cemeteries14.1$CRA, errors=cemeteries14.1$Error, normalised=FALSE, calCurves = Curve14.1)
cemeteries14.6=subset(cemeteries, Marine.protein==14.6)
cemeteriescal14.6=calibrate(x=cemeteries14.6$CRA, errors=cemeteries14.6$Error, normalised=FALSE, calCurves = Curve14.6)
cemeteries15.3=subset(cemeteries, Marine.protein==15.3)
cemeteriescal15.3=calibrate(x=cemeteries15.3$CRA, errors=cemeteries15.3$Error, normalised=FALSE, calCurves = Curve15.3)
cemeteries16.5=subset(cemeteries, Marine.protein==16.5)
cemeteriescal16.5=calibrate(x=cemeteries16.5$CRA, errors=cemeteries16.5$Error, normalised=FALSE, calCurves = Curve16.5)
cemeteries17.6=subset(cemeteries, Marine.protein==17.6)
cemeteriescal17.6=calibrate(x=cemeteries17.6$CRA, errors=cemeteries17.6$Error, normalised=FALSE, calCurves = Curve17.6)
cemeteries18.8=subset(cemeteries, Marine.protein==18.8)
cemeteriescal18.8=calibrate(x=cemeteries18.8$CRA, errors=cemeteries18.8$Error, normalised=FALSE, calCurves = Curve18.8)
cemeteries18.9=subset(cemeteries, Marine.protein==18.9)
cemeteriescal18.9=calibrate(x=cemeteries18.9$CRA, errors=cemeteries18.9$Error, normalised=FALSE, calCurves = Curve18.9)
cemeteries20=subset(cemeteries, Marine.protein==20)
cemeteriescal20=calibrate(x=cemeteries20$CRA, errors=cemeteries20$Error, normalised=FALSE, calCurves = Curve20)
cemeteries21.1=subset(cemeteries, Marine.protein==21.1)
cemeteriescal21.1=calibrate(x=cemeteries21.1$CRA, errors=cemeteries21.1$Error, normalised=FALSE, calCurves = Curve21.1)
cemeteries21.2=subset(cemeteries, Marine.protein==21.2)
cemeteriescal21.2=calibrate(x=cemeteries21.2$CRA, errors=cemeteries21.2$Error, normalised=FALSE, calCurves = Curve21.2)
cemeteries22.4=subset(cemeteries, Marine.protein==22.4)
cemeteriescal22.4=calibrate(x=cemeteries22.4$CRA, errors=cemeteries22.4$Error, normalised=FALSE, calCurves = Curve22.4)
cemeteries23.5=subset(cemeteries, Marine.protein==23.5)
cemeteriescal23.5=calibrate(x=cemeteries23.5$CRA, errors=cemeteries23.5$Error, normalised=FALSE, calCurves = Curve23.5)
cemeteries24.7=subset(cemeteries, Marine.protein==24.7)
cemeteriescal24.7=calibrate(x=cemeteries24.7$CRA, errors=cemeteries24.7$Error, normalised=FALSE, calCurves = Curve24.7)
cemeteries25.9=subset(cemeteries, Marine.protein==25.9)
cemeteriescal25.9=calibrate(x=cemeteries25.9$CRA, errors=cemeteries25.9$Error, normalised=FALSE, calCurves = Curve25.9)
cemeteries27.1=subset(cemeteries, Marine.protein==27.1)
cemeteriescal27.1=calibrate(x=cemeteries27.1$CRA, errors=cemeteries27.1$Error, normalised=FALSE, calCurves = Curve27.1)
cemeteries28.2=subset(cemeteries, Marine.protein==28.2)
cemeteriescal28.2=calibrate(x=cemeteries28.2$CRA, errors=cemeteries28.2$Error, normalised=FALSE, calCurves = Curve28.2)
cemeteries28.4=subset(cemeteries, Marine.protein==28.4)
cemeteriescal28.4=calibrate(x=cemeteries28.4$CRA, errors=cemeteries28.4$Error, normalised=FALSE, calCurves = Curve28.4)
cemeteries29.4=subset(cemeteries, Marine.protein==29.4)
cemeteriescal29.4=calibrate(x=cemeteries29.4$CRA, errors=cemeteries29.4$Error, normalised=FALSE, calCurves = Curve29.4)
cemeteries30.6=subset(cemeteries, Marine.protein==30.6)
cemeteriescal30.6=calibrate(x=cemeteries30.6$CRA, errors=cemeteries30.6$Error, normalised=FALSE, calCurves = Curve30.6)
cemeteries32=subset(cemeteries, Marine.protein==32)
cemeteriescal32=calibrate(x=cemeteries32$CRA, errors=cemeteries32$Error, normalised=FALSE, calCurves = Curve32)
cemeteries31.8=subset(cemeteries, Marine.protein==31.8)
cemeteriescal31.8=calibrate(x=cemeteries31.8$CRA, errors=cemeteries31.8$Error, normalised=FALSE, calCurves = Curve31.8)
cemeteries32.7=subset(cemeteries, Marine.protein==32.7)
cemeteriescal32.7=calibrate(x=cemeteries32.7$CRA, errors=cemeteries32.7$Error, normalised=FALSE, calCurves = Curve32.7)
cemeteries32.9=subset(cemeteries, Marine.protein==32.9)
cemeteriescal32.9=calibrate(x=cemeteries32.9$CRA, errors=cemeteries32.9$Error, normalised=FALSE, calCurves = Curve32.9)
cemeteries34.1=subset(cemeteries, Marine.protein==34.1)
cemeteriescal34.1=calibrate(x=cemeteries34.1$CRA, errors=cemeteries34.1$Error, normalised=FALSE, calCurves = Curve34.1)
cemeteries37.6=subset(cemeteries, Marine.protein==37.6)
cemeteriescal37.6=calibrate(x=cemeteries37.6$CRA, errors=cemeteries37.6$Error, normalised=FALSE, calCurves = Curve37.6)
cemeteriesterrestrial=subset(cemeteries, is.na(cemeteries$Marine.protein)| Marine.protein=="0")
cemeteriescalterrestrial=calibrate(x=cemeteriesterrestrial$CRA, errors=cemeteriesterrestrial$Error, normalised=FALSE)
combdatecemeteries=combine(cemeteriescal1.2, cemeteriescal2.4, cemeteriescal3.5, cemeteriescal4.7, cemeteriescal5.9, cemeteriescal6.4, cemeteriescal6.7, cemeteriescal7.1, cemeteriescal8.2, cemeteriescal8.3, cemeteriescal8.7, cemeteriescal9.4, cemeteriescal9.5, cemeteriescal9.6, cemeteriescal10.1, cemeteriescal10.3, cemeteriescal10.5, cemeteriescal10.6, cemeteriescal11.8, cemeteriescal12.2, cemeteriescal12.5, cemeteriescal12.9, cemeteriescal14.1, cemeteriescal14.6, cemeteriescal15.3, cemeteriescal16.5, cemeteriescal17.6, cemeteriescal18.8, cemeteriescal18.9, cemeteriescal20, cemeteriescal21.1, cemeteriescal21.2, cemeteriescal22.4, cemeteriescal23.5, cemeteriescal24.7, cemeteriescal25.9, cemeteriescal27.1, cemeteriescal28.2, cemeteriescal28.4, cemeteriescal29.4, cemeteriescal30.6, cemeteriescal31.8,cemeteriescal32, cemeteriescal32.7, cemeteriescal32.9, cemeteriescal34.1, cemeteriescal37.6, cemeteriescalterrestrial, fixIDs = TRUE)
binsc=binPrep(sites=cemeteries$SiteName, ages=cemeteries$CRA, h=50)
samplecem=sampleDates(combdatecemeteries, bins=binsc, nsim=1000, boot=TRUE)
kdecem30=ckde(samplecem, timeRange=c(1650,1050), bw=30)
plot(kdecem30, calendar="BCAD", ylim=c(0,0.002)) + title("Larger cemeteries, n=881")

#all graves without marine correction
callall=calibrate(x=alldates$CRA, errors=alldates$Error, normalised=FALSE)
bins=binPrep(sites=alldates$SiteName, ages=alldates$CRA, h=50)
sampleall=sampleDates(callall, bins=bins, nsim=1000, boot=TRUE)
kdeall30=ckde(sampleall, timeRange = c(1650,1050), bw=30)
plot(kdeall30, calendar="BCAD") +title("All radiocarbon-dated graves without marine correction, n=1086")

#furnished graves without marine correction
furnished=subset(alldates, Furnished.=="Furnished")
calfurnished=calibrate(x=furnished$CRA, errors=furnished$Error, normalised=FALSE)
binsf=binPrep(sites = furnished$SiteName, ages=furnished$CRA, h=50)
samplefurn=sampleDates(calfurnished, bins=binsf, nsim=1000, boot=TRUE)
kdefurnished30=ckde(samplefurn, timeRange=c(1650,1050), bw=30)
plot(kdefurnished30, calendar="BCAD") + title("Furnished graves without marine correction, n=353")

#unfurnished graves without marine correction
unfurnished=subset(alldates, Furnished.=="Unfurnished")
calunfurnished=calibrate(x=unfurnished$CRA, errors=unfurnished$Error, normalised=FALSE)
binsu=binPrep(sites=unfurnished$SiteName, ages=unfurnished$CRA, h=50)
sampleunfurnished=sampleDates(calunfurnished, bins=binsu, nsim=1000, boot=TRUE)
kdeunfurnished30=ckde(sampleunfurnished, timeRange=c(1650,1050), bw=30)
plot(kdeunfurnished30, calendar="BCAD") +title("Unfurnished graves without marine correction, n=715")


#Female furnished graves without marine correction
femfurnished=subset(furnished, Sex=="F" | Sex=="F?")
calfemfurnished=calibrate(x=femfurnished$CRA, errors=femfurnished$Error, normalised=FALSE)
binsff=binPrep(sites = femfurnished$SiteName, ages = femfurnished$CRA, h=50)
samplefemfurn=sampleDates(calfemfurnished, bins=binsff, nsim=1000, boot=TRUE)
kdefemfurnished30=ckde(samplefemfurn, timeRange=c(1650,1050), bw=30)
plot(kdefemfurnished30, calendar="BCAD", ylim=c(0,0.004)) + title("Female furnished graves without marine correction, n=162")


#female unfurnished graves without marine correction
femunfurnished=subset(unfurnished, Sex=="F" | Sex=="F?")
calfemunfurnished=calibrate(x=femunfurnished$CRA, errors=femunfurnished$Error, normalised=FALSE)
binsfu=binPrep(sites = femunfurnished$SiteName, ages = femunfurnished$CRA, h=50)
samplefemunfurnished=sampleDates(calfemunfurnished, bins=binsfu, nsim=1000, boot = TRUE)
kdefemunfurnished30=ckde(samplefemunfurnished, timeRange=c(1650,1050), bw=30)
plot(kdefemunfurnished30, calendar="BCAD", ylim=c(0,0.0025)) + title("Female unfurnished graves without marine correction, n=202")

#male furnished graves without marine correction
malfurnished=subset(furnished, Sex=="M" | Sex=="M?")
calmalfurnished=calibrate(x=malfurnished$CRA, errors=malfurnished$Error, normalised=FALSE)
binsmf=binPrep(sites = malfurnished$SiteName, ages = malfurnished$CRA, h=50)
samplemalfurn=sampleDates(calmalfurnished, bins=binsmf, nsim=1000, boot=TRUE)
kdemalfurnished30=ckde(samplemalfurn, timeRange=c(1650,1050), bw=30)
plot(kdemalfurnished30, calendar="BCAD", ylim=c(0.0045))+ title("Male furnished graves without marine correction, n=156")

#male unfurnished graves without marine correction
malunfurnished=subset(unfurnished, Sex=="M" | Sex=="M?")
calmalunfurnished=calibrate(x=malunfurnished$CRA, errors=malunfurnished$Error, normalised=FALSE)
binsmu=binPrep(sites = malunfurnished$SiteName, ages = malunfurnished$CRA, h=50)
samplemalunfurnished=sampleDates(calmalunfurnished, bins=binsmu, nsim=1000, boot=TRUE)
kdemalunfurnished30=ckde(samplemalunfurnished, timeRange=c(1650,1050), bw=30)
plot(kdemalunfurnished30, calendar="BCAD", ylim=c(0,0.0025))+ title("Male unfurnished graves without marine correction, n=269")

#Furnished burials in the north and west without marine correction
furnishednw=subset(furnished, Region=="NW")
calfnw=calibrate(x=furnishednw$CRA, errors=furnishednw$Error, normalised=FALSE)
binsfnw=binPrep(sites = furnishednw$SiteName, ages=furnishednw$CRA, h=50)
samplefurnnw=sampleDates(calfnw, bins=binsfnw, nsim=1000, boot=TRUE)
kdefurnishednw30=ckde(samplefurnnw, timeRange=c(1650,1050), bw=30)
plot(kdefurnishednw30, calendar="BCAD", ylim=c(0,0.003))+ title("Furnished graves in the north and west without marine correction, n=72")

#Unfurnished burials in the north and west without marine correction
unfurnishednw=subset(unfurnished, Region=="NW")
calunw=calibrate(x=unfurnishednw$CRA, errors=unfurnishednw$Error, normalised=FALSE)
binsunw=binPrep(sites = unfurnishednw$SiteName, ages=unfurnishednw$CRA, h=50)
sampleunfurnnw=sampleDates(calunw, bins=binsunw, nsim=1000, boot=TRUE)
kdeunfurnishednw30=ckde(sampleunfurnnw, timeRange=c(1650,1050), bw=30)
plot(kdeunfurnishednw30, calendar="BCAD") + title ("Unfurnished burials in the north and west without marine correction, n=348")

#Furnished burials in the south and east without marine correction
furnishedse=subset(furnished, Region=="SE")
calfse=calibrate(x=furnishedse$CRA, errors=furnishedse$Error, normalised=FALSE)
binsfse=binPrep(sites = furnishedse$SiteName, ages=furnishedse$CRA, h=50)
samplefurnse=sampleDates(calfse, bins=binsfse, nsim=1000, boot=TRUE)
kdefurnishedse30=ckde(samplefurnse, timeRange=c(1650,1050), bw=30)
plot(kdefurnishedse30, calendar="BCAD", ylim=c(0,0.004)) + title("Furnished graves in the south and east without marine correction, n=281")

#Unfurnished burials in the south and east
unfurnishedse=subset(unfurnished, Region=="SE")
caluse=calibrate(x=unfurnishedse$CRA, errors=unfurnishedse$Error, normalised=FALSE)
binsuse=binPrep(sites = unfurnishedse$SiteName, ages=unfurnishedse$CRA, h=50)
sampleunfurnse=sampleDates(caluse, bins=binsuse, nsim=1000, boot=TRUE)
kdeunfurnishedse30=ckde(sampleunfurnse, timeRange=c(1650,1050), bw=30)
plot(kdeunfurnishedse30, calendar="BCAD", ylim=c(0,0.0025)) + title ("Unfurnished burials in the south and east without marine correction, n=367")

#Richly furnished graves without marine correction
rich=subset(alldates, Rich.materials=="Y")
calrich=calibrate(x=rich$CRA, errors=rich$Error, normalised=FALSE)
binsr=binPrep(sites= rich$SiteName, ages= rich$CRA, h=50)
samplerichmarine=sampleDates(calrich, bins=binsr, nsim=1000, boot=TRUE)
kdemarinerich30=ckde(samplerichmarine, timeRange=c(1650,1050), bw=30)
plot(kdemarinerich30, calendar="BCAD", ylim=c(0,0.0065)) + title("Richly furnished graves without marine correction, n=49")

#Isolated burials without marine correction
isoburial=subset(alldates, Burial.Type=="Isolated burial")
calisoburial=calibrate(x=isoburial$CRA, errors=isoburial$Error, normalised=FALSE)
sampleisoburial=sampleDates(calisoburial, nsim=1000, boot=TRUE)
kdeisoburial=ckde(sampleisoburial, timeRange = c(1650, 1050), bw=30)
plot(kdeisoburial, calendar = "BCAD", ylim=c(0,0.003))+ title("Isolated graves without marine correction, n=72")

#Small groups of burials without marine correction
smallgroup=subset(alldates, Burial.Type=="Small group of burials")
calsmallgroup=calibrate(x=smallgroup$CRA, errors=smallgroup$Error, normalised=FALSE)
binssgc=binPrep(sites = smallgroup$SiteName, ages = smallgroup$CRA, h=50)
samplesmallgroup=sampleDates(calsmallgroup, bins=binssgc, nsim=1000, boot=TRUE)
kdesmallgroup=ckde(samplesmallgroup, timeRange = c(1650, 1050), bw=30)
plot(kdesmallgroup, calendar="BCAD", ylim=c(0,0.003))+ title("Small groups of burials without marine correction, n=81")

#large cemeteries without marine correction
cemeteries=subset(alldates, Burial.Type== "Unfurnished cemetery" | Burial.Type=="Furnished cemetery")
calcem=calibrate(x=cemeteries$CRA, errors=cemeteries$Error, normalised=FALSE)
binsc=binPrep(sites=cemeteries$SiteName, ages = cemeteries$CRA, h=50)
samplecem=sampleDates(calcem, bins=binsc, nsim=1000, boot=TRUE)
kdecem=ckde(samplecem, timeRange = c(1650,1050), bw=30)
plot(kdecem, calendar = "BCAD") + title ("Larger cemeteries without marine correction, n=881")

#Execution cemeteries
execution=subset(alldates, Burial.Type=="Execution cemetery")
execution4.7=subset(execution, Marine.protein==4.7)
executioncal4.7=calibrate(x=execution4.7$CRA, errors=execution4.7$Error, normalised=FALSE, calCurves = Curve4.7)
execution9.4=subset(execution, Marine.protein==9.4)
executioncal9.4=calibrate(x=execution9.4$CRA, errors=execution9.4$Error, normalised=FALSE, calCurves = Curve9.4)
execution11.8=subset(execution, Marine.protein==11.8)
executioncal11.8=calibrate(x=execution11.8$CRA, errors=execution11.8$Error, normalised=FALSE, calCurves = Curve11.8)
execution12.9=subset(execution, Marine.protein==12.9)
executioncal12.9=calibrate(x=execution12.9$CRA, errors=execution12.9$Error, normalised=FALSE, calCurves = Curve12.9)
execution14.1=subset(execution, Marine.protein==14.1)
executioncal14.1=calibrate(x=execution14.1$CRA, errors=execution14.1$Error, normalised=FALSE, calCurves = Curve14.1)
execution15.3=subset(execution, Marine.protein==15.3)
executioncal15.3=calibrate(x=execution15.3$CRA, errors=execution15.3$Error, normalised=FALSE, calCurves = Curve15.3)
execution16.5=subset(execution, Marine.protein==16.5)
executioncal16.5=calibrate(x=execution16.5$CRA, errors=execution16.5$Error, normalised=FALSE, calCurves = Curve16.5)
execution17.6=subset(execution, Marine.protein==17.6)
executioncal17.6=calibrate(x=execution17.6$CRA, errors=execution17.6$Error, normalised=FALSE, calCurves = Curve17.6)
execution18.8=subset(execution, Marine.protein==18.8)
executioncal18.8=calibrate(x=execution18.8$CRA, errors=execution18.8$Error, normalised=FALSE, calCurves = Curve18.8)
terrestrialexecution=subset(execution, is.na(execution$Marine.protein)| Marine.protein=="0")
executioncalterrestial=calibrate(x=terrestrialexecution$CRA, errors=terrestrialexecution$Error, normalised=FALSE)
combdateexecution=combine(executioncal4.7, executioncal9.4, executioncal11.8, executioncal12.9, executioncal14.1, executioncal15.3, executioncal16.5, executioncal17.6, executioncal18.8, executioncalterrestial, fixIDs = TRUE)
binsexecution= binPrep(sites = execution$SiteName, ages=execution$CRA, h=50)
sampleexecutionmarine=sampleDates(combdateexecution, bins=binsexecution, nsim=1000)
kdeexecution30marine=ckde(sampleexecutionmarine, timeRange=c(1650,1050), bw=30)
plot(kdeexecution30marine, calendar="BCAD", ylim=c(0,0.005)) + title("Execution cemeteries, n=39")

#KDE for furnished burials, without the dates from the ASCF
furnishedminusascf=subset(furnished, Notes!="From ASCF")
furnishedminusascf1.2=subset(furnishedminusascf, Marine.protein==1.2)
furnishedminusascfcal1.2=calibrate(x=furnishedminusascf1.2$CRA, errors=furnishedminusascf1.2$Error, normalised=FALSE, calCurves=Curve1.2)
furnishedminusascf2.4=subset(furnishedminusascf, Marine.protein==2.4)
furnishedminusascfcal2.4=calibrate(x=furnishedminusascf2.4$CRA, errors=furnishedminusascf2.4$Error, normalised=FALSE, calCurves = Curve2.4)
furnishedminusascf3.5=subset(furnishedminusascf, Marine.protein==3.5)
furnishedminusascfcal3.5=calibrate(x=furnishedminusascf3.5$CRA, errors=furnishedminusascf3.5$Error, normalised=FALSE, calCurves = Curve3.5)
furnishedminusascf4.7=subset(furnishedminusascf, Marine.protein==4.7)
furnishedminusascfcal4.7=calibrate(x=furnishedminusascf4.7$CRA, errors=furnishedminusascf4.7$Error, normalised=FALSE, calCurves = Curve4.7)
furnishedminusascf5.9=subset(furnishedminusascf, Marine.protein==5.9)
furnishedminusascfcal5.9=calibrate(x=furnishedminusascf5.9$CRA, errors=furnishedminusascf5.9$Error, normalised=FALSE, calCurves = Curve5.9)
furnishedminusascf6.7=subset(furnishedminusascf, Marine.protein==6.7)
furnishedminusascfcal6.7=calibrate(x=furnishedminusascf6.7$CRA, errors=furnishedminusascf6.7$Error, normalised=FALSE, calCurves = Curve6.7)
furnishedminusascf7.1=subset(furnishedminusascf, Marine.protein==7.1)
furnishedminusascfcal7.1=calibrate(x=furnishedminusascf7.1$CRA, errors=furnishedminusascf7.1$Error, normalised=FALSE, calCurves = Curve7.1)
furnishedminusascf8=subset(furnishedminusascf, Marine.protein==8)
furnishedminusascfcal8=calibrate(x=furnishedminusascf8$CRA, errors=furnishedminusascf8$Error, normalised=FALSE, calCurves = Curve8)
furnishedminusascf8.2=subset(furnishedminusascf, Marine.protein==8.2)
furnishedminusascfcal8.2=calibrate(x=furnishedminusascf8.2$CRA, errors=furnishedminusascf8.2$Error, normalised=FALSE, calCurves = Curve8.2)
furnishedminusascf8.3=subset(furnishedminusascf, Marine.protein==8.3)
furnishedminusascfcal8.3=calibrate(x=furnishedminusascf8.3$CRA, errors=furnishedminusascf8.3$Error, normalised=FALSE, calCurves = Curve8.3)
furnishedminusascf8.7=subset(furnishedminusascf, Marine.protein==8.7)
furnishedminusascfcal8.7=calibrate(x=furnishedminusascf8.7$CRA, errors=furnishedminusascf8.7$Error, normalised=FALSE, calCurves = Curve8.7)
furnishedminusascf9.4=subset(furnishedminusascf, Marine.protein==9.4)
furnishedminusascfcal9.4=calibrate(x=furnishedminusascf9.4$CRA, errors=furnishedminusascf9.4$Error, normalised=FALSE, calCurves = Curve9.4)
furnishedminusascf9.5=subset(furnishedminusascf, Marine.protein==9.5)
furnishedminusascfcal9.5=calibrate(x=furnishedminusascf9.5$CRA, errors=furnishedminusascf9.5$Error, normalised=FALSE, calCurves = Curve9.5)
furnishedminusascf9.6=subset(furnishedminusascf, Marine.protein==9.6)
furnishedminusascfcal9.6=calibrate(x=furnishedminusascf9.6$CRA, errors=furnishedminusascf9.6$Error, normalised=FALSE, calCurves = Curve9.6)
furnishedminusascf10.1=subset(furnishedminusascf, Marine.protein==10.1)
furnishedminusascfcal10.1=calibrate(x=furnishedminusascf10.1$CRA, errors=furnishedminusascf10.1$Error, normalised=FALSE, calCurves = Curve10.1)
furnishedminusascf10.6=subset(furnishedminusascf, Marine.protein==10.6)
furnishedminusascfcal10.6=calibrate(x=furnishedminusascf10.6$CRA, errors=furnishedminusascf10.6$Error, normalised=FALSE, calCurves = Curve10.6)
furnishedminusascf11.8=subset(furnishedminusascf, Marine.protein==11.8)
furnishedminusascfcal11.8=calibrate(x=furnishedminusascf11.8$CRA, errors=furnishedminusascf11.8$Error, normalised=FALSE, calCurves = Curve11.8)
furnishedminusascf12.9=subset(furnishedminusascf, Marine.protein==12.9)
furnishedminusascfcal12.9=calibrate(x=furnishedminusascf12.9$CRA, errors=furnishedminusascf12.9$Error, normalised=FALSE, calCurves = Curve12.9)
furnishedminusascf13.9=subset(furnishedminusascf, Marine.protein==13.9)
furnishedminusascfcal13.9=calibrate(x=furnishedminusascf13.9$CRA, errors=furnishedminusascf13.9$Error, normalised=FALSE, calCurves = Curve13.9)
furnishedminusascf14.1=subset(furnishedminusascf, Marine.protein==14.1)
furnishedminusascfcal14.1=calibrate(x=furnishedminusascf14.1$CRA, errors=furnishedminusascf14.1$Error, normalised=FALSE, calCurves = Curve14.1)
furnishedminusascf15.3=subset(furnishedminusascf, Marine.protein==15.3)
furnishedminusascfcal15.3=calibrate(x=furnishedminusascf15.3$CRA, errors=furnishedminusascf15.3$Error, normalised=FALSE, calCurves = Curve15.3)
furnishedminusascf16.5=subset(furnishedminusascf, Marine.protein==16.5)
furnishedminusascfcal16.5=calibrate(x=furnishedminusascf16.5$CRA, errors=furnishedminusascf16.5$Error, normalised=FALSE, calCurves = Curve16.5)
furnishedminusascf17.6=subset(furnishedminusascf, Marine.protein==17.6)
furnishedminusascfcal17.6=calibrate(x=furnishedminusascf17.6$CRA, errors=furnishedminusascf17.6$Error, normalised=FALSE, calCurves = Curve17.6)
furnishedminusascf18.8=subset(furnishedminusascf, Marine.protein==18.8)
furnishedminusascfcal18.8=calibrate(x=furnishedminusascf18.8$CRA, errors=furnishedminusascf18.8$Error, normalised=FALSE, calCurves = Curve18.8)
furnishedminusascf20=subset(furnishedminusascf, Marine.protein==20)
furnishedminusascfcal20=calibrate(x=furnishedminusascf20$CRA, errors=furnishedminusascf20$Error, normalised=FALSE, calCurves = Curve20)
furnishedminusascf21.2=subset(furnishedminusascf, Marine.protein==21.2)
furnishedminusascfcal21.2=calibrate(x=furnishedminusascf21.2$CRA, errors=furnishedminusascf21.2$Error, normalised=FALSE, calCurves = Curve21.2)
furnishedminusascf22.4=subset(furnishedminusascf, Marine.protein==22.4)
furnishedminusascfcal22.4=calibrate(x=furnishedminusascf22.4$CRA, errors=furnishedminusascf22.4$Error, normalised=FALSE, calCurves = Curve22.4)
furnishedminusascf23.5=subset(furnishedminusascf, Marine.protein==23.5)
furnishedminusascfcal23.5=calibrate(x=furnishedminusascf23.5$CRA, errors=furnishedminusascf23.5$Error, normalised=FALSE, calCurves = Curve23.5)
furnishedminusascf25.9=subset(furnishedminusascf, Marine.protein==25.9)
furnishedminusascfcal25.9=calibrate(x=furnishedminusascf25.9$CRA, errors=furnishedminusascf25.9$Error, normalised=FALSE, calCurves = Curve25.9)
furnishedminusascf27.1=subset(furnishedminusascf, Marine.protein==27.1)
furnishedminusascfcal27.1=calibrate(x=furnishedminusascf27.1$CRA, errors=furnishedminusascf27.1$Error, normalised=FALSE, calCurves = Curve27.1)
furnishedminusascf28.2=subset(furnishedminusascf, Marine.protein==28.2)
furnishedminusascfcal28.2=calibrate(x=furnishedminusascf28.2$CRA, errors=furnishedminusascf28.2$Error, normalised=FALSE, calCurves = Curve28.2)
furnishedminusascf28.4=subset(furnishedminusascf, Marine.protein==28.4)
furnishedminusascfcal28.4=calibrate(x=furnishedminusascf28.4$CRA, errors=furnishedminusascf28.4$Error, normalised=FALSE, calCurves = Curve28.4)
furnishedminusascf32.7=subset(furnishedminusascf, Marine.protein==32.7)
furnishedminusascfcal32.7=calibrate(x=furnishedminusascf32.7$CRA, errors=furnishedminusascf32.7$Error, normalised=FALSE, calCurves = Curve32.7)
furnishedminusascf32.9=subset(furnishedminusascf, Marine.protein==32.9)
furnishedminusascfcal32.9=calibrate(x=furnishedminusascf32.9$CRA, errors=furnishedminusascf32.9$Error, normalised=FALSE, calCurves = Curve32.9)
furnishedminusascf37.6=subset(furnishedminusascf, Marine.protein==37.6)
furnishedminusascfcal37.6=calibrate(x=furnishedminusascf37.6$CRA, errors=furnishedminusascf37.6$Error, normalised=FALSE, calCurves = Curve37.6)
furnminusascfterrestrial=subset(furnishedminusascf, is.na(furnishedminusascf$Marine.protein)| Marine.protein=="0")
furnishedminusascfcalterrestrial=calibrate(x=furnminusascfterrestrial$CRA, errors=furnminusascfterrestrial$Error, normalised=FALSE)
combfurnishedminusascfdate=combine(furnishedminusascfcal1.2, furnishedminusascfcal2.4, furnishedminusascfcal3.5, furnishedminusascfcal4.7, furnishedminusascfcal5.9, furnishedminusascfcal6.7, furnishedminusascfcal7.1, furnishedminusascfcal8, furnishedminusascfcal8.2, furnishedminusascfcal8.3, furnishedminusascfcal8.7, furnishedminusascfcal9.4, furnishedminusascfcal9.5, furnishedminusascfcal9.6, furnishedminusascfcal10.1, furnishedminusascfcal10.6, furnishedminusascfcal11.8, furnishedminusascfcal12.9, furnishedminusascfcal13.9, furnishedminusascfcal14.1, furnishedminusascfcal15.3, furnishedminusascfcal16.5, furnishedminusascfcal17.6, furnishedminusascfcal18.8, furnishedminusascfcal20, furnishedminusascfcal21.2, furnishedminusascfcal22.4, furnishedminusascfcal23.5, furnishedminusascfcal25.9, furnishedminusascfcal27.1, furnishedminusascfcal28.2, furnishedminusascfcal28.4, furnishedminusascfcal32.7, furnishedminusascfcal32.9, furnishedminusascfcal37.6, furnishedminusascfcalterrestrial, fixIDs = TRUE)
binsfmascf=binPrep(sites = furnishedminusascf$SiteName, ages=furnishedminusascf$CRA, h=50)
samplefurnminusascfmarine=sampleDates(combfurnishedminusascfdate, bins=binsfmascf, nsim=1000, boot=TRUE)
kdemarinefurnishedminusascf30=ckde(samplefurnminusascfmarine, timeRange=c(1650,1050), bw=30)
plot(kdemarinefurnishedminusascf30, calendar="BCAD", ylim=c(0,0.0035)) + title("Furnished graves without those dates in the ASCF, n=269")