require(smatr)
require(MASS)
require(vegan)

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
}


bootci.ext=function(temp) {
x.resample<-numeric(); xmean.resample<-numeric(); xmedian.resample<-numeric();x25.resample<-numeric();
x75.resample<-numeric();x99.resample<-numeric();x95.resample<-numeric();x05.resample<-numeric();iqr.resample<-numeric();range.resample<-numeric();
x975.resample<-numeric();x025.resample<-numeric();
for (i in 1:1000) {
for (j in 1:length(temp)) {x.resample[j]<-sample(temp,1,replace=T)}
xmean.resample[i]<-mean(x.resample, na.rm=T)
xmedian.resample[i]<-median(x.resample, na.rm=T)
x25.resample[i]<-quantile(x.resample, c(0.25), na.rm=T)
x75.resample[i]<-quantile(x.resample, c(0.75), na.rm=T)
iqr.resample[i]=x75.resample[i]-x25.resample[i]
x95.resample[i]<-quantile(x.resample, c(0.95), na.rm=T)
x99.resample[i]<-quantile(x.resample, c(0.99), na.rm=T)
x975.resample[i]<-quantile(x.resample, c(0.975), na.rm=T)
x05.resample[i]<-quantile(x.resample, c(0.05), na.rm=T)
x025.resample[i]<-quantile(x.resample, c(0.025), na.rm=T)
range.resample[i]=x975.resample[i]-x025.resample[i]
}
out=list(resample=xmean.resample, resamplemedian=xmedian.resample, resample25=x25.resample, resample75=x75.resample, resampleiqr=iqr.resample, 
resamplerange=range.resample, resample25=x25.resample, resample99=x99.resample,resample95=x95.resample,resample05=x05.resample)
out
}

##################################################################################################################################
#INDIVIDUAL-BASED SCORED DATASET
##################################################################################################################################
setwd("D:/Data/Adriatic_sea_project/Taph scoring")
taphonomic.data <- read.delim("Table S2-Taphonomic scores of Varicorbula-pyrite linings.txt", header=TRUE)
ages.for.clock=taphonomic.data[,"Age"]
specimen=taphonomic.data$"specimen_no"
total.acids.for.clock=taphonomic.data[,"X.total."]
taxon=taphonomic.data$taxon
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"
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"
oxides=taphonomic.data$"iron.oxide.coating"
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."

##################################################################################################################################
Corbula.cementation=cementation[taxon=="Corbula gibba"]
Corbula.pyrite=pyrite[taxon=="Corbula gibba"]
Corbula.periostracum=periostracum[taxon=="Corbula gibba"]
Corbula.conchiolin=conchiolin[taxon=="Corbula gibba"]
Corbula.int.conch=int.conch[taxon=="Corbula gibba"]
Corbula.discoloration=discoloration[taxon=="Corbula gibba"]
Corbula.ornamentation=ornamentation[taxon=="Corbula gibba"]
Corbula.ext.encrustation=ext.encrustation[taxon=="Corbula gibba"]
Corbula.int.encrustation=int.encrustation[taxon=="Corbula gibba"]
Corbula.core=core[taxon=="Corbula gibba"]
Corbula.fragmentation=fragmentation[taxon=="Corbula gibba"]
Corbula.stained=stained[taxon=="Corbula gibba"]
Corbula.bioerosion=bioerosion[taxon=="Corbula gibba"]
Corbula.dissolution=dissolution[taxon=="Corbula gibba"]
Corbula.oxides=oxides[taxon=="Corbula gibba"]
Corbula.depths=depths[taxon=="Corbula gibba"]
Corbula.10cm=d10cm[taxon=="Corbula gibba"]
Corbula.5cm=d5cm[taxon=="Corbula gibba"]
Corbula.ages.for.clock=ages.for.clock[taxon=="Corbula gibba"]
Corbula.total.acids.for.clock=total.acids.for.clock[taxon=="Corbula gibba"]

