#Tomašových, A., Kidwell, S.M., Alexander, C.R. and Kaufman, D.S., 2019. Millennial‐scale age offsets within fossil assemblages: Result of bioturbation below the taphonomic active zone and out‐of‐phase production. Paleoceanography and Paleoclimatology, 34(6), pp.954-977. #Age distributions are first scaled up to the total number of individuals in each increment (just extrapolating by resampling from the shape of distribution #from adjacent increments), and then those full distributions are divided by the survival (preservation) function to account for the disintegrated specimens. #The actual parameters that are used for lambda 1 (disintegration in the taphonomically active zone) and tau (rate at which lambda 1 declines to very low value or zero - #this can be actually burial rate) are in the script. The estimation of both lambda 1 and tau can be noisy especially when input of shells is not #constant in time - so shells from the whole southern California Bight were used in model fitting, not just from cores themselves (which would also work). #Postmortem age data from amino acid racemization calibrated with 14C can have large age errors - therefore, #when estimating uncertainties, there are also script parts where shell ages are drawn from lognormal distribution (for Nuculana) #and gamma distribution (Parvilucina), with or without bootstrap, and where rates of disintegration are also repeated with these uncertainties. #note that the model parameters in the survival function (r1, r2, alpha) are just tranformations of lambda1, lambda2, tau ####################################################################### install.packages("TruncatedDistributions", repos="http://R-Forge.R-project.org") require(TruncatedDistributions) require(truncdist) require(gridBase) require(grid) require(vegan) #MODEL WITH DISINTEGRATION DECLINING FROM LAMBDA 1 to LAMBDA 2 AT RATE TAU (SEE CONVERSION TO R1, R2 below) FitExpMix=function(X,niters=10000,a=2,b=0.001,rmean=mean(X)){ ## prior on beta=Beta(a,a); posterior = Beta(a+N1,a+N2) ## prior on r1 & r2 =Gamma(b,1/(rmean*b); posterior = Gamma(b+n,1/(rmean*b+mean(X1 or X2)*n)) N=length(X) d=function(x){ r1*exp(-r1*x)*beta+r2*exp(-r2*x)*(1-beta) } beta=0.5; r1=1; r2=1; best=c(beta,r1,r2); dbest=sum(log(d(X))) cluster=rbinom(N,1,beta); c=which(cluster==1) for(iter in 1:niters){ beta=rbeta(1,a+length(c),a+N-length(c)) r1=rgamma(1,b+length(c),rate=rmean*b+sum(X[c])); r2=rgamma(1,b+N-length(c),rate=rmean*b+sum(X[-c])) L1=beta*r1*exp(-r1*X) L2=(1-beta)*r2*exp(-r2*X) cluster=rbinom(N,1,L1/(L1+L2)); c=which(cluster==1) dnew=sum(log(d(X))); if(dnew>dbest){ dbest=dnew;best=c(beta,r1,r2) } } beta=best[1];r1=best[2];r2=best[3] if(r1>r2){beta=1-beta;r1temp=r1;r1=r2;r2=r1temp} out=list();out$beta=beta;out$r1=r1;out$r2=r2 ## convert parameters out$alpha=(beta*r1)/(beta*r1+(1-beta)*r2) out$tau=out$alpha*(r2-r1);out$lambda1=r2-out$tau;out$lambda2=r1 out$lik=dbest out } #MODEL WITH DISINTEGRATION DECLINING FROM LAMBDA 1 to LAMBDA 2 AT RATE TAU #TIME AND TIMEMAX REFER TO OFFSET AND ONSET IN PRODUCTION OneStateConstantTwiceTruncatedExpMix=function(theta){ lambda1=exp(theta[1]) lambda2=exp(theta[2]) tau=exp(theta[3]) N=length(ages) first01=(tau*(exp(-lambda2*TIME)-exp(-lambda2*TIMEMAX)))/lambda2 first02=((lambda1-lambda2)*(exp(-(tau+lambda1)*TIME)-exp(-(tau+lambda1)*TIMEMAX)))/(tau+lambda1) #first01=ifelse(first01>0,first01,10e-50) #first02=ifelse(first02>0,first02,10e-50) first=(1/(first01+first02)) second=(tau*exp(-lambda2*ages) + (lambda1-lambda2)*exp(-(tau+lambda1)*ages)) #second=ifelse(second>0,second,10e-50) prsecond=sum(log(second)) loglik=N*log(first)+prsecond if (lambda2>=(lambda1)) {loglik=-100000} -(loglik) } #MODEL WITH CONSTANT DISINTEGRATION #TIME AND TIMEMAX REFER TO OFFSET AND ONSET IN PRODUCTION OneStateConstantTwiceTruncatedExp=function(theta){ lambda1=exp(theta) N=length(ages) first01=0 first02=(lambda1*(exp(-lambda1*TIME)-exp(-lambda1*TIMEMAX)))/lambda1 first=(1/(first01+first02)) second=lambda1*exp(-lambda1*ages) prsecond=sum(log(second)) loglik=N*log(first)+prsecond -(loglik) } ###################################################################################################################### horiz.hist <- function(Data, breaks="Sturges", col="transparent", las=1, ylim=range(HBreaks), labelat=pretty(ylim), labels=labelat, border=par("fg"), ... ) {a <- hist(Data, plot=FALSE, breaks=breaks) HBreaks <- a$breaks HBreak1 <- a$breaks[1] hpos <<- function(Pos) (Pos-HBreak1)*(length(HBreaks)-1)/ diff(range(HBreaks)) barplot(a$counts, space=0, horiz=T, ylim=hpos(ylim), col=col, border=border,...) axis(2, at=hpos(labelat), labels=labels, las=las, ...) print("use hpos() to address y-coordinates") } ############################################################################################# setwd("D:/Data/Melville_2012/MS-Statistic_project") load("Parvilucina and Nuculana postmortem age data.Rdata") ############################################################################################# #ASSIGN SHELL AGES TO DECIMETER-SCALE UNITS (BEDS) - 3 UNITS IN EACH CORE, SEPARATELY FOR EACH SPECIES ############################################################################################# Nuculana.PVL10.50.bed1=c(Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths<23]) Nuculana.PVL10.50.bed2=c(Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths>23 & Nuculana.PVL1050.meanSdepths<35]) Nuculana.PVL10.50.bed3=c(Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths>35]) Nuculana.OC50.list1=split(Nuculana.OC50.ages,Nuculana.OC50.meanSdepths) Nuculana.OC.50.bed1=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths<23] Nuculana.OC.50.bed2=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths>23 & Nuculana.OC50.meanSdepths<38] Nuculana.OC.50.bed3=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths>38] Parvilucina.PVL10.50.bed1=c(Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths<23]) Parvilucina.PVL10.50.bed2=c(Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths>23 & Parvilucina.PVL1050.meanSdepths<35]) Parvilucina.PVL10.50.bed3=c(Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths>35]) Parvilucina.OC.50.bed1=Parvilucina.OC50.ages[Parvilucina.OC50.meanSdepths<23] Parvilucina.OC.50.bed2=Parvilucina.OC50.ages[Parvilucina.OC50.meanSdepths>23 & Parvilucina.OC50.meanSdepths<38] Parvilucina.OC.50.bed3=Parvilucina.OC50.ages[Parvilucina.OC50.meanSdepths>38] ############################################################################################################################################################ #SURFACE AGE-FREQUENCY DISTRIBUTIONS AT PVL10-50 (BED 1 - upper 24 cm), OC-50 (20 cm), #and all Van Veen grabs from the Southern California shelf shallower than 50 m (Tomasovych et al. 2016 Paleobiology age data) ############################################################################################################################################################ #cairo_pdf(file='plot.pdf', width=7, height=7) par(mfrow=c(2,3)) par(mar=c(4,3,1,0)) hist(Nuculana.PVL10.50.bed1[Nuculana.PVL10.50.bed1<200], breaks=seq(0,200,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="Palos Verdes Shelf-core top 24 cm", cex.main=0.9, xlab="") hist(Parvilucina.PVL10.50.bed1[Parvilucina.PVL10.50.bed1<200], add=T, breaks=seq(0,200,by=10), col=rgb(1,1,1,0.5)) legend(x=100, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1.3) hist(Nuculana.OC.50.bed1[Nuculana.OC.50.bed1<200], breaks=seq(0,200,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="San Pedro Shelf-core top 24 cm", cex.main=0.9, xlab="") hist(Parvilucina.OC.50.bed1[Parvilucina.OC.50.bed1<200], add=T, breaks=seq(0,200,by=10), col=rgb(1,1,1,0.5)) legend(x=100, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1.3) hist(Nuculana.shallower.50[Nuculana.shallower.50<200], breaks=seq(0,200,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="Southern California Bight-grabs", cex.main=0.9, xlab="") hist(Parvilucina.shallower.50[Parvilucina.shallower.50<200], add=T, breaks=seq(0,200,by=10), col=rgb(1,1,1,0.5)) legend(x=100, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1.3) #dev.off() ############################################################################################################################################################# #SURFACE AGE-FREQUENCY DISTRIBUTIONS AT PVL10-50 (BED 1 - upper 24 cm), OC-50 (20 cm), #and all Van Veen grabs from the Southern California shelf shallower than 50 m (Tomasovych et al. 2016 Paleobiology data) #FINER BINS and DIFFERENT DURATIONS ############################################################################################################################################################ par(mfrow=c(3,3)) par(mar=c(4,2,1,0)) hist(Nuculana.PVL10.50.bed1[Nuculana.PVL10.50.bed1<200], breaks=seq(0,200,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="Palos Verdes Shelf-core top 24 cm", cex.main=0.9, xlab="") hist(Parvilucina.PVL10.50.bed1[Parvilucina.PVL10.50.bed1<200], add=T, breaks=seq(0,200,by=10), col=rgb(1,1,1,0.5)) legend(x=100, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1.3) hist(Nuculana.OC.50.bed1[Nuculana.OC.50.bed1<200], breaks=seq(0,200,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="San Pedro Shelf-core top 24 cm", cex.main=0.9, xlab="") hist(Parvilucina.OC.50.bed1[Parvilucina.OC.50.bed1<200], add=T, breaks=seq(0,200,by=10), col=rgb(1,1,1,0.5)) legend(x=100, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1.3) hist(Nuculana.shallower.50[Nuculana.shallower.50<200], breaks=seq(0,200,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="Southern California Bight-grabs", cex.main=0.9, xlab="") hist(Parvilucina.shallower.50[Parvilucina.shallower.50<200], add=T, breaks=seq(0,200,by=10), col=rgb(1,1,1,0.5)) legend(x=100, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1.3) hist(Nuculana.PVL10.50.bed1[Nuculana.PVL10.50.bed1<1000], breaks=seq(0,1000,by=25), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,60), main="Palos Verdes Shelf-core top 24 cm", cex.main=0.9, xlab="") hist(Parvilucina.PVL10.50.bed1[Parvilucina.PVL10.50.bed1<1000], add=T, breaks=seq(0,1000,by=25), col=rgb(1,1,1,0.5)) hist(Nuculana.OC.50.bed1[Nuculana.OC.50.bed1<1000], breaks=seq(0,1000,by=25), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,60), main="San Pedro Shelf-core top 24 cm", cex.main=0.9, xlab="") hist(Parvilucina.OC.50.bed1[Parvilucina.OC.50.bed1<1000], add=T, breaks=seq(0,1000,by=25), col=rgb(1,1,1,0.5)) hist(Nuculana.shallower.50[Nuculana.shallower.50<1000], breaks=seq(0,1000,by=25), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,60), main="Southern California Bight-grabs", cex.main=0.9, xlab="") hist(Parvilucina.shallower.50[Parvilucina.shallower.50<1000], add=T, breaks=seq(0,1000,by=25), col=rgb(1,1,1,0.5)) hist(Nuculana.PVL10.50.bed1[Nuculana.PVL10.50.bed1<5000], breaks=seq(0,5000,by=200), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,80), main="Palos Verdes Shelf-core top 24 cm", cex.main=0.9, xlab="Postmortem age (y)") hist(Parvilucina.PVL10.50.bed1[Parvilucina.PVL10.50.bed1<5000], add=T, breaks=seq(0,5000,by=200), col=rgb(1,1,1,0.5)) hist(Nuculana.OC.50.bed1[Nuculana.OC.50.bed1<5000], breaks=seq(0,5000,by=200), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,80), main="San Pedro Shelf-core top 24 cm", cex.main=0.9, xlab="Postmortem age (y)") hist(Parvilucina.PVL10.50.bed2[Parvilucina.PVL10.50.bed2<5000], add=T, breaks=seq(0,5000,by=200), col=rgb(1,1,1,0.5)) hist(Nuculana.shallower.50[Nuculana.shallower.50<5000], breaks=seq(0,5000,by=200), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,80), main="Southern California Bight-grabs", cex.main=0.9, xlab="Postmortem age (y)") hist(Parvilucina.PVL10.50.bed3[Parvilucina.PVL10.50.bed3<5000], add=T, breaks=seq(0,5000,by=200), col=rgb(1,1,1,0.5)) ############################################################################################################################################################ #AGE DISTRIBUTIONS OF BOTH SPECIES IN BOTH CORES AT ~8-15-CM SCALE #DOWNCORE AGE DISTRIBUTIONS - NUCULANA ############################################################################################################################################################ Parvilucina.PVL1050.2.6=Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths==2 | Parvilucina.PVL1050.meanSdepths==6] Parvilucina.PVL1050.10.14=Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths==10 | Parvilucina.PVL1050.meanSdepths==14] Parvilucina.PVL1050.18.22=Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths==18 | Parvilucina.PVL1050.meanSdepths==22] Parvilucina.PVL1050.26.34=Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths==26 | Parvilucina.PVL1050.meanSdepths==30 | Parvilucina.PVL1050.meanSdepths==34] Parvilucina.PVL1050.38.41=Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths==38 | Parvilucina.PVL1050.meanSdepths==41.5] Parvilucina.PVL1050.45.49=Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths==45 | Parvilucina.PVL1050.meanSdepths==49] Nuculana.PVL1050.2.6=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==2 | Nuculana.PVL1050.meanSdepths==6] Nuculana.PVL1050.10.14=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==10 | Nuculana.PVL1050.meanSdepths==14] Nuculana.PVL1050.18.22=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==18 | Nuculana.PVL1050.meanSdepths==22] Nuculana.PVL1050.26.30=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==26 | Nuculana.PVL1050.meanSdepths==30] Nuculana.PVL1050.34.38=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==34 | Nuculana.PVL1050.meanSdepths==38] Nuculana.PVL1050.41=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==41.5] Nuculana.PVL1050.45=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==45] Nuculana.PVL1050.49=Nuculana.PVL1050.ages[Nuculana.PVL1050.meanSdepths==49] Nuculana.OC50.2.6=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==2 | Nuculana.OC50.meanSdepths==6] Nuculana.OC50.10.14=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==10 | Nuculana.OC50.meanSdepths==14] Nuculana.OC50.18.22=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==18 | Nuculana.OC50.meanSdepths==22.5] Nuculana.OC50.26.32=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==27.5 | Nuculana.OC50.meanSdepths==32] Nuculana.OC50.36.39=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==35.5 | Nuculana.OC50.meanSdepths==39] Nuculana.OC50.43.47=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==43 | Nuculana.OC50.meanSdepths==47] Nuculana.OC50.51.55=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==51 | Nuculana.OC50.meanSdepths==55] Nuculana.OC50.59.63=Nuculana.OC50.ages[Nuculana.OC50.meanSdepths==59 | Nuculana.OC50.meanSdepths==63] par(mfcol=c(8,3)) par(mar=c(2,4,2,1)) hist(Nuculana.PVL1050.2.6, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="0-8 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.2.6),lwd=2) hist(Nuculana.PVL1050.10.14, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="8-16 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.10.14),lwd=2) hist(Nuculana.PVL1050.18.22, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="16-24 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.18.22),lwd=2) hist(Nuculana.PVL1050.26.30, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="24-32 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.26.30),lwd=2) hist(Nuculana.PVL1050.34.38, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="32-40 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.34.38),lwd=2) hist(Nuculana.PVL1050.41, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="40-43 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.41),lwd=2) hist(Nuculana.PVL1050.45, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="43-47 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.45),lwd=2) hist(Nuculana.PVL1050.49, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="47-51 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.PVL1050.49),lwd=2) hist(Nuculana.OC50.2.6, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="0-8 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.2.6),lwd=2) hist(Nuculana.OC50.10.14, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,25), main="8-16 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.10.14),lwd=2) hist(Nuculana.OC50.18.22, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="16-24 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.18.22),lwd=2) hist(Nuculana.OC50.26.32, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="24-34 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.26.32),lwd=2) hist(Nuculana.OC50.36.39, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="34-41 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.36.39),lwd=2) hist(Nuculana.OC50.43.47, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="43-49 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.43.47),lwd=2) hist(Nuculana.OC50.51.55, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="49-57 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.51.55),lwd=2) hist(Nuculana.OC50.59.63, breaks=seq(0,20000,by=500), col="gray", ylim=c(0,15), main="57-65 cm", cex.main=0.9, xlab="") abline(v=median(Nuculana.OC50.59.63),lwd=2) ############################################################################################################################################################ #INITIALIZE PER-SPECIES PER-INCREMENT DISTRIBUTIONS FOR UNMIXING - LISTS #INTERPOLATE INCREMENTS WITH UNDATED SHELLS WITH RTNORM FUNCTION FOR SUBSURFACE INCREMENTS #REXP FUNCTION FOR SURFACE INCREMENTS IN THE UPPERMOST 20-24 CM ############################################################################################################################################################ Nuculana.PVL1050.list1=split(Nuculana.PVL1050.ages,Nuculana.PVL1050.meanSdepths) Parvilucina.PVL1050.list1=split(Parvilucina.PVL1050.ages,Parvilucina.PVL1050.meanSdepths) ############################################################################################################################################################ #items in the list are stratigraphically ordered, labelled according to the mean depth of the increment in centimeters ############################################################################################################################################################ Nuculana.PVL1050.list1=split(Nuculana.PVL1050.ages,Nuculana.PVL1050.meanSdepths) Nuculana.PVL1050.list1$"2"=Nuculana.PVL1050.list1$"2" Nuculana.PVL1050.list1$"6"=Nuculana.PVL1050.list1$"6" Nuculana.PVL1050.list1$"10"=Nuculana.PVL1050.list1$"10" Nuculana.PVL1050.list1$"14"=Nuculana.PVL1050.list1$"14" Nuculana.PVL1050.list1$"18"=Nuculana.PVL1050.list1$"18" Nuculana.PVL1050.list1$"22"=Nuculana.PVL1050.list1$"22" Nuculana.PVL1050.list1$"26"=Nuculana.PVL1050.list1$"26" Nuculana.PVL1050.list1$"30"=Nuculana.PVL1050.list1$"30" Nuculana.PVL1050.list1$"34"=Nuculana.PVL1050.list1$"34" Nuculana.PVL1050.list1$"38"=Nuculana.PVL1050.list1$"38" Nuculana.PVL1050.list1$"41.5"=Nuculana.PVL1050.list1$"41.5" Nuculana.PVL1050.list1$"45"=Nuculana.PVL1050.list1$"45" Nuculana.PVL1050.list1$"49"=Nuculana.PVL1050.list1$"49" ############################################################################################################################################################ #items in the list are stratigraphically ordered, labelled according to the mean depth of the increment in centimeters ############################################################################################################################################################ Nuculana.OC50.list1=split(Nuculana.OC50.ages,Nuculana.OC50.meanSdepths) Nuculana.OC50.list1$"2"=Nuculana.OC50.list1$"2" Nuculana.OC50.list1$"6"=Nuculana.OC50.list1$"6" Nuculana.OC50.list1$"10"=Nuculana.OC50.list1$"10" Nuculana.OC50.list1$"14"=Nuculana.OC50.list1$"14" Nuculana.OC50.list1$"18"=Nuculana.OC50.list1$"18" Nuculana.OC50.list1$"22.5"=Nuculana.OC50.list1$"22.5" Nuculana.OC50.list1$"27.5"=Nuculana.OC50.list1$"27.5" Nuculana.OC50.list1$"32"=Nuculana.OC50.list1$"32" Nuculana.OC50.list1$"35.5"=Nuculana.OC50.list1$"35.5" Nuculana.OC50.list1$"39"=Nuculana.OC50.list1$"39" Nuculana.OC50.list1$"43"=Nuculana.OC50.list1$"43" Nuculana.OC50.list1$"47"=Nuculana.OC50.list1$"47" Nuculana.OC50.list1$"51"=Nuculana.OC50.list1$"51" Nuculana.OC50.list1$"55"=Nuculana.OC50.list1$"55" Nuculana.OC50.list1$"59"=Nuculana.OC50.list1$"59" Nuculana.OC50.list1$"63"=Nuculana.OC50.list1$"63" Nuculana.OC50.list1$"67"=Nuculana.OC50.list1$"67" Nuculana.OC50.list1$"71"=Nuculana.OC50.list1$"71" Nuculana.OC50.list1$"75"=Nuculana.OC50.list1$"75" Nuculana.OC50.list1$"79"=Nuculana.OC50.list1$"79" Nuculana.OC50.list1$"83"=Nuculana.OC50.list1$"83" Nuculana.OC50.list1$"87"=Nuculana.OC50.list1$"87" ############################################################################################################################################################ #items in the list are stratigraphically ordered, labelled according to the mean depth of the increment in centimeters ############################################################################################################################################################ Parvilucina.PVL1050.list1=split(Parvilucina.PVL1050.ages,Parvilucina.PVL1050.meanSdepths) Parvilucina.PVL1050.list1$"2"=Parvilucina.PVL1050.list1$"2" Parvilucina.PVL1050.list1$"6"=Parvilucina.PVL1050.list1$"6" Parvilucina.PVL1050.list1$"10"=Parvilucina.PVL1050.list1$"10" Parvilucina.PVL1050.list1$"14"=Parvilucina.PVL1050.list1$"14" Parvilucina.PVL1050.list1$"18"=rexp(100,rate=1/mean(c(Parvilucina.PVL1050.list1$"14",Parvilucina.PVL1050.list1$"22"))) Parvilucina.PVL1050.list1$"18"=c(Parvilucina.PVL1050.list1$"14",Parvilucina.PVL1050.list1$"22") Parvilucina.PVL1050.list1$"22"=Parvilucina.PVL1050.list1$"22" Parvilucina.PVL1050.list1$"26"=Parvilucina.PVL1050.list1$"26" Parvilucina.PVL1050.list1$"30"=Parvilucina.PVL1050.list1$"30" Parvilucina.PVL1050.list1$"34"=Parvilucina.PVL1050.list1$"34" Parvilucina.PVL1050.list1$"38"=Parvilucina.PVL1050.list1$"38" Parvilucina.PVL1050.list1$"41.5"=Parvilucina.PVL1050.list1$"41.5" Parvilucina.PVL1050.list1$"45"=Parvilucina.PVL1050.list1$"45" Parvilucina.PVL1050.list1$"49"=Parvilucina.PVL1050.list1$"49" temp=names(Parvilucina.PVL1050.list1) temp=as.numeric(temp) ord=order(temp) Parvilucina.PVL1050.list1=Parvilucina.PVL1050.list1[ord] ############################################################################################################################################################ #items in the list are stratigraphically ordered, labelled according to the mean depth of the increment in centimeters ############################################################################################################################################################ Parvilucina.OC50.list1=split(Parvilucina.OC50.ages,Parvilucina.OC50.meanSdepths) Parvilucina.OC50.list1$"2"=Parvilucina.OC50.list1$"2" Parvilucina.OC50.list1$"6"=Parvilucina.OC50.list1$"6" Parvilucina.OC50.list1$"10"=Parvilucina.OC50.list1$"10" Parvilucina.OC50.list1$"14"=Parvilucina.OC50.list1$"14" Parvilucina.OC50.list1$"18"=Parvilucina.OC50.list1$"18" Parvilucina.OC50.list1$"22.5"=Parvilucina.OC50.list1$"22.5" Parvilucina.OC50.list1$"27.5"=Parvilucina.OC50.list1$"27.5" Parvilucina.OC50.list1$"32"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"27.5")), sd=sd(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"27.5")), lower=0, upper=14500) Parvilucina.OC50.list1$"35.5"=Parvilucina.OC50.list1$"35.5" Parvilucina.OC50.list1$"39"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"51")), sd=sd(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"51")), lower=0, upper=14500) Parvilucina.OC50.list1$"43"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"51")), sd=sd(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"51")), lower=0, upper=14500) Parvilucina.OC50.list1$"47"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"51")), sd=sd(c(Parvilucina.OC50.list1$"35.5",Parvilucina.OC50.list1$"51")), lower=0, upper=14500) Parvilucina.OC50.list1$"51"=Parvilucina.OC50.list1$"51" Parvilucina.OC50.list1$"55"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"51",Parvilucina.OC50.list1$"59")), sd=sd(c(Parvilucina.OC50.list1$"51",Parvilucina.OC50.list1$"59")), lower=0, upper=14500) Parvilucina.OC50.list1$"59"=Parvilucina.OC50.list1$"59" Parvilucina.OC50.list1$"63"=Parvilucina.OC50.list1$"63" Parvilucina.OC50.list1$"67"=Parvilucina.OC50.list1$"67" Parvilucina.OC50.list1$"71"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"67",Parvilucina.OC50.list1$"83")), sd=sd(c(Parvilucina.OC50.list1$"67",Parvilucina.OC50.list1$"83")), lower=0, upper=14500) Parvilucina.OC50.list1$"75"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"67",Parvilucina.OC50.list1$"83")), sd=sd(c(Parvilucina.OC50.list1$"67",Parvilucina.OC50.list1$"83")), lower=0, upper=14500) Parvilucina.OC50.list1$"79"=rtnorm(100,mean=mean(c(Parvilucina.OC50.list1$"67",Parvilucina.OC50.list1$"83")), sd=sd(c(Parvilucina.OC50.list1$"67",Parvilucina.OC50.list1$"83")), lower=0, upper=14500) Parvilucina.OC50.list1$"83"=Parvilucina.OC50.list1$"83" temp=names(Parvilucina.OC50.list1) temp=as.numeric(temp) ord=order(temp) Parvilucina.OC50.list1=Parvilucina.OC50.list1[ord] ############################################################################################################################################################ Parvilucina.PVL10.50.bed1=c(Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths<23]) Parvilucina.PVL10.50.bed2=c(Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths>23 & Parvilucina.PVL1050.meanSdepths<35]) Parvilucina.PVL10.50.bed3=c(Parvilucina.PVL1050.ages[Parvilucina.PVL1050.meanSdepths>35]) Parvilucina.OC.50.bed1=c(Parvilucina.OC50.ages[Parvilucina.OC50.meanSdepths<20]) Parvilucina.OC.50.bed2=c(Parvilucina.OC50.ages[Parvilucina.OC50.meanSdepths>20 & Parvilucina.OC50.meanSdepths<40]) Parvilucina.OC.50.bed3=c(Parvilucina.OC50.ages[Parvilucina.OC50.meanSdepths>40]) ############################################################################################################################################################ #SOUTHERN CALIFORNIA BIGHT #RATES OF DISINTEGRATION BASED ON SURFACE AFDS OF SURFACE SAMPLES (EXCLUDING OC50 and PVL10-50)- PARVILUCINA #FINAL ESTIMATES OF DISINTEGRATION ARE ALSO IN SUPPLEMENTARY TABLE 5 #FOR PARVILUCINA - THREE MODELS DIFFERING IN THE LACKING PRODUCTION FOR 0, 20 and 25 years #REFLECTING THE SHARP DECLINE IN ABUNDANCE OF PARVILUCINA in 1986 AT PALOS VERDES SHELF ############################################################################################################################################################ age=Parvilucina.surface.age hist(age[age<200], breaks=seq(0,300,by=10)) ages=age TIMEMAX=15000 TIME=25 Parvilucina.TIMEMIN=25 #model where disintegration rate does not decline with shell age start1=log(c(0.1,0.001,0.00001)) out=optim(par=start1, OneStateConstantTwiceTruncatedExpMix) out$value LIK=out$value k=5 AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) AICc fit2phasetauUNDERtruncation=exp(out$par[3]) fit2phaselambda2UNDERtruncation=exp(out$par[2]) fit2phaselambda1UNDERtruncation=exp(out$par[1]) r1=fit2phaselambda2UNDERtruncation r2=fit2phaselambda1UNDERtruncation+fit2phasetauUNDERtruncation alpha=fit2phasetauUNDERtruncation/(r2-r1) beta1=((1/r1)*alpha) beta2=((1/r1)*alpha)+((1/r2)*(1-alpha)) T.Parvilucina.beta.surface=beta1/beta2 T.Parvilucina.r1.surface=r1 T.Parvilucina.r2.surface=r2 fit2phasetauUNDERtruncation T.Parvilucina.r2.surface out$value LIK=out$value k=5 AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) AICc T.Parvilucina.lambda1=fit2phaselambda1UNDERtruncation T.Parvilucina.lambda2=fit2phaselambda2UNDERtruncation T.Parvilucina.tau=fit2phasetauUNDERtruncation T.Parvilucina.lambda1 T.Parvilucina.lambda2 T.Parvilucina.tau Parvilucina.onephaserate<-1/mean(age) lkexp=function(rho) {-length(age)*log(rho) + rho*sum(age)} like<--(sum(log(Parvilucina.onephaserate)-Parvilucina.onephaserate*age)) #negative loglikelihood for exponential AIC=2*1-(-2*like) AICc=AIC+((2*1*(1+1))/(length(age)-1-1)) AICc onephaseAIC<-AICc #RANDOM 2-PHASE EXPONENTIAL X=age out<-FitExpMix(X) Parvilucina.randomtwophaserate1.surface<-out$r1 Parvilucina.randomtwophaserate2.surface<-out$r2 Parvilucina.randomtwophasealpha.surface<-out$alpha Parvilucina.randomtwophasebeta.surface<-out$beta lik<--(out$lik) likrandom=lik rAIC=2*3-(-2*lik) randomAIC=rAIC+((2*3*(3+1))/(length(X)-3-1)) randomAIC Parvilucina.tau.surface<-Parvilucina.randomtwophasealpha.surface*(Parvilucina.randomtwophaserate2.surface-Parvilucina.randomtwophaserate1.surface) Parvilucina.lambda1.surface<-Parvilucina.randomtwophaserate1.surface*Parvilucina.randomtwophasealpha.surface+Parvilucina.randomtwophaserate2.surface*(1-Parvilucina.randomtwophasealpha.surface) Parvilucina.lambda2.surface<-Parvilucina.randomtwophaserate1.surface Parvilucina.lambda1.surface Parvilucina.tau.surface Parvilucina.lambda2.surface ############################################################################################################################### #SOUTHERN CALIFORNIA BIGHT #DISINTEGRATION RATES FOR NUCULANA - FROM SITES SHALLOWER THAN 50 M (EXCLUDING OC50 and PVL10-50) #THERE IS NO TRUNCATION FOR NUCULANA (TIME=0) ############################################################################################################################### age=Nuculana.shallower.50 ages=age hist(age[age<200], breaks=seq(0,300,by=10)) TIMEMAX=15000 TIMEMIN=0 TIME=TIMEMIN Nuculana.TIMEMIN=TIMEMIN ########################################################################################### #check a model with a single disintegration rate that does not change with shell age ########################################################################################### out1=optim(par=log(0.0001), OneStateConstantTwiceTruncatedExp) fit1phaseUNDERtruncation=exp(out1$par) lambda=exp(out1$par) lambda LIK=out1$value k=3 AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) AICc ########################################################################################### #a model with disintegration rate declining with shell age ########################################################################################### start1=log(c(0.002,0.0001,0.00001)) out=optim(par=start1, OneStateConstantTwiceTruncatedExpMix) fit2phasetauUNDERtruncation=exp(out$par[3]) fit2phaselambda2UNDERtruncation=exp(out$par[2]) fit2phaselambda1UNDERtruncation=exp(out$par[1]) r1=fit2phaselambda2UNDERtruncation r2=fit2phaselambda1UNDERtruncation+fit2phasetauUNDERtruncation alpha=fit2phasetauUNDERtruncation/(r2-r1) beta1=((1/r1)*alpha) beta2=((1/r1)*alpha)+((1/r2)*(1-alpha)) T.Nuculana.beta.surface=beta1/beta2 T.Nuculana.r1.surface=r1 T.Nuculana.r2.surface=r2 out$value LIK=out$value k=5 AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) AICc T.Nuculana.lambda1=fit2phaselambda1UNDERtruncation T.Nuculana.lambda2=fit2phaselambda2UNDERtruncation T.Nuculana.tau=fit2phasetauUNDERtruncation T.Nuculana.lambda1 T.Nuculana.lambda2 T.Nuculana.tau Nuculana.onephaserate<-1/mean(age) Nuculana.onephaserate lkexp=function(rho) {-length(age)*log(rho) + rho*sum(age)} like<--(sum(log(Nuculana.onephaserate)-Nuculana.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) Nuculana.randomtwophaserate1.surface<-out$r1 Nuculana.randomtwophaserate2.surface<-out$r2 Nuculana.randomtwophasealpha.surface<-out$alpha Nuculana.randomtwophasebeta.surface<-out$beta lik<--(out$lik) likrandom=lik rAIC=2*3-(-2*lik) randomAIC=rAIC+((2*3*(3+1))/(length(X)-3-1)) randomAIC Nuculana.tau.surface<-Nuculana.randomtwophasealpha.surface*(Nuculana.randomtwophaserate2.surface-Nuculana.randomtwophaserate1.surface) Nuculana.lambda1.surface<-Nuculana.randomtwophaserate1.surface*Nuculana.randomtwophasealpha.surface+Nuculana.randomtwophaserate2.surface*(1-Nuculana.randomtwophasealpha.surface) Nuculana.lambda2.surface<-Nuculana.randomtwophaserate1.surface Nuculana.lambda1.surface Nuculana.tau.surface Nuculana.lambda2.surface ############################################################################################################################### #FIGURES SHOWING THE FIT OF MODELS WITH DISINTEGRATION ESTIMATES #FOR SURFACE AGE DISTRIBUTIONS FROM SOUTHERN CALIFORNIA BIGHT (EXCLUDING OC50 and PVL10-50) ############################################################################################################################### ##################################################################### #show the subset of the fit - just to shells younger than 10000 years ##################################################################### par(mfcol=c(3,2)) par(mar=c(4,4,2,1)) out1=hist(Parvilucina.surface.age[Parvilucina.surface.age<10000], breaks=seq(0,10000,by=100), col="gray91", ylim=c(0,300), main="Parvilucina", xlab="") times1=out1$mids-out1$mids[1] predictSIMPLE1<-dexp(times1,rate=Parvilucina.onephaserate) lines(times1, (predictSIMPLE1/sum(predictSIMPLE1))*length(Parvilucina.surface.age), lty=2, col="gray51", lwd=1) predictRANDOM1<-Parvilucina.randomtwophaserate1.surface*exp(-Parvilucina.randomtwophaserate1.surface*times1)*Parvilucina.randomtwophasebeta.surface+ Parvilucina.randomtwophaserate2.surface*exp(-Parvilucina.randomtwophaserate2.surface*times1)*(1-Parvilucina.randomtwophasebeta.surface) lines(times1, predictRANDOM1/sum(predictRANDOM1)*length(Parvilucina.surface.age), lwd=1, col="gray",lty=1) predictTRUNCRANDOM1<-T.Parvilucina.r1.surface*exp(-T.Parvilucina.r1.surface*times1)*T.Parvilucina.beta.surface+T.Parvilucina.r2.surface*exp(-T.Parvilucina.r2.surface*times1)*(1-T.Parvilucina.beta.surface) lines(times1+Parvilucina.TIMEMIN, predictTRUNCRANDOM1/sum(predictTRUNCRANDOM1)*length(Parvilucina.surface.age[Parvilucina.surface.age>Parvilucina.TIMEMIN]), lwd=2, col="black",lty=1) ##################################################################### #show the subset of the fit - just to shells younger than 1000 years ##################################################################### out2=hist(Parvilucina.surface.age[Parvilucina.surface.age<1000], breaks=seq(0,1000,by=10), col="gray91", ylim=c(0,80), main="", xlab="") times2=out2$mids-out2$mids[1] predictSIMPLE2<-dexp(times2,rate=Parvilucina.onephaserate) lines(times2, (predictSIMPLE2/sum(predictSIMPLE2))*length(Parvilucina.surface.age), lty=2, col="gray51", lwd=3) predictRANDOM2<-Parvilucina.randomtwophaserate1.surface*exp(-Parvilucina.randomtwophaserate1.surface*times2)*Parvilucina.randomtwophasebeta.surface+ Parvilucina.randomtwophaserate2.surface*exp(-Parvilucina.randomtwophaserate2.surface*times2)*(1-Parvilucina.randomtwophasebeta.surface) lines(times2, predictRANDOM2/sum(predictRANDOM2)*length(Parvilucina.surface.age), lwd=3, col="gray",lty=1) predictTRUNCRANDOM2<-T.Parvilucina.r1.surface*exp(-T.Parvilucina.r1.surface*times2)*T.Parvilucina.beta.surface+T.Parvilucina.r2.surface*exp(-T.Parvilucina.r2.surface*times2)*(1-T.Parvilucina.beta.surface) age2=Parvilucina.surface.age[Parvilucina.surface.age>Parvilucina.TIMEMIN] lines(times2[1:100]+Parvilucina.TIMEMIN, predictTRUNCRANDOM2/sum(predictTRUNCRANDOM2)*length(age2), lwd=3, col="black",lty=1) ##################################################################### #show the subset of the fit - just to shells younger than 200 years ##################################################################### out3=hist(Parvilucina.surface.age[Parvilucina.surface.age<200], breaks=seq(0,200,by=10), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age (y)") times3=out3$mids-out3$mids[1] predictSIMPLE3<-dexp(times3,rate=Parvilucina.onephaserate) lines(times3, (predictSIMPLE3/sum(predictSIMPLE3))*length(Parvilucina.surface.age), lty=2, col="gray51", lwd=3) predictRANDOM3<-Parvilucina.randomtwophaserate1.surface*exp(-Parvilucina.randomtwophaserate1.surface*times3)*Parvilucina.randomtwophasebeta.surface+ Parvilucina.randomtwophaserate2.surface*exp(-Parvilucina.randomtwophaserate2.surface*times3)*(1-Parvilucina.randomtwophasebeta.surface) lines(times3, predictRANDOM3/sum(predictRANDOM3)*length(Parvilucina.surface.age), lwd=3, col="gray",lty=1) predictTRUNCRANDOM3<-T.Parvilucina.r1.surface*exp(-T.Parvilucina.r1.surface*times3)*T.Parvilucina.beta.surface+T.Parvilucina.r2.surface*exp(-T.Parvilucina.r2.surface*times3)*(1-T.Parvilucina.beta.surface) legend(x="topright", legend=c("Time-constant model","Rate declines with shell age","Rate declines with shell age-truncated"), lwd=3, lty=c(2,1,1), col=c("gray51","gray","black")) lines(times3[1:20]+Parvilucina.TIMEMIN, predictTRUNCRANDOM3[1:20]/sum(predictTRUNCRANDOM3[1:20])*length(age2), lwd=3, col="black",lty=1) ##################################################################### #show the subset of the fit - just to shells younger than 10000 years ##################################################################### out1=hist(Nuculana.shallower.50[Nuculana.shallower.50<10000], breaks=seq(0,10000,by=100), col="gray91", ylim=c(0,120), main="Nuculana", xlab="") times1=out1$mids-out1$mids[1] predictSIMPLE1<-dexp(times1,rate=Nuculana.onephaserate) lines(times1, (predictSIMPLE1/sum(predictSIMPLE1))*length(Nuculana.shallower.50), lty=1, col="gray51", lwd=1) predictRANDOM1<-Nuculana.randomtwophaserate1.surface*exp(-Nuculana.randomtwophaserate1.surface*times1)*Nuculana.randomtwophasebeta.surface+ Nuculana.randomtwophaserate2.surface*exp(-Nuculana.randomtwophaserate2.surface*times1)*(1-Nuculana.randomtwophasebeta.surface) lines(times1, predictRANDOM1/sum(predictRANDOM1)*length(Nuculana.shallower.50), lwd=2, col="black",lty=1) predictTRUNCRANDOM1<-T.Nuculana.r1.surface*exp(-T.Nuculana.r1.surface*times1)*T.Nuculana.beta.surface+T.Nuculana.r2.surface*exp(-T.Nuculana.r2.surface*times1)*(1-T.Nuculana.beta.surface) lines(times1, predictTRUNCRANDOM1/sum(predictTRUNCRANDOM1)*length(Nuculana.shallower.50), lwd=1, col="black",lty=1) ##################################################################### #show the subset of the fit - just to shells younger than 1000 years ##################################################################### out2=hist(Nuculana.shallower.50[Nuculana.shallower.50<1000], breaks=seq(0,1000,by=10), col="gray91", ylim=c(0,100), main="", xlab="") times2=out2$mids-out2$mids[1] predictSIMPLE2<-dexp(times2,rate=Nuculana.onephaserate) lines(times2, (predictSIMPLE2/sum(predictSIMPLE2))*length(Nuculana.shallower.50), lty=2, col="gray51", lwd=3) predictRANDOM2<-Nuculana.randomtwophaserate1.surface*exp(-Nuculana.randomtwophaserate1.surface*times2)*Nuculana.randomtwophasebeta.surface+ Nuculana.randomtwophaserate2.surface*exp(-Nuculana.randomtwophaserate2.surface*times2)*(1-Nuculana.randomtwophasebeta.surface) lines(times2, predictRANDOM2/sum(predictRANDOM2)*length(Nuculana.shallower.50), lwd=3, col="gray",lty=1) predictTRUNCRANDOM2<-T.Nuculana.r1.surface*exp(-T.Nuculana.r1.surface*times2)*T.Nuculana.beta.surface+T.Nuculana.r2.surface*exp(-T.Nuculana.r2.surface*times2)*(1-T.Nuculana.beta.surface) age2=Nuculana.shallower.50[Nuculana.shallower.50>Nuculana.TIMEMIN] lines(times2, predictTRUNCRANDOM2/sum(predictTRUNCRANDOM2)*length(age2), lwd=3, col="black",lty=1) ##################################################################### #show the subset of the fit - just to shells younger than 200 years ##################################################################### out3=hist(Nuculana.shallower.50[Nuculana.shallower.50<200], breaks=seq(0,200,by=10), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age(y)") times3=out3$mids-out3$mids[1] predictSIMPLE3<-dexp(times3,rate=Nuculana.onephaserate) lines(times3, (predictSIMPLE3/sum(predictSIMPLE3))*length(Nuculana.shallower.50), lty=2, col="gray51", lwd=3) predictRANDOM3<-Nuculana.randomtwophaserate1.surface*exp(-Nuculana.randomtwophaserate1.surface*times3)*Nuculana.randomtwophasebeta.surface+ Nuculana.randomtwophaserate2.surface*exp(-Nuculana.randomtwophaserate2.surface*times3)*(1-Nuculana.randomtwophasebeta.surface) lines(times3, predictRANDOM3/sum(predictRANDOM3)*length(Nuculana.shallower.50), lwd=3, col="gray",lty=1) predictTRUNCRANDOM3<-T.Nuculana.r1.surface*exp(-T.Nuculana.r1.surface*times3)*T.Nuculana.beta.surface+T.Nuculana.r2.surface*exp(-T.Nuculana.r2.surface*times3)*(1-T.Nuculana.beta.surface) lines(times3, predictTRUNCRANDOM3/sum(predictTRUNCRANDOM3)*length(age2), lwd=3, col="black",lty=1) ############################################################################################################################### #NUCULANA - UNCERTAINTY IN DISINTEGRATION ESTIMATION DUE TO CALIBRATION ERRORS #SOUTHERN CALIFORNIA BIGHT (EXCLUDING OC50 and PVL10-50) ############################################################################################################################### out3=hist(Nuculana.shallower.50[Nuculana.shallower.50<200], breaks=seq(0,200,by=10), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age(y)") times3=out3$mids-out3$mids[1] out4=hist(Nuculana.shallower.50[Nuculana.shallower.50<10000], breaks=seq(0,10000,by=100), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age(y)") times4=out4$mids-out4$mids[1] #Nuculana boot.T.Nuculana.lambda1=numeric() boot.T.Nuculana.lambda2=numeric() boot.T.Nuculana.tau=numeric() boot.T.Nuculana.alpha=numeric() boot.T.Nuculana.beta=numeric() boot.T.Nuculana.r2=numeric() boot.T.Nuculana.r1=numeric() predictRANDOM.Nuculana.sim=array(NA, dim=c(100,length(times3))) predictRANDOM.Nuculana.sim.10000=array(NA, dim=c(100,length(times4))) for (i in 1:100) { sample.age=Nuculana.shallower.50 #sample shell ages from the lognormal distribution (shell ages represent medians of posterior distributions, not means), so the mean of the lognormal just corresponds to the #log-transformed median age estimate, using Nuculana.ln.d from the Bayesian calibration #this step may be skipped bootages=rtlnorm(length(sample.age),log(sample.age),sdlog=sqrt(exp(Nuculana.ln.d)),a=0,b=20000) bootages=sample(bootages,length(bootages),replace=T) #just bootstrapping shell ages #bootages=sample(sample.age, length(sample.age), replace=T) hist(bootages[bootages<5000], breaks=seq(0,5000,by=50)) ages=bootages TIMEMAX=15000 TIMEMIN=0 TIME=TIMEMIN start1=log(c(0.02,0.0001,0.00001)) out=optim(par=start1, OneStateConstantTwiceTruncatedExpMix) out$value; LIK=out$value; k=5 fit2phasetauUNDERtruncation=exp(out$par[3]) fit2phaselambda2UNDERtruncation=exp(out$par[2]) fit2phaselambda1UNDERtruncation=exp(out$par[1]) r1=fit2phaselambda2UNDERtruncation r2=fit2phaselambda1UNDERtruncation+fit2phasetauUNDERtruncation boot.T.Nuculana.r2[i]=r2 boot.T.Nuculana.lambda1[i]=exp(out$par[1]) boot.T.Nuculana.tau[i]=exp(out$par[3]) boot.T.Nuculana.alpha[i]=fit2phasetauUNDERtruncation/(r2-r1) boot.T.Nuculana.lambda2[i]=r1 beta1=((1/r1)*boot.T.Nuculana.alpha[i]) beta2=((1/r1)*boot.T.Nuculana.alpha[i])+((1/r2)*(1-boot.T.Nuculana.alpha[i])) beta=beta1/beta2 boot.T.Nuculana.r1[i]=boot.T.Nuculana.lambda2[i] boot.T.Nuculana.beta[i]=beta predictRANDOM.Nuculana.sim[i,]=r1*exp(-r1*times3)*beta+r2*exp(-r2*times3)*(1-beta) predictRANDOM.Nuculana.sim.10000[i,]=r1*exp(-r1*times4)*beta+r2*exp(-r2*times4)*(1-beta) } predict.Nuculana.RANDOM.mean=apply(predictRANDOM.Nuculana.sim,2,mean) predict.Nuculana.RANDOM.lci=apply(predictRANDOM.Nuculana.sim,2,quantile,0.025) predict.Nuculana.RANDOM.uci=apply(predictRANDOM.Nuculana.sim,2,quantile,0.975) predict.Nuculana.RANDOM.mean.10000=apply(predictRANDOM.Nuculana.sim.10000,2,mean) predict.Nuculana.RANDOM.lci.10000=apply(predictRANDOM.Nuculana.sim.10000,2,quantile,0.025) predict.Nuculana.RANDOM.uci.10000=apply(predictRANDOM.Nuculana.sim.10000,2,quantile,0.975) quantile(boot.T.Nuculana.lambda1, c(0.025,0.5,0.975)) quantile(boot.T.Nuculana.tau, c(0.025,0.5,0.975)) quantile(boot.T.Nuculana.lambda2, c(0.025,0.5,0.975)) out3=hist(Nuculana.shallower.50[Nuculana.shallower.50<10000], breaks=seq(0,10000,by=100), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age(y)") times3=out3$mids-out3$mids[1] #predictSIMPLE3<-dexp(times3,rate=Nuculana.onephaserate) #lines(times3, (predictSIMPLE3/sum(predictSIMPLE3))*length(Nuculana.shallower.50), lty=2, col="gray51", lwd=3) lines(times3, predict.Nuculana.RANDOM.mean.10000/sum(predict.Nuculana.RANDOM.mean.10000)*length(Nuculana.shallower.50), lwd=2, col="gray",lty=1) lines(times3, predict.Nuculana.RANDOM.lci.10000/sum(predict.Nuculana.RANDOM.mean.10000)*length(Nuculana.shallower.50), lwd=1, col="black",lty=2) lines(times3, predict.Nuculana.RANDOM.uci.10000/sum(predict.Nuculana.RANDOM.mean.10000)*length(Nuculana.shallower.50), lwd=1, col="black",lty=2) predict.Nuculana.RANDOM3<-Nuculana.randomtwophaserate1.surface*exp(-Nuculana.randomtwophaserate1.surface*times3)*Nuculana.randomtwophasebeta.surface+ Nuculana.randomtwophaserate2.surface*exp(-Nuculana.randomtwophaserate2.surface*times3)*(1-Nuculana.randomtwophasebeta.surface) out3=hist(Nuculana.shallower.50[Nuculana.shallower.50<200], breaks=seq(0,200,by=10), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age(y)") times3=out3$mids-out3$mids[1] predictSIMPLE3<-dexp(times3,rate=Nuculana.onephaserate) lines(times3, (predictSIMPLE3/sum(predictSIMPLE3))*length(Nuculana.shallower.50), lty=2, col="gray51", lwd=3) lines(times3, predict.Nuculana.RANDOM.mean/sum(predict.Nuculana.RANDOM.mean)*length(Nuculana.shallower.50), lwd=2, col="gray",lty=1) lines(times3, predict.Nuculana.RANDOM.lci/sum(predict.Nuculana.RANDOM.mean)*length(Nuculana.shallower.50), lwd=1, col="black",lty=2) lines(times3, predict.Nuculana.RANDOM.uci/sum(predict.Nuculana.RANDOM.mean)*length(Nuculana.shallower.50), lwd=1, col="black",lty=2) ############################################################################################################################### #ESTIMATE DISINTEGRATION IN THE MIXED LAYER (TAPHONOMIC ACTIVE ZONE) OF PARVILUCINA, ASSUMING A SIGNIFICANT DECLINE #IN PRODUCTION OVER THE PAST 20 YEARS, 25 YEARS or 0 YEARS - variable TIME BELOW #PARVILUCINA UNCERTAINTY IN PARAMETER ESTIMATION DUE TO CALIBRATION ERRORS #SOUTHERN CALIFORNIA BIGHT (EXCLUDING OC50 and PVL10-50) ############################################################################################################################### par(mfrow=c(2,3)) boot.T.Parvilucina.lambda2=numeric() boot.T.Parvilucina.lambda1=numeric() boot.T.Parvilucina.tau=numeric() boot.T.Parvilucina.alpha=numeric() boot.T.Parvilucina.r2=numeric() boot.T.Parvilucina.r1=numeric() boot.T.Parvilucina.beta=numeric() boot.T.Parvilucina.lambda1.03=numeric() boot.T.Parvilucina.tau.03=numeric() boot.T.Parvilucina.lambda1.02=numeric() boot.T.Parvilucina.tau.02=numeric() boot.T.Parvilucina.lambda1.01=numeric() boot.T.Parvilucina.tau.01=numeric() boot.T.Parvilucina.tau.03[i]=exp(out$par[3]) boot.T.Parvilucina.lambda2.01=numeric() boot.T.Parvilucina.lambda2.02=numeric() boot.T.Parvilucina.lambda2.03=numeric() predictRANDOM01.Parvilucina.sim=array(NA, dim=c(1000,length(times3))) predictRANDOM02.Parvilucina.sim=array(NA, dim=c(1000,length(times3))) predictRANDOM03.Parvilucina.sim=array(NA, dim=c(1000,length(times3))) predictRANDOM01.Parvilucina.sim.10000=array(NA, dim=c(1000,length(times4))) predictRANDOM02.Parvilucina.sim.10000=array(NA, dim=c(1000,length(times4))) predictRANDOM03.Parvilucina.sim.10000=array(NA, dim=c(1000,length(times4))) minimum=numeric() for (i in 1:1000) { sample.age=Parvilucina.surface.age bootages=sample(sample.age, length(sample.age), replace=T) hist(bootages[bootages<500], breaks=seq(0,500,by=5)) ages=bootages TIMEMAX=10000 TIME=25 minimum[i]=TIME start1=log(0.01) out=optim(par=start1,OneStateConstantTwiceTruncatedExp) out$value; LIK=out$value; k=3 exp(out$par[1]) AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) Parvilucina.l1.T25=exp(out$par[1]) Parvilucina.l1.T25.AIC=AICc start1=log(c(0.5,0.001,0.0001)) out=optim(par=start1, OneStateConstantTwiceTruncatedExpMix) out$value; LIK=out$value; k=5 AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) Parvilucina.T25.AIC=AICc fit2phasetauUNDERtruncation=exp(out$par[3]) fit2phaselambda2UNDERtruncation=exp(out$par[2]) fit2phaselambda1UNDERtruncation=exp(out$par[1]) r1=fit2phaselambda2UNDERtruncation r2=fit2phaselambda1UNDERtruncation+fit2phasetauUNDERtruncation boot.T.Parvilucina.r2[i]=r2 boot.T.Parvilucina.lambda1[i]=exp(out$par[1]) boot.T.Parvilucina.tau[i]=exp(out$par[3]) boot.T.Parvilucina.alpha[i]=fit2phasetauUNDERtruncation/(r2-r1) boot.T.Parvilucina.lambda2[i]=r1 beta1=((1/r1)*boot.T.Parvilucina.alpha[i]) beta2=((1/r1)*boot.T.Parvilucina.alpha[i])+((1/r2)*(1-boot.T.Parvilucina.alpha[i])) beta=beta1/beta2 boot.T.Parvilucina.r1[i]=boot.T.Parvilucina.lambda2[i] boot.T.Parvilucina.beta[i]=beta predictRANDOM03.Parvilucina.sim[i,]=r1*exp(-r1*times3)*beta+r2*exp(-r2*times3)*(1-beta) predictRANDOM03.Parvilucina.sim.10000[i,]=r1*exp(-r1*times4)*beta+r2*exp(-r2*times4)*(1-beta) boot.T.Parvilucina.lambda1.03[i]=exp(out$par[1]) boot.T.Parvilucina.tau.03[i]=exp(out$par[3]) boot.T.Parvilucina.lambda2.03[i]=exp(out$par[2]) TIME=20 minimum[i]=TIME start1=log(0.01) out=optim(par=start1,OneStateConstantTwiceTruncatedExp) out$value; LIK=out$value; k=3 exp(out$par[1]) AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) Parvilucina.l1.T20.AIC=AICc Parvilucina.l1.T20=exp(out$par[1]) start1=log(c(0.5,0.001,0.0001)) out=optim(par=start1, OneStateConstantTwiceTruncatedExpMix) out$value; LIK=out$value; k=5 AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) Parvilucina.T20.AIC=AICc fit2phasetauUNDERtruncation=exp(out$par[3]) fit2phaselambda2UNDERtruncation=exp(out$par[2]) fit2phaselambda1UNDERtruncation=exp(out$par[1]) r1=fit2phaselambda2UNDERtruncation r2=fit2phaselambda1UNDERtruncation+fit2phasetauUNDERtruncation boot.T.Parvilucina.r2[i]=r2 boot.T.Parvilucina.lambda1[i]=exp(out$par[1]) boot.T.Parvilucina.tau[i]=exp(out$par[3]) boot.T.Parvilucina.alpha[i]=fit2phasetauUNDERtruncation/(r2-r1) boot.T.Parvilucina.lambda2[i]=r1 beta1=((1/r1)*boot.T.Parvilucina.alpha[i]) beta2=((1/r1)*boot.T.Parvilucina.alpha[i])+((1/r2)*(1-boot.T.Parvilucina.alpha[i])) beta=beta1/beta2 boot.T.Parvilucina.r1[i]=boot.T.Parvilucina.lambda2[i] boot.T.Parvilucina.beta[i]=beta predictRANDOM02.Parvilucina.sim[i,]=r1*exp(-r1*times3)*beta+r2*exp(-r2*times3)*(1-beta) predictRANDOM02.Parvilucina.sim.10000[i,]=r1*exp(-r1*times4)*beta+r2*exp(-r2*times4)*(1-beta) boot.T.Parvilucina.lambda1.02[i]=exp(out$par[1]) boot.T.Parvilucina.tau.02[i]=exp(out$par[3]) boot.T.Parvilucina.lambda2.02[i]=exp(out$par[2]) TIME=0 start1=log(0.01) out=optim(par=start1,OneStateConstantTwiceTruncatedExp) out$value; LIK=out$value; k=3 exp(out$par[1]) AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) Parvilucina.l1.T0.AIC=AICc Parvilucina.l1.T0=exp(out$par[1]) start1=log(c(0.5,0.001,0.0001)) out=optim(par=start1, OneStateConstantTwiceTruncatedExpMix) out$value; LIK=out$value; k=5 AIC=k*1-(-2*LIK) AICc=AIC+((2*k*(k+1))/(length(age)-k-1)) Parvilucina.T0.AIC=AICc fit2phasetauUNDERtruncation=exp(out$par[3]) fit2phaselambda2UNDERtruncation=exp(out$par[2]) fit2phaselambda1UNDERtruncation=exp(out$par[1]) r1=fit2phaselambda2UNDERtruncation r2=fit2phaselambda1UNDERtruncation+fit2phasetauUNDERtruncation boot.T.Parvilucina.r2[i]=r2 boot.T.Parvilucina.lambda1[i]=exp(out$par[1]) boot.T.Parvilucina.tau[i]=exp(out$par[3]) boot.T.Parvilucina.alpha[i]=fit2phasetauUNDERtruncation/(r2-r1) boot.T.Parvilucina.lambda2[i]=r1 beta1=((1/r1)*boot.T.Parvilucina.alpha[i]) beta2=((1/r1)*boot.T.Parvilucina.alpha[i])+((1/r2)*(1-boot.T.Parvilucina.alpha[i])) beta=beta1/beta2 boot.T.Parvilucina.r1[i]=boot.T.Parvilucina.lambda2[i] boot.T.Parvilucina.beta[i]=beta predictRANDOM01.Parvilucina.sim[i,]=r1*exp(-r1*times3)*beta+r2*exp(-r2*times3)*(1-beta) predictRANDOM01.Parvilucina.sim.10000[i,]=r1*exp(-r1*times4)*beta+r2*exp(-r2*times4)*(1-beta) boot.T.Parvilucina.lambda1.01[i]=exp(out$par[1]) boot.T.Parvilucina.tau.01[i]=exp(out$par[3]) boot.T.Parvilucina.lambda2.01[i]=exp(out$par[2]) } ######################################################################### #CHECK PARAMETER ESTIMATES - USED IN UNMIXING ######################################################################### Parvilucina.l1.T20 Parvilucina.l1.T20.AIC Parvilucina.l1.T25 Parvilucina.l1.T25.AIC Parvilucina.l1.T0 Parvilucina.l1.T0.AIC #no truncation quantile(boot.T.Parvilucina.lambda1.01, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.tau.01, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.lambda2.01, c(0.025,0.5,0.975)) #truncation=20 years quantile(boot.T.Parvilucina.lambda1.02, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.tau.02, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.lambda2.02, c(0.025,0.5,0.975)) #truncation=25 years quantile(boot.T.Parvilucina.lambda1.03, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.tau.03, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.lambda2.03, c(0.025,0.5,0.975)) predict.Parvilucina.RANDOM01.mean=apply(predictRANDOM01.Parvilucina.sim,2,mean) predict.Parvilucina.RANDOM01.lci=apply(predictRANDOM01.Parvilucina.sim,2,quantile,0.025) predict.Parvilucina.RANDOM01.uci=apply(predictRANDOM01.Parvilucina.sim,2,quantile,0.975) predict.Parvilucina.RANDOM02.mean=apply(predictRANDOM02.Parvilucina.sim,2,mean) predict.Parvilucina.RANDOM02.lci=apply(predictRANDOM02.Parvilucina.sim,2,quantile,0.025) predict.Parvilucina.RANDOM02.uci=apply(predictRANDOM02.Parvilucina.sim,2,quantile,0.975) predict.Parvilucina.RANDOM03.mean=apply(predictRANDOM03.Parvilucina.sim,2,mean, na.rm=T) predict.Parvilucina.RANDOM03.lci=apply(predictRANDOM03.Parvilucina.sim,2,quantile,0.025, na.rm=T) predict.Parvilucina.RANDOM03.uci=apply(predictRANDOM03.Parvilucina.sim,2,quantile,0.975, na.rm=T) out3=hist(Parvilucina.surface.age[Parvilucina.surface.age<200], breaks=seq(0,200,by=10), col="gray91", ylim=c(0,80), main="", xlab="", ylab="") times3=out3$mids-out3$mids[1] lines(times3, predict.Parvilucina.RANDOM01.mean/sum(predict.Parvilucina.RANDOM01.mean)*length(Parvilucina.surface.age), lwd=2, col="gray",lty=1) lines(times3, predict.Parvilucina.RANDOM01.lci/sum(predict.Parvilucina.RANDOM01.mean)*length(Parvilucina.surface.age), lwd=1, col="black",lty=2) lines(times3, predict.Parvilucina.RANDOM01.uci/sum(predict.Parvilucina.RANDOM01.mean)*length(Parvilucina.surface.age), lwd=1, col="black",lty=2) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.mean[3:length(times3)]/sum(predict.Parvilucina.RANDOM02.mean)*length(Parvilucina.surface.age), lwd=2, col="black",lty=1) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.lci[3:length(times3)]/sum(predict.Parvilucina.RANDOM02.mean)*length(Parvilucina.surface.age), lwd=1, col="black",lty=2) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.uci[3:length(times3)]/sum(predict.Parvilucina.RANDOM02.mean)*length(Parvilucina.surface.age), lwd=1, col="black",lty=2) out3=hist(Nuculana.shallower.50[Nuculana.shallower.50<200], breaks=seq(0,200,by=10), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age(y)") times3=out3$mids-out3$mids[1] predictSIMPLE3<-dexp(times3,rate=Nuculana.onephaserate) lines(times3, (predictSIMPLE3/sum(predictSIMPLE3))*length(Nuculana.shallower.50), lty=2, col="gray51", lwd=3) lines(times3, predict.Nuculana.RANDOM.mean/sum(predict.Nuculana.RANDOM.mean)*length(Nuculana.shallower.50), lwd=2, col="gray",lty=1) lines(times3, predict.Nuculana.RANDOM.lci/sum(predict.Nuculana.RANDOM.mean)*length(Nuculana.shallower.50), lwd=1, col="black",lty=2) lines(times3, predict.Nuculana.RANDOM.uci/sum(predict.Nuculana.RANDOM.mean)*length(Nuculana.shallower.50), lwd=1, col="black",lty=2) ########################################################################################################################################## #FIGURE 08 MAIN TEXT IN PP2019 - SURFACE AFDS - SCB AND PVL10-50 AND OC-50 - 3 ROWS TIMES 3 COLUMNS ########################################################################################################################################## #cairo_pdf(file='plot.pdf', width=7, height=7) par(mfrow=c(3,3)) par(mar=c(2,2,1,0)) hist(Nuculana.PVL10.50.bed1, breaks=seq(0,15000,by=250), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="PVL10-50-N. taphria", cex.main=0.9, xlab="", cex.axis=1, cex.lab=1,ylab="Number of specimens") hist(Parvilucina.PVL10.50.bed1, breaks=seq(0,15000,by=250), col=rgb(1,1,1,0.5), ylim=c(0,80), main="PVL10-50-P. tenuisculpta", cex.main=0.9, xlab="", cex.axis=1) hist(Nuculana.PVL10.50.bed1[Nuculana.PVL10.50.bed1<210], breaks=seq(0,210,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="18-20th c.", cex.main=0.9, xlab="", cex.axis=1) hist(Parvilucina.PVL10.50.bed1[Parvilucina.PVL10.50.bed1<210], add=T, breaks=seq(0,210,by=10), col=rgb(1,1,1,0.5)) legend(x=75, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1) hist(Nuculana.OC.50.bed1, breaks=seq(0,15000,by=250), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="OC-50-N. taphria", cex.main=0.9, cex.axis=1, xlab="Postmortem age (y)", cex.lab=1, ylab="Number of specimens") hist(Parvilucina.OC.50.bed1, breaks=seq(0,15000,by=250), col=rgb(1,1,1,0.5), ylim=c(0,80), main="OC-50-P. tenuisculpta", cex.main=0.9, cex.axis=1, xlab="Postmortem age (y)", cex.lab=1) hist(Nuculana.OC.50.bed1[Nuculana.OC.50.bed1<210], breaks=seq(0,210,by=10), col=rgb(190/255,190/255,190/255,0.5), ylim=c(0,40), main="18-20th c.", cex.main=0.9, cex.axis=1, xlab="Postmortem age (y)", cex.lab=1) hist(Parvilucina.OC.50.bed1[Parvilucina.OC.50.bed1<210], add=T, breaks=seq(0,210,by=10), col=rgb(1,1,1,0.5)) legend(x=75, y=35, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1) out4=hist(Nuculana.surface.age[Nuculana.surface.age<15000], breaks=seq(0,15000,by=250), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age (y)") times4=out4$mids-out4$mids[1] vp <- baseViewports() pushViewport(vp$inner,vp$figure,vp$plot) pushViewport(viewport(x=0.25,y=0.99,width=0.7,height=0.7,just=c("left","top"))) par(plt=gridPLT(),new=T) out3=hist(Nuculana.surface.age[Nuculana.surface.age<200], breaks=seq(0,200,by=10), col="gray91", main="SCB-N. taphria", cex.main=0.9, ylab="",xlab="", freq=F) times3=out3$mids-out3$mids[1] lines(times3, predict.Nuculana.RANDOM.mean, lwd=2, col="gray",lty=1) lines(times3, predict.Nuculana.RANDOM.lci, lwd=1, col="black",lty=2) lines(times3, predict.Nuculana.RANDOM.uci, lwd=1, col="black",lty=2) popViewport(4) par(mar=c(2,2,1,0)) out4=hist(Parvilucina.surface.age[Parvilucina.surface.age<15000], breaks=seq(0,15000,by=250), col="gray91", ylim=c(0,200), main="", xlab="Postmortem age (y)") times4=out4$mids-out4$mids[1] vp <- baseViewports() pushViewport(vp$inner,vp$figure,vp$plot) pushViewport(viewport(x=0.25,y=0.99,width=0.7,height=0.7,just=c("left","top"))) par(plt=gridPLT(),new=T) out3=hist(Parvilucina.surface.age[Parvilucina.surface.age<200], breaks=seq(0,200,by=10), col="gray91", main="SCB-P. tenuisculpta", cex.main=0.9, ylab="",xlab="", frame.plot=F, freq=FALSE) times3=out3$mids-out3$mids[1] #RANDOM01 = Tmin is 0 years lines(times3, predict.Parvilucina.RANDOM01.mean, lwd=2, col="gray",lty=1) lines(times3, predict.Parvilucina.RANDOM01.lci, lwd=1, col="black",lty=2) lines(times3, predict.Parvilucina.RANDOM01.uci, lwd=1, col="black",lty=2) #RANDOM02 = Tmin is 20 years #lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.mean[3:length(times3)], lwd=2, col="black",lty=1) #lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.lci[3:length(times3)], lwd=1, col="black",lty=2) #lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.uci[3:length(times3)], lwd=1, col="black",lty=2) #RANDOM03 = Tmin is 25 years lines(times3[3:length(times3)], predict.Parvilucina.RANDOM03.mean[3:length(times3)], lwd=2, col="black",lty=1) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM03.lci[3:length(times3)], lwd=1, col="black",lty=2) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM03.uci[3:length(times3)], lwd=1, col="black",lty=2) popViewport(4) par(mar=c(2,2,1,0)) hist(Parvilucina.surface.age[Parvilucina.surface.age<210], breaks=seq(0,210,by=10), col=rgb(1,1,1,0.5), ylim=c(0,100), main="18-20th c.", cex.main=0.9, xlab="", cex.axis=1) hist(Nuculana.surface.age[Nuculana.surface.age<210], add=T, breaks=seq(0,210,by=10), col=rgb(190/255,190/255,190/255,0.5),xlab="",main="") legend(x=75, y=85, legend=c("Parvilucina","Nuculana"), pch=22, pt.bg=c("white","gray"), bty="n", cex=1) #dev.off() ########################################################################################################################################## #SUMMARIZING FITS FOR PARVILUCINA AND NUCULANA #BASED ON VAN VEEN GRAB SAMPLES - ~UPPER 20 CM in the TAPHONOMIC ACTIVE ZONE ########################################################################################################################################## par(mfrow=c(2,2)) par(mar=c(4,4,2,1)) out4=hist(Parvilucina.surface.age[Parvilucina.surface.age<10000], breaks=seq(0,10000,by=100), col="gray91", ylim=c(0,200), main="", xlab="Postmortem age (y)") times4=out4$mids-out4$mids[1] vp <- baseViewports() pushViewport(vp$inner,vp$figure,vp$plot) pushViewport(viewport(x=0.25,y=0.99,width=0.7,height=0.7,just=c("left","top"))) par(plt=gridPLT(),new=T) out3=hist(Parvilucina.surface.age[Parvilucina.surface.age<200], breaks=seq(0,200,by=10), col="gray91", main="", ylab="",xlab="", frame.plot=F, freq=FALSE) times3=out3$mids-out3$mids[1] lines(times3, predict.Parvilucina.RANDOM01.mean, lwd=2, col="gray",lty=1) lines(times3, predict.Parvilucina.RANDOM01.lci, lwd=1, col="black",lty=2) lines(times3, predict.Parvilucina.RANDOM01.uci, lwd=1, col="black",lty=2) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.mean[3:length(times3)], lwd=2, col="black",lty=1) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.lci[3:length(times3)], lwd=1, col="black",lty=2) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM02.uci[3:length(times3)], lwd=1, col="black",lty=2) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM03.mean[3:length(times3)], lwd=2, col="black",lty=1) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM03.lci[3:length(times3)], lwd=1, col="black",lty=2) lines(times3[3:length(times3)], predict.Parvilucina.RANDOM03.uci[3:length(times3)], lwd=1, col="black",lty=2) popViewport(4) par(mar=c(4,4,2,1)) out4=hist(Nuculana.surface.age[Nuculana.surface.age<10000], breaks=seq(0,10000,by=100), col="gray91", ylim=c(0,80), main="", xlab="Postmortem age (y)") times4=out4$mids-out4$mids[1] vp <- baseViewports() pushViewport(vp$inner,vp$figure,vp$plot) pushViewport(viewport(x=0.25,y=0.99,width=0.7,height=0.7,just=c("left","top"))) par(plt=gridPLT(),new=T) out3=hist(Nuculana.surface.age[Nuculana.surface.age<200], breaks=seq(0,200,by=10), col="gray91", main="", ylab="",xlab="", frame.plot=F, freq=F) times3=out3$mids-out3$mids[1] lines(times3, predict.Nuculana.RANDOM.mean, lwd=2, col="gray",lty=1) lines(times3, predict.Nuculana.RANDOM.lci, lwd=1, col="black",lty=2) lines(times3, predict.Nuculana.RANDOM.uci, lwd=1, col="black",lty=2) popViewport(4) par(mar=c(4,4,2,1)) ####################################################################################################### #UNMIXING RECONSTRUCTION OF PARVILUCINA TRAJECTORY (DISINTEGRATION NOT ACCOUNTED FOR) ####################################################################################################### ##################################################################### #absolute Parvilucina total abundance at PVL10-50 ##################################################################### totalsize=c(49,105,194,200,171,135,378,84,126,147,441,588,588) FACTOR=exp(Parvilucina.ln.d) Full.Parvilucina.PVL1050.bootout.250=array(NA, dim=c(100,120)) Full.Parvilucina.PVL1050.bootout.10=array(NA, dim=c(100,3000)) for (j in 1:100) { temp.bootages=numeric() for (i in 1:length(Parvilucina.PVL1050.list1)) { temp=unlist(Parvilucina.PVL1050.list1[i]) sample.age=sample(temp,totalsize[i],replace=T) temp.bootages=append(temp.bootages,sample.age,length(temp.bootages)) } #sample shell ages from the gamma distribution of age errors picked by the Bayesian model (shell ages represent medians of posterior distributions, not means), #using Parvilucina.ln.d from the Bayesian calibration #this step may be skipped - not affecting results strongly, shells with ages < 100 years are just bootstrapped, #age errors were overestimated too much with gamma distribution for young shells (not for old shells) for (x in 1:length(temp.bootages)) { temp=ifelse(temp.bootages[x]>100,rtrunc(1,spec="gamma", shape=temp.bootages[x]/FACTOR,scale=FACTOR,a=0, b=20000),temp.bootages[x]) temp.bootages[x]=temp } out1=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=250)) out3=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=10)) Full.Parvilucina.PVL1050.bootout.250[j,]=out1$counts Full.Parvilucina.PVL1050.bootout.10[j,]=out3$counts } ##################################################################### #absolute Parvilucina total abundance at OC50 ##################################################################### #2 6 10 14 18 22.5 27.5 32 35.5 39 43 47 51 55 59 63 67 71 75 79 83 totalsize=c(79,121,98,145,88,403,62,186,310,62,31,31,62,31,62,155,31,31,93,62,93) Full.Parvilucina.OC50.bootout.250=array(NA, dim=c(100,120)) Full.Parvilucina.OC50.bootout.10=array(NA, dim=c(100,3000)) for (j in 1:100) { temp.bootages=numeric() for (i in 1:13) { temp=unlist(Parvilucina.OC50.list1[i]) sample.age=sample(temp,totalsize[i],replace=T) temp.bootages=append(temp.bootages,sample.age,length(temp.bootages)) } #sample shell ages from the gamma distribution of age errors picked by the Bayesian model (shell ages represent medians of posterior distributions, not means), #using Parvilucina.ln.d from the Bayesian calibration #this step may be skipped - not affecting results strongly for (x in 1:length(temp.bootages)) { temp=ifelse(temp.bootages[x]>100,rtrunc(1,spec="gamma", shape=temp.bootages[x]/FACTOR,scale=FACTOR,a=0, b=30000),temp.bootages[x]) temp.bootages[x]=temp } out1=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=250)) out3=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=10)) Full.Parvilucina.OC50.bootout.250[j,]=out1$counts Full.Parvilucina.OC50.bootout.10[j,]=out3$counts } totalproduction.OC50=totalsize Parvilucina.PVL1050.10.means=apply(Full.Parvilucina.PVL1050.bootout.10,MARGIN=2,mean) Parvilucina.OC50.10.means=apply(Full.Parvilucina.OC50.bootout.10,MARGIN=2,mean) scale.10=seq(0,29990,by=10) Parvilucina.PVL1050.250.means=apply(Full.Parvilucina.PVL1050.bootout.250,MARGIN=2,mean) Parvilucina.PVL1050.250.025=apply(Full.Parvilucina.PVL1050.bootout.250,MARGIN=2,quantile, 0.025) Parvilucina.PVL1050.250.975=apply(Full.Parvilucina.PVL1050.bootout.250,MARGIN=2,quantile, 0.975) names(Parvilucina.PVL1050.250.means)=seq(0,29750,by=250) Parvilucina.OC50.250.means=apply(Full.Parvilucina.OC50.bootout.250,MARGIN=2,mean) Parvilucina.OC50.250.025=apply(Full.Parvilucina.OC50.bootout.250,MARGIN=2,quantile, 0.025) Parvilucina.OC50.250.975=apply(Full.Parvilucina.OC50.bootout.250,MARGIN=2,quantile, 0.975) names(Parvilucina.OC50.250.means)=seq(0,29750,by=250) scale.250=seq(0,29750,by=250) ###################################################################################################################### #CORRECTING FOR SHELL DISINTEGRATON IN PARVILUCINA ###################################################################################################################### quantile(boot.T.Parvilucina.lambda1.01, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.lambda1.02, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.lambda1.03, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.tau.01, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.tau.02, c(0.025,0.5,0.975)) quantile(boot.T.Parvilucina.tau.03, c(0.025,0.5,0.975)) #final values used in reconstructions - see Supplement Table 5 in Tomasovych et al. 2019 PP #the crucial value is lambda 1 and tau, lambda2 is set to zero (or very small value) #tau - sequestration rate, the estimates of tau are variable - between 0.003 and ~0.0001 #here tau is set to 1/mean time to burial below the TAZ (about 2000 years for both cores) #this should be mean time needed to get shells below the TAZ (20 cm) #used the same value for both Parvilucina and Nuculana tau<-0.0005 #lambda 2 is set to zero for both Parvilucina and Nuculana - this assumes zero disintegration below the TAZ, even when the models give nonzero values) #this parameter also incoporates loss of shells by burial from the last increment (i.e., beyond the last core increment) #e.g., if core has ~10,000 year-duration, this parameter would be on the order of 1/10000, so should be negligible lambda2=0 lambda1=0.042 #estimate for Parvilucina with truncation in production at 20 years lambda1.lci=0.035 lambda1.uci=0.052 #transform lambda1, lambda2 and tau to r1, r2, and alpha r1=lambda2 r2=lambda1+tau alpha=tau/(tau + lambda1 - lambda2) r2.lci=lambda1.lci+tau alpha.uci=tau/(tau + lambda1.uci - lambda2) r2.uci=lambda1.uci+tau alpha.lci=tau/(tau + lambda1.lci - lambda2) ###################################################################################################################### #THESE ARE SURVIVAL FUNCTIONS FOR EACH SITE AND DIFFERENT BIN SIZES (10, 250 YEARS) ###################################################################################################################### invsurv.PVL10.50.250=(alpha*exp(-r1*scale.250)+(1-alpha)*exp(-r2*scale.250)) invsurv.PVL10.50.250=(alpha*exp(-r1*scale.250)+(1-alpha)*exp(-r2*scale.250)) invsurv.OC.50.250=(alpha*exp(-r1*scale.250)+(1-alpha)*exp(-r2*scale.250)) invsurv.PVL10.50.10=(alpha*exp(-r1*scale.10)+(1-alpha)*exp(-r2*scale.10)) invsurv.PVL10.50.10.lci=(alpha.lci*exp(-r1*scale.10)+(1-alpha.lci)*exp(-r2.lci*scale.10)) invsurv.PVL10.50.10.uci=(alpha.uci*exp(-r1*scale.10)+(1-alpha.uci)*exp(-r2.uci*scale.10)) invsurv.OC.50.10=(alpha*exp(-r1*scale.10)+(1-alpha)*exp(-r2*scale.10)) invsurv.PVL10.50.250.lci=(alpha.lci*exp(-r1*scale.250)+(1-alpha.lci)*exp(-r2.lci*scale.250)) invsurv.OC.50.250.lci=(alpha.lci*exp(-r1*scale.250)+(1-alpha.lci)*exp(-r2.lci*scale.250)) invsurv.OC.50.10.lci=(alpha.lci*exp(-r1*scale.10)+(1-alpha.lci)*exp(-r2.lci*scale.10)) invsurv.PVL10.50.250.uci=(alpha.uci*exp(-r1*scale.250)+(1-alpha.uci)*exp(-r2.uci*scale.250)) invsurv.OC.50.250.uci=(alpha.uci*exp(-r1*scale.250)+(1-alpha.uci)*exp(-r2.uci*scale.250)) invsurv.OC.50.10.uci=(alpha.uci*exp(-r1*scale.10)+(1-alpha.uci)*exp(-r2.uci*scale.10)) ###################################################################################################################### #DIVIDE BY SURVIVAL FUNCTIONS - RECONSTRUCT ABUNDANCES ###################################################################################################################### Corr.Full.Parvilucina.PVL1050.bootout.10=t(Full.Parvilucina.PVL1050.bootout.10)/invsurv.PVL10.50.10 Corr.Full.Parvilucina.OC50.bootout.10=t(Full.Parvilucina.OC50.bootout.10)/invsurv.OC.50.10 Corr.Full.Parvilucina.PVL1050.bootout.250=t(Full.Parvilucina.PVL1050.bootout.250)/invsurv.PVL10.50.250 Corr.Full.Parvilucina.OC50.bootout.250=t(Full.Parvilucina.OC50.bootout.250)/invsurv.OC.50.250 Corr.Parvilucina.PVL1050.10.means=apply(Corr.Full.Parvilucina.PVL1050.bootout.10,MARGIN=1,mean) Corr.Parvilucina.OC50.10.means=apply(Corr.Full.Parvilucina.OC50.bootout.10,MARGIN=1,mean) Corr.Full.Parvilucina.PVL1050.bootout.10.lci=t(Full.Parvilucina.PVL1050.bootout.10)/invsurv.PVL10.50.10.lci Corr.Full.Parvilucina.OC50.bootout.10.lci=t(Full.Parvilucina.OC50.bootout.10)/invsurv.OC.50.10.lci Corr.Full.Parvilucina.PVL1050.bootout.10.uci=t(Full.Parvilucina.PVL1050.bootout.10)/invsurv.PVL10.50.10.uci Corr.Full.Parvilucina.OC50.bootout.10.uci=t(Full.Parvilucina.OC50.bootout.10)/invsurv.OC.50.10.uci Corr.Parvilucina.PVL1050.10.means.lower=apply(Corr.Full.Parvilucina.PVL1050.bootout.10.lci,MARGIN=1,mean) Corr.Parvilucina.PVL1050.10.means.upper=apply(Corr.Full.Parvilucina.PVL1050.bootout.10.uci,MARGIN=1,mean) Corr.Parvilucina.OC50.10.means.lower=apply(Corr.Full.Parvilucina.OC50.bootout.10.lci,MARGIN=1,mean) Corr.Parvilucina.OC50.10.means.upper=apply(Corr.Full.Parvilucina.OC50.bootout.10.uci,MARGIN=1,mean) Corr.Parvilucina.PVL1050.250.means=apply(Corr.Full.Parvilucina.PVL1050.bootout.250,MARGIN=1,mean) Corr.Parvilucina.OC50.250.means=apply(Corr.Full.Parvilucina.OC50.bootout.250,MARGIN=1,mean) ####################################################################################################### #RECONSTRUCTION OF NUCULANA TRAJECTORY (DISINTEGRATION NOT ACCOUNTED FOR) ####################################################################################################### ##################################################################### #absolute abundance of Nuculana taphria in increments at PVL10-50 ##################################################################### totalsize=c(93,137,184,252,278,228,1008,1386,1029,1134,966,1386,1008) Full.Nuculana.PVL1050.bootout.250=array(NA, dim=c(100,120)) Full.Nuculana.PVL1050.bootout.10=array(NA, dim=c(100,3000)) for (j in 1:100) { temp.bootages=numeric() for (i in 1:length(Nuculana.PVL1050.list1)) { temp=unlist(Nuculana.PVL1050.list1[i]) sample.age=sample(temp,totalsize[i],replace=T) temp.bootages=append(temp.bootages,sample.age,length(temp.bootages)) } #accounting for age error calibration (lognormal errors for Nuculana) for (x in 1:length(temp.bootages)) { temp.bootages[x]=rtlnorm(1,log(temp.bootages[x]),sdlog=sqrt(exp(Nuculana.ln.d)),a=0,b=20000) } out1=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=250)) out3=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=10)) Full.Nuculana.PVL1050.bootout.250[j,]=out1$counts Full.Nuculana.PVL1050.bootout.10[j,]=out3$counts } ##################################################################### #absolute abundance of Nuculana taphria in increments at OC-50 ##################################################################### # 2 6 10 14 18 22.5 27.5 32 35.5 39 43 47 51 55 59 63 67 71 75 79 83 totalsize=c(30,77,79,100,89,620,1426,992,837,992,465,651,496,899,1767,961,930,1147,1178,1395,1240,341,403,341,31,341) Full.Nuculana.OC50.bootout.250=array(NA, dim=c(100,120)) Full.Nuculana.OC50.bootout.10=array(NA, dim=c(100,3000)) for (j in 1:100) { temp.bootages=numeric() for (i in 1:length(totalsize)) { temp=unlist(Nuculana.OC50.list1[i]) sample.age=sample(temp,totalsize[i],replace=T) temp.bootages=append(temp.bootages,sample.age,length(temp.bootages)) } #accounting for age error calibration (lognormal errors for Nuculana) for (x in 1:length(temp.bootages)) { temp.bootages[x]=rtlnorm(1,log(temp.bootages[x]),sdlog=sqrt(exp(Nuculana.ln.d)),a=0,b=20000) } out1=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=250)) out3=hist(temp.bootages[temp.bootages<30000], breaks=seq(0,30000,by=10)) Full.Nuculana.OC50.bootout.250[j,]=out1$counts Full.Nuculana.OC50.bootout.10[j,]=out3$counts } Nuculana.PVL1050.10.means=apply(Full.Nuculana.PVL1050.bootout.10,MARGIN=2,mean) Nuculana.OC50.10.means=apply(Full.Nuculana.OC50.bootout.10,MARGIN=2,mean) scale.10=seq(0,29990,by=10) Nuculana.PVL1050.250.means=apply(Full.Nuculana.PVL1050.bootout.250,MARGIN=2,mean) Nuculana.PVL1050.250.025=apply(Full.Nuculana.PVL1050.bootout.250,MARGIN=2,quantile, 0.025) Nuculana.PVL1050.250.975=apply(Full.Nuculana.PVL1050.bootout.250,MARGIN=2,quantile, 0.975) names(Nuculana.PVL1050.250.means)=seq(0,29750,by=250) scale.250=seq(0,29750,by=250) Nuculana.OC50.250.means=apply(Full.Nuculana.OC50.bootout.250,MARGIN=2,mean) Nuculana.OC50.250.025=apply(Full.Nuculana.OC50.bootout.250,MARGIN=2,quantile, 0.025) Nuculana.OC50.250.975=apply(Full.Nuculana.OC50.bootout.250,MARGIN=2,quantile, 0.975) names(Nuculana.OC50.250.means)=seq(0,29750,by=250) scale.250=seq(0,29750,by=250) ###################################################################################################################### #CORRECTING FOR SHELL DISINTEGRATION IN NUCULANA ###################################################################################################################### T.Nuculana.lambda1.lci=quantile(boot.T.Nuculana.lambda1, 0.025) T.Nuculana.lambda1.uci=quantile(boot.T.Nuculana.lambda1, 0.975) T.Nuculana.lambda1.median=quantile(boot.T.Nuculana.lambda1, 0.5) T.Nuculana.tau.lci=quantile(boot.T.Nuculana.tau, 0.025) T.Nuculana.tau.median=quantile(boot.T.Nuculana.tau, 0.5) T.Nuculana.tau.uci=quantile(boot.T.Nuculana.tau, 0.975) T.Nuculana.r2.lci=quantile(boot.T.Nuculana.r2, 0.025) T.Nuculana.r2.median=quantile(boot.T.Nuculana.r2, 0.5) T.Nuculana.r2.uci=quantile(boot.T.Nuculana.r2, 0.975) #final values used in reconstructions - see Supplement Table 5 in Tomasovych et al. 2019 PP #the crucial value is lambda 1 and tau, lambda2 is set to zero (or very small value) #tau - sequestration rate, the estimates of tau are variable - between 0.003 and ~0.0001 #here tau is set to 1/mean time to burial below the TAZ (about 2000 years for both cores) #this should be mean time needed to get shells below the TAZ (20 cm) #used the same value for both Parvilucina and Nuculana tau<-0.0005 #lambda 2 is set to zero for both Parvilucina and Nuculana - this assumes zero disintegration below the TAZ, even when the models give nonzero values) #this parameter also incoporates loss of shells by burial from the last increment (i.e., beyond the last core increment) #e.g., if core has ~10,000 year-duration, this parameter would be on the order of 1/10000, so should be negligible lambda2=0 lambda1=T.Nuculana.r2.median #0.047 from Supplementary Table 5 lambda1.uci=T.Nuculana.r2.uci lambda1.lci=T.Nuculana.r2.lci r2=lambda1+tau alpha=tau/(tau + lambda1 - lambda2) r2.lci=lambda1.lci+tau alpha.uci=tau/(tau + lambda1.uci - lambda2) r2.uci=lambda1.uci+tau alpha.lci=tau/(tau + lambda1.lci - lambda2) ###################################################################################################################### #THESE ARE SURVIVAL FUNCTIONS FOR EACH SITE AND DIFFERENT BIN SIZES (10, 250 YEARS) ###################################################################################################################### invsurv.PVL10.50.250=(alpha*exp(-r1*scale.250)+(1-alpha)*exp(-r2*scale.250)) invsurv.OC.50.250=(alpha*exp(-r1*scale.250)+(1-alpha)*exp(-r2*scale.250)) invsurv.PVL10.50.10=(alpha*exp(-r1*scale.10)+(1-alpha)*exp(-r2*scale.10)) invsurv.OC.50.10=(alpha*exp(-r1*scale.10)+(1-alpha)*exp(-r2*scale.10)) invsurv.PVL10.50.250.lci=(alpha.lci*exp(-r1*scale.250)+(1-alpha.lci)*exp(-r2.lci*scale.250)) invsurv.OC.50.250.lci=(alpha.lci*exp(-r1*scale.250)+(1-alpha.lci)*exp(-r2.lci*scale.250)) invsurv.PVL10.50.10.lci=(alpha.lci*exp(-r1*scale.10)+(1-alpha.lci)*exp(-r2.lci*scale.10)) invsurv.OC.50.10.lci=(alpha.lci*exp(-r1*scale.10)+(1-alpha.lci)*exp(-r2.lci*scale.10)) invsurv.PVL10.50.250.uci=(alpha.uci*exp(-r1*scale.250)+(1-alpha.uci)*exp(-r2.uci*scale.250)) invsurv.OC.50.250.uci=(alpha.uci*exp(-r1*scale.250)+(1-alpha.uci)*exp(-r2.uci*scale.250)) invsurv.PVL10.50.10.uci=(alpha.uci*exp(-r1*scale.10)+(1-alpha.uci)*exp(-r2.uci*scale.10)) invsurv.OC.50.10.uci=(alpha.uci*exp(-r1*scale.10)+(1-alpha.uci)*exp(-r2.uci*scale.10)) ###################################################################################################################### #DIVIDE BY SURVIVAL FUNCTIONS - RECONSTRUCT ABUNDANCES ###################################################################################################################### Corr.Full.Nuculana.PVL1050.bootout.250=t(Full.Nuculana.PVL1050.bootout.250)/invsurv.PVL10.50.250 Corr.Full.Nuculana.OC50.bootout.250=t(Full.Nuculana.OC50.bootout.250)/invsurv.OC.50.250 Corr.Full.Nuculana.PVL1050.bootout.10=t(Full.Nuculana.PVL1050.bootout.10)/invsurv.PVL10.50.10 Corr.Full.Nuculana.OC50.bootout.10=t(Full.Nuculana.OC50.bootout.10)/invsurv.OC.50.10 Corr.Nuculana.PVL1050.10.means=apply(Corr.Full.Nuculana.PVL1050.bootout.10,MARGIN=1,mean) Corr.Nuculana.OC50.10.means=apply(Corr.Full.Nuculana.OC50.bootout.10,MARGIN=1,mean) Corr.Full.Nuculana.PVL1050.bootout.10.lci=t(Full.Nuculana.PVL1050.bootout.10)/invsurv.PVL10.50.10.lci Corr.Full.Nuculana.OC50.bootout.10.lci=t(Full.Nuculana.OC50.bootout.10)/invsurv.OC.50.10.lci Corr.Full.Nuculana.PVL1050.bootout.10.uci=t(Full.Nuculana.PVL1050.bootout.10)/invsurv.PVL10.50.10.uci Corr.Full.Nuculana.OC50.bootout.10.uci=t(Full.Nuculana.OC50.bootout.10)/invsurv.OC.50.10.uci Corr.Nuculana.PVL1050.10.means.lower=apply(Corr.Full.Nuculana.PVL1050.bootout.10.lci,MARGIN=1,mean) Corr.Nuculana.PVL1050.10.means.upper=apply(Corr.Full.Nuculana.PVL1050.bootout.10.uci,MARGIN=1,mean) Corr.Nuculana.OC50.10.means.lower=apply(Corr.Full.Nuculana.OC50.bootout.10.lci,MARGIN=1,mean) Corr.Nuculana.OC50.10.means.upper=apply(Corr.Full.Nuculana.OC50.bootout.10.uci,MARGIN=1,mean) Corr.Nuculana.PVL1050.250.means=apply(Corr.Full.Nuculana.PVL1050.bootout.250,MARGIN=1,mean) Corr.Nuculana.OC50.250.means=apply(Corr.Full.Nuculana.OC50.bootout.250,MARGIN=1,mean) ###################################################################################################################### #FIGURE PARVILUCINA RECONSTRUCTED TRAJECTORY ###################################################################################################################### par(mfcol=c(2,3)) par(mar=c(4,4,2,1)) horiz.hist(-Parvilucina.PVL1050.ages[Parvilucina.PVL1050.ages<15000], breaks=seq(-15000,0,by=250), xlim=c(0,100), col="gray", main="", cex.main=0.9, xlab="Abundance") horiz.hist(-Parvilucina.OC50.ages[Parvilucina.OC50.ages<15000], breaks=seq(-15000,0,by=250), xlim=c(0,100), col="gray", main="", cex.main=0.9, xlab="Abundance") barplot(Parvilucina.PVL1050.250.means[60:1], xlim=c(0,800), las=2, xaxt="n", horiz=TRUE, ylab="", names.arg="", xlab="Extrapolated abundance") cc=barplot(Parvilucina.PVL1050.250.means[1:60],horiz=TRUE, plot=F) axis(2,at=cc[c(1,20,40,60),1],labels=-rev(c(1,20,40,60))*250, las=2) axis(1,at=c(0,400,800),labels=c(0,400,800)) barplot(Parvilucina.OC50.250.means[60:1], xlim=c(0,800), las=2, xaxt="n", horiz=TRUE, ylab="", names.arg="", xlab="Extrapolated abundance") cc=barplot(Parvilucina.OC50.250.means[1:60],horiz=TRUE, plot=F) axis(2,at=cc[c(1,20,40,60),1],labels=-rev(c(1,20,40,60))*250, las=2) axis(1,at=c(0,400,800),labels=c(0,400,800)) plot(Corr.Parvilucina.PVL1050.10.means[1:1500], -c(1:1500)*10, type="l", ylab="Time (y)", frame.plot=F, xlab="Reconstructed abundance/10 y", xlim=c(0,4000), ylim=c(-15000,0)) lines(Corr.Parvilucina.PVL1050.10.means.lower[1:1500], -c(1:1500)*10,lty=1, col="gray", lwd=2) lines(Corr.Parvilucina.PVL1050.10.means.upper[1:1500],-c(1:1500)*10,lty=1, col="gray", lwd=2) plot(Corr.Parvilucina.OC50.10.means[1:1500], -c(1:1500)*10,type="l", ylab="Time (y)", frame.plot=F, xlab="Reconstructed abundance/10 y", xlim=c(0,4000), ylim=c(-15000,0)) lines(Corr.Parvilucina.OC50.10.means.lower[1:1500], -c(1:1500)*10,lty=1, col="gray", lwd=2) lines(Corr.Parvilucina.OC50.10.means.upper[1:1500], -c(1:1500)*10,lty=1, col="gray", lwd=2) ###################################################################################################################### #FIGURE NUCULANA RECONSTRUCTED TRAJECTORY ###################################################################################################################### par(mfcol=c(2,3)) par(mar=c(4,4,2,1)) horiz.hist(-Nuculana.PVL1050.ages[Nuculana.PVL1050.ages<15000], breaks=seq(-15000,0,by=250), xlim=c(0,40), col="gray", main="", cex.main=0.9, xlab="Abundance") horiz.hist(-Nuculana.OC50.ages[Nuculana.OC50.ages<15000], breaks=seq(-15000,0,by=250), xlim=c(0,40), col="gray", main="", cex.main=0.9, xlab="Abundance") barplot(Nuculana.PVL1050.250.means[60:1], xlim=c(0,300), las=2, xaxt="n", horiz=TRUE, ylab="", names.arg="", xlab="Extrapolated abundance") cc=barplot(Nuculana.PVL1050.250.means[1:60],horiz=TRUE, plot=F) axis(2,at=cc[c(1,20,40,60),1],labels=-rev(c(1,20,40,60))*250, las=2) axis(1,at=c(0,100,200,300),labels=c(0,100,200,300)) barplot(Nuculana.OC50.250.means[60:1], xlim=c(0,300), las=2, xaxt="n", horiz=TRUE, ylab="", names.arg="", xlab="Extrapolated abundance") cc=barplot(Nuculana.OC50.250.means[1:60],horiz=TRUE, plot=F) axis(2,at=cc[c(1,20,40,60),1],labels=-rev(c(1,20,40,60))*250, las=2) axis(1,at=c(0,100,200,300),labels=c(0,100,200,300)) plot(Corr.Nuculana.PVL1050.10.means[1:1500], -c(1:1500)*10, type="l", ylab="", xlab="Reconstructed abundance/10 y", xlim=c(0,4000), ylim=c(-15000,0)) lines(Corr.Nuculana.PVL1050.10.means.lower[1:1500],-c(1:1500)*10, lty=1, col="gray", lwd=2) lines(Corr.Nuculana.PVL1050.10.means.upper[1:1500],-c(1:1500)*10, lty=1, col="gray", lwd=2) plot(Corr.Nuculana.OC50.10.means[1:1500],-c(1:1500)*10, type="l",ylab="", xlab="Reconstructed abundance/10 y", xlim=c(0,4000), ylim=c(-15000,0)) lines(Corr.Nuculana.OC50.10.means.lower[1:1500],-c(1:1500)*10, lty=1, col="gray", lwd=2) lines(Corr.Nuculana.OC50.10.means.upper[1:1500], -c(1:1500)*10, lty=1, col="gray", lwd=2) ###################################################################################################################### #PARVILUCINA AND NUCULANA RECONSTRUCTED TRAJECTORY - VERTICAL TIME ###################################################################################################################### par(mfrow=c(1,2)) par(mar=c(15,4,2,1)) plot(Corr.Nuculana.PVL1050.10.means[1:100], -c(1:100)*10, type="l", lwd=2, ylab="Time (since 2012)", xlim=c(0,2500), frame.plot=F, xlab="Reconstructed abundance/10 y", main="PVL10-50", cex.main=0.9) lines(Corr.Nuculana.PVL1050.10.means.upper[1:100],-c(1:100)*10, lty=2, lwd=1) lines(Corr.Nuculana.PVL1050.10.means.lower[1:100],-c(1:100)*10, lty=2, lwd=1) lines(Corr.Parvilucina.PVL1050.10.means[1:100], -c(1:100)*10, col="gray", lwd=2) lines(Corr.Parvilucina.PVL1050.10.means.lower[1:100],-c(1:100)*10, lty=1, col="gray", lwd=1) lines(Corr.Parvilucina.PVL1050.10.means.upper[1:100],-c(1:100)*10, lty=1, col="gray", lwd=1) plot(Corr.Nuculana.OC50.10.means[1:100],-c(1:100)*10, type="l", lwd=2, ylab="", xlim=c(0,2500), frame.plot=F, xlab="Reconstructed abundance/10 y", main="OC-50", cex.main=0.9) lines(Corr.Nuculana.OC50.10.means.upper[1:100],-c(1:100)*10, lty=2, lwd=1) lines(Corr.Nuculana.OC50.10.means.lower[1:100], -c(1:100)*10, lty=2, lwd=1) lines(Corr.Parvilucina.OC50.10.means[1:100],-c(1:100)*10, col="gray", lwd=2) lines(Corr.Parvilucina.OC50.10.means.lower[1:100], -c(1:100)*10, lty=2, col="gray", lwd=1) lines(Corr.Parvilucina.OC50.10.means.upper[1:100], -c(1:100)*10, lty=2, col="gray", lwd=1) legend(x="bottomright", lwd=c(2,2), lty=c(1,1),legend=c("N. taphria","Parvilucina"), cex=1, bty="n", col=c("black","gray")) ###################################################################################################################### #FIGURE 11 PARVILUCINA AND NUCULANA RECONSTRUCTED TRAJECTORY - HORIZONTAL TIME ###################################################################################################################### par(mfrow=c(2,2)) par(mar=c(6,2,2,1)) plot(-c(1:100)*10, Corr.Nuculana.PVL1050.10.means[1:100], type="l", lwd=2, xlab="Calendar year (AD)", ylim=c(0,2000), frame.plot=F, ylab="Reconstructed abundance/10 y", main="PVL10-50", cex.main=0.9, xaxt="n", cex.axis=0.9) axis(1, at=seq(-820,-20,by=200), labels=seq(1200,2000,by=200), cex.axis=0.9) lines(-c(1:100)*10, Corr.Nuculana.PVL1050.10.means.upper[1:100], lty=2, lwd=1) lines(-c(1:100)*10, Corr.Nuculana.PVL1050.10.means.lower[1:100], lty=2, lwd=1) lines(-c(1:100)*10, Corr.Parvilucina.PVL1050.10.means[1:100], col="gray", lwd=2) lines(-c(1:100)*10, Corr.Parvilucina.PVL1050.10.means.lower[1:100], lty=1, col="gray", lwd=1) lines(-c(1:100)*10, Corr.Parvilucina.PVL1050.10.means.upper[1:100], lty=1, col="gray", lwd=1) plot(-c(1:100)*10, Corr.Nuculana.OC50.10.means[1:100], type="l", lwd=2, xlab="Calendar year (AD)", ylim=c(0,2000), frame.plot=F, ylab="", main="OC-50", cex.main=0.9, xaxt="n", cex.axis=0.9) axis(1, at=seq(-820,-20,by=200), labels=seq(1200,2000,by=200), cex.axis=0.9) lines(-c(1:100)*10, Corr.Nuculana.OC50.10.means.upper[1:100], lty=2, lwd=1) lines(-c(1:100)*10, Corr.Nuculana.OC50.10.means.lower[1:100], lty=2, lwd=1) lines(-c(1:100)*10, Corr.Parvilucina.OC50.10.means[1:100], col="gray", lwd=2) lines(-c(1:100)*10, Corr.Parvilucina.OC50.10.means.lower[1:100], lty=2, col="gray", lwd=1) lines(-c(1:100)*10, Corr.Parvilucina.OC50.10.means.upper[1:100], lty=2, col="gray", lwd=1) legend(x="topright", lwd=c(2,2), lty=c(1,1),legend=c("N. taphria","Parvilucina"), cex=1, bty="n", col=c("black","gray"))