inst/doc/Bchron.R

In Bchron: Radiocarbon Dating, Age-Depth Modelling, Relative Sea Level Rate Estimation, and Non-Parametric Phase Modelling

## ----global_options, include=FALSE--------------------------------------------
knitr::opts_chunk$set(fig.width=7, fig.height = 5, fig.align = 'center',
                      warning=FALSE, message=FALSE)

## ----eval=FALSE---------------------------------------------------------------
#  install.packages('Bchron')

## -----------------------------------------------------------------------------
library(Bchron)

## ----fig.align='center',fig.width=6,fig.height=5------------------------------
ages1 = BchronCalibrate(ages=11553,
                        ageSds=230,
                        calCurves='intcal20',
                        ids='Ox-123456')
summary(ages1)
plot(ages1)

## ----results='hide'-----------------------------------------------------------
ages2 = BchronCalibrate(ages=c(3445,11553,7456), 
                        ageSds=c(50,230,110), 
                        calCurves=c('intcal20','intcal20','shcal20'))
summary(ages2)
plot(ages2)

## ---- eval=FALSE--------------------------------------------------------------
#  plot(ages2, date = 'Date3')
#  plot(ages2, date = 1:2)

## ----fig.align='center',fig.width=6,fig.height=5------------------------------
ages3 = BchronCalibrate(ages=c(3445,11000), 
                        ageSds=c(50,200), 
                        positions=c(100,150), 
                        calCurves=c('intcal20','normal'))
summary(ages3)
plot(ages3)

## -----------------------------------------------------------------------------
library(ggplot2)
plot(ages3) +
  labs(x = 'Age (years BP)',
       y = 'Depth (cm)',
       title = 'Two dates at different depths')

## -----------------------------------------------------------------------------
plot(ages3, ageScale = 'bc', scaleReverse = FALSE) +
  labs(x = 'Age (years BC/AD)',
       y = 'Depth (cm)',
       title = 'Two dates at different depths')

## -----------------------------------------------------------------------------
plot(ages1, includeCal = TRUE)

## -----------------------------------------------------------------------------
ages4 = BchronCalibrate(ages=c(3445,11553,7456), 
                        ageSds=c(50,230,110), 
                        calCurves=rep('intcal20',3))
plot(ages4, includeCal = TRUE, fillCol = 'orange')

## -----------------------------------------------------------------------------
# First create age samples for each date
age_samples = sampleAges(ages3)
# Now summarise them with quantile - this gives a 95% credible interval
apply(age_samples, 2, quantile, prob=c(0.025,0.975))

## -----------------------------------------------------------------------------
apply(age_samples, 2, quantile, prob=c(0.5))

## -----------------------------------------------------------------------------
age_diff = age_samples[,2] - age_samples[,1]
qplot(age_diff, geom = 'histogram', bins = 30,
     main = 'Age difference between 3445+/-50 and 11553+/-230',
     ylab = 'Frequency', xlab = 'Age difference')

## -----------------------------------------------------------------------------
data(Glendalough)
print(Glendalough)

## ----results='hide'-----------------------------------------------------------
GlenOut = with(Glendalough, 
               Bchronology(ages=ages,
                           ageSds=ageSds, 
                           calCurves=calCurves,
                           positions=position, 
                           positionThicknesses=thickness,
                           ids=id, 
                           predictPositions=seq(0,1500,by=10)))

## ----eval=FALSE---------------------------------------------------------------
#  help(Bchronology)

## ----eval=FALSE---------------------------------------------------------------
#  summary(GlenOut)

## -----------------------------------------------------------------------------
summary(GlenOut, type='convergence')
summary(GlenOut, type='outliers')

## ----fig.align='center',fig.width=6,fig.height=5------------------------------
plot(GlenOut)

## ----results='hide'-----------------------------------------------------------
predictAges = predict(GlenOut, 
                      newPositions = c(150,725,1500), 
                      newPositionThicknesses=c(5,0,20))
predictAges = predict(GlenOut, 
                      newPositions = seq(0,1500,by=10))