Corbula.M44.cementation=cementation[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.pyrite=pyrite[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.discoloration=discoloration[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.ornamentation=ornamentation[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.conchiolin=conchiolin[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.int.conch=int.conch[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.ext.encrustation=ext.encrustation[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.int.encrustation=int.encrustation[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.core=core[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.fragmentation=fragmentation[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.stained=stained[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.bioerosion=bioerosion[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.dissolution=dissolution[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.depths=depths[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.10cm=d10cm[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.5cm=d5cm[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.finecm=dfinecm[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.periostracum=periostracum[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.int.conch=int.conch[taxon=="Corbula gibba" & core=="Brijuni 44"]
Corbula.M44.oxides=oxides[taxon=="Corbula gibba" & core=="Brijuni 44"]

Corbula.Piran.cementation=cementation[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.pyrite=pyrite[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.periostracum=periostracum[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.conchiolin=conchiolin[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.int.conch=int.conch[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.discoloration=discoloration[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.ornamentation=ornamentation[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.ext.encrustation=ext.encrustation[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.int.encrustation=int.encrustation[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.core=core[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.fragmentation=fragmentation[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.stained=stained[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.bioerosion=bioerosion[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.dissolution=dissolution[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.oxides=oxides[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.depths=depths[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.10cm=d10cm[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.5cm=d5cm[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.finecm=dfinecm[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.ages.for.clock=ages.for.clock[taxon=="Corbula gibba" & core=="PiranII M53"]
Corbula.Piran.total.acids.for.clock=total.acids.for.clock[taxon=="Corbula gibba" & core=="PiranII M53"]

Corbula.Panzano.cementation=ornamentation[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.pyrite=pyrite[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.periostracum=periostracum[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.conchiolin=conchiolin[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.int.conch=int.conch[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.discoloration=discoloration[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.ornamentation=ornamentation[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.ext.encrustation=ext.encrustation[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.int.encrustation=int.encrustation[taxon=="Corbula gibba" & core=="Pan M283"]
Corbula.Panzano.core=core[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.fragmentation=fragmentation[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.stained=stained[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.bioerosion=bioerosion[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.dissolution=dissolution[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.oxides=oxides[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.depths=depths[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.10cm=d10cm[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.5cm=d5cm[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.finecm=dfinecm[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.ages.for.clock=ages.for.clock[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.total.acids.for.clock=total.acids.for.clock[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.M28.pyrite=pyrite[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.M28.depths=depths[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.M28.10cm=d10cm[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.M28.5cm=d5cm[taxon=="Corbula gibba" & core=="Pan M28"]
Corbula.Panzano.M28.finecm=dfinecm[taxon=="Corbula gibba" & core=="Pan M28"]

Corbula.Po3.cementation=cementation[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.pyrite=pyrite[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.periostracum=periostracum[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.conchiolin=conchiolin[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.int.conch=int.conch[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.discoloration=discoloration[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.ornamentation=ornamentation[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.ext.encrustation=ext.encrustation[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.int.encrustation=int.encrustation[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.core=core[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.fragmentation=fragmentation[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.stained=stained[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.bioerosion=bioerosion[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.dissolution=dissolution[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.oxides=oxides[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.depths=depths[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.10cm=d10cm[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.5cm=d5cm[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.finecm=dfinecm[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.ages.for.clock=ages.for.clock[taxon=="Corbula gibba" & core=="Po3 M13"]
Corbula.Po3.total.acids.for.clock=total.acids.for.clock[taxon=="Corbula gibba" & core=="Po3 M13"]

Corbula.Po4.cementation=cementation[taxon=="Corbula gibba" & core=="Po3 M21"]
Corbula.Po4.pyrite=pyrite[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.periostracum=periostracum[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.conchiolin=conchiolin[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.int.conch=int.conch[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.discoloration=discoloration[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.ornamentation=ornamentation[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.ext.encrustation=ext.encrustation[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.int.encrustation=int.encrustation[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.core=core[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.fragmentation=fragmentation[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.stained=stained[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.bioerosion=bioerosion[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.dissolution=dissolution[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.oxides=oxides[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.depths=depths[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.10cm=d10cm[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.5cm=d5cm[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.finecm=dfinecm[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.ages.for.clock=ages.for.clock[taxon=="Corbula gibba" & core=="Po4 M21"]
Corbula.Po4.total.acids.for.clock=total.acids.for.clock[taxon=="Corbula gibba" & core=="Po4 M21"]
###########################################################################################################
#AMINO ACIDS
###########################################################################################################
par(mfrow=c(2,2))
plot(Corbula.Po4.ages.for.clock, Corbula.Po4.total.acids.for.clock, pch=16, log="xy", frame=F)
plot(Corbula.Po3.ages.for.clock, Corbula.Po3.total.acids.for.clock, pch=16, log="xy", frame=F)
plot(Corbula.Panzano.ages.for.clock, Corbula.Panzano.total.acids.for.clock, pch=16, log="xy", frame=F)
plot(Corbula.Piran.ages.for.clock, Corbula.Piran.total.acids.for.clock, pch=16, log="xy", frame=F)

Corbula.Po4.ages.10y=cut(Corbula.Po4.ages.for.clock, breaks=seq(0,15000,by=10), labels=F)*10
Corbula.Po3.ages.10y=cut(Corbula.Po3.ages.for.clock, breaks=seq(0,15000,by=10), labels=F)*10
Corbula.Po4.ages.25y=cut(Corbula.Po4.ages.for.clock, breaks=seq(0,15000,by=25), labels=F)*25
Corbula.Po3.ages.25y=cut(Corbula.Po3.ages.for.clock, breaks=seq(0,15000,by=25), labels=F)*25
Corbula.Panzano.ages.10y=cut(Corbula.Panzano.ages.for.clock, breaks=seq(0,15000,by=10), labels=F)*10
Corbula.Panzano.ages.25y=cut(Corbula.Panzano.ages.for.clock, breaks=seq(0,15000,by=25), labels=F)*25
Corbula.Piran.ages.10y=cut(Corbula.Piran.ages.for.clock, breaks=seq(0,15000,by=10), labels=F)*10
Corbula.Piran.ages.25y=cut(Corbula.Piran.ages.for.clock, breaks=seq(0,15000,by=25), labels=F)*25
Corbula.Po4.ages.20y=cut(Corbula.Po4.ages.for.clock, breaks=seq(0,15000,by=20), labels=F)*20
Corbula.Po3.ages.20y=cut(Corbula.Po3.ages.for.clock, breaks=seq(0,15000,by=20), labels=F)*20
Corbula.Panzano.ages.20y=cut(Corbula.Panzano.ages.for.clock, breaks=seq(0,15000,by=20), labels=F)*20
Corbula.Piran.ages.20y=cut(Corbula.Piran.ages.for.clock, breaks=seq(0,15000,by=20), labels=F)*20
Corbula.POS.ages.20y=cut(Corbula.POS.ages.for.clock, breaks=seq(0,15000,by=20), labels=F)*20
Corbula.POS.ages.10y=cut(Corbula.POS.ages.for.clock, breaks=seq(0,15000,by=10), labels=F)*10
Corbula.POS.ages.25y=cut(Corbula.POS.ages.for.clock, breaks=seq(0,15000,by=25), labels=F)*25
Corbula.Panzano.ages.50y=cut(Corbula.Panzano.ages.for.clock, breaks=seq(0,15000,by=50), labels=F)*50

#############################################################################################################################################
#FINAL FIGURE 11
#############################################################################################################################################
par(mfrow=c(3,3))
par(mar=c(4,4,2,1))
boxplot(split(Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==0], Corbula.Po3.ages.10y[Corbula.Po3.pyrite==0])[1:12], 
at=as.numeric(names(split(Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==0], Corbula.Po3.ages.10y[Corbula.Po3.pyrite==0])))[1:12], 
range=0, boxwex=3, xlab="Age (years, based on Asp D/L)", main="Po 3", cex.main=1, col="white", log="y", frame=F, ylim=c(1000,100000), ylab="Total concentration of amino acids")
boxplot(split(Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==1], Corbula.Po3.ages.10y[Corbula.Po3.pyrite==1]+5)[1:12],  
at=as.numeric(names(split(Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==1], Corbula.Po3.ages.10y[Corbula.Po3.pyrite==1]+5)))[1:12], 
names=rep("",12), range=0, boxwex=3, col="gray", add=T, frame=F, ylim=c(1000,100000), yaxt="n", xaxt="n")
legend(x="bottomright", pch=22, pt.bg=c("white","gray"), legend=c("without pyrite","with pyrite"), bty="n", cex=1.2)

boxplot(split(Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==0], Corbula.Po4.ages.10y[Corbula.Po4.pyrite==0])[1:12], 
at=as.numeric(names(split(Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==0], Corbula.Po4.ages.10y[Corbula.Po4.pyrite==0])))[1:12], 
range=0, boxwex=3, xlab="Age (years, based on Asp D/L)", main="Po 4", cex.main=1, col="white", log="y", frame=F, ylim=c(1000,100000), ylab="")
boxplot(split(Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==1], Corbula.Po4.ages.10y[Corbula.Po4.pyrite==1]+5)[1:12],  
at=as.numeric(names(split(Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==1], Corbula.Po4.ages.10y[Corbula.Po4.pyrite==1]+5)))[1:12], 
names=rep("",12), range=0, boxwex=3, col="gray", add=T, frame=F, ylim=c(1000,100000), yaxt="n", xaxt="n")

boxplot(split(Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==0], Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==0])[1:9], 
at=as.numeric(names(split(Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==0], Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==0])))[1:9], 
range=0, boxwex=10, xlab="Age (years, based on Asp D/L)", main="Panzano", cex.main=1, col="white", log="y", frame=F, ylim=c(1000,100000), ylab="")
boxplot(split(Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==1], Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==1]+5)[1:9],  
at=as.numeric(names(split(Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==1], Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==1]+10)))[1:9], 
names=rep("",9), range=0, boxwex=10, col="gray", add=T, frame=F, ylim=c(1000,100000), yaxt="n", xaxt="n")

plot(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==0], Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==0],
xlab="Age (years, based on Asp D/L)", main="", cex.main=1, bg="white", pch=21, log="xy", frame=F, ylim=c(1000,100000), ylab="Total concentration of amino acids", cex=1.2)
points(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==1], Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==1],  bg="gray81", pch=21, cex=1.2)
out1=lm(log(Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==0])~log(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==0]))
out2=lm(log(Corbula.Po3.total.acids.for.clock[Corbula.Po3.pyrite==1])~log(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==1]))

temp.predict=predict.lm(out1)
temp.predict.se=predict.lm(out1, se.fit=T)$se.fit
temp.ages=Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==0]
temp.ages=temp.ages[!is.na(temp.ages)]
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]), lwd=2)
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]+2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2)
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]-2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2)
temp.predict=predict.lm(out2)
temp.predict.se=predict.lm(out2, se.fit=T)$se.fit
temp.ages=Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==1]
temp.ages=temp.ages[!is.na(temp.ages)]
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]), lwd=3, col="gray51")
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]+2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2, col="gray41")
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]-2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2, col="gray41")

plot(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==0], Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==0],
xlab="Age (years, based on Asp D/L)", main="", cex.main=1, bg="white", pch=21, log="xy", frame=F, ylim=c(1000,100000), ylab="", cex=1.2)
points(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==1], Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==1],  bg="gray81", pch=21, cex=1.2)
out1=lm(log(Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==0])~log(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==0]))
out2=lm(log(Corbula.Po4.total.acids.for.clock[Corbula.Po4.pyrite==1])~log(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==1]))
legend(x="bottomleft", col=c("black","gray51"), lwd=c(2,3), legend=c("without pyrite","with pyrite"), bty="n")

temp.predict=predict.lm(out1)
temp.predict.se=predict.lm(out1, se.fit=T)$se.fit
temp.ages=Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==0]
temp.ages=temp.ages[!is.na(temp.ages)]
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]), lwd=2)
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]+2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2)
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]-2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2)
temp.predict=predict.lm(out2)
temp.predict.se=predict.lm(out2, se.fit=T)$se.fit
temp.ages=Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==1]
temp.ages=temp.ages[!is.na(temp.ages)]
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]), lwd=3, col="gray51")
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]+2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2, col="gray51")
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]-2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2, col="gray51")

plot(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0], Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==0], xlim=c(1,500), 
xlab="Age (years, based on Asp D/L)", main="", cex.main=1, bg="white", pch=21, log="xy", frame=F, ylim=c(1000,100000), ylab="", cex=1.2)
points(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1], Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==1],  bg="gray81", pch=21, cex=1.2)
out1=lm(log(Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==0])~log(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0]))
out2=lm(log(Corbula.Panzano.total.acids.for.clock[Corbula.Panzano.pyrite==1])~log(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1]))

temp.predict=predict.lm(out1)
temp.predict.se=predict.lm(out1, se.fit=T)$se.fit
temp.ages=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0]
temp.ages=temp.ages[!is.na(temp.ages)]
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]), lwd=2)
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]+2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2)
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]-2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2)
temp.predict=predict.lm(out2)
temp.predict.se=predict.lm(out2, se.fit=T)$se.fit
temp.ages=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1]
temp.ages=temp.ages[!is.na(temp.ages)]

lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]), lwd=3, col="gray51")
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]+2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2, col="gray51")
lines(temp.ages[order(temp.predict)], exp(temp.predict[order(temp.predict)]-2*temp.predict.se[order(temp.predict)]), lwd=2, lty=2, col="gray51")

###############################
#PER-STATION PERIOSTRACUM LOSS
###############################
mean.conchiolin.Po3.0=lapply(split(Corbula.Po3.conchiolin[Corbula.Po3.pyrite==0],Corbula.Po3.ages.10y[Corbula.Po3.pyrite==0]), mean)
length.conchiolin.Po3.0=tapply(Corbula.Po3.conchiolin[Corbula.Po3.pyrite==0],Corbula.Po3.ages.10y[Corbula.Po3.pyrite==0], length)
out=tapply(Corbula.Po3.conchiolin[Corbula.Po3.pyrite==0],Corbula.Po3.ages.10y[Corbula.Po3.pyrite==0], bootci)
LCI.0=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI.0=unlist(lapply(out,quantile,c(0.975), na.rm=T))
mean.conchiolin.Po3.1=lapply(split(Corbula.Po3.conchiolin[Corbula.Po3.pyrite==1],Corbula.Po3.ages.10y[Corbula.Po3.pyrite==1]), mean)
length.conchiolin.Po3.1=tapply(Corbula.Po3.conchiolin[Corbula.Po3.pyrite==1],Corbula.Po3.ages.10y[Corbula.Po3.pyrite==1], length)
out=tapply(Corbula.Po3.conchiolin[Corbula.Po3.pyrite==1],Corbula.Po3.ages.10y[Corbula.Po3.pyrite==1], bootci)
LCI.1=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI.1=unlist(lapply(out,quantile,c(0.975), na.rm=T))
select.0=which(length.conchiolin.Po3.0>2)
select.1=which(length.conchiolin.Po3.1>2)
plot(as.numeric(names(mean.conchiolin.Po3.0)), mean.conchiolin.Po3.0, type="n", pch=16, frame=F, ylim=c(0,1), xlim=c(0, 150), xlab="Age (years)", ylab="Frequency of valves with periostracum", cex=1.4)
lines(as.numeric(names(mean.conchiolin.Po3.0))[select.0], mean.conchiolin.Po3.0[select.0])
lines(as.numeric(names(mean.conchiolin.Po3.1))[select.1]+3, mean.conchiolin.Po3.1[select.1])
segments(x0 = as.numeric(names(mean.conchiolin.Po3.0))[select.0], x1 = as.numeric(names(mean.conchiolin.Po3.0))[select.0], y0 = LCI.0[select.0], y1 = UCI.0[select.0])
segments(x0 = as.numeric(names(mean.conchiolin.Po3.1))[select.1]+3, x1 = as.numeric(names(mean.conchiolin.Po3.1))[select.1]+3, y0 = LCI.1[select.1], y1 = UCI.1[select.1])
points(as.numeric(names(mean.conchiolin.Po3.0)[select.0]), mean.conchiolin.Po3.0[select.0], pch=21, bg="black", cex=1.2)
points(as.numeric(names(mean.conchiolin.Po3.1))[select.1]+3, mean.conchiolin.Po3.1[select.1], pch=21, bg="gray", cex=1.2)

mean.conchiolin.Po4.0=lapply(split(Corbula.Po4.conchiolin[Corbula.Po4.pyrite==0],Corbula.Po4.ages.10y[Corbula.Po4.pyrite==0]), mean)
out=tapply(Corbula.Po4.conchiolin[Corbula.Po4.pyrite==0],Corbula.Po4.ages.10y[Corbula.Po4.pyrite==0], bootci)
LCI.0=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI.0=unlist(lapply(out,quantile,c(0.975), na.rm=T))
mean.conchiolin.Po4.1=lapply(split(Corbula.Po4.conchiolin[Corbula.Po4.pyrite==1],Corbula.Po4.ages.10y[Corbula.Po4.pyrite==1]), mean)
out=tapply(Corbula.Po4.conchiolin[Corbula.Po4.pyrite==1],Corbula.Po4.ages.10y[Corbula.Po4.pyrite==1], bootci)
LCI.1=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI.1=unlist(lapply(out,quantile,c(0.975), na.rm=T))
length.conchiolin.Po4.0=tapply(Corbula.Po4.conchiolin[Corbula.Po4.pyrite==0],Corbula.Po4.ages.10y[Corbula.Po4.pyrite==0], length)
length.conchiolin.Po4.1=tapply(Corbula.Po4.conchiolin[Corbula.Po4.pyrite==1],Corbula.Po4.ages.10y[Corbula.Po4.pyrite==1], length)
select.0=which(length.conchiolin.Po4.0>2)
select.1=which(length.conchiolin.Po4.1>2)
plot(as.numeric(names(mean.conchiolin.Po4.0)), mean.conchiolin.Po4.0, type="n", pch=16, frame=F, ylim=c(0,1), xlim=c(0, 150), xlab="Age (years)", ylab="", cex=1.4)
segments(x0 = as.numeric(names(mean.conchiolin.Po4.0))[select.0], x1 = as.numeric(names(mean.conchiolin.Po4.0))[select.0], y0 = LCI.0[select.0], y1 = UCI.0[select.0])
segments(x0 = as.numeric(names(mean.conchiolin.Po4.1))[select.1]+3, x1 = as.numeric(names(mean.conchiolin.Po4.1))[select.1]+3, y0 = LCI.1[select.1], y1 = UCI.1[select.1])
lines(as.numeric(names(mean.conchiolin.Po4.0))[select.0], mean.conchiolin.Po4.0[select.0])
lines(as.numeric(names(mean.conchiolin.Po4.1))[select.1]+3, mean.conchiolin.Po4.1[select.1])
points(as.numeric(names(mean.conchiolin.Po4.0))[select.0], mean.conchiolin.Po4.0[select.0], pch=21, bg="black", cex=1.2)
points(as.numeric(names(mean.conchiolin.Po4.1))[select.1]+3, mean.conchiolin.Po4.1[select.1], pch=21, bg="gray", cex=1.2)

mean.conchiolin.Panzano.0=lapply(split(Corbula.Panzano.conchiolin[Corbula.Panzano.pyrite==0],Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==0]), mean)
out=tapply(Corbula.Panzano.conchiolin[Corbula.Panzano.pyrite==0],Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==0], bootci)
LCI.0=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI.0=unlist(lapply(out,quantile,c(0.975), na.rm=T))
mean.conchiolin.Panzano.1=lapply(split(Corbula.Panzano.conchiolin[Corbula.Panzano.pyrite==1],Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==1]), mean)
out=tapply(Corbula.Panzano.conchiolin[Corbula.Panzano.pyrite==1],Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==1], bootci)
LCI.1=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI.1=unlist(lapply(out,quantile,c(0.975), na.rm=T))
length.conchiolin.Panzano.0=tapply(Corbula.Panzano.conchiolin[Corbula.Panzano.pyrite==0],Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==0], length)
length.conchiolin.Panzano.1=tapply(Corbula.Panzano.conchiolin[Corbula.Panzano.pyrite==1],Corbula.Panzano.ages.50y[Corbula.Panzano.pyrite==1], length)
select.0=which(length.conchiolin.Panzano.0>2)
select.1=which(length.conchiolin.Panzano.1>2)
plot(as.numeric(names(mean.conchiolin.Panzano.0)), mean.conchiolin.Panzano.0, pch=16, type="n", frame=F, ylim=c(0,1), xlim=c(0, 500), xlab="Age (years)", ylab="", cex=1.4)
segments(x0 = as.numeric(names(mean.conchiolin.Panzano.0))[select.0], x1 = as.numeric(names(mean.conchiolin.Panzano.0))[select.0], y0 = LCI.0[select.0], y1 = UCI.0[select.0])
segments(x0 = as.numeric(names(mean.conchiolin.Panzano.1))[select.1]+5, x1 = as.numeric(names(mean.conchiolin.Panzano.1))[select.1]+5, y0 = LCI.1[select.1], y1 = UCI.1[select.1])
lines(as.numeric(names(mean.conchiolin.Panzano.0))[select.0], mean.conchiolin.Panzano.0[select.0])
lines(as.numeric(names(mean.conchiolin.Panzano.1))[select.1]+5, mean.conchiolin.Panzano.1[select.1])
points(as.numeric(names(mean.conchiolin.Panzano.0))[select.0], mean.conchiolin.Panzano.0[select.0], pch=21, bg="black", cex=1.2)
points(as.numeric(names(mean.conchiolin.Panzano.1))[select.1]+5, mean.conchiolin.Panzano.1[select.1], pch=21, bg="gray", cex=1.2)
####################################################################################################
#END OF FIGURE 11
####################################################################################################

###########################################################################################################
#ASSEMBLAGE-LEVEL INCREMENT-SPECIFIC ALTERATION MEANS BASED ON DATED SHELLS
###########################################################################################################
Corbula.Po3.N.inc=tapply(Corbula.Po3.pyrite[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], length)
Corbula.Po3.pyrite.5.inc=tapply(Corbula.Po3.pyrite[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.5cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.pyrite.inc=tapply(Corbula.Po3.pyrite[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.conchiolin.inc=tapply(Corbula.Po3.conchiolin[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.discoloration.inc=tapply(Corbula.Po3.discoloration[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.ornamentation.inc=tapply(Corbula.Po3.ornamentation[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.ext.encrustation.inc=tapply(Corbula.Po3.ext.encrustation[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.int.encrustation.inc=tapply(Corbula.Po3.int.encrustation[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.fragmentation.inc=tapply(Corbula.Po3.fragmentation[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.stained.inc=tapply(Corbula.Po3.stained[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.bioerosion.inc=tapply(Corbula.Po3.bioerosion[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.dissolution.inc=tapply(Corbula.Po3.dissolution[!is.na(Corbula.Po3.ages.for.clock)],Corbula.Po3.10cm[!is.na(Corbula.Po3.ages.for.clock)], mean, na.rm=T)
Corbula.Po3.inc.scores=cbind(Corbula.Po3.discoloration.inc, Corbula.Po3.ornamentation.inc, Corbula.Po3.ext.encrustation.inc,Corbula.Po3.int.encrustation.inc,
Corbula.Po3.fragmentation.inc, Corbula.Po3.stained.inc, Corbula.Po3.bioerosion.inc, Corbula.Po3.dissolution.inc)
Corbula.Po3.inc.scores.with.pyrite=cbind(Corbula.Po3.discoloration.inc, Corbula.Po3.ornamentation.inc, Corbula.Po3.ext.encrustation.inc,Corbula.Po3.int.encrustation.inc,
Corbula.Po3.fragmentation.inc, Corbula.Po3.stained.inc, Corbula.Po3.bioerosion.inc, Corbula.Po3.dissolution.inc, Corbula.Po3.pyrite.inc)

Corbula.Po4.N.inc=tapply(Corbula.Po4.pyrite[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], length)
Corbula.Po4.5.pyrite.inc=tapply(Corbula.Po4.pyrite[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.5cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.pyrite.inc=tapply(Corbula.Po4.pyrite[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.conchiolin.inc=tapply(Corbula.Po4.conchiolin[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.discoloration.inc=tapply(Corbula.Po4.discoloration[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.ornamentation.inc=tapply(Corbula.Po4.ornamentation[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.ext.encrustation.inc=tapply(Corbula.Po4.ext.encrustation[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.int.encrustation.inc=tapply(Corbula.Po4.int.encrustation[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.fragmentation.inc=tapply(Corbula.Po4.fragmentation[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.stained.inc=tapply(Corbula.Po4.stained[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.bioerosion.inc=tapply(Corbula.Po4.bioerosion[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.dissolution.inc=tapply(Corbula.Po4.dissolution[!is.na(Corbula.Po4.ages.for.clock)],Corbula.Po4.10cm[!is.na(Corbula.Po4.ages.for.clock)], mean, na.rm=T)
Corbula.Po4.inc.scores=cbind(Corbula.Po4.discoloration.inc, Corbula.Po4.ornamentation.inc, Corbula.Po4.ext.encrustation.inc,Corbula.Po4.int.encrustation.inc,
Corbula.Po4.fragmentation.inc, Corbula.Po4.stained.inc, Corbula.Po4.bioerosion.inc, Corbula.Po4.dissolution.inc)
Corbula.Po4.inc.scores.with.pyrite=cbind(Corbula.Po4.discoloration.inc, Corbula.Po4.ornamentation.inc, Corbula.Po4.ext.encrustation.inc,Corbula.Po4.int.encrustation.inc,
Corbula.Po4.fragmentation.inc, Corbula.Po4.stained.inc, Corbula.Po4.bioerosion.inc, Corbula.Po4.dissolution.inc, Corbula.Po4.pyrite.inc)

Corbula.Panzano.N.inc=tapply(Corbula.Panzano.pyrite[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], length)
Corbula.Panzano.pyrite.5.inc=tapply(Corbula.Panzano.pyrite[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.5cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.pyrite.inc=tapply(Corbula.Panzano.pyrite[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.conchiolin.inc=tapply(Corbula.Panzano.conchiolin[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.discoloration.inc=tapply(Corbula.Panzano.discoloration[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.ornamentation.inc=tapply(Corbula.Panzano.ornamentation[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.ext.encrustation.inc=tapply(Corbula.Panzano.ext.encrustation[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.int.encrustation.inc=tapply(Corbula.Panzano.int.encrustation[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.fragmentation.inc=tapply(Corbula.Panzano.fragmentation[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.stained.inc=tapply(Corbula.Panzano.stained[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.bioerosion.inc=tapply(Corbula.Panzano.bioerosion[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.dissolution.inc=tapply(Corbula.Panzano.dissolution[!is.na(Corbula.Panzano.ages.for.clock)],Corbula.Panzano.10cm[!is.na(Corbula.Panzano.ages.for.clock)], mean, na.rm=T)
Corbula.Panzano.inc.scores=cbind(Corbula.Panzano.discoloration.inc, Corbula.Panzano.ornamentation.inc, Corbula.Panzano.ext.encrustation.inc,Corbula.Panzano.int.encrustation.inc,
Corbula.Panzano.fragmentation.inc, Corbula.Panzano.stained.inc, Corbula.Panzano.bioerosion.inc, Corbula.Panzano.dissolution.inc)
Corbula.Panzano.inc.scores.with.pyrite=cbind(Corbula.Panzano.discoloration.inc, Corbula.Panzano.ornamentation.inc, Corbula.Panzano.ext.encrustation.inc,Corbula.Panzano.int.encrustation.inc,
Corbula.Panzano.fragmentation.inc, Corbula.Panzano.stained.inc, Corbula.Panzano.bioerosion.inc, Corbula.Panzano.dissolution.inc, Corbula.Panzano.pyrite.inc)

Corbula.Piran.N.inc=tapply(Corbula.Piran.pyrite[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], length)
Corbula.Piran.pyrite.5.inc=tapply(Corbula.Piran.pyrite[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.5cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.pyrite.inc=tapply(Corbula.Piran.pyrite[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.conchiolin.inc=tapply(Corbula.Piran.conchiolin[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.discoloration.inc=tapply(Corbula.Piran.discoloration[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.ornamentation.inc=tapply(Corbula.Piran.ornamentation[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.ext.encrustation.inc=tapply(Corbula.Piran.ext.encrustation[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.int.encrustation.inc=tapply(Corbula.Piran.int.encrustation[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.fragmentation.inc=tapply(Corbula.Piran.fragmentation[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.stained.inc=tapply(Corbula.Piran.stained[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.bioerosion.inc=tapply(Corbula.Piran.bioerosion[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.dissolution.inc=tapply(Corbula.Piran.dissolution[!is.na(Corbula.Piran.ages.for.clock)],Corbula.Piran.10cm[!is.na(Corbula.Piran.ages.for.clock)], mean, na.rm=T)
Corbula.Piran.inc.scores=cbind(Corbula.Piran.discoloration.inc, Corbula.Piran.ornamentation.inc, Corbula.Piran.ext.encrustation.inc,Corbula.Piran.int.encrustation.inc,
Corbula.Piran.fragmentation.inc, Corbula.Piran.stained.inc, Corbula.Piran.bioerosion.inc, Corbula.Piran.dissolution.inc)
Corbula.Piran.inc.scores.with.pyrite=cbind(Corbula.Piran.discoloration.inc, Corbula.Piran.ornamentation.inc, Corbula.Piran.ext.encrustation.inc,Corbula.Piran.int.encrustation.inc,
Corbula.Piran.fragmentation.inc, Corbula.Piran.stained.inc, Corbula.Piran.bioerosion.inc, Corbula.Piran.dissolution.inc, Corbula.Piran.pyrite.inc)

Corbula.Brijuni.N.inc=tapply(Corbula.M44.pyrite,Corbula.M44.10cm, length)
Corbula.Brijuni.pyrite.5.inc=tapply(Corbula.M44.pyrite,Corbula.M44.5cm, mean, na.rm=T)
Corbula.Brijuni.pyrite.inc=tapply(Corbula.M44.pyrite,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.conchiolin.inc=tapply(Corbula.M44.conchiolin,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.discoloration.inc=tapply(Corbula.M44.discoloration,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.ornamentation.inc=tapply(Corbula.M44.ornamentation,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.ext.encrustation.inc=tapply(Corbula.M44.ext.encrustation,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.int.encrustation.inc=tapply(Corbula.M44.int.encrustation,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.fragmentation.inc=tapply(Corbula.M44.fragmentation,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.stained.inc=tapply(Corbula.M44.stained,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.bioerosion.inc=tapply(Corbula.M44.bioerosion,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.dissolution.inc=tapply(Corbula.M44.dissolution,Corbula.M44.10cm, mean, na.rm=T)
Corbula.Brijuni.inc.scores=cbind(Corbula.Brijuni.discoloration.inc, Corbula.Brijuni.ornamentation.inc, Corbula.Brijuni.ext.encrustation.inc,Corbula.Brijuni.int.encrustation.inc,
Corbula.Brijuni.fragmentation.inc, Corbula.Brijuni.stained.inc, Corbula.Brijuni.bioerosion.inc, Corbula.Brijuni.dissolution.inc)
Corbula.Brijuni.inc.scores.with.pyrite=cbind(Corbula.Brijuni.discoloration.inc, Corbula.Brijuni.ornamentation.inc, Corbula.Brijuni.ext.encrustation.inc,Corbula.Brijuni.int.encrustation.inc,
Corbula.Brijuni.fragmentation.inc, Corbula.Brijuni.stained.inc, Corbula.Brijuni.bioerosion.inc, Corbula.Brijuni.dissolution.inc, Corbula.Brijuni.pyrite.inc)

##########################################################################################
#ADDITIONAL SPECIMENS FOR PYRITIZATION STUDY
##########################################################################################
#PO 3 - M12 - 2 cm and 5 cm depth 
Corbula.Po3.M12.add.fine.depth.for.pyrite=c(rep(2,0),rep(4,0),rep(6,1),rep(8,3),rep(10,4),rep(12,4),rep(14,7),rep(16,17),rep(18,5),rep(20,5),rep(25,9),rep(30,6),rep(35,12),rep(40,3),rep(45,7),rep(50,10),rep(55,17),rep(60,29),rep(65,12),rep(70,18),rep(75,11),rep(80,44),rep(85,20),rep(90,21),rep(95,9),rep(100,7),rep(105,8),rep(110,4),rep(115,0),rep(120,7),rep(125,6),rep(130,0),rep(135,1),rep(140,1),rep(145,3),rep(150,1))
Corbula.Po3.M12.add.for.pyrite=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,1),rep(1,2)),c(rep(0,4),rep(1,0)),c(rep(0,4),rep(1,0)),c(rep(0,7),rep(1,0)),c(rep(0,8),rep(1,9)),c(rep(0,5),rep(1,0)),c(rep(0,3),rep(1,2)),c(rep(0,6),rep(1,3)),c(rep(0,1),rep(1,5)),c(rep(0,4),rep(1,8)),c(rep(0,1),rep(1,2)),c(rep(0,2),rep(1,5)),c(rep(0,5),rep(1,5)),c(rep(0,9),rep(1,8)),c(rep(0,5),rep(1,24)),c(rep(0,5),rep(1,7)),c(rep(0,10),rep(1,8)),c(rep(0,6),rep(1,5)),c(rep(0,20),rep(1,24)),c(rep(0,8),rep(1,12)),c(rep(0,11),rep(1,10)),c(rep(0,7),rep(1,2)),c(rep(0,5),rep(1,2)),c(rep(0,6),rep(1,2)),c(rep(0,4),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,6)),c(rep(0,6),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,3),rep(1,0)),c(rep(0,1),rep(1,0)))
Corbula.Po3.M12.5cm.depth.for.pyrite=c(rep(4,0),rep(8,4),rep(12,8),rep(16,24),rep(20,10),rep(25,9),rep(30,6),rep(35,12),rep(40,3),rep(45,7),rep(50,10),rep(55,17),rep(60,29),rep(65,12),rep(70,18),rep(75,11),rep(80,44),rep(85,20),rep(90,21),rep(95,9),rep(100,7),rep(105,8),rep(110,4),rep(115,0),rep(120,7),rep(125,6),rep(130,0),rep(135,1),rep(140,1),rep(145,3),rep(150,1))
Corbula.Po3.M12.5cm.for.pyrite=c(c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,2)),c(rep(0,8),rep(1,0)),c(rep(0,15),rep(1,9)),c(rep(0,8),rep(1,2)),c(rep(0,6),rep(1,3)),c(rep(0,1),rep(1,5)),c(rep(0,4),rep(1,8)),c(rep(0,1),rep(1,2)),c(rep(0,2),rep(1,5)),c(rep(0,5),rep(1,5)),c(rep(0,9),rep(1,8)),c(rep(0,5),rep(1,24)),c(rep(0,5),rep(1,7)),c(rep(0,10),rep(1,8)),c(rep(0,6),rep(1,5)),c(rep(0,20),rep(1,24)),c(rep(0,8),rep(1,12)),c(rep(0,11),rep(1,10)),c(rep(0,7),rep(1,2)),c(rep(0,5),rep(1,2)),c(rep(0,6),rep(1,2)),c(rep(0,4),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,6)),c(rep(0,6),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,3),rep(1,0)),c(rep(0,1),rep(1,0)))

#Po 3 - M14 - 2 cm and 5 cm depth 
Corbula.Po3.M14.add.fine.depth.for.pyrite=c(rep(2,0),rep(4,1),rep(6,4),rep(8,5),rep(10,11),rep(12,6),rep(14,8),rep(16,17),rep(18,17),rep(20,12),rep(25,23),rep(30,14),rep(35,32),rep(40,17),rep(45,48),rep(50,25),rep(55,25),rep(60,22),rep(65,10),rep(70,8),rep(75,36),rep(80,9),rep(85,15),rep(90,29),rep(95,16),rep(100,17),rep(105,2),rep(110,12),rep(115,4),rep(120,0),rep(125,2),rep(130,1),rep(135,6),rep(140,0),rep(145,4),rep(150,1))
Corbula.Po3.M14.add.for.pyrite=c(c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,3),rep(1,1)),c(rep(0,4),rep(1,1)),c(rep(0,10),rep(1,1)),c(rep(0,6),rep(1,0)),c(rep(0,7),rep(1,1)),c(rep(0,12),rep(1,5)),c(rep(0,9),rep(1,8)),c(rep(0,9),rep(1,3)),c(rep(0,13),rep(1,10)),c(rep(0,8),rep(1,6)),c(rep(0,13),rep(1,19)),c(rep(0,9),rep(1,8)),c(rep(0,30),rep(1,18)),c(rep(0,14),rep(1,11)),c(rep(0,10),rep(1,15)),c(rep(0,9),rep(1,13)),c(rep(0,5),rep(1,5)),c(rep(0,3),rep(1,5)),c(rep(0,18),rep(1,18)),c(rep(0,3),rep(1,6)),c(rep(0,6),rep(1,9)),c(rep(0,22),rep(1,7)),c(rep(0,13),rep(1,3)),c(rep(0,14),rep(1,3)),c(rep(0,2),rep(1,0)),c(rep(0,11),rep(1,1)),c(rep(0,2),rep(1,2)),c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,6),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,4),rep(1,0)),c(rep(0,1),rep(1,0)))
Corbula.Po3.M14.5cm.depth.for.pyrite=c(rep(4,1),rep(8,9),rep(12,17),rep(16,25),rep(20,29),rep(25,23),rep(30,14),rep(35,32),rep(40,17),rep(45,48),rep(50,25),rep(55,25),rep(60,22),rep(65,10),rep(70,8),rep(75,36),rep(80,9),rep(85,15),rep(90,29),rep(95,16),rep(100,17),rep(105,2),rep(110,12),rep(115,4),rep(120,0),rep(125,2),rep(130,1),rep(135,6),rep(140,0),rep(145,4),rep(150,1))
Corbula.Po3.M14.5cm.for.pyrite=c(c(rep(0,1),rep(1,0)),c(rep(0,7),rep(1,2)),c(rep(0,16),rep(1,1)),c(rep(0,19),rep(1,6)),c(rep(0,18),rep(1,11)),c(rep(0,13),rep(1,10)),c(rep(0,8),rep(1,6)),c(rep(0,13),rep(1,19)),c(rep(0,9),rep(1,8)),c(rep(0,30),rep(1,18)),c(rep(0,14),rep(1,11)),c(rep(0,10),rep(1,15)),c(rep(0,9),rep(1,13)),c(rep(0,5),rep(1,5)),c(rep(0,3),rep(1,5)),c(rep(0,18),rep(1,18)),c(rep(0,3),rep(1,6)),c(rep(0,6),rep(1,9)),c(rep(0,22),rep(1,7)),c(rep(0,13),rep(1,3)),c(rep(0,14),rep(1,3)),c(rep(0,2),rep(1,0)),c(rep(0,11),rep(1,1)),c(rep(0,2),rep(1,2)),c(rep(0,0),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,6),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,4),rep(1,0)),c(rep(0,1),rep(1,0)))

#2 cm and 5 cm depth 
Corbula.Po3.M13.add.fine.depth.for.pyrite=c(rep(2,0),rep(4,0),rep(6,0),rep(8,0),rep(10,1),rep(12,0),rep(14,0),rep(16,0),rep(18,3),rep(20,3),rep(25,17),rep(30,56),rep(35,82),rep(40,25),rep(45,0),rep(50,31),rep(55,30),rep(60,35),rep(65,47),rep(70,24),rep(75,81),rep(80,24),rep(85,0),rep(90,11),rep(95,16),rep(100,17),rep(105,17),rep(110,1),rep(115,0),rep(120,13),rep(125,10),rep(130,1),rep(135,0),rep(140,0),rep(145,4),rep(150,4),rep(155,4))
Corbula.Po3.M13.add.for.pyrite=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,0),rep(1,0)),c(rep(0,3),rep(1,0)),c(rep(0,2),rep(1,1)),c(rep(0,13),rep(1,4)),c(rep(0,48),rep(1,8)),c(rep(0,66),rep(1,16)),c(rep(0,15),rep(1,10)),c(rep(0,0),rep(1,0)),c(rep(0,14),rep(1,17)),c(rep(0,13),rep(1,17)),c(rep(0,25),rep(1,10)),c(rep(0,24),rep(1,23)),c(rep(0,15),rep(1,9)),c(rep(0,60),rep(1,21)),c(rep(0,21),rep(1,3)),c(rep(0,0),rep(1,0)),c(rep(0,10),rep(1,1)),c(rep(0,14),rep(1,2)),c(rep(0,16),rep(1,1)),c(rep(0,14),rep(1,3)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,12),rep(1,1)),c(rep(0,7),rep(1,3)),c(rep(0,0),rep(1,1)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,4),rep(1,0)),c(rep(0,3),rep(1,1)),c(rep(0,2),rep(1,2)))
Corbula.Po3.M13.5cm.depth.for.pyrite=c(rep(4,0),rep(8,0),rep(12,1),rep(16,0),rep(20,6),rep(25,17),rep(30,56),rep(35,82),rep(40,25),rep(45,0),rep(50,31),rep(55,30),rep(60,35),rep(65,47),rep(70,24),rep(75,81),rep(80,24),rep(85,0),rep(90,11),rep(95,16),rep(100,17),rep(105,17),rep(110,1),rep(115,0),rep(120,13),rep(125,10),rep(130,1),rep(135,0),rep(140,0),rep(145,4),rep(150,4),rep(155,4))
Corbula.Po3.M13.5cm.for.pyrite=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,5),rep(1,1)),c(rep(0,13),rep(1,4)),c(rep(0,48),rep(1,8)),c(rep(0,66),rep(1,16)),c(rep(0,15),rep(1,10)),c(rep(0,0),rep(1,0)),c(rep(0,14),rep(1,17)),c(rep(0,13),rep(1,17)),c(rep(0,25),rep(1,10)),c(rep(0,24),rep(1,23)),c(rep(0,15),rep(1,9)),c(rep(0,60),rep(1,21)),c(rep(0,21),rep(1,3)),c(rep(0,0),rep(1,0)),c(rep(0,10),rep(1,1)),c(rep(0,14),rep(1,2)),c(rep(0,16),rep(1,1)),c(rep(0,14),rep(1,3)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,12),rep(1,1)),c(rep(0,7),rep(1,3)),c(rep(0,0),rep(1,1)),c(rep(0,0),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,4),rep(1,0)),c(rep(0,3),rep(1,1)),c(rep(0,2),rep(1,2)))

#2 cm and 5 cm depth Po4 M20
Corbula.Po4.M20.add.fine.depth.for.pyrite=c(rep(2,1),rep(4,1),rep(6,2),rep(8,22),rep(10,1),rep(12,0),rep(14,1),rep(16,9),rep(18,8),rep(20,9),rep(25,24),rep(30,20),rep(35,40),rep(40,43),rep(45,16),rep(50,68),rep(55,66),rep(60,8),rep(65,78),rep(70,48),rep(75,38),rep(80,104),rep(85,57),rep(90,17),rep(95,23),rep(100,10),rep(105,24),rep(110,4),rep(115,13),rep(120,15),rep(125,10),rep(130,12),rep(135,6),rep(140,2),rep(145,7),rep(150,1))
Corbula.Po4.M20.add.for.pyrite=c(c(rep(0,1),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,2),rep(1,0)),c(rep(0,19),rep(1,3)),c(rep(0,1),rep(1,0)),c(rep(0,0),rep(1,0)),c(rep(0,1),rep(1,0)),c(rep(0,9),rep(1,0)),c(rep(0,8),rep(1,0)),c(rep(0,8),rep(1,1)),c(rep(0,17),rep(1,7)),c(rep(0,14),rep(1,6)),c(rep(0,28),rep(1,12)),c(rep(0,33),rep(1,10)),c(rep(0,12),rep(1,4)),c(rep(0,53),rep(1,15)),c(rep(0,49),rep(1,17)),c(rep(0,4),rep(1,4)),c(rep(0,61),rep(1,17)),c(rep(0,37),rep(1,11)),c(rep(0,28),rep(1,10)),c(rep(0,70),rep(1,34)),c(rep(0,45),rep(1,12)),c(rep(0,13),rep(1,4)),c(rep(0,21),rep(1,2)),c(rep(0,8),rep(1,2)),c(rep(0,23),rep(1,1)),c(rep(0,4),rep(1,0)),c(rep(0,12),rep(1,1)),c(rep(0,12),rep(1,3)),c(rep(0,8),rep(1,2)),c(rep(0,10),rep(1,2)),c(rep(0,6),rep(1,0)),c(rep(0,1),rep(1,1)),c(rep(0,4),rep(1,3)),c(rep(0,1),rep(1,0)))
Corbula.Po4.M20.5cm.depth.for.pyrite=c(rep(4,2),rep(8,24),rep(12,1),rep(16,10),rep(20,17),rep(25,24),rep(30,20),rep(35,40),rep(40,43),rep(45,16),rep(50,68),rep(55,66),rep(60,8),rep(65,78),rep(70,48),rep(75,38),rep(80,104),rep(85,57),rep(90,17),rep(95,23),rep(100,10),rep(105,24),rep(110,4),rep(115,13),rep(120,15),rep(125,10),rep(130,12),rep(135,6),rep(140,2),rep(145,7),rep(150,1))
Corbula.Po4.M20.5cm.for.pyrite=c(c(rep(0,2),rep(1,0)),c(rep(0,21),rep(1,3)),c(rep(0,1),rep(1,0)),c(rep(0,10),rep(1,0)),c(rep(0,16),rep(1,1)),c(rep(0,17),rep(1,7)),c(rep(0,14),rep(1,6)),c(rep(0,28),rep(1,12)),c(rep(0,33),rep(1,10)),c(rep(0,12),rep(1,4)),c(rep(0,53),rep(1,15)),c(rep(0,49),rep(1,17)),c(rep(0,4),rep(1,4)),c(rep(0,61),rep(1,17)),c(rep(0,37),rep(1,11)),c(rep(0,28),rep(1,10)),c(rep(0,70),rep(1,34)),c(rep(0,45),rep(1,12)),c(rep(0,13),rep(1,4)),c(rep(0,21),rep(1,2)),c(rep(0,8),rep(1,2)),c(rep(0,23),rep(1,1)),c(rep(0,4),rep(1,0)),c(rep(0,12),rep(1,1)),c(rep(0,12),rep(1,3)),c(rep(0,8),rep(1,2)),c(rep(0,10),rep(1,2)),c(rep(0,6),rep(1,0)),c(rep(0,1),rep(1,1)),c(rep(0,4),rep(1,3)),c(rep(0,1),rep(1,0)))

#2 cm and 5 cm depth Po4 M21
Corbula.Po4.M21.add.fine.depth.for.pyrite=c(rep(2,0),rep(4,0),rep(6,0),rep(8,2),rep(10,0),rep(12,0),rep(14,0),rep(16,0),rep(18,0),rep(20,4),rep(25,57),rep(30,111),rep(35,28),rep(40,31),rep(45,11),rep(50,36),rep(55,135),rep(60,19),rep(65,14),rep(70,30),rep(75,54),rep(80,12),rep(85,12),rep(90,6),rep(95,6),rep(100,12),rep(105,8),rep(110,6),rep(115,10),rep(120,16),rep(125,11),rep(130,1),rep(135,0),rep(140,0),rep(145,2),rep(150,2))
Corbula.Po4.M21.add.for.pyrite=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,2),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,4),rep(1,0)),c(rep(0,40),rep(1,17)),c(rep(0,96),rep(1,15)),c(rep(0,20),rep(1,8)),c(rep(0,22),rep(1,9)),c(rep(0,4),rep(1,7)),c(rep(0,28),rep(1,8)),c(rep(0,86),rep(1,49)),c(rep(0,13),rep(1,6)),c(rep(0,7),rep(1,7)),c(rep(0,27),rep(1,3)),c(rep(0,28),rep(1,26)),c(rep(0,10),rep(1,2)),c(rep(0,10),rep(1,2)),c(rep(0,5),rep(1,1)),c(rep(0,5),rep(1,1)),c(rep(0,11),rep(1,1)),c(rep(0,8),rep(1,0)),c(rep(0,5),rep(1,1)),c(rep(0,8),rep(1,2)),c(rep(0,13),rep(1,3)),c(rep(0,9),rep(1,2)),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,2),rep(1,0)))
Corbula.Po4.M21.5cm.depth.for.pyrite=c(rep(4,0),rep(8,2),rep(12,0),rep(16,0),rep(20,4),rep(25,57),rep(30,111),rep(35,28),rep(40,31),rep(45,11),rep(50,36),rep(55,135),rep(60,19),rep(65,14),rep(70,30),rep(75,54),rep(80,12),rep(85,12),rep(90,6),rep(95,6),rep(100,12),rep(105,8),rep(110,6),rep(115,10),rep(120,16),rep(125,11),rep(130,1),rep(135,0),rep(140,0),rep(145,2),rep(150,2))
Corbula.Po4.M21.5cm.for.pyrite=c(c(rep(0,0),rep(1,0)),c(rep(0,2),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,40),rep(1,17)),c(rep(0,96),rep(1,15)),c(rep(0,20),rep(1,8)),c(rep(0,22),rep(1,9)),c(rep(0,4),rep(1,7)),c(rep(0,28),rep(1,8)),c(rep(0,86),rep(1,49)),c(rep(0,13),rep(1,6)),c(rep(0,7),rep(1,7)),c(rep(0,27),rep(1,3)),c(rep(0,28),rep(1,26)),c(rep(0,10),rep(1,2)),c(rep(0,10),rep(1,2)),c(rep(0,5),rep(1,1)),c(rep(0,5),rep(1,1)),c(rep(0,11),rep(1,1)),c(rep(0,8),rep(1,0)),c(rep(0,5),rep(1,1)),c(rep(0,8),rep(1,2)),c(rep(0,13),rep(1,3)),c(rep(0,9),rep(1,2)),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,2),rep(1,0)))

#Panzano M28 - 2 cm and 5 cm depth 
Corbula.Panzano.M28.add.fine.depth.for.pyrite=c(rep(2,2),rep(4,12),rep(6,58),rep(8,23),rep(10,47),rep(12,23),rep(14,51),rep(16,47),rep(18,4),rep(20,18),rep(25,41),rep(30,74),rep(35,37),rep(40,102),rep(45,207),rep(50,42),rep(55,40),rep(60,32),rep(65,10),rep(70,12),rep(75,30),rep(80,28),rep(85,13),rep(90,2),rep(95,7),rep(100,31),rep(105,38),rep(110,45),rep(115,14),rep(120,38),rep(125,36),rep(130,18),rep(135,68),rep(140,45),rep(145,26),rep(150,5))
Corbula.Panzano.M28.add.for.pyrite=c(c(rep(0,2),rep(1,0)),c(rep(0,12),rep(1,0)),c(rep(0,38),rep(1,20)),c(rep(0,8),rep(1,15)),c(rep(0,19),rep(1,28)),c(rep(0,13),rep(1,10)),c(rep(0,29),rep(1,22)),c(rep(0,23),rep(1,24)),c(rep(0,3),rep(1,1)),c(rep(0,9),rep(1,9)),c(rep(0,22),rep(1,19)),c(rep(0,46),rep(1,28)),c(rep(0,19),rep(1,18)),c(rep(0,58),rep(1,44)),c(rep(0,137),rep(1,70)),c(rep(0,32),rep(1,10)),c(rep(0,31),rep(1,9)),c(rep(0,19),rep(1,13)),c(rep(0,7),rep(1,3)),c(rep(0,8),rep(1,4)),c(rep(0,19),rep(1,11)),c(rep(0,18),rep(1,10)),c(rep(0,8),rep(1,5)),c(rep(0,2),rep(1,0)),c(rep(0,5),rep(1,2)),c(rep(0,22),rep(1,9)),c(rep(0,29),rep(1,9)),c(rep(0,30),rep(1,15)),c(rep(0,11),rep(1,3)),c(rep(0,26),rep(1,12)),c(rep(0,25),rep(1,11)),c(rep(0,12),rep(1,6)),c(rep(0,48),rep(1,20)),c(rep(0,29),rep(1,16)),c(rep(0,17),rep(1,9)),c(rep(0,4),rep(1,1)))
Corbula.Panzano.M28.5cm.depth.for.pyrite=c(rep(4,14),rep(8,81),rep(12,70),rep(16,98),rep(20,22),rep(25,41),rep(30,74),rep(35,37),rep(40,102),rep(45,207),rep(50,42),rep(55,40),rep(60,32),rep(65,10),rep(70,12),rep(75,30),rep(80,28),rep(85,13),rep(90,2),rep(95,7),rep(100,31),rep(105,38),rep(110,45),rep(115,14),rep(120,38),rep(125,36),rep(130,18),rep(135,68),rep(140,45),rep(145,26),rep(150,5))
Corbula.Panzano.M28.5cm.for.pyrite=c(c(rep(0,14),rep(1,0)),c(rep(0,46),rep(1,35)),c(rep(0,32),rep(1,38)),c(rep(0,52),rep(1,46)),c(rep(0,12),rep(1,10)),c(rep(0,22),rep(1,19)),c(rep(0,46),rep(1,28)),c(rep(0,19),rep(1,18)),c(rep(0,58),rep(1,44)),c(rep(0,137),rep(1,70)),c(rep(0,32),rep(1,10)),c(rep(0,31),rep(1,9)),c(rep(0,19),rep(1,13)),c(rep(0,7),rep(1,3)),c(rep(0,8),rep(1,4)),c(rep(0,19),rep(1,11)),c(rep(0,18),rep(1,10)),c(rep(0,8),rep(1,5)),c(rep(0,2),rep(1,0)),c(rep(0,5),rep(1,2)),c(rep(0,22),rep(1,9)),c(rep(0,29),rep(1,9)),c(rep(0,30),rep(1,15)),c(rep(0,11),rep(1,3)),c(rep(0,26),rep(1,12)),c(rep(0,25),rep(1,11)),c(rep(0,12),rep(1,6)),c(rep(0,48),rep(1,20)),c(rep(0,29),rep(1,16)),c(rep(0,17),rep(1,9)),c(rep(0,4),rep(1,1)))

#pyrite PanzanoM29
Corbula.Panzano.M29.add.for.pyrite=c(c(rep(0,42),rep(1,6)),c(rep(0,28),rep(1,3)),c(rep(0,41),rep(1,7)),c(rep(0,17),rep(1,12)),c(rep(0,18),rep(1,10)),c(rep(0,24),rep(1,19)),c(rep(0,45),rep(1,24)),c(rep(0,38),rep(1,29)),c(rep(0,32),rep(1,15)),c(rep(0,17),rep(1,15)),c(rep(0,36),rep(1,26)),c(rep(0,35),rep(1,25)),c(rep(0,40),rep(1,17)),c(rep(0,32),rep(1,19)),c(rep(0,30),rep(1,8)),c(rep(0,31),rep(1,13)),c(rep(0,24),rep(1,6)),c(rep(0,18),rep(1,6)),c(rep(0,18),rep(1,21)),c(rep(0,15),rep(1,7)),c(rep(0,9),rep(1,7)),c(rep(0,15),rep(1,7)),c(rep(0,3),rep(1,2)),c(rep(0,13),rep(1,8)),c(rep(0,10),rep(1,6)),c(rep(0,21),rep(1,10)),c(rep(0,45),rep(1,27)),c(rep(0,28),rep(1,12)),c(rep(0,23),rep(1,7)),c(rep(0,19),rep(1,11)),c(rep(0,29),rep(1,10)),c(rep(0,24),rep(1,14)),c(rep(0,45),rep(1,18)),c(rep(0,27),rep(1,14)),c(rep(0,22),rep(1,3)),c(rep(0,23),rep(1,14)))
Corbula.Panzano.M29.add.fine.depth.for.pyrite=c(rep(2,48),rep(4,31),rep(6,48),rep(8,29),rep(10,28),rep(12,43),rep(14,69),rep(16,67),rep(18,47),rep(20,32),rep(25,62),rep(30,60),rep(35,57),rep(40,51),rep(45,38),rep(50,44),rep(55,30),rep(60,24),rep(65,39),rep(70,22),rep(75,16),rep(80,22),rep(85,5),rep(90,21),rep(95,16),rep(100,31),rep(105,72),rep(110,40),rep(115,30),rep(120,30),rep(125,39),rep(130,38),rep(135,63),rep(140,41),rep(145,25),rep(150,37))
Corbula.Panzano.M29.5cm.for.pyrite=c(c(rep(0,70),rep(1,9)),c(rep(0,58),rep(1,19)),c(rep(0,42),rep(1,29)),c(rep(0,83),rep(1,53)),c(rep(0,49),rep(1,30)),c(rep(0,36),rep(1,26)),c(rep(0,35),rep(1,25)),c(rep(0,40),rep(1,17)),c(rep(0,32),rep(1,19)),c(rep(0,30),rep(1,8)),c(rep(0,31),rep(1,13)),c(rep(0,24),rep(1,6)),c(rep(0,18),rep(1,6)),c(rep(0,18),rep(1,21)),c(rep(0,15),rep(1,7)),c(rep(0,9),rep(1,7)),c(rep(0,15),rep(1,7)),c(rep(0,3),rep(1,2)),c(rep(0,13),rep(1,8)),c(rep(0,10),rep(1,6)),c(rep(0,21),rep(1,10)),c(rep(0,45),rep(1,27)),c(rep(0,28),rep(1,12)),c(rep(0,23),rep(1,7)),c(rep(0,19),rep(1,11)),c(rep(0,29),rep(1,10)),c(rep(0,24),rep(1,14)),c(rep(0,45),rep(1,18)),c(rep(0,27),rep(1,14)),c(rep(0,22),rep(1,3)),c(rep(0,23),rep(1,14)))
Corbula.Panzano.M29.5cm.depth.for.pyrite=c(rep(4,79),rep(8,77),rep(12,71),rep(16,136),rep(20,79),rep(25,62),rep(30,60),rep(35,57),rep(40,51),rep(45,38),rep(50,44),rep(55,30),rep(60,24),rep(65,39),rep(70,22),rep(75,16),rep(80,22),rep(85,5),rep(90,21),rep(95,16),rep(100,31),rep(105,72),rep(110,40),rep(115,30),rep(120,30),rep(125,39),rep(130,38),rep(135,63),rep(140,41),rep(145,25),rep(150,37))

########################################
#PYRITE PERCENTAGE IN CORBULA
########################################
all.Corbula.Panzano.M29.finecm=c(Corbula.Panzano.M29.add.fine.depth.for.pyrite)
all.Corbula.Panzano.M29.pyrite=c(Corbula.Panzano.M29.add.for.pyrite)
all.Corbula.Panzano.M28.finecm=c(Corbula.Panzano.finecm,Corbula.Panzano.M28.add.fine.depth.for.pyrite)
all.Corbula.Panzano.M28.pyrite=c(Corbula.Panzano.pyrite,Corbula.Panzano.M28.add.for.pyrite)
all.Corbula.Po3.finecm=c(Corbula.Po3.finecm,Corbula.Po3.M13.add.fine.depth.for.pyrite,Corbula.Po3.M14.add.fine.depth.for.pyrite,Corbula.Po3.M12.add.fine.depth.for.pyrite)
all.Corbula.Po3.pyrite=c(Corbula.Po3.pyrite,Corbula.Po3.M13.add.for.pyrite,Corbula.Po3.M14.add.for.pyrite,Corbula.Po3.M12.add.for.pyrite)
all.Corbula.Po4.finecm=c(Corbula.Po4.finecm,Corbula.Po4.M20.add.fine.depth.for.pyrite,Corbula.Po4.M21.add.fine.depth.for.pyrite)
all.Corbula.Po4.pyrite=c(Corbula.Po4.pyrite,Corbula.Po4.M20.add.for.pyrite, Corbula.Po4.M21.add.for.pyrite)
all.Corbula.Panzano.M29.5cm=c(Corbula.Panzano.M29.5cm.depth.for.pyrite)
all.Corbula.Panzano.M29.5cm.pyrite=c(Corbula.Panzano.M29.5cm.for.pyrite)
all.Corbula.Panzano.M28.5cm=c(Corbula.Panzano.5cm,Corbula.Panzano.M28.5cm.depth.for.pyrite)
all.Corbula.Panzano.M28.5cm.pyrite=c(Corbula.Panzano.pyrite,Corbula.Panzano.M28.5cm.for.pyrite)
all.Corbula.Po3.5cm=c(Corbula.Po3.5cm,Corbula.Po3.M13.5cm.depth.for.pyrite,Corbula.Po3.M14.5cm.depth.for.pyrite,Corbula.Po3.M12.5cm.depth.for.pyrite)
all.Corbula.Po3.5cm.pyrite=c(Corbula.Po3.pyrite,Corbula.Po3.M13.5cm.for.pyrite,Corbula.Po3.M14.5cm.for.pyrite,Corbula.Po3.M12.5cm.for.pyrite)
all.Corbula.Po4.5cm=c(Corbula.Po4.5cm,Corbula.Po4.M20.5cm.depth.for.pyrite,Corbula.Po4.M21.5cm.depth.for.pyrite)
all.Corbula.Po4.5cm.pyrite=c(Corbula.Po4.pyrite,Corbula.Po4.M20.5cm.for.pyrite, Corbula.Po4.M21.5cm.for.pyrite)
all.Corbula.Brijuni.5cm=c(Corbula.M44.5cm)
all.Corbula.Brijuni.5cm.pyrite=c(Corbula.M44.pyrite)
all.Corbula.Piran2.5cm=c(Corbula.Piran.5cm)
all.Corbula.Piran2.5cm.pyrite=c(Corbula.Piran.pyrite)
all.Corbula.Panzano.5cm=c(Corbula.Panzano.5cm,Corbula.Panzano.M28.5cm.depth.for.pyrite,Corbula.Panzano.M29.5cm.depth.for.pyrite)
all.Corbula.Panzano.5cm.pyrite=c(Corbula.Panzano.pyrite,Corbula.Panzano.M28.5cm.for.pyrite,Corbula.Panzano.M29.5cm.for.pyrite)

#sample sizes for per-increment percentages
pyrite.N.Piran2.5cm=tapply(all.Corbula.Piran2.5cm.pyrite,all.Corbula.Piran2.5cm, length)
pyrite.N.Brijuni.5cm=tapply(all.Corbula.Brijuni.5cm.pyrite,all.Corbula.Brijuni.5cm, length)
pyrite.N.Po3.5cm=tapply(all.Corbula.Po3.5cm.pyrite,all.Corbula.Po3.5cm, length)
pyrite.N.Po4.5cm=tapply(all.Corbula.Po4.5cm.pyrite,all.Corbula.Po4.5cm, length)
pyrite.N.Panzano.M28.5cm=tapply(all.Corbula.Panzano.M28.5cm.pyrite,all.Corbula.Panzano.M28.5cm, length)
pyrite.N.Panzano.M29.5cm=tapply(all.Corbula.Panzano.M29.5cm.pyrite,all.Corbula.Panzano.M29.5cm, length)
pyrite.N.Panzano.5cm=tapply(all.Corbula.Panzano.5cm.pyrite,all.Corbula.Panzano.5cm, length)

#depths for pyrite percentages
cor.depths.Po3.5cm=as.numeric(names(tapply(all.Corbula.Po3.5cm.pyrite,all.Corbula.Po3.5cm, mean)))
cor.depths.Po4.5cm=as.numeric(names(tapply(all.Corbula.Po4.5cm.pyrite,all.Corbula.Po4.5cm, mean)))
cor.depths.Panzano.M28.5cm=as.numeric(names(tapply(all.Corbula.Panzano.M28.5cm.pyrite,all.Corbula.Panzano.M28.5cm, mean)))
cor.depths.Panzano.M29.5cm=as.numeric(names(tapply(all.Corbula.Panzano.M29.5cm.pyrite,all.Corbula.Panzano.M29.5cm, mean)))
cor.depths.Brijuni.5cm=as.numeric(names(tapply(all.Corbula.Brijuni.5cm.pyrite,all.Corbula.Brijuni.5cm, mean)))
cor.depths.Piran2.5cm=as.numeric(names(tapply(all.Corbula.Piran2.5cm.pyrite,all.Corbula.Piran2.5cm, mean)))
cor.depths.Panzano.5cm=as.numeric(names(tapply(all.Corbula.Panzano.5cm.pyrite,all.Corbula.Panzano.5cm, mean)))

##########################################################################################
#COVARIATES FOR PYRITE PRESERVATION
##########################################################################################
#match TOC levels with Corbula levels
#not all increments at Piran and Brijuni were scored for pyrite in Corbula
#for Brijuni
temp1=c(4,8,12,16,20,seq(25,155,by=5))
temp2=c(4,8,12,16,20,seq(30,150,by=10),155)
Brijuni.ID=match(temp2, temp1)
#for Piran
temp1=c(4,8,12,16,20,seq(25,155,by=5))
temp2=c(4,12,30,50,70,90,105,130,150)
Piran2.ID=match(temp2, temp1)

##############################################################################
#CORBULA PERCENTAGES
##############################################################################
Corbula.percent.Po3=c(0.1666667, 0.4482759, 0.5423729, 0.5316456, 0.5657895, 0.6543210, 0.5822785, 0.6785714, 0.5000000, 0.5140187, 0.5066667, 0.6349206, 0.6125000, 0.5000000, 0.5972222, 0.6915888, 0.5000000, 0.6545455, 0.3978495, 0.4150943, 0.4343434, 0.3392857, 0.3333333, 0.2142857, 0.2653061, 0.2444444, 0.1250000, 0.3142857, 0.2727273, 0.2580645, 0.2272727, 0.1333333)
Corbula.percent.Po4=c(0.3750000, 0.5633803, 0.4444444, 0.5555556, 0.4693878, 0.6557377, 0.6463415, 0.6475410, 0.5090909, 0.5294118, 0.6375000, 0.5432692, 0.3733333, 0.5076142, 0.5213675, 0.5639098, 0.5659341, 0.4607843, 0.2542373, 0.3402778, 0.3773585, 0.2906977, 0.4142857, 0.3275862, 0.3636364, 0.3409091, 0.2448980, 0.2708333, 0.1500000, 0.1875000, 0.1363636, 0.3684211) 
Corbula.percent.Panzano.M28=c(0.17129630,  0.27237354,  0.28372093,  0.16615385,  0.14227642,  0.08992806,  0.13636364,  0.13194444,  0.18237082,  0.31593407, 0.20000000,  0.16216216, 0.13636364,  0.15231788,  0.22950820,  0.15315315,  0.18000000,  0.10256410,  0.07258065, 0.10559006,0.12318841,  0.15432099,  0.18055556,  0.24193548,  0.25225225,  0.21359223,  0.27551020,  0.26060606,  0.12935323,  0.11016949, 0.14583333)
Corbula.percent.Panzano.M29=c(0.26436782,  0.20796460,  0.24074074,  0.25000000,  0.13056380,  0.12639405,  0.13718412,  0.18848168,  0.16000000,  0.12371134, 0.12558140,  0.13422819,  0.12380952,  0.19230769,  0.27272727,  0.18518519,  0.25925926,  0.05660377,  0.11111111,  0.08333333, 0.19298246,  0.21808511,  0.14379085,  0.18446602,  0.16964286,  0.23000000,  0.29577465,  0.28448276,  0.17910448,  0.13000000, 0.16788321)
Corbula.percent.Piran2=c(0.08706625, 0.07420495, 0.06207882, 0.05806452,0.06095041, 0.04655674,0.07175573, 0.07046980,0.08407517,0.08320726,0.08186687,0.08533654,0.11649215,0.11608222,0.14942529, 0.16566265,0.16094987,0.12903226,0.12852665,0.13058419,0.14347826,0.17508418,0.16172507,0.20416667,0.21428571,0.14487633,0.21010638,0.13938053,0.16018307,0.14006515,0.19591837,0.17391304) 
Corbula.percent.Brijuni=c(0.032258065,0.056250000,0.046979866,0.050156740,0.029629630,0.013651877,0.027842227,0.025316456,0.025553663,0.045908184,0.033755274,0.033600000,0.033395176,0.031446541,0.027737226,0.034768212,0.030395137,0.030588235,0.016309888,0.024504084,0.014023732,0.026796590,0.009922822,0.020227560,0.012141280,0.020300088,0.007513148,0.015205725,0.011398964,0.024509804,0.021052632,0.022950820,0.008264463) 
all.Corbula.percent=c(Corbula.percent.Po3[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20],Corbula.percent.Po4[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20],Corbula.percent.Panzano.M28[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>6],Corbula.percent.Panzano.M29[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>6],Corbula.percent.Piran2[Piran2.ID][pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>8],Corbula.percent.Brijuni[Brijuni.ID][pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20])
all.Corbula.percent.for.pyrite=all.Corbula.percent

##########################
#MEAN SIZE OF CORBULA
##########################
mean.size.Brijuni=c(4.538029,4.232404,4.361662,4.895595,4.782005,3.761450,3.935959,4.161395,3.802079,3.620235,3.060489,3.280399,3.889179,3.668294,3.101922,3.730864,3.338956,3.337209,3.530185,4.175014,3.763468,3.642849,4.057574,3.880884,3.556666,3.939638,4.163679,4.937091,4.013974,4.196133,5.667017,3.630503)
mean.size.Po3=c(3.725907,4.771816,4.794191,6.021267,6.126549,6.904323,7.933230,9.142016,6.535067,6.989911,8.577965,8.731000,7.432732,8.376845,9.025681,7.089800,5.817659,6.087129,5.308381,3.865567,4.310520,3.617360,3.289417,2.938245,6.054400,5.168333,3.944000,3.147706,3.024557,2.371333,5.377500,2.870000)
mean.size.Po4=c(6.123876,5.253910,6.742328,6.268744,5.195384,7.301230,5.836574,8.428187,8.441882,9.344864,8.423565,7.394548,7.044735,6.431091,5.016429,3.784605,4.428767,3.117654,3.379500,3.401499,4.120286,3.314429,3.706968,3.781750,3.988181,3.949000,2.620024,3.193636,2.786000,3.352500,2.215783,4.188477)
mean.size.Panzano.M28=c(4.278590,5.530471,5.724320,5.591535,3.569322,6.459666,5.439440,4.408201,3.653616,3.364687,3.342619,2.908874,3.835034,3.848876,3.868926,4.062944,3.743138,3.667982,4.022748,4.056716,2.990933,3.400659,3.358591,3.341236,2.965961,3.462900,3.473914,3.588521,4.379905,3.759660,3.205527)
mean.size.Panzano.M29=c(4.627855,5.303088,8.050561,7.484742,4.960588,4.298009,5.045937,3.951358,4.130932,3.437200,4.913526,3.219603,3.940383,3.844825,3.715164,4.086737,3.472852,4.653793,3.865996,4.613065,3.334466,3.751045,4.498967,3.172502,3.909852,4.126252,4.126252,3.287393,4.077851,4.157685,4.938948)
mean.size.Piran2=c(3.543109,3.761236,3.640591,4.168434,3.560846,3.651127,3.443406,3.729228,3.449448,4.175324,3.424421,3.571916,3.701388,3.264221,3.307315,3.367719,3.297307,3.460639,3.182368,3.741585,3.435437,4.173260,3.754391,4.757026,4.049901,4.178691,4.750240,4.547396,4.636377,3.425787,3.376536)
all.Corbula.mean.size=c(mean.size.Po3[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20],mean.size.Po4[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20],
mean.size.Panzano.M28[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>6],mean.size.Panzano.M29[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>6],
mean.size.Brijuni[Brijuni.ID][pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20],mean.size.Piran2[Piran2.ID][pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>8])

##########################
#MEDIAN SIZE OF CORBULA
##########################
Po4.q50=c(4.640000, 3.910000, 7.342500, 7.490000, 4.551000, 7.995000, 4.981000, 9.995000, 9.930000, 9.700000, 8.880000, 7.870000,7.318000, 6.745000, 3.668500, 3.651000, 4.108000, 2.836000, 3.047000, 3.278811, 3.095000, 3.177500, 3.390000, 2.974500,2.850000, 3.284500, 2.239000, 2.848237, 2.414000, 3.352500, 2.215783, 3.390000)
Po3.q50=c(3.719792,  4.078353,  4.195500,  5.631213,  5.889000,  7.338000,  9.460000, 10.042000,  6.115000,  7.038000,  9.290000,10.110000,  8.35600,  9.627568,  9.557500,  7.930250,  5.609000,  5.037500,  4.123000,  4.208000,  3.705000,  2.910000,3.108726,  2.963000,  5.793000,  4.949000,  2.000000,  3.013793,  2.700000,  2.329000,  5.377500,  2.000000)
Panzano.M28.q50=c(3.508693, 4.126082, 3.998956, 3.750000, 3.137773, 3.890000, 3.707821, 3.120000, 3.237944, 3.288002, 2.945055, 2.685209, 2.987109, 2.930000, 2.750000, 3.290000, 3.328042, 3.059414, 3.588044, 3.460000, 2.757761, 2.710000, 2.970000, 2.938874, 2.619956, 2.900000, 2.769847, 3.107587, 2.975053, 3.160000, 2.795000)
Panzano.M29.q50=c(3.800000, 4.200000, 9.751886, 8.748146, 3.468136, 2.800000, 3.719792, 3.100000, 3.599723, 2.993132, 3.863345, 2.721492, 3.140000, 3.005190, 3.405074, 3.901420, 3.468136, 4.653793, 3.004707, 3.140000, 3.070000, 2.907373, 2.975053, 3.047367, 2.900000, 3.204414, 3.201930, 3.131897, 2.987109, 3.100000, 3.600000)
Brijuni.M44.q50=c(4.538029, 3.970896, 3.552085, 4.269162, 5.200000, 3.222408, 3.100000, 3.280965, 3.528107, 2.950000, 2.972836, 3.080000, 3.348095, 3.450000, 2.866473, 3.158884, 2.902681, 3.082354, 2.985000, 3.561995, 3.550000, 3.336083, 4.221486, 3.675000, 3.731761, 3.645000, 3.705000, 4.450000, 3.468136, 4.099970, 5.750000, 3.167769)
Piran2.q50=c(3.120000, 3.312055, 3.209543, 3.360105, 3.300000, 3.326028, 3.282012, 3.408796, 3.135832, 3.570000, 3.070000, 3.220000, 3.197845, 2.950000, 2.930000, 2.878544,       NA, 3.004248, 2.940000, 2.827800, 3.360105, 2.775000, 3.460000, 3.251961, 3.970896, 3.324070, 3.700179, 3.704656, 3.560000, 3.802848, 3.190000, 3.504122)
all.Corbula.median.size=c(Po3.q50[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20],Po4.q50[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20],
Panzano.M28.q50[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>6],Panzano.M29.q50[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>6],
Brijuni.M44.q50[Brijuni.ID][pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20],Piran2.q50[Piran2.ID][pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>8])
labels.Corbula=c(rep("Po 3", length(Po3.q50[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20])),rep("Po 4", length(Po4.q50[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20])),
rep("Panzano M28", length(Panzano.M28.q50[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>6])),rep("Panzano M29", length(Panzano.M29.q50[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>6])),
rep("Brijuni", length(Brijuni.M44.q50[Brijuni.ID][pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20])),rep("Piran", length(Piran2.q50[Piran2.ID][pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>8])))

##########################
#MAXIMUM SIZE OF CORBULA
##########################
Po4.q95=c(13.718011, 11.462500,  9.765000,  9.982137,  8.370000, 12.452900, 11.556000, 12.225500, 12.026000, 11.905000, 12.420000,11.735000, 11.038750, 10.118500,  9.561850,  5.802950,  6.639800,  4.773250,  5.008200,  5.397279,  8.340100,  5.391100,5.815944,  6.887800,  8.036403,  5.998500,  4.152076,  4.613340,  3.912500,  3.589650,  2.322178,  8.005518)
Po3.q95=c(5.167042, 11.289509,  8.974500, 10.236650,  9.718366, 12.459000, 11.960000, 12.938000, 11.277600, 11.936628, 12.163550,12.570000, 11.230000, 11.380000, 11.710000, 10.310000,  8.743750, 11.070000, 11.148000,  4.661500,  8.244136,  7.373400,4.887450,  3.766373,  9.384400,  5.642000,  3.000000,  4.204500,  4.790876,  2.676400,  5.566950,  2.000000)
Panzano.M28.q95=c(8.904630, 10.920000, 11.797216, 12.604243,  5.883500, 12.522937, 11.243318, 10.021922,  6.549909,  4.720000,  5.617607,3.850000,  7.877623,  8.529382,  8.669644, 10.206697,  6.975000,  6.272000,  8.704713,  8.308285,  3.723170, 8.300000,5.815686,  6.431401,  4.614087,  5.930000,  6.805374,  8.580620,  9.552884,  6.236437,  4.661849)
Panzano.M29.q95=c(9.840000, 10.350000, 11.570745, 11.500000, 10.508984, 10.859889, 10.600000,  7.244815,  7.108000,  6.335000,  9.235326,7.021380,  9.085367,  7.442283,  6.009138,  6.302743,  5.200000,  6.045379,  8.499409,  9.170000,  4.700000,  8.984336,9.320577,  4.991748,  7.174000,  7.530900,  7.765840,  4.889348,  8.976331,  8.409628,  9.700000)
Piran2.q95=c(6.629500, 7.140000, 6.582500, 8.003152, 5.820000, 6.140000, 5.600000, 6.663041, 5.486500, 6.770000, 6.735000, 6.500000,7.291158, 5.491773, 5.988604, 5.800764,       NA, 5.733054, 5.600000, 5.039500, 5.525265, 8.632500, 9.559876, 6.862873,9.715727, 7.943093, 8.336393, 9.764370, 9.827000, 9.386000, 5.478274, 4.471595)
Brijuni.M44.q95=c(5.403803, 6.786496, 8.383945, 8.673009, 7.856491, 5.792215, 7.229052, 8.100000, 5.750000, 6.450000, 3.897500, 5.819765,7.150025, 6.160000, 4.895660, 7.600000, 6.364255, 5.439295, 5.574069, 8.162261, 5.814436, 5.445861, 5.685987, 5.671820,5.299989, 6.533727, 6.434494, 7.284608, 6.356045, 7.370000, 9.006709, 4.955000)
all.Corbula.maximum.size=c(Po3.q95[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20],Po4.q95[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20],Panzano.M28.q95[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>6],Panzano.M29.q95[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>6],Brijuni.M44.q95[Brijuni.ID][pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20],Piran2.q95[Piran2.ID][pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>8])

####################################################
#INCREMENTS - PYRITE PERCENTAGE IN CORBULA
####################################################
#per-increment percentages
pyrite.percent.Piran2.5cm=tapply(all.Corbula.Piran2.5cm.pyrite,all.Corbula.Piran2.5cm, mean, na.rm=T)
pyrite.percent.Brijuni.5cm=tapply(all.Corbula.Brijuni.5cm.pyrite,all.Corbula.Brijuni.5cm, mean, na.rm=T)
pyrite.percent.Po3.5cm=tapply(all.Corbula.Po3.5cm.pyrite,all.Corbula.Po3.5cm, mean, na.rm=T)
pyrite.percent.Po4.5cm=tapply(all.Corbula.Po4.5cm.pyrite,all.Corbula.Po4.5cm, mean, na.rm=T)
pyrite.percent.Panzano.M28.5cm=tapply(all.Corbula.Panzano.M28.5cm.pyrite,all.Corbula.Panzano.M28.5cm, mean, na.rm=T)
pyrite.percent.Panzano.M29.5cm=tapply(all.Corbula.Panzano.M29.5cm.pyrite,all.Corbula.Panzano.M29.5cm, mean, na.rm=T)
pyrite.percent.Panzano.5cm=tapply(all.Corbula.Panzano.5cm.pyrite,all.Corbula.Panzano.5cm, mean, na.rm=T)
all.pyrite.percent.for.Corbula=c(pyrite.percent.Po3.5cm[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20], pyrite.percent.Po4.5cm[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20], 
pyrite.percent.Panzano.M28.5cm[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>6],pyrite.percent.Panzano.M29.5cm[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>6],
pyrite.percent.Brijuni.5cm[pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20], pyrite.percent.Piran2.5cm[pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>8])


########################################################################################################################################
#ALL SITES - PYRITE FREQUENCY VERSUS CORBULA
########################################################################################################################################
options(scipen = 999)
par(mfrow=c(2,2))
par(mar=c(4,2,2,1))
plot(all.Corbula.percent, all.pyrite.percent.for.Corbula, pch=16, xlim=c(0,1),ylim=c(0,1),frame=F, xlab="Proportion of Corbula (relative to all mollusks)", ylab="Proportion of pyrite-lined shells")
points(all.Corbula.percent[labels.Corbula=="Po 3"], all.pyrite.percent.for.Corbula[labels.Corbula=="Po 3"], pch=21, bg="black")
points(all.Corbula.percent[labels.Corbula=="Po 4"], all.pyrite.percent.for.Corbula[labels.Corbula=="Po 4"], pch=22, bg="black")
points(all.Corbula.percent[labels.Corbula=="Panzano M28"], all.pyrite.percent.for.Corbula[labels.Corbula=="Panzano M28"], pch=21, bg="gray")
points(all.Corbula.percent[labels.Corbula=="Panzano M29"], all.pyrite.percent.for.Corbula[labels.Corbula=="Panzano M29"], pch=22, bg="gray")
points(all.Corbula.percent[labels.Corbula=="Piran"], all.pyrite.percent.for.Corbula[labels.Corbula=="Piran"], pch=22, bg="white")
points(all.Corbula.percent[labels.Corbula=="Brijuni"], all.pyrite.percent.for.Corbula[labels.Corbula=="Brijuni"], pch=21, bg="white")
legend(x="topleft", ncol=2, legend=c("Panzano 28","Panzano 29","Brijuni","Piran","Po 3","Po 4"), 
pt.bg=c("gray","gray","white","white","black","black"),pch=c(21,22,21,22,21,22), bty="n")

plot(all.Corbula.maximum.size, all.pyrite.percent.for.Corbula, pch=16, xlim=c(2,15), frame=F, ylim=c(0,1), xlab="Corbula maximum length (cm) (per ", ylab="Proportion of pyrite-lined shells")
points(all.Corbula.maximum.size[labels.Corbula=="Po 3"], all.pyrite.percent.for.Corbula[labels.Corbula=="Po 3"], pch=21, bg="black")
points(all.Corbula.maximum.size[labels.Corbula=="Po 4"], all.pyrite.percent.for.Corbula[labels.Corbula=="Po 4"], pch=22, bg="black")
points(all.Corbula.maximum.size[labels.Corbula=="Panzano M28"], all.pyrite.percent.for.Corbula[labels.Corbula=="Panzano M28"], pch=21, bg="gray")
points(all.Corbula.maximum.size[labels.Corbula=="Panzano M29"], all.pyrite.percent.for.Corbula[labels.Corbula=="Panzano M29"], pch=22, bg="gray")
points(all.Corbula.maximum.size[labels.Corbula=="Piran"], all.pyrite.percent.for.Corbula[labels.Corbula=="Piran"], pch=22, bg="white")
points(all.Corbula.maximum.size[labels.Corbula=="Brijuni"], all.pyrite.percent.for.Corbula[labels.Corbula=="Brijuni"], pch=21, bg="white")
########################################################################################################################################
#PER-CORE RELATONSHIP BETWEEN PERCENT CORBULA AND PYRITE FREQUENCY
########################################################################################################################################
par(mfrow=c(2,3))
plot(Corbula.percent.Po3[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20], pyrite.percent.Po3.5cm[pyrite.N.Po3.5cm>10 & cor.depths.Po3.5cm>20], pch=16, frame=F, xlim=c(0,1), ylim=c(0,1),xlab="Corbula proportion", ylab="Proportion of pyrite-lined shells")
plot(Corbula.percent.Po4[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20], pyrite.percent.Po4.5cm[pyrite.N.Po4.5cm>10 & cor.depths.Po4.5cm>20], pch=16, frame=F, xlim=c(0,1), ylim=c(0,1),xlab="Corbula proportion", ylab="Proportion of pyrite-lined shells")
plot(Corbula.percent.Panzano.M28[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>8], pyrite.percent.Panzano.M28.5cm[pyrite.N.Panzano.M28.5cm>10 & cor.depths.Panzano.M28.5cm>8], pch=16, frame=F, xlim=c(0,1), ylim=c(0,1),xlab="Corbula proportion", ylab="Proportion of pyrite-lined shells")
plot(Corbula.percent.Panzano.M29[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>8], pyrite.percent.Panzano.M29.5cm[pyrite.N.Panzano.M29.5cm>10 & cor.depths.Panzano.M29.5cm>8], pch=16, frame=F, xlim=c(0,1), ylim=c(0,1),xlab="Corbula proportion", ylab="Proportion of pyrite-lined shells")
plot(Corbula.percent.Piran2[Piran2.ID][pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>20], pyrite.percent.Piran2.5cm[pyrite.N.Piran2.5cm>10 & cor.depths.Piran2.5cm>20], pch=16, frame=F, xlim=c(0,1), ylim=c(0,1),xlab="Corbula proportion", ylab="Proportion of pyrite-lined shells")
plot(Corbula.percent.Brijuni[Brijuni.ID][pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20], pyrite.percent.Brijuni.5cm[pyrite.N.Brijuni.5cm>10 & cor.depths.Brijuni.5cm>20], pch=16, frame=F, xlim=c(0,1), ylim=c(0,1),xlab="Corbula proportion", ylab="Proportion of pyrite-lined shells")



##########################################################################################
#ADDITIONAL FIGURE - GEOGRAPHIC GRADIENT IN BIOEROSION
##########################################################################################
par(mfrow=c(1,5))
par(mar=c(20,1,1,0))
depths=as.numeric(names(tapply(Corbula.Po3.bioerosion,Corbula.Po3.finecm, mean)))
plot(tapply(Corbula.Po3.pyrite,Corbula.Po3.finecm, mean, na.rm=T),-depths, xlim=c(0,1), type="n", xlab="Bioerosion", main="Po 3",ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, 
ylim=c(-155,0))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po3.bioerosion,Corbula.Po3.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po3.bioerosion,Corbula.Po3.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
lines(tapply(Corbula.Po3.bioerosion,Corbula.Po3.finecm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Po4.bioerosion,Corbula.Po4.finecm, mean)))
plot(tapply(Corbula.Po4.bioerosion,Corbula.Po4.finecm, mean, na.rm=T),-depths, xlim=c(0,1), xlab="", main="Po 4", ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, 
ylim=c(-155,0))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po4.bioerosion,Corbula.Po4.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po4.bioerosion,Corbula.Po4.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
lines(tapply(Corbula.Po4.bioerosion,Corbula.Po4.finecm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.finecm, mean)))
plot(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.finecm, mean, na.rm=T),-depths, xlim=c(0,1), 
xlab="Bioerosion", main="Panzano", ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, ylim=c(-155,0))
abline(v=0.5, lty=2)
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray91", border=NA) 
out=tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
lines(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.finecm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Piran.bioerosion,Corbula.Piran.finecm, mean)))
plot(tapply(Corbula.Piran.bioerosion,Corbula.Piran.finecm, mean, na.rm=T),-depths, xlim=c(0,1), 
xlab="Bioerosion", main="Piran", ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, ylim=c(-155,0))
abline(v=0.5, lty=2)
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray91", border=NA) 
out=tapply(Corbula.Piran.bioerosion,Corbula.Piran.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Piran.bioerosion,Corbula.Piran.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
lines(tapply(Corbula.Piran.bioerosion,Corbula.Piran.finecm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.M44.bioerosion,Corbula.M44.finecm, mean)))
plot(tapply(Corbula.M44.bioerosion,Corbula.M44.finecm, mean, na.rm=T),-depths, xlim=c(0,1), 
xlab="Bioerosion", main="M44", ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, ylim=c(-155,0))
abline(v=0.5, lty=2)
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray91", border=NA) 
out=tapply(Corbula.M44.bioerosion,Corbula.M44.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.M44.bioerosion,Corbula.M44.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
lines(tapply(Corbula.M44.bioerosion,Corbula.M44.finecm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

##########################################################################################
#FIGURE 13 SUBSET - GEOGRAPHY OF PYRITE PRESERVATION TRENDS 
##########################################################################################
MINIMUM=10
##########################
#PYRITE AT PO3
##########################
par(mfrow=c(1,5))
par(mar=c(20,1,1,0))
samplesizes=tapply(Corbula.Po3.pyrite,Corbula.Po3.finecm, length)
depths=as.numeric(names(tapply(Corbula.Po3.pyrite,Corbula.Po3.finecm, mean)))
plot(tapply(Corbula.Po3.pyrite,Corbula.Po3.finecm, mean, na.rm=T),-depths, xlim=c(0,1), type="n",  xlab="Shell-lining pyrite", main="Po 3",ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, ylim=c(-155,0))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray61", border=NA) 
out=tapply(Corbula.Po3.pyrite,Corbula.Po3.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
#points(tapply(Corbula.Po3.pyrite,Corbula.Po3.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
#lines(tapply(Corbula.Po3.pyrite,Corbula.Po3.finecm, mean, na.rm=T),-depths)
#segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
out=tapply(Corbula.Po3.pyrite[Corbula.Po3.bioerosion<1],Corbula.Po3.finecm[Corbula.Po3.bioerosion<1], bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
means.not.bored=unlist(lapply(out,mean, na.rm=T))
points(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], pch=21, bg="black", cex=1.4)
lines(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], col="black")

samplesizes=tapply(all.Corbula.Po3.pyrite,all.Corbula.Po3.finecm, length)
all.depths=as.numeric(names(tapply(all.Corbula.Po3.pyrite,all.Corbula.Po3.finecm, mean)))
out=tapply(all.Corbula.Po3.pyrite,all.Corbula.Po3.finecm, 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[samplesizes>10],x1=UCI[samplesizes>MINIMUM],y0=-all.depths[samplesizes>MINIMUM], y1=-all.depths[samplesizes>MINIMUM], col="gray")
points(tapply(all.Corbula.Po3.pyrite,all.Corbula.Po3.finecm, mean, na.rm=T)[samplesizes>MINIMUM],-all.depths[samplesizes>MINIMUM], pch=21, bg="gray91", cex=1.4)
lines(tapply(all.Corbula.Po3.pyrite,all.Corbula.Po3.finecm, mean, na.rm=T)[samplesizes>MINIMUM],-all.depths[samplesizes>MINIMUM], col="gray91")
abline(v=0.5, lty=2)
text(0.8, -120, label="Early 20th c.", cex=1.2, srt=90)
text(0.8, -40, label="Late 20th c.", cex=1.2, srt=90)
text(0.8, -10, label="Mixed l.", cex=1.2, srt=90)

##########################
#PYRITE AT PO4
##########################
samplesizes=tapply(Corbula.Po4.pyrite,Corbula.Po4.finecm, length)
depths=as.numeric(names(tapply(Corbula.Po4.pyrite,Corbula.Po4.finecm, mean)))
plot(tapply(Corbula.Po4.pyrite,Corbula.Po4.finecm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", type="n",  xlab="Shell-lining pyrite", main="Po 4", ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, ylim=c(-155,0))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray61", border=NA) 
out=tapply(Corbula.Po4.pyrite,Corbula.Po4.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
#points(tapply(Corbula.Po4.pyrite,Corbula.Po4.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
#lines(tapply(Corbula.Po4.pyrite,Corbula.Po4.finecm, mean, na.rm=T),-depths)
#segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
out=tapply(Corbula.Po4.pyrite[Corbula.Po4.bioerosion<1],Corbula.Po4.finecm[Corbula.Po4.bioerosion<1], bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
means.not.bored=unlist(lapply(out,mean, na.rm=T))
points(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], pch=21, bg="black", cex=1.4)
lines(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], col="black")

samplesizes=tapply(all.Corbula.Po4.pyrite,all.Corbula.Po4.finecm, length)
all.depths=as.numeric(names(tapply(all.Corbula.Po4.pyrite,all.Corbula.Po4.finecm, mean)))
out=tapply(all.Corbula.Po4.pyrite,all.Corbula.Po4.finecm, 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[samplesizes>MINIMUM],x1=UCI[samplesizes>MINIMUM],y0=-all.depths[samplesizes>MINIMUM], y1=-all.depths[samplesizes>MINIMUM], col="gray")
points(tapply(all.Corbula.Po4.pyrite,all.Corbula.Po4.finecm, mean, na.rm=T)[samplesizes>MINIMUM],-all.depths[samplesizes>MINIMUM], pch=21, bg="gray91", cex=1.4)
lines(tapply(all.Corbula.Po4.pyrite,all.Corbula.Po4.finecm, mean, na.rm=T)[samplesizes>MINIMUM],-all.depths[samplesizes>MINIMUM], col="gray91")
abline(v=0.5, lty=2)
text(0.8, -120, label="Early 20th c.", cex=1.2, srt=90)
text(0.8, -40, label="Late 20th c.", cex=1.2, srt=90)
text(0.8, -10, label="Mixed l.", cex=1.2, srt=90)

##########################
#PYRITE AT PANZANO
##########################
samplesizes=tapply(Corbula.Panzano.M28.pyrite,Corbula.Panzano.M28.finecm, length)
depths=as.numeric(names(tapply(Corbula.Panzano.M28.pyrite,Corbula.Panzano.M28.finecm, mean)))
plot(tapply(Corbula.Panzano.M28.pyrite,Corbula.Panzano.M28.finecm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n",  xlab="Shell-lining pyrite", main="Panzano", ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, ylim=c(-155,0))
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray61", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Panzano.M28.pyrite,Corbula.Panzano.M28.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
#points(tapply(Corbula.Panzano.M28.pyrite,Corbula.Panzano.M28.finecm, mean, na.rm=T),-depths, pch=16, cex=1.4)
#lines(tapply(Corbula.Panzano.M28.pyrite,Corbula.Panzano.M28.finecm, mean, na.rm=T),-depths)
#segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
out=tapply(Corbula.Panzano.M28.pyrite[Corbula.Panzano.bioerosion<1],Corbula.Panzano.M28.finecm[Corbula.Panzano.bioerosion<1], bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
means.not.bored=unlist(lapply(out,mean, na.rm=T))
points(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], pch=21, bg="black", cex=1.4)
lines(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], col="black")

means.Panzano=tapply(c(all.Corbula.Panzano.M28.pyrite,all.Corbula.Panzano.M29.pyrite),c(all.Corbula.Panzano.M28.finecm,all.Corbula.Panzano.M29.finecm), mean)
samplesizes=tapply(c(all.Corbula.Panzano.M28.pyrite,all.Corbula.Panzano.M29.pyrite),c(all.Corbula.Panzano.M28.finecm,all.Corbula.Panzano.M29.finecm), length)
all.depths=as.numeric(names(tapply(c(all.Corbula.Panzano.M28.pyrite,all.Corbula.Panzano.M29.pyrite),c(all.Corbula.Panzano.M28.finecm,all.Corbula.Panzano.M29.finecm), mean)))
out=tapply(c(all.Corbula.Panzano.M28.pyrite,all.Corbula.Panzano.M29.pyrite),c(all.Corbula.Panzano.M28.finecm,all.Corbula.Panzano.M29.finecm), 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[samplesizes>MINIMUM],x1=UCI[samplesizes>MINIMUM],y0=-all.depths[samplesizes>MINIMUM], y1=-all.depths[samplesizes>MINIMUM], col="gray")
points(means.Panzano[samplesizes>MINIMUM],-all.depths[samplesizes>MINIMUM], pch=21, bg="gray91", cex=1.4)
lines(means.Panzano[samplesizes>MINIMUM],-all.depths[samplesizes>MINIMUM], col="gray91")
abline(v=0.5, lty=2)
text(0.8, -80, label="Highstand phase-Pre-20th c.", cex=1.2, srt=90)
text(0.8, -25, label="E. 20th c.", cex=1.2, srt=90)
text(0.8, -3, label="ML", cex=1.2, srt=90)

##########################
#PYRITE AT PIRAN
##########################
samplesizes=tapply(Corbula.Piran.pyrite,Corbula.Piran.finecm, length)
depths=as.numeric(names(tapply(Corbula.Piran.pyrite,Corbula.Piran.finecm, mean)))
plot(tapply(Corbula.Piran.pyrite,Corbula.Piran.finecm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n",  xlab="Shell-lining pyrite", main="Piran M53", ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, 
ylim=c(-155,0))
rect(xleft=0, xright=1, ytop=0, ybottom=-8, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-8, ybottom=-16, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-16, ybottom=-35, col="gray81", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-70, col="white", border=NA) 
out=tapply(Corbula.Piran.pyrite,Corbula.Piran.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Piran.pyrite,Corbula.Piran.finecm, mean, na.rm=T)[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM],  pch=21, bg="gray91", cex=1.4)
lines(tapply(Corbula.Piran.pyrite,Corbula.Piran.finecm, mean, na.rm=T)[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], col="gray91")
segments(x0=LCI[samplesizes>MINIMUM],x1=UCI[samplesizes>MINIMUM],y0=-depths[samplesizes>MINIMUM], y1=-depths[samplesizes>MINIMUM], col="gray91")

out=tapply(Corbula.Piran.pyrite[Corbula.Piran.bioerosion<1],Corbula.Piran.finecm[Corbula.Piran.bioerosion<1], bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
means.not.bored=unlist(lapply(out,mean, na.rm=T))
points(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], pch=21, bg="black", cex=1.4)
lines(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], col="black")
abline(v=0.5, lty=2)
legend(x="bottomright", pch=21, cex=1.2, pt.bg=c("gray","black"), 
legend=c("All","Non-bored"), bg="white") 
mean(tapply(Corbula.Piran.pyrite,Corbula.Piran.finecm, mean, na.rm=T))
mean(tapply(Corbula.Piran.pyrite[Corbula.Piran.bioerosion<1],Corbula.Piran.finecm[Corbula.Piran.bioerosion<1], mean, na.rm=T))
text(0.8, -80, label="Trangressive phase", cex=1.1, srt=90)
text(0.8, -28, label="MFZ", cex=1.1, srt=90)
text(0.8, -12, label="HST", cex=1.1, srt=0)
text(0.8, -5, label="ML", cex=1.1, srt=0)

##########################
#PYRITE AT BRIJUNI
##########################
samplesizes=tapply(Corbula.M44.pyrite,Corbula.M44.finecm, length)
depths=as.numeric(names(tapply(Corbula.M44.pyrite,Corbula.M44.finecm, mean)))
plot(tapply(Corbula.M44.pyrite,Corbula.M44.finecm, mean, na.rm=T),-depths, type="n", xlim=c(0,1),yaxt="n",  xlab="Shell-lining pyrite", main="Brijuni M44",ylab="Sediment depth (cm)", pch=16, cex=1.4, frame=F, ylim=c(-155,0))
rect(xleft=0, xright=1, ytop=0, ybottom=-20, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-20, ybottom=-90, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-90, ybottom=-120, col="gray81", border=NA) 
rect(xleft=0, xright=1, ytop=-120, ybottom=-155, col="white", border=NA) 
out=tapply(Corbula.M44.pyrite,Corbula.M44.finecm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
means=unlist(lapply(out,mean, na.rm=T))
points(means[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], pch=21, bg="gray91", cex=1.4)
lines(means[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], col="gray91")
segments(x0=LCI[samplesizes>MINIMUM],x1=UCI[samplesizes>MINIMUM],y0=-depths[samplesizes>MINIMUM], y1=-depths[samplesizes>MINIMUM], col="gray91")

out=tapply(Corbula.M44.pyrite[Corbula.M44.bioerosion<1],Corbula.M44.finecm[Corbula.M44.bioerosion<1], bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
means.not.bored=unlist(lapply(out,mean, na.rm=T))
points(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], pch=21, bg="black", cex=1.4)
lines(means.not.bored[samplesizes>MINIMUM],-depths[samplesizes>MINIMUM], col="black")
abline(v=0.5, lty=2)
text(0.8, -140, label="TST", cex=1.1, srt=90)
text(0.8, -105, label="MFZ", cex=1.1, srt=90)
text(0.8, -50, label="HST", cex=1.1, srt=90)
text(0.8, -10, label="ML", cex=1.1, srt=0)

mean(tapply(Corbula.M44.pyrite,Corbula.M44.finecm, mean, na.rm=T))
mean(tapply(Corbula.M44.pyrite[Corbula.M44.bioerosion<1],Corbula.M44.finecm[Corbula.M44.bioerosion<1], mean, na.rm=T))
##############################################################################
#END OF FIGURE 13 SUBSET
##############################################################################



##########################################################################################
#ARTICULATION
##########################################################################################
d.art=c(seq(4,20,by=4),seq(25,150,by=5))
art.Po3=c(0,2,0,6,3,9,6,8,4,8,4,6,7,3,5,2,2,5,2,0,0,0,0,2,0,1,0,0,1,0,0)
art.Po4=c(5,6,4,7,3,12,5,8,5,8,7,14,2,10,4,4,6,2,0,0,1,1,1,1,1,0,0,0,2,1,0)
art.Panzano=c(9,9,4,7,8,3,4,1,6,2,5,1,0,2,1,1,0,0,0,1,2,1,2,3,1,1,1,4,2,0,3)
art.Brijuni=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0)
art.Piran2=c(1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,1,0,0,2,1,0,2,0,0,1,0,1,1,0,0)

art.Po3.MNI=c(3,14,30,41,41,54,49,78,32,56,40,43,49,43,45,74,46,38,38,22,43,19,17,12,13,12,4,11,9,8,5)
art.Po4.MNI=c(14,32,9,19,20,67,68,55,60,41,73,88,20,82,50,57,84,44,12,31,15,23,18,15,16,11,11,8,5,7,2)
art.Panzano.MNI=c(83,114, 100, 127, 79, 59, 86, 74, 96, 139,71, 50, 34, 43, 40, 27, 32, 11, 21, 26, 39, 66, 48, 49, 47, 45, 48, 76, 50, 26, 37)
art.Brijuni.MNI=c(7,11,13,10,9,12,12,15,23,16,21,18,15,19,21,30,26,16,21,13,22,9,16,11,23,10,17,11,20,10,7)
art.Piran2.MNI=c(138,168,169,45,59,48,94,42,85,55,107,71,89,96,65,55,61,52,41,38,33,52,60,49,51,41,79,63,70,43,48)

disart.Po3=art.Po3.MNI-art.Po3
disart.Po4=art.Po4.MNI-art.Po4
disart.Panzano=art.Panzano.MNI-art.Panzano
disart.Brijuni=art.Brijuni.MNI-art.Brijuni
disart.Piran2=art.Piran2.MNI-art.Piran2

Po3.art.inc=art.Po3/art.Po3.MNI
Po4.art.inc=art.Po4/art.Po4.MNI
Panzano.art.inc=art.Panzano/art.Panzano.MNI
Brijuni.art.inc=art.Brijuni/art.Brijuni.MNI
Piran2.art.inc=art.Piran2/art.Piran2.MNI

temp.layer.0=numeric();temp.score.0=numeric()
for (i in 1:length(disart.Po3)) {
temp.layer.0=append(temp.layer.0,rep(i, times=disart.Po3[i]), length(temp.layer.0))
temp.score.0=append(temp.score.0,rep(0, times=disart.Po3[i]), length(temp.score.0))
}
temp.layer.1=numeric();temp.score.1=numeric(); 
for (i in 1:length(art.Po3)) {
temp.layer.1=append(temp.layer.1,rep(i, times=art.Po3[i]), length(temp.layer.1))
temp.score.1=append(temp.score.1,rep(1, times=art.Po3[i]), length(temp.score.1))
}
Po3.temp.layer=c(temp.layer.0,temp.layer.1); Po3.temp.score=c(temp.score.0,temp.score.1)

temp.layer.0=numeric();temp.score.0=numeric()
for (i in 1:length(disart.Po4)) {
temp.layer.0=append(temp.layer.0,rep(i, times=disart.Po4[i]), length(temp.layer.0))
temp.score.0=append(temp.score.0,rep(0, times=disart.Po4[i]), length(temp.score.0))
}
temp.layer.1=numeric();temp.score.1=numeric(); 
for (i in 1:length(art.Po4)) {
temp.layer.1=append(temp.layer.1,rep(i, times=art.Po4[i]), length(temp.layer.1))
temp.score.1=append(temp.score.1,rep(1, times=art.Po4[i]), length(temp.score.1))
}
Po4.temp.layer=c(temp.layer.0,temp.layer.1); Po4.temp.score=c(temp.score.0,temp.score.1)

temp.layer.0=numeric();temp.score.0=numeric()
for (i in 1:length(disart.Panzano)) {
temp.layer.0=append(temp.layer.0,rep(i, times=disart.Panzano[i]), length(temp.layer.0))
temp.score.0=append(temp.score.0,rep(0, times=disart.Panzano[i]), length(temp.score.0))
}
temp.layer.1=numeric();temp.score.1=numeric(); 
for (i in 1:length(art.Panzano)) {
temp.layer.1=append(temp.layer.1,rep(i, times=art.Panzano[i]), length(temp.layer.1))
temp.score.1=append(temp.score.1,rep(1, times=art.Panzano[i]), length(temp.score.1))
}
Panzano.temp.layer=c(temp.layer.0,temp.layer.1); Panzano.temp.score=c(temp.score.0,temp.score.1)

temp.layer.0=numeric();temp.score.0=numeric()
for (i in 1:length(disart.Brijuni)) {
temp.layer.0=append(temp.layer.0,rep(i, times=disart.Brijuni[i]), length(temp.layer.0))
temp.score.0=append(temp.score.0,rep(0, times=disart.Brijuni[i]), length(temp.score.0))
}
temp.layer.1=numeric();temp.score.1=numeric(); 
for (i in 1:length(art.Brijuni)) {
temp.layer.1=append(temp.layer.1,rep(i, times=art.Brijuni[i]), length(temp.layer.1))
temp.score.1=append(temp.score.1,rep(1, times=art.Brijuni[i]), length(temp.score.1))
}
Brijuni.temp.layer=c(temp.layer.0,temp.layer.1); Brijuni.temp.score=c(temp.score.0,temp.score.1)

temp.layer.0=numeric();temp.score.0=numeric()
for (i in 1:length(disart.Piran2)) {
temp.layer.0=append(temp.layer.0,rep(i, times=disart.Piran2[i]), length(temp.layer.0))
temp.score.0=append(temp.score.0,rep(0, times=disart.Piran2[i]), length(temp.score.0))
}
temp.layer.1=numeric();temp.score.1=numeric(); 
for (i in 1:length(art.Piran2)) {
temp.layer.1=append(temp.layer.1,rep(i, times=art.Piran2[i]), length(temp.layer.1))
temp.score.1=append(temp.score.1,rep(1, times=art.Piran2[i]), length(temp.score.1))
}
Piran2.temp.layer=c(temp.layer.0,temp.layer.1); Piran2.temp.score=c(temp.score.0,temp.score.1)
##########################################################################################
all.pyrite.percent.temp=c(pyrite.percent.Po3.5cm[6:31][art.Po3.MNI>10],
pyrite.percent.Po4.5cm[6:31][art.Po4.MNI>10],
pyrite.percent.Panzano.5cm[2:31][art.Panzano.MNI>10],
pyrite.percent.Piran2.5cm[2:31][art.Piran2.MNI>10],
pyrite.percent.Brijuni.5cm[6:31][art.Brijuni.MNI>10])

all.pyrite.percent.temp=c(pyrite.percent.Po3.5cm[6:31][art.Po3.MNI>10],pyrite.percent.Po4.5cm[6:31][art.Po4.MNI>10],pyrite.percent.Panzano.5cm[2:31][pyrite.N.Panzano.5cm>10],
pyrite.percent.Piran2.5cm[2:9],pyrite.percent.Brijuni.5cm[6:18])

all.art.percent.temp=c(Po3.art.inc[6:31][art.Po3.MNI>10],Po4.art.inc[6:31][art.Po4.MNI>10],Panzano.art.inc[2:31][art.Panzano.MNI>10],
Piran2.art.inc[c(3,7,11,15,19,22,27,31)],Brijuni.art.inc[c(7,9,11,13,15,17,19,21,23,25,27,29,31)])

core.labels.for.art=c(rep("Po 3",length(pyrite.percent.Po3.5cm[6:31][art.Po3.MNI>10])),rep("Po 4",length(pyrite.percent.Po4.5cm[6:31][art.Po4.MNI>10])),
rep("Panzano",length(pyrite.percent.Panzano.5cm[2:31][art.Panzano.MNI>10])),
rep("Piran",length(pyrite.percent.Piran2.5cm[2:9])),rep("Brijuni",length(pyrite.percent.Brijuni.5cm[6:18])))

par(mfrow=c(2,2))
plot(pyrite.percent.Po3.5cm[6:31][art.Po3.MNI>10], Po3.art.inc[6:31][art.Po3.MNI>10], pch=16, frame=F, ylim=c(0,0.5), xlab="Pyrite-lined shells", ylab="Articulated shells")
plot(pyrite.percent.Po4.5cm[6:31][art.Po4.MNI>10], Po4.art.inc[6:31][art.Po4.MNI>10], pch=16, frame=F, ylim=c(0,0.5), xlab="Pyrite-lined shells", ylab="Articulated shells")
plot(pyrite.percent.Panzano.5cm[2:31][pyrite.N.Panzano.5cm>10], Panzano.art.inc[2:31][pyrite.N.Panzano.5cm>10], pch=16, frame=F, ylim=c(0,0.5), xlab="Pyrite-lined shells", ylab="Articulated shells")
plot(all.pyrite.percent.temp,all.art.percent.temp, pch=16, frame=F, ylim=c(0,0.5), xlab="Pyrite-lined shells", ylab="Articulated shells")

cor.test(pyrite.percent.Po3.5cm[6:31][art.Po3.MNI>10], Po3.art.inc[6:31][art.Po3.MNI>10], method="s")
cor.test(pyrite.percent.Po4.5cm[6:31][art.Po4.MNI>10], Po4.art.inc[6:31][art.Po4.MNI>10], method="s")
cor.test(pyrite.percent.Panzano.5cm[2:31][pyrite.N.Panzano.5cm>10], Panzano.art.inc[2:31][pyrite.N.Panzano.5cm>10], method="s")
cor.test(all.pyrite.percent.temp,all.art.percent.temp, method="s")

##########################################################################################
#STRATIGRAPHIC TRENDS IN AGE-DATED SHELLS
##########################################################################################
#FIGURE 14 - TRENDS IN PERIOSTRACUM IN SEDIMENT CORES - PO AND PANZANO
##########################################################################################
par(mfrow=c(4,5))
par(mar=c(4,1,0,0))
depths=as.numeric(names(tapply(Corbula.Po3.conchiolin,Corbula.Po3.5cm, mean)))
plot(tapply(Corbula.Po3.conchiolin,Corbula.Po3.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))

rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po3.conchiolin,Corbula.Po3.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po3.conchiolin,Corbula.Po3.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Po3.conchiolin,Corbula.Po3.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
text(x=0, y=-30, labels="Late 20th c.", pos=4)
text(x=0, y=-145, labels="Early 20th c.", pos=4)

depths=as.numeric(names(tapply(Corbula.Po4.conchiolin,Corbula.Po4.5cm, mean)))
plot(tapply(Corbula.Po4.conchiolin,Corbula.Po4.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n",xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po4.conchiolin,Corbula.Po4.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po4.conchiolin,Corbula.Po4.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Po4.conchiolin,Corbula.Po4.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
text(x=0, y=-30, labels="Late 20th c.", pos=4)
text(x=0, y=-145, labels="Early 20th c.", pos=4)

depths=as.numeric(names(tapply(Corbula.Panzano.conchiolin,Corbula.Panzano.5cm, mean)))
plot(tapply(Corbula.Panzano.conchiolin,Corbula.Panzano.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray91", border=NA) 
out=tapply(Corbula.Panzano.conchiolin,Corbula.Panzano.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Panzano.conchiolin,Corbula.Panzano.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Panzano.conchiolin,Corbula.Panzano.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
text(x=0, y=-10, labels="Late 20th c.", pos=4)
text(x=0, y=-30, labels="Early 20th c.", pos=4)

depths=as.numeric(names(tapply(Corbula.Piran.conchiolin,Corbula.Piran.5cm, mean)))
plot(tapply(Corbula.Piran.conchiolin,Corbula.Piran.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-8, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-8, ybottom=-14, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray91", border=NA) 
out=tapply(Corbula.Piran.conchiolin,Corbula.Piran.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Piran.conchiolin,Corbula.Piran.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Piran.conchiolin,Corbula.Piran.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
text(x=0, y=-10, labels="HST", pos=4)
text(x=0, y=-25, labels="MFZ", pos=4)
text(x=0, y=-100, labels="TST", pos=4)

depths=as.numeric(names(tapply(Corbula.M44.conchiolin,Corbula.M44.5cm, mean)))
plot(tapply(Corbula.M44.conchiolin,Corbula.M44.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-20, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-20, ybottom=-90, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-90, ybottom=-120, col="gray91", border=NA) 
out=tapply(Corbula.M44.conchiolin,Corbula.M44.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.M44.conchiolin,Corbula.M44.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.M44.conchiolin,Corbula.M44.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
text(x=0, y=-55, labels="HST", pos=4)
text(x=0, y=-105, labels="MFZ", pos=4)
text(x=0, y=-140, labels="TST", pos=4)

######################################
#DISSOLUTION - PO AND PANZANO
######################################
depths=as.numeric(names(tapply(Corbula.Po3.dissolution,Corbula.Po3.5cm, mean)))
plot(tapply(Corbula.Po3.dissolution,Corbula.Po3.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po3.dissolution,Corbula.Po3.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po3.dissolution,Corbula.Po3.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Po3.dissolution,Corbula.Po3.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Po4.dissolution,Corbula.Po4.5cm, mean)))
plot(tapply(Corbula.Po4.dissolution,Corbula.Po4.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po4.dissolution,Corbula.Po4.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po4.dissolution,Corbula.Po4.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Po4.dissolution,Corbula.Po4.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Panzano.dissolution,Corbula.Panzano.5cm, mean)))
plot(tapply(Corbula.Panzano.dissolution,Corbula.Panzano.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray91", border=NA) 
out=tapply(Corbula.Panzano.dissolution,Corbula.Panzano.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Panzano.dissolution,Corbula.Panzano.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Panzano.dissolution,Corbula.Panzano.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Piran.dissolution,Corbula.Piran.5cm, mean)))
plot(tapply(Corbula.Piran.dissolution,Corbula.Piran.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-8, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-8, ybottom=-14, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray91", border=NA) 
out=tapply(Corbula.Piran.dissolution,Corbula.Piran.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Piran.dissolution,Corbula.Piran.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Piran.dissolution,Corbula.Piran.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.M44.dissolution,Corbula.M44.5cm, mean)))
plot(tapply(Corbula.M44.dissolution,Corbula.M44.5cm, mean, na.rm=T),-depths, xlim=c(0,1), 
yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-20, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-20, ybottom=-90, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-90, ybottom=-120, col="gray91", border=NA) 
out=tapply(Corbula.M44.dissolution,Corbula.M44.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.M44.dissolution,Corbula.M44.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.M44.dissolution,Corbula.M44.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

######################################
#ARTICULATION - PO AND PANZANO
######################################
#par(mfrow=c(3,4))
#par(mar=c(4,2,1,0))
depths=d.art
plot(art.Po3[art.Po3.MNI>10]/art.Po3.MNI[art.Po3.MNI>10],-depths[art.Po3.MNI>10], xlim=c(0,0.5),  yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Po3.temp.score,Po3.temp.layer, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(art.Po3[art.Po3.MNI>10]/art.Po3.MNI[art.Po3.MNI>10],-depths[art.Po3.MNI>10], pch=16)
lines(art.Po3[art.Po3.MNI>10]/art.Po3.MNI[art.Po3.MNI>10],-depths[art.Po3.MNI>10])
segments(x0=LCI[art.Po3.MNI>10],x1=UCI[art.Po3.MNI>10],y0=-depths[art.Po3.MNI>10], y1=-depths[art.Po3.MNI>10])

plot(art.Po4[art.Po4.MNI>10]/art.Po4.MNI[art.Po4.MNI>10],-depths[art.Po4.MNI>10], xlim=c(0,0.5), yaxt="n",  xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Po4.temp.score,Po4.temp.layer, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(art.Po4[art.Po4.MNI>10]/art.Po4.MNI[art.Po4.MNI>10],-depths[art.Po4.MNI>10], pch=16)
lines(art.Po4[art.Po4.MNI>10]/art.Po4.MNI[art.Po4.MNI>10],-depths[art.Po4.MNI>10])
segments(x0=LCI[art.Po4.MNI>10],x1=UCI[art.Po4.MNI>10],y0=-depths[art.Po4.MNI>10], y1=-depths[art.Po4.MNI>10])

plot(art.Panzano/art.Panzano.MNI,-depths, xlim=c(0,0.5), yaxt="n",  xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray91", border=NA) 
out=tapply(Panzano.temp.score,Panzano.temp.layer, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(art.Panzano/art.Panzano.MNI,-depths, pch=16)
lines(art.Panzano/art.Panzano.MNI,-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

plot(art.Piran2/art.Piran2.MNI,-depths, xlim=c(0,0.5), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
rect(xleft=0, xright=1, ytop=0, ybottom=-8, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-8, ybottom=-14, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray91", border=NA) 
out=tapply(Piran2.temp.score,Piran2.temp.layer, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(art.Piran2/art.Piran2.MNI,-depths, pch=16)
lines(art.Piran2/art.Piran2.MNI,-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
axis(2, at=seq(-150,0,by=50), labels=rep("",4))

plot(art.Brijuni/art.Brijuni.MNI,-depths, xlim=c(0,0.5), yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
rect(xleft=0, xright=1, ytop=0, ybottom=-20, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-20, ybottom=-90, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-90, ybottom=-120, col="gray91", border=NA) 
out=tapply(Brijuni.temp.score,Brijuni.temp.layer, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(art.Brijuni/art.Brijuni.MNI,-depths, pch=16)
lines(art.Brijuni/art.Brijuni.MNI,-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
axis(2, at=seq(-150,0,by=50), labels=rep("",4))

######################################
#BIOEROSION - PO AND PANZANO
######################################
depths=as.numeric(names(tapply(Corbula.Po3.bioerosion,Corbula.Po3.5cm, mean)))
plot(tapply(Corbula.Po3.bioerosion,Corbula.Po3.5cm, mean, na.rm=T),-depths, xlim=c(0,1), 
yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po3.bioerosion,Corbula.Po3.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po3.bioerosion,Corbula.Po3.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Po3.bioerosion,Corbula.Po3.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Po4.bioerosion,Corbula.Po4.5cm, mean)))
plot(tapply(Corbula.Po4.bioerosion,Corbula.Po4.5cm, mean, na.rm=T),-depths, xlim=c(0,1), 
yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-25, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-25, ybottom=-85, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-85, ybottom=-155, col="gray71", border=NA) 
out=tapply(Corbula.Po4.bioerosion,Corbula.Po4.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Po4.bioerosion,Corbula.Po4.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Po4.bioerosion,Corbula.Po4.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.5cm, mean)))
plot(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.5cm, mean, na.rm=T),-depths, xlim=c(0,1), yaxt="n", 
xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-6, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-6, ybottom=-14, col="gray41", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-35, ybottom=-155, col="gray91", border=NA) 
out=tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.Piran.bioerosion,Corbula.Piran.5cm, mean)))
plot(tapply(Corbula.Piran.bioerosion,Corbula.Piran.5cm, mean, na.rm=T),-depths, xlim=c(0,1), 
yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-8, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-8, ybottom=-14, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-14, ybottom=-35, col="gray91", border=NA) 
out=tapply(Corbula.Piran.bioerosion,Corbula.Piran.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.Piran.bioerosion,Corbula.Piran.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.Piran.bioerosion,Corbula.Piran.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)

depths=as.numeric(names(tapply(Corbula.M44.bioerosion,Corbula.M44.5cm, mean)))
plot(tapply(Corbula.M44.bioerosion,Corbula.M44.5cm, mean, na.rm=T),-depths, xlim=c(0,1), 
yaxt="n", xlab="", ylab="Sediment depth (cm)", pch=16, cex=1, frame=F, ylim=c(-155,0), main="")
axis(2, at=seq(-150,0,by=50), labels=rep("",4))
rect(xleft=0, xright=1, ytop=0, ybottom=-20, col="gray31", border=NA) 
rect(xleft=0, xright=1, ytop=-20, ybottom=-90, col="gray71", border=NA) 
rect(xleft=0, xright=1, ytop=-90, ybottom=-120, col="gray91", border=NA) 
out=tapply(Corbula.M44.bioerosion,Corbula.M44.5cm, bootci)
LCI=unlist(lapply(out,quantile,c(0.025), na.rm=T))
UCI=unlist(lapply(out,quantile,c(0.975), na.rm=T))
points(tapply(Corbula.M44.bioerosion,Corbula.M44.5cm, mean, na.rm=T),-depths, pch=16)
lines(tapply(Corbula.M44.bioerosion,Corbula.M44.5cm, mean, na.rm=T),-depths)
segments(x0=LCI,x1=UCI,y0=-depths, y1=-depths)
##############################################################
#END OF FIGURE 14
##############################################################

###########################################################################################################
#PER-INCREMENT AGE DISTRIBUTIONS
###########################################################################################################
Corbula.Po.ages.for.clock=c(Corbula.Po3.ages.for.clock, Corbula.Po4.ages.for.clock)
Corbula.Po.pyrite=c(Corbula.Po3.pyrite, Corbula.Po4.pyrite)
Corbula.Po.dissolution=c(Corbula.Po3.dissolution, Corbula.Po4.dissolution)
Corbula.Po.periostracum=c(Corbula.Po3.periostracum, Corbula.Po4.periostracum)
Corbula.Po.5cm=c(Corbula.Po3.5cm, Corbula.Po4.5cm)
Corbula.Po.10cm=c(Corbula.Po3.10cm, Corbula.Po4.10cm)

CUTOFF=15
CUTOFF=5

###########################################################################################################
#RATES OF BURIAL
###########################################################################################################
source("D:/Data/Codes/function_random-time_2-phase_exponential.txt")
par(mfrow=c(3,2))
Po3.Po4.ages.without.pyrite=c(Corbula.Po3.ages.for.clock[Corbula.Po3.5cm<21 & Corbula.Po3.pyrite==0],Corbula.Po4.ages.for.clock[Corbula.Po4.5cm<21 & Corbula.Po4.pyrite==0])
Po3.Po4.ages.with.pyrite=c(Corbula.Po3.ages.for.clock[Corbula.Po3.5cm<21 & Corbula.Po3.pyrite==1],Corbula.Po4.ages.for.clock[Corbula.Po4.5cm<21 & Corbula.Po4.pyrite==1])

1/mean(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm < 21)], na.rm=T)
1/mean(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm < 21)], na.rm=T)
1/mean(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm < 8)], na.rm=T)
1/mean(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm < 8)], na.rm=T)

age=Po3.Po4.ages.without.pyrite
age=age[!is.na(age)]

onephaserate<-1/mean(age)
lkexp=function(rho) {-length(age)*log(rho) + rho*sum(age)}
like<--(sum(log(onephaserate)-onephaserate*age))	#negative loglikelihood for exponential
AIC=2*1-(-2*like)
AICc=AIC+((2*1*(1+1))/(length(age)-1-1))
onephaseAIC<-AICc
#RANDOM 2-PHASE EXPONENTIAL
X=age
out<-FitExpMix(X)
randomtwophaserate1<-out$r1
randomtwophaserate2<-out$r2
randomtwophasealpha<-out$alpha
randomtwophasebeta<-out$beta
lik<--(out$lik)
likrandom=lik
rAIC=2*3-(-2*lik)
randomAIC=rAIC+((2*3*(3+1))/(length(X)-3-1))
tau<-randomtwophasealpha*(randomtwophaserate2-randomtwophaserate1)
lambda1<-randomtwophaserate1*randomtwophasealpha+randomtwophaserate2*(1-randomtwophasealpha)
lambda2<-randomtwophaserate1

onephaseAIC
randomAIC
onephaserate
lambda1


###########################################################################################################
#FIGURE 10 - UPPER PART - MODEL-FITTED AGE DISTRIBUTIONS
###########################################################################################################
expAIC=numeric()
onephaserate=numeric()
onephaserate.lci=numeric()
onephaserate.uci=numeric()
randomAIC=numeric()
tau=numeric()
lambda1=numeric()
lambda2=numeric()
like=numeric()
likrandom=numeric()
size=numeric()
par(mfcol=c(4,3))
par(mar=c(2,3,2,1))


for (i in 1:7) {
if (i==1) {age=Corbula.Po3.ages.for.clock[Corbula.Po3.5cm<21 & Corbula.Po3.pyrite==0]; MAX=20; COL="gray81"; MAIN="Po 3-without pyrite"; thickness=20}
if (i==2) {age=Corbula.Po3.ages.for.clock[Corbula.Po3.5cm<21 & Corbula.Po3.pyrite==1]; MAX=20; COL="gray51"; MAIN="Po 3-with pyrite"; thickness=20}
if (i==3) {age=Corbula.Po4.ages.for.clock[Corbula.Po4.5cm<21 & Corbula.Po4.pyrite==0]; MAX=40; COL="gray81"; MAIN="Po 4-without pyrite"; thickness=20}
if (i==4) {age=Corbula.Po4.ages.for.clock[Corbula.Po4.5cm<21 & Corbula.Po4.pyrite==1]; MAX=20; COL="gray51"; MAIN="Po 4-with pyrite"; thickness=20}
if (i==5) {age=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm < 8)]; MAX=20; COL="gray81"; MAIN="Panzano-without pyrite"; thickness=6}
if (i==6) {age=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm < 8)]; MAX=20; COL="gray51"; MAIN="Panzano-with pyrite"; thickness=6}
if (i==7) {age=Corbula.Piran.ages.for.clock[Corbula.Piran.5cm<13 & Corbula.Piran.pyrite==0]; MAX=20; COL="gray81"; MAIN="Piran-without pyrite"; thickness=8}

age=age[!is.na(age)]
size[i]=length(age)
onephaserate[i]<-1/mean(age)
like[i]<--(sum(log(onephaserate[i])-onephaserate[i]*age))
sd=fitdistr(age, densfun="exponential")$sd
onephaserate.lci[i]<-onephaserate[i]-2*sd
onephaserate.uci[i]<-onephaserate[i]+2*sd
	
AIC=2*1-(-2*like[i])
AICc=AIC+((2*1*(1+1))/(length(age)-1-1))
expAIC[i]<-AICc
onephaseAIC<-AICc

#RANDOM 2-PHASE EXPONENTIAL
X=age
source("D:/Data/Codes/function_random-time_2-phase_exponential.txt")
out<-FitExpMix(X)
randomtwophaserate1<-out$r1
randomtwophaserate2<-out$r2
randomtwophasealpha<-out$alpha
randomtwophasebeta<-out$beta
lik<--(out$lik)
likrandom[i]=lik
rAIC=2*3-(-2*lik)
randomAIC[i]=rAIC+((2*3*(3+1))/(length(age)-3-1))
tau[i]<-randomtwophasealpha*(randomtwophaserate2-randomtwophaserate1)
lambda1[i]<-randomtwophaserate1*randomtwophasealpha+randomtwophaserate2*(1-randomtwophasealpha)
lambda2[i]<-randomtwophaserate1
if (i < 7) {
BIN=5
hist(age, breaks=seq(0,210,by=BIN), xlab="Postmortem age (years)", ylab="Frequency",ylim=c(0,MAX), xlim=c(0,100), col=COL, main=MAIN, cex.main=0.9, bty="u")
times=seq(0,2000, by=BIN)+(BIN/2)
lambda=onephaserate[i]
predictSIMPLE<-dexp(times,rate=lambda)
lines(times, (predictSIMPLE/sum(predictSIMPLE))*length(age), lty=1, col="black", lwd=3)
if (i==3) {XLAB1=35;XLAB2=30;XLAB3=25}
if (i==1 | i==2 | i==4 | i==5 | i==6) {XLAB1=15;XLAB2=12.5;XLAB3=10;XLAB4=5}
text(50,XLAB1,labels=paste("Loss rate from ML =","",round(onephaserate[i], digits=3)))
#text(50,XLAB2,labels=paste("Time to loss =","",round(1/onephaserate[i], digits=1),"","years"))
text(50,XLAB2,labels=paste("Burial =","",round(thickness/(1/onephaserate[i]), digits=1),"","cm/y"))
text(50,XLAB3,labels=paste("(",round(thickness/(1/onephaserate.lci[i]), digits=2),"-",round(thickness/(1/onephaserate.uci[i]), digits=2),""," cm/y)",sep="" ))
if (i==2 | i==4) {plot.new(); plot.new()}
}
}


labels=c("Po 3-without pyrite","Po 3-with pyrite","Po 4-without pyrite", "Po 4-with pyrite","Panzano-without pyrite", "Panzano-with pyrite","Piran-without pyrite")
out=cbind(size,onephaserate, lambda1, lambda2, tau, like, likrandom, expAIC,  randomAIC)
rownames(out)=labels
colnames(out)=c("n","lambda", "lambda1", "lambda2", "tau", "Dis. model-neg. log-likelihood","Seq. model-neg. log-likelihood", "Disintegration model-expAICc",   "Sequestratiom model-AICc")
out
write.table(out, file="Fitting AFD of valves with and without pyrite.txt")

par(mfrow=c(4,3))
par(mar=c(2,3,2,1))

temp=Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==0]
hist(temp[temp<200], xlab="Postmortem age (y)", main="Po3-without pyrite", col="gray71", breaks=seq(0,200,5), ylim=c(0,60), cex.main=0.9)
abline(v=median(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==0], na.rm=T), lwd=2)
median(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==0], na.rm=T)

temp=Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==0]
hist(temp[temp<200], xlab="Postmortem age (y)", main="Po4-without pyrite", col="gray71", breaks=seq(0,200,5), ylim=c(0,60), cex.main=0.9)
abline(v=median(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==0], na.rm=T), lwd=2)
median(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==0], na.rm=T)

temp=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0]
hist(temp[temp<500], xlab="Postmortem age (y)", main="Panzano-without pyrite", col="gray71",border="gray71", breaks=seq(0,500,5), ylim=c(0,20), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0], na.rm=T), lwd=2)
median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0], na.rm=T)

hist(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==1], xlab="Postmortem age (y)", main="Po3-with pyrite",col="gray31", breaks=seq(0,200,5), ylim=c(0,60), cex.main=0.9)
abline(v=median(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==1], na.rm=T), lwd=2)
median(Corbula.Po3.ages.for.clock[Corbula.Po3.pyrite==1], na.rm=T)

hist(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==1], xlab="Postmortem age (y)", main="Po4-with pyrite",col="gray31", breaks=seq(0,200,5), ylim=c(0,60), cex.main=0.9)
abline(v=median(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==1], na.rm=T), lwd=2)
median(Corbula.Po4.ages.for.clock[Corbula.Po4.pyrite==1], na.rm=T)

temp=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1]
hist(temp[temp<500], xlab="Postmortem age (y)", main="Panzano-with pyrite",col="gray31", border="gray31", breaks=seq(0,500,5), ylim=c(0,20), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1], na.rm=T), lwd=2)
median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1], na.rm=T)
###########################################################################################################
#PIRAN AGE DISTRIBUTIONS
###########################################################################################################
temp=Corbula.Piran.ages.for.clock[Corbula.Piran.pyrite==0 & (Corbula.Piran.5cm < 13)]
hist(temp[temp<6000], xlab="Postmortem age (y)", main="Piran-without pyrite", col="gray71", breaks=seq(0,6000,250), ylim=c(0,60), cex.main=0.9)
abline(v=median(Corbula.Piran.ages.for.clock[Corbula.Piran.pyrite==0 & (Corbula.Piran.5cm < 13)], na.rm=T), lwd=2)
median(Corbula.Piran.ages.for.clock[Corbula.Piran.pyrite==0], na.rm=T)

temp=Corbula.Piran.ages.for.clock[Corbula.Piran.pyrite==1 & (Corbula.Piran.5cm < 13)]
hist(temp[temp<6000], xlab="Postmortem age (y)", main="Piran-with pyrite", col="gray71", breaks=seq(0,6000,250), ylim=c(0,60), cex.main=0.9)
abline(v=median(Corbula.Piran.ages.for.clock[Corbula.Piran.pyrite==1 & (Corbula.Piran.5cm < 13)], na.rm=T), lwd=2)
median(Corbula.Piran.ages.for.clock[Corbula.Piran.pyrite==1], na.rm=T)

###########################################################################################################
#DOWNCORE CHANGES IN AGE DISTRIBUTIONS WITH AND WITHOUT PYRITE AT PO AND PANZANO
###########################################################################################################
par(mfrow=c(5,4))
par(mar=c(4,2,2,1))
hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm < 21)], xlab="", main="Po-top 20 cm-without pyrite", col="gray71", breaks=seq(0,150,5), ylim=c(0,50), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm < 21)], na.rm=T), lwd=2)
hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm < 21)], xlab="", main="Po-top 20 cm-with pyrite",col="gray31", breaks=seq(0,150,5), ylim=c(0,50), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm < 21)], na.rm=T), lwd=2)
#plot.new()
hist(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm < 15)], xlab="", main="Panzano-top 14 cm-without pyrite", col="gray71", breaks=seq(0,500,20), ylim=c(0,50), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm < 15)], na.rm=T), lwd=2)
hist(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm < 15)], xlab="", main="Panzano-top 14 cm-without pyrite", col="gray31", breaks=seq(0,500,20), ylim=c(0,50), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm < 15)], na.rm=T), lwd=2)

hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 20 & Corbula.Po.5cm < 51)], xlab="", main="Po-20-50 cm", col="gray71", breaks=seq(0,150,5), ylim=c(0,20), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 20 & Corbula.Po.5cm < 51)], na.rm=T), lwd=2)
hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 20 & Corbula.Po.5cm < 51)], xlab="", main="Po-20-50 cm",col="gray31", breaks=seq(0,150,5), ylim=c(0,20), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 20 & Corbula.Po.5cm < 51)], na.rm=T), lwd=2)
#plot.new()
hist(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 14 & Corbula.Panzano.5cm < 37)], xlab="", main="Panzano-14-35 cm", col="gray71", breaks=seq(0,500,20), ylim=c(0,20), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 14 & Corbula.Panzano.5cm < 37)], na.rm=T), lwd=2)
hist(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 14 & Corbula.Panzano.5cm < 37)], xlab="", main="Panzano-14-35 cm", col="gray31", breaks=seq(0,500,20), ylim=c(0,20), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 14 & Corbula.Panzano.5cm < 37)], na.rm=T), lwd=2)
#
hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 50 & Corbula.Po.5cm < 90)], xlab="", main="Po-50-90 cm", col="gray71", breaks=seq(0,150,5), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 50 & Corbula.Po.5cm < 90)], na.rm=T), lwd=2)
hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 50 & Corbula.Po.5cm < 90)], xlab="", main="Po-50-90 cm",col="gray31", breaks=seq(0,150,5), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 50 & Corbula.Po.5cm < 90)], na.rm=T), lwd=2)
#plot.new()
hist(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 35 & Corbula.Panzano.5cm < 86)], xlab="", main="Panzano-35-85 cm", col="gray71", breaks=seq(0,500,20), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 35 & Corbula.Panzano.5cm < 86)], na.rm=T), lwd=2)
hist(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 35 & Corbula.Panzano.5cm < 86)], xlab="", main="Panzano-35-85 cm", col="gray31", breaks=seq(0,500,20), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 35 & Corbula.Panzano.5cm < 86)], na.rm=T), lwd=2)

