## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # library(DEoptimR) # N <- 1500 # # # generate 5 sets of random calendar dates under the toy model. # set.seed(882) # cal1 <- simulateCalendarDates(model = toy, N) # set.seed(884) # cal2 <- simulateCalendarDates(model = toy, N) # set.seed(886) # cal3 <- simulateCalendarDates(model = toy, N) # set.seed(888) # cal4 <- simulateCalendarDates(model = toy, N) # set.seed(890) # cal5 <- simulateCalendarDates(model = toy, N) # # # Convert to 14C dates. # age1 <- uncalibrateCalendarDates(cal1, shcal20) # age2 <- uncalibrateCalendarDates(cal2, shcal20) # age3 <- uncalibrateCalendarDates(cal3, shcal20) # age4 <- uncalibrateCalendarDates(cal4, shcal20) # age5 <- uncalibrateCalendarDates(cal5, shcal20) # # # construct data frames. One date per phase. # data1 <- data.frame(age = age1, sd = 25, phase = 1:N, datingType = '14C') # data2 <- data.frame(age = age2, sd = 25, phase = 1:N, datingType = '14C') # data3 <- data.frame(age = age3, sd = 25, phase = 1:N, datingType = '14C') # data4 <- data.frame(age = age4, sd = 25, phase = 1:N, datingType = '14C') # data5 <- data.frame(age = age5, sd = 25, phase = 1:N, datingType = '14C') # # # Calibrate each phase, taking care to restrict to the modelled date range # CalArray <- makeCalArray(shcal20, calrange = range(toy$year), inc = 5) # # PD1 <- phaseCalibrator(data1, CalArray, remove.external = TRUE) # PD2 <- phaseCalibrator(data2, CalArray, remove.external = TRUE) # PD3 <- phaseCalibrator(data3, CalArray, remove.external = TRUE) # PD4 <- phaseCalibrator(data4, CalArray, remove.external = TRUE) # PD5 <- phaseCalibrator(data5, CalArray, remove.external = TRUE) # # # Generate SPD of each dataset # SPD1 <- summedCalibrator(data1, CalArray, normalise='full') # SPD2 <- summedCalibrator(data2, CalArray, normalise='full') # SPD3 <- summedCalibrator(data3, CalArray, normalise='full') # SPD4 <- summedCalibrator(data4, CalArray, normalise='full') # SPD5 <- summedCalibrator(data5, CalArray, normalise='full') # # # 3-CPL parameter search # lower <- rep(0,5) # upper <- rep(1,5) # fn <- objectiveFunction # best1 <- JDEoptim(lower, upper, fn, PDarray=PD1, type='CPL',trace=T,NP=100) # best2 <- JDEoptim(lower, upper, fn, PDarray=PD2, type='CPL',trace=T,NP=100) # best3 <- JDEoptim(lower, upper, fn, PDarray=PD3, type='CPL',trace=T,NP=100) # best4 <- JDEoptim(lower, upper, fn, PDarray=PD4, type='CPL',trace=T,NP=100) # best5 <- JDEoptim(lower, upper, fn, PDarray=PD5, type='CPL',trace=T,NP=100) # # #save results, for separate plotting # save(best1,best2,best3,best4,best5,SPD1,SPD2,SPD3,SPD4,SPD5, file='results.RData',version=2) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # load('results.RData') # oldpar <- par(no.readonly = TRUE) # pdf('Fig1.pdf',height=4,width=10) # par(mar=c(2,4,0.1,2)) # plot(NULL, xlim=c(7500,5500), ylim=c(0,0.0011), xlab='', ylab='', xaxs='i',cex.axis=0.7, bty='n',las=1) # axis(1,at=6400,labels='calBP',tick=F) # axis(2,at=-0.00005,labels='PD',tick=F, las=1) # lwd1 <- 1 # lwd2 <- 2 # lwd3 <- 3 # legend(x=6000, y = 0.0011, bty='n', cex=0.7, # legend=c('True (toy) population', # 'SPD 1', # 'SPD 2', # 'SPD 3', # 'SPD 4', # 'SPD 5', # 'Pop model 1', # 'Pop model 2', # 'Pop model 3', # 'Pop model 4', # 'Pop model 5'), # lwd=c(lwd3,rep(lwd1,5),rep(lwd2,5)), # col=c(1,2:6,2:6) # ) # # years <- as.numeric(row.names(SPD1)) # # # plot SPDs # lines(years,SPD1[,1],col=2, lwd=lwd1) # lines(years,SPD2[,1],col=3, lwd=lwd1) # lines(years,SPD3[,1],col=4, lwd=lwd1) # lines(years,SPD4[,1],col=5, lwd=lwd1) # lines(years,SPD5[,1],col=6, lwd=lwd1) # # # convert parameters to model pdfs # mod.1 <- convertPars(pars=best1$par, years=years, type='CPL') # mod.2 <- convertPars(pars=best2$par, years=years, type='CPL') # mod.3 <- convertPars(pars=best3$par, years=years, type='CPL') # mod.4 <- convertPars(pars=best4$par, years=years, type='CPL') # mod.5 <- convertPars(pars=best5$par, years=years, type='CPL') # # lines(mod.1$year,mod.1$pdf,col=2,lwd=lwd2) # lines(mod.2$year,mod.2$pdf,col=3,lwd=lwd2) # lines(mod.3$year,mod.3$pdf,col=4,lwd=lwd2) # lines(mod.4$year,mod.4$pdf,col=5,lwd=lwd2) # lines(mod.5$year,mod.5$pdf,col=6,lwd=lwd2) # # # plot true toy model # lines(toy$year, toy$pdf, lwd=lwd3) # # dev.off() # par(oldpar) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # library(DEoptimR) # set.seed(888) # N <- c(6,20,60,180,360,540) # names <- c('sample1','sample2','sample3','sample4','sample5','sample6') # # # generate 6 sets of random calendar dates under the toy model. # cal1 <- simulateCalendarDates(model = toy, N[1]) # cal2 <- simulateCalendarDates(model = toy, N[2]) # cal3 <- simulateCalendarDates(model = toy, N[3]) # cal4 <- simulateCalendarDates(model = toy, N[4]) # cal5 <- simulateCalendarDates(model = toy, N[5]) # cal6 <- simulateCalendarDates(model = toy, N[6]) # # # Convert to 14C dates. # age1 <- uncalibrateCalendarDates(cal1, shcal20) # age2 <- uncalibrateCalendarDates(cal2, shcal20) # age3 <- uncalibrateCalendarDates(cal3, shcal20) # age4 <- uncalibrateCalendarDates(cal4, shcal20) # age5 <- uncalibrateCalendarDates(cal5, shcal20) # age6 <- uncalibrateCalendarDates(cal6, shcal20) # # # construct data frames. One date per phase. # data1 <- data.frame(age = age1, sd = 25, phase = 1:N[1], datingType = '14C') # data2 <- data.frame(age = age2, sd = 25, phase = 1:N[2], datingType = '14C') # data3 <- data.frame(age = age3, sd = 25, phase = 1:N[3], datingType = '14C') # data4 <- data.frame(age = age4, sd = 25, phase = 1:N[4], datingType = '14C') # data5 <- data.frame(age = age5, sd = 25, phase = 1:N[5], datingType = '14C') # data6 <- data.frame(age = age6, sd = 25, phase = 1:N[6], datingType = '14C') # # # narrow domain of the model to the range of data, # # since absence of evidence in periods well outside the data should # # not be interpreted as evidence of absence. # # Only required when sample sizes are extremely small. # # Otherwise the data domain is constrained by the model date range. # r1 <- estimateDataDomain(data1, shcal20) # # # narrower range for extremely small samples # CalArray1 <- makeCalArray(shcal20, calrange = c( max(r1[1],5500) , min(r1[2],7500) ), inc = 5) # CalArray <- makeCalArray(shcal20, calrange = range(toy$year), inc = 5) # # # Calibrate each phase # PD1 <- phaseCalibrator(data1, CalArray1, remove.external = TRUE) # PD2 <- phaseCalibrator(data2, CalArray, remove.external = TRUE) # PD3 <- phaseCalibrator(data3, CalArray, remove.external = TRUE) # PD4 <- phaseCalibrator(data4, CalArray, remove.external = TRUE) # PD5 <- phaseCalibrator(data5, CalArray, remove.external = TRUE) # PD6 <- phaseCalibrator(data6, CalArray, remove.external = TRUE) # PD <- list(PD1, PD2, PD3, PD4, PD5, PD6); names(PD) <- names # # # Generate SPD of each dataset # SPD1 <- summedCalibrator(data1, CalArray, normalise='full') # SPD2 <- summedCalibrator(data2, CalArray, normalise='full') # SPD3 <- summedCalibrator(data3, CalArray, normalise='full') # SPD4 <- summedCalibrator(data4, CalArray, normalise='full') # SPD5 <- summedCalibrator(data5, CalArray, normalise='full') # SPD6 <- summedCalibrator(data6, CalArray, normalise='full') # SPD <- list(SPD1, SPD2, SPD3, SPD4, SPD5, SPD6); names(SPD) <- names # # # Uniform model: No parameters. # # Log Likelihood calculated directly using objectiveFunction, without a search required. # unif1.loglik <- -objectiveFunction(pars = NULL, PDarray = PD1, type = 'uniform') # unif2.loglik <- -objectiveFunction(pars = NULL, PDarray = PD2, type = 'uniform') # unif3.loglik <- -objectiveFunction(pars = NULL, PDarray = PD3, type = 'uniform') # unif4.loglik <- -objectiveFunction(pars = NULL, PDarray = PD4, type = 'uniform') # unif5.loglik <- -objectiveFunction(pars = NULL, PDarray = PD5, type = 'uniform') # unif6.loglik <- -objectiveFunction(pars = NULL, PDarray = PD6, type = 'uniform') # uniform <- list(unif1.loglik, unif2.loglik, unif3.loglik, unif4.loglik, unif5.loglik, unif6.loglik) # names(uniform) <- names # # # Best 1-CPL model. Parameters and log likelihood found using search # lower <- rep(0,1) # upper <- rep(1,1) # fn <- objectiveFunction # best1 <- JDEoptim(lower, upper, fn, PDarray=PD1, type='CPL',trace=T,NP=20) # best2 <- JDEoptim(lower, upper, fn, PDarray=PD2, type='CPL',trace=T,NP=20) # best3 <- JDEoptim(lower, upper, fn, PDarray=PD3, type='CPL',trace=T,NP=20) # best4 <- JDEoptim(lower, upper, fn, PDarray=PD4, type='CPL',trace=T,NP=20) # best5 <- JDEoptim(lower, upper, fn, PDarray=PD5, type='CPL',trace=T,NP=20) # best6 <- JDEoptim(lower, upper, fn, PDarray=PD6, type='CPL',trace=T,NP=20) # CPL1 <- list(best1, best2, best3, best4, best5, best6); names(CPL1) <- names # # # Best 2-CPL model. Parameters and log likelihood found using search # lower <- rep(0,3) # upper <- rep(1,3) # fn <- objectiveFunction # best1 <- JDEoptim(lower, upper, fn, PDarray=PD1, type='CPL',trace=T,NP=60) # best2 <- JDEoptim(lower, upper, fn, PDarray=PD2, type='CPL',trace=T,NP=60) # best3 <- JDEoptim(lower, upper, fn, PDarray=PD3, type='CPL',trace=T,NP=60) # best4 <- JDEoptim(lower, upper, fn, PDarray=PD4, type='CPL',trace=T,NP=60) # best5 <- JDEoptim(lower, upper, fn, PDarray=PD5, type='CPL',trace=T,NP=60) # best6 <- JDEoptim(lower, upper, fn, PDarray=PD6, type='CPL',trace=T,NP=60) # CPL2 <- list(best1, best2, best3, best4, best5, best6); names(CPL2) <- names # # # Best 3-CPL model. Parameters and log likelihood found using search # lower <- rep(0,5) # upper <- rep(1,5) # fn <- objectiveFunction # best1 <- JDEoptim(lower, upper, fn, PDarray=PD1, PDarray=PD1, type='CPL',trace=T,NP=100) # best2 <- JDEoptim(lower, upper, fn, PDarray=PD1, PDarray=PD2, type='CPL',trace=T,NP=100) # best3 <- JDEoptim(lower, upper, fn, PDarray=PD1, PDarray=PD3, type='CPL',trace=T,NP=100) # best4 <- JDEoptim(lower, upper, fn, PDarray=PD1, PDarray=PD4, type='CPL',trace=T,NP=100) # best5 <- JDEoptim(lower, upper, fn, PDarray=PD1, PDarray=PD5, type='CPL',trace=T,NP=100) # best6 <- JDEoptim(lower, upper, fn, PDarray=PD1, PDarray=PD6, type='CPL',trace=T,NP=100) # CPL3 <- list(best1, best2, best3, best4, best5, best6); names(CPL3) <- names # # # Best 4-CPL model. Parameters and log likelihood found using search # lower <- rep(0,7) # upper <- rep(1,7) # fn <- objectiveFunction # best1 <- JDEoptim(lower, upper, fn, PDarray=PD1, type='CPL',trace=T,NP=140) # best2 <- JDEoptim(lower, upper, fn, PDarray=PD2, type='CPL',trace=T,NP=140) # best3 <- JDEoptim(lower, upper, fn, PDarray=PD3, type='CPL',trace=T,NP=140) # best4 <- JDEoptim(lower, upper, fn, PDarray=PD4, type='CPL',trace=T,NP=140) # best5 <- JDEoptim(lower, upper, fn, PDarray=PD5, type='CPL',trace=T,NP=140) # best6 <- JDEoptim(lower, upper, fn, PDarray=PD6, type='CPL',trace=T,NP=140) # CPL4 <- list(best1, best2, best3, best4, best5, best6); names(CPL4) <- names # # # Best 5-CPL model. Parameters and log likelihood found using search # lower <- rep(0,9) # upper <- rep(1,9) # fn <- objectiveFunction # best1 <- JDEoptim(lower, upper, fn, PDarray=PD1, type='CPL',trace=T,NP=180) # best2 <- JDEoptim(lower, upper, fn, PDarray=PD2, type='CPL',trace=T,NP=180) # best3 <- JDEoptim(lower, upper, fn, PDarray=PD3, type='CPL',trace=T,NP=180) # best4 <- JDEoptim(lower, upper, fn, PDarray=PD4, type='CPL',trace=T,NP=180) # best5 <- JDEoptim(lower, upper, fn, PDarray=PD5, type='CPL',trace=T,NP=180) # best6 <- JDEoptim(lower, upper, fn, PDarray=PD6, type='CPL',trace=T,NP=180) # CPL5 <- list(best1, best2, best3, best4, best5, best6); names(CPL5) <- names # # # save results, for separate plotting # save(SPD, PD, uniform, CPL1, CPL2, CPL3, CPL4, CPL5, file='results.RData',version=2) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # load('results.RData') # # # Calculate BICs for all six sample sizes and all six models # BIC <- as.data.frame(matrix(,6,6)) # row.names(BIC) <- c('uniform','1-CPL','2-CPL','3-CPL','4-CPL','5-CPL') # for(s in 1:6){ # # # extract log likelihoods for each model # loglik <- c(uniform[[s]], # -CPL1[[s]]$value, # -CPL2[[s]]$value, # -CPL3[[s]]$value, # -CPL4[[s]]$value, # -CPL5[[s]]$value) # # # extract effective sample sizes for each model # N <- c(rep(ncol(PD[[s]]),6)) # # # number of parameters for each model # K <- c(0, 1, 3, 5, 7, 9) # # # calculate BIC for each model # BIC[,s] <- log(N)*K - 2*loglik # # # store effective sample size # names(BIC)[s] <- paste('N',N[1],sep='=') # } # # # Show all BICs for all sample sizes and models # print(BIC) ## ---- eval = FALSE------------------------------------------------------------ # oldpar <- par(no.readonly = TRUE) # # # Fig 2 plot # pdf('Fig2.pdf',height=6,width=13) # layout(mat=matrix(1:14, 2, 7, byrow = F),widths=c(0.3,rep(1,6)), heights=c(1,1.5),respect=T) # # # plot two blanks first # par(mar=c(5,4,1.5,0),las=2) # ymax <- 0.0032 # plot(NULL, xlim=c(0,1),ylim=c(0,1),main='', xlab='',ylab='',bty='n',xaxt='n',yaxt='n') # mtext(side=2, at=0.5,text='BIC',las=0,line=1) # plot(NULL, xlim=c(0,1),ylim=sqrt(c(0,ymax)),main='', xlab='',ylab='',bty='n',xaxt='n',yaxt='n') # axis(side=2, at=sqrt(seq(0,ymax,by=0.001)), labels=round(seq(0,ymax,by=0.001),4),las=1) # mtext(side=2, at=sqrt(0.00025),text='PD',las=0,line=0.8,cex=1) # abline(h=sqrt(seq(0,ymax,by=0.001)),col='grey') # # for(n in 1:6){ # # # extract the best model (lowest BIC) # BICs <- BIC[,n] # best <- which(BICs==min(BICs)) # # # convert parameters to model # if(best==1){ # type <- 'uniform' # pars <- NULL # } # if(best!=1)type <- 'CPL' # if(best==2)pars <- CPL1[[n]]$par # if(best==3)pars <- CPL2[[n]]$par # if(best==4)pars <- CPL3[[n]]$par # if(best==5)pars <- CPL4[[n]]$par # if(best==6)pars <- CPL5[[n]]$par # # spd.years <- as.numeric(row.names(SPD[[n]])) # spd.pdf <- SPD[[n]][,1] # mod.years <- as.numeric(row.names(PD[[n]])) # model <- convertPars(pars, mod.years, type) # # # # plot # red <- 'firebrick' # col <- rep('grey35',6); col[best] <- red # ymin <- min(BIC)-diff(range(BIC))*0.15 # par(mar=c(5,3,1.5,1),las=2) # plot(BICs,xlab='',ylab='',xaxt='n',pch=20,cex=3,col=col, main='') # axis(side=1, at=1:6, labels=c('Uniform','1-piece','2-piece','3-piece','4-piece','5-piece')) # # par(mar=c(5,1,1.5,1),las=2) # plot(NULL,type='l', # xlab='Cal Yrs BP', # ylab='',yaxt='n', # col='steelblue', # main=paste('N =',ncol(PD[[n]])), # ylim=sqrt(c(0,ymax)), # xlim=c(7500,5500)) # abline(h=sqrt(seq(0,ymax,by=0.001)),col='grey') # polygon(c(min(spd.years),spd.years,max(spd.years)),sqrt(c(0,spd.pdf,0)),col='steelblue',border=NA) # # lwd=3 # lines(toy$year,sqrt(toy$pdf),lwd=lwd) # lines(model$year, sqrt(model$pdf), lwd=lwd, col=red) # } # dev.off() # par(oldpar) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # set.seed(999) # # best exponential parameter previously found using ML search for Fig 5. # summary <- SPDsimulationTest(data=SAAD, # calcurve=shcal20, # calrange=c(2500,14000), # pars=-0.0001674152, # type='exp', # N=20000) # save(summary, file='results.RData',version=2) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # load('results.RData') # oldpar <- par(no.readonly = TRUE) # pdf('Fig4.pdf',height=4,width=10) # par(mar=c(2,4,0.1,0.1)) # plotSimulationSummary(summary, legend.x=11500,legend.y=0.0003) # axis(side=1, at=2500,labels='calBP',tick=F) # dev.off() # par(oldpar) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # library(DEoptimR) # # # Generate SPD # SPD <- summedPhaseCalibrator(data=SAAD, calcurve=shcal20, calrange = c(2500,14000)) # # # Calibrate each phase # CalArray <- makeCalArray(calcurve=shcal20, calrange = c(2500,14000)) # PD <- phaseCalibrator(data=SAAD, CalArray, remove.external = TRUE) # # # Best exponential model. Parameter and log likelihood found using seach # exp <- JDEoptim(lower=-0.01, upper=0.01, fn=objectiveFunction, PDarray=PD, type='exp', trace=T, NP=20) # # # Best CPL models. Parameters and log likelihood found using seach # fn <- objectiveFunction # CPL1 <- JDEoptim(lower=rep(0,1), upper=rep(1,1), fn, PDarray=PD, type='CPL',trace=T,NP=20) # CPL2 <- JDEoptim(lower=rep(0,3), upper=rep(1,3), fn, PDarray=PD, type='CPL',trace=T,NP=60) # CPL3 <- JDEoptim(lower=rep(0,5), upper=rep(1,5), fn, PDarray=PD, type='CPL',trace=T,NP=100) # CPL4 <- JDEoptim(lower=rep(0,7), upper=rep(1,7), fn, PDarray=PD, type='CPL',trace=T,NP=140) # CPL5 <- JDEoptim(lower=rep(0,9), upper=rep(1,9), fn, PDarray=PD, type='CPL',trace=T,NP=180) # CPL6 <- JDEoptim(lower=rep(0,11),upper=rep(1,11),fn, PDarray=PD, type='CPL',trace=T,NP=220) # # # save results, for separate plotting # save(SPD, PD, exp, CPL1, CPL2, CPL3, CPL4, CPL5, CPL6, file='results.RData',version=2) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # load('results.RData') # # # Calculate BICs for all six models # # name of each model # model <- c('exponential','1-CPL','2-CPL','3-CPL','4-CPL','5-CPL','6-CPL') # # # extract log likelihoods for each model # loglik <- c(-exp$value, -CPL1$value, -CPL2$value, -CPL3$value, -CPL4$value, -CPL5$value, -CPL6$value) # # # extract effective sample sizes # N <- c(rep(ncol(PD),7)) # # # number of parameters for each model # K <- c(1, 1, 3, 5, 7, 9, 11) # # # calculate BIC for each model # BICs <- log(N)*K - 2*loglik # # # convert best 3-CPL parameters into model pdf # best <- convertPars(pars=CPL3$par, years=c(2500:14000), type='CPL') ## ---- eval = FALSE------------------------------------------------------------ # oldpar <- par(no.readonly = TRUE) # pdf('Fig5.pdf',height=4,width=10) # par(mfrow=c(1,2)) # # # model comparison # par(mar=c(6,6,2,0.1)) # red <- 'firebrick' # blue <- 'steelblue' # col <- rep('grey35',7); col[which(BICs==min(BICs))] <- red # plot(BICs,xlab='',ylab='',xaxt='n', pch=20,cex=2,col=col,main='',las=1,cex.axis=0.7) # labels <- c('exponential','1-CPL','2-CPL','3-CPL','4-CPL','5-CPL','6-CPL') # axis(side=1, at=1:7, las=2, labels=labels, cex.axis=0.9) # mtext(side=2, at=mean(BICs),text='BIC',las=0,line=3) # # # best fitting CPL # years <- as.numeric(row.names(SPD)) # plot(NULL,xlim=rev(range(years)), ylim=range(SPD), # type='l',xlab='kyr cal BP',xaxt='n', ylab='',las=1,cex.axis=0.7) # axis(1,at=seq(14000,3000,by=-1000), labels=seq(14,3,by=-1),cex.axis=0.9) # mtext(side=2, at=max(SPD[,1])/2,text='PD',las=0,line=3.5,cex=1) # polygon(c(min(years),years,max(years)),c(0,SPD[,1],0),col=blue,border=NA) # lines(best$year,best$pdf,col=red,lwd=3) # # legend(x=14000,y=0.0003,lwd=c(5,3),col=c(blue,red),bty='n',legend=c('SPD','3-CPL')) # dev.off() # par(oldpar) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # library(DEoptimR) # # # Calibrate each phase # CalArray <- makeCalArray(calcurve=shcal20, calrange = c(2500,14000)) # PD <- phaseCalibrator(data=SAAD, CalArray, remove.external = TRUE) # # # arbitrary starting parameters # chain <- mcmc(PDarray=PD, startPars=rep(0.5,5), type='CPL', N=100000, burn=2000, thin=5, jumps=0.025) # # # find ML parameters # best.pars <- JDEoptim(lower=rep(0,5),upper=rep(1,5),fn=objectiveFunction,PDarray=PD,type='CPL',trace=T,NP=100)$par # # # save results, for separate plotting # save(chain, best.pars, file='results.RData',version=2) ## ---- eval = FALSE------------------------------------------------------------ # library('ADMUR') # library('scales') # load('results.RData') # # # Convert Maximum Likelihood parameters to hinge coordinates # ML <- CPLparsToHinges(best.pars,years=c(2500,14000)) # # # Convert MCMC chain of parameters to hinge coordinates # hinges <- CPLparsToHinges(chain$res, years=c(2500,14000)) # # # check the acceptance ratio is sensible (c. 0.2 to 0.5) # chain$acceptance.ratio # # # Eyeball the entire chain, before burn-in and thinning # for(n in 1:5)plot(chain$all.pars[,n], type='l', ylim=c(0,1)) # # # Generate CI for Fig 7 # N <- nrow(hinges) # years <- 2500:14000 # Y <- length(years) # pdf.matrix <- matrix(,N,Y) # for(n in 1:N){ # yr <- c('yr1','yr2','yr3','yr4') # pdf <- c('pdf1','pdf2','pdf3','pdf4') # pdf.matrix[n,] <- approx(x=hinges[n,yr],y=hinges[n,pdf],xout=years, ties='ordered')$y # } # CI <- matrix(,Y,6) # for(y in 1:Y)CI[y,] <- quantile(pdf.matrix[,y],prob=c(0.025,0.125,0.25,0.75,0.875,0.975)) ## ---- eval = FALSE------------------------------------------------------------ # oldpar <- par(no.readonly = TRUE) # pdf('Fig6.pdf',height=5,width=11) # par(mfrow=c(2,3)) # lwd <- 3 # red='firebrick' # grey='grey65' # breaks.yr <- seq(14000,2000,length.out=80) # breaks.pdf <- seq(0,0.0003,length.out=80) # xlab.yr <- 'yrs BP' # xlab.pdf <-'PD' # names <- c('Date of Hinge B','Date of Hinge C','PD of Hinge A','PD fo Hinge B','PD of Hinge C','PD of Hinge D') # hist(hinges$yr3, breaks=breaks.yr, col=grey, border=NA, main=names[1], xlab=xlab.yr) # abline(v = ML$year[3], col=red, lwd=lwd) # hist(hinges$yr2, breaks=breaks.yr, col=grey, border=NA, main=names[2], xlab=xlab.yr) # abline(v = ML$year[2], col=red, lwd=lwd) # hist(hinges$pdf4, breaks=breaks.pdf, col=grey, border=NA, main=names[3], xlab=xlab.pdf) # abline(v = ML$pdf[4], col=red, lwd=lwd) # hist(hinges$pdf3, breaks=breaks.pdf, col=grey, border=NA, main=names[5], xlab=xlab.pdf) # abline(v = ML$pdf[3], col=red, lwd=lwd) # hist(hinges$pdf2, breaks=breaks.pdf, col=grey, border=NA, main=names[4], xlab=xlab.pdf) # abline(v = ML$pdf[2], col=red, lwd=lwd) # hist(hinges$pdf1, breaks=breaks.pdf, col=grey, border=NA, main=names[6], xlab=xlab.pdf) # abline(v = ML$pdf[1], col=red, lwd=lwd) # dev.off() # par(oldpar) ## ---- eval = FALSE------------------------------------------------------------ # oldpar <- par(no.readonly = TRUE) # pdf('Fig7.pdf',height=5,width=12) # grey1 <- 'grey90' # grey2 <- 'grey70' # grey3 <- 'grey50' # red <- 'firebrick' # par(mfrow=c(1,2),las=0) # plot(NULL,xlim=c(14000,2500),ylim=c(0,0.00025),xlab='kyr cal BP',xaxt='n', ylab='PD', las=1, cex.axis=0.7) # set.seed(888) # S <- sample(1:N,size=1000) # for(n in 1:1000){ # lines(x=hinges[S[n],c('yr1','yr2','yr3','yr4')],y=hinges[S[n],c('pdf1','pdf2','pdf3','pdf4')],col=alpha('black',0.05)) # } # lines(ML$year, ML$pdf,col='firebrick',lwd=2) # axis(1,at=seq(14000,3000,by=-1000), labels=seq(14,3,by=-1)) # text(x=ML$year, y=ML$pdf + c(-0.00002,-0.00002,0.00002,0.00002), labels=rev(c('A','B','C','D'))) # legend(legend=c('Maximum Likelihood model PDF','Model PDF sampled from joint posterior parameters'), # x = 6000,y = 0.00024,cex = 0.7,bty = 'n',border = NA, xjust = 1, lwd=c(2,1), col=c(red,grey3)) # # plot(NULL,xlim=c(14000,2500),ylim=c(0,0.00025),xlab='kyr cal BP',xaxt='n', ylab='PD', las=1, cex.axis=0.7) # polygon(x=c(years,rev(years)),c(CI[,1],rev(CI[,6])),col=grey1,border=F) # polygon(x=c(years,rev(years)),c(CI[,2],rev(CI[,5])),col=grey2,border=F) # polygon(x=c(years,rev(years)),c(CI[,3],rev(CI[,4])),col=grey3,border=F) # a <- 0.05 # cex <- 0.2 # points(hinges$yr1,hinges$pdf1,pch=20,col=alpha(red,alpha=a),cex=cex) # points(hinges$yr2,hinges$pdf2,pch=20,col=alpha(red,alpha=a),cex=cex) # points(hinges$yr3,hinges$pdf3,pch=20,col=alpha(red,alpha=a),cex=cex) # points(hinges$yr4,hinges$pdf4,pch=20,col=alpha(red,alpha=a),cex=cex) # axis(1,at=seq(14000,3000,by=-1000), labels=seq(14,3,by=-1)) # legend(legend=c('Joint posterior parameters','50% CI of model PDF','75% CI of model PDF','95% CI of model PDF'), # x = 10000,y = 0.00024,cex = 0.7,bty = 'n',border = NA, xjust = 1, # pch = c(16,NA,NA,NA), # col = c(red,NA,NA,NA), # fill = c(NA,grey3,grey2,grey1), # x.intersp = c(1.5,1,1,1)) # dev.off() # par(oldpar) ## ---- eval = FALSE------------------------------------------------------------ # #---------------------------------------------------------------------------------------------- # # dates (H = hinge) # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # H.A.date <- ML$year[4] # H.B.date <- round(ML$year[3]) # H.C.date <- round(ML$year[2]) # H.D.date <- ML$year[1] # # H.B.date.CI <- round(quantile(hinges$yr3,prob=c(0.025,0.975))) # H.C.date.CI <- round(quantile(hinges$yr2,prob=c(0.025,0.975))) # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # # gradients (P = phase or piece) # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # P.1.gradient <- (ML$pdf[3] - ML$pdf[4]) / (ML$year[4] - ML$year[3]) # P.2.gradient <- (ML$pdf[2] - ML$pdf[3]) / (ML$year[3] - ML$year[2]) # P.3.gradient <- (ML$pdf[1] - ML$pdf[2]) / (ML$year[2] - ML$year[1]) # # P.1.gradient.mcmc <- (hinges$pdf3 - hinges$pdf4) / (hinges$yr4 - hinges$yr3) # P.2.gradient.mcmc <- (hinges$pdf2 - hinges$pdf3) / (hinges$yr3 - hinges$yr2) # P.3.gradient.mcmc <- (hinges$pdf1 - hinges$pdf2) / (hinges$yr2 - hinges$yr1) # # P.1.gradient.CI <- quantile(P.1.gradient.mcmc,prob=c(0.025,0.975)) # P.2.gradient.CI <- quantile(P.2.gradient.mcmc,prob=c(0.025,0.975)) # P.3.gradient.CI <- quantile(P.3.gradient.mcmc,prob=c(0.025,0.975)) # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # # relative growth rate per generation (P = phase or piece) # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # P.1.growth <- round(relativeRate(x=c(ML$year[3],ML$year[4]), y=c(ML$pdf[3],ML$pdf[4]) ),2) # P.2.growth <- round(relativeRate(x=c(ML$year[2],ML$year[3]), y=c(ML$pdf[2],ML$pdf[3]) ),2) # P.3.growth <- round(relativeRate(x=c(ML$year[1],ML$year[2]), y=c(ML$pdf[1],ML$pdf[2]) ),2) # # P.1.growth.mcmc <- relativeRate(x=hinges[,c('yr3','yr4')], y=hinges[,c('pdf3','pdf4')] ) # P.2.growth.mcmc <- relativeRate(x=hinges[,c('yr2','yr3')], y=hinges[,c('pdf2','pdf3')] ) # P.3.growth.mcmc <- relativeRate(x=hinges[,c('yr1','yr2')], y=hinges[,c('pdf1','pdf2')] ) # # P.1.growth.CI <- round(quantile(P.1.growth.mcmc,prob=c(0.025,0.975)),2) # P.2.growth.CI <- round(quantile(P.2.growth.mcmc,prob=c(0.025,0.975)),2) # P.3.growth.CI <- round(quantile(P.3.growth.mcmc,prob=c(0.025,0.975)),2) # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # # summary # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # headings <- c('Linear phase between hinges', # 'Start yrs BP (95% CI)','End yrs BP (95% CI)', # 'Gradient (x 10^-9 per year)(95% CI)', # 'Relative growth rate per 25 yr generation (95% CI)') # # all.dates <- c(H.A.date, # paste(H.B.date,' (',H.B.date.CI[2],' to ',H.B.date.CI[1],')',sep=''), # paste(H.C.date,' (',H.C.date.CI[2],' to ',H.C.date.CI[1],')',sep=''), # H.D.date) # # all.gradients <- round(c(P.1.gradient, P.2.gradient, P.3.gradient) / 1e-09, 1) # all.gradients.lower <- round(c(P.1.gradient.CI[1], P.2.gradient.CI[1], P.3.gradient.CI[1]) / 1e-09, 1) # all.gradients.upper <- round(c(P.1.gradient.CI[2], P.2.gradient.CI[2], P.3.gradient.CI[2]) / 1e-09, 1) # # col.1 <- c('1 (A-B)', # '2 (B-C)', # '3 (C-D)') # # col.2 <- all.dates[1:3] # # col.3 <- all.dates[2:4] # # col.4 <- c(paste(all.gradients[1],' (',all.gradients.lower[1],' to ',all.gradients.upper[1],')',sep=''), # paste(all.gradients[2],' (',all.gradients.lower[2],' to ',all.gradients.upper[2],')',sep=''), # paste(all.gradients[3],' (',all.gradients.lower[3],' to ',all.gradients.upper[3],')',sep='')) # # col.5 <- c(paste(P.1.growth,'%',' (',P.1.growth.CI[1],' to ',P.1.growth.CI[2],')',sep=''), # paste(P.2.growth,'%',' (',P.2.growth.CI[1],' to ',P.2.growth.CI[2],')',sep=''), # paste(P.3.growth,'%',' (',P.3.growth.CI[1],' to ',P.3.growth.CI[2],')',sep='')) # # res <- cbind(col.1,col.2,col.3,col.4,col.5); colnames(res) <- headings # write.csv(res, 'Table 2.csv', row.names=F) ## ---- eval = TRUE, echo = FALSE----------------------------------------------- tb2 <- read.csv(file='Table2.csv') print(tb2) ## ---- eval = FALSE------------------------------------------------------------ # library(ADMUR) # library(DEoptimR) # # # generate a set of random calendar dates under the toy model. # set.seed(888) # cal <- simulateCalendarDates(model = toy, 1500) # # # Convert to 14C dates. # age <- uncalibrateCalendarDates(cal, shcal20) # # # construct data frame. One date per phase. # data <- data.frame(age = age, sd = 25, phase = 1:1500, datingType = '14C') # # # Calibrate each phase # CalArray <- makeCalArray(shcal20, calrange = range(toy$year), inc = 5) # PD <- phaseCalibrator(data, CalArray, remove.external = TRUE) # # # Generate SPD # SPD <- summedCalibrator(data, CalArray) # # # Uniform model: No parameters. # # Log Likelihood calculated directly using objectiveFunction, without a search required. # unif.loglik <- -objectiveFunction(pars = NULL, PDarray = PD, type = 'uniform') # # # Best CPL models. Parameters and log likelihood found using seach # fn <- objectiveFunction # CPL1 <- JDEoptim(lower=rep(0,1), upper=rep(1,1), fn, PDarray=PD, type='CPL',trace=T,NP=20) # CPL2 <- JDEoptim(lower=rep(0,3), upper=rep(1,3), fn, PDarray=PD, type='CPL',trace=T,NP=60) # CPL3 <- JDEoptim(lower=rep(0,5), upper=rep(1,5), fn, PDarray=PD, type='CPL',trace=T,NP=100) # CPL4 <- JDEoptim(lower=rep(0,7), upper=rep(1,7), fn, PDarray=PD, type='CPL',trace=T,NP=140) # CPL5 <- JDEoptim(lower=rep(0,9), upper=rep(1,9), fn, PDarray=PD, type='CPL',trace=T,NP=180) # # # save results, for separate plotting # save(SPD, PD, unif.loglik, CPL1, CPL2, CPL3, CPL4, CPL5, file='results.RData',version=2) ## ---- eval = FALSE------------------------------------------------------------ # load('results.RData') # # # Calculate BICs for all six models # # name of each model # model <- c('uniform','1-CPL','2-CPL','3-CPL','4-CPL','5-CPL') # # # extract log likelihoods for each model # loglik <- c(unif.loglik, -CPL1$value, -CPL2$value, -CPL3$value, -CPL4$value, -CPL5$value) # # # extract effective sample sizes # N <- c(rep(ncol(PD),6)) # # # number of parameters for each model # K <- c(0, 1, 3, 5, 7, 9) # # # calculate BIC for each model # BIC <- log(N)*K - 2*loglik # # table <- data.frame(Model=model, Parameters=K, MaxLogLikelihood=loglik, BIC=BIC) # names(table) <- c('model','parameter','maximum log likelihood','BIC') # # print(table) # write.csv(table,file='Table 1.csv', row.names=F) ## ---- eval = TRUE, echo = FALSE----------------------------------------------- tb1 <- read.csv(file='Table1.csv') print(tb1)