## ---- results ='hide'---------------------------------------------------------
acc_rate = summary(GlenOut, type = 'acc_rate',
                   probs=c(0.25, 0.5, 0.75))

## ---- eval=FALSE--------------------------------------------------------------
#  acc_rate_2 = acc_rate %>%
#    gather(key = Percentile,
#           value = value,
#           -age_grid)
#  ggplot(acc_rate_2, aes(x = age_grid,
#                       y = value,
#                       colour = Percentile)) +
#    scale_colour_brewer(palette = 1) +
#    geom_line() + theme_bw() +
#    scale_x_reverse() +
#    labs(y = 'cm per year',
#         x = 'Age (cal years BP)')

## ---- eval=FALSE--------------------------------------------------------------
#  sed_rate = summary(GlenOut, type = 'sed_rate', useExisting = FALSE,
#                     probs=c(0.25, 0.5, 0.75))

## ---- eval=FALSE--------------------------------------------------------------
#  sed_rate_2 = sed_rate %>%
#    gather(key = Percentile,
#           value = value,
#           -position_grid)
#  ggplot(sed_rate_2, aes(x = position_grid,
#                       y = value,
#                       colour = Percentile)) +
#    scale_colour_brewer(palette = 1) +
#    geom_line() + theme_bw() +
#    scale_x_reverse() +
#    labs(y = 'Years per cm',
#         x = 'Depth (cm)')

## -----------------------------------------------------------------------------
summary(GlenOut, type = 'max_var')

## -----------------------------------------------------------------------------
max_var = summary(GlenOut, type = 'max_var', numPos = 1)
plot(GlenOut) + 
  labs(title = "Glendalough",
       x = 'Age (cal years BP)',
       y = 'Depth (cm)') + 
  geom_hline(yintercept = max_var)

## ---- eval=FALSE--------------------------------------------------------------
#  dateInfluence(GlenOut,
#                whichDate = 'Beta-100901',
#                measure = 'absMedianDiff')

## ---- eval=FALSE--------------------------------------------------------------
#  dateInfluence(GlenOut,
#                whichDate = 'all',
#                measure = 'absMedianDiff')

## ---- eval=FALSE--------------------------------------------------------------
#  dateInfluence(GlenOut,
#                whichDate = 'all',
#                measure = 'KL')

## -----------------------------------------------------------------------------
data(TestChronData)
data(TestRSLData)

## ----messages=FALSE, results='hide', eval=FALSE-------------------------------
#  RSLchron = with(TestChronData,
#                  Bchronology(ages = ages,
#                              ageSds = ageSds,
#                              positions = position,
#                              positionThicknesses = thickness,
#                              ids = id,
#                              calCurves = calCurves,
#                              predictPositions = TestRSLData$Depth))
#  RSLrun = with(TestRSLData,
#                BchronRSL(RSLchron,
#                          RSLmean = RSL,
#                          RSLsd = Sigma,
#                          degree = 3))

## ---- eval=FALSE--------------------------------------------------------------
#  summary(RSLrun, type = 'RSL', age_grid = seq(0, 2000, by  = 250))
#  plot(RSLrun, type = 'RSL') + ggtitle('Relative sea level plot')
#  plot(RSLrun, type = 'rate') + ggtitle('Rate of RSL change') +
#    ylab('Rate (mm per year)')

## ----results='hide', eval=FALSE-----------------------------------------------
#  data(Sluggan)
#  SlugDens = BchronDensity(ages=Sluggan$ages,
#                           ageSds=Sluggan$ageSds,
#                           calCurves=Sluggan$calCurves)

## ---- eval=FALSE--------------------------------------------------------------
#  summary(SlugDens, prob = 0.95)

## ---- eval=FALSE--------------------------------------------------------------
#  plot(SlugDens,xlab='Age (cal years BP)', xaxp=c(0, 16000, 16))

## ----eval=FALSE---------------------------------------------------------------
#  SlugDensFast = BchronDensityFast(ages=Sluggan$ages,
#                                   ageSds=Sluggan$ageSds,
#                                   calCurves=Sluggan$calCurves)