temp=Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 90 & Corbula.Po.5cm < 121)]
hist(temp[temp<150], xlab="Postmortem age (y)", main="Po-90-120 cm", col="gray71", breaks=seq(0,150,5), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 90 & Corbula.Po.5cm < 121)], na.rm=T), lwd=2)
hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 90 & Corbula.Po.5cm < 121)], xlab="", main="Po-90-120 cm",col="gray31", breaks=seq(0,150,5), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 90 & Corbula.Po.5cm < 121)], na.rm=T), lwd=2)
#plot.new()
temp=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 85 & Corbula.Panzano.5cm < 131)]
hist(temp[temp<500], xlab="", main="Panzano-85-130 cm", col="gray71", breaks=seq(0,500,20), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 85 & Corbula.Panzano.5cm < 131)], na.rm=T), lwd=2)
hist(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 85 & Corbula.Panzano.5cm < 131)], xlab="", main="Panzano-85-130 cm", col="gray31", breaks=seq(0,500,20), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 85 & Corbula.Panzano.5cm < 131)], na.rm=T), lwd=2)

temp=Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 120)]
hist(temp[temp<150], xlab="Postmortem age (y)", main="Po-120-150 cm", col="gray71", breaks=seq(0,150,5), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==0 & (Corbula.Po.5cm > 120)], na.rm=T), lwd=2)
hist(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 120)], xlab="Postmortem age (y)", main="Po-120-150 cm",col="gray31", breaks=seq(0,150,5), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Po.ages.for.clock[Corbula.Po.pyrite==1 & (Corbula.Po.5cm > 120)], na.rm=T), lwd=2)
#plot.new()
temp=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 130)]
hist(temp[temp<500], xlab="Postmortem age (y)", main="Panzano-130-150 cm", col="gray71", breaks=seq(0,500,20), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==0 & (Corbula.Panzano.5cm > 130)], na.rm=T), lwd=2)
temp=Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 130)]
hist(temp[temp<500], xlab="Postmortem age (y)", main="Panzano-130-150 cm", col="gray31", breaks=seq(0,500,20), ylim=c(0,15), cex.main=0.9)
abline(v=median(Corbula.Panzano.ages.for.clock[Corbula.Panzano.pyrite==1 & (Corbula.Panzano.5cm > 130)], na.rm=T), lwd=2)


