##################################################################################### #two packages that need to be installed require(vegan) require(rsq) #function for confidence intervales bootci=function(temp) { x.resample<-numeric(); xmean.resample<-numeric() for (i in 1:100) { for (j in 1:length(temp)) {x.resample[j]<-sample(temp,1,replace=T)} xmean.resample[i]<-mean(x.resample, na.rm=T) } out=list(resample=xmean.resample) out } ##################################################################################### #Tomašových, Adam; Gallmetzer, Ivo; Haselmair, Alexandra; Zuschin, Martin; 2021: #Inferring time averaging and hiatus durations in the stratigraphic record of high-frequency depositional sequences. Sedimentology ##################################################################################### #INDIVIDUAL-BASED SCORED DATASET - BASED ON AGE-DATED SHELLS ONLY - FIGURES 10 and 11 #ANALYSES BASED ON ALL SHELLS - SEE SCRIPT BELOW ##################################################################################### setwd("D:/Data/Adriatic_Sea_project/MS-Brijuni/Revision II/Tables final") taphonomic.data <- read.delim("Table S1.txt", header=TRUE) ages.for.clock=taphonomic.data[,"Age"] specimen=taphonomic.data$"specimen_no" taxon=taphonomic.data$taxon periostracum=taphonomic.data$"Periostracum" periostracum=ifelse(periostracum>0.5, 1,0) ext.conch=taphonomic.data$"Conchiolin.external" int.conch=taphonomic.data$"Conchiolin.internal" conchiolin=ifelse(ext.conch==0 | int.conch==0, 0, 1) periostracum=ext.conch pyrite=taphonomic.data$"pyrite" discoloration=taphonomic.data$"External.discoloration" ornamentation=taphonomic.data$"External.ornamentation" ext.encrustation=taphonomic.data$"External.encrustation" ext.worn.encrustation=taphonomic.data$"Worn.ext.encrustation" ext.unworn.encrustation=taphonomic.data$"Unworn.ext.encrustation" int.worn.encrustation=taphonomic.data$"Worn.int.encrustation" int.unworn.encrustation=taphonomic.data$"Unworn.int.encrustation" int.encrustation=taphonomic.data$"Internal.encrustation" core=taphonomic.data$"Station" completeness=taphonomic.data$"Valve.completeness"/100 fragmentation=1-completeness stained=taphonomic.data$"Staining" bioerosion=taphonomic.data$"Internal.microbioerosion" dissolution=taphonomic.data$"Internal.dissolution" cementation=taphonomic.data$"Internal.cementation" depths=as.numeric(as.vector(taphonomic.data$"max..depth..cm.")) d10cm=as.numeric(as.vector(taphonomic.data$"X10.cm.group..cm.")) d5cm=as.numeric(as.vector(taphonomic.data$"X5.cm.group..cm.")) dfinecm=as.numeric(as.vector(taphonomic.data$"max..depth..cm.")) ages=ages.for.clock length=taphonomic.data$"length..mm." height=taphonomic.data$"height..mm." Timoclea.length=length[taxon=="Timoclea ovata"] Timoclea.height=height[taxon=="Timoclea ovata"] Timoclea.discoloration=discoloration[taxon=="Timoclea ovata"] Timoclea.ornamentation=ornamentation[taxon=="Timoclea ovata"] Timoclea.ext.encrustation=ext.encrustation[taxon=="Timoclea ovata"] Timoclea.ext.worn.encrustation=ext.worn.encrustation[taxon=="Timoclea ovata"] Timoclea.ext.unworn.encrustation=ext.unworn.encrustation[taxon=="Timoclea ovata"] Timoclea.int.encrustation=int.encrustation[taxon=="Timoclea ovata"] Timoclea.int.worn.encrustation=int.worn.encrustation[taxon=="Timoclea ovata"] Timoclea.int.unworn.encrustation=int.unworn.encrustation[taxon=="Timoclea ovata"] temp1=Timoclea.ext.unworn.encrustation; temp1=ifelse(Timoclea.ext.unworn.encrustation==0 & Timoclea.int.unworn.encrustation==0,0,temp1); temp1=ifelse(Timoclea.ext.unworn.encrustation==1 | Timoclea.int.unworn.encrustation==1,1,temp1) Timoclea.unworn.encrustation=temp1 temp1=Timoclea.ext.worn.encrustation; temp1=ifelse(Timoclea.ext.worn.encrustation==0 & Timoclea.int.worn.encrustation==0,0,temp1); temp1=ifelse(Timoclea.ext.worn.encrustation==1 | Timoclea.int.worn.encrustation==1,1,temp1) Timoclea.worn.encrustation=temp1 Timoclea.cementation=cementation[taxon=="Timoclea ovata"] Timoclea.int.worn.encrustation=int.worn.encrustation[taxon=="Timoclea ovata"] Timoclea.int.unworn.encrustation=int.unworn.encrustation[taxon=="Timoclea ovata"] Timoclea.int.encrustation=int.encrustation[taxon=="Timoclea ovata"] Timoclea.core=core[taxon=="Timoclea ovata"] Timoclea.fragmentation=fragmentation[taxon=="Timoclea ovata"] Timoclea.stained=stained[taxon=="Timoclea ovata"] Timoclea.bioerosion=bioerosion[taxon=="Timoclea ovata"] Timoclea.dissolution=dissolution[taxon=="Timoclea ovata"] Timoclea.depths=depths[taxon=="Timoclea ovata"] Timoclea.10cm=d10cm[taxon=="Timoclea ovata"] Timoclea.5cm=d5cm[taxon=="Timoclea ovata"] Timoclea.ages.for.clock=ages.for.clock[taxon=="Timoclea ovata"] Timoclea.specimen=specimen[taxon=="Timoclea ovata"] #EIGHT ALTERATION VARIABLES, COMBINE INTERNAL AND EXTERNAL ENCRUSTATION Timoclea.scores=cbind(Timoclea.discoloration,Timoclea.ornamentation,Timoclea.unworn.encrustation,Timoclea.worn.encrustation, Timoclea.stained,Timoclea.bioerosion,Timoclea.dissolution,Timoclea.cementation) labels=c("discoloration","ornament loss","fresh encrusters","worn encrusters","stained","bioerosion","dissolution","cement") Timoclea.dated.depth.10cm=Timoclea.10cm[!is.na(Timoclea.ages.for.clock)] Timoclea.dated.core=Timoclea.core[!is.na(Timoclea.ages.for.clock)] Timoclea.dated.age=Timoclea.ages.for.clock[!is.na(Timoclea.ages.for.clock)] Timoclea.dated.specimen=Timoclea.specimen[!is.na(Timoclea.ages.for.clock)] Timoclea.ind.scores=Timoclea.scores[!is.na(Timoclea.ages.for.clock),] Timoclea.ind.ages=ages.for.clock[!is.na(Timoclea.ages.for.clock)] colnames(Timoclea.ind.scores)=labels ########################################################################################## #TIMOCLEA - INDIVIDUAL SCORES BY CORE ########################################################################################## Timoclea.M44.pyrite=pyrite[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.cementation=cementation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.discoloration=discoloration[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.ornamentation=ornamentation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.ext.encrustation=ext.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.ext.worn.encrustation=ext.worn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.ext.unworn.encrustation=ext.unworn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.int.worn.encrustation=int.worn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.int.unworn.encrustation=int.unworn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.int.encrustation=int.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.core=core[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.fragmentation=fragmentation[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.stained=stained[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.bioerosion=bioerosion[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.dissolution=dissolution[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.depths=depths[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.10cm=d10cm[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.5cm=d5cm[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.ages.for.clock=ages.for.clock[taxon=="Timoclea ovata" & core=="Brijuni 44"] Timoclea.M44.total.encrustation=ifelse(int.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"] ==1 | ext.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 44"]==1,1,0) Timoclea.M44.worn.encrustation=Timoclea.worn.encrustation[Timoclea.core=="Brijuni 44"] Timoclea.M44.unworn.encrustation=Timoclea.unworn.encrustation[Timoclea.core=="Brijuni 44"] Timoclea.M40.pyrite=pyrite[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.cementation=cementation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.discoloration=discoloration[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.ornamentation=ornamentation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.ext.encrustation=ext.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.ext.worn.encrustation=ext.worn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.ext.unworn.encrustation=ext.unworn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.int.worn.encrustation=int.worn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.int.unworn.encrustation=int.unworn.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.int.encrustation=int.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.core=core[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.fragmentation=fragmentation[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.stained=stained[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.bioerosion=bioerosion[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.dissolution=dissolution[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.depths=depths[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.10cm=d10cm[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.5cm=d5cm[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.ages.for.clock=ages.for.clock[taxon=="Timoclea ovata" & core=="Brijuni 40"] Timoclea.M40.total.encrustation=ifelse(int.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"] ==1 | ext.encrustation[taxon=="Timoclea ovata" & core=="Brijuni 40"]==1,1,0) Timoclea.M40.worn.encrustation=Timoclea.worn.encrustation[Timoclea.core=="Brijuni 40"] Timoclea.M40.unworn.encrustation=Timoclea.unworn.encrustation[Timoclea.core=="Brijuni 40"] ########################################################################################################### #MULTIVARIATE - INDIRECT ORDINATIONS - PRINCIPAL COORDINATE ANALYSES #DIRECT ORDINATONS - REDUNDANCY ANALYSES AND CONSTRAINED ANALYSES OF PRINCIPAL COORDINATES ########################################################################################################### #M44 CORE - STRATIGRAPHIC UNITS IN ORDINATIONS AND EFFECT OF POSTMORTEM AGE ON TIMOCLEA TAPHONOMY ########################################################################################################### Timoclea.ind.rda.0.20cm.M44=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm<21 & Timoclea.dated.core=="Brijuni 44",]~log(Timoclea.dated.age[Timoclea.dated.depth.10cm<21 & Timoclea.dated.core=="Brijuni 44"])) Timoclea.ind.rda.20.90cm.M44=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm>20 & Timoclea.dated.depth.10cm<90 & Timoclea.dated.core=="Brijuni 44",]~log(Timoclea.dated.age[Timoclea.dated.depth.10cm>20 & Timoclea.dated.depth.10cm<90 & Timoclea.dated.core=="Brijuni 44"])) Timoclea.ind.rda.90.120cm.M44=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm>89 & Timoclea.dated.depth.10cm<121 & Timoclea.dated.core=="Brijuni 44",]~log(Timoclea.dated.age[Timoclea.dated.depth.10cm>89 & Timoclea.dated.depth.10cm<121 & Timoclea.dated.core=="Brijuni 44"])) Timoclea.ind.rda.120.160cm.M44=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm>120 & Timoclea.dated.core=="Brijuni 44",]~log(Timoclea.dated.age[Timoclea.dated.depth.10cm>120 & Timoclea.dated.core=="Brijuni 44"])) anova(Timoclea.ind.rda.0.20cm.M44) anova(Timoclea.ind.rda.20.90cm.M44) anova(Timoclea.ind.rda.90.120cm.M44) anova(Timoclea.ind.rda.120.160cm.M44) ############################################################################################## #M40 CORE - STRATIGRAPHIC UNITS IN ORDINATIONS AND EFFECT OF POSTMORTEM AGE ON TIMOCLEA TAPHONOMY ############################################################################################## Timoclea.ind.rda.0.20cm.M40=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm<21 & Timoclea.dated.core=="Brijuni 40",]~Timoclea.dated.age[Timoclea.dated.depth.10cm<21 & Timoclea.dated.core=="Brijuni 40"]) Timoclea.ind.rda.20.40cm.M40=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm>20 & Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40",]~Timoclea.dated.age[Timoclea.dated.depth.10cm>20 & Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40"]) Timoclea.ind.rda.40.80cm.M40=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm>40 & Timoclea.dated.depth.10cm<81 & Timoclea.dated.core=="Brijuni 40",]~Timoclea.dated.age[Timoclea.dated.depth.10cm>40 & Timoclea.dated.depth.10cm<81 & Timoclea.dated.core=="Brijuni 40"]) Timoclea.ind.rda.80.120cm.M40=rda(Timoclea.ind.scores[Timoclea.dated.depth.10cm>81 & Timoclea.dated.core=="Brijuni 40",]~Timoclea.dated.age[Timoclea.dated.depth.10cm>81 & Timoclea.dated.core=="Brijuni 40"]) anova(Timoclea.ind.rda.0.20cm.M40) anova(Timoclea.ind.rda.20.40cm.M40) anova(Timoclea.ind.rda.40.80cm.M40) anova(Timoclea.ind.rda.80.120cm.M40) ###################################################################################################################################################### #WHOLE-CORE CONSTRAINED ANALYSIS OF PRINCIPAL COORDINATES, PARTIALLING OUT THE DEPTH EFFECT ON TIMOCLEA TAPHONOMY ###################################################################################################################################################### Timoclea.ind.cap.all.M44.part=capscale(Timoclea.ind.scores[Timoclea.dated.core=="Brijuni 44",]~log(Timoclea.dated.age[Timoclea.dated.core=="Brijuni 44"]) +Condition(Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"]), distance = "manhattan") anova(Timoclea.ind.cap.all.M44.part) Timoclea.ind.cap.all.M40.part=capscale(Timoclea.ind.scores[Timoclea.dated.core=="Brijuni 40",]~log(Timoclea.dated.age[Timoclea.dated.core=="Brijuni 40"]) +Condition(Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 40"]), distance = "manhattan") anova(Timoclea.ind.cap.all.M40.part) Timoclea.ind.cap.0.40cm.M44=capscale(Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44",]~Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44"], distance="manhattan") Timoclea.ind.cap.0.40cm.M40=capscale(Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40",]~Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40"], distance="manhattan") Timoclea.ind.pco.0.40cm.M44=cmdscale(vegdist(Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44",], method="manhattan"), eig=T) Timoclea.ind.pco.0.40cm.M40=cmdscale(vegdist(Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40",], method="manhattan"), eig=T) VARIABLE="stained" par(mfrow=c(3,2)) par(mar=c(4,4,2,1)) plot(scores(Timoclea.ind.cap.all.M44.part)$sites[,1:2], xlab="CAP axis 1", ylab="PCO axis 1", frame=F, pch=16, type="p", main="Whole core-M44", cex.main=0.9) points(scores(Timoclea.ind.cap.all.M44.part)$sites[Timoclea.ind.scores[Timoclea.dated.core=="Brijuni 44",VARIABLE]==0,1:2],pch=21, bg="white") points(scores(Timoclea.ind.cap.all.M44.part)$sites[Timoclea.ind.scores[Timoclea.dated.core=="Brijuni 44",VARIABLE]==1,1:2],pch=21, bg="black") ordisurf(x=scores(Timoclea.ind.cap.all.M44.part)$sites, y=Timoclea.dated.age[Timoclea.dated.core=="Brijuni 44"], family = "quasipoisson", lty=2, add=T, col="gray31", levels=seq(0,10000,by=500), labcex=0.8, knots=5) plot(scores(Timoclea.ind.cap.all.M40.part)$sites[,1:2], xlab="CAP axis 1", ylab="PCO axis 1", frame=F, pch=16, type="p", main="Whole core-M40", cex.main=0.9) points(scores(Timoclea.ind.cap.all.M40.part)$sites[Timoclea.ind.scores[Timoclea.dated.core=="Brijuni 40",VARIABLE]==0,1:2],pch=21, bg="white") points(scores(Timoclea.ind.cap.all.M40.part)$sites[Timoclea.ind.scores[Timoclea.dated.core=="Brijuni 40",VARIABLE]==1,1:2],pch=21, bg="black") ordisurf(x=scores(Timoclea.ind.cap.all.M40.part)$sites, y=Timoclea.dated.age[Timoclea.dated.core=="Brijuni 40"], family = "quasipoisson", lty=2, add=T, col="gray31", levels=seq(0,10000,by=500), labcex=0.8, knots=5) ###################################################################################################################### #TOP 40 CM ONLY - CONSTRAINED ANALYSIS OF PRINCIPAL COORDINATES, REMOVING SUBSURFACE INCREMENTS TO AVOID THE DEPTH EFFECT ###################################################################################################################### VARIABLE="stained" plot(scores(Timoclea.ind.cap.0.40cm.M44)$sites[,1:2], xlab="CAP axis 1", ylab="PCO axis 1", frame=F, pch=16, type="p", main="M44-top 40 cm", cex.main=0.9) points(scores(Timoclea.ind.cap.0.40cm.M44)$sites[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44",VARIABLE]==0,1:2],pch=21, bg="white") points(scores(Timoclea.ind.cap.0.40cm.M44)$sites[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44",VARIABLE]==1,1:2],pch=21, bg="black") ordisurf(x=scores(Timoclea.ind.cap.0.40cm.M44)$sites, y=Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44"], family = "quasipoisson", lty=2, add=T, col="gray31", levels=seq(0,10000,by=500), labcex=0.8, knots=5) plot(scores(Timoclea.ind.cap.0.40cm.M40)$sites[,1:2],xlab="CAP axis 1", ylab="PCO axis 1", frame=F, pch=16, type="p", main="M40-top 40 cm", cex.main=0.9) points(scores(Timoclea.ind.cap.0.40cm.M40)$sites[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40",VARIABLE]==0,1:2],pch=21, bg="white") points(scores(Timoclea.ind.cap.0.40cm.M40)$sites[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40",VARIABLE]==1,1:2],pch=21, bg="black") ordisurf(x=scores(Timoclea.ind.cap.0.40cm.M40)$sites, y=Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40"], family = "quasipoisson", lty=2, add=T, col="gray31", levels=seq(0,10000,by=500), labcex=0.8, knots=5) ########################################################################### #TOP 40 CM - PRINCIPAL COORDINATE ANALYSIS ########################################################################### plot(Timoclea.ind.pco.0.40cm.M44$points[,1:2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, type="p", main="Top 40 cm-M44", cex.main=0.9) points(Timoclea.ind.pco.0.40cm.M44$points[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44",VARIABLE]==0,1:2],pch=21, bg="white", cex=1.4) points(Timoclea.ind.pco.0.40cm.M44$points[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44",VARIABLE]==1,1:2],pch=21, bg="black", cex=1.4) ordisurf(x=Timoclea.ind.pco.0.40cm.M44$points, y=Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44"], family = "quasipoisson", lty=2, add=T, col="gray31", levels=seq(0,10000,by=500), labcex=0.8, knots=5) plot(Timoclea.ind.pco.0.40cm.M40$points[,1:2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, type="p", main="Top 40 cm-M40", cex.main=0.9) points(Timoclea.ind.pco.0.40cm.M40$points[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40",VARIABLE]==0,1:2],pch=21, bg="white", cex=1.4) points(Timoclea.ind.pco.0.40cm.M40$points[Timoclea.ind.scores[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40",VARIABLE]==1,1:2],pch=21, bg="black", cex=1.4) ordisurf(x=Timoclea.ind.pco.0.40cm.M40$points, y=Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40"], family = "quasipoisson", lty=2, add=T, col="gray31", levels=seq(0,10000,by=500), labcex=0.8, knots=5) legend(x="bottomleft", legend=c("Non-stained","Stained"), pt.bg=c("white","black"), pch=21, bty="n") ###################################################################################################################################################### #AGGREGATE SHELLS INTO 500 YEAR AGE COHORTS FOR TOP 40 CM - AGGREGATE-LEVEL MULTIVARIATE TAPHONOMIC CLOCK ###################################################################################################################################################### temp.M44=cut(Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44"], breaks=seq(0,12000,by=500), labels=F)*500 temp.M40=cut(Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40"], breaks=seq(0,12000,by=500), labels=F)*500 par(mfrow=c(3,2)) par(mar=c(4,4,2,1)) plot(Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 44"],scores(Timoclea.ind.cap.0.40cm.M44)$sites[,1], main="Not aggregated-M44", xlab="Age (y)", ylab="Preservation axis CAP 1",pch=16, frame=F) plot(Timoclea.dated.age[Timoclea.dated.depth.10cm<41 & Timoclea.dated.core=="Brijuni 40"],scores(Timoclea.ind.cap.0.40cm.M40)$sites[,1], main="Not aggregated-M40", , xlab="Age (y)", ylab="Preservation axis CAP 1",pch=16, frame=F) mean.M44.CAP1.500.y=tapply(scores(Timoclea.ind.cap.0.40cm.M44)$sites[,1],temp.M44, mean) boot.M44.CAP1.500.y=tapply(scores(Timoclea.ind.cap.0.40cm.M44)$sites[,1],temp.M44, bootci) LCI.M44.CAP1.500.y=unlist(lapply(boot.M44.CAP1.500.y, quantile, 0.025)) UCI.M44.CAP1.500.y=unlist(lapply(boot.M44.CAP1.500.y, quantile, 0.975)) mean.M40.CAP1.500.y=tapply(scores(Timoclea.ind.cap.0.40cm.M40)$sites[,1],temp.M40, mean) boot.M40.CAP1.500.y=tapply(scores(Timoclea.ind.cap.0.40cm.M40)$sites[,1],temp.M40, bootci) LCI.M40.CAP1.500.y=unlist(lapply(boot.M40.CAP1.500.y, quantile, 0.025)) UCI.M40.CAP1.500.y=unlist(lapply(boot.M40.CAP1.500.y, quantile, 0.975)) plot(as.numeric(names(mean.M44.CAP1.500.y)),mean.M44.CAP1.500.y, pch=16, ylim=c(-1.5,3.5), xlim=c(0,5000), frame=F, xlab="Postmortem age (y)", ylab="Preservation axis CAP 1", main="aggregated-M44") segments(x0=as.numeric(names(mean.M44.CAP1.500.y)), x1=as.numeric(names(mean.M44.CAP1.500.y)), y0 = LCI.M44.CAP1.500.y, y1 = UCI.M44.CAP1.500.y, lty=1) plot(as.numeric(names(mean.M40.CAP1.500.y)),mean.M40.CAP1.500.y, ylim=c(-2.5,4), xlim=c(0,5000), pch=16, frame=F, xlab="Postmortem age (y)", ylab="Preservation axis CAP 1", main="aggregated-M40") segments(x0=as.numeric(names(mean.M40.CAP1.500.y)), x1=as.numeric(names(mean.M40.CAP1.500.y)), y0 = LCI.M40.CAP1.500.y, y1 = UCI.M40.CAP1.500.y, lty=1) cor.test(as.numeric(names(mean.M44.CAP1.500.y)),mean.M44.CAP1.500.y) cor.test(as.numeric(names(mean.M40.CAP1.500.y)),mean.M40.CAP1.500.y) ############################################################################# #SUPPLEMENTARY FIGURE 1 - PRINCIPAL COORDINATE ANALYSES FOR STRATIGRAPHIC UNITS ############################################################################# par(mfrow=c(2,2)) par(mar=c(4,4,2,1)) select=c(which(Timoclea.dated.depth.10cm<21 & Timoclea.dated.core=="Brijuni 44"),which(Timoclea.dated.depth.10cm<21 & Timoclea.dated.core=="Brijuni 40")) select.1=select Timoclea.ind.distance.top.M40.M44=vegdist(Timoclea.ind.scores[select,],method="manhattan") Timoclea.ind.pco.top.M40.M44=cmdscale(Timoclea.ind.distance.top.M40.M44, eig=T) plot(Timoclea.ind.pco.top.M40.M44$points[,1],Timoclea.ind.pco.top.M40.M44$points[,2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, type="p", main="Surface (top 20 cm)", cex.main=0.9, xlim=c(-4,7), ylim=c(-3,3)) envi=envfit(Timoclea.ind.pco.top.M40.M44$points, Timoclea.ind.scores[select,]) plot(envi, add=T, cex=1, col="gray31", arrow.mul=c(3,3)) points(Timoclea.ind.pco.top.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,1],Timoclea.ind.pco.top.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,2], pch=21, bg="gray") points(Timoclea.ind.pco.top.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,1],Timoclea.ind.pco.top.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,2], pch=21, bg="black") select=c(which(Timoclea.dated.depth.10cm>20 & Timoclea.dated.depth.10cm<90 & Timoclea.dated.core=="Brijuni 44"), which(Timoclea.dated.depth.10cm>20 & Timoclea.dated.depth.10cm <41 & Timoclea.dated.core=="Brijuni 40")) select.2=select Timoclea.ind.distance.HST.M40.M44=vegdist(Timoclea.ind.scores[select,],method="euclidean") Timoclea.ind.pco.HST.M40.M44=cmdscale(Timoclea.ind.distance.HST.M40.M44, eig=T) plot(Timoclea.ind.pco.HST.M40.M44$points[,1],Timoclea.ind.pco.HST.M40.M44$points[,2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, type="p", main="HST", cex.main=0.9, xlim=c(-1.5,1.5), ylim=c(-1.5,1.5)) envi=envfit(Timoclea.ind.pco.HST.M40.M44$points, Timoclea.ind.scores[select,]) plot(envi, add=T, cex=1, col="gray31", arrow.mul=c(1,1)) points(Timoclea.ind.pco.HST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,1],Timoclea.ind.pco.HST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,2], pch=21, bg="gray") points(Timoclea.ind.pco.HST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,1],Timoclea.ind.pco.HST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,2], pch=21, bg="black") select=c(which(Timoclea.dated.depth.10cm>89 & Timoclea.dated.depth.10cm<121 & Timoclea.dated.core=="Brijuni 44"), which(Timoclea.dated.depth.10cm>41 & Timoclea.dated.depth.10cm <71 & Timoclea.dated.core=="Brijuni 40")) select.3=select Timoclea.ind.distance.MFZ.M40.M44=vegdist(Timoclea.ind.scores[select,],method="euclidean") Timoclea.ind.pco.MFZ.M40.M44=cmdscale(Timoclea.ind.distance.MFZ.M40.M44, eig=T) plot(Timoclea.ind.pco.MFZ.M40.M44$points[,1],Timoclea.ind.pco.MFZ.M40.M44$points[,2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, type="p", main="MFZ", cex.main=0.9, ylim=c(1.5,-1), xlim=c(2,-2)) envi=envfit(Timoclea.ind.pco.MFZ.M40.M44$points, Timoclea.ind.scores[select,]) plot(envi, add=T, cex=1, col="gray31", arrow.mul=c(1,1)) points(Timoclea.ind.pco.MFZ.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,1],Timoclea.ind.pco.MFZ.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,2], pch=21, bg="gray") points(Timoclea.ind.pco.MFZ.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,1],Timoclea.ind.pco.MFZ.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,2], pch=21, bg="black") select=c(which(Timoclea.dated.depth.10cm>121 & Timoclea.dated.core=="Brijuni 44"), which(Timoclea.dated.depth.10cm>71 & Timoclea.dated.core=="Brijuni 40")) select.4=select Timoclea.ind.distance.TST.M40.M44=vegdist(Timoclea.ind.scores[select,],method="euclidean") Timoclea.ind.pco.TST.M40.M44=cmdscale(Timoclea.ind.distance.TST.M40.M44, eig=T) plot(Timoclea.ind.pco.TST.M40.M44$points[,1],Timoclea.ind.pco.TST.M40.M44$points[,2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, type="p", main="TST", cex.main=0.9, xlim=c(-2,2), ylim=c(-1.5,1.5)) envi=envfit(Timoclea.ind.pco.TST.M40.M44$points, Timoclea.ind.scores[select,]) plot(envi, add=T, cex=1, col="gray31", arrow.mul=c(1,1)) points(Timoclea.ind.pco.TST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,1],Timoclea.ind.pco.TST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==0,2], pch=21, bg="gray") points(Timoclea.ind.pco.TST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,1],Timoclea.ind.pco.TST.M40.M44$points[Timoclea.ind.scores[select,"stained"]==1,2], pch=21, bg="black") ########################################################################################################### #CORRELATIONS AMONG ALTERATION VARIABLES AT THE LEVEL OF INDIVIDUAL VALVES ########################################################################################################### mydata=cbind(Timoclea.ind.pco.top.M40.M44$points[,2], Timoclea.ind.scores[select.1,]) cor(mydata, method="s")[2:9,1] mydata=cbind(-Timoclea.ind.pco.HST.M40.M44$points[,2], Timoclea.ind.scores[select.2,]) cor(mydata, method="s")[2:9,1] mydata=cbind(Timoclea.ind.pco.MFZ.M40.M44$points[,2], Timoclea.ind.scores[select.3,]) cor(mydata, method="s")[2:9,1] mydata=cbind(-Timoclea.ind.pco.TST.M40.M44$points[,2], Timoclea.ind.scores[select.4,]) cor(mydata, method="s")[2:9,1] ########################################################################################################### #INDIVIDUAL-LEVEL TAPHONOMIC CLOCK WITH GENERALIZED LINEAR MODELS ########################################################################################################### age.variable=Timoclea.M44.ages.for.clock[!is.na(Timoclea.M44.ages.for.clock)] depth.variable=Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)] out=glm(Timoclea.M44.ornamentation[!is.na(Timoclea.M44.ages.for.clock)]~age.variable+depth.variable, family="quasibinomial") out=glm(Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)]~age.variable+depth.variable, family="quasibinomial") out=glm(Timoclea.M44.ext.encrustation[!is.na(Timoclea.M44.ages.for.clock)]~age.variable+depth.variable, family="quasibinomial") out=glm(Timoclea.M44.total.encrustation[!is.na(Timoclea.M44.ages.for.clock)]~age.variable+depth.variable, family="quasibinomial") out=glm(Timoclea.M44.bioerosion[!is.na(Timoclea.M44.ages.for.clock)]~age.variable+depth.variable, family="quasibinomial") out=glm(Timoclea.M44.discoloration[!is.na(Timoclea.M44.ages.for.clock)]~age.variable+depth.variable, family="quasibinomial") out=glm(Timoclea.M44.fragmentation[!is.na(Timoclea.M44.ages.for.clock)]~age.variable+depth.variable, family="quasibinomial") out=glm(Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)]~log(age.variable)+depth.variable, family="quasibinomial") rsq(out) boxplot(split(Timoclea.M44.ages.for.clock[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)])) boxplot(split(Timoclea.M44.ages.for.clock[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.total.encrustation[!is.na(Timoclea.M44.ages.for.clock)])) age.variable=log(Timoclea.M40.ages.for.clock) out=glm(Timoclea.M40.ornamentation~Timoclea.M40.ages.for.clock+Timoclea.M40.10cm, family="quasibinomial") out=glm(Timoclea.M40.stained~Timoclea.M40.ages.for.clock+Timoclea.M40.10cm, family="quasibinomial") out=glm(Timoclea.M40.ext.encrustation~Timoclea.M40.ages.for.clock+Timoclea.M40.10cm, family="quasibinomial") out=glm(Timoclea.M40.total.encrustation~Timoclea.M40.ages.for.clock+Timoclea.M40.10cm, family="quasibinomial") out=glm(Timoclea.M40.bioerosion~Timoclea.M40.ages.for.clock+Timoclea.M40.10cm, family="quasibinomial") out=glm(Timoclea.M40.discoloration~Timoclea.M40.ages.for.clock+Timoclea.M40.10cm, family="quasibinomial") out=glm(Timoclea.M40.fragmentation~Timoclea.M40.ages.for.clock+Timoclea.M40.10cm, family="quasibinomial") ############################################################################################################################################## #MANUSCRIPT FIGURE 9 - UPPER PART - AGE DISTRIBUTIONS IN BOXPLOTS FOR ALTERED AND UNALTERED VALVES IN THE UPPER 20 CM ############################################################################################################################################## LIMIT=26 #LIMIT=51 par(mfrow=c(3,2)) par(mar=c(4,4,2,1)) out1=boxplot(c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm20 & Timoclea.M44.5cm<91 & Timoclea.M44.discoloration==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>20 & Timoclea.M40.5cm<91 & Timoclea.M40.discoloration==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>20 & Timoclea.M44.5cm<91 &Timoclea.M44.unworn.encrustation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>20 & Timoclea.M40.5cm<91 & Timoclea.M40.unworn.encrustation==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>20 & Timoclea.M44.5cm<91 &Timoclea.M44.bioerosion==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>20 & Timoclea.M40.5cm<91 & Timoclea.M40.bioerosion==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>20 & Timoclea.M44.5cm<91 &Timoclea.M44.dissolution==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>20 & Timoclea.M40.5cm<91 & Timoclea.M40.dissolution==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>20 & Timoclea.M44.5cm<91 &Timoclea.M44.ornamentation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>20 & Timoclea.M40.5cm<91 & Timoclea.M40.ornamentation==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>20 & Timoclea.M44.5cm<91 &Timoclea.M44.stained==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>20 & Timoclea.M40.5cm<91 & Timoclea.M40.stained==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>20 & Timoclea.M44.5cm<91 &Timoclea.M44.cementation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>20 & Timoclea.M40.5cm<91 & Timoclea.M40.cementation==1]), las=2, range=0, names=c("Discoloration","Encrustation","Bioerosion","Dissolution","Wear","Staining","Cement"), col="gray", ylab="Postmortem age (yrs)", frame=F) boxplot(c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>90 & Timoclea.M44.5cm<121 & Timoclea.M44.discoloration==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>90 & Timoclea.M40.5cm<121 & Timoclea.M40.discoloration==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>90 & Timoclea.M44.5cm<121 &Timoclea.M44.unworn.encrustation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>90 & Timoclea.M40.5cm<121 & Timoclea.M40.unworn.encrustation==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>90 & Timoclea.M44.5cm<121 &Timoclea.M44.bioerosion==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>90 & Timoclea.M40.5cm<121 & Timoclea.M40.bioerosion==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>90 & Timoclea.M44.5cm<121 &Timoclea.M44.dissolution==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>90 & Timoclea.M40.5cm<121 & Timoclea.M40.dissolution==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>90 & Timoclea.M44.5cm<121 &Timoclea.M44.ornamentation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>90 & Timoclea.M40.5cm<121 & Timoclea.M40.ornamentation==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>90 & Timoclea.M44.5cm<121 &Timoclea.M44.stained==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>90 & Timoclea.M40.5cm<121 & Timoclea.M40.stained==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>90 & Timoclea.M44.5cm<121 &Timoclea.M44.cementation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>90 & Timoclea.M40.5cm<121 & Timoclea.M40.cementation==1]), las=2, range=0, names=c("Discoloration","Encrustation","Bioerosion","Dissolution","Wear","Staining","Cement"), col="gray", ylab="Postmortem age (yrs)", frame=F) boxplot(c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>120 & Timoclea.M44.discoloration==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>120 & Timoclea.M40.discoloration==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>120 &Timoclea.M44.unworn.encrustation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>120 & Timoclea.M40.unworn.encrustation==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>120 &Timoclea.M44.bioerosion==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>120 & Timoclea.M40.bioerosion==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>120 &Timoclea.M44.dissolution==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>120 & Timoclea.M40.dissolution==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>120 &Timoclea.M44.ornamentation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>120 & Timoclea.M40.ornamentation==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>120 &Timoclea.M44.stained==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>120 & Timoclea.M40.stained==1]), c(Timoclea.M44.ages.for.clock[Timoclea.M44.5cm>120 &Timoclea.M44.cementation==1],Timoclea.M40.ages.for.clock[Timoclea.M40.5cm>120 & Timoclea.M40.cementation==1]), las=2, range=0, names=c("Discoloration","Encrustation","Bioerosion","Dissolution","Wear","Staining","Cement"), col="gray", ylab="Postmortem age (yrs)", frame=F) ############################################################################################################################################## #AGE DISTRIBUTIONS FOR COMBINATIONS OF ALTERATION VARIABLES - NOT IN THE MANUSCRIPT ############################################################################################################################################## par(mfrow=c(3,3)) par(mar=c(4,4,2,1)) hist(c(Timoclea.M44.ages.for.clock,Timoclea.M44.ages.for.clock), breaks=seq(0,12000,by=250), col="gray", main="All valves", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.ornamentation==1],Timoclea.M44.ages.for.clock[Timoclea.M40.ornamentation==1]), breaks=seq(0,12000,by=250), col="gray", main="All worn valves ", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.stained==1],Timoclea.M44.ages.for.clock[Timoclea.M40.stained==1]), breaks=seq(0,12000,by=250), col="gray", main="All stained valves ", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.stained==1 & Timoclea.M44.unworn.encrustation==1],Timoclea.M44.ages.for.clock[Timoclea.M40.stained==1 & Timoclea.M40.unworn.encrustation==1]), breaks=seq(0,12000,by=250), col="gray", main="Staining with fresh encrusters", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.stained==1 & Timoclea.M44.worn.encrustation==1],Timoclea.M44.ages.for.clock[Timoclea.M40.stained==1 & Timoclea.M40.worn.encrustation==1]), breaks=seq(0,12000,by=250), col="gray", main="Staining with worn encrusters ", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.stained==1 & Timoclea.M44.total.encrustation==0],Timoclea.M44.ages.for.clock[Timoclea.M40.stained==1 & Timoclea.M40.total.encrustation==0]), breaks=seq(0,12000,by=250), col="gray", main="Staining without encrusters ", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.stained==0 & Timoclea.M44.total.encrustation==1],Timoclea.M44.ages.for.clock[Timoclea.M40.stained==0 & Timoclea.M40.total.encrustation==1]), breaks=seq(0,12000,by=250), col="gray", main="Not-Stained with encrusters ", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.ornamentation==1 & Timoclea.M44.total.encrustation==0],Timoclea.M44.ages.for.clock[Timoclea.M40.ornamentation==1 & Timoclea.M40.total.encrustation==0]), breaks=seq(0,12000,by=250), col="gray", main="Worn without encrusters ", xlab="Postmortem age (yrs)") hist(c(Timoclea.M44.ages.for.clock[Timoclea.M44.ornamentation==1 & Timoclea.M44.total.encrustation==1],Timoclea.M44.ages.for.clock[Timoclea.M40.ornamentation==1 & Timoclea.M40.total.encrustation==1]), breaks=seq(0,12000,by=250), col="gray", main="Worn with encrusters ", xlab="Postmortem age (yrs)") ########################################################################################################### #RESPONSES OF TAPHONOMIC ALTERATION TO POSTMORTEM AGE - INDIVIDUAL-LEVEL CLOCK, ESTIMATED WITH GLM ########################################################################################################### #M44 - GENERALIZED LINEAR MODELS ########################################################################################################### #PARTITIONING DOES NOT WORK BECAUSE OF COLLINEARITY BETWEEN AGE AND SEDIMENT DEPTH age.variable.M44=Timoclea.M44.ages.for.clock[!is.na(Timoclea.M44.ages.for.clock)] age.ord.M44=age.variable.M44[order(age.variable.M44)] depth.variable.M44=Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)] stained.variable.M44=Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)] ext.encrusted.variable.M44=Timoclea.M44.ext.encrustation[!is.na(Timoclea.M44.ages.for.clock)] int.encrusted.variable.M44=Timoclea.M44.int.encrustation[!is.na(Timoclea.M44.ages.for.clock)] total.encrusted.variable.M44=ifelse(ext.encrusted.variable.M44 == 1 | int.encrusted.variable.M44 == 1, 1, 0) ornamentation.variable.M44=Timoclea.M44.ornamentation[!is.na(Timoclea.M44.ages.for.clock)] bioerosion.variable.M44=Timoclea.M44.bioerosion[!is.na(Timoclea.M44.ages.for.clock)] cementation.variable.M44=Timoclea.M44.cementation[!is.na(Timoclea.M44.ages.for.clock)] dissolution.variable.M44=Timoclea.M44.dissolution[!is.na(Timoclea.M44.ages.for.clock)] discoloration.variable.M44=Timoclea.M44.discoloration[!is.na(Timoclea.M44.ages.for.clock)] worn.encrusted.variable.M44=Timoclea.M44.worn.encrustation[!is.na(Timoclea.M44.ages.for.clock)] unworn.encrusted.variable.M44=Timoclea.M44.unworn.encrustation[!is.na(Timoclea.M44.ages.for.clock)] depth.ord.M44=depth.variable.M44[order(age.variable.M44)] stained.ord.M44=stained.variable.M44[order(age.variable.M44)] ext.encrusted.ord.M44=ext.encrusted.variable.M44[order(age.variable.M44)] total.encrusted.ord.M44=total.encrusted.variable.M44[order(age.variable.M44)] ornamentation.ord.M44=ornamentation.variable.M44[order(age.variable.M44)] bioerosion.ord.M44=bioerosion.variable.M44[order(age.variable.M44)] cementation.ord.M44=cementation.variable.M44[order(age.variable.M44)] dissolution.ord.M44=dissolution.variable.M44[order(age.variable.M44)] discoloration.ord.M44=discoloration.variable.M44[order(age.variable.M44)] worn.encrusted.ord.M44=worn.encrusted.variable.M44[order(age.variable.M44)] unworn.encrusted.ord.M44=unworn.encrusted.variable.M44[order(age.variable.M44)] sum.ord.M44=discoloration.ord.M44+dissolution.ord.M44+cementation.ord.M44+bioerosion.ord.M44+ornamentation.ord.M44+stained.ord.M44+worn.encrusted.ord.M44+unworn.encrusted.ord.M44 sum.ord.M44=discoloration.ord.M44+dissolution.ord.M44+cementation.ord.M44+bioerosion.ord.M44+ornamentation.ord.M44+stained.ord.M44+total.encrusted.ord.M44 options(scipen = 999) output.glm.M44=array(NA, dim=c(7,12)) rownames(output.glm.M44)=c("M44-staining","M44-encrustation","M44-bioerosion", "M44-ornamentation loss","M44-dissolution","M44-cementation","M44-discoloration") colnames(output.glm.M44)=c("0-20 cm","p","R2","20-90 cm","p","R2","90-120 cm","p","R2", "120-160 cm", "p","R2") out=glm(stained.ord.M44[depth.ord.M44<21]~age.ord.M44[depth.ord.M44<21], family="quasibinomial") pred.stained.M44.0.20cm=predict.glm(out, predict="response") output.glm.M44[1,1:3]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(total.encrusted.ord.M44[depth.ord.M44<21]~age.ord.M44[depth.ord.M44<21], family="quasibinomial") pred.total.encrusted.M44.0.20cm=predict.glm(out, predict="response") output.glm.M44[2,1:3]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(bioerosion.ord.M44[depth.ord.M44<21]~age.ord.M44[depth.ord.M44<21], family="quasibinomial") pred.bioerosion.M44.0.20cm=predict.glm(out, predict="response") output.glm.M44[3,1:3]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(ornamentation.ord.M44[depth.ord.M44<21]~age.ord.M44[depth.ord.M44<21], family="quasibinomial") pred.ornamentation.M44.0.20cm=predict.glm(out, predict="response") output.glm.M44[4,1:3]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(dissolution.ord.M44[depth.ord.M44<21]~age.ord.M44[depth.ord.M44<21], family="quasibinomial") pred.dissolution.M44.0.20cm=predict.glm(out, predict="response") output.glm.M44[5,1:3]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(cementation.ord.M44[depth.ord.M44<21]~age.ord.M44[depth.ord.M44<21], family="quasibinomial") pred.cementation.M44.0.20cm=predict.glm(out, predict="response") output.glm.M44[6,1:3]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(discoloration.ord.M44[depth.ord.M44<21]~age.ord.M44[depth.ord.M44<21], family="quasibinomial") pred.discoloration.M44.0.20cm=predict.glm(out, predict="response") output.glm.M44[7,1:3]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(stained.ord.M44[depth.ord.M44>20 & depth.ord.M44<89]~age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89], family="quasibinomial") pred.stained.M44.20.90cm=predict.glm(out, predict="response") output.glm.M44[1,4:6]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(total.encrusted.ord.M44[depth.ord.M44>20 & depth.ord.M44<89]~age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89], family="quasibinomial") pred.total.encrusted.M44.20.90cm=predict.glm(out, predict="response") output.glm.M44[2,4:6]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(bioerosion.ord.M44[depth.ord.M44>20 & depth.ord.M44<89]~age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89], family="quasibinomial") pred.bioerosion.M44.20.90cm=predict.glm(out, predict="response") output.glm.M44[3,4:6]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(ornamentation.ord.M44[depth.ord.M44>20 & depth.ord.M44<89]~age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89], family="quasibinomial") pred.ornamentation.M44.20.90cm=predict.glm(out, predict="response") output.glm.M44[4,4:6]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(dissolution.ord.M44[depth.ord.M44>20 & depth.ord.M44<89]~age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89], family="quasibinomial") pred.dissolution.M44.20.90cm=predict.glm(out, predict="response") output.glm.M44[5,4:6]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(cementation.ord.M44[depth.ord.M44>20 & depth.ord.M44<89]~age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89], family="quasibinomial") pred.cementation.M44.20.90cm=predict.glm(out, predict="response") output.glm.M44[6,4:6]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(discoloration.ord.M44[depth.ord.M44>20 & depth.ord.M44<89]~age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89], family="quasibinomial") pred.discoloration.M44.20.90cm=predict.glm(out, predict="response") output.glm.M44[7,4:6]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(stained.ord.M44[depth.ord.M44>89 & depth.ord.M44<121]~age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121], family="quasibinomial") pred.stained.M44.90.120cm=predict.glm(out, predict="response") output.glm.M44[1,7:9]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(total.encrusted.ord.M44[depth.ord.M44>89 & depth.ord.M44<121]~age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121], family="quasibinomial") pred.total.encrusted.M44.90.120cm=predict.glm(out, predict="response") output.glm.M44[2,7:9]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(bioerosion.ord.M44[depth.ord.M44>89 & depth.ord.M44<121]~age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121], family="quasibinomial") pred.bioerosion.M44.90.120cm=predict.glm(out, predict="response") output.glm.M44[3,7:9]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(ornamentation.ord.M44[depth.ord.M44>89 & depth.ord.M44<121]~age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121], family="quasibinomial") pred.ornamentation.M44.90.120cm=predict.glm(out, predict="response") output.glm.M44[4,7:9]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(dissolution.ord.M44[depth.ord.M44>89 & depth.ord.M44<121]~age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121], family="quasibinomial") pred.dissolution.M44.90.120cm=predict.glm(out, predict="response") output.glm.M44[5,7:9]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(cementation.ord.M44[depth.ord.M44>89 & depth.ord.M44<121]~age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121], family="quasibinomial") pred.cementation.M44.90.120cm=predict.glm(out, predict="response") output.glm.M44[6,7:9]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(discoloration.ord.M44[depth.ord.M44>89 & depth.ord.M44<121]~age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121], family="quasibinomial") pred.discoloration.M44.90.120cm=predict.glm(out, predict="response") output.glm.M44[7,7:9]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(stained.ord.M44[depth.ord.M44>121]~age.ord.M44[depth.ord.M44>121], family="quasibinomial") pred.stained.M44.120.160cm=predict.glm(out, predict="response") output.glm.M44[1,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(total.encrusted.ord.M44[depth.ord.M44>121]~age.ord.M44[depth.ord.M44>121], family="quasibinomial") pred.total.encrusted.M44.120.160cm=predict.glm(out, predict="response") output.glm.M44[2,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(bioerosion.ord.M44[depth.ord.M44>121]~age.ord.M44[depth.ord.M44>121], family="quasibinomial") pred.bioerosion.M44.120.160cm=predict.glm(out, predict="response") output.glm.M44[3,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(ornamentation.ord.M44[depth.ord.M44>121]~age.ord.M44[depth.ord.M44>121], family="quasibinomial") pred.ornamentation.M44.120.160cm=predict.glm(out, predict="response") output.glm.M44[4,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(dissolution.ord.M44[depth.ord.M44>121]~age.ord.M44[depth.ord.M44>121], family="quasibinomial") pred.dissolution.M44.120.160cm=predict.glm(out, predict="response") output.glm.M44[5,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(cementation.ord.M44[depth.ord.M44>121]~age.ord.M44[depth.ord.M44>121], family="quasibinomial") pred.cementation.M44.120.160cm=predict.glm(out, predict="response") output.glm.M44[6,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(discoloration.ord.M44[depth.ord.M44>121]~age.ord.M44[depth.ord.M44>121], family="quasibinomial") pred.discoloration.M44.120.160cm=predict.glm(out, predict="response") output.glm.M44[7,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) output.glm.M44[,c(2,5,8,9)]=round(output.glm.M44[,c(2,5,8,9)], digits=3) output.glm.M44[,c(1,4,7,10)]=round(output.glm.M44[,c(1,4,7,10)], digits=5) output.glm.M44[,c(3,6,9,12)]=round(output.glm.M44[,c(3,6,9,12)], digits=2) ########################################################################################################### #M40 - GENERALIZED LINEAR MODELS ########################################################################################################### age.variable.M40=Timoclea.M40.ages.for.clock[!is.na(Timoclea.M40.ages.for.clock)] age.ord.M40=age.variable.M40[order(age.variable.M40)] depth.variable.M40=Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)] stained.variable.M40=Timoclea.M40.stained[!is.na(Timoclea.M40.ages.for.clock)] ext.encrusted.variable.M40=Timoclea.M40.ext.encrustation[!is.na(Timoclea.M40.ages.for.clock)] int.encrusted.variable.M40=Timoclea.M40.int.encrustation[!is.na(Timoclea.M40.ages.for.clock)] total.encrusted.variable.M40=ifelse(ext.encrusted.variable.M40 == 1 | int.encrusted.variable.M40 == 1, 1, 0) ornamentation.variable.M40=Timoclea.M40.ornamentation[!is.na(Timoclea.M40.ages.for.clock)] bioerosion.variable.M40=Timoclea.M40.bioerosion[!is.na(Timoclea.M40.ages.for.clock)] cementation.variable.M40=Timoclea.M40.cementation[!is.na(Timoclea.M40.ages.for.clock)] dissolution.variable.M40=Timoclea.M40.dissolution[!is.na(Timoclea.M40.ages.for.clock)] discoloration.variable.M40=Timoclea.M40.discoloration[!is.na(Timoclea.M40.ages.for.clock)] worn.encrusted.variable.M40=Timoclea.M40.worn.encrustation[!is.na(Timoclea.M40.ages.for.clock)] unworn.encrusted.variable.M40=Timoclea.M40.unworn.encrustation[!is.na(Timoclea.M40.ages.for.clock)] depth.ord.M40=depth.variable.M40[order(age.variable.M40)] stained.ord.M40=stained.variable.M40[order(age.variable.M40)] ext.encrusted.ord.M40=ext.encrusted.variable.M40[order(age.variable.M40)] total.encrusted.ord.M40=total.encrusted.variable.M40[order(age.variable.M40)] ornamentation.ord.M40=ornamentation.variable.M40[order(age.variable.M40)] bioerosion.ord.M40=bioerosion.variable.M40[order(age.variable.M40)] cementation.ord.M40=cementation.variable.M40[order(age.variable.M40)] dissolution.ord.M40=dissolution.variable.M40[order(age.variable.M40)] discoloration.ord.M40=discoloration.variable.M40[order(age.variable.M40)] worn.encrusted.ord.M40=worn.encrusted.variable.M40[order(age.variable.M40)] unworn.encrusted.ord.M40=unworn.encrusted.variable.M40[order(age.variable.M40)] sum.ord.M40=discoloration.ord.M40+dissolution.ord.M40+cementation.ord.M40+bioerosion.ord.M40+ornamentation.ord.M40+stained.ord.M40+worn.encrusted.ord.M40+unworn.encrusted.ord.M40 sum.ord.M40=discoloration.ord.M40+dissolution.ord.M40+cementation.ord.M40+bioerosion.ord.M40+ornamentation.ord.M40+stained.ord.M40+total.encrusted.ord.M40 options(scipen = 999) output.glm.M40=array(NA, dim=c(7,12)) rownames(output.glm.M40)=c("M40-staining","M40-encrustation","M40-bioerosion", "M40-ornamentation loss","M40-dissolution","M40-cementation","M40-discoloration") colnames(output.glm.M40)=c("0-20 cm","p","R2","20-45 cm","p","R2","45-70 cm","p","R2", "70-120 cm", "p","R2") out=glm(stained.ord.M40[depth.ord.M40<21]~age.ord.M40[depth.ord.M40<21], family="quasibinomial") pred.stained.M40.0.20cm=predict.glm(out, predict="response") output.glm.M40[1,1:3]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(total.encrusted.ord.M40[depth.ord.M40<21]~age.ord.M40[depth.ord.M40<21], family="quasibinomial") pred.total.encrusted.M40.0.20cm=predict.glm(out, predict="response") output.glm.M40[2,1:3]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(bioerosion.ord.M40[depth.ord.M40<21]~age.ord.M40[depth.ord.M40<21], family="quasibinomial") pred.bioerosion.M40.0.20cm=predict.glm(out, predict="response") output.glm.M40[3,1:3]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(ornamentation.ord.M40[depth.ord.M40<21]~age.ord.M40[depth.ord.M40<21], family="quasibinomial") pred.ornamentation.M40.0.20cm=predict.glm(out, predict="response") output.glm.M40[4,1:3]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(dissolution.ord.M40[depth.ord.M40<21]~age.ord.M40[depth.ord.M40<21], family="quasibinomial") pred.dissolution.M40.0.20cm=predict.glm(out, predict="response") output.glm.M40[5,1:3]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(cementation.ord.M40[depth.ord.M40<21]~age.ord.M40[depth.ord.M40<21], family="quasibinomial") pred.cementation.M40.0.20cm=predict.glm(out, predict="response") output.glm.M40[6,1:3]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(discoloration.ord.M40[depth.ord.M40<21]~age.ord.M40[depth.ord.M40<21], family="quasibinomial") pred.discoloration.M40.0.20cm=predict.glm(out, predict="response") output.glm.M40[7,1:3]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(stained.ord.M40[depth.ord.M40>20 & depth.ord.M40<89]~age.ord.M40[depth.ord.M40>20 & depth.ord.M40<89], family="quasibinomial") pred.stained.M40.20.90cm=predict.glm(out, predict="response") output.glm.M40[1,4:6]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(total.encrusted.ord.M40[depth.ord.M40>20 & depth.ord.M40<89]~age.ord.M40[depth.ord.M40>20 & depth.ord.M40<89], family="quasibinomial") pred.total.encrusted.M40.20.90cm=predict.glm(out, predict="response") output.glm.M40[2,4:6]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(bioerosion.ord.M40[depth.ord.M40>20 & depth.ord.M40<89]~age.ord.M40[depth.ord.M40>20 & depth.ord.M40<89], family="quasibinomial") pred.bioerosion.M40.20.90cm=predict.glm(out, predict="response") output.glm.M40[3,4:6]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(ornamentation.ord.M40[depth.ord.M40>20 & depth.ord.M40<89]~age.ord.M40[depth.ord.M40>20 & depth.ord.M40<89], family="quasibinomial") pred.ornamentation.M40.20.90cm=predict.glm(out, predict="response") output.glm.M40[4,4:6]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(dissolution.ord.M40[depth.ord.M40>20 & depth.ord.M40<89]~age.ord.M40[depth.ord.M40>20 & depth.ord.M40<89], family="quasibinomial") pred.dissolution.M40.20.90cm=predict.glm(out, predict="response") output.glm.M40[5,4:6]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(cementation.ord.M40[depth.ord.M40>20 & depth.ord.M40<89]~age.ord.M40[depth.ord.M40>20 & depth.ord.M40<89], family="quasibinomial") pred.cementation.M40.20.90cm=predict.glm(out, predict="response") output.glm.M40[6,4:6]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(discoloration.ord.M40[depth.ord.M40>20 & depth.ord.M40<89]~age.ord.M40[depth.ord.M40>20 & depth.ord.M40<89], family="quasibinomial") pred.discoloration.M40.20.90cm=predict.glm(out, predict="response") output.glm.M40[7,4:6]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(stained.ord.M40[depth.ord.M40>89 & depth.ord.M40<121]~age.ord.M40[depth.ord.M40>89 & depth.ord.M40<121], family="quasibinomial") pred.stained.M40.90.120cm=predict.glm(out, predict="response") output.glm.M40[1,7:9]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(total.encrusted.ord.M40[depth.ord.M40>89 & depth.ord.M40<121]~age.ord.M40[depth.ord.M40>89 & depth.ord.M40<121], family="quasibinomial") pred.total.encrusted.M40.90.120cm=predict.glm(out, predict="response") output.glm.M40[2,7:9]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(bioerosion.ord.M40[depth.ord.M40>89 & depth.ord.M40<121]~age.ord.M40[depth.ord.M40>89 & depth.ord.M40<121], family="quasibinomial") pred.bioerosion.M40.90.120cm=predict.glm(out, predict="response") output.glm.M40[3,7:9]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(ornamentation.ord.M40[depth.ord.M40>89 & depth.ord.M40<121]~age.ord.M40[depth.ord.M40>89 & depth.ord.M40<121], family="quasibinomial") pred.ornamentation.M40.90.120cm=predict.glm(out, predict="response") output.glm.M40[4,7:9]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(dissolution.ord.M40[depth.ord.M40>89 & depth.ord.M40<121]~age.ord.M40[depth.ord.M40>89 & depth.ord.M40<121], family="quasibinomial") pred.dissolution.M40.90.120cm=predict.glm(out, predict="response") output.glm.M40[5,7:9]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(cementation.ord.M40[depth.ord.M40>89 & depth.ord.M40<121]~age.ord.M40[depth.ord.M40>89 & depth.ord.M40<121], family="quasibinomial") pred.cementation.M40.90.120cm=predict.glm(out, predict="response") output.glm.M40[6,7:9]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(discoloration.ord.M40[depth.ord.M40>89 & depth.ord.M40<121]~age.ord.M40[depth.ord.M40>89 & depth.ord.M40<121], family="quasibinomial") pred.discoloration.M40.90.120cm=predict.glm(out, predict="response") output.glm.M40[7,7:9]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(stained.ord.M40[depth.ord.M40>71]~age.ord.M40[depth.ord.M40>71], family="quasibinomial") pred.stained.M40.120.160cm=predict.glm(out, predict="response") output.glm.M40[1,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(total.encrusted.ord.M40[depth.ord.M40>71]~age.ord.M40[depth.ord.M40>71], family="quasibinomial") pred.total.encrusted.M40.120.160cm=predict.glm(out, predict="response") output.glm.M40[2,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(bioerosion.ord.M40[depth.ord.M40>71]~age.ord.M40[depth.ord.M40>71], family="quasibinomial") pred.bioerosion.M40.120.160cm=predict.glm(out, predict="response") output.glm.M40[3,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(ornamentation.ord.M40[depth.ord.M40>71]~age.ord.M40[depth.ord.M40>71], family="quasibinomial") pred.ornamentation.M40.120.160cm=predict.glm(out, predict="response") output.glm.M40[4,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(dissolution.ord.M40[depth.ord.M40>71]~age.ord.M40[depth.ord.M40>71], family="quasibinomial") pred.dissolution.M40.120.160cm=predict.glm(out, predict="response") output.glm.M40[5,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(cementation.ord.M40[depth.ord.M40>71]~age.ord.M40[depth.ord.M40>71], family="quasibinomial") pred.cementation.M40.120.160cm=predict.glm(out, predict="response") output.glm.M40[6,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) out=glm(discoloration.ord.M40[depth.ord.M40>71]~age.ord.M40[depth.ord.M40>71], family="quasibinomial") pred.discoloration.M40.120.160cm=predict.glm(out, predict="response") output.glm.M40[7,10:12]=c(coef(summary(out))[2,c(1,4)], rsq(out)) output.glm.M40[,c(2,5,8,9)]=round(output.glm.M40[,c(2,5,8,9)], digits=3) output.glm.M40[,c(1,4,7,10)]=round(output.glm.M40[,c(1,4,7,10)], digits=5) output.glm.M40[,c(3,6,9,12)]=round(output.glm.M40[,c(3,6,9,12)], digits=2) glm.output=rbind(output.glm.M44,output.glm.M40) write.table(glm.output, file="GLM individual-level output.txt") ########################################################################################################### #PLOT GLM RESPONSES FOR M44 ########################################################################################################### par(cex=1.4) par(mfrow=c(2,2)) #exp(x)/(1+exp(x)) plot(age.ord.M44[depth.ord.M44<21],exp(pred.stained.M44.0.20cm)/(1+exp(pred.stained.M44.0.20cm)), ylim=c(0,1), type="n", frame=F, ylab="Proportion", xlab="Time (years)") lines(age.ord.M44[depth.ord.M44<21],exp(pred.stained.M44.0.20cm)/(1+exp(pred.stained.M44.0.20cm)), lwd=2) lines(age.ord.M44[depth.ord.M44<21],exp(pred.total.encrusted.M44.0.20cm)/(1+exp(pred.total.encrusted.M44.0.20cm)), lwd=2, col="black", lty=2) lines(age.ord.M44[depth.ord.M44<21],exp(pred.bioerosion.M44.0.20cm)/(1+exp(pred.bioerosion.M44.0.20cm)), col="black", lwd=2, lty=3) lines(age.ord.M44[depth.ord.M44<21],exp(pred.ornamentation.M44.0.20cm)/(1+exp(pred.ornamentation.M44.0.20cm)), col="gray71", lwd=2) lines(age.ord.M44[depth.ord.M44<21],exp(pred.cementation.M44.0.20cm)/(1+exp(pred.cementation.M44.0.20cm)), col="gray71", lwd=2, lty=2) lines(age.ord.M44[depth.ord.M44<21],exp(pred.dissolution.M44.0.20cm)/(1+exp(pred.dissolution.M44.0.20cm)), col="gray71", lwd=2, lty=3) legend(x="bottomright", col=c("black","black","black","gray71","gray71","gray71"),lwd=2, lty=c(1,2,3,1,2,3), legend=c("staining","encr","bioer","ornam","cement","diss"), bty="n", ncol=2) plot(age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89],exp(pred.stained.M44.20.90cm)/(1+exp(pred.stained.M44.20.90cm)), ylim=c(0,1), type="n", frame=F, ylab="Proportion", xlab="Time (years)") lines(age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89],exp(pred.stained.M44.20.90cm)/(1+exp(pred.stained.M44.20.90cm)), lwd=2) lines(age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89],exp(pred.total.encrusted.M44.20.90cm)/(1+exp(pred.total.encrusted.M44.20.90cm)), lwd=2, col="black", lty=2) lines(age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89],exp(pred.bioerosion.M44.20.90cm)/(1+exp(pred.bioerosion.M44.20.90cm)),col="black", lwd=2, lty=3) lines(age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89],exp(pred.ornamentation.M44.20.90cm)/(1+exp(pred.ornamentation.M44.20.90cm)), col="gray71", lwd=2) lines(age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89],exp(pred.cementation.M44.20.90cm)/(1+exp(pred.cementation.M44.20.90cm)), col="gray71", lwd=2, lty=2) lines(age.ord.M44[depth.ord.M44>20 & depth.ord.M44<89],exp(pred.dissolution.M44.20.90cm)/(1+exp(pred.dissolution.M44.20.90cm)), col="gray71", lwd=2, lty=3) plot(age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121],exp(pred.stained.M44.90.120cm)/(1+exp(pred.stained.M44.90.120cm)), ylim=c(0,1), type="n", frame=F, ylab="Proportion", xlab="Time (years)") lines(age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121],exp(pred.stained.M44.90.120cm)/(1+exp(pred.stained.M44.90.120cm)), lwd=2) lines(age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121],exp(pred.total.encrusted.M44.90.120cm)/(1+exp(pred.total.encrusted.M44.90.120cm)), lwd=2, col="black", lty=2) lines(age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121],exp(pred.bioerosion.M44.90.120cm)/(1+exp(pred.bioerosion.M44.90.120cm)), col="black", lwd=2, lty=3) lines(age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121],exp(pred.ornamentation.M44.90.120cm)/(1+exp(pred.ornamentation.M44.90.120cm)), col="gray71", lwd=2) lines(age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121],exp(pred.cementation.M44.90.120cm)/(1+exp(pred.cementation.M44.90.120cm)), col="gray71", lwd=2, lty=2) lines(age.ord.M44[depth.ord.M44>89 & depth.ord.M44<121],exp(pred.dissolution.M44.90.120cm)/(1+exp(pred.dissolution.M44.90.120cm)), col="gray71", lwd=2, lty=3) plot(age.ord.M44[depth.ord.M44>121],exp(pred.stained.M44.120.160cm)/(1+exp(pred.stained.M44.120.160cm)), ylim=c(0,1), type="n", frame=F, ylab="Proportion", xlab="Time (years)") lines(age.ord.M44[depth.ord.M44>121],exp(pred.stained.M44.120.160cm)/(1+exp(pred.stained.M44.120.160cm)), lwd=2) lines(age.ord.M44[depth.ord.M44>121],exp(pred.total.encrusted.M44.120.160cm)/(1+exp(pred.total.encrusted.M44.120.160cm)), lwd=2, col="black", lty=2) lines(age.ord.M44[depth.ord.M44>121],exp(pred.bioerosion.M44.120.160cm)/(1+exp(pred.bioerosion.M44.120.160cm)), col="black", lwd=2, lty=3) lines(age.ord.M44[depth.ord.M44>121],exp(pred.ornamentation.M44.120.160cm)/(1+exp(pred.ornamentation.M44.120.160cm)), col="gray71", lwd=2) lines(age.ord.M44[depth.ord.M44>121],exp(pred.cementation.M44.120.160cm)/(1+exp(pred.cementation.M44.120.160cm)), col="gray71", lwd=2, lty=2) lines(age.ord.M44[depth.ord.M44>121],exp(pred.dissolution.M44.120.160cm)/(1+exp(pred.dissolution.M44.120.160cm)), col="gray71", lwd=2, lty=3) ############################################################################################### #AGGREGATING AGE COHORTS - FIRST 1000 YEARS TO 100 YEAR COHORTS, THEN 500 YEARS ############################################################################################### TOP.SURFACE=91 COHORT=500 BREAKS=c(0,100,200,300,400,500,600,700,800,900,seq(1000,13000,by=COHORT)) temp.M44=cut(age.ord.M44[depth.ord.M44CUTOFF],mean.M44.bioerosion.500.y[N.M44.stained.500.y>CUTOFF], pch=16, ylim=c(0,1), frame=F, xlab="Postmortem age (y)", ylab="Proportion of altered", xlim=c(0,5000)) lines(as.numeric(names(mean.M44.bioerosion.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.bioerosion.500.y[N.M44.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M44.bioerosion.500.y))[N.M44.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M44.bioerosion.500.y))[N.M44.stained.500.y>CUTOFF], y0 = LCI.M44.bioerosion.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.bioerosion.500.y[N.M44.stained.500.y>CUTOFF], lty=1) legend(x="bottomright", pch=21, pt.bg=c("black","gray71"), legend=c("Bioerosion","Dissolution"), bty="n", cex=1) lines(as.numeric(names(mean.M44.dissolution.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.dissolution.500.y[N.M44.stained.500.y>CUTOFF], col="black") segments(x0=as.numeric(names(mean.M44.dissolution.500.y))[N.M44.stained.500.y>CUTOFF]+50, x1=as.numeric(names(mean.M44.dissolution.500.y))[N.M44.stained.500.y>CUTOFF]+50, y0 = LCI.M44.dissolution.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.dissolution.500.y[N.M44.stained.500.y>CUTOFF], lty=1) points(as.numeric(names(mean.M44.bioerosion.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.bioerosion.500.y[N.M44.stained.500.y>CUTOFF], pch=16, bg="gray71", cex=1.4) points(as.numeric(names(mean.M44.dissolution.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.dissolution.500.y[N.M44.stained.500.y>CUTOFF], pch=21, bg="gray71", cex=1.4) corr=cor.test(as.numeric(names(mean.M44.bioerosion.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.bioerosion.500.y[N.M44.stained.500.y>CUTOFF]) out=glm(mean.M44.bioerosion.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.bioerosion.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") output.M44[7,1:2]=c(corr$estimate,corr$p.value) output.M44[7,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(mean.M44.dissolution.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.dissolution.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M44.dissolution.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.dissolution.500.y[N.M44.stained.500.y>CUTOFF]) output.M44[8,1:2]=c(corr$estimate,corr$p.value) output.M44[8,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) plot(as.numeric(names(mean.M40.bioerosion.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF], ylim=c(0,1), pch=16, frame=F, xlab="Postmortem age (y)",ylab="", xlim=c(0,5000)) lines(as.numeric(names(mean.M40.bioerosion.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M40.bioerosion.500.y))[N.M40.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M40.bioerosion.500.y))[N.M40.stained.500.y>CUTOFF], y0 = LCI.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M40.dissolution.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.dissolution.500.y[N.M40.stained.500.y>CUTOFF], col="black") segments(x0=as.numeric(names(mean.M40.dissolution.500.y[N.M40.stained.500.y>CUTOFF]))+50, x1=as.numeric(names(mean.M40.dissolution.500.y[N.M40.stained.500.y>CUTOFF]))+50, y0 = LCI.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF], lty=1) points(as.numeric(names(mean.M40.bioerosion.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF], pch=16, cex=1.4) points(as.numeric(names(mean.M40.dissolution.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.dissolution.500.y[N.M40.stained.500.y>CUTOFF], pch=21, bg="gray71", cex=1.4) out=glm(mean.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.bioerosion.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M40.bioerosion.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.bioerosion.500.y[N.M40.stained.500.y>CUTOFF]) output.M40[7,1:2]=c(corr$estimate,corr$p.value) output.M40[7,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(mean.M40.dissolution.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.dissolution.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M40.dissolution.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.dissolution.500.y[N.M40.stained.500.y>CUTOFF]) output.M40[8,1:2]=c(corr$estimate,corr$p.value) output.M40[8,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) ################################################ #AVERAGE ACROSS 500 YEARS TOTAL ENCRUSTATION ################################################ plot(as.numeric(names(mean.M44.total.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.total.encrusted.500.y[N.M44.stained.500.y>CUTOFF], type="n", pch=16, ylim=c(0,1), frame=F, xlab="Postmortem age (y)", ylab="Proportion of altered", xlim=c(0,5000)) lines(as.numeric(names(mean.M44.total.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.total.encrusted.500.y[N.M44.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M44.total.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M44.total.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], y0 = LCI.M44.total.encrusted.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.total.encrusted.500.y[N.M44.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M44.ornamentation.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.ornamentation.500.y[N.M44.stained.500.y>CUTOFF],bg="gray71") segments(x0=as.numeric(names(mean.M44.ornamentation.500.y))[N.M44.stained.500.y>CUTOFF]+50, x1=as.numeric(names(mean.M44.ornamentation.500.y))[N.M44.stained.500.y>CUTOFF]+50, y0 = LCI.M44.ornamentation.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.ornamentation.500.y[N.M44.stained.500.y>CUTOFF], lty=1) corr=cor.test(as.numeric(names(mean.M44.total.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.total.encrusted.500.y[N.M44.stained.500.y>CUTOFF]) out=glm(mean.M44.total.encrusted.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.total.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") output.M44[3,1:2]=c(corr$estimate,corr$p.value) output.M44[3,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) legend(x="topleft", pch=c(21,8,21,22), pt.bg=c("gray71","black","white","black"), legend=c("Ornam. loss","All encrusters", "Unworn encrusters","Worn encrusters"), bty="n", cex=1, ncol=2) out=glm(mean.M44.ornamentation.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.ornamentation.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M44.ornamentation.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.ornamentation.500.y[N.M44.stained.500.y>CUTOFF]) output.M44[6,1:2]=c(corr$estimate,corr$p.value) output.M44[6,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) corr=cor.test(as.numeric(names(mean.M44.unworn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.unworn.encrusted.500.y[N.M44.stained.500.y>CUTOFF]) out=glm(mean.M44.unworn.encrusted.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.unworn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") output.M44[4,1:2]=c(corr$estimate,corr$p.value) output.M44[4,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(mean.M44.worn.encrusted.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.worn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M44.worn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.worn.encrusted.500.y[N.M44.stained.500.y>CUTOFF]) output.M44[5,1:2]=c(corr$estimate,corr$p.value) output.M44[5,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) lines(as.numeric(names(mean.M44.unworn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.unworn.encrusted.500.y[N.M44.stained.500.y>CUTOFF],bg="gray71") segments(x0=as.numeric(names(mean.M44.unworn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M44.unworn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], y0 = LCI.M44.unworn.encrusted.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.unworn.encrusted.500.y[N.M44.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M44.worn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.worn.encrusted.500.y[N.M44.stained.500.y>CUTOFF],bg="gray71") segments(x0=as.numeric(names(mean.M44.worn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M44.worn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF], y0 = LCI.M44.worn.encrusted.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.worn.encrusted.500.y[N.M44.stained.500.y>CUTOFF], lty=1) points(as.numeric(names(mean.M44.total.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.total.encrusted.500.y[N.M44.stained.500.y>CUTOFF], pch=8, cex=1.4) points(as.numeric(names(mean.M44.ornamentation.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.ornamentation.500.y[N.M44.stained.500.y>CUTOFF], pch=21, bg="gray71", cex=1.4) points(as.numeric(names(mean.M44.unworn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.unworn.encrusted.500.y[N.M44.stained.500.y>CUTOFF], pch=21, bg="white", cex=1.5) points(as.numeric(names(mean.M44.worn.encrusted.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.worn.encrusted.500.y[N.M44.stained.500.y>CUTOFF], pch=22, bg="black", cex=1.5) plot(as.numeric(names(mean.M40.total.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.total.encrusted.500.y[N.M40.stained.500.y>CUTOFF], ylim=c(0,1), type="n", pch=16, frame=F, xlab="Postmortem age (y)", ylab="", xlim=c(0,5000)) lines(as.numeric(names(mean.M40.total.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.total.encrusted.500.y[N.M40.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M40.total.encrusted.500.y))[N.M40.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M40.total.encrusted.500.y))[N.M40.stained.500.y>CUTOFF], y0 = LCI.M40.total.encrusted.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.total.encrusted.500.y[N.M40.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M40.unworn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.unworn.encrusted.500.y[N.M40.stained.500.y>CUTOFF],bg="gray71") segments(x0=as.numeric(names(mean.M40.unworn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF]+50, x1=as.numeric(names(mean.M40.unworn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF]+50, y0 = LCI.M40.unworn.encrusted.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.unworn.encrusted.500.y[N.M40.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M40.worn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.worn.encrusted.500.y[N.M40.stained.500.y>CUTOFF],bg="gray71") segments(x0=as.numeric(names(mean.M40.worn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M40.worn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF], y0 = LCI.M40.worn.encrusted.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.worn.encrusted.500.y[N.M40.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M40.ornamentation.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.ornamentation.500.y[N.M40.stained.500.y>CUTOFF],bg="gray71") segments(x0=as.numeric(names(mean.M40.ornamentation.500.y))[N.M40.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M40.ornamentation.500.y))[N.M40.stained.500.y>CUTOFF], y0 = LCI.M40.ornamentation.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.ornamentation.500.y[N.M40.stained.500.y>CUTOFF], lty=1) points(as.numeric(names(mean.M40.total.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.total.encrusted.500.y[N.M40.stained.500.y>CUTOFF], pch=8, cex=1.4) points(as.numeric(names(mean.M40.ornamentation.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.ornamentation.500.y[N.M40.stained.500.y>CUTOFF], pch=21, bg="gray71", cex=1.4) points(as.numeric(names(mean.M40.unworn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.unworn.encrusted.500.y[N.M40.stained.500.y>CUTOFF], pch=21, bg="white", cex=1.5) points(as.numeric(names(mean.M40.worn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.worn.encrusted.500.y[N.M40.stained.500.y>CUTOFF], pch=22, bg="black", cex=1.5) corr=cor.test(as.numeric(names(mean.M40.total.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.total.encrusted.500.y[N.M40.stained.500.y>CUTOFF]) out=glm(mean.M40.total.encrusted.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.total.encrusted.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") output.M40[3,1:2]=c(corr$estimate,corr$p.value) output.M40[3,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(mean.M40.ornamentation.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.ornamentation.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M40.ornamentation.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.ornamentation.500.y[N.M40.stained.500.y>CUTOFF]) output.M40[6,1:2]=c(corr$estimate,corr$p.value) output.M40[6,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) corr=cor.test(as.numeric(names(mean.M40.unworn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.unworn.encrusted.500.y[N.M40.stained.500.y>CUTOFF]) out=glm(mean.M40.unworn.encrusted.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.unworn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") output.M40[4,1:2]=c(corr$estimate,corr$p.value) output.M40[4,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(mean.M40.worn.encrusted.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.worn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M40.worn.encrusted.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.worn.encrusted.500.y[N.M40.stained.500.y>CUTOFF]) output.M40[5,1:2]=c(corr$estimate,corr$p.value) output.M40[5,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) ##################################################################################################### #AVERAGE ACROSS 500 YEARS STAINING ##################################################################################################### plot(as.numeric(names(mean.M44.stained.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.stained.500.y[N.M44.stained.500.y>CUTOFF], pch=16, cex=1.2, ylim=c(0,1), frame=F, xlab="Postmortem age (y)", ylab="Proportion of altered", xlim=c(0,5000)) lines(as.numeric(names(mean.M44.stained.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.stained.500.y[N.M44.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M44.stained.500.y))[N.M44.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M44.stained.500.y))[N.M44.stained.500.y>CUTOFF], y0 = LCI.M44.stained.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.stained.500.y[N.M44.stained.500.y>CUTOFF], lty=1) legend(x="topleft", pch=21, pt.bg=c("black","gray71","white"), legend=c("Staining","Cementation","Discoloration"), bty="n", cex=1) corr=cor.test(as.numeric(names(mean.M44.stained.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.stained.500.y[N.M44.stained.500.y>CUTOFF]) out=glm(mean.M44.stained.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.stained.500.y)[N.M44.stained.500.y>CUTOFF]), family="quasibinomial") output.M44[1,1:2]=c(corr$estimate,corr$p.value) output.M44[1,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) lines(as.numeric(names(mean.M44.cementation.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.cementation.500.y[N.M44.stained.500.y>CUTOFF], col="black") segments(x0=as.numeric(names(mean.M44.cementation.500.y))[N.M44.stained.500.y>CUTOFF]+50, x1=as.numeric(names(mean.M44.cementation.500.y))[N.M44.stained.500.y>CUTOFF]+50, y0 = LCI.M44.cementation.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.cementation.500.y[N.M44.stained.500.y>CUTOFF], lty=1) corr=cor.test(as.numeric(names(mean.M44.cementation.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.cementation.500.y[N.M44.stained.500.y>CUTOFF]) out=glm(mean.M44.cementation.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.cementation.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") output.M44[2,1:2]=c(corr$estimate,corr$p.value) output.M44[2,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) lines(as.numeric(names(mean.M44.discoloration.500.y))[N.M44.stained.500.y>CUTOFF]+100,mean.M44.discoloration.500.y[N.M44.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M44.discoloration.500.y))[N.M44.stained.500.y>CUTOFF]+100, x1=as.numeric(names(mean.M44.discoloration.500.y))[N.M44.stained.500.y>CUTOFF]+100, y0 = LCI.M44.discoloration.500.y[N.M44.stained.500.y>CUTOFF], y1 = UCI.M44.discoloration.500.y[N.M44.stained.500.y>CUTOFF], lty=1) points(as.numeric(names(mean.M44.discoloration.500.y))[N.M44.stained.500.y>CUTOFF]+100,mean.M44.discoloration.500.y[N.M44.stained.500.y>CUTOFF], pch=21, bg="white", cex=1.4) points(as.numeric(names(mean.M44.stained.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.stained.500.y[N.M44.stained.500.y>CUTOFF], pch=16, cex=1.4) points(as.numeric(names(mean.M44.cementation.500.y))[N.M44.stained.500.y>CUTOFF]+50,mean.M44.cementation.500.y[N.M44.stained.500.y>CUTOFF], pch=21, bg="gray71", cex=1.4) plot(as.numeric(names(mean.M40.stained.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.stained.500.y[N.M40.stained.500.y>CUTOFF], ylim=c(0,1), pch=16, cex=1.2, frame=F, xlab="Postmortem age (y)", ylab="", xlim=c(0,5000)) lines(as.numeric(names(mean.M40.stained.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.stained.500.y[N.M40.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M40.stained.500.y))[N.M40.stained.500.y>CUTOFF], x1=as.numeric(names(mean.M40.stained.500.y))[N.M40.stained.500.y>CUTOFF], y0 = LCI.M40.stained.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.stained.500.y[N.M40.stained.500.y>CUTOFF], lty=1) corr=cor.test(as.numeric(names(mean.M44.discoloration.500.y))[N.M44.stained.500.y>CUTOFF],mean.M44.discoloration.500.y[N.M44.stained.500.y>CUTOFF]) out=glm(mean.M44.discoloration.500.y[N.M44.stained.500.y>CUTOFF]~as.numeric(names(mean.M44.discoloration.500.y))[N.M44.stained.500.y>CUTOFF], family="quasibinomial") output.M44[9,1:2]=c(corr$estimate,corr$p.value) output.M44[9,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) corr=cor.test(as.numeric(names(mean.M40.stained.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.stained.500.y[N.M40.stained.500.y>CUTOFF]) out=glm(mean.M40.stained.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.stained.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") output.M40[1,1:2]=c(corr$estimate,corr$p.value) output.M40[1,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) lines(as.numeric(names(mean.M40.cementation.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.cementation.500.y[N.M40.stained.500.y>CUTOFF], col="black") segments(x0=as.numeric(names(mean.M40.cementation.500.y))[N.M40.stained.500.y>CUTOFF]+50, x1=as.numeric(names(mean.M40.cementation.500.y))[N.M40.stained.500.y>CUTOFF]+50, y0 = LCI.M40.cementation.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.cementation.500.y[N.M40.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M40.discoloration.500.y))[N.M40.stained.500.y>CUTOFF]+100,mean.M40.discoloration.500.y[N.M40.stained.500.y>CUTOFF]) segments(x0=as.numeric(names(mean.M40.discoloration.500.y))[N.M40.stained.500.y>CUTOFF]+100, x1=as.numeric(names(mean.M40.discoloration.500.y))[N.M40.stained.500.y>CUTOFF]+100, y0 = LCI.M40.discoloration.500.y[N.M40.stained.500.y>CUTOFF], y1 = UCI.M40.discoloration.500.y[N.M40.stained.500.y>CUTOFF], lty=1) lines(as.numeric(names(mean.M40.discoloration.500.y))[N.M40.stained.500.y>CUTOFF]+100,mean.M40.discoloration.500.y[N.M40.stained.500.y>CUTOFF], col="black") points(as.numeric(names(mean.M40.discoloration.500.y))[N.M40.stained.500.y>CUTOFF]+100,mean.M40.discoloration.500.y[N.M40.stained.500.y>CUTOFF], pch=21, bg="white", cex=1.4) points(as.numeric(names(mean.M40.stained.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.stained.500.y[N.M40.stained.500.y>CUTOFF], pch=16, cex=1.4) points(as.numeric(names(mean.M40.cementation.500.y))[N.M40.stained.500.y>CUTOFF]+50,mean.M40.cementation.500.y[N.M40.stained.500.y>CUTOFF], pch=21, bg="gray71", cex=1.4) out=glm(mean.M40.cementation.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.cementation.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M40.cementation.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.cementation.500.y[N.M40.stained.500.y>CUTOFF]) output.M40[2,1:2]=c(corr$estimate,corr$p.value) output.M40[2,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) out=glm(mean.M40.discoloration.500.y[N.M40.stained.500.y>CUTOFF]~as.numeric(names(mean.M40.discoloration.500.y))[N.M40.stained.500.y>CUTOFF], family="quasibinomial") corr=cor.test(as.numeric(names(mean.M40.discoloration.500.y))[N.M40.stained.500.y>CUTOFF],mean.M40.discoloration.500.y[N.M40.stained.500.y>CUTOFF]) output.M40[9,1:2]=c(corr$estimate,corr$p.value) output.M40[9,3:5]=c(coef(summary(out))[2,c(1,4)],rsq(out)) output.M44[,c(1,3,5)]=c(round(output.M44[,1], digits=2),round(output.M44[,3], digits=5),round(output.M44[,5], digits=2)) output.M40[,c(1,3,5)]=c(round(output.M40[,1], digits=2),round(output.M40[,3], digits=5),round(output.M40[,5], digits=2)) output.M44[,c(2,4)]=c(round(output.M44[,2], digits=4),round(output.M44[,4], digits=4)) output.M40[,c(2,4)]=c(round(output.M40[,2], digits=4),round(output.M40[,4], digits=4)) output.M44 output.M40 temp=rbind(output.M44,output.M40) write.table(temp, file="Aggregate output.txt") ############################################# #END OF FIGURE WITH AGGREGATED COHORTS ############################################# ############################################# #SUM OF ALTERATIONS ############################################# plot(as.numeric(names(mean.M44.sum.500.y)),mean.M44.sum.500.y, cex=1.4, pch=16, ylim=c(0,8), xlim=c(0,5000), frame=F, xlab="Postmortem age (y)", ylab="Sum of altered states") lines(as.numeric(names(mean.M44.sum.500.y)),mean.M44.sum.500.y) segments(x0=as.numeric(names(mean.M44.sum.500.y)), x1=as.numeric(names(mean.M44.sum.500.y)), y0 = LCI.M44.sum.500.y, y1 = UCI.M44.sum.500.y, lty=1) plot(as.numeric(names(mean.M40.sum.500.y)),mean.M40.sum.500.y, cex=1.4, ylim=c(0,8), pch=16, xlim=c(0,5000), frame=F, xlab="Postmortem age (y)", ylab="") lines(as.numeric(names(mean.M40.sum.500.y)),mean.M40.sum.500.y) segments(x0=as.numeric(names(mean.M40.sum.500.y)), x1=as.numeric(names(mean.M40.sum.500.y)), y0 = LCI.M40.sum.500.y, y1 = UCI.M40.sum.500.y, lty=1) ########################################################################################################### #ESTIMATE MEANS OF ALTERATION PER ASSEMBLAGE BASED ON DATED SHELLS TO COMPARE WITH IQR (TIME AVERAGING) ########################################################################################################### Timoclea.M44.ext.sum=Timoclea.M44.unworn.encrustation+Timoclea.M44.stained+Timoclea.M44.ornamentation+Timoclea.M44.bioerosion+Timoclea.M44.worn.encrustation+Timoclea.M44.discoloration+Timoclea.M44.cementation+Timoclea.M44.dissolution Timoclea.M40.ext.sum=Timoclea.M40.unworn.encrustation+Timoclea.M40.stained+Timoclea.M40.ornamentation+Timoclea.M40.bioerosion+Timoclea.M40.worn.encrustation+Timoclea.M40.discoloration+Timoclea.M40.cementation+Timoclea.M40.dissolution M44.diss.dated=tapply(Timoclea.M44.dissolution[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M44.dissolution[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], bootci) M44.diss.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M44.diss.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M44.enc.dated=tapply(Timoclea.M44.total.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M44.total.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], bootci) M44.enc.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M44.enc.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M44.ext.enc.dated=tapply(Timoclea.M44.ext.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) M44.int.enc.dated=tapply(Timoclea.M44.int.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) M44.orn.dated=tapply(Timoclea.M44.ornamentation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M44.ornamentation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], bootci) M44.orn.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M44.orn.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M44.stain.dated=tapply(Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], bootci) M44.stain.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M44.stain.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M44.disc.dated=tapply(Timoclea.M44.discoloration[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) M44.cement.dated=tapply(Timoclea.M44.cementation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M44.cementation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], bootci) M44.cement.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M44.cement.dated.UCI=unlist(lapply(temp, quantile, 0.975)) Timoclea.M44.st.cem=ifelse(Timoclea.M44.stained | Timoclea.M44.cementation >0,1,0) M44.stain.cement.dated=tapply(Timoclea.M44.st.cem[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) M44.worn.enc.dated=tapply(Timoclea.M44.worn.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) M44.unworn.enc.dated=tapply(Timoclea.M44.unworn.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) M44.bio.dated=tapply(Timoclea.M44.bioerosion[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M44.bioerosion[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], bootci) M44.bio.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M44.bio.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M44.ext.sum.dated=tapply(Timoclea.M44.ext.sum[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M44.ext.sum[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.5cm[!is.na(Timoclea.M44.ages.for.clock)], bootci) M44.ext.sum.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M44.ext.sum.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M40.diss.dated=tapply(Timoclea.M40.dissolution[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M40.dissolution[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], bootci) M40.diss.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M40.diss.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M40.enc.dated=tapply(Timoclea.M40.total.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M40.total.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], bootci) M40.enc.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M40.enc.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M40.ext.enc.dated=tapply(Timoclea.M40.ext.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) M40.int.enc.dated=tapply(Timoclea.M40.int.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) M40.orn.dated=tapply(Timoclea.M40.ornamentation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M40.ornamentation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], bootci) M40.orn.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M40.orn.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M40.stain.dated=tapply(Timoclea.M40.stained[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M40.stained[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], bootci) M40.stain.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M40.stain.dated.UCI=unlist(lapply(temp, quantile, 0.975)) Timoclea.M40.st.cem=ifelse(Timoclea.M40.stained | Timoclea.M40.cementation >0,1,0) M40.stain.cement.dated=tapply(Timoclea.M40.st.cem[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) M40.disc.dated=tapply(Timoclea.M40.discoloration[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) M40.cement.dated=tapply(Timoclea.M40.cementation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M40.cementation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], bootci) M40.cement.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M40.cement.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M40.worn.enc.dated=tapply(Timoclea.M40.worn.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) M40.unworn.enc.dated=tapply(Timoclea.M40.unworn.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) M40.bio.dated=tapply(Timoclea.M40.bioerosion[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M40.bioerosion[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], bootci) M40.bio.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M40.bio.dated.UCI=unlist(lapply(temp, quantile, 0.975)) M40.ext.sum.dated=tapply(Timoclea.M40.ext.sum[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) temp=tapply(Timoclea.M40.ext.sum[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.5cm[!is.na(Timoclea.M40.ages.for.clock)], bootci) M40.ext.sum.dated.LCI=unlist(lapply(temp, quantile, 0.025)) M40.ext.sum.dated.UCI=unlist(lapply(temp, quantile, 0.975)) ############################################################ #TIMOCLEA TIME AVERAGING (IQR IN YEARS) ############################################################ #M44 IQR M44.raw.IQR=c(280,644,1017,1787,1418,2966,1894,1248,1188,1607,1542,1902,2117) M44.raw.IQR.LCI=c(6,96,457,305,197,996,1038,755,567,1116,705,1294,1386) M44.raw.IQR.UCI=c(1298,454,3001,2036,2409,1460,4907,5945,7033,3603) M44.corr.IQR=c(239,473,692,1409,1020,2258,1076,NA,NA,NA,NA,NA,NA) M44.corr.IQR.LCI=c(201,410,651,1338,925,2148,993,NA,NA,NA,NA,NA,NA) M44.corr.IQR.UCI=c(263,526,729,1474,1103,2358,1155,NA,NA,NA,NA,NA,NA) M44.group=c("SML","SML","SML","SML","SML","HST","HST","HST","HST","MFZ","MFZ","TST","TST") #M40 IQR M40.raw.IQR=c(45,415,296,655,638,547,2072,4881,5527,1856) M40.raw.IQR.LCI=c(3,0,5,154,0,55,964,1412,2796,789) M40.raw.IQR.UCI=c(1298,454,3001,2036,2409,1460,4907,5945,7033,3603) M40.corr.IQR=c(NA,364,140,433,171,154,1116,3632,4152,NA) M40.corr.IQR.LCI=c(NA,323,28,356,10,10,964,3508,3986,NA) M40.corr.IQR.UCI=c(NA,395,226,494,397,295,1261,3756,4316,NA) M40.group=c("SML","SML","SML","SML","SML","HST","HST","MFZ","MFZ","TST") ################################################################################ #MANUSCRIPT FIGURE 11 #BIVARIATE PLOTS FOR PER-INCREMENT PROPORTIONS IN DATED SPECIMENS ################################################################################ #cairo_pdf(file='plot.pdf', width=7, height=7) par(mfrow=c(4,4)) par(mar=c(3,3,1,0)) ord=cbind(c(M44.enc.dated,M40.enc.dated), c(M44.bio.dated,M40.bio.dated)) plot(c(M44.enc.dated,M40.enc.dated), c(M44.bio.dated,M40.bio.dated), pch=16, log="", xlab="", ylab="", frame=F, xlim=c(0,1), ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Encrustation", side = 1, line = 2.5, cex=0.7) mtext(text="Bioerosion", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.enc.dated,M40.enc.dated), x1=c(M44.enc.dated,M40.enc.dated), y0=c(M44.bio.dated.LCI,M40.bio.dated.LCI), y1=c(M44.bio.dated.UCI,M40.bio.dated.UCI), lty=1, col="gray") segments(x0=c(M44.enc.dated.LCI,M40.enc.dated.LCI), x1=c(M44.enc.dated.UCI,M40.enc.dated.UCI), y0=c(M44.bio.dated,M40.bio.dated), y1=c(M44.bio.dated,M40.bio.dated), lty=1, col="gray") points(c(M44.enc.dated,M40.enc.dated), c(M44.bio.dated,M40.bio.dated), pch=16) ord=cbind(c(M44.orn.dated,M40.orn.dated), c(M44.bio.dated,M40.bio.dated)) plot(c(M44.orn.dated,M40.orn.dated), c(M44.bio.dated,M40.bio.dated), pch=16, log="", xlab="", ylab="", frame=F, xlim=c(0,1), ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Ornamentation loss", side = 1, line = 2.5, cex=0.7) mtext(text="Bioerosion", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.orn.dated,M40.orn.dated), x1=c(M44.orn.dated,M40.orn.dated), y0=c(M44.bio.dated.LCI,M40.bio.dated.LCI), y1=c(M44.bio.dated.UCI,M40.bio.dated.UCI), lty=1, col="gray") segments(x0=c(M44.orn.dated.LCI,M40.orn.dated.LCI), x1=c(M44.orn.dated.UCI,M40.orn.dated.UCI), y0=c(M44.bio.dated,M40.bio.dated), y1=c(M44.bio.dated,M40.bio.dated), lty=1, col="gray") points(c(M44.orn.dated,M40.orn.dated), c(M44.bio.dated,M40.bio.dated), pch=16) ord=cbind(c(M44.stain.dated,M40.stain.dated), c(M44.bio.dated,M40.bio.dated)) plot(c(M44.stain.dated,M40.stain.dated), c(M44.bio.dated,M40.bio.dated), pch=16, log="", xlab="", ylab="", frame=F, xlim=c(0,1), ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Staining", side = 1, line = 2.5, cex=0.7) mtext(text="Bioerosion", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.stain.dated,M40.stain.dated), x1=c(M44.stain.dated,M40.stain.dated), y0=c(M44.bio.dated.LCI,M40.bio.dated.LCI), y1=c(M44.bio.dated.UCI,M40.bio.dated.UCI), lty=1, col="gray") segments(x0=c(M44.stain.dated.LCI,M40.stain.dated.LCI), x1=c(M44.stain.dated.UCI,M40.stain.dated.UCI), y0=c(M44.bio.dated,M40.bio.dated), y1=c(M44.bio.dated,M40.bio.dated), lty=1, col="gray") points(c(M44.stain.dated,M40.stain.dated), c(M44.bio.dated,M40.bio.dated), pch=16) ord=cbind(c(M44.cement.dated,M40.cement.dated), c(M44.bio.dated,M40.bio.dated)) plot(c(M44.cement.dated,M40.cement.dated), c(M44.bio.dated,M40.bio.dated), pch=16, log="", xlab="", ylab="", frame=F, xlim=c(0,1), ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Micrite envelopes", side = 1, line = 2.5, cex=0.7) mtext(text="Bioerosion", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.cement.dated,M40.cement.dated), x1=c(M44.cement.dated,M40.cement.dated), y0=c(M44.bio.dated.LCI,M40.bio.dated.LCI), y1=c(M44.bio.dated.UCI,M40.bio.dated.UCI), lty=1, col="gray") segments(x0=c(M44.cement.dated.LCI,M40.cement.dated.LCI), x1=c(M44.cement.dated.UCI,M40.cement.dated.UCI), y0=c(M44.bio.dated,M40.bio.dated), y1=c(M44.bio.dated,M40.bio.dated), lty=1, col="gray") points(c(M44.cement.dated,M40.cement.dated), c(M44.bio.dated,M40.bio.dated), pch=16) ord=cbind(c(M44.enc.dated,M40.enc.dated), c(M44.stain.dated,M40.stain.dated)) plot(c(M44.enc.dated,M40.enc.dated), c(M44.stain.dated,M40.stain.dated), pch=16, log="", xlab="", ylab="", frame=F, xlim=c(0,1), ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Encrustation", side = 1, line = 2.5, cex=0.7) mtext(text="Staining", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.enc.dated,M40.enc.dated), x1=c(M44.enc.dated,M40.enc.dated), y0=c(M44.stain.dated.LCI,M40.stain.dated.LCI), y1=c(M44.stain.dated.UCI,M40.stain.dated.UCI), lty=1, col="gray") segments(x0=c(M44.enc.dated.LCI,M40.enc.dated.LCI), x1=c(M44.enc.dated.UCI,M40.enc.dated.UCI), y0=c(M44.stain.dated,M40.stain.dated), y1=c(M44.stain.dated,M40.stain.dated), lty=1, col="gray") points(c(M44.enc.dated,M40.enc.dated), c(M44.stain.dated,M40.stain.dated), pch=16) ord=cbind(c(M44.orn.dated,M40.orn.dated), c(M44.stain.dated,M40.stain.dated)) plot(c(M44.orn.dated,M40.orn.dated), c(M44.stain.dated,M40.stain.dated), pch=16, log="", xlab="", ylab="", frame=F, xlim=c(0,1), ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Ornamentation loss", side = 1, line = 2.5, cex=0.7) mtext(text="Staining", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.orn.dated,M40.orn.dated), x1=c(M44.orn.dated,M40.orn.dated), y0=c(M44.stain.dated.LCI,M40.stain.dated.LCI), y1=c(M44.stain.dated.UCI,M40.stain.dated.UCI), lty=1, col="gray") segments(x0=c(M44.orn.dated.LCI,M40.orn.dated.LCI), x1=c(M44.orn.dated.UCI,M40.orn.dated.UCI), y0=c(M44.stain.dated,M40.stain.dated), y1=c(M44.stain.dated,M40.stain.dated), lty=1, col="gray") points(c(M44.orn.dated,M40.orn.dated), c(M44.stain.dated,M40.stain.dated), pch=16) ord=cbind(c(M44.orn.dated,M40.orn.dated), c(M44.cement.dated,M40.cement.dated)) plot(c(M44.orn.dated,M40.orn.dated), c(M44.cement.dated,M40.cement.dated), pch=16, log="", xlab="", ylab="", frame=F, xlim=c(0,1), ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Ornamentation loss", side = 1, line = 2.5, cex=0.7) mtext(text="Micrite envelopes", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.orn.dated,M40.orn.dated), x1=c(M44.orn.dated,M40.orn.dated), y0=c(M44.cement.dated.LCI,M40.cement.dated.LCI), y1=c(M44.cement.dated.UCI,M40.cement.dated.UCI), lty=1, col="gray") segments(x0=c(M44.orn.dated.LCI,M40.orn.dated.LCI), x1=c(M44.orn.dated.UCI,M40.orn.dated.UCI), y0=c(M44.cement.dated,M40.cement.dated), y1=c(M44.cement.dated,M40.cement.dated), lty=1, col="gray") points(c(M44.orn.dated,M40.orn.dated), c(M44.cement.dated,M40.cement.dated), pch=16) ord=cbind(c(M44.cement.dated,M40.cement.dated), c(M44.stain.dated,M40.stain.dated)) plot(c(M44.cement.dated,M40.cement.dated), c(M44.stain.dated,M40.stain.dated), xlim=c(0,1), pch=16, log="", xlab="", ylab="", frame=F, ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Micrite envelopes", side = 1, line = 2.5, cex=0.7) mtext(text="Staining", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.cement.dated,M40.cement.dated), x1=c(M44.cement.dated,M40.cement.dated), y0=c(M44.stain.dated.LCI,M40.stain.dated.LCI), y1=c(M44.stain.dated.UCI,M40.stain.dated.UCI), lty=1, col="gray") segments(x0=c(M44.cement.dated.LCI,M40.cement.dated.LCI), x1=c(M44.cement.dated.UCI,M40.cement.dated.UCI), y0=c(M44.stain.dated,M40.stain.dated), y1=c(M44.stain.dated,M40.stain.dated), lty=1, col="gray") points(c(M44.cement.dated,M40.cement.dated), c(M44.stain.dated,M40.stain.dated), pch=16) ord=cbind(c(M44.raw.IQR,M40.raw.IQR), c(M44.bio.dated,M40.bio.dated)) plot(c(M44.raw.IQR,M40.raw.IQR), c(M44.bio.dated,M40.bio.dated), pch=16, log="", xlim=c(10,7000), xlab="", ylab="", frame=F, ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Raw time averaging (y)", side = 1, line = 2.5, cex=0.7) mtext(text="Bioerosion", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.raw.IQR,M40.raw.IQR), x1=c(M44.raw.IQR,M40.raw.IQR), y0=c(M44.bio.dated.LCI,M40.bio.dated.LCI), y1=c(M44.bio.dated.UCI,M40.bio.dated.UCI), lty=1, col="gray") segments(x0=c(M44.raw.IQR.LCI,M40.raw.IQR.LCI), x1=c(M44.raw.IQR.UCI,M40.raw.IQR.UCI), y0=c(M44.bio.dated,M40.bio.dated), y1=c(M44.bio.dated,M40.bio.dated), lty=1, col="gray") ordihull(ord, group=c(M44.group,M40.group), display = "sites", draw = "polygon", col=c("gray81","gray51","white","gray21")) points(c(M44.raw.IQR,M40.raw.IQR), c(M44.bio.dated,M40.bio.dated), pch=16) points(c(M44.raw.IQR), c(M44.bio.dated), pch=16, cex=1.2) points(c(M40.raw.IQR), c(M40.bio.dated), pch=21, bg="white", cex=1.2) ord=cbind(c(M44.raw.IQR,M40.raw.IQR), c(M44.enc.dated,M40.enc.dated)) plot(c(M44.raw.IQR,M40.raw.IQR), c(M44.enc.dated,M40.enc.dated), pch=16, log="", xlim=c(10,7000), xlab="", ylab="", frame=F, ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Raw time averaging (y)", side = 1, line = 2.5, cex=0.7) mtext(text="Encrustation", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.raw.IQR,M40.raw.IQR), x1=c(M44.raw.IQR,M40.raw.IQR), y0=c(M44.enc.dated.LCI,M40.enc.dated.LCI), y1=c(M44.enc.dated.UCI,M40.enc.dated.UCI), lty=1, col="gray") segments(x0=c(M44.raw.IQR.LCI,M40.raw.IQR.LCI), x1=c(M44.raw.IQR.UCI,M40.raw.IQR.UCI), y0=c(M44.enc.dated,M40.enc.dated), y1=c(M44.enc.dated,M40.enc.dated), lty=1, col="gray") ordihull(ord, group=c(M44.group,M40.group), display = "sites", draw = "polygon", col=c("gray81","gray51","white","gray21")) points(c(M44.raw.IQR,M40.raw.IQR), c(M44.enc.dated,M40.enc.dated), pch=16) points(c(M44.raw.IQR), c(M44.enc.dated), pch=16, cex=1.2) points(c(M40.raw.IQR), c(M40.enc.dated), pch=21, bg="white", cex=1.2) ord=cbind(c(M44.raw.IQR,M40.raw.IQR), c(M44.orn.dated,M40.orn.dated)) plot(c(M44.raw.IQR,M40.raw.IQR), c(M44.orn.dated,M40.orn.dated), pch=16, log="", xlim=c(10,7000), xlab="", ylab="", frame=F, ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Raw time averaging (y)", side = 1, line = 2.5, cex=0.7) mtext(text="Ornamentation loss", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.raw.IQR,M40.raw.IQR), x1=c(M44.raw.IQR,M40.raw.IQR), y0=c(M44.orn.dated.LCI,M40.orn.dated.LCI), y1=c(M44.orn.dated.UCI,M40.orn.dated.UCI), lty=1, col="gray") segments(x0=c(M44.raw.IQR.LCI,M40.raw.IQR.LCI), x1=c(M44.raw.IQR.UCI,M40.raw.IQR.UCI), y0=c(M44.orn.dated,M40.orn.dated), y1=c(M44.orn.dated,M40.orn.dated), lty=1, col="gray") ordihull(ord, group=c(M44.group,M40.group), display = "sites", draw = "polygon", col=c("gray81","gray51","white","gray21")) points(c(M44.raw.IQR,M40.raw.IQR), c(M44.orn.dated,M40.orn.dated), pch=16) points(c(M44.raw.IQR), c(M44.orn.dated), pch=16, cex=1.2) points(c(M40.raw.IQR), c(M40.orn.dated), pch=21, bg="white", cex=1.2) ord=cbind(c(M44.raw.IQR,M40.raw.IQR), c(M44.stain.dated,M40.stain.dated)) plot(c(M44.raw.IQR,M40.raw.IQR), c(M44.stain.dated,M40.stain.dated), pch=16, log="", xlim=c(10,7000), xlab="", ylab="", frame=F, ylim=c(0,1), cex.axis=0.9, cex.lab=0.9) mtext(text="Raw time averaging (y)", side = 1, line = 2.5, cex=0.7) mtext(text="Staining", side = 2, line = 2.25, cex=0.7) segments(x0=c(M44.raw.IQR,M40.raw.IQR), x1=c(M44.raw.IQR,M40.raw.IQR), y0=c(M44.stain.dated.LCI,M40.stain.dated.LCI), y1=c(M44.stain.dated.UCI,M40.stain.dated.UCI), lty=1, col="gray") segments(x0=c(M44.raw.IQR.LCI,M40.raw.IQR.LCI), x1=c(M44.raw.IQR.UCI,M40.raw.IQR.UCI), y0=c(M44.stain.dated,M40.stain.dated), y1=c(M44.stain.dated,M40.stain.dated), lty=1, col="gray") ordihull(ord, group=c(M44.group,M40.group), display = "sites", draw = "polygon", col=c("gray81","gray51","white","gray21")) points(c(M44.raw.IQR,M40.raw.IQR), c(M44.stain.dated,M40.stain.dated), pch=16) points(c(M44.raw.IQR), c(M44.stain.dated), pch=16, cex=1.2) points(c(M40.raw.IQR), c(M40.stain.dated), pch=21, bg="white", cex=1.2) legend(x="topleft", pch=22, pt.bg=c("white","gray81","gray51","gray21"), legend=c("ML","HST","MFZ","TST"), cex=1.2, bty="n", ncol=2) #dev.off() ######################################################################################################### #END OF ASSEMBLAGE-LEVEL TAPHONOMIC CLOCK FIGURE ######################################################################################################### ############################################################################################### #CORRELATIONS FOR ASSEMBLAGE-LEVEL VARIABLES ############################################################################################### out=cor.test(c(M44.bio.dated,M40.bio.dated), c(M44.stain.dated,M40.stain.dated), method="s") out=cor.test(c(M44.bio.dated,M40.bio.dated), c(M44.enc.dated,M40.enc.dated), method="s") out=cor.test(c(M44.bio.dated,M40.bio.dated), c(M44.diss.dated,M40.diss.dated), method="s") out=cor.test(c(M44.bio.dated,M40.bio.dated), c(M44.cement.dated,M40.cement.dated), method="s") out=cor.test(c(M44.bio.dated,M40.bio.dated), c(M44.orn.dated,M40.orn.dated), method="s") out=cor.test(c(M44.cement.dated,M40.cement.dated), c(M44.stain.dated,M40.stain.dated), method="s") out=cor.test(c(M44.cement.dated,M40.cement.dated), c(M44.diss.dated,M40.diss.dated), method="s") out=cor.test(c(M44.cement.dated,M40.cement.dated), c(M44.enc.dated,M40.enc.dated), method="s") out=cor.test(c(M44.stain.dated,M40.stain.dated), c(M44.orn.dated,M40.orn.dated), method="s") out=cor.test(c(M44.cement.dated,M40.cement.dated), c(M44.orn.dated,M40.orn.dated), method="s") out=cor.test(c(M44.stain.dated,M40.stain.dated),c(M44.enc.dated,M40.enc.dated), method="s") ###################################################### #FIGURE BOXPLOT WITH SYSTEM TRACT VS AVERAGING ###################################################### boxplot(split(c(M44.raw.IQR,M40.raw.IQR), c(M44.group,M40.group))[c(3,1,2,4)], col="gray", frame=F, ylab="Interquartile age range (y)") ###################################################### #TAPHONOMIC CLOCK FOR SUM OF ALTERATIONS# ###################################################### #cairo_pdf(file='plot.pdf', width=7, height=7) par(mfrow=c(2,2)) par(mar=c(4,4,2,1)) ord=cbind(c(M44.raw.IQR,M40.raw.IQR), c(M44.ext.sum.dated,M40.ext.sum.dated)) plot(c(M44.raw.IQR,M40.raw.IQR), c(M44.ext.sum.dated,M40.ext.sum.dated), pch=16, log="", xlim=c(10,7000), xlab="Raw time averaging (y)", ylab="Sum of alteration signatures", frame=F, ylim=c(0,6)) segments(x0=c(M44.raw.IQR,M40.raw.IQR), x1=c(M44.raw.IQR,M40.raw.IQR), y0=c(M44.ext.sum.dated.LCI,M40.ext.sum.dated.LCI), y1=c(M44.ext.sum.dated.UCI,M40.ext.sum.dated.UCI), lty=1, col="gray") segments(x0=c(M44.raw.IQR.LCI,M40.raw.IQR.LCI), x1=c(M44.raw.IQR.UCI,M40.raw.IQR.UCI), y0=c(M44.ext.sum.dated,M40.ext.sum.dated), y1=c(M44.ext.sum.dated,M40.ext.sum.dated), lty=1, col="gray") ordihull(ord, group=c(M44.group,M40.group), display = "sites", draw = "polygon", col=c("gray81","gray51","white","gray21")) points(c(M44.raw.IQR,M40.raw.IQR), c(M44.ext.sum.dated,M40.ext.sum.dated), pch=16) points(c(M44.raw.IQR), c(M44.ext.sum.dated), pch=16, cex=1.2) points(c(M40.raw.IQR), c(M40.ext.sum.dated), pch=21, bg="white", cex=1.2) ord=cbind(c(M44.corr.IQR,M40.corr.IQR), c(M44.ext.sum.dated,M40.ext.sum.dated)) plot(c(M44.corr.IQR,M40.corr.IQR), c(M44.ext.sum.dated,M40.ext.sum.dated), pch=16, log="", xlim=c(10,7000), xlab="Corrected time averaging (y)", ylab="Sum of alteration signatures", frame=F, ylim=c(0,6)) segments(x0=c(M44.corr.IQR,M40.corr.IQR), x1=c(M44.corr.IQR,M40.corr.IQR), y0=c(M44.ext.sum.dated.LCI,M40.ext.sum.dated.LCI), y1=c(M44.ext.sum.dated.UCI,M40.ext.sum.dated.UCI), lty=1, col="gray") segments(x0=c(M44.corr.IQR.LCI,M40.corr.IQR.LCI), x1=c(M44.corr.IQR.UCI,M40.corr.IQR.UCI), y0=c(M44.ext.sum.dated,M40.ext.sum.dated), y1=c(M44.ext.sum.dated,M40.ext.sum.dated), lty=1, col="gray") ordihull(ord, group=c(M44.group,M40.group), display = "sites", draw = "polygon",col=c("gray81","gray51","white","gray21")) points(c(M44.corr.IQR,M40.corr.IQR), c(M44.ext.sum.dated,M40.ext.sum.dated), pch=16) points(c(M44.corr.IQR), c(M44.ext.sum.dated), pch=16, cex=1.2) points(c(M40.corr.IQR), c(M40.ext.sum.dated), pch=21, bg="white", cex=1.2) legend(x="bottomright", pch=21, pt.bg=c("black","white"),legend=c("M44","M40"), bty="n", cex=1.2) #dev.off() ########################################################################################################### #OUTPUT FOR ASSEMBLAGE-LEVEL CLOCK ########################################################################################################### par(mfrow=c(2,4)) par(mar=c(4,4,2,1)) METHOD="s" options(scipen = 999) output.ass.clock=array(NA, dim=c(10,2)) rownames(output.ass.clock)=c("Total encrustation","Fresh encrusters","Worn encrusters","Bioerosion","Dissolution", "Ornamentation loss","Staining","Cementation","Discoloration","Sum") colnames(output.ass.clock)=c("Spearman r", "p") out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.enc.dated,M40.enc.dated), method=METHOD) output.ass.clock[1,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.unworn.enc.dated,M40.unworn.enc.dated), method=METHOD) output.ass.clock[2,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.worn.enc.dated,M40.worn.enc.dated), method=METHOD) output.ass.clock[3,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.bio.dated,M40.bio.dated), method=METHOD) output.ass.clock[4,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.diss.dated,M40.diss.dated), method=METHOD) output.ass.clock[5,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.orn.dated,M40.orn.dated), method=METHOD) output.ass.clock[6,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.stain.dated,M40.stain.dated), method=METHOD) output.ass.clock[7,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.cement.dated,M40.cement.dated), method=METHOD) output.ass.clock[8,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.disc.dated,M40.disc.dated), method=METHOD) output.ass.clock[9,1:2]=c(out$estimate,out$p.value) out=cor.test(c(M44.raw.IQR,M40.raw.IQR), c(M44.ext.sum.dated,M40.ext.sum.dated), method=METHOD) output.ass.clock[10,1:2]=c(out$estimate,out$p.value) output.ass.clock=round(output.ass.clock, digits=3) output.ass.clock write.table(output.ass.clock, file="Assemblage clock.txt") ######################################################################################################################## #TAPHONOMIC CLOCK WITH DEPTH AT M44 ######################################################################################################################## par(mfrow=c(2,3)) par(mar=c(4,2,2,1)) boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.discoloration==0], Timoclea.M44.10cm[Timoclea.M44.discoloration==0]), at=-as.numeric(names(split(Timoclea.M44.discoloration[Timoclea.M44.discoloration==0], Timoclea.M44.10cm[Timoclea.M44.discoloration==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Discoloration") boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.discoloration==1], Timoclea.M44.10cm[Timoclea.M44.discoloration==1]), at=-as.numeric(names(split(Timoclea.M44.discoloration[Timoclea.M44.discoloration==1], Timoclea.M44.10cm[Timoclea.M44.discoloration==1])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.ornamentation==0], Timoclea.M44.10cm[Timoclea.M44.ornamentation==0]), at=-as.numeric(names(split(Timoclea.M44.ornamentation[Timoclea.M44.ornamentation==0], Timoclea.M44.10cm[Timoclea.M44.ornamentation==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Ornamentation loss") boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.ornamentation>0], Timoclea.M44.10cm[Timoclea.M44.ornamentation>0]), at=-as.numeric(names(split(Timoclea.M44.ornamentation[Timoclea.M44.ornamentation>0], Timoclea.M44.10cm[Timoclea.M44.ornamentation>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.bioerosion==0], Timoclea.M44.10cm[Timoclea.M44.bioerosion==0]), at=-as.numeric(names(split(Timoclea.M44.bioerosion[Timoclea.M44.bioerosion==0], Timoclea.M44.10cm[Timoclea.M44.bioerosion==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Bioerosion") boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.bioerosion>0], Timoclea.M44.10cm[Timoclea.M44.bioerosion>0]), at=-as.numeric(names(split(Timoclea.M44.bioerosion[Timoclea.M44.bioerosion>0], Timoclea.M44.10cm[Timoclea.M44.bioerosion>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.dissolution==0], Timoclea.M44.10cm[Timoclea.M44.dissolution==0]), at=-as.numeric(names(split(Timoclea.M44.dissolution[Timoclea.M44.dissolution==0], Timoclea.M44.10cm[Timoclea.M44.dissolution==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Dissolution") boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.dissolution>0], Timoclea.M44.10cm[Timoclea.M44.dissolution>0]), at=-as.numeric(names(split(Timoclea.M44.dissolution[Timoclea.M44.dissolution>0], Timoclea.M44.10cm[Timoclea.M44.dissolution>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.stained==0], Timoclea.M44.10cm[Timoclea.M44.stained==0]), at=-as.numeric(names(split(Timoclea.M44.bioerosion[Timoclea.M44.stained==0], Timoclea.M44.10cm[Timoclea.M44.stained==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Staining") boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.stained>0], Timoclea.M44.10cm[Timoclea.M44.stained>0]), at=-as.numeric(names(split(Timoclea.M44.bioerosion[Timoclea.M44.stained>0], Timoclea.M44.10cm[Timoclea.M44.stained>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.ext.encrustation==0], Timoclea.M44.10cm[Timoclea.M44.ext.encrustation==0]), at=-as.numeric(names(split(Timoclea.M44.bioerosion[Timoclea.M44.ext.encrustation==0], Timoclea.M44.10cm[Timoclea.M44.ext.encrustation==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Ext encrustation") boxplot(split(Timoclea.M44.ages.for.clock[Timoclea.M44.ext.encrustation>0], Timoclea.M44.10cm[Timoclea.M44.ext.encrustation>0]), at=-as.numeric(names(split(Timoclea.M44.bioerosion[Timoclea.M44.ext.encrustation>0], Timoclea.M44.10cm[Timoclea.M44.ext.encrustation>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) #M40 par(mfrow=c(2,3)) par(mar=c(4,2,2,1)) boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.discoloration==0], Timoclea.M40.10cm[Timoclea.M40.discoloration==0]), at=-as.numeric(names(split(Timoclea.M40.discoloration[Timoclea.M40.discoloration==0], Timoclea.M40.10cm[Timoclea.M40.discoloration==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Discoloration") boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.discoloration==1], Timoclea.M40.10cm[Timoclea.M40.discoloration==1]), at=-as.numeric(names(split(Timoclea.M40.discoloration[Timoclea.M40.discoloration==1], Timoclea.M40.10cm[Timoclea.M40.discoloration==1])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.ornamentation==0], Timoclea.M40.10cm[Timoclea.M40.ornamentation==0]), at=-as.numeric(names(split(Timoclea.M40.ornamentation[Timoclea.M40.ornamentation==0], Timoclea.M40.10cm[Timoclea.M40.ornamentation==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Ornamentation loss") boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.ornamentation>0], Timoclea.M40.10cm[Timoclea.M40.ornamentation>0]), at=-as.numeric(names(split(Timoclea.M40.ornamentation[Timoclea.M40.ornamentation>0], Timoclea.M40.10cm[Timoclea.M40.ornamentation>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.bioerosion==0], Timoclea.M40.10cm[Timoclea.M40.bioerosion==0]), at=-as.numeric(names(split(Timoclea.M40.bioerosion[Timoclea.M40.bioerosion==0], Timoclea.M40.10cm[Timoclea.M40.bioerosion==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Bioerosion") boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.bioerosion>0], Timoclea.M40.10cm[Timoclea.M40.bioerosion>0]), at=-as.numeric(names(split(Timoclea.M40.bioerosion[Timoclea.M40.bioerosion>0], Timoclea.M40.10cm[Timoclea.M40.bioerosion>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.dissolution==0], Timoclea.M40.10cm[Timoclea.M40.dissolution==0]), at=-as.numeric(names(split(Timoclea.M40.dissolution[Timoclea.M40.dissolution==0], Timoclea.M40.10cm[Timoclea.M40.dissolution==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Dissolution") boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.dissolution>0], Timoclea.M40.10cm[Timoclea.M40.dissolution>0]), at=-as.numeric(names(split(Timoclea.M40.dissolution[Timoclea.M40.dissolution>0], Timoclea.M40.10cm[Timoclea.M40.dissolution>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.stained==0], Timoclea.M40.10cm[Timoclea.M40.stained==0]), at=-as.numeric(names(split(Timoclea.M40.bioerosion[Timoclea.M40.stained==0], Timoclea.M40.10cm[Timoclea.M40.stained==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Staining") boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.stained>0], Timoclea.M40.10cm[Timoclea.M40.stained>0]), at=-as.numeric(names(split(Timoclea.M40.bioerosion[Timoclea.M40.stained>0], Timoclea.M40.10cm[Timoclea.M40.stained>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.ext.encrustation==0], Timoclea.M40.10cm[Timoclea.M40.ext.encrustation==0]), at=-as.numeric(names(split(Timoclea.M40.bioerosion[Timoclea.M40.ext.encrustation==0], Timoclea.M40.10cm[Timoclea.M40.ext.encrustation==0]))),range=0, col="white", boxwex=3, horizontal=T, frame=F, main="Ext encrustation") boxplot(split(Timoclea.M40.ages.for.clock[Timoclea.M40.ext.encrustation>0], Timoclea.M40.10cm[Timoclea.M40.ext.encrustation>0]), at=-as.numeric(names(split(Timoclea.M40.bioerosion[Timoclea.M40.ext.encrustation>0], Timoclea.M40.10cm[Timoclea.M40.ext.encrustation>0])))-3,range=0, col="gray", boxwex=3, add=T, horizontal=T, yaxt="n", xaxt="n", frame=F) ########################################################################################################### ########################################################################################################### #COMPREHENSIVE DATASET - SCORING ALL INDIVIDUALS OF CORBULA AND TIMOCLEA, SOME ANALYSES JUST FOR FIVE VARIABLES ########################################################################################################### ########################################################################################################### #here, taphonomic scoring is based on preservation of non-dated shells entered directly into vectors below - #these scores are then combined with scores in Table S1 as computed above ################ #BRIJUNI M44 ################ #2 cm and 5 cm depth Corbula at M44 - each increment is represented by a pair - the first number is the count of unaltered valves #the second number is the count of altered valves Corbula.Brijuni.M44.add.for.stained=c(c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,8),rep(1,0)),c(rep(0,5),rep(1,0)),c(rep(0,11),rep(1,1)),c(rep(0,13),rep(1,0)),c(rep(0,10),rep(1,0)),c(rep(0,13),rep(1,2)),c(rep(0,2),rep(1,0)),c(rep(0,8),rep(1,1)),c(rep(0,6),rep(1,1)),c(rep(0,18),rep(1,1)),c(rep(0,14),rep(1,4)),c(rep(0,28),rep(1,1)),c(rep(0,31),rep(1,4)),c(rep(0,19),rep(1,3)),c(rep(0,23),rep(1,7)),c(rep(0,23),rep(1,6)),c(rep(0,17),rep(1,8)),c(rep(0,18),rep(1,11)),c(rep(0,21),rep(1,15)),c(rep(0,34),rep(1,17)),c(rep(0,34),rep(1,18)), c(rep(0,20),rep(1,14)),c(rep(0,21),rep(1,9)),c(rep(0,12),rep(1,11)),c(rep(0,16),rep(1,12)),c(rep(0,8),rep(1,8)),c(rep(0,11),rep(1,11)),c(rep(0,14),rep(1,5)),c(rep(0,14),rep(1,15)),c(rep(0,12),rep(1,6)),c(rep(0,13),rep(1,18)),c(rep(0,7),rep(1,8)),c(rep(0,12),rep(1,12)),c(rep(0,4),rep(1,10)),c(rep(0,5),rep(1,6))) Corbula.Brijuni.M44.add.fine.depth.for.stained=c(rep(2,0),rep(4,2),rep(6,8),rep(8,5),rep(10,12),rep(12,13),rep(14,10),rep(16,15),rep(18,2),rep(20,9),rep(25,7),rep(30,19),rep(35,18),rep(40,29),rep(45,35),rep(50,22),rep(55,30),rep(60,29),rep(65,25),rep(70,29), rep(75,36),rep(80,51),rep(85,52), rep(90,34),rep(95,30),rep(100,23),rep(105,28),rep(110,16),rep(115,22),rep(120,19),rep(125,29), rep(130,18),rep(135,31),rep(140,15),rep(145,24),rep(150,14),rep(155,11)) Corbula.Brijuni.M44.add.coarse.depth.for.stained=ceiling(Corbula.Brijuni.M44.add.fine.depth.for.stained/10)*10 Corbula.Brijuni.M44.add.for.encrusted=c(c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,7),rep(1,1)),c(rep(0,5),rep(1,0)),c(rep(0,11),rep(1,1)),c(rep(0,13),rep(1,0)),c(rep(0,9),rep(1,1)),c(rep(0,12),rep(1,3)),c(rep(0,2),rep(1,0)),c(rep(0,8),rep(1,1)),c(rep(0,7),rep(1,0)),c(rep(0,16),rep(1,3)),c(rep(0,16),rep(1,2)),c(rep(0,27),rep(1,2)),c(rep(0,32),rep(1,3)),c(rep(0,19),rep(1,3)),c(rep(0,20),rep(1,6)),c(rep(0,21),rep(1,8)),c(rep(0,19),rep(1,6)),c(rep(0,21),rep(1,8)),c(rep(0,32),rep(1,6)),c(rep(0,44),rep(1,11)),c(rep(0,40),rep(1,9)),c(rep(0,28),rep(1,4)),c(rep(0,23),rep(1,7)),c(rep(0,14),rep(1,7)),c(rep(0,22),rep(1,7)),c(rep(0,10),rep(1,4)),c(rep(0,17),rep(1,6)),c(rep(0,14),rep(1,5)),c(rep(0,18),rep(1,11)),c(rep(0,13),rep(1,5)),c(rep(0,17),rep(1,14)),c(rep(0,10),rep(1,6)),c(rep(0,16),rep(1,8)),c(rep(0,4),rep(1,10)),c(rep(0,8),rep(1,3))) Corbula.Brijuni.M44.add.fine.depth.for.encrusted=c(rep(2,0),rep(4,2),rep(6,8),rep(8,5),rep(10,12),rep(12,13),rep(14,10),rep(16,15),rep(18,2),rep(20,9),rep(25,7),rep(30,19),rep(35,18),rep(40,29),rep(45,35),rep(50,22),rep(55,26),rep(60,29),rep(65,25),rep(70,29), rep(75,38),rep(80,55),rep(85,49), rep(90,32),rep(95,30),rep(100,21),rep(105,29),rep(110,14),rep(115,23),rep(120,19),rep(125,29), rep(130,18),rep(135,31),rep(140,16),rep(145,24),rep(150,14),rep(155,11)) Corbula.Brijuni.M44.add.coarse.depth.for.encrusted=ceiling(Corbula.Brijuni.M44.add.fine.depth.for.encrusted/10)*10 Corbula.Brijuni.M44.add.for.disornamented=c(c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,7),rep(1,1)),c(rep(0,5),rep(1,0)),c(rep(0,9),rep(1,3)),c(rep(0,11),rep(1,2)),c(rep(0,9),rep(1,1)),c(rep(0,10),rep(1,5)),c(rep(0,2),rep(1,0)),c(rep(0,7),rep(1,2)),c(rep(0,7),rep(1,0)),c(rep(0,16),rep(1,3)),c(rep(0,13),rep(1,5)),c(rep(0,26),rep(1,3)),c(rep(0,31),rep(1,4)),c(rep(0,16),rep(1,6)),c(rep(0,14),rep(1,12)),c(rep(0,17),rep(1,12)),c(rep(0,14),rep(1,11)),c(rep(0,17),rep(1,12)), c(rep(0,22),rep(1,16)),c(rep(0,38),rep(1,17)),c(rep(0,30),rep(1,19)),c(rep(0,18),rep(1,14)),c(rep(0,18),rep(1,12)),c(rep(0,12),rep(1,9)),c(rep(0,17),rep(1,12)),c(rep(0,5),rep(1,9)),c(rep(0,13),rep(1,10)),c(rep(0,11),rep(1,8)),c(rep(0,14),rep(1,15)),c(rep(0,9),rep(1,9)),c(rep(0,12),rep(1,19)),c(rep(0,9),rep(1,7)),c(rep(0,11),rep(1,13)),c(rep(0,3),rep(1,11)),c(rep(0,5),rep(1,6))) Corbula.Brijuni.M44.add.coarse.depth.for.disornamented=Corbula.Brijuni.M44.add.coarse.depth.for.encrusted #the first number is the count of unworn encrusters, the second is worn encruster Corbula.Brijuni.M44.add.for.enc.preserved=c(c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,3),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,3),rep(1,0)),c(rep(0,1),rep(1,1)),c(rep(0,1),rep(1,1)),c(rep(0,3),rep(1,0)),c(rep(0,2),rep(1,1)),c(rep(0,5),rep(1,1)),c(rep(0,3),rep(1,5)),c(rep(0,3),rep(1,3)),c(rep(0,6),rep(1,2)),c(rep(0,3),rep(1,3)),c(rep(0,7),rep(1,4)),c(rep(0,4),rep(1,5)),c(rep(0,1),rep(1,3)),c(rep(0,3),rep(1,4)),c(rep(0,4),rep(1,3)),c(rep(0,4),rep(1,3)),c(rep(0,1),rep(1,3)),c(rep(0,3),rep(1,3)),c(rep(0,3),rep(1,2)),c(rep(0,4),rep(1,7)),c(rep(0,1),rep(1,4)),c(rep(0,3),rep(1,11)),c(rep(0,1),rep(1,5)),c(rep(0,1),rep(1,7)),c(rep(0,3),rep(1,7)),c(rep(0,1),rep(1,2))) Corbula.Brijuni.M44.add.fine.depth.for.enc.preserved=c(rep(2,0),rep(4,0),rep(6,1),rep(8,0),rep(10,1),rep(12,0),rep(14,1),rep(16,3),rep(18,0),rep(20,1),rep(25,0),rep(30,3),rep(35,2),rep(40,2),rep(45,3),rep(50,3),rep(55,6),rep(60,8),rep(65,6),rep(70,8), rep(75,6),rep(80,11),rep(85,9),rep(90,4),rep(95,7),rep(100,7),rep(105,7),rep(110,4),rep(115,6),rep(120,5),rep(125,11), rep(130,5),rep(135,14),rep(140,6),rep(145,8),rep(150,10),rep(155,3)) Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved=ceiling(Corbula.Brijuni.M44.add.fine.depth.for.enc.preserved/10)*10 #the first number is the count of stained encrusters, the second is unstained encruster, both on stained shells Corbula.Brijuni.M44.add.for.stained.enc.preserved=c(c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,2),rep(1,1)),c(rep(0,1),rep(1,0)),c(rep(0,2),rep(1,1)),c(rep(0,4),rep(1,1)),c(rep(0,4),rep(1,0)),c(rep(0,3),rep(1,3)),c(rep(0,3),rep(1,3)),c(rep(0,6),rep(1,3)),c(rep(0,3),rep(1,3)),c(rep(0,4),rep(1,1)),c(rep(0,5),rep(1,1)),c(rep(0,1),rep(1,4)),c(rep(0,5),rep(1,1)),c(rep(0,3),rep(1,0)),c(rep(0,5),rep(1,2)),c(rep(0,1),rep(1,1)),c(rep(0,4),rep(1,4)),c(rep(0,2),rep(1,1)),c(rep(0,3),rep(1,4)),c(rep(0,1),rep(1,2)),c(rep(0,3),rep(1,6)),c(rep(0,3),rep(1,5)),c(rep(0,1),rep(1,1)),c(rep(0,0),rep(1,1))) Corbula.Brijuni.M44.add.fine.depth.for.stained.enc.preserved=c(rep(2,0),rep(4,0),rep(6,0),rep(8,0),rep(10,1),rep(12,0),rep(14,0),rep(16,2),rep(18,0),rep(20,1),rep(25,0),rep(30,1),rep(35,2),rep(40,1),rep(45,3), rep(50,1),rep(55,3),rep(60,5),rep(65,4),rep(70,6),rep(75,6),rep(80,9),rep(85,6),rep(90,5),rep(95,6),rep(100,5),rep(105,6),rep(110,3),rep(115,7),rep(120,2),rep(125,8), rep(130,3),rep(135,7),rep(140,3),rep(145,9),rep(150,8),rep(155,2),rep(160,1)) Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved=ceiling(Corbula.Brijuni.M44.add.fine.depth.for.stained.enc.preserved/10)*10 #Timoclea - the first number is the count of unworn, the second is worn encruster Timoclea.Brijuni.M44.add.for.enc.preserved=c(c(rep(0,0),rep(1,2)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,1)),c(rep(0,4),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,3)),c(rep(0,1),rep(1,1)),c(rep(0,7),rep(1,1)),c(rep(0,8),rep(1,2)),c(rep(0,7),rep(1,3)),c(rep(0,9),rep(1,4)),c(rep(0,5),rep(1,3)),c(rep(0,11),rep(1,2)),c(rep(0,9),rep(1,8)),c(rep(0,20),rep(1,8)),c(rep(0,17),rep(1,12)),c(rep(0,18),rep(1,12)),c(rep(0,27),rep(1,21)),c(rep(0,35),rep(1,22)),c(rep(0,29),rep(1,16)),c(rep(0,32),rep(1,29)),c(rep(0,34),rep(1,26)),c(rep(0,30),rep(1,26)),c(rep(0,33),rep(1,22)),c(rep(0,35),rep(1,27)),c(rep(0,20),rep(1,18)),c(rep(0,15),rep(1,16)),c(rep(0,32),rep(1,33)),c(rep(0,31),rep(1,24)),c(rep(0,38),rep(1,34)),c(rep(0,28),rep(1,31)),c(rep(0,31),rep(1,39)),c(rep(0,10),rep(1,18)),c(rep(0,14),rep(1,19))) Timoclea.Brijuni.M44.add.fine.depth.for.enc.preserved=c(rep(2,2),rep(4,0),rep(6,0),rep(8,1),rep(10,0),rep(12,3),rep(14,4),rep(16,1),rep(18,3),rep(20,2),rep(25,8),rep(30,10),rep(35,10),rep(40,13),rep(45,8),rep(50,13),rep(55,17),rep(60,28),rep(65,29),rep(70,30), rep(75,48),rep(80,57),rep(85,45),rep(90,61),rep(95,60),rep(100,56),rep(105,55),rep(110,62),rep(115,38),rep(120,31),rep(125,65), rep(130,55),rep(135,72),rep(140,59),rep(145,70),rep(150,28),rep(155,33)) Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved=ceiling(Timoclea.Brijuni.M44.add.fine.depth.for.enc.preserved/10)*10 #Timoclea - the first number is the count of stained encrusters, the second is unstained encruster, both on stained shells Timoclea.Brijuni.M44.add.for.stained.enc.preserved=c(c(rep(0,0),rep(1,2)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,1)),c(rep(0,2),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,6),rep(1,2)),c(rep(0,8),rep(1,2)),c(rep(0,8),rep(1,5)),c(rep(0,8),rep(1,3)),c(rep(0,18),rep(1,7)),c(rep(0,20),rep(1,9)),c(rep(0,12),rep(1,3)),c(rep(0,25),rep(1,16)),c(rep(0,27),rep(1,15)),c(rep(0,22),rep(1,13)),c(rep(0,24),rep(1,10)),c(rep(0,33),rep(1,15)),c(rep(0,11),rep(1,8)),c(rep(0,14),rep(1,7)),c(rep(0,16),rep(1,15)),c(rep(0,9),rep(1,11)),c(rep(0,15),rep(1,16)),c(rep(0,14),rep(1,13)),c(rep(0,15),rep(1,18)),c(rep(0,9),rep(1,10)),c(rep(0,6),rep(1,9)),c(rep(0,2),rep(1,3))) Timoclea.Brijuni.M44.add.fine.depth.for.stained.enc.preserved=c(rep(2,2),rep(4,0),rep(6,0),rep(8,1),rep(10,2),rep(12,0),rep(14,1),rep(16,1),rep(18,2),rep(20,1),rep(25,1),rep(30,0),rep(35,1),rep(40,1),rep(45,1),rep(50,0),rep(55,8),rep(60,10),rep(65,13),rep(70,11),rep(75,25),rep(80,29),rep(85,15),rep(90,41),rep(95,42),rep(100,35),rep(105,34),rep(110,48),rep(115,19),rep(120,21),rep(125,31), rep(130,20),rep(135,31),rep(140,27),rep(145,33),rep(150,19),rep(155,15),rep(160,5)) Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved=ceiling(Timoclea.Brijuni.M44.add.fine.depth.for.stained.enc.preserved/10)*10 ################ #BRIJUNI M40 ################ #counts of stained, encrusted, and disornamented valves in 2 cm and 5 cm increments, from top to bottom Corbula.Brijuni.M40.add.fine.depth.for.stained=c(rep(2,0),rep(4,1),rep(6,0),rep(8,0),rep(10,0),rep(12,4),rep(14,0),rep(16,7),rep(18,0),rep(20,9),rep(25,17),rep(35,12),rep(45,19),rep(55,44),rep(60,28),rep(65,76),rep(75,61),rep(80,26),rep(85,69),rep(90,59),rep(95,46),rep(100,45),rep(105,9),rep(110,44),rep(115,10)) Corbula.Brijuni.M40.add.coarse.depth.for.stained=ceiling(Corbula.Brijuni.M40.add.fine.depth.for.stained/10)*10 Corbula.Brijuni.M40.add.for.stained=c(c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,4),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,7),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,8),rep(1,1)),c(rep(0,13),rep(1,4)),c(rep(0,8),rep(1,4)),c(rep(0,9),rep(1,10)),c(rep(0,15),rep(1,29)),c(rep(0,11),rep(1,17)),c(rep(0,19),rep(1,57)),c(rep(0,11),rep(1,50)),c(rep(0,7),rep(1,19)),c(rep(0,24),rep(1,45)),c(rep(0,10),rep(1,49)),c(rep(0,12),rep(1,34)),c(rep(0,7),rep(1,38)),c(rep(0,1),rep(1,8)),c(rep(0,12),rep(1,32)),c(rep(0,0),rep(1,10))) Corbula.Brijuni.M40.add.fine.depth.for.encrusted=c(rep(2,0),rep(4,1),rep(6,0),rep(8,0),rep(10,0),rep(12,4),rep(14,0),rep(16,7),rep(18,0),rep(20,9),rep(25,17),rep(35,12),rep(45,19),rep(55,44),rep(60,28),rep(65,76),rep(75,61),rep(80,24),rep(85,67),rep(90,59),rep(95,46),rep(100,40),rep(105,9),rep(110,44),rep(115,16)) Corbula.Brijuni.M40.add.coarse.depth.for.encrusted=ceiling(Corbula.Brijuni.M40.add.fine.depth.for.encrusted/10)*10 Corbula.Brijuni.M40.add.coarse.depth.for.disornamented=Corbula.Brijuni.M40.add.coarse.depth.for.encrusted Corbula.Brijuni.M40.add.for.encrusted=c(c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,4),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,7),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,8),rep(1,1)),c(rep(0,16),rep(1,1)),c(rep(0,8),rep(1,4)),c(rep(0,14),rep(1,5)),c(rep(0,32),rep(1,12)),c(rep(0,17),rep(1,11)),c(rep(0,46),rep(1,30)),c(rep(0,34),rep(1,27)),c(rep(0,13),rep(1,11)),c(rep(0,48),rep(1,19)),c(rep(0,34),rep(1,25)),c(rep(0,29),rep(1,17)),c(rep(0,28),rep(1,12)),c(rep(0,6),rep(1,3)),c(rep(0,33),rep(1,11)),c(rep(0,8),rep(1,8))) Corbula.Brijuni.M40.add.for.disornamented=c(c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,3),rep(1,1)),c(rep(0,0),rep(1,0)),c(rep(0,7),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,6),rep(1,3)),c(rep(0,13),rep(1,4)),c(rep(0,6),rep(1,6)),c(rep(0,9),rep(1,10)),c(rep(0,18),rep(1,26)),c(rep(0,11),rep(1,17)),c(rep(0,22),rep(1,54)),c(rep(0,19),rep(1,42)),c(rep(0,9),rep(1,15)),c(rep(0,34),rep(1,33)),c(rep(0,23),rep(1,36)),c(rep(0,20),rep(1,26)),c(rep(0,17),rep(1,23)),c(rep(0,2),rep(1,7)),c(rep(0,21),rep(1,23)),c(rep(0,4),rep(1,12))) #the first number is the count of unworn encrusters, the second is the count of worn encrusters Timoclea.Brijuni.M40.add.for.enc.preserved=c(c(rep(0,2),rep(1,1)),c(rep(0,2),rep(1,2)),c(rep(0,0),rep(1,2)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,4),rep(1,2)),c(rep(0,1),rep(1,0)),c(rep(0,1),rep(1,2)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,5),rep(1,2)),c(rep(0,7),rep(1,8)),c(rep(0,10),rep(1,7)),c(rep(0,13),rep(1,14)),c(rep(0,10),rep(1,22)),c(rep(0,8),rep(1,19)),c(rep(0,17),rep(1,26)),c(rep(0,18),rep(1,18)),c(rep(0,16),rep(1,7)),c(rep(0,14),rep(1,10)),c(rep(0,11),rep(1,4)),c(rep(0,19),rep(1,4)),c(rep(0,6),rep(1,3)),c(rep(0,4),rep(1,1)),c(rep(0,6),rep(1,3)),c(rep(0,2),rep(1,1))) Timoclea.Brijuni.M40.add.fine.depth.for.enc.preserved=c(rep(2,3),rep(4,4),rep(6,2),rep(8,0),rep(10,0),rep(12,6),rep(14,1),rep(16,3),rep(18,0),rep(20,1),rep(25,7),rep(35,15),rep(45,17),rep(50,27),rep(55,32),rep(60,27),rep(65,43),rep(75,36),rep(80,23),rep(85,24),rep(90,15),rep(95,23),rep(100,9),rep(105,5),rep(110,9),rep(115,3)) Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved=ceiling(Timoclea.Brijuni.M40.add.fine.depth.for.enc.preserved/10)*10 #the first is unworn, the second is worn encruster Timoclea.Brijuni.M40.add.for.stained.enc.preserved=c(c(rep(0,1),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,4),rep(1,2)),c(rep(0,10),rep(1,4)),c(rep(0,13),rep(1,3)),c(rep(0,14),rep(1,9)),c(rep(0,23),rep(1,8)),c(rep(0,18),rep(1,6)),c(rep(0,17),rep(1,9)),c(rep(0,9),rep(1,10)),c(rep(0,3),rep(1,7)),c(rep(0,9),rep(1,9)),c(rep(0,5),rep(1,4)),c(rep(0,8),rep(1,7)),c(rep(0,3),rep(1,3)),c(rep(0,1),rep(1,0)),c(rep(0,1),rep(1,5)),c(rep(0,1),rep(1,1))) Timoclea.Brijuni.M40.add.fine.depth.for.stained.enc.preserved=c(rep(2,1),rep(4,1),rep(6,0),rep(8,0),rep(10,0),rep(12,2),rep(14,0),rep(16,2),rep(18,0),rep(20,1),rep(25,6),rep(35,14),rep(45,16),rep(50,23),rep(55,31),rep(60,24),rep(65,26),rep(75,19),rep(80,10),rep(85,18),rep(90,9),rep(95,15),rep(100,6),rep(105,1),rep(110,6),rep(115,2)) Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved=ceiling(Timoclea.Brijuni.M40.add.fine.depth.for.stained.enc.preserved/10)*10 #the first number is the count of unworn encrusters, the second is the count of worn encrusters Corbula.Brijuni.M40.add.for.enc.preserved=c(c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,1)),c(rep(0,0),rep(1,1)),c(rep(0,6),rep(1,1)),c(rep(0,1),rep(1,2)),c(rep(0,4),rep(1,6)),c(rep(0,4),rep(1,7)),c(rep(0,8),rep(1,22)),c(rep(0,3),rep(1,21)),c(rep(0,5),rep(1,4)),c(rep(0,7),rep(1,12)),c(rep(0,6),rep(1,17)),c(rep(0,6),rep(1,10)),c(rep(0,7),rep(1,6)),c(rep(0,0),rep(1,5)),c(rep(0,2),rep(1,8)),c(rep(0,1),rep(1,4))) Corbula.Brijuni.M40.add.fine.depth.for.enc.preserved=c(rep(2,0),rep(4,0),rep(6,0),rep(8,0),rep(10,0),rep(12,0),rep(14,0),rep(16,0),rep(18,0),rep(20,1),rep(25,1),rep(35,7),rep(45,3),rep(55,10),rep(60,11),rep(65,30),rep(75,24),rep(80,9),rep(85,19),rep(90,23),rep(95,16),rep(100,13),rep(105,5),rep(110,10),rep(115,5)) Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved=ceiling(Corbula.Brijuni.M40.add.fine.depth.for.enc.preserved/10)*10 #the first number is the count of stained encrusters, the second number is the count of non-stained encruster on stained shells Corbula.Brijuni.M40.add.for.stained.enc.preserved=c(c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,3),rep(1,0)),c(rep(0,2),rep(1,2)),c(rep(0,2),rep(1,2)),c(rep(0,7),rep(1,2)),c(rep(0,5),rep(1,3)),c(rep(0,12),rep(1,7)),c(rep(0,11),rep(1,5)),c(rep(0,5),rep(1,4)),c(rep(0,9),rep(1,8)),c(rep(0,10),rep(1,6)),c(rep(0,9),rep(1,6)),c(rep(0,5),rep(1,8)),c(rep(0,3),rep(1,0)),c(rep(0,7),rep(1,3)),c(rep(0,1),rep(1,1))) Corbula.Brijuni.M40.add.fine.depth.for.stained.enc.preserved=c(rep(2,0),rep(4,0),rep(6,0),rep(8,0),rep(10,0),rep(12,0),rep(14,0),rep(16,0),rep(18,0),rep(20,1),rep(25,3),rep(35,4),rep(45,4),rep(55,9),rep(60,8),rep(65,19),rep(75,16),rep(80,9),rep(85,17),rep(90,16),rep(95,15),rep(100,13),rep(105,3),rep(110,10),rep(115,2)) Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved=ceiling(Corbula.Brijuni.M40.add.fine.depth.for.stained.enc.preserved/10)*10 ########################################################################################## #PER-INCREMENT FREQUENCIES OF CORBULA ########################################################################################## Corbula.M40.ext.encrustation.inc=tapply(Corbula.Brijuni.M40.add.for.encrusted,Corbula.Brijuni.M40.add.coarse.depth.for.encrusted, mean, na.rm=T) Corbula.M40.stained.inc=tapply(Corbula.Brijuni.M40.add.for.stained,Corbula.Brijuni.M40.add.coarse.depth.for.stained, mean, na.rm=T) Corbula.M40.disornamented.inc=tapply(Corbula.Brijuni.M40.add.for.disornamented,Corbula.Brijuni.M40.add.coarse.depth.for.disornamented, mean, na.rm=T) Corbula.M44.ext.encrustation.inc=tapply(Corbula.Brijuni.M44.add.for.encrusted,Corbula.Brijuni.M44.add.coarse.depth.for.encrusted, mean, na.rm=T) Corbula.M44.stained.inc=tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T) Corbula.M44.disornamented.inc=tapply(Corbula.Brijuni.M44.add.for.disornamented,Corbula.Brijuni.M44.add.coarse.depth.for.disornamented, mean, na.rm=T) ########################################################################################## #DEFINE SYSTEMS TRACTS FOR PER-INCREMENT ANALYSES ########################################################################################## TO.hold.M44.ST=Timoclea.10cm[Timoclea.core=="Brijuni 44"] TO.hold.M44.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 44"]<21,"ML",TO.hold.M44.ST) TO.hold.M44.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 44"]<90 & Timoclea.10cm[Timoclea.core=="Brijuni 44"]>20,"HST",TO.hold.M44.ST) TO.hold.M44.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 44"]>120,"TST",TO.hold.M44.ST) TO.hold.M44.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 44"]>85 & Timoclea.10cm[Timoclea.core=="Brijuni 44"]<125,"MFZ",TO.hold.M44.ST) TO.enc.preserved.hold.M44.ST=Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved TO.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved<21,"HML",TO.enc.preserved.hold.M44.ST) TO.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved<90 & Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved>20,"HST",TO.enc.preserved.hold.M44.ST) TO.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved>120,"TST",TO.enc.preserved.hold.M44.ST) TO.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved>85 & Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved<125,"MFZ",TO.enc.preserved.hold.M44.ST) TO.stained.enc.preserved.hold.M44.ST=Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved TO.stained.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved<21,"ML",TO.stained.enc.preserved.hold.M44.ST) TO.stained.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved<90 & Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved>20,"HST",TO.stained.enc.preserved.hold.M44.ST) TO.stained.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved>120,"TST",TO.stained.enc.preserved.hold.M44.ST) TO.stained.enc.preserved.hold.M44.ST=ifelse(Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved>85 & Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved<125,"MFZ",TO.stained.enc.preserved.hold.M44.ST) TO.hold.M40.ST=Timoclea.10cm[Timoclea.core=="Brijuni 40"] TO.hold.M40.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 40"]<21,"ML",TO.hold.M44.ST) TO.hold.M40.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 40"]<45 & Timoclea.10cm[Timoclea.core=="Brijuni 40"]>20,"HST",TO.hold.M44.ST) TO.hold.M40.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 40"]>70,"TST",TO.hold.M44.ST) TO.hold.M40.ST=ifelse(Timoclea.10cm[Timoclea.core=="Brijuni 40"]>40 & Timoclea.10cm[Timoclea.core=="Brijuni 40"]<75,"MFZ",TO.hold.M40.ST) TO.enc.preserved.hold.M40.ST=Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved TO.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved<21,"ML",TO.enc.preserved.hold.M40.ST) TO.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved<45 & Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved>20,"HST",TO.enc.preserved.hold.M40.ST) TO.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved>70,"TST",TO.enc.preserved.hold.M40.ST) TO.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved>40 & Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved<75,"MFZ",TO.enc.preserved.hold.M40.ST) TO.stained.enc.preserved.hold.M40.ST=Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved TO.stained.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved<21,"ML",TO.stained.enc.preserved.hold.M40.ST) TO.stained.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved<45 & Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved>20,"HST",TO.stained.enc.preserved.hold.M40.ST) TO.stained.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved>70,"TST",TO.stained.enc.preserved.hold.M40.ST) TO.stained.enc.preserved.hold.M40.ST=ifelse(Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved>40 & Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved<75,"MFZ",TO.stained.enc.preserved.hold.M40.ST) CG.hold.M44.ST=Corbula.Brijuni.M44.add.coarse.depth.for.encrusted CG.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.encrusted<21,"ML",CG.hold.M44.ST) CG.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.encrusted<90 & Corbula.Brijuni.M44.add.coarse.depth.for.encrusted>20,"HST",CG.hold.M44.ST) CG.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.encrusted>120,"TST",CG.hold.M44.ST) CG.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.encrusted>85 & Corbula.Brijuni.M44.add.coarse.depth.for.encrusted<125,"MFZ",CG.hold.M44.ST) CG.stain.hold.M44.ST=Corbula.Brijuni.M44.add.coarse.depth.for.stained CG.stain.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained<21,"ML",CG.stain.hold.M44.ST) CG.stain.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained<90 & Corbula.Brijuni.M44.add.coarse.depth.for.stained>20,"HST",CG.stain.hold.M44.ST) CG.stain.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained>120,"TST",CG.stain.hold.M44.ST) CG.stain.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained>85 & Corbula.Brijuni.M44.add.coarse.depth.for.stained<125,"MFZ",CG.stain.hold.M44.ST) CG.enc.preserved.hold.M44.ST=Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved CG.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved<21,"ML",CG.enc.preserved.hold.M44.ST) CG.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved<90 & Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved>20,"HST",CG.enc.preserved.hold.M44.ST) CG.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved>120,"TST",CG.enc.preserved.hold.M44.ST) CG.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved>85 & Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved<125,"MFZ",CG.enc.preserved.hold.M44.ST) CG.stained.enc.preserved.hold.M44.ST=Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved CG.stained.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved<21,"ML",CG.stained.enc.preserved.hold.M44.ST) CG.stained.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved<90 & Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved>20,"HST",CG.stained.enc.preserved.hold.M44.ST) CG.stained.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved>120,"TST",CG.stained.enc.preserved.hold.M44.ST) CG.stained.enc.preserved.hold.M44.ST=ifelse(Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved>85 & Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved<125,"MFZ",CG.stained.enc.preserved.hold.M44.ST) CG.hold.M40.ST=Corbula.Brijuni.M44.add.coarse.depth.for.encrusted CG.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.encrusted<21,"ML",CG.hold.M40.ST) CG.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.encrusted<45 & Corbula.Brijuni.M40.add.coarse.depth.for.encrusted>20,"HST",CG.hold.M40.ST) CG.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.encrusted>70,"TST",CG.hold.M40.ST) CG.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.encrusted>40 & Corbula.Brijuni.M40.add.coarse.depth.for.encrusted<75,"MFZ",CG.hold.M40.ST) CG.stain.hold.M40.ST=Corbula.Brijuni.M40.add.coarse.depth.for.stained CG.stain.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained<21,"ML",CG.stain.hold.M40.ST) CG.stain.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained<45 & Corbula.Brijuni.M40.add.coarse.depth.for.stained>20,"HST",CG.stain.hold.M40.ST) CG.stain.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained>70,"TST",CG.stain.hold.M40.ST) CG.stain.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained>40 & Corbula.Brijuni.M40.add.coarse.depth.for.stained<75,"MFZ",CG.stain.hold.M40.ST) CG.enc.preserved.hold.M40.ST=Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved CG.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved<21,"ML",CG.enc.preserved.hold.M40.ST) CG.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved<45 & Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved>20,"HST",CG.enc.preserved.hold.M40.ST) CG.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved>70,"TST",CG.enc.preserved.hold.M40.ST) CG.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved>40 & Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved<75,"MFZ",CG.enc.preserved.hold.M40.ST) CG.stained.enc.preserved.hold.M40.ST=Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved CG.stained.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved<21,"ML",CG.stained.enc.preserved.hold.M40.ST) CG.stained.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved<45 & Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved>20,"HST",CG.stained.enc.preserved.hold.M40.ST) CG.stained.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved>70,"TST",CG.stained.enc.preserved.hold.M40.ST) CG.stained.enc.preserved.hold.M40.ST=ifelse(Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved>40 & Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved<75,"MFZ",CG.stained.enc.preserved.hold.M40.ST) TO.dated.M44.ST=Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"] TO.dated.M44.ST=ifelse(Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"]<21,"ML",TO.dated.M44.ST) TO.dated.M44.ST=ifelse(Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"]<90 & Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"]>20,"HST",TO.dated.M44.ST) TO.dated.M44.ST=ifelse(Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"]>120,"TST",TO.dated.M44.ST) TO.dated.M44.ST=ifelse(Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"]>85 & Timoclea.dated.depth.10cm[Timoclea.dated.core=="Brijuni 44"]<125,"MFZ",TO.dated.M44.ST) ########################################################################################################### #TAPHONOMIC TRENDS BASED ON ALL SCORED SHELLS ########################################################################################################### depths.Corbula.M44=as.numeric(names(tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean))) depths.Timoclea.M44=as.numeric(names(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean))) depths.dated=as.numeric(names(tapply(Timoclea.M44.ext.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean))) depths.Corbula.M40=as.numeric(names(tapply(Corbula.Brijuni.M40.add.for.stained,Corbula.Brijuni.M40.add.coarse.depth.for.stained, mean))) depths.Timoclea.M40=as.numeric(names(tapply(Timoclea.discoloration[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean))) ####################################################################################################################################################### #BOTH SPECIES POOLED ####################################################################################################################################################### par(mfrow=c(2,4)) par(mar=c(4,2,2,1)) plot(tapply(c(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.for.encrusted),c(Timoclea.10cm[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.coarse.depth.for.encrusted), mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Encrustation", ylab="Sediment depth (cm)", type="b", pch=16, cex=1.4) plot(tapply(c(Timoclea.Brijuni.M44.add.for.enc.preserved, Corbula.Brijuni.M44.add.for.enc.preserved),c(Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved, Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved), mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Worn encrusters", ylab="Sediment depth (cm)", type="b", pch=16, cex=1.4) plot(tapply(c(Timoclea.stained[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.for.stained),c(Timoclea.10cm[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.coarse.depth.for.stained), mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Staining", ylab="Sediment depth (cm)", type="b", pch=16, cex=1.4) plot(tapply(c(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.for.disornamented),c(Timoclea.10cm[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.coarse.depth.for.disornamented), mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Ornamentation loss", ylab="Sediment depth (cm)", type="b", pch=16, cex=1.4) ####################################################################################################################################################### #MANUSCRIPT FIGURE 13- STRATIGRAPHIC DIFFERENCES IN PER-UNIT SPECIES ALTERATION BETWEEN UNITS (SYSTEMS TRACTS) ####################################################################################################################################################### ################### #TIMOCLEA AT M44 ################### par(mfrow=c(2,4)) par(mar=c(4,2,2,1)) plot(c(1:4)-0.1,tapply(Corbula.Brijuni.M44.add.for.encrusted,CG.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)], type="n", xlab="", xaxt="n", ylim=c(0,1), pch=21, bg="black", frame=F, xlim=c(0.5,4.5), ylab="Proportion", main="Timoclea-M44", cex.main=1) axis(1, at=1:4, labels=c("TST","MFZ","HST","ML")) segments(x0=c(1:4)-0.1, x1=c(1:4)-0.1, y0=unlist(lapply(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)-0.05, x1=c(1:4)-0.05, y0=unlist(lapply(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4), x1=c(1:4), y0=unlist(lapply(tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.05, x1=c(1:4)+0.05, y0=unlist(lapply(tapply(Timoclea.Brijuni.M44.add.for.enc.preserved,TO.enc.preserved.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.Brijuni.M44.add.for.enc.preserved,TO.enc.preserved.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.1, x1=c(1:4)+0.1, y0=unlist(lapply(tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) lines(c(1:4)-0.1,tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)-0.05,tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0,tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.05,tapply(Timoclea.Brijuni.M44.add.for.enc.preserved,TO.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.1,tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) points(c(1:4)-0.1,tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="black", cex=1.4) points(c(1:4)-0.05,tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="gray", cex=1.4) points(c(1:4)+0,tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],TO.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="gray", cex=1.4) points(c(1:4)+0.05,tapply(Timoclea.Brijuni.M44.add.for.enc.preserved,TO.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="white", cex=1.4) points(c(1:4)+0.1,tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="white", cex=1.4) legend(x="topleft", pch=c(21,21,22), pt.bg=c("black","gray","gray"), bty="n", legend=c("Encrustation","Ornamentation loss","Staining")) legend(x="bottomright", pch=c(21,22), pt.bg=c("white","white"), bty="n", legend=c("Worn encrusters","Non-stained enc.")) ################### #CORBULA AT M44 ################### plot(c(1:4)-0.1,tapply(Corbula.Brijuni.M44.add.for.encrusted,CG.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)], type="n", xlab="", xaxt="n", ylim=c(0,1), pch=21, bg="black", frame=F, xlim=c(0.5,4.5), ylab="Proportion", main="Corbula-M44", cex.main=1) axis(1, at=1:4, labels=c("TST","MFZ","HST","ML")) segments(x0=c(1:4)-0.1, x1=c(1:4)-0.1, y0=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.encrusted,CG.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.encrusted,CG.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)-0.05, x1=c(1:4)-0.05, y0=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.disornamented,CG.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.disornamented,CG.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4), x1=c(1:4), y0=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.stained,CG.stain.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.stained,CG.stain.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.05, x1=c(1:4)+0.05, y0=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.enc.preserved,CG.enc.preserved.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.enc.preserved,CG.enc.preserved.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.1, x1=c(1:4)+0.1, y0=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M44.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M44.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) lines(c(1:4)-0.1,tapply(Corbula.Brijuni.M44.add.for.encrusted,CG.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)-0.05,tapply(Corbula.Brijuni.M44.add.for.disornamented,CG.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0,tapply(Corbula.Brijuni.M44.add.for.stained,CG.stain.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.05,tapply(Corbula.Brijuni.M44.add.for.enc.preserved,CG.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.1,tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)]) points(c(1:4)-0.1,tapply(Corbula.Brijuni.M44.add.for.encrusted,CG.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="black", cex=1.4) points(c(1:4)-0.05,tapply(Corbula.Brijuni.M44.add.for.disornamented,CG.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="gray", cex=1.4) points(c(1:4)+0,tapply(Corbula.Brijuni.M44.add.for.stained,CG.stain.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="gray", cex=1.4) points(c(1:4)+0.05,tapply(Corbula.Brijuni.M44.add.for.enc.preserved,CG.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="white", cex=1.4) points(c(1:4)+0.1,tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M44.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="white", cex=1.4) ################### #TIMOCLEA AT M40 ################### plot(c(1:4)-0.1,tapply(Corbula.Brijuni.M40.add.for.encrusted,CG.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)], type="n", xlab="", xaxt="n", ylim=c(0,1), pch=21, bg="black", frame=F, xlim=c(0.5,4.5), ylab="Proportion", main="Timoclea-M40", cex.main=1) axis(1, at=1:4, labels=c("TST","MFZ","HST","ML")) segments(x0=c(1:4)-0.1, x1=c(1:4)-0.1, y0=unlist(lapply(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)-0.05, x1=c(1:4)-0.05, y0=unlist(lapply(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4), x1=c(1:4), y0=unlist(lapply(tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.05, x1=c(1:4)+0.05, y0=unlist(lapply(tapply(Timoclea.Brijuni.M40.add.for.enc.preserved,TO.enc.preserved.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.Brijuni.M40.add.for.enc.preserved,TO.enc.preserved.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.1, x1=c(1:4)+0.1, y0=unlist(lapply(tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) lines(c(1:4)-0.1,tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)-0.05,tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0,tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.05,tapply(Timoclea.Brijuni.M40.add.for.enc.preserved,TO.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.1,tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) points(c(1:4)-0.1,tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="black", cex=1.4) points(c(1:4)-0.05,tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="gray", cex=1.4) points(c(1:4)+0,tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],TO.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="gray", cex=1.4) points(c(1:4)+0.05,tapply(Timoclea.Brijuni.M40.add.for.enc.preserved,TO.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="white", cex=1.4) points(c(1:4)+0.1,tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,TO.stained.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="white", cex=1.4) ################### #CORBULA AT M40 ################### plot(c(1:4)-0.1,tapply(Corbula.Brijuni.M40.add.for.encrusted,CG.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)], type="n", xlab="", xaxt="n", ylim=c(0,1), pch=21, bg="black", frame=F, xlim=c(0.5,4.5), ylab="Proportion", main="Corbula-M40", cex.main=1) axis(1, at=1:4, labels=c("TST","MFZ","HST","ML")) segments(x0=c(1:4)-0.1, x1=c(1:4)-0.1, y0=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.encrusted,CG.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.encrusted,CG.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)-0.05, x1=c(1:4)-0.05, y0=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.disornamented,CG.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.disornamented,CG.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4), x1=c(1:4), y0=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.stained,CG.stain.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.stained,CG.stain.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.05, x1=c(1:4)+0.05, y0=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.enc.preserved,CG.enc.preserved.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.enc.preserved,CG.enc.preserved.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) segments(x0=c(1:4)+0.1, x1=c(1:4)+0.1, y0=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M40.ST, bootci), quantile, 0.025))[c(4,2,1,3)],y1=unlist(lapply(tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M40.ST, bootci), quantile, 0.975))[c(4,2,1,3)]) lines(c(1:4)-0.1,tapply(Corbula.Brijuni.M40.add.for.encrusted,CG.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)-0.05,tapply(Corbula.Brijuni.M40.add.for.disornamented,CG.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0,tapply(Corbula.Brijuni.M40.add.for.stained,CG.stain.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.05,tapply(Corbula.Brijuni.M40.add.for.enc.preserved,CG.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) lines(c(1:4)+0.1,tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)]) points(c(1:4)-0.1,tapply(Corbula.Brijuni.M40.add.for.encrusted,CG.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="black", cex=1.4) points(c(1:4)-0.05,tapply(Corbula.Brijuni.M40.add.for.disornamented,CG.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="gray", cex=1.4) points(c(1:4)+0,tapply(Corbula.Brijuni.M40.add.for.stained,CG.stain.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="gray", cex=1.4) points(c(1:4)+0.05,tapply(Corbula.Brijuni.M40.add.for.enc.preserved,CG.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=21, bg="white", cex=1.4) points(c(1:4)+0.1,tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,CG.stained.enc.preserved.hold.M40.ST, mean, na.rm=T)[c(4,2,1,3)],pch=22, bg="white", cex=1.4) ####################################################################################################################################################### #MANUSCRIPT FIGURE 12 - STRATIGRAPHIC TRENDS IN PER-INCREMENT ALTERATION FOR ALL SHELLS ####################################################################################################################################################### #BRIJUNI M44 par(mfrow=c(2,5)) par(mar=c(4,2,2,0)) ########################### #EXTERNAL ENCRUSTATION ########################### plot(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0), xlim=c(0,1), frame=F, xlab="Encrustation", ylab="Sediment depth (cm)", pch=16, cex=1.4) rect(xleft=0,xright=1,ytop=-120, ybottom=-160, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-90, ybottom=-120, border=NA, col="gray61") lines(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44) points(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44,pch=16, cex=1.2) out=tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44) out=tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 44" & Timoclea.length>5], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44) lines(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 44" & Timoclea.length>5], mean, na.rm=T),-depths.Timoclea.M44, col="black") points(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 44" & Timoclea.length>5], mean, na.rm=T),-depths.Timoclea.M44, pch=21, bg="white", cex=1.2) out=tapply(Corbula.Brijuni.M44.add.for.encrusted,Corbula.Brijuni.M44.add.coarse.depth.for.encrusted, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44, col="gray") lines(tapply(Corbula.Brijuni.M44.add.for.encrusted,Corbula.Brijuni.M44.add.coarse.depth.for.encrusted, mean, na.rm=T),-depths.Corbula.M44, col="gray") points(tapply(Corbula.Brijuni.M44.add.for.encrusted,Corbula.Brijuni.M44.add.coarse.depth.for.encrusted, mean, na.rm=T),-depths.Corbula.M44, bg="gray", pch=21, cex=1.2) ########################### #LOSS OF ORNAMENTATION ########################### plot(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Worn shells", ylab="Sediment depth (cm)", pch=16, cex=1.4) rect(xleft=0,xright=1,ytop=-120, ybottom=-160, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-90, ybottom=-120, border=NA, col="gray61") points(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44, pch=16, cex=1.2) lines(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44) out=tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44) out=tapply(Corbula.Brijuni.M44.add.for.disornamented,Corbula.Brijuni.M44.add.coarse.depth.for.disornamented, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44, col="gray") lines(tapply(Corbula.Brijuni.M44.add.for.disornamented,Corbula.Brijuni.M44.add.coarse.depth.for.disornamented, mean, na.rm=T),-depths.Corbula.M44, col="gray") points(tapply(Corbula.Brijuni.M44.add.for.disornamented,Corbula.Brijuni.M44.add.coarse.depth.for.disornamented, mean, na.rm=T),-depths.Corbula.M44, bg="gray", pch=21, cex=1.2) ########################### #PERCENT OF WORN ENCRUSTERS ########################### plot(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Worn encrusters", ylab="Sediment depth (cm)", pch=16, cex=1.4, type="n") rect(xleft=0,xright=1,ytop=-120, ybottom=-160, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-90, ybottom=-120, border=NA, col="gray61") out=tapply(Timoclea.Brijuni.M44.add.for.enc.preserved, Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44, col="black") out=tapply(Corbula.Brijuni.M44.add.for.enc.preserved, Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Corbula.M44, y1=-depths.Corbula.M44, col="gray") lines(tapply(Timoclea.Brijuni.M44.add.for.enc.preserved, Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved, mean, rm=T),-depths.Timoclea.M44, col="black") points(tapply(Timoclea.Brijuni.M44.add.for.enc.preserved, Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved, mean, na.rm=T),-depths.Timoclea.M44, bg="black", pch=21, cex=1.2) lines(tapply(Corbula.Brijuni.M44.add.for.enc.preserved, Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved, mean, rm=T),-depths.Corbula.M44, col="gray") points(tapply(Corbula.Brijuni.M44.add.for.enc.preserved, Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved, mean, na.rm=T),-depths.Corbula.M44, bg="gray", pch=21, cex=1.2) ########################### #STAINING ########################### plot(tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Shell staining", ylab="Sediment depth (cm)", pch=16, cex=1.4) rect(xleft=0,xright=1,ytop=-120, ybottom=-160, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-90, ybottom=-120, border=NA, col="gray61") points(tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44, pch=16, cex=1.2) out=tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44) lines(tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44) #lines(tapply(Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T),-depths.dated) #points(tapply(Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T),-depths.dated, bg="white", pch=21, cex=1.2) out=tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44, col="gray") lines(tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T),-depths.Corbula.M44, col="gray") points(tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T),-depths.Corbula.M44, bg="gray", pch=21, cex=1.2) legend(x="topright", pt.bg=c("gray","black"), pch=21, legend=c("Corbula","Timoclea"), cex=1.2) ########################################## #NON-STAINED ENCRUSTERS ON STAINED SHELLS ########################################## plot(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T),-depths.Timoclea.M44, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Non-stained encrusters", ylab="Sediment depth (cm)", pch=16, cex=1.4, type="n") rect(xleft=0,xright=1,ytop=-120, ybottom=-160, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-90, ybottom=-120, border=NA, col="gray61") out=tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M44, y1=-depths.Timoclea.M44, col="black") lines(tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-depths.Timoclea.M44, col="black") points(tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-depths.Timoclea.M44, bg="black", pch=21, cex=1.4) out=tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Corbula.M40, y1=-depths.Corbula.M40, col="gray") lines(tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-depths.Corbula.M44, col="gray") points(tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-depths.Corbula.M44, bg="gray", pch=21, cex=1.4) ####################################################################################################################################################### #BRIJUNI M40 ####################################################################################################################################################### depths=as.numeric(names(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean))) plot(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, ylim=c(-160,0),xlim=c(0,1), frame=F, xlab="Encrustation", ylab="Sediment depth (cm)", pch=16, cex=1.4) rect(xleft=0,xright=1,ytop=-75, ybottom=-120, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-45, ybottom=-75, border=NA, col="gray61") points(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, pch=16, cex=1.4) lines(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40) out=tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M40, y1=-depths.Timoclea.M40) out=tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 40" & Timoclea.length>5], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M40, y1=-depths.Timoclea.M40) lines(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 40" & Timoclea.length>5], mean, na.rm=T),-depths.Timoclea.M40[1:11], col="black") points(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 40" & Timoclea.length>5], mean, na.rm=T),-depths.Timoclea.M40[1:11], pch=21, bg="white", cex=1.4) depths.Corbula.M40.encrusted=as.numeric(names(tapply(Corbula.Brijuni.M40.add.for.encrusted,Corbula.Brijuni.M40.add.coarse.depth.for.encrusted, mean))) out=tapply(Corbula.Brijuni.M40.add.for.encrusted,Corbula.Brijuni.M40.add.coarse.depth.for.encrusted, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Corbula.M40, y1=-depths.Corbula.M40, col="gray") lines(tapply(Corbula.Brijuni.M40.add.for.encrusted,Corbula.Brijuni.M40.add.coarse.depth.for.encrusted, mean, na.rm=T),-depths.Corbula.M40.encrusted, col="gray") points(tapply(Corbula.Brijuni.M40.add.for.encrusted,Corbula.Brijuni.M40.add.coarse.depth.for.encrusted, mean, na.rm=T),-depths.Corbula.M40.encrusted, bg="gray", pch=21, cex=1.4) ########################### #LOSS OF ORNAMENTATION ########################### plot(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, ylim=c(-160,0), xlim=c(0,1), frame=F, xlab="Worn shells", ylab="Sediment depth (cm)", pch=16, cex=1.4) rect(xleft=0,xright=1,ytop=-75, ybottom=-120, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-45, ybottom=-75, border=NA, col="gray61") points(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, pch=16, cex=1.4) lines(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40) out=tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M40, y1=-depths.Timoclea.M40) out=tapply(Corbula.Brijuni.M40.add.for.disornamented,Corbula.Brijuni.M40.add.coarse.depth.for.disornamented, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Corbula.M40, y1=-depths.Corbula.M40, col="gray") lines(tapply(Corbula.Brijuni.M40.add.for.disornamented,Corbula.Brijuni.M40.add.coarse.depth.for.disornamented, mean, na.rm=T),-depths.Corbula.M40, col="gray") points(tapply(Corbula.Brijuni.M40.add.for.disornamented,Corbula.Brijuni.M40.add.coarse.depth.for.disornamented, mean, na.rm=T),-depths.Corbula.M40, bg="gray", pch=21, cex=1.4) ########################### #PERCENT OF WORN ENCRUSTERS ########################### plot(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, ylim=c(-160,0), xlim=c(0,1), frame=F, xlab="Worn encrusters", ylab="Sediment depth (cm)", pch=16, cex=1.4, type="n") rect(xleft=0,xright=1,ytop=-75, ybottom=-120, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-45, ybottom=-75, border=NA, col="gray61") out=tapply(Timoclea.Brijuni.M40.add.for.enc.preserved, Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M40, y1=-depths.Timoclea.M40, col="black") temp=as.numeric(names(tapply(Corbula.Brijuni.M40.add.for.enc.preserved, Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved, mean, rm=T))) out=tapply(Corbula.Brijuni.M40.add.for.enc.preserved, Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-temp, y1=-temp, col="gray") lines(tapply(Timoclea.Brijuni.M40.add.for.enc.preserved, Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved, mean, rm=T),-depths.Timoclea.M40, col="black") points(tapply(Timoclea.Brijuni.M40.add.for.enc.preserved, Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved, mean, na.rm=T),-depths.Timoclea.M40, bg="black", pch=21, cex=1.4) lines(tapply(Corbula.Brijuni.M40.add.for.enc.preserved, Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved, mean, rm=T),-temp, col="gray") points(tapply(Corbula.Brijuni.M40.add.for.enc.preserved, Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved, mean, na.rm=T),-temp, bg="gray", pch=21, cex=1.4) ########################### #STAINING ########################### plot(tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, ylim=c(-160,0), xlim=c(0,1), frame=F, xlab="Shell staining", ylab="Sediment depth (cm)", pch=16, cex=1.4) rect(xleft=0,xright=1,ytop=-75, ybottom=-120, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-45, ybottom=-75, border=NA, col="gray61") points(tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, pch=16, cex=1.4) lines(tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40) out=tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M40, y1=-depths.Timoclea.M40) out=tapply(Corbula.Brijuni.M40.add.for.stained,Corbula.Brijuni.M40.add.coarse.depth.for.stained, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Corbula.M40, y1=-depths.Corbula.M40, col="gray") lines(tapply(Corbula.Brijuni.M40.add.for.stained,Corbula.Brijuni.M40.add.coarse.depth.for.stained, mean, na.rm=T),-depths.Corbula.M40, col="gray") points(tapply(Corbula.Brijuni.M40.add.for.stained,Corbula.Brijuni.M40.add.coarse.depth.for.stained, mean, na.rm=T),-depths.Corbula.M40, bg="gray", pch=21, cex=1.4) ########################################## #NON-STAINED ENCRUSTERS ON STAINED SHELLS ########################################## plot(tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T),-depths.Timoclea.M40, ylim=c(-160,0), xlim=c(0,1), frame=F, xlab="Non-stained encrusters", ylab="Sediment depth (cm)", pch=16, cex=1.4, type="n") rect(xleft=0,xright=1,ytop=-75, ybottom=-120, border=NA, col="gray31") rect(xleft=0,xright=1,ytop=-45, ybottom=-75, border=NA, col="gray61") out=tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-depths.Timoclea.M40, y1=-depths.Timoclea.M40, col="black") lines(tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-depths.Timoclea.M40, col="black") points(tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-depths.Timoclea.M40, bg="black", pch=21, cex=1.4) temp=as.numeric(names(tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved, Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, mean, rm=T))) out=tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, bootci) LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T)) UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T)) segments(x0=LCI,x1=UCI,y0=-temp, y1=-temp, col="gray") lines(tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-temp, col="gray") points(tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),-temp, bg="gray", pch=21, cex=1.4) ################################################################################################################################## #MANUSCRIPT FIGURE 13 LOWER PART - MULTIVARIATE ASSEMBLAGE-LEVEL ORDINATION ################################################################################################################################## #cairo_pdf(file='plot.pdf', width=7, height=7) par(mfcol=c(2,2)) par(mar=c(4,4,2,1)) multi=array(NA, dim=c(28,5)) depth.labels=c(-depths.Timoclea.M44, -depths.Timoclea.M40) core.labels=c(rep("M44",length(depths.Timoclea.M44)), rep("M40", length(depths.Timoclea.M40))) rownames(multi)=depth.labels colnames(multi)=c("encrusted","worn","worn encrusters","stained","non-stained encrusters") multi[,1]=c(tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T), tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T)) multi[,2]=c(tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T), tapply(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T)) multi[,3]=c(tapply(Timoclea.Brijuni.M44.add.for.enc.preserved, Timoclea.Brijuni.M44.add.coarse.depth.for.enc.preserved, mean, rm=T), tapply(Timoclea.Brijuni.M40.add.for.enc.preserved, Timoclea.Brijuni.M40.add.coarse.depth.for.enc.preserved, mean, rm=T)) multi[,4]=c(tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T), tapply(Timoclea.stained[Timoclea.core=="Brijuni 40"],Timoclea.10cm[Timoclea.core=="Brijuni 40"], mean, na.rm=T)) multi[,5]=c(tapply(Timoclea.Brijuni.M44.add.for.stained.enc.preserved,Timoclea.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T), tapply(Timoclea.Brijuni.M40.add.for.stained.enc.preserved,Timoclea.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T)) core.labels=c(rep("M44",length(depths.Timoclea.M44)-2), rep("M40", length(depths.Timoclea.M40))) ST.labels=c(rep("HST",6),rep("MFZ",4),rep("TST",4), rep("HST",4),rep("MFZ",3),rep("TST",5)) ST.labels=c(rep("ML",2),rep("HST",4),rep("MFZ",4),rep("TST",4), rep("ML",2),rep("HST",2),rep("MFZ",3),rep("TST",5)) assemblage.pco=cmdscale(vegdist(multi[3:28,], method="manhattan"), eig=T) assemblage.pco$eig[1]/sum(assemblage.pco$eig[assemblage.pco$eig>0]) assemblage.pco$eig[2]/sum(assemblage.pco$eig[assemblage.pco$eig>0]) #assemblage.pco=isoMDS(vegdist(multi[3:28,], method="manhattan"), k=2) plot(assemblage.pco$points[,1],assemblage.pco$points[,2], ylim=c(-0.8,0.6), xlim=c(-0.8,0.6), xlab="NMDS axis 1", ylab="NMDS axis 2", frame=F, pch=16, asp=1, main="Timoclea ovata", cex.main=0.9) text(assemblage.pco$points[,1],assemblage.pco$points[,2], labels=depth.labels[3:28], cex=0.8, pos=3) #for (i in 1:13) {arrows(x0=assemblage.pco$points[i,1], y0=assemblage.pco$points[i,2], x1 = assemblage.pco$points[i+1,1], y1 = assemblage.pco$points[i+1,2], length = 0.25)} #for (i in 15:25) {arrows(x0=assemblage.pco$points[i,1], y0=assemblage.pco$points[i,2], x1 = assemblage.pco$points[i+1,1], y1 = assemblage.pco$points[i+1,2], length = 0.25)} points(assemblage.pco$points[core.labels=="M44",1],assemblage.pco$points[core.labels=="M44",2], pch=21, bg="black", cex=1.4) points(assemblage.pco$points[core.labels=="M40",1],assemblage.pco$points[core.labels=="M40",2], pch=21, bg="gray", cex=1.4) ordihull(ord=assemblage.pco, group=ST.labels, draw = "polygon", col=c("gray21","gray51","white","gray81")) plot(assemblage.pco$points[,1],assemblage.pco$points[,2], ylim=c(-0.8,0.6), xlim=c(-0.8,0.6), xlab="NMDS axis 1", ylab="NMDS axis 2", frame=F, pch=16, asp=1, main="Timoclea ovata", cex.main=0.9, type="n") envi=envfit(assemblage.pco$points, multi[3:28,]) plot(envi, add=T, cex=1, col="gray31", lwd=2, arrow.mul=0.5) envi Timoclea.PERMANOVA=adonis(vegdist(multi[3:28,], method="manhattan")~ST.labels) multi=array(NA, dim=c(28,5)) multi[,1]=c(tapply(Corbula.Brijuni.M44.add.for.encrusted,Corbula.Brijuni.M44.add.coarse.depth.for.encrusted, mean, na.rm=T), tapply(Corbula.Brijuni.M40.add.for.encrusted,Corbula.Brijuni.M40.add.coarse.depth.for.encrusted, mean, na.rm=T)) multi[,2]=c(tapply(Corbula.Brijuni.M44.add.for.disornamented,Corbula.Brijuni.M44.add.coarse.depth.for.disornamented, mean, na.rm=T), tapply(Corbula.Brijuni.M40.add.for.disornamented,Corbula.Brijuni.M40.add.coarse.depth.for.disornamented, mean, na.rm=T)) multi[,3]=c(tapply(Corbula.Brijuni.M44.add.for.enc.preserved, Corbula.Brijuni.M44.add.coarse.depth.for.enc.preserved, mean, rm=T),0, tapply(Corbula.Brijuni.M40.add.for.enc.preserved, Corbula.Brijuni.M40.add.coarse.depth.for.enc.preserved, mean, rm=T)) multi[,4]=c(tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T), tapply(Corbula.Brijuni.M40.add.for.stained,Corbula.Brijuni.M40.add.coarse.depth.for.stained, mean, na.rm=T)) multi[,5]=c(tapply(Corbula.Brijuni.M44.add.for.stained.enc.preserved,Corbula.Brijuni.M44.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T),0, tapply(Corbula.Brijuni.M40.add.for.stained.enc.preserved,Corbula.Brijuni.M40.add.coarse.depth.for.stained.enc.preserved, mean, na.rm=T)) rownames(multi)=depth.labels colnames(multi)=c("encrusted","worn","worn encrusters","stained","non-stained encrusters") assemblage.pco=cmdscale(vegdist(multi[3:28,], method="manhattan"), eig=T) assemblage.pco$eig[1]/sum(assemblage.pco$eig[assemblage.pco$eig>0]) assemblage.pco$eig[2]/sum(assemblage.pco$eig[assemblage.pco$eig>0]) plot(assemblage.pco$points[,1],assemblage.pco$points[,2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, asp=1, main="Corbula gibba", cex.main=0.9) text(assemblage.pco$points[,1],assemblage.pco$points[,2], labels=depth.labels[3:28], cex=0.8, pos=3) #for (i in 1:13) {arrows(x0=assemblage.pco$points[i,1], y0=assemblage.pco$points[i,2], x1 = assemblage.pco$points[i+1,1], y1 = assemblage.pco$points[i+1,2], length = 0.25)} #for (i in 15:25) {arrows(x0=assemblage.pco$points[i,1], y0=assemblage.pco$points[i,2], x1 = assemblage.pco$points[i+1,1], y1 = assemblage.pco$points[i+1,2], length = 0.25)} points(assemblage.pco$points[core.labels=="M44",1],assemblage.pco$points[core.labels=="M44",2], pch=21, bg="black", cex=1.4) points(assemblage.pco$points[core.labels=="M40",1],assemblage.pco$points[core.labels=="M40",2], pch=21, bg="gray", cex=1.4) ordihull(ord=assemblage.pco, group=ST.labels, draw = "polygon", col=c("gray21","gray51","white","gray81")) legend(x="bottomleft", legend=c("M44","M40"), pt.bg=c("black","gray"), pch=21, bty="n", cex=1.4) plot(assemblage.pco$points[,1],assemblage.pco$points[,2], xlab="PCO axis 1", ylab="PCO axis 2", frame=F, pch=16, asp=1, main="Corbula gibba", cex.main=0.9, type="n") envi=envfit(assemblage.pco$points, multi[3:28,]) plot(envi, add=T, cex=1, col="gray31", lwd=2, arrow.mul=0.5) Corbula.PERMANOVA=adonis(vegdist(multi[3:28,], method="manhattan")~ST.labels) #dev.off() Timoclea.PERMANOVA Corbula.PERMANOVA ########################################################################################################### #TOTAL ABUNDANCE OF ALL MOLLUSCS PER 10 CM INCREMENT - PROXY FOR SHELLINESS ########################################################################################################### #increments=10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 Brijuni.M44.10cm.total.ab=c(489, 587, 717, 1062, 979, 1193, 1189, 1680, 1919, 1973, 1936, 1912, 2832, 2412, 1472, 471) ########################################################################################################### #PERCENT GRAVEL - ANOTHER PROXY FOR SHELLINESS ########################################################################################################### Brijuni.M44.10cm.gravel=c(9.34, 6.52, 11.80, 14.30, 16.65, 17.10, 16.45, 23.80, 28.65, 34.35, 49.90, 51.60, 49.85, 50.60, 20.10, 1.70) Brijuni.M44.10cm.depth=c(6, 16, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155) ################################################################################################################################## #CORBULA AND TIMOCLEA POOLED - RELATIONS BETWEEN THE PERCENT SHELL GRAVEL AND TOTAL ABUNDANCE OF MOLLUSCS VS ALTERATION ################################################################################################################################## mean.encrustation.M44.2.species=tapply(c(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.for.encrusted),c(Timoclea.10cm[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.coarse.depth.for.encrusted), mean, na.rm=T)[-16] mean.staining.M44.2.species=tapply(c(Timoclea.stained[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.for.stained),c(Timoclea.10cm[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.coarse.depth.for.stained), mean, na.rm=T)[-16] mean.ornamentation.M44.2.species=tapply(c(Timoclea.ornamentation[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.for.disornamented),c(Timoclea.10cm[Timoclea.core=="Brijuni 44"],Corbula.Brijuni.M44.add.coarse.depth.for.disornamented), mean, na.rm=T)[-16] mean.encrustation.M40.2.species=tapply(c(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 40"],Corbula.Brijuni.M40.add.for.encrusted),c(Timoclea.10cm[Timoclea.core=="Brijuni 40"],Corbula.Brijuni.M40.add.coarse.depth.for.encrusted), mean, na.rm=T) mean.staining.M40.2.species=tapply(c(Timoclea.stained[Timoclea.core=="Brijuni 40"],Corbula.Brijuni.M40.add.for.stained),c(Timoclea.10cm[Timoclea.core=="Brijuni 40"],Corbula.Brijuni.M40.add.coarse.depth.for.stained), mean, na.rm=T) mean.ornamentation.M40.2.species=tapply(c(Timoclea.ornamentation[Timoclea.core=="Brijuni 40"],Corbula.Brijuni.M40.add.for.disornamented),c(Timoclea.10cm[Timoclea.core=="Brijuni 40"],Corbula.Brijuni.M40.add.coarse.depth.for.disornamented), mean, na.rm=T) par(mfcol=c(2,2)) par(mar=c(4,4,2,1)) plot(Brijuni.M44.10cm.gravel[-16],mean.encrustation.M44.2.species, frame=F, ylab="Encrustation", xlab="Percent shell gravel", pch=16, cex=1.4, ylim=c(0,0.5)) plot(Brijuni.M44.10cm.total.ab[-16],mean.encrustation.M44.2.species, frame=F, ylab="Encrustation", xlab="Total abundance", pch=16, cex=1.4, ylim=c(0,0.5)) plot(Brijuni.M44.10cm.gravel[-16],mean.staining.M44.2.species, frame=F, ylab="Staining", xlab="Percent shell gravel", pch=16, cex=1.4, ylim=c(0,0.5)) plot(Brijuni.M44.10cm.total.ab[-16],mean.staining.M44.2.species, frame=F, ylab="Staining", xlab="Total abundance", pch=16, cex=1.4, ylim=c(0,0.5)) #cor.test(Brijuni.M44.10cm.gravel[-16],mean.encrustation.M44.2.species, method="s") #cor.test(Brijuni.M44.10cm.total.ab[-16],mean.encrustation.M44.2.species, method="s") #cor.test(Brijuni.M44.10cm.gravel[-16],mean.staining.M44.2.species, method="s") #cor.test(Brijuni.M44.10cm.total.ab[-16],mean.staining.M44.2.species, method="s") ################################################################################################################################## #TAPHONOMIC FEEDBACK CORBULA AND TIMOCLEA SEPARATELY ################################################################################################################################## par(mfrow=c(2,2)) plot(Brijuni.M44.10cm.gravel,-Brijuni.M44.10cm.depth, pch=16, frame=F, xlab="Shell gravel") lines(Brijuni.M44.10cm.gravel,-Brijuni.M44.10cm.depth) plot(Brijuni.M44.10cm.total.ab,-Brijuni.M44.10cm.depth, pch=16, frame=F, xlab="Total abundance") lines(Brijuni.M44.10cm.total.ab,-Brijuni.M44.10cm.depth) #ENCRUSTATION par(mfrow=c(3,4)) par(mar=c(4,2,2,1)) plot(Brijuni.M44.10cm.gravel[-16], tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T)[-16],frame=F, ylab="Timoclea encrustation", xlab="Percent shell gravel", pch=16, cex=1.4, ylim=c(0,1), main="Timoclea", cex.main=0.9) points(Brijuni.M44.10cm.gravel[-16], tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 44" & Timoclea.length>5], mean, na.rm=T)[-16],pch=21, bg="gray", cex=1.4) plot(Brijuni.M44.10cm.total.ab[-16], tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T)[-16], frame=F, ylab="Timoclea encrustation", xlab="Total molluscan abundance", pch=16, cex=1.4, ylim=c(0,1), main="Timoclea", cex.main=0.9) points(Brijuni.M44.10cm.total.ab[-16], tapply(Timoclea.ext.encrustation[Timoclea.core=="Brijuni 44" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 44" & Timoclea.length>5], mean, na.rm=T)[-16],pch=21, bg="gray", cex=1.4) plot(Brijuni.M44.10cm.gravel[-16], tapply(Corbula.Brijuni.M44.add.for.encrusted,Corbula.Brijuni.M44.add.coarse.depth.for.encrusted, mean, na.rm=T)[-16], frame=F, ylab="Corbula encrustation", xlab="Percent shell gravel", pch=16, cex=1.4, ylim=c(0,1), main="Corbula", cex.main=0.9) plot(Brijuni.M44.10cm.total.ab[-16], tapply(Corbula.Brijuni.M44.add.for.encrusted,Corbula.Brijuni.M44.add.coarse.depth.for.encrusted, mean, na.rm=T)[-16], frame=F, ylab="Corbula encrustation", xlab="Total molluscan abundance", pch=16, cex=1.4, ylim=c(0,1), main="Corbula", cex.main=0.9) #STAINING plot(Brijuni.M44.10cm.gravel[-16], tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T)[-16], frame=F, ylab="Timoclea staining", xlab="Percent shell gravel", pch=16, cex=1.4, ylim=c(0,1)) points(Brijuni.M44.10cm.gravel[-16], tapply(Timoclea.stained[Timoclea.core=="Brijuni 44" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 44" & Timoclea.length>5], mean, na.rm=T)[-16],pch=21, bg="gray", cex=1.4) legend(x="topright",legend=c("> 1 mm","> 5 mm"), pt.bg=c("black","gray"), pch=21, bty="n", cex=1.4) plot(Brijuni.M44.10cm.total.ab[-16], tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T)[-16], frame=F, ylab="Timoclea staining", xlab="Total molluscan abundance", pch=16, cex=1.4, ylim=c(0,1)) points(Brijuni.M44.10cm.total.ab[-16], tapply(Timoclea.stained[Timoclea.core=="Brijuni 44" & Timoclea.length>5],Timoclea.10cm[Timoclea.core=="Brijuni 44" & Timoclea.length>5], mean, na.rm=T)[-16],pch=21, bg="gray", cex=1.4) plot(Brijuni.M44.10cm.gravel[-16], tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T)[-16], frame=F, ylab="Corbula staining", xlab="Percent shell gravel", pch=16, cex=1.4, ylim=c(0,1)) plot(Brijuni.M44.10cm.total.ab[-16], tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T)[-16], frame=F, ylab="Corbula staining", xlab="Total molluscan abundance", pch=16, cex=1.4, ylim=c(0,1)) #cor.test(Brijuni.M44.10cm.gravel[-16], tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T)[-16],method="s") #cor.test(Brijuni.M44.10cm.total.ab[-16], tapply(Timoclea.stained[Timoclea.core=="Brijuni 44"],Timoclea.10cm[Timoclea.core=="Brijuni 44"], mean, na.rm=T)[-16],method="s") #cor.test(Brijuni.M44.10cm.gravel[-16], tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T)[-16],method="s") #cor.test(Brijuni.M44.10cm.total.ab[-16], tapply(Corbula.Brijuni.M44.add.for.stained,Corbula.Brijuni.M44.add.coarse.depth.for.stained, mean, na.rm=T)[-16],method="s") ########################################################################### #PER-INCREMENT ALTERATION IN TIMOCLEA ########################################################################### Timoclea.M44.N.inc=tapply(Timoclea.M44.pyrite[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], length) Timoclea.M44.pyrite.inc=tapply(Timoclea.M44.pyrite[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.discoloration.inc=tapply(Timoclea.M44.discoloration[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.ornamentation.inc=tapply(Timoclea.M44.ornamentation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.ext.encrustation.inc=tapply(Timoclea.M44.ext.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.int.encrustation.inc=tapply(Timoclea.M44.int.encrustation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.fragmentation.inc=tapply(Timoclea.M44.fragmentation[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.stained.inc=tapply(Timoclea.M44.stained[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.bioerosion.inc=tapply(Timoclea.M44.bioerosion[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.dissolution.inc=tapply(Timoclea.M44.dissolution[!is.na(Timoclea.M44.ages.for.clock)],Timoclea.M44.10cm[!is.na(Timoclea.M44.ages.for.clock)], mean, na.rm=T) Timoclea.M44.inc.scores=cbind(Timoclea.M44.discoloration.inc, Timoclea.M44.ornamentation.inc, Timoclea.M44.ext.encrustation.inc,Timoclea.M44.int.encrustation.inc, Timoclea.M44.fragmentation.inc, Timoclea.M44.stained.inc, Timoclea.M44.bioerosion.inc, Timoclea.M44.dissolution.inc) Timoclea.M44.inc.scores.with.pyrite=cbind(Timoclea.M44.discoloration.inc, Timoclea.M44.ornamentation.inc, Timoclea.M44.ext.encrustation.inc,Timoclea.M44.int.encrustation.inc, Timoclea.M44.fragmentation.inc, Timoclea.M44.stained.inc, Timoclea.M44.bioerosion.inc, Timoclea.M44.dissolution.inc, Timoclea.M44.pyrite.inc) Timoclea.M40.N.inc=tapply(Timoclea.M40.pyrite[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], length) Timoclea.M40.pyrite.inc=tapply(Timoclea.M40.pyrite[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.discoloration.inc=tapply(Timoclea.M40.discoloration[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.ornamentation.inc=tapply(Timoclea.M40.ornamentation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.ext.encrustation.inc=tapply(Timoclea.M40.ext.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.int.encrustation.inc=tapply(Timoclea.M40.int.encrustation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.fragmentation.inc=tapply(Timoclea.M40.fragmentation[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.stained.inc=tapply(Timoclea.M40.stained[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.bioerosion.inc=tapply(Timoclea.M40.bioerosion[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.dissolution.inc=tapply(Timoclea.M40.dissolution[!is.na(Timoclea.M40.ages.for.clock)],Timoclea.M40.10cm[!is.na(Timoclea.M40.ages.for.clock)], mean, na.rm=T) Timoclea.M40.inc.scores=cbind(Timoclea.M40.discoloration.inc, Timoclea.M40.ornamentation.inc, Timoclea.M40.ext.encrustation.inc,Timoclea.M40.int.encrustation.inc, Timoclea.M40.fragmentation.inc, Timoclea.M40.stained.inc, Timoclea.M40.bioerosion.inc, Timoclea.M40.dissolution.inc) Timoclea.M40.inc.scores.with.pyrite=cbind(Timoclea.M40.discoloration.inc, Timoclea.M40.ornamentation.inc, Timoclea.M40.ext.encrustation.inc,Timoclea.M40.int.encrustation.inc, Timoclea.M40.fragmentation.inc, Timoclea.M40.stained.inc, Timoclea.M40.bioerosion.inc, Timoclea.M40.dissolution.inc, Timoclea.M40.pyrite.inc) ################################################################################################################################## #PAIRWISE PER-ASSEMBLAGE RELATIONS BASED ON ALL BRIJUNI SHELLS ################################################################################################################################## par(mfcol=c(2,3)) plot(c(Corbula.M44.ext.encrustation.inc,Corbula.M44.ext.encrustation.inc),c(Corbula.M44.stained.inc,Corbula.M44.stained.inc), pch=16, xlab="Corbula ext. encrustation", ylab="Corbula staining", frame=F, xlim=c(0,1), ylim=c(0,1)) lines(c(0,1),c(0,1), lty=2) plot(c(Timoclea.M44.ext.encrustation.inc,Timoclea.M40.ext.encrustation.inc),c(Timoclea.M44.stained.inc,Timoclea.M40.stained.inc), pch=16, xlab="Timoclea ext. encrustation", ylab="Timoclea staining", frame=F, xlim=c(0,1), ylim=c(0,1)) lines(c(0,1),c(0,1), lty=2) points(Timoclea.M44.ext.encrustation.inc,Timoclea.M44.stained.inc, pch=21, bg="white") points(Timoclea.M40.ext.encrustation.inc,Timoclea.M40.stained.inc, pch=16) plot(c(Corbula.M44.ext.encrustation.inc,Corbula.M40.ext.encrustation.inc),c(Corbula.M44.disornamented.inc,Corbula.M40.disornamented.inc), pch=16, xlab="Corbula ext. encrustation", ylab="Corbula ornamentation loss", frame=F, xlim=c(0,1), ylim=c(0,1)) lines(c(0,1),c(0,1), lty=2) plot(c(Timoclea.M44.ext.encrustation.inc,Timoclea.M40.ext.encrustation.inc),c(Timoclea.M44.ornamentation.inc,Timoclea.M40.ornamentation.inc), pch=16, xlab="Timoclea ext. encrustation", ylab="Timoclea ornamentation loss", frame=F, xlim=c(0,1), ylim=c(0,1)) lines(c(0,1),c(0,1), lty=2) plot(c(Corbula.M44.stained.inc,Corbula.M40.stained.inc),c(Corbula.M44.disornamented.inc,Corbula.M40.disornamented.inc), pch=16, xlab="Corbula staining", ylab="Corbula ornamentation loss", frame=F, xlim=c(0,1), ylim=c(0,1)) lines(c(0,1),c(0,1), lty=2) plot(c(Timoclea.M44.stained.inc,Timoclea.M40.stained.inc),c(Timoclea.M44.ornamentation.inc,Timoclea.M40.ornamentation.inc), pch=16, xlab="Timoclea staining", ylab="Timoclea ornamentation loss", frame=F, xlim=c(0,1), ylim=c(0,1)) lines(c(0,1),c(0,1), lty=2)