## ---- eval = FALSE------------------------------------------------------------
#  # Load in the calibration curve with:
#  intcal09 = read.table(system.file('extdata/intcal09.14c',
#                                    package = 'Bchron'), sep=',')
#  
#  # Run createCalCurve
#  createCalCurve(name='intcal09',calAges=intcal09[,1],
#                 uncalAges=intcal09[,2],oneSigma=intcal09[,3])

## ---- eval = FALSE------------------------------------------------------------
#  file.copy('intcal09.rda', system.file('data',package='Bchron'))

## ---- eval = FALSE------------------------------------------------------------
#  age_09 = BchronCalibrate(age=15500,ageSds=150,calCurves = 'intcal09',
#                           ids='IntCal09')
#  age_20 = BchronCalibrate(age=15500,ageSds=150,calCurves = 'intcal20',
#                           ids = 'Intcal20')
#  library(ggplot2)
#  plot(age_09) +
#    geom_line(data = as.data.frame(age_20$Intcal20),
#              aes(x = ageGrid, y = densities), col = 'red') +
#    ggtitle('Intcal09 vs Intcal20')

## ---- results = 'hide'--------------------------------------------------------
unCal1 = unCalibrate(2350, type = 'ages')

## -----------------------------------------------------------------------------
print(unCal1)

## ---- results = 'hide'--------------------------------------------------------
unCal2 = unCalibrate(calAge = c(2350, 4750, 11440),
                     type = 'ages',
                     calCurve = 'intcal20')

## -----------------------------------------------------------------------------
print(unCal2)

## ---- results = 'hide'--------------------------------------------------------
ageRange = seq(8000, 9000, by = 5)
c14Ages = unCalibrate(ageRange,
                      type = 'ages')
load(system.file('data/intcal20.rda', package = 'Bchron'))
ggplot(intcal20, aes(x = V1, y = V2)) + 
  geom_line() + 
  theme_bw() + 
  scale_x_continuous(limits = range(ageRange)) + 
  scale_y_continuous(limits = range(c14Ages)) + 
  geom_rug(data = data.frame(c14Ages), aes(y = c14Ages), 
           sides = 'l', inherit.aes = F) + 
  geom_rug(data = data.frame(ageRange), aes(x = ageRange), 
           sides = 'b', inherit.aes = F) + 
  labs(x = 'Calendar years BP', y = '14C years BP')

## -----------------------------------------------------------------------------
calAge = BchronCalibrate(ages = 11255,
                         ageSds = 25,
                         calCurves = 'intcal20')
calSampleAges = sampleAges(calAge)

## ---- results = 'hide'--------------------------------------------------------
unCal = unCalibrate(calSampleAges,
            type = 'samples')

## -----------------------------------------------------------------------------
print(unCal)

## ---- results = 'hide'--------------------------------------------------------
data(Glendalough)
GlenOut = with(Glendalough, 
               Bchronology(ages=ages,
                           ageSds=ageSds, 
                           calCurves=calCurves,
                           positions=position, 
                           positionThicknesses=thickness,
                           ids=id, 
                           predictPositions=seq(0,1500,by=10)))

## -----------------------------------------------------------------------------
library(ggplot2)
plot(GlenOut) +
labs(title = "Glendalough",
     xlab = "Age (cal years BP)",
     ylab = "Depth (cm)")

## -----------------------------------------------------------------------------
new_cal = BchronCalibrate(ages = 7000, 
                          ageSds = 40, 
                          calCurve = 'intcal20')

## -----------------------------------------------------------------------------
library(ggridges)
plot(GlenOut) +
  geom_ridgeline(data = as.data.frame(new_cal$Date1), 
                 aes(x = ageGrid, 
                     y = 600,
                     height = densities*10000, # Note the 10000 came from trial and error
                     group = 'New date',
                     ),
                 fill = 'grey',
                 colour = 'black') +
  annotate("text", x = 9000, y = 570, label = "New date")