#########################################
#DEPTH BOXPLOTS WITHIN AGE COHORTS
#########################################
Corbula.Po3.10y.cohorts=cut(Corbula.Po3.ages.for.clock, breaks=seq(0,10000,by=10), labels=F)*10
Corbula.Po4.10y.cohorts=cut(Corbula.Po4.ages.for.clock, breaks=seq(0,10000,by=10), labels=F)*10
Corbula.Po.10y.cohorts=cut(Corbula.Po.ages.for.clock, breaks=seq(0,10000,by=10), labels=F)*10
Corbula.Panzano.25y.cohorts=cut(Corbula.Panzano.ages.for.clock, breaks=seq(0,10000,by=25), labels=F)*25
Corbula.Panzano.50y.cohorts=cut(Corbula.Panzano.ages.for.clock, breaks=seq(0,10000,by=50), labels=F)*50

#########################################
#FIGURE 10 - LOWER PART - SEDIMENT DEPTHS OF AGE COHORTS - PO3 and PO4 SEPARATED
#########################################
out1=lapply(split(-Corbula.Po3.5cm[Corbula.Po3.pyrite==0], Corbula.Po3.10y.cohorts[Corbula.Po3.pyrite==0])[1:5], function (x) quantile(bootci.ext(x)$resample, c(0.025,0.5,0.975)))
out2=lapply(split(-Corbula.Po3.5cm[Corbula.Po3.pyrite==1], Corbula.Po3.10y.cohorts[Corbula.Po3.pyrite==1])[1:5], function (x) quantile(bootci.ext(x)$resample, c(0.025,0.5,0.975)))
out3=lapply(split(-Corbula.Po4.5cm[Corbula.Po4.pyrite==0], Corbula.Po4.10y.cohorts[Corbula.Po4.pyrite==0])[1:5], function (x) quantile(bootci.ext(x)$resample, c(0.025,0.5,0.975)))
out4=lapply(split(-Corbula.Po4.5cm[Corbula.Po4.pyrite==1], Corbula.Po4.10y.cohorts[Corbula.Po4.pyrite==1])[1:5], function (x) quantile(bootci.ext(x)$resample, c(0.025,0.5,0.975)))
out5=lapply(split(-Corbula.Panzano.5cm[Corbula.Panzano.pyrite==0], Corbula.Panzano.50y.cohorts[Corbula.Panzano.pyrite==0])[1:5], function (x) quantile(bootci.ext(x)$resample, c(0.025,0.5,0.975)))
out6=lapply(split(-Corbula.Panzano.5cm[Corbula.Panzano.pyrite==1], Corbula.Panzano.50y.cohorts[Corbula.Panzano.pyrite==1])[1:5], function (x) quantile(bootci.ext(x)$resample, c(0.025,0.5,0.975)))

par(mfrow=c(3,3))
par(mar=c(4,4,2,1))
plot(1:10, c(out1$"10"[2],out2$"10"[2],out1$"20"[2],out2$"20"[2],out1$"30"[2],out2$"30"[2],out1$"40"[2],out2$"40"[2],out1$"50"[2],out2$"50"[2]), 
xlab="", ylim=c(-100,0), pch=21, bg=rep(c("black","gray"), 5),
ylab="Depth (cm)", xlim=c(0,11), xaxt="n", frame=F, main="Po prodelta (Po 3)", cex.main=0.9)
axis(1, at=c(1.5,3.5,5.5,7.5,9.5), labels=c("0-10 y","10-20 y","20-30 y","30-40 y","40-50 y"))
abline(v=c(2.5,4.5,6.5,8.5), lty=2)
segments(x0=1:10, x1=1:10,y0=c(out1$"10"[1],out2$"10"[1],out1$"20"[1],out2$"20"[1],out1$"30"[1],out2$"30"[1],out1$"40"[1],out2$"40"[1],out1$"50"[1],out2$"50"[1]),
y1=c(out1$"10"[3],out2$"10"[3],out1$"20"[3],out2$"20"[3],out1$"30"[3],out2$"30"[3],out1$"40"[3],out2$"40"[3],out1$"50"[3],out2$"50"[3]))
points(1:10, c(out1$"10"[2],out2$"10"[2],out1$"20"[2],out2$"20"[2],out1$"30"[2],out2$"30"[2],out1$"40"[2],out2$"40"[2],out1$"50"[2],out2$"50"[2]), pch=21, bg=rep(c("black","gray"), 5))

plot(1:10, c(out3$"10"[2],out4$"10"[2],out3$"20"[2],out4$"20"[2],out3$"30"[2],out4$"30"[2],out3$"40"[2],out4$"40"[2],out3$"50"[2],out4$"50"[2]), 
xlab="", ylim=c(-100,0), pch=21, bg=rep(c("black","gray"), 5),
ylab="Depth (cm)", xlim=c(0,11), xaxt="n", frame=F, main="Po prodelta (Po 4)", cex.main=0.9)
axis(1, at=c(1.5,3.5,5.5,7.5,9.5), labels=c("0-10 y","10-20 y","20-30 y","30-40 y","40-50 y"))
abline(v=c(2.5,4.5,6.5,8.5), lty=2)
segments(x0=1:10, x1=1:10,y0=c(out3$"10"[1],out4$"10"[1],out3$"20"[1],out4$"20"[1],out3$"30"[1],out4$"30"[1],out3$"40"[1],out4$"40"[1],out3$"50"[1],out4$"50"[1]),
y1=c(out3$"10"[3],out4$"10"[3],out3$"20"[3],out4$"20"[3],out3$"30"[3],out4$"30"[3],out3$"40"[3],out4$"40"[3],out3$"50"[3],out4$"50"[3]))
points(1:10, c(out3$"10"[2],out4$"10"[2],out3$"20"[2],out4$"20"[2],out3$"30"[2],out4$"30"[2],out3$"40"[2],out4$"40"[2],out3$"50"[2],out4$"50"[2]), pch=21, bg=rep(c("black","gray"), 5))

plot(1:10, c(out5$"50"[2],out6$"50"[2],out5$"100"[2],out6$"100"[2],out5$"150"[2],out6$"150"[2],out5$"200"[2],out6$"200"[2],out5$"250"[2],out6$"250"[2]), 
 xlab="", ylim=c(-100,0), pch=21, bg=rep(c("black","gray"), 5),
ylab="Depth (cm)", xlim=c(0.5,10.5), xaxt="n", frame=F, main="Panzano core", cex.main=0.9)
axis(1, at=c(1.5,3.5,5.5,7.5,9.5), labels=c("0-50 y","50-100 y","100-150 y","150-200 y","200-250 y"))
abline(v=c(2.5,4.5,6.5,8.5), lty=2)
segments(x0=1:10, x1=1:10,y0=c(out5$"50"[1],out6$"50"[1],out5$"100"[1],out6$"100"[1],out5$"150"[1],out6$"150"[1],out5$"200"[1],out6$"200"[1],out5$"250"[1],out6$"250"[1]),
y1=c(out5$"50"[3],out6$"50"[3],out5$"100"[3],out6$"100"[3],out5$"150"[3],out6$"150"[3],out5$"200"[3],out6$"200"[3],out5$"250"[3],out6$"250"[3]))
legend(x="topright", legend=c("without pyrite","with pyrite"), pch=21, pt.bg=c("black","gray"), bg="white")
points(1:10, c(out5$"50"[2],out6$"50"[2],out5$"100"[2],out6$"100"[2],out5$"150"[2],out6$"150"[2],out5$"200"[2],out6$"200"[2],out5$"250"[2],out6$"250"[2]), 
pch=21, bg=rep(c("black","gray"), 5))

###########################################################################################################
#PER-INCREMENT FREQUENCIES OF ALTERATION
###########################################################################################################
#INCLUDING SURFACE MIXED LAYER
#########################################################################################################################

bioerosion.Panzano.inc=tapply(Corbula.Panzano.bioerosion,Corbula.Panzano.5cm, mean, na.rm=T)
bioerosion.Po4.inc=tapply(Corbula.Po4.bioerosion,Corbula.Po4.5cm, mean, na.rm=T)
bioerosion.Po3.inc=tapply(Corbula.Po3.bioerosion,Corbula.Po3.5cm, mean, na.rm=T)
bioerosion.Piran.inc=tapply(Corbula.Piran.bioerosion,Corbula.Piran.5cm, mean, na.rm=T)
bioerosion.Brijuni.inc=tapply(Corbula.M44.bioerosion,Corbula.M44.5cm, mean, na.rm=T)

ext.encrustation.Panzano.inc=tapply(Corbula.Panzano.ext.encrustation,Corbula.Panzano.5cm, mean, na.rm=T)
ext.encrustation.Po4.inc=tapply(Corbula.Po4.ext.encrustation,Corbula.Po4.5cm, mean, na.rm=T)
ext.encrustation.Po3.inc=tapply(Corbula.Po3.ext.encrustation,Corbula.Po3.5cm, mean, na.rm=T)
ext.encrustation.Piran.inc=tapply(Corbula.Piran.ext.encrustation,Corbula.Piran.5cm, mean, na.rm=T)
ext.encrustation.Brijuni.inc=tapply(Corbula.M44.ext.encrustation,Corbula.M44.5cm, mean, na.rm=T)

dissolution.Panzano.inc=tapply(Corbula.Panzano.dissolution,Corbula.Panzano.5cm, mean, na.rm=T)
dissolution.Po4.inc=tapply(Corbula.Po4.dissolution,Corbula.Po4.5cm, mean, na.rm=T)
dissolution.Po3.inc=tapply(Corbula.Po3.dissolution,Corbula.Po3.5cm, mean, na.rm=T)
dissolution.Piran.inc=tapply(Corbula.Piran.dissolution,Corbula.Piran.5cm, mean, na.rm=T)
dissolution.Brijuni.inc=tapply(Corbula.M44.dissolution,Corbula.M44.5cm, mean, na.rm=T)

pyrite.Panzano.inc=tapply(Corbula.Panzano.pyrite,Corbula.Panzano.5cm, mean, na.rm=T)
pyrite.Po4.inc=tapply(Corbula.Po4.pyrite,Corbula.Po4.5cm, mean, na.rm=T)
pyrite.Po3.inc=tapply(Corbula.Po3.pyrite,Corbula.Po3.5cm, mean, na.rm=T)
pyrite.Piran.inc=tapply(Corbula.Piran.pyrite,Corbula.Piran.5cm, mean, na.rm=T)
pyrite.Brijuni.inc=tapply(Corbula.M44.pyrite,Corbula.M44.5cm, mean, na.rm=T)

periostracum.Panzano.inc=tapply(Corbula.Panzano.periostracum,Corbula.Panzano.5cm, mean, na.rm=T)
periostracum.Po4.inc=tapply(Corbula.Po4.periostracum,Corbula.Po4.5cm, mean, na.rm=T)
periostracum.Po3.inc=tapply(Corbula.Po3.periostracum,Corbula.Po3.5cm, mean, na.rm=T)
periostracum.Piran.inc=tapply(Corbula.Piran.periostracum,Corbula.Piran.5cm, mean, na.rm=T)
periostracum.Brijuni.inc=tapply(Corbula.M44.periostracum,Corbula.M44.5cm, mean, na.rm=T)

ornamentation.Panzano.inc=tapply(Corbula.Panzano.ornamentation,Corbula.Panzano.5cm, mean, na.rm=T)
ornamentation.Po4.inc=tapply(Corbula.Po4.ornamentation,Corbula.Po4.5cm, mean, na.rm=T)
ornamentation.Po3.inc=tapply(Corbula.Po3.ornamentation,Corbula.Po3.5cm, mean, na.rm=T)
ornamentation.Piran.inc=tapply(Corbula.Piran.ornamentation,Corbula.Piran.5cm, mean, na.rm=T)
ornamentation.Brijuni.inc=tapply(Corbula.M44.ornamentation,Corbula.M44.5cm, mean, na.rm=T)

ornamentation.inc=c(ornamentation.Po3.inc, ornamentation.Po4.inc,ornamentation.Panzano.inc, ornamentation.Piran.inc, ornamentation.Brijuni.inc)
bioerosion.inc=c(bioerosion.Po3.inc, bioerosion.Po4.inc,bioerosion.Panzano.inc, bioerosion.Piran.inc, bioerosion.Brijuni.inc)
ext.encrustation.inc=c(ext.encrustation.Po3.inc, ext.encrustation.Po4.inc,ext.encrustation.Panzano.inc, ext.encrustation.Piran.inc, ext.encrustation.Brijuni.inc)
dissolution.inc=c(dissolution.Po3.inc, dissolution.Po4.inc,dissolution.Panzano.inc, dissolution.Piran.inc, dissolution.Brijuni.inc)
pyrite.inc=c(pyrite.Po3.inc,pyrite.Po4.inc,pyrite.Panzano.inc,pyrite.Piran.inc,pyrite.Brijuni.inc)
periostracum.inc=c(periostracum.Po3.inc,periostracum.Po4.inc,periostracum.Panzano.inc,periostracum.Piran.inc,periostracum.Brijuni.inc)

#########################################################################################################################
#EXCLUDING SURFACE MIXED LAYER
#########################################################################################################################
bioerosion.Panzano.inc=tapply(Corbula.Panzano.bioerosion[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
bioerosion.Po4.inc=tapply(Corbula.Po4.bioerosion[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
bioerosion.Po3.inc=tapply(Corbula.Po3.bioerosion[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
bioerosion.Piran.inc=tapply(Corbula.Piran.bioerosion[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
bioerosion.Brijuni.inc=tapply(Corbula.M44.bioerosion[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

ext.encrustation.Panzano.inc=tapply(Corbula.Panzano.ext.encrustation[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
ext.encrustation.Po4.inc=tapply(Corbula.Po4.ext.encrustation[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
ext.encrustation.Po3.inc=tapply(Corbula.Po3.ext.encrustation[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
ext.encrustation.Piran.inc=tapply(Corbula.Piran.ext.encrustation[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
ext.encrustation.Brijuni.inc=tapply(Corbula.M44.ext.encrustation[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

dissolution.Panzano.inc=tapply(Corbula.Panzano.dissolution[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
dissolution.Po4.inc=tapply(Corbula.Po4.dissolution[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
dissolution.Po3.inc=tapply(Corbula.Po3.dissolution[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
dissolution.Piran.inc=tapply(Corbula.Piran.dissolution[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8],mean, na.rm=T)
dissolution.Brijuni.inc=tapply(Corbula.M44.dissolution[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20],mean, na.rm=T)

pyrite.Panzano.inc=tapply(Corbula.Panzano.pyrite[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
pyrite.Po4.inc=tapply(Corbula.Po4.pyrite[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
pyrite.Po3.inc=tapply(Corbula.Po3.pyrite[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
pyrite.Piran.inc=tapply(Corbula.Piran.pyrite[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
pyrite.Brijuni.inc=tapply(Corbula.M44.pyrite[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

periostracum.Panzano.inc=tapply(Corbula.Panzano.periostracum[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
periostracum.Po4.inc=tapply(Corbula.Po4.periostracum[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
periostracum.Po3.inc=tapply(Corbula.Po3.periostracum[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
periostracum.Piran.inc=tapply(Corbula.Piran.periostracum[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
periostracum.Brijuni.inc=tapply(Corbula.M44.periostracum[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

int.conch.Panzano.inc=tapply(Corbula.Panzano.int.conch[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
int.conch.Po4.inc=tapply(Corbula.Po4.int.conch[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
int.conch.Po3.inc=tapply(Corbula.Po3.int.conch[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
int.conch.Piran.inc=tapply(Corbula.Piran.int.conch[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
int.conch.Brijuni.inc=tapply(Corbula.M44.int.conch[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

ornamentation.Panzano.inc=tapply(Corbula.Panzano.ornamentation[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
ornamentation.Po4.inc=tapply(Corbula.Po4.ornamentation[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
ornamentation.Po3.inc=tapply(Corbula.Po3.ornamentation[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
ornamentation.Piran.inc=tapply(Corbula.Piran.ornamentation[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
ornamentation.Brijuni.inc=tapply(Corbula.M44.ornamentation[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

staining.Panzano.inc=tapply(Corbula.Panzano.stained[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
staining.Po4.inc=tapply(Corbula.Po4.stained[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
staining.Po3.inc=tapply(Corbula.Po3.stained[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
staining.Piran.inc=tapply(Corbula.Piran.stained[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
staining.Brijuni.inc=tapply(Corbula.M44.stained[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

conchiolin.Panzano.inc=tapply(Corbula.Panzano.conchiolin[Corbula.Panzano.5cm>6],Corbula.Panzano.5cm[Corbula.Panzano.5cm>6], mean, na.rm=T)
conchiolin.Po4.inc=tapply(Corbula.Po4.conchiolin[Corbula.Po4.5cm>20],Corbula.Po4.5cm[Corbula.Po4.5cm>20], mean, na.rm=T)
conchiolin.Po3.inc=tapply(Corbula.Po3.conchiolin[Corbula.Po3.5cm>20],Corbula.Po3.5cm[Corbula.Po3.5cm>20], mean, na.rm=T)
conchiolin.Piran.inc=tapply(Corbula.Piran.conchiolin[Corbula.Piran.5cm>8],Corbula.Piran.5cm[Corbula.Piran.5cm>8], mean, na.rm=T)
conchiolin.Brijuni.inc=tapply(Corbula.M44.conchiolin[Corbula.M44.5cm>20],Corbula.M44.5cm[Corbula.M44.5cm>20], mean, na.rm=T)

#internal conchiolin and bioerosion is just scored for dated shells, so the number of levels is not ~31
int.conch.percent.Piran2.5cm=tapply(Corbula.Piran.int.conch,Corbula.Piran.5cm, mean, na.rm=T)
int.conch.percent.Po3.5cm=tapply(Corbula.Po3.int.conch,Corbula.Po3.5cm, mean, na.rm=T)
int.conch.percent.Po4.5cm=tapply(Corbula.Po4.int.conch,Corbula.Po4.5cm, mean, na.rm=T)
int.conch.percent.Panzano.5cm=tapply(Corbula.Panzano.int.conch,Corbula.Panzano.5cm, mean, na.rm=T)

#match TOC levels with Corbula levels
#not all increments at Piran and Brijuni were scored for pyrite in Corbula
#for Brijuni
#temp1=c(4,8,12,16,20,seq(25,155,by=5))
#temp2=c(4,8,12,16,20,seq(30,150,by=10),160)
temp1=c(4,8,12,16,20,seq(25,155,by=5))
temp2=c(4,8,12,16,20,seq(30,150,by=10),160)
Brijuni.ID=match(temp2, temp1)
#for Piran
temp1=c(4,8,12,16,20,seq(25,150,by=5))
temp2=c(4,12,30,50,70,90,105,130,150)
Piran2.ID=match(temp2, temp1)
#for Po3
temp1=c(4,8,12,16,20,seq(25,150,by=5))
temp2=c(4,8,12,16,20,25,45,70,85,100,110,115,135,140)
Po3.ID=match(temp2, temp1)
#for Po4
temp1=c(4,8,12,16,20,seq(25,150,by=5))
temp2=c(4,8,12,16,20,35,50,65,80,95,110,130,135,150)
Po4.ID=match(temp2, temp1)
#for Panzano
temp1=c(4,8,12,16,20,seq(25,150,by=5))
temp2=c(4,8,12,20,35,50,65,70,90,95,115,130,150)
Panzano.ID=match(temp2, temp1)

art.Po3.MNI
art.Po4.MNI
pyrite.N.Panzano.5cm
Brijuni.art.inc.shortseries=Brijuni.art.inc[Brijuni.ID[6:length(Brijuni.ID)]]
Piran2.art.inc.shortseries=Piran2.art.inc[Piran2.ID[2:length(Piran2.ID)]]
Po3.art.inc.shortseries=Po3.art.inc[Po3.ID[6:length(Po3.ID)]]
Po4.art.inc.shortseries=Po4.art.inc[Po4.ID[6:length(Po4.ID)]]
Panzano.art.inc.shortseries=Panzano.art.inc[Panzano.ID[2:length(Panzano.ID)]]

Po3.depths.no.ML=c(25,45,70,85,100,110,115,135,140)
Po4.depths.no.ML=c(35,50,65,80,95,110,130,135,150)
Panzano.depths.no.ML=c(8,12,20,35,50,65,70,90,95,115,130,150)
Piran.depths.no.ML=c(12,30,50,70,90,105,130,150)
Brijuni.depths.no.ML=c(seq(30,150,by=10),160)

staining.inc=c(staining.Po3.inc, staining.Po4.inc,staining.Panzano.inc,staining.Piran.inc,staining.Brijuni.inc)
ornamentation.inc=c(ornamentation.Po3.inc, ornamentation.Po4.inc,ornamentation.Panzano.inc, ornamentation.Piran.inc, ornamentation.Brijuni.inc)
bioerosion.inc=c(bioerosion.Po3.inc, bioerosion.Po4.inc,bioerosion.Panzano.inc, bioerosion.Piran.inc, bioerosion.Brijuni.inc)
ext.encrustation.inc=c(ext.encrustation.Po3.inc, ext.encrustation.Po4.inc,ext.encrustation.Panzano.inc, ext.encrustation.Piran.inc, ext.encrustation.Brijuni.inc)
dissolution.inc=c(dissolution.Po3.inc, dissolution.Po4.inc,dissolution.Panzano.inc, dissolution.Piran.inc, dissolution.Brijuni.inc)
pyrite.inc=c(pyrite.Po3.inc,pyrite.Po4.inc,pyrite.Panzano.inc,pyrite.Piran.inc, pyrite.Brijuni.inc)
periostracum.inc=c(periostracum.Po3.inc,periostracum.Po4.inc,periostracum.Panzano.inc,periostracum.Piran.inc,periostracum.Brijuni.inc)
int.conch.inc=c(int.conch.Po3.inc,int.conch.Po4.inc,int.conch.Panzano.inc,int.conch.Piran.inc,rep(1, length(periostracum.Brijuni.inc)))
conchiolin.inc=c(conchiolin.Po3.inc,conchiolin.Po4.inc,conchiolin.Panzano.inc,conchiolin.Piran.inc,rep(1, length(periostracum.Brijuni.inc)))
core.labels=c(rep("Po", length(pyrite.Po3.inc)), rep("Po", length(pyrite.Po4.inc)), rep("Panzano", length(pyrite.Panzano.inc)),rep("Piran", length(pyrite.Piran.inc)),rep("Brijuni", length(pyrite.Brijuni.inc)))
all.art.for.shortseries=c(Po3.art.inc.shortseries,Po4.art.inc.shortseries,Panzano.art.inc.shortseries,Piran2.art.inc.shortseries,Brijuni.art.inc.shortseries)
core.depths.noML=c(Po3.depths.no.ML,Po4.depths.no.ML,Panzano.depths.no.ML,Piran.depths.no.ML,Brijuni.depths.no.ML)

##################################################################################################################
#FINAL FIGURE 12
##################################################################################################################
par(mfrow=c(4,3))
par(mar=c(4,4,1,0))

plot(all.pyrite.percent.temp,1-all.art.percent.temp, pch=16, frame=F, xlim=c(0,1), ylim=c(0.5,1), 
xlab="Pyrite-lined valves", ylab="Disarticulation")
points(all.pyrite.percent.temp[core.labels.for.art=="Po 3"], 1-all.art.percent.temp[core.labels.for.art=="Po 3"], pch=16)
points(all.pyrite.percent.temp[core.labels.for.art=="Po 4"], 1-all.art.percent.temp[core.labels.for.art=="Po 4"], pch=16)
points(all.pyrite.percent.temp[core.labels.for.art=="Panzano"], 1-all.art.percent.temp[core.labels.for.art=="Panzano"], pch=21, bg="gray")
points(all.pyrite.percent.temp[core.labels.for.art=="Piran"], 1-all.art.percent.temp[core.labels.for.art=="Piran"], pch=21, bg="white")
points(all.pyrite.percent.temp[core.labels.for.art=="Brijuni"], 1-all.art.percent.temp[core.labels.for.art=="Brijuni"], pch=22, bg="white")

plot(conchiolin.inc, pyrite.inc, pch=16, bg="black", frame=F, xlab="Conchiolin loss", ylab="Pyrite-lined valves", ylim=c(0,1), xlim=c(0,1))
points(conchiolin.inc[core.labels=="Po"], pyrite.inc[core.labels=="Po"], pch=16)
points(conchiolin.inc[core.labels=="Panzano"], pyrite.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(conchiolin.inc[core.labels=="Piran"], pyrite.inc[core.labels=="Piran"], pch=21, bg="white")
points(conchiolin.inc[core.labels=="Brijuni"], pyrite.inc[core.labels=="Brijuni"], pch=22, bg="white")

plot(bioerosion.inc, pyrite.inc, pch=16, bg="black", frame=F, xlab="Bioerosion", ylab="Pyrite-lined valves", ylim=c(0,1), xlim=c(0,1))
points(bioerosion.inc[core.labels=="Po"], pyrite.inc[core.labels=="Po"], pch=16)
points(bioerosion.inc[core.labels=="Panzano"], pyrite.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(bioerosion.inc[core.labels=="Piran"], pyrite.inc[core.labels=="Piran"], pch=21, bg="white")
points(bioerosion.inc[core.labels=="Brijuni"], pyrite.inc[core.labels=="Brijuni"], pch=22, bg="white")

plot(staining.inc, pyrite.inc, pch=16, bg="black", frame=F, xlab="Staining", ylab="Pyrite-lined valves", ylim=c(0,1), xlim=c(0,1))
points(staining.inc[core.labels=="Po"], pyrite.inc[core.labels=="Po"], pch=16)
points(staining.inc[core.labels=="Panzano"], pyrite.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(staining.inc[core.labels=="Piran"], pyrite.inc[core.labels=="Piran"], pch=21, bg="white")
points(staining.inc[core.labels=="Brijuni"], pyrite.inc[core.labels=="Brijuni"], pch=22, bg="white")
#legend(x="topright", legend=c("Po","Panzano","Piran","Brijuni"), pt.bg=c("black","gray","white","white"), pch=c(21,21,21,22))


plot(ornamentation.inc, pyrite.inc, pch=16, bg="black", frame=F, xlab="Ornamentation loss", ylab="Pyrite-lined valves", ylim=c(0,1), xlim=c(0,1))
points(ornamentation.inc[core.labels=="Po"], pyrite.inc[core.labels=="Po"], pch=16)
points(ornamentation.inc[core.labels=="Panzano"], pyrite.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(ornamentation.inc[core.labels=="Piran"], pyrite.inc[core.labels=="Piran"], pch=21, bg="white")
points(ornamentation.inc[core.labels=="Brijuni"], pyrite.inc[core.labels=="Brijuni"], pch=22, bg="white")

plot(ext.encrustation.inc, pyrite.inc, pch=16, bg="black", frame=F, xlab="Encrustation", ylab="Pyrite-lined valves", ylim=c(0,1), xlim=c(0,1))
points(ext.encrustation.inc[core.labels=="Po"], pyrite.inc[core.labels=="Po"], pch=16)
points(ext.encrustation.inc[core.labels=="Panzano"], pyrite.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(ext.encrustation.inc[core.labels=="Piran"], pyrite.inc[core.labels=="Piran"], pch=21, bg="white")
points(ext.encrustation.inc[core.labels=="Brijuni"], pyrite.inc[core.labels=="Brijuni"], pch=22, bg="white")
#legend(x="topright", legend=c("Po","Panzano","Piran","Brijuni"), pt.bg=c("black","gray","white","white"), pch=c(21,21,21,22))

plot(dissolution.inc, pyrite.inc, pch=16, bg="black", frame=F, xlab="Dissolution", ylab="Pyrite-lined valves", ylim=c(0,1), xlim=c(0,1))
points(dissolution.inc[core.labels=="Po"], pyrite.inc[core.labels=="Po"], pch=16)
points(dissolution.inc[core.labels=="Panzano"], pyrite.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(dissolution.inc[core.labels=="Piran"], pyrite.inc[core.labels=="Piran"], pch=21, bg="white")
points(dissolution.inc[core.labels=="Brijuni"], pyrite.inc[core.labels=="Brijuni"], pch=22, bg="white")
#legend(x="topleft", legend=c("Po","Panzano","Piran","Brijuni"), pt.bg=c("black","gray","white","white"), pch=c(21,21,21,22))


plot(periostracum.inc, 1-all.art.for.shortseries, pch=16, bg="black", frame=F, ylab="Disarticulation", xlab="Periostracum loss", ylim=c(0.5,1), xlim=c(0,1))
points(periostracum.inc[core.labels=="Po 3"], 1-all.art.for.shortseries[core.labels=="Po 3"], pch=16)
points(periostracum.inc[core.labels=="Po 4"], 1-all.art.for.shortseries[core.labels=="Po 4"], pch=16)
points(periostracum.inc[core.labels=="Panzano"], 1-all.art.for.shortseries[core.labels=="Panzano"], pch=21, bg="gray")
points(periostracum.inc[core.labels=="Piran"], 1-all.art.for.shortseries[core.labels=="Piran"], pch=21, bg="white")
points(periostracum.inc[core.labels=="Brijuni"], 1-all.art.for.shortseries[core.labels=="Brijuni"], pch=22, bg="white")

plot(bioerosion.inc, 1-all.art.for.shortseries, pch=16, bg="black", frame=F, ylab="Disarticulation", xlab="Bioerosion", ylim=c(0.5,1), xlim=c(0,1))
points(bioerosion.inc[core.labels=="Po 3"], 1-all.art.for.shortseries[core.labels=="Po 3"], pch=16)
points(bioerosion.inc[core.labels=="Po 4"], 1-all.art.for.shortseries[core.labels=="Po 4"], pch=16)
points(bioerosion.inc[core.labels=="Panzano"], 1-all.art.for.shortseries[core.labels=="Panzano"], pch=21, bg="gray")
points(bioerosion.inc[core.labels=="Piran"], 1-all.art.for.shortseries[core.labels=="Piran"], pch=21, bg="white")
points(bioerosion.inc[core.labels=="Brijuni"], 1-all.art.for.shortseries[core.labels=="Brijuni"], pch=22, bg="white")
legend(x="bottomright", legend=c("Po","Panzano","Piran","Brijuni"), pt.bg=c("black","gray","white","white"), pch=c(21,21,21,22))

plot(ext.encrustation.inc, bioerosion.inc, pch=16, bg="black", frame=F, xlab="Encrustation", ylab="Bioerosion", ylim=c(0,1), xlim=c(0,1))
points(ext.encrustation.inc[core.labels=="Po"], bioerosion.inc[core.labels=="Po"], pch=16)
points(ext.encrustation.inc[core.labels=="Panzano"], bioerosion.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(ext.encrustation.inc[core.labels=="Piran"], bioerosion.inc[core.labels=="Piran"], pch=21, bg="white")
points(ext.encrustation.inc[core.labels=="Brijuni"], bioerosion.inc[core.labels=="Brijuni"], pch=22, bg="white")

plot(staining.inc, bioerosion.inc, pch=16, bg="black", frame=F, xlab="Staining", ylab="Bioerosion", ylim=c(0,1), xlim=c(0,1))
points(staining.inc[core.labels=="Po"], bioerosion.inc[core.labels=="Po"], pch=16)
points(staining.inc[core.labels=="Panzano"], bioerosion.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(staining.inc[core.labels=="Piran"], bioerosion.inc[core.labels=="Piran"], pch=21, bg="white")
points(staining.inc[core.labels=="Brijuni"], bioerosion.inc[core.labels=="Brijuni"], pch=22, bg="white")

plot(ornamentation.inc, bioerosion.inc, pch=16, bg="black", frame=F, xlab="Loss of ornamentation", ylab="Bioerosion", ylim=c(0,1), xlim=c(0,1))
points(ornamentation.inc[core.labels=="Po"], bioerosion.inc[core.labels=="Po"], pch=16)
points(ornamentation.inc[core.labels=="Panzano"], bioerosion.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(ornamentation.inc[core.labels=="Piran"], bioerosion.inc[core.labels=="Piran"], pch=21, bg="white")
points(ornamentation.inc[core.labels=="Brijuni"], bioerosion.inc[core.labels=="Brijuni"], pch=22, bg="white")

###################################33
#ADDITIONAL PLOTS
###################################33
plot(staining.inc, ext.encrustation.inc, pch=16, bg="black", frame=F, xlab="Staining", ylab="Encrustation", ylim=c(0,1), xlim=c(0,1))
points(staining.inc[core.labels=="Po"], ext.encrustation.inc[core.labels=="Po"], pch=16)
points(staining.inc[core.labels=="Panzano"], ext.encrustation.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(staining.inc[core.labels=="Piran"], ext.encrustation.inc[core.labels=="Piran"], pch=21, bg="white")
points(staining.inc[core.labels=="Brijuni"], ext.encrustation.inc[core.labels=="Brijuni"], pch=22, bg="white")

plot(ornamentation.inc, ext.encrustation.inc, pch=16, bg="black", frame=F, xlab="Ornamentation loss", ylab="Encrustation", ylim=c(0,1), xlim=c(0,1))
points(ornamentation.inc[core.labels=="Po"], ext.encrustation.inc[core.labels=="Po"], pch=16)
points(ornamentation.inc[core.labels=="Panzano"], ext.encrustation.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(ornamentation.inc[core.labels=="Piran"], ext.encrustation.inc[core.labels=="Piran"], pch=21, bg="white")
points(ornamentation.inc[core.labels=="Brijuni"], ext.encrustation.inc[core.labels=="Brijuni"], pch=22, bg="white")

plot(ornamentation.inc, staining.inc, pch=16, bg="black", frame=F, xlab="Ornamentation loss", ylab="Staining", ylim=c(0,1), xlim=c(0,1))
points(ornamentation.inc[core.labels=="Po"], staining.inc[core.labels=="Po"], pch=16)
points(ornamentation.inc[core.labels=="Panzano"], staining.inc[core.labels=="Panzano"], pch=21, bg="gray")
points(ornamentation.inc[core.labels=="Piran"], staining.inc[core.labels=="Piran"], pch=21, bg="white")
points(ornamentation.inc[core.labels=="Brijuni"], staining.inc[core.labels=="Brijuni"], pch=22, bg="white")


############################################################################################################################################################################################################
#FIGURE 13 UPPER PART - MULTIVARIATE ANALYSES
############################################################################################################################################################################################################
boxplot(split(pyrite.inc, core.labels), col="gray", frame=F)


par(mfrow=c(2,3))
par(mar=c(4,4,2,0))
out=cbind(pyrite.inc, conchiolin.inc, ext.encrustation.inc, bioerosion.inc, ornamentation.inc, staining.inc, dissolution.inc)
pco=cmdscale(vegdist(out, method="manhattan"), eig=T)
plot(pco$points[,1], pco$points[,2], xlim=c(-2.5, 4.5), ylim=c(-3.5,3), pch=16, xlab="PCO 1 (81%)", ylab="PCO 2 (12%)", frame=F, cex=1.4)
points(pco$points[core.labels=="Piran",1], pco$points[core.labels=="Piran",2], pch=21, bg="white", cex=1.4)
points(pco$points[core.labels=="Brijuni",1], pco$points[core.labels=="Brijuni",2], pch=22, bg="white", cex=1.4)
points(pco$points[core.labels=="Po",1], pco$points[core.labels=="Po",2], pch=21, bg="black", cex=1.4)
points(pco$points[core.labels=="Panzano",1], pco$points[core.labels=="Panzano",2], pch=21, bg="gray", cex=1.4)
pco$eig/sum(pco$eig)
fit=envfit(pco, out, perm = 999)
envfit.cors=scores(fit, "vectors")
plot(fit, p.max = 0.05, col = "gray31", lwd=2, labels=c("pyrite", "conchiolin loss", "encrustation", "bioerosion", "ornamentation loss", "staining", "dissolution"))

barplot(envfit.cors[,1], ylim=c(-1, 1), cex=1, cex.names=0.9, las=2,  names.arg =c("pyrite",  "conchiolin loss", "encrustation", "bioerosion", "ornamentation loss", "staining", "dissolution"))

cbind(core.depths.noML,pco$points[,1])

out=cbind(pyrite.inc, conchiolin.inc, ext.encrustation.inc, bioerosion.inc, ornamentation.inc, staining.inc, dissolution.inc)
out=out[which(core.labels=="Po" | core.labels=="Panzano"),]
new.depths=core.depths.noML[which(core.labels=="Po" | core.labels=="Panzano")]
new.labels=core.labels[which(core.labels=="Po" | core.labels=="Panzano")]
new.time=c(rep("Po-AST",4),rep("Po-HST",5),rep("Po-AST",4),rep("Po-HST",5),rep("Panzano-AST",4),rep("Panzano-HST",8))
pco=cmdscale(vegdist(out, method="manhattan"), eig=T)

plot(pco$points[,1], pco$points[,2], xlim=c(-1.25, 1.25), ylim=c(-1.25,1.25), pch=16, xlab="PCO 1 (62%)", ylab="PCO 2 (18%)", frame=F, cex=1.4)
pco$eig/sum(pco$eig)
ordihull(pco, groups=new.time, display = "sites", draw = c("polygon"), col=c("gray","white","gray","white"))
ordisurf(x=pco, y=out[,"pyrite.inc"], family="quasibinomial", add=T, col="black")
points(pco$points[new.time=="Po-AST",1], pco$points[new.time=="Po-AST",2], pch=21, bg="black", cex=1.4)
points(pco$points[new.time=="Panzano-AST",1], pco$points[new.time=="Panzano-AST",2], pch=22, bg="black", cex=1.4)
points(pco$points[new.time=="Po-HST",1], pco$points[new.time=="Po-HST",2], pch=21, bg="white", cex=1.4)
points(pco$points[new.time=="Panzano-HST",1], pco$points[new.time=="Panzano-HST",2], pch=22, bg="white", cex=1.4)