,
Enrico Crema [aut, cre] ,
R. Kyle Bocinsky [ctb],
Martin Hinz [ctb],
Philip Riris [ctb],
Fabio Silva [ctb]| Title: | Calibration and Analysis of Radiocarbon Dates |
|---|---|
| Description: | Enables the calibration and analysis of radiocarbon dates, often but not exclusively for the purposes of archaeological research. It includes functions not only for basic calibration, uncalibration, and plotting of one or more dates, but also a statistical framework for building demographic and related longitudinal inferences from aggregate radiocarbon date lists, including: Monte-Carlo simulation test (Timpson et al 2014 <doi:10.1016/j.jas.2014.08.011>), random mark permutation test (Crema et al 2016 <doi:10.1371/journal.pone.0154809>) and spatial permutation tests (Crema, Bevan, and Shennan 2017 <doi:10.1016/j.jas.2017.09.007>). |
| Authors: | Andrew Bevan [aut] ,
Enrico Crema [aut, cre] ,
R. Kyle Bocinsky [ctb],
Martin Hinz [ctb],
Philip Riris [ctb],
Fabio Silva [ctb] |
| Maintainer: | Enrico Crema <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.5.1 |
| Built: | 2024-11-09 05:34:24 UTC |
| Source: | https://github.com/ahb108/rcarbon |
Convert other calibrated date formats to an rcarbon CalDates object.
as.CalDates(x)as.CalDates(x)
x |
One or more calibrated dated to convert (currently only BchronCalibratedDates and oxcAARCalibratedDatesList objects are supported) |
A CalDates object
## Not run: library(Bchron) library(oxcAAR) quickSetupOxcal() dates <- data.frame(CRA=c(3200,2100,1900), Error=c(35,40,50)) bcaldates <- BchronCalibrate(ages=dates$CRA, ageSds=dates$Error, calCurves=rep("intcal13", nrow(dates))) rcaldates <- rcarbon::calibrate(dates$CRA, dates$Error, calCurves=rep("intcal13")) ocaldates <- oxcalCalibrate(c(3200,2100,1900),c(35,40,50),c("a","b","c")) ## Convert to rcarbon format caldates.b <- as.CalDates(bcaldates) caldates.o <- as.CalDates(ocaldates) ## Comparison plot plot(rcaldates$grids[[2]]$calBP,rcaldates$grids[[2]]$PrDens, type="l", col="green", xlim=c(2300,1900)) lines(caldates.b$grids[[2]]$calBP,caldates.b$grids[[2]]$PrDens, col="red") lines(caldates.o$grids[[2]]$calBP,caldates.o$grids[[2]]$PrDens, col="blue") legend("topright", legend=c("rcarbon","Bchron","OxCal"), col=c("green","red","blue"), lwd=2) ## End(Not run)## Not run: library(Bchron) library(oxcAAR) quickSetupOxcal() dates <- data.frame(CRA=c(3200,2100,1900), Error=c(35,40,50)) bcaldates <- BchronCalibrate(ages=dates$CRA, ageSds=dates$Error, calCurves=rep("intcal13", nrow(dates))) rcaldates <- rcarbon::calibrate(dates$CRA, dates$Error, calCurves=rep("intcal13")) ocaldates <- oxcalCalibrate(c(3200,2100,1900),c(35,40,50),c("a","b","c")) ## Convert to rcarbon format caldates.b <- as.CalDates(bcaldates) caldates.o <- as.CalDates(ocaldates) ## Comparison plot plot(rcaldates$grids[[2]]$calBP,rcaldates$grids[[2]]$PrDens, type="l", col="green", xlim=c(2300,1900)) lines(caldates.b$grids[[2]]$calBP,caldates.b$grids[[2]]$PrDens, col="red") lines(caldates.o$grids[[2]]$calBP,caldates.o$grids[[2]]$PrDens, col="blue") legend("topright", legend=c("rcarbon","Bchron","OxCal"), col=c("green","red","blue"), lwd=2) ## End(Not run)
Tries to coerce any two-column matrix or data.frame to a calibrated probability distribution (an object of class "CalGrid") for use by the rcarbon package.
as.CalGrid(x)as.CalGrid(x)
x |
A two-column |
A CalGrid class object of probabilities or summed probabilities per calendar year BP.
df <- data.frame(calBP=5000:2000,PrDens=runif(length(5000:2000))) mycalgrid <- as.CalGrid(df) plot(mycalgrid)df <- data.frame(calBP=5000:2000,PrDens=runif(length(5000:2000))) mycalgrid <- as.CalGrid(df) plot(mycalgrid)
Tries to coerce any two-column matrix or data.frame to a uncalibrated probability distribution (an object of class "UncalGrid") for use by the rcarbon package.
as.UncalGrid(x)as.UncalGrid(x)
x |
A two-column |
A CalGrid class object of probabilities or summed probabilities per CRA.
df <- data.frame(CRA=5000:2000,PrDens=runif(length(5000:2000))) mycalgrid <- as.UncalGrid(df)df <- data.frame(CRA=5000:2000,PrDens=runif(length(5000:2000))) mycalgrid <- as.UncalGrid(df)
Plot the median values of multiple calibrated radiocarbon dates or bins in a barcode-like strip.
barCodes( x, yrng = c(0, 0.03), width = 20, col = rgb(0, 0, 0, 25, maxColorValue = 255), border = NA, ... )barCodes( x, yrng = c(0, 0.03), width = 20, col = rgb(0, 0, 0, 25, maxColorValue = 255), border = NA, ... )
x |
A vector containing median values obtained from |
yrng |
y-axis range of the bars. |
width |
width of the bars (optional) |
col |
color of the bars |
border |
the color to draw the border. Use border = NA to omit borders. |
... |
Additional arguments affecting the plot |
## Not run: #Load EUROEVOL Data data(euroevol) #Subset Danish Dates denmark <- subset(euroevol,Country=="Denmark") #Calibrate and Bin denmarkDates <- calibrate(x=denmark$C14Age,errors=denmark$C14SD) denmarkBins <- binPrep(sites=denmark$SiteID,ages=denmark$C14Age,h=200) #200 years bin size #Compute median date for each bin bm <- binMed(x=denmarkDates,bins=denmarkBins) #Compute median date for each date dm <- medCal(denmarkDates) #Compute SPD denmarkSPD <- spd(x=denmarkDates,bins=denmarkBins,timeRange=c(10000,4000)) #Plot SPD and barCodes of median dates plot(denmarkSPD,runm=200) barCodes(dm,yrng=c(0,0.01)) #Plot SPD and barCodes of median bins in BC/AD plot(denmarkSPD,runm=200,calendar="BCAD") barCodes(BPtoBCAD(bm),yrng=c(0,0.01)) ## End(Not run)## Not run: #Load EUROEVOL Data data(euroevol) #Subset Danish Dates denmark <- subset(euroevol,Country=="Denmark") #Calibrate and Bin denmarkDates <- calibrate(x=denmark$C14Age,errors=denmark$C14SD) denmarkBins <- binPrep(sites=denmark$SiteID,ages=denmark$C14Age,h=200) #200 years bin size #Compute median date for each bin bm <- binMed(x=denmarkDates,bins=denmarkBins) #Compute median date for each date dm <- medCal(denmarkDates) #Compute SPD denmarkSPD <- spd(x=denmarkDates,bins=denmarkBins,timeRange=c(10000,4000)) #Plot SPD and barCodes of median dates plot(denmarkSPD,runm=200) barCodes(dm,yrng=c(0,0.01)) #Plot SPD and barCodes of median bins in BC/AD plot(denmarkSPD,runm=200,calendar="BCAD") barCodes(BPtoBCAD(bm),yrng=c(0,0.01)) ## End(Not run)
Converts BC/AD dates to BP format while handling the absence of 'year 0'
BCADtoBP(x)BCADtoBP(x)
x |
A numerical vector (currently only basic checks that these numbers are in a sensible range). |
A vector with BP dates.
BCADtoBP(-1268)BCADtoBP(-1268)
Function for generating a vector of median calibrated dates for each each bin.
binMed(x, bins, verbose = TRUE)binMed(x, bins, verbose = TRUE)
x |
A |
bins |
vector containing the bin names associated with each radiocarbon date. Can be generated using |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
A vector of median dates in cal BP
## Not run: #Load EUROEVOL Data data(euroevol) #Subset Danish Dates denmark <- subset(euroevol,Country=="Denmark") #Calibrate and Bin denmarkDates <- calibrate(x=denmark$C14Age,errors=denmark$C14SD) denmarkBins <- binPrep(sites=denmark$SiteID,ages=denmark$C14Age,h=200) #200 years bin size #Compute median date for each bin binMed(x=denmarkDates,bins=denmarkBins) ## End(Not run)## Not run: #Load EUROEVOL Data data(euroevol) #Subset Danish Dates denmark <- subset(euroevol,Country=="Denmark") #Calibrate and Bin denmarkDates <- calibrate(x=denmark$C14Age,errors=denmark$C14SD) denmarkBins <- binPrep(sites=denmark$SiteID,ages=denmark$C14Age,h=200) #200 years bin size #Compute median date for each bin binMed(x=denmarkDates,bins=denmarkBins) ## End(Not run)
Prepare a set of bins for controlling the aggregation of radiocarbon dates
known to be from the same phase of same archaeological site (for use with spd). Used in cases where there is a concern that unusually high levels of sampling for radiocarbon at a given site or in a given site phase will impede comparison between sites or phases.
binPrep(sites, ages, h, method = "complete")binPrep(sites, ages, h, method = "complete")
sites |
a vector of character strings (or number to coerce to character) of all sites or site phases. If character strings are used these should not contain underscores (see also below) |
ages |
a vector of uncalibrated conventional radiocarbon ages or a |
h |
a single numeric value passed to |
method |
the agglomeration method to be used, passed on to |
If ages is a CalDates class object, median dates are used for the clustering.
A vector of character strings with the same length of the object supplied for the argument ages identifying intra-site or intra-phase grouping, for use with spd.The character strings effectively provide a "name" for each "phase" within a "site", using sequential integers after an underscore. For example if a site named "S001" had four dates grouped into two bins with two dates each, the resulting vector would be "S001_1", "S001_1", "S001_2", and "S001_2".
spd for generating SPD; binsense for sensitivity analysis pertaining the choice of the parameter h.
Visually explores how choosing different values for h in the binPrep function affects the shape of the SPD.
binsense( x, y, h, timeRange, calendar = "BP", binning = "CRA", raw = F, verbose = T, legend = T, ... )binsense( x, y, h, timeRange, calendar = "BP", binning = "CRA", raw = F, verbose = T, legend = T, ... )
x |
A |
y |
A vector containing the locations ID (e.g. site ID) of each calibrated date to be used for the binning process. |
h |
A vector of numbers containing values for the |
timeRange |
A vector of length 2 indicating the start and end date of the analysis in cal BP. |
calendar |
Either |
binning |
Either |
raw |
A logical variable indicating whether all SPDs should be returned or not. Default is FALSE. |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
legend |
A logical variable indicating whether the legend should be displayed. Default is TRUE |
... |
Additional arguments to be passed to the |
## Not run: data(euroevol) #subset Danish dates denmark=subset(euroevol,Country=="Denmark") denmarkDates=calibrate(x=denmark$C14Age,errors=denmark$C14SD) binsense(x=denmarkDates,y=denmark$SiteID,h=seq(0,200,20),timeRange=c(10000,4000),runm=200) ## End(Not run)## Not run: data(euroevol) #subset Danish dates denmark=subset(euroevol,Country=="Denmark") denmarkDates=calibrate(x=denmark$C14Age,errors=denmark$C14SD) binsense(x=denmarkDates,y=denmark$SiteID,h=seq(0,200,20),timeRange=c(10000,4000),runm=200) ## End(Not run)
Converts calibrated BP dates to BC/AD dates, omitting ‘year 0’
BPtoBCAD(x)BPtoBCAD(x)
x |
A numerical vector (currently only basic checks that these numbers are in a sensible range). |
A vector with BC/BCE dates expressed as negative numbers and AD/CE dates as positive ones.
BPtoBCAD(4200)BPtoBCAD(4200)
Function for calibrating radiocarbon dates or uncalibrated SPDs.
Function for calibrating a SPD generated by summing uncalibrated dates.
calibrate(x, ...) ## Default S3 method: calibrate( x, errors, ids = NA, dateDetails = NA, calCurves = "intcal20", resOffsets = 0, resErrors = 0, timeRange = c(55000, 0), normalised = TRUE, F14C = FALSE, calMatrix = FALSE, eps = 1e-05, ncores = 1, verbose = TRUE, ... ) ## S3 method for class 'UncalGrid' calibrate( x, errors = 0, calCurves = "intcal20", timeRange = c(55000, 0), eps = 1e-05, type = "fast", datenormalised = FALSE, spdnormalised = FALSE, verbose = TRUE, ... )calibrate(x, ...) ## Default S3 method: calibrate( x, errors, ids = NA, dateDetails = NA, calCurves = "intcal20", resOffsets = 0, resErrors = 0, timeRange = c(55000, 0), normalised = TRUE, F14C = FALSE, calMatrix = FALSE, eps = 1e-05, ncores = 1, verbose = TRUE, ... ) ## S3 method for class 'UncalGrid' calibrate( x, errors = 0, calCurves = "intcal20", timeRange = c(55000, 0), eps = 1e-05, type = "fast", datenormalised = FALSE, spdnormalised = FALSE, verbose = TRUE, ... )
x |
A vector of uncalibrated radiocarbon ages or a UncalGrid class object. |
... |
ignored |
errors |
A vector of standard deviations corresponding to each estimated radiocarbon age. |
ids |
An optional vector of IDs for each date. |
dateDetails |
An optional vector of details for each date which will be returned in the output metadata. |
calCurves |
Either a character string naming a calibration curve already provided with the rcarbon package (currently 'intcal20','intcal13','intcal13nhpine16','shcal20','shcal13','shcal13shkauri16',”marine13','marine20' and 'normal'(i.e. no calibration) are possible; default is 'intcal20') or a custom calibration curve as three-column matrix or data.frame (calibrated year BP, uncalibrated age bp, standard deviation). Different existing curves can be specified per dated sample, but only one custom curve can be provided for all dates. |
resOffsets |
A vector of offset values for any marine reservoir effect (default is no offset). |
resErrors |
A vector of offset value errors for any marine reservoir effect (default is no offset). |
timeRange |
Earliest and latest data to calibrate for, in calendar years. Posterior probabilities beyond this range will be excluded (the default is sensible in most cases). |
normalised |
A logical variable indicating whether the calibration should be normalised or not. Default is TRUE. |
F14C |
A logical variable indicating whether calibration should be carried out in F14C space or not. Default is FALSE. |
calMatrix |
a logical variable indicating whether the age grid should be limited to probabilities higher than |
eps |
Cut-off value for density calculation. Default is 1e-5. |
ncores |
Number of cores/workers used for parallel execution. Default is 1 (>1 requires doSnow package). |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
type |
Currently, two options are available. With "fast", a look-up is conducted for each calendar year that picks up the amplitude of the nearest CRA age. With "full", each CRA is calibrated individually and then summed (much slower). |
datenormalised |
Controls for calibrated dates with probability mass outside the timerange of analysis. If set to TRUE the total probability mass within the time-span of analysis is normalised to sum to unity. Should be set to FALSE when the parameter |
spdnormalised |
A logical variable indicating whether the total probability mass of the SPD is normalised to sum to unity. |
This function computes one or more calibrated radiocarbon ages using the method described in Bronk Ramsey 2008 (see also Parnell 2017). It is possible to specify different calibration curves or reservoir offsets individually for each date, and control whether the resulting calibrated distribution is normalised to 1 under-the-curve or not. Calculations can also be executed in parallel to reduce computing time. The function was modified from the BchronCalibrate function in the Bchron package developed by A.Parnell (see references below).
An object of class CalDates with the following elements:
metadata A data.frame containing relevant information regarding each radiocarbon date and the parameter used in the calibration process.
grids A list of calGrid class objects, containing the posterior probabilities for each calendar year. The most memory-efficient way to store calibrated dates, as only years with non-zero probability are stored, but aggregation methods such as spd() may then take longer to extract and combine multiple dates. NA when the parameter calMatrix is set to TRUE.
calMatrix A matrix of probability values, one row per calendar year in timeRange and one column per date. By storing all possible years, not just those with non-zero probability, this approach takes more memory, but speeds up spd() and is suggested whenever the latter is to be used. NA when the parameter calMatrix is set to FALSE.
Bronk Ramsey, C. 2008. Radiocarbon dating: revolutions in understanding, Archaeometry 50.2: 249-75. DOI: https://doi.org/10.1111/j.1475-4754.2008.00394.x
Parnell, A. 2017. Bchron: Radiocarbon Dating, Age-Depth Modelling, Relative Sea Level Rate Estimation, and Non-Parametric Phase Modelling, R package: https://CRAN.R-project.org/package=Bchron
x1 <- calibrate(x=4000, errors=30) plot(x1) summary(x1) # Example with a Marine Date, using a DeltaR of 300 and a DeltaR error of 30 x2 <- calibrate(x=4000, errors=30, calCurves='marine20', resOffsets=300, resErrors=30) plot(x2)x1 <- calibrate(x=4000, errors=30) plot(x1) summary(x1) # Example with a Marine Date, using a DeltaR of 300 and a DeltaR error of 30 x2 <- calibrate(x=4000, errors=30, calCurves='marine20', resOffsets=300, resErrors=30) plot(x2)
Computes a Composite Kernel Density Estimate (CKDE) from multiple sets of randomly sampled calendar dates.
ckde(x, timeRange, bw, normalised = FALSE)ckde(x, timeRange, bw, normalised = FALSE)
x |
A |
timeRange |
A vector of length 2 indicating the start and end date of the analysis in cal BP. |
bw |
Kernel bandwidth to be used. |
normalised |
A logical variable indicating whether the contribution of individual dates should be equal (TRUE), or weighted based on the area under the curve of non-normalised calibration (FALSE). Default is TRUE. |
The function computes Kernel Density Estimates using randomly sampled calendar dates contained in a simdates class object (generated using the sampledates() function). The output contains nsim KDEs, where nsim is the argument used in simulate.dates(). The resulting object can be plotted to visualise a CKDE (cf Brown 2017), and if boot was set to TRUE in sampleDates its bootstrapped variant (cf McLaughlin 2018 for a similar analysis). The shape of the CKDE is comparable to an SPD generated from non-normalised dates when the argument normalised is set to FALSE.
An object of class ckdeSPD with the following elements
timeRange The timeRange setting used.
res.matrix A matrix containing the KDE values with rows representing calendar dates.
Brown, W. A. 2017. The past and future of growth rate estimation in demographic temporal frequency analysis: Biodemographic interpretability and the ascendance of dynamic growth models. Journal of Archaeological Science, 80: 96–108. DOI: https://doi.org/10.1016/j.jas.2017.02.003
McLaughlin, T. R. 2018. On Applications of Space–Time Modelling with Open-Source 14C Age Calibration. Journal of Archaeological Method and Theory. DOI https://doi.org/10.1007/s10816-018-9381-3
## Not run: data(emedyd) x = calibrate(x=emedyd$CRA, errors=emedyd$Error,normalised=FALSE) bins = binPrep(sites=emedyd$SiteName, ages=emedyd$CRA,h=50) s = sampleDates(x,bins=bins,nsim=100,boot=FALSE) ckdeNorm = ckde(s,timeRange=c(16000,9000),bw=100,normalised=TRUE) plot(ckdeNorm,type='multiline',calendar='BCAD') ## End(Not run)## Not run: data(emedyd) x = calibrate(x=emedyd$CRA, errors=emedyd$Error,normalised=FALSE) bins = binPrep(sites=emedyd$SiteName, ages=emedyd$CRA,h=50) s = sampleDates(x,bins=bins,nsim=100,boot=FALSE) ckdeNorm = ckde(s,timeRange=c(16000,9000),bw=100,normalised=TRUE) plot(ckdeNorm,type='multiline',calendar='BCAD') ## End(Not run)
Combine multiple CalDates Class Objects into one.
combine(..., fixIDs = FALSE)combine(..., fixIDs = FALSE)
... |
|
fixIDs |
logical. If set to TRUE, each date is given a new ID based on sequential integer. Default is FALSE |
An object of class CalDates
x1 = calibrate(c(2000,3400),c(20,20),ids=1:2) x2 = calibrate(c(4000,3000),c(30,30),calCurves=c('intcal20','marine20'), resOffsets=c(0,30),resErrors=c(0,20),ids=3:4) mcurve <- mixCurves('intcal20','marine20',p=0.7,resOffsets=300,resErrors=20) x3 = calibrate(5300,20,calCurves=mcurve,ids=5) x = combine(x1,x2,x3) ## x$metadatax1 = calibrate(c(2000,3400),c(20,20),ids=1:2) x2 = calibrate(c(4000,3000),c(30,30),calCurves=c('intcal20','marine20'), resOffsets=c(0,30),resErrors=c(0,20),ids=3:4) mcurve <- mixCurves('intcal20','marine20',p=0.7,resOffsets=300,resErrors=20) x3 = calibrate(5300,20,calCurves=mcurve,ids=5) x = combine(x1,x2,x3) ## x$metadata
Radiocarbon dates (n=1915) and site coordinates (n=201) from a paper considering the relationship between human activity in the eastern Mediterranean/Middle East and early Holocene climate change, including the Younger Dryas.
emedydemedyd
A data.frame with the following variables:
LabIDLaboratory ID assigned to each radiocarbon date (where known)
CRARadiocarbon age in 14C years BP
ErrorRadiocarbon age error
MaterialMaterial of the dated sample
SpeciesSpecies of the dated sample (where identified)
SiteNameName of the site from which the sample has been recovered
CountryCountry where the sampling site is located
LongitudeLongitude of the sampling site in decimal degrees
LatitudeLatitude of the sampling site in decimal degrees
RegionOne of three analytical regions (1=southern Levant, 2=Northern Levant, 3= South-central Anatolia
Palmisano, A., Bevan, A. and S. Shennan 2017. Data and code for demographic trends in the paper "Human responses and non-responses to climatic variations during the Last Glacial-Interglacial transition in the eastern Mediterranean", UCL Discovery Archive 1570274. doi:10.14324/000.ds.1570274.
Roberts, N., Woodbridge, J., Bevan, A., Palmisano, A., Shennan, S. and E. Asouti 2017. Human responses and non-responses to climatic variations during the Last Glacial-Interglacial transition in the eastern Mediterranean. Quaternary Science Reviews, 184, 47-67. doi:10.1016/j.quascirev.2017.09.011.
## Not run: data(emedyd) northernlevant <- emedyd[emedyd$Region=="2",] bins <- binPrep(northernlevant$SiteName, northernlevant$CRA, h=50) x <- calibrate(northernlevant$CRA, northernlevant$Error, normalised=FALSE) spd.northernlevant <- spd(x, bins=bins, runm=50, timeRange=c(17000,8000)) plot(spd.northernlevant) ## End(Not run)## Not run: data(emedyd) northernlevant <- emedyd[emedyd$Region=="2",] bins <- binPrep(northernlevant$SiteName, northernlevant$CRA, h=50) x <- calibrate(northernlevant$CRA, northernlevant$Error, normalised=FALSE) spd.northernlevant <- spd(x, bins=bins, runm=50, timeRange=c(17000,8000)) plot(spd.northernlevant) ## End(Not run)
Radiocarbon dates (n=14,053) and site coordinates (n=4,213) from the EUROEVOL project database. Sites without radiocarbon dates (n=544), phase-codes, and other data have been omitted.
euroevoleuroevol
A data.frame with the following variables:
C14IDID of each radiocarbon date
C14AgeRadiocarbon age in 14C years BP
C14SDRadiocarbon age error
LabCodeLabcode of the radiocarbon date
MaterialMaterial of the dated sample
SiteIDID of the site from which the sample has been recovered
LatitudeLatitude of the sampling site in decimal degrees
LongitudeLongitude of the sampling site in decimal degrees
CountryCountry where the sampling site is located
Manning, K., Colledge, S., Crema, E., Shennan, S., Timpson, A., 2016. The Cultural Evolution of Neolithic Europe. EUROEVOL Dataset 1: Sites, Phases and Radiocarbon Data. Journal of Open Archaeology Data 5. doi:10.5334/joad.40
Shennan, S., Downey, S.S., Timpson, A., Edinborough, K., Colledge, S., Kerig, T., Manning, K., Thomas, M.G., 2013. Regional population collapse followed initial agriculture booms in mid-Holocene Europe. Nature Communications 4, ncomms3486. doi:10.1038/ncomms3486
Timpson, A., Colledge, S., Crema, E., Edinborough, K., Kerig, T., Manning, K., Thomas, M.G., Shennan, S., 2014. Reconstructing regional population fluctuations in the European Neolithic using radiocarbon dates: a new case-study using an improved method. Journal of Archaeological Science 52, 549-557. doi:10.1016/j.jas.2014.08.011
## Not run: data(euroevol) Ireland <- subset(euroevol,Country=="Ireland") bins <- binPrep(Ireland$SiteID,Ireland$C14Age,h=200) x <- calibrate(Ireland$C14Age,Ireland$C14SD) spd.ireland <- spd(x,bins=bins,runm=200,timeRange=c(8000,4000)) plot(spd.ireland) ## End(Not run)## Not run: data(euroevol) Ireland <- subset(euroevol,Country=="Ireland") bins <- binPrep(Ireland$SiteID,Ireland$C14Age,h=200) x <- calibrate(Ireland$C14Age,Ireland$C14SD) spd.ireland <- spd(x,bins=bins,runm=200,timeRange=c(8000,4000)) plot(spd.ireland) ## End(Not run)
Radiocarbon dates (n=2,324) and site coordinates (n=652) from England and Wales collected from the EUROEVOL project database. See euroevol for more details regarding the source data.
ewdatesewdates
A data.frame with the following variables:
C14IDID of each radiocarbon date
C14AgeRadiocarbon age in 14C years BP
C14SDRadiocarbon age error
LabCodeLabcode of the radiocarbon date
MaterialMaterial of the dated sample
SiteIDID of the site from which the sample has been recovered
EastingsEasting coordinates of the sampling site in meters (OSGB 1936 epsg:27700)
NorthingsNorthing coordinates of the sampling site in meters (OSGB 1936 epsg:27700)
Manning, K., Colledge, S., Crema, E., Shennan, S., Timpson, A., 2016. The Cultural Evolution of Neolithic Europe. EUROEVOL Dataset 1: Sites, Phases and Radiocarbon Data. Journal of Open Archaeology Data 5. doi:10.5334/joad.40
An owin class polygonal window of England and Wales.
ewowinewowin
An owin class object.
South,A.2011.rworldmap: A New R package for Mapping Global Data. The R Journal Vol. 3/1 : 35-43.
## Not run: data(ewowin) # Obtained from rworldmap: library(maptools) library(rgeos) library(rworldmap) bng <- CRS("+init=epsg:27700") sbrit <- getMap(resolution="high") sbrit <- spTransform(sbrit[sbrit$GEOUNIT =="United Kingdom",], bng) tmp <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(130000,130000,310000, 425000,700000,700000,130000), c(0,300000,560000,598000,300000,0,0)))), "1")), 1:1, proj4string=bng) sbrit <- gIntersection(sbrit,tmp) ewowin <- as.owin(sbrit) ## End(Not run)## Not run: data(ewowin) # Obtained from rworldmap: library(maptools) library(rgeos) library(rworldmap) bng <- CRS("+init=epsg:27700") sbrit <- getMap(resolution="high") sbrit <- spTransform(sbrit[sbrit$GEOUNIT =="United Kingdom",], bng) tmp <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(130000,130000,310000, 425000,700000,700000,130000), c(0,300000,560000,598000,300000,0,0)))), "1")), 1:1, proj4string=bng) sbrit <- gIntersection(sbrit,tmp) ewowin <- as.owin(sbrit) ## End(Not run)
Rescale a numeric vector to a specified minimum and maximum.
gaussW(x, mean, sd, type = "weights")gaussW(x, mean, sd, type = "weights")
x |
A numeric vector or an object of class CalDates. |
mean |
A single numeric value indicating the value to centre the Gaussian kernel on. |
sd |
A single numeric value indicating the standard deviation of the Gaussian kernel to be used. |
type |
The type of output to produce: currently either "weighted" (for a simple total weight value for each date) or "raw" (a list of reweighted calibrated radiocarbon probabilities for each calibrated date). |
A numeric vector of weights (or optionally a list of reweighted calibrated radiocarbon probabilities).
## Example weighting fo a set of dates versus a focal date of 5950 calBP years <- seq(6500, 5500, -10) plot(cbind(years, gaussW(years, 5950, 50))) ## Example weighting of three calibrated dates versus a focal date of 5950 calBP dates <- calibrate(c(5280, 5180, 5080), c(30,30,30), normalised=FALSE) gaussW(dates, 5950, 50) ## Or the same with raw output dateswt <- gaussW(dates, 5950, 50, type="raw") head(dateswt[[1]])## Example weighting fo a set of dates versus a focal date of 5950 calBP years <- seq(6500, 5500, -10) plot(cbind(years, gaussW(years, 5950, 50))) ## Example weighting of three calibrated dates versus a focal date of 5950 calBP dates <- calibrate(c(5280, 5180, 5080), c(30,30,30), normalised=FALSE) gaussW(dates, 5950, 50) ## Or the same with raw output dateswt <- gaussW(dates, 5950, 50, type="raw") head(dateswt[[1]])
Computes the highest posterior density interval (HPDI) for calibrated dates
hpdi( x, credMass = 0.95, asList = TRUE, calendar = "BP", sep = "|", pdigits = 3 )hpdi( x, credMass = 0.95, asList = TRUE, calendar = "BP", sep = "|", pdigits = 3 )
x |
A |
credMass |
Interval probability mass |
asList |
Logical variable determining whether the output should be a list (TRUE) or a character vector (FALSE). Default is TRUE |
calendar |
Whether daes should be reported in cal BP ( |
sep |
character string to separate date ranges. Default is '|'. Ignored when |
pdigits |
indicate the number of decimal places for reporting probabilities. Default is 3. Ignored when |
A list of matrices with HPDI ranges and associated probabilities in cal BP. Note that the sum of the probability mass will be only approximately equal to the argument credMass.
x <- calibrate(c(2000,3050,2950),c(35,20,20)) hpdi(x) hpdi(x,asList=FALSE,calendar='BCAD')x <- calibrate(c(2000,3050,2950),c(35,20,20)) hpdi(x) hpdi(x,asList=FALSE,calendar='BCAD')
Function for generating a vector median calibrated dates from a CalDates class object.
medCal(x)medCal(x)
x |
A |
A vector of median dates in cal BP
x <- calibrate(c(3050,2950),c(20,20)) medCal(x)x <- calibrate(c(3050,2950),c(20,20)) medCal(x)
Function for generating mixed calibration curves (e.g. intcal20/marine20, intcal20/shcal20) with user-defined proportions.
mixCurves( calCurve1 = "intcal20", calCurve2 = "marine20", p = 0.5, resOffsets = 0, resErrors = 0 )mixCurves( calCurve1 = "intcal20", calCurve2 = "marine20", p = 0.5, resOffsets = 0, resErrors = 0 )
calCurve1 |
Name of the first curve, chosen from 'intcal20','intcal13','intcal13nhpine16','shcal20','shcal13','shcal13shkauri16',”marine13','marine20'. Default is 'intcal20'. |
calCurve2 |
Name of the second curve, from the same list. Default is 'marine20'. |
p |
The proportion of contribution of the first curve. Default is 1. |
resOffsets |
Offset value for the marine reservoir effect to be applied to calCurve2. Default is 0. |
resErrors |
Error of the marine reservoir effect offset to be applied to calCurve2. Default is 0. |
The function is based on the mix.calibrationcurves function of the clam package.
A three-column matrix containing calibrated year BP, uncalibrated age bp, and standard deviation. To be used as custom calibration curve for the calibrate function.
Blaauw, M. and Christen, J.A.. 2011. Flexible paleoclimate age-depth models using an autoregressive gamma process. Bayesian Analysis, 6, 457-474. Blaauw, M. 2018. clam: Classical Age-Depth Modelling of Cores from Deposits. R package version 2.3.1. https://CRAN.R-project.org/packacge=clam Marsh, E.J., Bruno, M.C., Fritz, S.C., Baker, P., Capriles, J.M. and Hastorf, C.A., 2018. IntCal, SHCal, or a Mixed Curve? Choosing a 14 C Calibration Curve for Archaeological and Paleoenvironmental Records from Tropical South America. Radiocarbon, 60(3), pp.925-940.
myCurve <- mixCurves('intcal20','marine20',p=0.7,resOffsets=300,resErrors=20) x <- calibrate(4000,30,calCurves=myCurve)myCurve <- mixCurves('intcal20','marine20',p=0.7,resOffsets=300,resErrors=20) x <- calibrate(4000,30,calCurves=myCurve)
Comparison of an observed summed radiocarbon date distribution (aka SPD) with simulated outcomes from a theoretical model.
modelTest( x, errors, nsim, bins = NA, runm = NA, timeRange = NA, backsight = 50, changexpr = expression((t1/t0)^(1/d) - 1), gridclip = TRUE, raw = FALSE, model = c("exponential"), method = c("uncalsample"), predgrid = NA, normalised = NA, datenormalised = NA, spdnormalised = FALSE, ncores = 1, fitonly = FALSE, a = 0, b = 0, edgeSize = 500, verbose = TRUE )modelTest( x, errors, nsim, bins = NA, runm = NA, timeRange = NA, backsight = 50, changexpr = expression((t1/t0)^(1/d) - 1), gridclip = TRUE, raw = FALSE, model = c("exponential"), method = c("uncalsample"), predgrid = NA, normalised = NA, datenormalised = NA, spdnormalised = FALSE, ncores = 1, fitonly = FALSE, a = 0, b = 0, edgeSize = 500, verbose = TRUE )
x |
A |
errors |
A vector of errors corresponding to each radiocarbon age |
nsim |
Number of simulations |
bins |
A vector indicating which bin each radiocarbon date is assigned to. |
runm |
A number indicating the window size of the moving average to smooth both observed and simulated SPDs. If set to |
timeRange |
A vector of length 2 indicating the start and end date of the analysis in cal BP. The fitting process is applied considering the SPD within the interval defined by this parameter. If no values are supplied the earliest and latest median calibrated dates of the observed data will be used. |
backsight |
A single numeric value defining the distance in time between the focal year and the backsight year for computing the rate of change. Default is 10. |
changexpr |
An expression for calculating the rate of change in SPD between the focal year and a backsight year. Available input options are t1 (the SPD for the focal year), t0 (the SPD for the backsight year), d (the distance between t0 and t1), and any other standard constants and mathematical operators. A sensible default is provided. |
gridclip |
Whether the sampling of random dates is constrained to the observed range (TRUE) or not (FALSE). Default is TRUE. |
raw |
A logical variable indicating whether all permuted SPDs should be returned or not. Default is FALSE. |
model |
A vector indicating the model to be fitted. Currently the acceptable options are |
method |
Method for the creation of random dates from the fitted model. Either |
predgrid |
A data.frame containing calendar years (column |
normalised |
Whether the simulated dates should be normalised or not. Default based on whether x is normalised or not. |
datenormalised |
Argument kept for backward compatibility with previous versions. |
spdnormalised |
A logical variable indicating whether the total probability mass of the SPD is normalised to sum to unity for both observed and simulated data. |
ncores |
Number of cores used for for parallel execution. Default is 1. |
fitonly |
A logical variable. If set to TRUE, only the the model fitting is executed and returned. Default is FALSE. |
a |
Starter value for the exponential fit with the |
b |
Starter value for the exponential fit with the |
edgeSize |
Controls edge effect by expanding the fitted model beyond the range defined by |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
The function implements a Monte-Carlo test for comparing a theoretical or fitted statistical model to an observed summed radiocarbon date distribution (aka SPD) and associated rates of changes. A variety of theoretical expectations can be compared to the observed distribution by setting the model argument, for example to fit basic 'uniform' (the mean of the SPD), 'linear' (fitted using the lm function) or model='exponential' models (fitted using the nls function). Models are fitted to the period spanned by timeRange although x can contain dates outside this range to mitigate possible edge effects (see also bracket). Notice that estimating growth model parameters directly from SPD can be biased, resulting into a null hypothesis that is not necessarily representative (see Crema 2022). It is also possible for the user to provide a model of their own by setting model='custom' and then supplying a two-column data.frame to predgrid. The function generates nsim theoretical SPDs from the fitted model via Monte-Carlo simulation, this is then used to define a 95% critical envelope for each calendar year. The observed SPD is then compared against the simulation envelope; local departures from the model are defined as instances where the observed SPD is outside such an envelope, while an estimate of the global significance of the observed SPD is also computed by comparing the total areas of observed and simulated SPDs that fall outside the simulation envelope. The theoretical SPDs can be generated using two different sampling approaches defined by the parameter method. If method is set to 'uncalsample' each date is drawn after the fitted model is backcalibrated as a whole and adjusted for a baseline expectation; if it is set to 'calsample' samples are drawn from the fitted model in calendar year then individually back calibrated and recalibrated (the approach of Timpson et al. 2014). For each simulation, both approaches produces samples, with equal to the number of bins or number of dates (when bins are not defined). Differences between these two approaches are particularly evident at dates coincident with steeper portions of the calibration curve. If more than one type of calibration curve is associated with the observed dates, at each Monte-Carlo iteration, the function randomly assigns each bin to one of the calibration curves with probability based on the proportion of dates within the bin associated to the specific curves. For example, if a bin is composed of four dates and three are calibrated with 'intcal20' the probability of that particular bin being assigned to 'intcal20' is 0.75.
An object of class SpdModelTest with the following elements
result A four column data.frame containing the observed probability density (column PrDens) and the lower and the upper values of the simulation envelope (columns lo and hi) for each calendar year column calBP
result.roc A four column data.frame containing the observed rates of change (column roc) and the lower and the upper values of the simulation envelope (columns lo.roc and hi.roc) for the mid point between two chronological blocks calBP
sim A matrix containing the simulation results of the summed probabilities. Available only when raw is set to TRUE
sim.roc A matrix containing the simulation results of the rate of change of summed probabilities. Available only when raw is set to TRUE
pval A numeric vector containing the p-value of the global significance test for the summed probabilities
pval.roc A numeric vector containing the p-value of the global significance test for the rates of change
fit A data.frame containing the probability densities of the fitted model for each calendar year within the time range of analysis
fitobject Fitted model. Not available when model is 'custom'
n Number of radiocarbon dates.
nbinsNumber of bins.
nsimNumber of Monte-Carlo simulations.
backsightBacksight size.
Windows users might receive a memory allocation error with larger time span of analysis (defined by the parameter timeRange). This can be avoided by increasing the memory limit with the memory.limit function.
Users experiencing a Error: cannot allocate vector of size ... error message can increase the memory size using the Sys.setenv(), for example: Sys.setenv("R_MAX_VSIZE" = 16e9).
The function currently supports only dates calibrated with 'intcal20',intcal13','intcal13nhpine16','shcal20','shcal13','shcal13shkauri16', 'marine20', and 'marine13'.
Crema, E.R. (2022). Statistical inference of prehistoric demography from frequency distributions of radiocarbon dates: a review and a guide for the perplexed. Journal of Archaeological Method and Theory. doi:10.1007/s10816-022-09559-5
Timpson, A., Colledge, S., Crema, E., Edinborough, K., Kerig, T., Manning, K., Thomas, M.G., Shennan, S., (2014). Reconstructing regional population fluctuations in the European Neolithic using radiocarbon dates: a new case-study using an improved method. Journal of Archaeological Science, 52, 549-557. doi:10.1016/j.jas.2014.08.011
## Example with Younger Dryas period Near East, including site bins ## Not run: data(emedyd) caldates <- calibrate(x=emedyd$CRA, errors=emedyd$Error, normalised=FALSE) bins <- binPrep(sites=emedyd$SiteName, ages=emedyd$CRA, h=50) nsim=5 #toy example expnull <- modelTest(caldates, errors=emedyd$Error, bins=bins, nsim=nsim, runm=50, timeRange=c(16000,9000), model="exponential", datenormalised=FALSE) plot(expnull, xlim=c(16000,9000)) round(expnull$pval,4) #p-value summary(expnull) ## End(Not run)## Example with Younger Dryas period Near East, including site bins ## Not run: data(emedyd) caldates <- calibrate(x=emedyd$CRA, errors=emedyd$Error, normalised=FALSE) bins <- binPrep(sites=emedyd$SiteName, ages=emedyd$CRA, h=50) nsim=5 #toy example expnull <- modelTest(caldates, errors=emedyd$Error, bins=bins, nsim=nsim, runm=50, timeRange=c(16000,9000), model="exponential", datenormalised=FALSE) plot(expnull, xlim=c(16000,9000)) round(expnull$pval,4) #p-value summary(expnull) ## End(Not run)
Plot multiple radiocarbon dates.
multiplot( x, type = "d", calendar = "BP", HPD = FALSE, credMass = 0.95, decreasing = NULL, label = TRUE, label.pos = 0.5, label.offset = 0, xlim = NULL, xlab = NA, ylab = NA, col.fill = "grey50", col.fill2 = "grey82", col.line = "black", lwd = 1, cex.id = 1, cex.lab = 1, cex.axis = 1, ydisp = FALSE, gapFactor = 0.2, rescale = FALSE )multiplot( x, type = "d", calendar = "BP", HPD = FALSE, credMass = 0.95, decreasing = NULL, label = TRUE, label.pos = 0.5, label.offset = 0, xlim = NULL, xlab = NA, ylab = NA, col.fill = "grey50", col.fill2 = "grey82", col.line = "black", lwd = 1, cex.id = 1, cex.lab = 1, cex.axis = 1, ydisp = FALSE, gapFactor = 0.2, rescale = FALSE )
x |
A CalDates class object with length > 1. |
type |
Whether the calibrated dates are displayed as distributions ( |
calendar |
Either |
HPD |
Logical value indicating whether intervals of higher posterior density should be displayed. Default is FALSE. |
credMass |
A numerical value indicating the size of the higher posterior density interval. Default is 0.95. |
decreasing |
Whether dates should be plotted with decreasing order of median calibrated date (i.e. old to new; TRUE) or increasing order (i.e. new to old; FALSE). If set to NULL the dates plotted in the supplied order. Default is NULL |
label |
Whether the ID of each date should be displayed. Default is TRUE. |
label.pos |
Relative position of the label in relation to the calibrated distribution expressed in quantiles. Default is 0.5 (median). |
label.offset |
Horrizontal offset of label position in number of years. Default is 0. |
xlim |
the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter |
xlab |
(optional) Label for the x axis. If unspecified the default setting will be applied ("Year BP" or "Year BC/AD"). |
ylab |
(optional) Label for the y axis. |
col.fill |
A vector of primary fill color for the calibrated date distribution. Default is 'grey50'. |
col.fill2 |
A vector of secondary secondary colour fill color for the calibrated date distribution, used for regions outside the higher posterior interval. Ignored when |
col.line |
A vector of line color (ignored when |
lwd |
Line width (ignored when |
cex.id |
The magnification to be used the date labels relative to the current setting of cex. Default is adjusted to 1. |
cex.lab |
The magnification to be used for x and y labels relative to the current setting of cex. Default is adjusted to 1. |
cex.axis |
The magnification to be used for axis annotation relative to the current setting of cex. Default is adjusted to 1. |
ydisp |
Whether the y axis should be displayed. Ignored when |
gapFactor |
Defines spacing between calibrated dates (when |
rescale |
Whether probability distributions should be rescaled (applicable only when |
data("emedyd") tellaswad = subset(emedyd,SiteName=='Tell Aswad') x = calibrate(tellaswad$CRA,tellaswad$Error,ids=tellaswad$LabID) multiplot(x,HPD=TRUE,decreasing=TRUE,label=FALSE,gapFactor = 0.1) multiplot(x,type='b',calendar='BCAD',cex.id = 0.5,lwd=2,gapFactor = 0.5)data("emedyd") tellaswad = subset(emedyd,SiteName=='Tell Aswad') x = calibrate(tellaswad$CRA,tellaswad$Error,ids=tellaswad$LabID) multiplot(x,HPD=TRUE,decreasing=TRUE,label=FALSE,gapFactor = 0.1) multiplot(x,type='b',calendar='BCAD',cex.id = 0.5,lwd=2,gapFactor = 0.5)
Test for evaluating the difference in the summed probability values associated with two points in time.
p2pTest(x, p1 = NA, p2 = NA, interactive = TRUE, plot = FALSE)p2pTest(x, p1 = NA, p2 = NA, interactive = TRUE, plot = FALSE)
x |
result of |
p1 |
calendar year (in BP) of start point. |
p2 |
calendar year (in BP) of end point. |
interactive |
if set to TRUE enables an interactive selection of p1 and p2 from a graphical display of the SPD. Disabled when |
plot |
if set to TRUE the function plots the location of p1 and p2 on the SPD. Default is FALSE. |
The function compares observed differences in the summed probability values associated with two points in time against a distribution of expected values under the null hypothesis defined with the modelTest function. The two points can be specified manually (assigning BP dates to the arguments p1 and p2) or interactively (clicking on a SPD plot). Note that modelTest should be executed setting the argument raw to TRUE (default is FALSE).
A list with: the BP dates for the two points and the p-value obtained from a two-sided test.
Edinborough, K., Porcic, M., Martindale, A., Brown, T.J., Supernant, K., Ames, K.M., (2017). Radiocarbon test for demographic events in written and oral history. PNAS 201713012. doi:10.1073/pnas.1713012114
## Example with Younger Dryas period Near East, including site bins ## Not run: data(emedyd) caldates <- calibrate(x=emedyd$CRA, errors=emedyd$Error, normalised=FALSE) bins <- binPrep(sites=emedyd$SiteName, ages=emedyd$CRA, h=50) nsim=10 #toy example expnull <- modelTest(caldates, errors=emedyd$Error, bins=bins, nsim=nsim, runm=50, timeRange=c(16000,9000), model="exponential", datenormalised=FALSE, raw=TRUE) p2pTest(x=expnull,p1=13000,p2=12500) #non-interactive mode p2pTest(x=expnull) #interactive mode ## End(Not run)## Example with Younger Dryas period Near East, including site bins ## Not run: data(emedyd) caldates <- calibrate(x=emedyd$CRA, errors=emedyd$Error, normalised=FALSE) bins <- binPrep(sites=emedyd$SiteName, ages=emedyd$CRA, h=50) nsim=10 #toy example expnull <- modelTest(caldates, errors=emedyd$Error, bins=bins, nsim=nsim, runm=50, timeRange=c(16000,9000), model="exponential", datenormalised=FALSE, raw=TRUE) p2pTest(x=expnull,p1=13000,p2=12500) #non-interactive mode p2pTest(x=expnull) #interactive mode ## End(Not run)
Global and local significance test for comparing shapes of multiple SPDs using random permutations.
permTest( x, marks, timeRange, backsight = 10, changexpr = expression((t1/t0)^(1/d) - 1), nsim, bins = NA, runm = NA, datenormalised = FALSE, spdnormalised = FALSE, raw = FALSE, verbose = TRUE )permTest( x, marks, timeRange, backsight = 10, changexpr = expression((t1/t0)^(1/d) - 1), nsim, bins = NA, runm = NA, datenormalised = FALSE, spdnormalised = FALSE, raw = FALSE, verbose = TRUE )
x |
A |
marks |
A numerical or character vector containing the marks associated to each radiocarbon date. |
timeRange |
A vector of length 2 indicating the start and end date of the analysis in cal BP. |
backsight |
A single numeric value defining the distance in time between the focal year and the backsight year for computing the rate of change. Default is 10. |
changexpr |
An expression for calculating the rate of change in SPD between the focal year and a backsight year. Available input options are t1 (the SPD for the focal year), t0 (the SPD for the backsight year), d (the distance between t0 and t1), and any other standard constants and mathematical operators. A sensible default is provided. |
nsim |
Number of random permutations |
bins |
A vector indicating which bin each radiocarbon date is assigned to. |
runm |
A number indicating the window size of the moving average to smooth the SPD. If set to |
datenormalised |
If set to TRUE the total probability mass of each calibrated date will be made to sum to unity (the default in most radiocarbon calibration software). This argument will only have an effect if the dates in |
spdnormalised |
A logical variable indicating whether the total probability mass of the SPD is normalised to sum to unity. |
raw |
A logical variable indicating whether all permuted SPDs should be returned or not. Default is FALSE. |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
The function generates a distribution of expected SPDs by randomly shuffling the marks assigned to each bin (see spd for details on binning). The resulting distribution of probabilities for each mark (i.e. group of dates) for each calendar year is z-transformed, and a 95% simulation envelope is computed. Local significant departures are defined as instances where the observed SPD (which is also z-transformed) is outside such envelope. A global significance is also computed by comparing the total "area" outside the simulation envelope in the observed and simulated data.
An object of class SpdPermTest with the following elements
observed A list containing data.frames with the summed probability (column PrDens for each calendar year (column calBP for each mark/group
envelope A list containing matrices with the lower and upper bound values of the simulation envelope for each mark/group
pValueList A list of p-value associated with each mark/group
Crema, E.R., Habu, J., Kobayashi, K., Madella, M., (2016). Summed Probability Distribution of 14 C Dates Suggests Regional Divergences in the Population Dynamics of the Jomon Period in Eastern Japan. PLOS ONE 11, e0154809. doi:10.1371/journal.pone.0154809
## compare demographic trajectories in Netherlands and Denmark ## Not run: data(euroevol) nld.dnk = subset(euroevol,Country=="Netherlands"|Country=="Denmark") bins = binPrep(nld.dnk$SiteID,nld.dnk$C14Age,h=200) dates = calibrate(nld.dnk$C14Age,nld.dnk$C14SD,normalised=FALSE) res = permTest(dates,marks=as.character(nld.dnk$Country),nsim=1000, bins=bins,runm=200,timeRange=c(10000,4000)) round(res$pValueList,4) #extract p-values summary(res) par(mfrow=c(2,1)) plot(res,focalm="Netherlands",main="Netherlands") plot(res,focalm="Denmark",main="Denmark") ## End(Not run)## compare demographic trajectories in Netherlands and Denmark ## Not run: data(euroevol) nld.dnk = subset(euroevol,Country=="Netherlands"|Country=="Denmark") bins = binPrep(nld.dnk$SiteID,nld.dnk$C14Age,h=200) dates = calibrate(nld.dnk$C14Age,nld.dnk$C14SD,normalised=FALSE) res = permTest(dates,marks=as.character(nld.dnk$Country),nsim=1000, bins=bins,runm=200,timeRange=c(10000,4000)) round(res$pValueList,4) #extract p-values summary(res) par(mfrow=c(2,1)) plot(res,focalm="Netherlands",main="Netherlands") plot(res,focalm="Denmark",main="Denmark") ## End(Not run)
Plot calibrated radiocarbon dates.
## S3 method for class 'CalDates' plot( x, ind = 1, label = NA, calendar = "BP", type = "standard", xlab = NA, ylab = NA, axis4 = TRUE, HPD = FALSE, credMass = 0.95, customCalCurve = NA, add = FALSE, col = "grey50", col2 = "grey82", cex.axis = 0.75, cex.xylab = 0.75, cex.label = 0.75, ... )## S3 method for class 'CalDates' plot( x, ind = 1, label = NA, calendar = "BP", type = "standard", xlab = NA, ylab = NA, axis4 = TRUE, HPD = FALSE, credMass = 0.95, customCalCurve = NA, add = FALSE, col = "grey50", col2 = "grey82", cex.axis = 0.75, cex.xylab = 0.75, cex.label = 0.75, ... )
x |
|
ind |
Number indicating the index value of the calibrated radiocarbon date to be displayed. Default is 1. |
label |
(optional) Character vector to be shown on the top-right corner of the display window. |
calendar |
Either |
type |
Either |
xlab |
(optional) Label for the x axis. If unspecified the default setting will be applied ("Year BP" or "Year BC/AD"). |
ylab |
(optional) Label for the y axis. If unspecified the default setting will be applied ("Radiocarbon Age"). |
axis4 |
Logical value indicating whether an axis of probabilities values should be displayed. Default is TRUE. |
HPD |
Logical value indicating whether intervals of higher posterior density should be displayed. Default is FALSE. |
credMass |
A numerical value indicating the size of the higher posterior density interval. Default is 0.95. |
customCalCurve |
A three column data.frame or matrix that allows you to pass and plot a custom calibration curve if you used one during calibration. You can currently only provide one such custom curve which is used for all dates. |
add |
if set to |
col |
The primary fill color for the calibrated date distribution. |
col2 |
The secondary colour fill color for the calibrated date distribution, used for regions outside the higher posterior interval. Ignored when |
cex.axis |
The magnification to be used for axis annotation relative to the current setting of cex. Default is adjusted to 0.75. |
cex.xylab |
The magnification to be used for x and y labels relative to the current setting of cex. Default is adjusted to 0.75. |
cex.label |
The magnification to be used for the label on the top right corner defined by the argument |
... |
Additional arguments affecting the plot. |
x <- calibrate(x=c(3402,3490,4042),errors=c(20,20,30)) plot(x) #display the first date plot(x,2) #displays the second date plot(x,3, calendar="BCAD", HPD=TRUE) #display in BC/AD with higher posterior density intervalx <- calibrate(x=c(3402,3490,4042),errors=c(20,20,30)) plot(x) #display the first date plot(x,2) #displays the second date plot(x,3, calendar="BCAD", HPD=TRUE) #display in BC/AD with higher posterior density interval
Plot a summed radiocarbon probability distribution. This is a basic function for plotting SPDs that have been constructed manually or by calibrating a summed or otherwise irregular CRA grid. In most instances, it is sensible to use plot.CalSPD instead.
## S3 method for class 'CalGrid' plot( x, runm = NA, calendar = "BP", fill.p = "grey50", border.p = NA, xlim = NA, ylim = NA, cex.lab = 0.75, cex.axis = cex.lab, mar = c(4, 4, 1, 3), add = FALSE, ... )## S3 method for class 'CalGrid' plot( x, runm = NA, calendar = "BP", fill.p = "grey50", border.p = NA, xlim = NA, ylim = NA, cex.lab = 0.75, cex.axis = cex.lab, mar = c(4, 4, 1, 3), add = FALSE, ... )
x |
A "CalGrid" class object of summed probabilities per calendar year BP. |
runm |
A number indicating the window size of the moving average to smooth the SPD. If set to |
calendar |
Either |
fill.p |
Fill colour of the polygon depicting the summed probability distribution. |
border.p |
Border colour of the polygon depicting the summed probability distribution. |
xlim |
the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter |
ylim |
Adjust y axis limits (otherwise sensible default). |
cex.lab |
Size of label text. |
cex.axis |
Size of axis text. |
mar |
Adjust margins around plot. |
add |
Whether or not the new graphic should be added to an existing plot. |
... |
Additional arguments affecting the plot |
data(euroevol) mycaldates <- calibrate(euroevol[1:10,"C14Age"], euroevol[1:10,"C14SD"], normalised=FALSE) myspd <- spd(mycaldates, timeRange=c(8000,2000)) plot(myspd) #ordinary plot using \code{plot.CalSPD} plot(myspd$grid) #working plot using the internal CalGrid objectdata(euroevol) mycaldates <- calibrate(euroevol[1:10,"C14Age"], euroevol[1:10,"C14SD"], normalised=FALSE) myspd <- spd(mycaldates, timeRange=c(8000,2000)) plot(myspd) #ordinary plot using \code{plot.CalSPD} plot(myspd$grid) #working plot using the internal CalGrid object
Plot a summed probability distribution (SPD) of radiocarbon dates
## S3 method for class 'CalSPD' plot( x, runm = NA, calendar = "BP", type = "standard", xlim = NA, ylim = NA, ylab = "Summed Probability", spdnormalised = FALSE, rescale = FALSE, fill.p = "grey75", border.p = NA, xaxt = "s", yaxt = "s", add = FALSE, cex.axis = 1, ... )## S3 method for class 'CalSPD' plot( x, runm = NA, calendar = "BP", type = "standard", xlim = NA, ylim = NA, ylab = "Summed Probability", spdnormalised = FALSE, rescale = FALSE, fill.p = "grey75", border.p = NA, xaxt = "s", yaxt = "s", add = FALSE, cex.axis = 1, ... )
x |
A |
runm |
A number indicating the window size of the moving average to smooth the SPD. If set to |
calendar |
Either |
type |
Either |
xlim |
the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter |
ylim |
the y limits of the plot. |
ylab |
(optional) Label for the y axis. If unspecified the default setting will be applied ("Summed Probability") |
spdnormalised |
A logical variable indicating whether the total probability mass of the SPD is normalised to sum to unity. |
rescale |
A logical variable indicating whether the SPD should be rescaled to range 0 to 1. |
fill.p |
Fill colour for the SPD |
border.p |
Border colour for the SPD |
xaxt |
Whether the x-axis tick marks should be displayed ( |
yaxt |
Whether the y-axis tick marks should be displayed ( |
add |
Whether or not the new graphic should be added to an existing plot. |
cex.axis |
The magnification to be used for axis annotation relative to the current setting of cex. Default is 1. |
... |
Additional arguments affecting the plot |
## Not run: data(emedyd) levant <- emedyd[emedyd$Region=="1"|emedyd$Region=="2",] bins <- binPrep(levant$SiteName, levant$CRA, h=50) x <- calibrate(levant$CRA, levant$Error, normalised=FALSE) spd.levant <- spd(x, bins=bins, timeRange=c(17000,8000)) spd.northernlevant <- spd(x[levant$Region=="2"], bins=bins[levant$Region=="2"], timeRange=c(17000,8000)) plot(spd.levant, runm=50, xlim=c(16000,9000)) plot(spd.northernlevant, runm=50, add=TRUE, fill.p="black") legend("topleft", legend=c("All Levant dates","Northern Levant only"), fill=c("grey75","black"), border=NA) plot(spd.levant, runm=50, xlim=c(16000,9000), type="simple") plot(spd.northernlevant, runm=50, col="red", type="simple", add=TRUE) ## End(Not run)## Not run: data(emedyd) levant <- emedyd[emedyd$Region=="1"|emedyd$Region=="2",] bins <- binPrep(levant$SiteName, levant$CRA, h=50) x <- calibrate(levant$CRA, levant$Error, normalised=FALSE) spd.levant <- spd(x, bins=bins, timeRange=c(17000,8000)) spd.northernlevant <- spd(x[levant$Region=="2"], bins=bins[levant$Region=="2"], timeRange=c(17000,8000)) plot(spd.levant, runm=50, xlim=c(16000,9000)) plot(spd.northernlevant, runm=50, add=TRUE, fill.p="black") legend("topleft", legend=c("All Levant dates","Northern Levant only"), fill=c("grey75","black"), border=NA) plot(spd.levant, runm=50, xlim=c(16000,9000), type="simple") plot(spd.northernlevant, runm=50, col="red", type="simple", add=TRUE) ## End(Not run)
Plots a Composite Kernel Density Estimate of sampled radiocarbon dates.
## S3 method for class 'compositeKDE' plot( x, calendar = "BP", type = "envelope", ylim = NA, xlim = NA, fill.col = "lightgrey", interval = 0.95, line.col = "black", line.type = 2, multiline.alpha = NA, multiline.col = "black", ... )## S3 method for class 'compositeKDE' plot( x, calendar = "BP", type = "envelope", ylim = NA, xlim = NA, fill.col = "lightgrey", interval = 0.95, line.col = "black", line.type = 2, multiline.alpha = NA, multiline.col = "black", ... )
x |
A |
calendar |
Either |
type |
Either |
ylim |
the y limits of the plot. |
xlim |
the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter |
fill.col |
Envelope color when |
interval |
Quantile interval for the envelope. Default is 0.95. |
line.col |
Line color when |
line.type |
Line type when |
multiline.alpha |
Alpha level for line transparency when |
multiline.col |
Line color when |
... |
Additional arguments affecting the plot |
Visualise a compositeKDE class object. If type is set 'envelope' an envelope of the percentile interval defined by the parameter interval is shown along with the mean KDE. If type is set 'multiline' all KDEs are shown.
ckde;
Displays local growth rates, p-values, and q-values retrieved from a spatialTest class object.
## S3 method for class 'spatialTest' plot( x, index = 1, option, breakRange, breakLength = 7, rd = 5, baseSize = 0.5, plim = 0.05, qlim = 0.05, legend = FALSE, legSize = 1, location = "bottomright", ... )## S3 method for class 'spatialTest' plot( x, index = 1, option, breakRange, breakLength = 7, rd = 5, baseSize = 0.5, plim = 0.05, qlim = 0.05, legend = FALSE, legSize = 1, location = "bottomright", ... )
x |
A |
index |
A numerical value indicating which transition to display. Ignored when |
option |
Either " |
breakRange |
A vector of length 2 defining the minimum and maximum values of growth rate to be displayed in the legend. If set to NA its computed from data range (default). |
breakLength |
A numerical vector defining the number of breaks for growth rates to be displayed in the legend. |
rd |
Number of decimal places of the growth rate to be displayed in the Legend |
baseSize |
Numerical value giving the amount by which points should be magnified relative to the default settings in R. Default is 0.5 |
plim |
Threshold value for the p-values. Default is 0.05. |
qlim |
Threshold value for the q-values. Default is 0.05. |
legend |
Logical values specifying whether the legend should be displayed or not. Default is FALSE. |
legSize |
Numerical value giving the amount by which points should be magnified relative to the default settings in R for the Legend. Default is 1. |
location |
A single keyword from the list "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center" to specify the location of the Legend. Default is "bottomright". |
... |
Graphical parameters to be passed to methods. |
The function displays a distribution map of local growth rates (when option="raw"), q- and p-values (when option="test"), and the associated legends (when option="rawlegend" or option="testlegend").
The function visualises the observed summed probability distribution of radiocarbon dates along with a simulation envelope for the null model and regions of positive and negative deviation.
## S3 method for class 'SpdModelTest' plot( x, calendar = "BP", type = "spd", ylim = NA, xlim = NA, col.obs = "black", lwd.obs = 0.5, xaxs = "i", yaxs = "i", bbty = "f", bbtyRet = TRUE, drawaxes = TRUE, ... )## S3 method for class 'SpdModelTest' plot( x, calendar = "BP", type = "spd", ylim = NA, xlim = NA, col.obs = "black", lwd.obs = 0.5, xaxs = "i", yaxs = "i", bbty = "f", bbtyRet = TRUE, drawaxes = TRUE, ... )
x |
A |
calendar |
Either |
type |
Whether to display SPDs ('spd') or rates of change ('roc'). Default is 'spd'. |
ylim |
the y limits of the plot. |
xlim |
the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter |
col.obs |
Line colour for the observed SPD. |
lwd.obs |
Line width for the observed SPD. |
xaxs |
The style of x-axis interval calculation (see |
yaxs |
The style of y-axis interval calculation (see |
bbty |
Display options; one between |
bbtyRet |
Whether details on the intervals of positive and negative deviations are returned. Default is TRUE. |
drawaxes |
A logical value determining whether the axes should be displayed or not. Default is TRUE. |
... |
Additional arguments affecting the plot |
The argument bbty controls the display options of the Monte-Carlo Test. Default settings (bbty='f') displays the observed SPD (solid black line), the simulation envelope of the fitted model (shaded grey polygon) and regions of significance positive (red semi-transparent rectangle) and negative (blue semi-transparent rectangle) deviation. The option bbty='b' removes the regions of positive/negative deviations, whilst the option bbty='n' displays the simulation envelope on existing plot.
Visualises the observed SPD along with the simulation envelope generated from permTest, with regions of positive and negative deviations highlighted in red and blue.
## S3 method for class 'SpdPermTest' plot( x, focalm = 1, type = "spd", calendar = "BP", xlim = NA, ylim = NA, col.obs = "black", col.env = rgb(0, 0, 0, 0.2), lwd.obs = 0.5, xaxs = "i", yaxs = "i", bbty = "f", drawaxes = TRUE, ... )## S3 method for class 'SpdPermTest' plot( x, focalm = 1, type = "spd", calendar = "BP", xlim = NA, ylim = NA, col.obs = "black", col.env = rgb(0, 0, 0, 0.2), lwd.obs = 0.5, xaxs = "i", yaxs = "i", bbty = "f", drawaxes = TRUE, ... )
x |
A |
focalm |
Value specifying the name of the focal mark (group) to be plotted. |
type |
Whether to display SPDs ('spd') or rates of change ('roc'). Default is 'spd'. |
calendar |
Either |
xlim |
the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter |
ylim |
the y limits of the plot. |
col.obs |
Line colour for the observed SPD |
col.env |
Colour for the simulation envelope |
lwd.obs |
Line width for the observed SPD |
xaxs |
The style of x-axis interval calculation (see |
yaxs |
The style of y-axis interval calculation (see |
bbty |
Display options; one between |
drawaxes |
A logical value determining whether the axes should be displayed or not. Default is TRUE. |
... |
Additional arguments affecting the plot |
spdRC class objectsPlot rates of change between time-blocks
## S3 method for class 'spdRC' plot( x, calendar = "BP", col.obs = "black", lwd.obs = 0.5, xaxs = "i", yaxs = "i", xlim = NA, ... )## S3 method for class 'spdRC' plot( x, calendar = "BP", col.obs = "black", lwd.obs = 0.5, xaxs = "i", yaxs = "i", xlim = NA, ... )
x |
|
calendar |
Either |
col.obs |
Line colour for the observed SPD |
lwd.obs |
Line width for the observed SPD |
xaxs |
The style of x-axis interval calculation (see |
yaxs |
The style of y-axis interval calculation (see |
xlim |
Limits for the x axis. |
... |
Additional arguments affecting the plot. |
Visualises multiple SPDs grouped as a stackCalSPD object.
## S3 method for class 'stackCalSPD' plot( x, type = "stacked", calendar = "BP", spdnormalised = FALSE, rescale = FALSE, runm = NA, xlim = NA, ylim = NA, xaxt = "s", yaxt = "s", gapFactor = 0.2, col.fill = NA, col.line = NA, lwd.obs = 1, lty.obs = 1, cex.lab = 1, cex.axis = cex.lab, legend = TRUE, legend.arg = NULL, ylab = NA, ymargin = 1.1, rnd = 2, ... )## S3 method for class 'stackCalSPD' plot( x, type = "stacked", calendar = "BP", spdnormalised = FALSE, rescale = FALSE, runm = NA, xlim = NA, ylim = NA, xaxt = "s", yaxt = "s", gapFactor = 0.2, col.fill = NA, col.line = NA, lwd.obs = 1, lty.obs = 1, cex.lab = 1, cex.axis = cex.lab, legend = TRUE, legend.arg = NULL, ylab = NA, ymargin = 1.1, rnd = 2, ... )
x |
A |
type |
How to display the SPDs.Current options are |
calendar |
Either |
spdnormalised |
A logical variable indicating whether the total probability mass of the SPDs are normalised to sum to unity. Default is FALSE. |
rescale |
A logical variable indicating whether the summed probabilities values should be rescaled to range 0 to 1. Default is FALSE.Notice that this is different from setting |
runm |
A number indicating the window size of the moving average to smooth the SPD. If set to |
xlim |
the x limits of the plot. In BP or in BC/AD depending on the choice of the parameter |
ylim |
the y limits of the plot. |
xaxt |
Whether the x-axis tick marks should be displayed ( |
yaxt |
Whether the y-axis tick marks should be displayed ( |
gapFactor |
Defines spacing between SPDs as proportion of the y-axis range for multipanel plots. Default is 0.2. |
col.fill |
Vector of fill color for the observed SPDs. The default color scheme is based on the Dark2 pallette of RColorBrewer package. |
col.line |
Line colour for the observed SPDs.The default color scheme is based on the Dark2 palette of RColorBrewer package. |
lwd.obs |
Line width for the observed SPDs. Default is 1. |
lty.obs |
Line type for the observed SPDs. Default is 1. |
cex.lab |
The magnification to be used for x and y labels relative to the current setting of cex. Default is adjusted to 1. |
cex.axis |
The magnification to be used for axis annotation relative to the current setting of cex. Default is adjusted to 1. |
legend |
Whether legend needs to be displayed. Item names will be retrieved from the values supplied in the argument |
legend.arg |
list of additional arguments to pass to |
ylab |
a title for the y axis |
ymargin |
multiplier for the maximum value on ylim range. Default is 1.1. |
rnd |
integer indicating the number of decimal places to be displayed in the y-axis for when |
... |
Additional arguments affecting the plot. |
The display order of the SPDs is given by the factor levels of the user-supplied group argument in the stackspd() function.
Erich Neuwirth (2014). RColorBrewer: ColorBrewer Palettes. R package version 1.1-2. https://CRAN.R-project.org/package=RColorBrewer.
## Not run: data(emedyd) x = calibrate(x=emedyd$CRA, errors=emedyd$Error,normalised=FALSE) bins = binPrep(sites=emedyd$SiteName, ages=emedyd$CRA,h=100) res = stackspd(x=x,timeRange=c(16000,8000),bins=bins,group=emedyd$Region) plot(res,type='stacked') plot(res,type='lines') plot(res,type='proportion') plot(res,type='multipanel') ## End(Not run)## Not run: data(emedyd) x = calibrate(x=emedyd$CRA, errors=emedyd$Error,normalised=FALSE) bins = binPrep(sites=emedyd$SiteName, ages=emedyd$CRA,h=100) res = stackspd(x=x,timeRange=c(16000,8000),bins=bins,group=emedyd$Region) plot(res,type='stacked') plot(res,type='lines') plot(res,type='proportion') plot(res,type='multipanel') ## End(Not run)
Plotting function for single or multiple maps of the spatio-temporal intensity of radiocarbon dates at given focal years.
## S3 method for class 'stKde' plot( x, focalyears = NULL, type = "focal", plotdir = NULL, imnames = "byyear", zlim = NULL, box = FALSE, main = "auto", col.mask = "grey75", ncolours = 256, ramptype = "raw", imdim = c(10, 10), mars = c(0.5, 0.5, 2.5, 2), ribbon = TRUE, ribargs = list(), cex.ribbon = 0.5, withpts = "n", cex.main = 0.8, pch.pts = 19, col.pts = "grey50", cex.pts = 0.1, tidydevs = TRUE, verbose = TRUE, ... )## S3 method for class 'stKde' plot( x, focalyears = NULL, type = "focal", plotdir = NULL, imnames = "byyear", zlim = NULL, box = FALSE, main = "auto", col.mask = "grey75", ncolours = 256, ramptype = "raw", imdim = c(10, 10), mars = c(0.5, 0.5, 2.5, 2), ribbon = TRUE, ribargs = list(), cex.ribbon = 0.5, withpts = "n", cex.main = 0.8, pch.pts = 19, col.pts = "grey50", cex.pts = 0.1, tidydevs = TRUE, verbose = TRUE, ... )
x |
An object of class stKde. |
focalyears |
A vector of numeric values for focal years, in calBP, that will be timesteps at which date intensity maps will be plotted. |
type |
A single character string stipulating which type of plot to create. Current options are "nonfocal", "mask", "focal", "proportion", "change" and "all". |
plotdir |
Optional output directory for plots. If NULL, then only a single plot is made to the current device. If a valid output direcotry is provided then one or more timeslices maps are saved as png files (e.g. as source images for an animation). |
imnames |
The format of the output files if output is as png files to a directory. The current two options are "basic" (labelling the images in basic sequence as preferred by animation software such as ffmpeg) or "byyear" (labelling the images by calBP year). |
zlim |
Numeric vector of length=2 which controlls the maximum or minimum of the colour ramp. |
box |
Logical. Plot a border around the map or not. |
main |
Single character string specifying a main title. "auto" implies internal default titles are used. |
col.mask |
The colour used to depict any areas that are being masked out. |
ncolours |
The maximum number of colours to use in the colour ramp. |
ramptype |
What kind of treatment for the colour ramp. Current options are "raw" (do not try to standardise the ramps across timeslices),"std" (standardise each plot, typically by capturing the first 3sd in the main colour ramp and then outliers beyond that with the extreme colours of the ramp),"unl" (do not standardise but plot generalised high/low ramp labels) and "mmx" (scale to the minimum and maximum values across all timeslices). |
imdim |
Height and width of the plot to png in cm. |
mars |
Margins of the plot to png. |
ribbon |
Whether to plot the colour ramp legend or not. |
ribargs |
Whether to plot the colour ramp legend or not. |
cex.ribbon |
The size of the ribbon font. |
withpts |
Plot with the original sample locations shown (current options are "y" and "n"). |
cex.main |
The size of the title font. |
pch.pts |
The symbols used for the plotted points. |
col.pts |
The colours used for the plotted points. |
cex.pts |
The size used for the plotted points. |
tidydevs |
Logical for whether to clean up any open grpahics devices or not (default is TRUE). |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
... |
Additional arguments affecting the plot. |
This function plots to a screen device if a single focal year is stipulated and no output directory. Or if an output directory is stipulated, it plots one or more focal years to png files, with some basic formatting options and optional cross-year standardisation of the colour ramps (with a view to them being stills-for-video). For even more control of plotting, call this function one year at a time in a loop.
## Not run: ## Example with a subset of English and Welsh dates from the Euroevol dataset data(ewdates) data(ewowin) dir.create(file.path("im"), showWarnings=FALSE) x <- calibrate(x=ewdates$C14Age, errors=ewdates$C14SD, normalised=FALSE) bins1 <- binPrep(sites=ewdates$SiteID, ages=ewdates$C14Age, h=50) stkde1 <- stkde(x=x, coords=ewdates[,c("Eastings", "Northings")], win=ewowin, sbw=40000, cellres=2000, focalyears=seq(6500, 5000, -100), tbw=50, bins=bins1, backsight=200, outdir="im") ## Plot just the proportion surface for just 5900 calBP plot(stkde1, 5900, type="proportion") ## Plot an example of all four basic outputs for just 5900 calBP dev.new(height=2.5, width=8) par(mar=c(0.5, 0.5, 2.5, 2)) plot(stkde1, 5900, type="all") ## Plot standardised change surfaces to a sub-directory called ## /png for all timeslices and save to file. dir.create(file.path("png"), showWarnings=FALSE) plot(stkde1, seq(6500,5000,-100), type="change", ramptype="std", withpts=TRUE, plotdir="png") ## Plot all four summary surfaces in one image, saving them to a sub-directory call 'pngall', ## and with the output of the change map standardised to a common ramp ## (but leaving the focal and proportion maps unstandardised with simple ramp labelling) dir.create(file.path("pngall"), showWarnings=FALSE) plot(stkde1, seq(6500,5000,-100), type="all", ramptype=c("unl","unl","std"), imdim=cm(c(2.5,8)), withpts=TRUE, plotdir="pngall") ## End(Not run)## Not run: ## Example with a subset of English and Welsh dates from the Euroevol dataset data(ewdates) data(ewowin) dir.create(file.path("im"), showWarnings=FALSE) x <- calibrate(x=ewdates$C14Age, errors=ewdates$C14SD, normalised=FALSE) bins1 <- binPrep(sites=ewdates$SiteID, ages=ewdates$C14Age, h=50) stkde1 <- stkde(x=x, coords=ewdates[,c("Eastings", "Northings")], win=ewowin, sbw=40000, cellres=2000, focalyears=seq(6500, 5000, -100), tbw=50, bins=bins1, backsight=200, outdir="im") ## Plot just the proportion surface for just 5900 calBP plot(stkde1, 5900, type="proportion") ## Plot an example of all four basic outputs for just 5900 calBP dev.new(height=2.5, width=8) par(mar=c(0.5, 0.5, 2.5, 2)) plot(stkde1, 5900, type="all") ## Plot standardised change surfaces to a sub-directory called ## /png for all timeslices and save to file. dir.create(file.path("png"), showWarnings=FALSE) plot(stkde1, seq(6500,5000,-100), type="change", ramptype="std", withpts=TRUE, plotdir="png") ## Plot all four summary surfaces in one image, saving them to a sub-directory call 'pngall', ## and with the output of the change map standardised to a common ramp ## (but leaving the focal and proportion maps unstandardised with simple ramp labelling) dir.create(file.path("pngall"), showWarnings=FALSE) plot(stkde1, seq(6500,5000,-100), type="all", ramptype=c("unl","unl","std"), imdim=cm(c(2.5,8)), withpts=TRUE, plotdir="pngall") ## End(Not run)
Computes a combined weighted mean and error and carries out a statistical test for internal consistency.
poolDates(x, errors, id = NULL, F14C = TRUE)poolDates(x, errors, id = NULL, F14C = TRUE)
x |
A vector of uncalibrated radiocarbon ages |
errors |
A vector of standard deviations corresponding to each estimated radiocarbon age |
id |
A vector of event/object identifiers to be matched with radiocarbon ages. If not supplied all radiocarbon ages are assumed to be from the same event/object. |
F14C |
Whether calculations are carried out in F14C space or not. Default is TRUE. |
The function calculates combined weighted error mean and standard error for each set of radiocarbon ages associated with the same event or object and computes a significance test for evaluating internal consistency following Ward and Wilson's method (1978). This is equivalent to OxCal's R_Combine routine.
A data.frame containing the weighted radiocarbon ages and errors for each event and the associated T-values and P-values.
Ward, G. K., & Wilson, S. R. (1978). Procedures for Comparing and Combining Radiocarbon Age Determinations: A Critique. Archaeometry, 20(1), 19–31. https://doi.org/10.1111/j.1475-4754.1978.tb00208.x
x = c(4300,4330,5600,5603,5620) errors = c(20,30,30,30,45) id = c(1,1,2,2,2) poolDates(x,errors,id)x = c(4300,4330,5600,5603,5620) errors = c(20,30,30,30,45) id = c(1,1,2,2,2) poolDates(x,errors,id)
Function for generating a vector quantile calibrated dates from a CalDates class object.
qCal(x, p = 0.5)qCal(x, p = 0.5)
x |
A |
p |
A numeric value of probability. Default is 0.5 (median). |
A vector of quantile dates in cal BP
x <- calibrate(c(3050,2950),c(20,20)) qCal(x,p=0.2)x <- calibrate(c(3050,2950),c(20,20)) qCal(x,p=0.2)
The rcarbon package handles the calibration and analysis of radiocarbon, often but not exclusively for the purposes of archaeological research. It includes functions not only for basic calibration, uncalibration and plotting of one or more dates, but also a statistical framework for building demographic and related longitudinal inferences from aggregate radiocarbon date lists.
Core functions in the rcarbon package can be grouped as follows:
calibrate and uncalibrate enable the calibration and back-calibration for a variety of curves.
spd generates a summed probability distribution (SPD) of radiocarbon dates; binPrep can be used to define clusters of radiocarbon dates associated with the same context/phase
modelTest compares the observed SPD against a variety of theoretical models (most typically an exponential curve) using the Monte-Carlo approach; p2pTest compares observed differences in SPD between two user-specified points in time against differences expected from a theoretical model; permTest compares two or more SPDs and test for the null hypothesis that all sets are derived from the same population; sptest identifies, for defined intervals, locations with significantly higher or lower growth rate in the SPD compared to the pan-regional trend in the data
Up-to-date development version, bug-reports, and further information concerning the rcarbon package can be found on GitHub (https://github.com/ahb108/rcarbon). To see the preferred citation for the package, type citation("rcarbon").
The rcarbon is developed and maintained by Andrew Bevan and Enrico Crema
See individual functions for references.
Rescale a numeric vector to a specified minimum and maximum.
reScale(x, type = "simple", to = c(0, 1), na.rm = TRUE)reScale(x, type = "simple", to = c(0, 1), na.rm = TRUE)
x |
numeric vector to smooth. |
type |
what kind of rescaling to perform. Current options are 'simple' (default) and 'normal' which produces a z-score and 'custom' for which the 'to' argument must be specified. |
to |
numeric vector of length 2 specifying the minimum and maximum value to perform a linear rescale between (default is 0 and 1) |
na.rm |
Set to TRUE,this removes NAs before rescaling. |
A numeric vector of rescaled values.
reScale(15:200)reScale(15:200)
Calculate a running mean from a numeric vector.
runMean(x, n, edge = "NA")runMean(x, n, edge = "NA")
x |
numeric vector to smooth. |
n |
the size of the window in which to smooth. |
edge |
How to treat edge cases where a full window is unavailable. Current options are 'NA' to fill with NAs or 'fill' to fill with original values |
A numeric vector of smoothed values.
x <- rnorm(1000) y <- c(1:1000) plot(y,x, type="l") lines(runMean(x,50), col="red")x <- rnorm(1000) y <- c(1:1000) plot(y,x, type="l") lines(runMean(x,50), col="red")
Randomly samples calendar dates from each calibrated date or bin.
sampleDates(x, bins = NA, nsim, boot = FALSE, verbose = TRUE)sampleDates(x, bins = NA, nsim, boot = FALSE, verbose = TRUE)
x |
A 'CalDates' class object. |
bins |
A vector containing the bin names associated with each radiocarbon date. If set to NA, binning is not carried out. |
nsim |
Number of sampling repetitions. |
boot |
A logical value indicating whether bootstrapping is carried out (see details below). Default is FALSE. |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
The function randomly samples calendar dates based from calibrated probability distributions. When the bins argument is supplied a single calendar date is sampled from each bin. When boot=TRUE, dates (or bins) are randomly sampled with replacement before calendar dates are sampled.
An object of class simdates with the following elements
sdates A matrix containing the randomly sampled calendar dates, with rows containing each of the nsim repetitions.
weight A vector (or matrix in when boot=TRUE) containing the total area under the curve of each date, normalised to sum to unity. Notice this will be identical for all dates if the calibration is carried out with the argument normalised set to TRUE.
Smooth a numeric vector using a Gaussian window
smoothGauss(x, alpha, window = 0.1)smoothGauss(x, alpha, window = 0.1)
x |
numeric vector of values to smooth. |
alpha |
numeric value controlling the size of the gaussian smoothing window. Proportional to the standard deviation of the Gaussian smoothing kernel where sd=(N-1)/(2*alpha) with N being the length of the input vector. |
window |
a fraction between 0 and 1 representing the proportion of the input vector to include in the moving window. |
Adapted from smth.gaussian in the smoother package.
Hamilton, N. (2015). smoother: Functions Relating to the Smoothing of Numerical Data, R package version 1.1, https://CRAN.R-project.org/package=smoother
smoothGauss(runif(200),alpha=5)smoothGauss(runif(200),alpha=5)
The function generates Summed probability distributions (SPD) of radiocarbon dates, with optional binning routine for controlling inter-site or inter-phase variation in sampling intensity.
spd( x, timeRange, bins = NA, datenormalised = FALSE, spdnormalised = FALSE, runm = NA, verbose = TRUE, edgeSize = 500 )spd( x, timeRange, bins = NA, datenormalised = FALSE, spdnormalised = FALSE, runm = NA, verbose = TRUE, edgeSize = 500 )
x |
A |
timeRange |
A vector of length 2 indicating the start and end date of the analysis in cal BP. |
bins |
A vector containing the bin names associated with each radiocarbon date. If set to NA, binning is not carried out. |
datenormalised |
Controls for calibrated dates with probability mass outside the timerange of analysis. If set to TRUE the total probability mass within the time-span of analysis is normalised to sum to unity. Should be set to FALSE when the parameter |
spdnormalised |
A logical variable indicating whether the total probability mass of the SPD is normalised to sum to unity. |
runm |
A number indicating the window size of the moving average to smooth the SPD. If set to |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
edgeSize |
Extra margin in C14 Age time to handle edge effect when |
The binning routine consists of computing summed probability distribution of all dates associated to a given bin, divided by the number of contributing dates. This controls for any striking differences in sampling intensity, and ensures that each site phase is equally contributing to the final SPD (see Timpson et al 2014 for details). Bins can be generated using the binPrep, whilst the sensitivity to parameter choice can be explored with the binsense function.
An object of class CalSPD with the following elements
metadata A data.frame containing relevant information regarding the parameters used to create the SPD as well as sample size and number of bins
grid A CalGrid class object containing the summed probability associated to each calendar year between timeRange[1] and timeRange[2]
Timpson, A., et al, (2014). Reconstructing regional population fluctuations in the European Neolithic using radiocarbon dates: a new case-study using an improved method. Journal of Archaeological Science 52: 549-557. DOI:10.1016/j.jas.2014.08.011
calibrate for calibrating radiocarbon dates; binPrep for preparing bins.
Function for computing rates of change between abutting user-defined time-blocks.
spd2rc( spd, breaks = NULL, backsight = NULL, changexpr = expression((t1/t0)^(1/d) - 1) )spd2rc( spd, breaks = NULL, backsight = NULL, changexpr = expression((t1/t0)^(1/d) - 1) )
spd |
Summed Probability Distribution obtained using the |
breaks |
A vector giving the breakpoints between the time-blocks. |
backsight |
A single numeric value defining the distance in time between the focal year and the backsight year for computing the rate of change. |
changexpr |
An expression defining how the rate of change is calculated, where |
When the argument breaks is supplied the function aggregates the summed probability within each time-block and compared them across abutting blocks using the expression defined by changexpr. When the argument backsight is provided he expression is based on the comparison between the summed probability of each year and the associated backsight year.
An object of class spdRC.
## Not run: data(emedyd) caldates <- calibrate(x=emedyd$CRA, errors=emedyd$Error, normalised=FALSE, calMatrix=TRUE) bins <- binPrep(sites=emedyd$SiteName, ages=emedyd$CRA, h=50) emedyd.spd <- spd(caldates,bins,timeRange=c(16000,9000),runm=100) emedyd.gg <- spd2rc(emedyd.spd,breaks=seq(16000,9000,-1000)) emedyd.gg2 <- spd2rc(emedyd.spd,backsight=10) plot(emedyd.gg) plot(emedyd.gg2) ## End(Not run)## Not run: data(emedyd) caldates <- calibrate(x=emedyd$CRA, errors=emedyd$Error, normalised=FALSE, calMatrix=TRUE) bins <- binPrep(sites=emedyd$SiteName, ages=emedyd$CRA, h=50) emedyd.spd <- spd(caldates,bins,timeRange=c(16000,9000),runm=100) emedyd.gg <- spd2rc(emedyd.spd,breaks=seq(16000,9000,-1000)) emedyd.gg2 <- spd2rc(emedyd.spd,backsight=10) plot(emedyd.gg) plot(emedyd.gg2) ## End(Not run)
Function for mapping the spatial intensity of radiocarbon dates for a given geographical region and focal year.
spkde( x, coords, sbw, focalyear, tbw, win, cellres, bins = NA, backsight = NA, nsim = NULL, maskthresh = 0, changexpr = expression((t1 - t0)/tk), raw = FALSE, spjitter = TRUE, amount = NULL, verbose = TRUE, ... )spkde( x, coords, sbw, focalyear, tbw, win, cellres, bins = NA, backsight = NA, nsim = NULL, maskthresh = 0, changexpr = expression((t1 - t0)/tk), raw = FALSE, spjitter = TRUE, amount = NULL, verbose = TRUE, ... )
x |
An object of class CalDates with calibrated radiocarbon ages. |
coords |
A two column matrix of geographical coordinates from a a projected coordinate system (no checks are made for this) and with the same number of rows as length(x). |
sbw |
A single numeric value for the spatial bandwidth to be applied around each raster cell, expressed as the standard deviation of a continuous Gaussian kernel (passed as the sigma argument to density.ppp()). |
focalyear |
A single numeric value for the focal year for the intensity map. |
tbw |
A single numeric value for the temporal bandwidth to be applied around each focal year, expressed as the standard deviation of a continuous Gaussian kernel. |
win |
The bounding polygon for the mapping (must be an object of class 'owin', see the spatstat package) |
cellres |
The cell or pixel resolution of the output raster maps. |
bins |
A vector of labels corresponding to site names, ids, bins or phases (same length as x) |
backsight |
A single numeric value (which will be coerced to be positive) that specifies a comparison timestep in the past for a mapping of temporal change. |
nsim |
How many bootstrap simulations to run (default is none). |
maskthresh |
A single numeric value for a lower-bound cut-off for all maps, based on a minimum required spatial intensity of all dates in x. |
changexpr |
An expression for calculating the change in spatial intensity between the focal year and a backsight year (as defined via the backsight argument). Available input options are t1 (the spatial intensity for the focal year), t0 (the spatial intensity for the backsight year) and tk (the overall spatial intensity for all dates irrespective of year), plus any other standard constants and mathematical operators. A sensible default is provided. |
raw |
Whether to output the raw simulations (if nsim is set) or just the summaries (the latter is default). |
spjitter |
Whether noise is applied to the spatial coordinates or not. Default is TRUE. |
amount |
Amount of jitter applied to the spatial coordinates when |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
... |
ignored or passed to internal functions. |
This function is not really intended for general use, but rather as an internal function for stkde(). Most applications should use the latter, but spkde has been exported and made externally available, both because this function retains the result in memory (in contrast to stkde) and with a view to possible addition of bootstrap methods in the future. Some function arguments therefore remain experimental. The function computes timeslice maps of the spatio-temporal kernel intensity of radiocarbon dates across a geographic region for a specific focal year. The user specifies the arbitrary size of both the spatial and the temporal Gaussian kernels that will be used to summarise radiocarbon date intensity per grid cell.
A list object of class spKde with the following elements:
A series of list items storing some of the input parameters such as the focal year, sbw, tbw, backsight, maskthresh.
nonfocal An im object mapping the basic spatial intensity of all dates, without reference to a focal year.
focal An im object mapping the spatial intensity of dates for the focal year (i.e. weighted by how much each dates probability distribution overlaps with a Gaussian kernel centred on the focal year with a standard deviation of tbw).
proportion An im object mapping the proportional intensity of dates for the focal year (i.e. the focal surface divided by the nonfocal surface).
change An im object mapping the amount of change between the intensity of dates for the focal year and a chosen backsight year (i.e. as defined by changexpr).
## Not run: ## Example for the focal year 5600 calBP (also with site binning), ## using a subset of English and Welsh dates from the Euroevol dataset data(ewdates) data(ewowin) x <- calibrate(x=ewdates$C14Age, errors=ewdates$C14SD, normalised=FALSE) bins1 <- binPrep(sites=ewdates$SiteID, ages=ewdates$C14Age, h=50) spkde1 <- spkde(x=x, coords=ewdates[,c("Eastings", "Northings")], win=ewowin, sbw=40000, cellres=2000, focalyear=5600, tbw=50, bins=bins1, backsight=200,amount=1) plot(spkde1$focal) plot(spkde1$proportion) ## End(Not run)## Not run: ## Example for the focal year 5600 calBP (also with site binning), ## using a subset of English and Welsh dates from the Euroevol dataset data(ewdates) data(ewowin) x <- calibrate(x=ewdates$C14Age, errors=ewdates$C14SD, normalised=FALSE) bins1 <- binPrep(sites=ewdates$SiteID, ages=ewdates$C14Age, h=50) spkde1 <- spkde(x=x, coords=ewdates[,c("Eastings", "Northings")], win=ewowin, sbw=40000, cellres=2000, focalyear=5600, tbw=50, bins=bins1, backsight=200,amount=1) plot(spkde1$focal) plot(spkde1$proportion) ## End(Not run)
This function carries out local spatial permutation tests of summed radiocarbon probability distributions in order to detect local deviations in growth rates (Crema et al 2017).
sptest( calDates, timeRange, bins, locations, locations.id.col = NULL, breaks, h, kernel = "gaussian", rate = expression((t2/t1)^(1/d) - 1), nsim = 1000, runm = NA, permute = "locations", ncores = 1, datenormalised = FALSE, verbose = TRUE, raw = FALSE )sptest( calDates, timeRange, bins, locations, locations.id.col = NULL, breaks, h, kernel = "gaussian", rate = expression((t2/t1)^(1/d) - 1), nsim = 1000, runm = NA, permute = "locations", ncores = 1, datenormalised = FALSE, verbose = TRUE, raw = FALSE )
calDates |
A |
timeRange |
A vector of length 2 indicating the start and end date of the analysis in cal BP |
bins |
A vector indicating which bin each radiocarbon date is assigned to. Must have the same length as the number of radiocarbon dates. Can be created using the |
locations |
A |
locations.id.col |
Column name containing the first part of the supplied bin names. If missing rownames are used. |
breaks |
A vector of break points for defining the temporal slices. |
h |
distance parameter for calculating spatial weights. |
kernel |
indicates the type of weighting function, either 'fixed' or 'gaussian'. When set to "fixed" weight is equal to 1 within a radius of |
rate |
An expression defining how the rate of change is calculated, where |
nsim |
The total number of simulations. Default is 1000. |
runm |
The window size of the moving window average. Must be set to |
permute |
Indicates whether the permutations should be based on the |
ncores |
Number of cores used for for parallel execution. Default is 1. |
datenormalised |
A logical variable indicating whether the probability mass of each date within |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
raw |
A logical variable indicating whether permuted sets of geometric growth rates for each location should be returned. Default is FALSE. |
The function consists of the following seven steps: 1) generate a weight matrix based on inter-site distance; 2) for each location (e.g. a site) generate a local SPD of radiocarbon dates, weighting the contribution of dates from neighbouring sites using a weight scheme computed in step 1; 3) define temporal slices (using breaks as break values), then compute the total probability mass within each slice; 4) compute the rate of change between abutting temporal slices by using the formula: ; 5) randomise the location of individual bins or the entire sequence of bins associated with a given location and carry out steps 2 to 4; 6) repeat step 5 nsim times and generate, for each location, a distribution of growth rates under the null hypothesis (i.e. spatial independence); 7) compare, for each location, the observed growth rate with the distribution under the null hypothesis and compute the p-values; and 8) compute the false-discovery rate for each location.
A spatialTest class object
Crema, E.R., Bevan, A., Shennan, S. (2017). Spatio-temporal approaches to archaeological radiocarbon dates. Journal of Archaeological Science, 87, 1-9.
permTest for a non-spatial permutation test; plot.spatialTest for plotting; spd2rc for computing geometric growth rates.
## Reproduce Crema et al 2017 ## ## Not run: data(euroevol) #load data ## Subset only for ca 8000 to 5000 Cal BP (7500-4000 14C Age) euroevol2=subset(euroevol,C14Age<=7500&C14Age>=4000) ## define chronological breaks breaks=seq(8000,5000,-500) ## Create a sf class object library(sf) sites.df = unique(data.frame(SiteID=euroevol2$SiteID, Longitude=euroevol2$Longitude,Latitude=euroevol2$Latitude)) sites.sf = st_as_sf(sites.df,coords=c('Longitude','Latitude'),crs=4326) ## Calibration and binning bins=binPrep(sites=euroevol2$SiteID,ages=euroevol2$C14Age,h=200) calDates=calibrate(x=euroevol2$C14Age,errors=euroevol2$C14SD,normalised=FALSE) ## Main Analysis (over 2 cores; requires doSnow package) ## NOTE: the number of simulations should be ideally larger ## to ensure more accurate and precise estimates of the p/q-values. res.locations=sptest(calDates,timeRange=c(8000,5000),bins=bins, locations=sites.sf,locations.id.col="SiteID",h=100,breaks=breaks,ncores=2,nsim=100, permute="locations",datenormalised=FALSE) ## Plot results library(rnaturalearth) win <- ne_countries(continent = 'europe',scale=10,returnclass='sf') #retrieve coordinate limits xrange <- st_bbox(sites.sf)[c(1,3)] yrange <- st_bbox(sites.sf)[c(2,4)] par(mfrow=c(1,2)) par(mar=c(0.1,0.1,0,0.5)) plot(sf::st_geometry(win),col="antiquewhite3",border="antiquewhite3",xlim=xrange,ylim=yrange) plot(res.locations,index=4,add=TRUE,legend=TRUE,option="raw",breakRange=c(-0.005,0.005)) plot(sf::st_geometry(win),col="antiquewhite3",border="antiquewhite3",xlim=xrange,ylim=yrange) plot(res.locations,index=4,add=TRUE,legend=TRUE,option="test") ## End(Not run)## Reproduce Crema et al 2017 ## ## Not run: data(euroevol) #load data ## Subset only for ca 8000 to 5000 Cal BP (7500-4000 14C Age) euroevol2=subset(euroevol,C14Age<=7500&C14Age>=4000) ## define chronological breaks breaks=seq(8000,5000,-500) ## Create a sf class object library(sf) sites.df = unique(data.frame(SiteID=euroevol2$SiteID, Longitude=euroevol2$Longitude,Latitude=euroevol2$Latitude)) sites.sf = st_as_sf(sites.df,coords=c('Longitude','Latitude'),crs=4326) ## Calibration and binning bins=binPrep(sites=euroevol2$SiteID,ages=euroevol2$C14Age,h=200) calDates=calibrate(x=euroevol2$C14Age,errors=euroevol2$C14SD,normalised=FALSE) ## Main Analysis (over 2 cores; requires doSnow package) ## NOTE: the number of simulations should be ideally larger ## to ensure more accurate and precise estimates of the p/q-values. res.locations=sptest(calDates,timeRange=c(8000,5000),bins=bins, locations=sites.sf,locations.id.col="SiteID",h=100,breaks=breaks,ncores=2,nsim=100, permute="locations",datenormalised=FALSE) ## Plot results library(rnaturalearth) win <- ne_countries(continent = 'europe',scale=10,returnclass='sf') #retrieve coordinate limits xrange <- st_bbox(sites.sf)[c(1,3)] yrange <- st_bbox(sites.sf)[c(2,4)] par(mfrow=c(1,2)) par(mar=c(0.1,0.1,0,0.5)) plot(sf::st_geometry(win),col="antiquewhite3",border="antiquewhite3",xlim=xrange,ylim=yrange) plot(res.locations,index=4,add=TRUE,legend=TRUE,option="raw",breakRange=c(-0.005,0.005)) plot(sf::st_geometry(win),col="antiquewhite3",border="antiquewhite3",xlim=xrange,ylim=yrange) plot(res.locations,index=4,add=TRUE,legend=TRUE,option="test") ## End(Not run)
Generates and combines multiple SPDs based on a user defined grouping.
stackspd( x, timeRange, bins = NA, group = NULL, datenormalised = FALSE, runm = NA, verbose = TRUE, edgeSize = 500 )stackspd( x, timeRange, bins = NA, group = NULL, datenormalised = FALSE, runm = NA, verbose = TRUE, edgeSize = 500 )
x |
A |
timeRange |
A vector of length 2 indicating the start and end date of the analysis in cal BP. |
bins |
A vector containing the bin names associated with each radiocarbon date. If set to NA, binning is not carried out. |
group |
A character or factor vector containing the grouping variable. |
datenormalised |
Controls for calibrated dates with probability mass outside the timerange of analysis. If set to TRUE the total probability mass within the time-span of analysis is normalised to sum to unity. Should be set to FALSE when the parameter |
runm |
A number indicating the window size of the moving average to smooth the SPD. If set to |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
edgeSize |
Controls edge effect by expanding the fitted model beyond the range defined by |
An object of class stackCalSPD
## Not run: data(emedyd) x = calibrate(x=emedyd$CRA, errors=emedyd$Error,normalised=FALSE) bins = binPrep(sites=emedyd$SiteName, ages=emedyd$CRA,h=50) res = stackspd(x=x,timeRange=c(16000,8000),bins=bins,group=emedyd$Region) ## End(Not run)## Not run: data(emedyd) x = calibrate(x=emedyd$CRA, errors=emedyd$Error,normalised=FALSE) bins = binPrep(sites=emedyd$SiteName, ages=emedyd$CRA,h=50) res = stackspd(x=x,timeRange=c(16000,8000),bins=bins,group=emedyd$Region) ## End(Not run)
Function for mapping the spatio-temporal intensity of radiocarbon dates for a given geographical region in one or more tim.
stkde( x, coords, sbw, focalyears, tbw, win, cellres, outdir = ".", bins = NA, backsight = NA, maskthresh = 0, changexpr = expression((t1 - t0)/tk), spjitter = TRUE, amount = NULL, verbose = TRUE, ... )stkde( x, coords, sbw, focalyears, tbw, win, cellres, outdir = ".", bins = NA, backsight = NA, maskthresh = 0, changexpr = expression((t1 - t0)/tk), spjitter = TRUE, amount = NULL, verbose = TRUE, ... )
x |
An object of class CalDates with calibrated radiocarbon ages. |
coords |
A two column matrix of geographical coordinates from a a projected coordinate system (no checks are made for this) and with the same number of rows as length(x). |
sbw |
A single numeric value for the spatial bandwidth to be applied around each raster cell, expressed as the standard deviation of a continuous Gaussian kernel (passed as the sigma argument to density.ppp()). |
focalyears |
A vector of numeric values for focal years, in calBP, that will be timesteps at which date intensity maps will be produced. |
tbw |
A single numeric value for the temporal bandwidth to be applied around each focal year, expressed as the standard deviation of a continuous Gaussian kernel. |
win |
The bounding polygon for the mapping (must be an object of class 'owin', see the spatstat package) |
cellres |
The cell or pixel resolution of the output raster maps. |
outdir |
The output directory for timeslice maps and data that are saved to file. |
bins |
A vector of labels corresponding to site names, ids, bins or phases (same length as x) |
backsight |
A single numeric value (which will be coerced to be positive) that specifies a comparison timestep in the past for a mapping of temporal change. |
maskthresh |
A single numeric value for a lower-bound cut-off for all maps, based on a minimum required spatial intensity of all dates in x. |
changexpr |
An expression for calculating the change in spatial intensity between the focal year and a backsight year (as defined via the backsight argument). Available input options are t1 (the spatial intensity for the focal year), t0 (the spatial intensity for the backsight year) and tk (the overall spatial intensity for all dates irrespective of year), plus any other standard constants and mathematical operators. A sensible default is provided. |
spjitter |
Whether noise is applied to the spatial coordinates or not. Default is TRUE. |
amount |
Amount of jitter applied to the spatial coordinates when |
verbose |
A logical variable indicating whether extra information on progress should be reported. Default is TRUE. |
... |
ignored or passed to internal functions. |
This function computes one or more timeslice maps of the spatio-temporal kernel intensity of radiocarbon dates across a geographic region and for a specific focal year. The user specifies the arbitrary sizes of both the spatial and the temporal Gaussian kernels that will be used to summarise radiocarbon date intensity per grid cell per timestep. The results allow standardisation of colour ramps, etc. across timesteps and are amenable to plotting individually via plot.stKde and/or for output to png for animation.
A list object of class stKde with the following elements:
A series of list items storing some of the input parameters such as the focalyears, sbw, tbw, backsight, maskthresh
nonfocal An im object mapping the basic spatial intensity of all dates, without reference to a focal year.
impaths A character vector of the paths to the individual timeslices stored on file. Maps are not stored in memory (see spkde() for further details of what is stored).
stats A list of data.frames offering summary statistics on each of the different types of output surface across all timeslices. Used primarily to allow consistent colour ramps across time-slices.
ppp The ppp object for all dates and the observation window.
## Not run: ## Example with a subset of English and Welsh dates from the Euroevol dataset data(ewdates) data(ewowin) x <- calibrate(x=ewdates$C14Age, errors=ewdates$C14SD, normalised=FALSE) ## Create centennial timeslices (also with site binning) bins1 <- binPrep(sites=ewdates$SiteID, ages=ewdates$C14Age, h=50) stkde1 <- stkde(x=x, coords=ewdates[,c("Eastings", "Northings")], win=ewowin, sbw=40000, cellres=2000, focalyears=seq(6500, 5000, -100), tbw=50, bins=bins1, backsight=200, outdir="im",amount=1) ## Plot an example of all four basic outputs for 5900 calBP dev.new(height=2.5, width=8) par(mar=c(0.5, 0.5, 2.5, 2)) plot(stkde1, 5900, type="all") ## End(Not run)## Not run: ## Example with a subset of English and Welsh dates from the Euroevol dataset data(ewdates) data(ewowin) x <- calibrate(x=ewdates$C14Age, errors=ewdates$C14SD, normalised=FALSE) ## Create centennial timeslices (also with site binning) bins1 <- binPrep(sites=ewdates$SiteID, ages=ewdates$C14Age, h=50) stkde1 <- stkde(x=x, coords=ewdates[,c("Eastings", "Northings")], win=ewowin, sbw=40000, cellres=2000, focalyears=seq(6500, 5000, -100), tbw=50, bins=bins1, backsight=200, outdir="im",amount=1) ## Plot an example of all four basic outputs for 5900 calBP dev.new(height=2.5, width=8) par(mar=c(0.5, 0.5, 2.5, 2)) plot(stkde1, 5900, type="all") ## End(Not run)
Subsets calibrated dates (CalDates class object) based on Logical expressions of time intervals.
## S3 method for class 'CalDates' subset(x, s, p, ...)## S3 method for class 'CalDates' subset(x, s, p, ...)
x |
A CalDates class object |
s |
Logical expression indicating dates to keep. The expression should include the term |
p |
Probability mass meeting the condition defined by |
... |
Further arguments to be passed to or from other methods (ignored). |
The function subsets CalDates class objects by identifying all dates that have a probability mass larger than p for a user defined logical expression of temporal interval containing the term BP, where BP refers to radiocarbon date. See examples for further detailes
A CalDates class object.
## Generate some calibrated dates x = calibrate(c(12100,5410,5320,3320),errors=c(20,20,30,30)) ## Subsets all dates that have a probability mass above 0.8 before 10000 BP x2 = subset(x,BP>10000,p=0.8) ## Subsets all dates that have a probability mass above 0.5 between 6000 and 6300 BP x3 = subset(x,BP>6000&BP<6300,p=0.5)## Generate some calibrated dates x = calibrate(c(12100,5410,5320,3320),errors=c(20,20,30,30)) ## Subsets all dates that have a probability mass above 0.8 before 10000 BP x2 = subset(x,BP>10000,p=0.8) ## Subsets all dates that have a probability mass above 0.5 between 6000 and 6300 BP x3 = subset(x,BP>6000&BP<6300,p=0.5)
CalDates class objectReturns summary statistics of calibrated dates.
## S3 method for class 'CalDates' summary(object, prob = NA, calendar = "BP", ...)## S3 method for class 'CalDates' summary(object, prob = NA, calendar = "BP", ...)
object |
A |
prob |
A vector containing probabilities for the higher posterior density interval. Default is |
calendar |
Whether the summary statistics should be computed in cal BP ( |
... |
further arguments passed to or from other methods. |
A data.frame class object containing the ID of each date, along with the median date and one and two sigma (or a user specified probability) higher posterior density ranges.
SpdModelTest class objectsummary method for class "SpdModelTest"
## S3 method for class 'SpdModelTest' summary(object, type = "spd", ...)## S3 method for class 'SpdModelTest' summary(object, type = "spd", ...)
object |
A |
type |
Specifies whether the summary should be based on SPD ('spd') or associated rates of change ('roc'). Default is 'spd'. |
... |
Ignored |
The summary function returns metadata (number of radiocarbon dates, bins, and simulations), the p-value of the global significance test, and the chronological interval of local positive and negative deviations from the simulation envelope.
SpdPermTest class objectsummary method for class "SpdPermTest"
## S3 method for class 'SpdPermTest' summary(object, type = "spd", ...)## S3 method for class 'SpdPermTest' summary(object, type = "spd", ...)
object |
A |
type |
Specifies whether the summary should be based on SPD ('spd') or associated rates of change ('roc'). Default is 'spd'. |
... |
Ignored |
The summary function returns metadata (number of radiocarbon dates, bins, and simulations), the p-value of the global significance test, and the chronological interval of local positive and negative deviations from the simulation envelope for each set.
Function to select a subset of uncalibrated radiocarbon dates up to a maximum sample size per site, bin or phase.
thinDates(ages, errors, bins, size, thresh = 0.5, method = "random", seed = NA)thinDates(ages, errors, bins, size, thresh = 0.5, method = "random", seed = NA)
ages |
A vector of uncalibrated radiocarbon ages |
errors |
A vector of uncalibrated radiocarbon errors (same length as ages) |
bins |
A vector of labels corresponding to site names, ids, bins or phases (same length as ages) |
size |
A single integer specifying the maximum number of desired dates for each label stated bin. |
thresh |
A single numeric value between 0 and 1 specifying the approximate proportion (after rounding) of the resulting sample that will be chosen according to lowest date errors. At the extremes, O produces a simple random sample whereas 1 selects the sample dates with the lowest errors. Ignored if method="random". |
method |
The method to be applied where "random" simple selects a random sample, whereas "splitsample", picks some proportion (see thresh) of the sample to minimise errors, and randomly samples the rest. At present, these are the only two options. |
seed |
Allows setting of a random seed to ensure reproducibility. |
A numeric vector of the row indices corresponding to those of the input data.
data(euroevol) foursites <- euroevol[euroevol$SiteID %in% c("S2072","S4380","S6139","S9222"),] table(as.character(foursites$SiteID)) ## Thin so each site has 10 dates each max, with random selection thinInds<- thinDates(ages=foursites$C14Age, errors=foursites$C14SD, bins=foursites$SiteID, size=10, method="random", seed=123) tdates <- foursites[thinInds,] tdates ## Same but choose the first 60% (i.e. 6 dates) from the lowest errors ## and then fill in the rest randomly. thinInds<- thinDates(ages=foursites$C14Age, errors=foursites$C14SD, bins=foursites$SiteID, size=10, method="splitsample", thresh=0.6, seed=123) tdates1 <- foursites[thinInds,] tdates1data(euroevol) foursites <- euroevol[euroevol$SiteID %in% c("S2072","S4380","S6139","S9222"),] table(as.character(foursites$SiteID)) ## Thin so each site has 10 dates each max, with random selection thinInds<- thinDates(ages=foursites$C14Age, errors=foursites$C14SD, bins=foursites$SiteID, size=10, method="random", seed=123) tdates <- foursites[thinInds,] tdates ## Same but choose the first 60% (i.e. 6 dates) from the lowest errors ## and then fill in the rest randomly. thinInds<- thinDates(ages=foursites$C14Age, errors=foursites$C14SD, bins=foursites$SiteID, size=10, method="splitsample", thresh=0.6, seed=123) tdates1 <- foursites[thinInds,] tdates1
Apply taphonomic corrections or other transformations to an SPD.
transformSPD( x, correction = expression(PrDens/(5.726442 * 10^6 * (CalBP + 2176.4)^-1.3925309)) )transformSPD( x, correction = expression(PrDens/(5.726442 * 10^6 * (CalBP + 2176.4)^-1.3925309)) )
x |
An object of class |
correction |
An expression for transforming the SPD. Available input terms include: CalBP, the vector of |
An object of the same class as x
Surovell, T.A., Finley, J.B., Smith, G.M., Brantingham, P.J., Kelly, R., 2009. Correcting temporal frequency distributions for taphonomic bias. Journal of Archaeological Science 36, 1715–1724.
## Not run: data(emedyd) region1 = subset(emedyd,Region==1) x = calibrate(x=region1$CRA, errors=region1$Error,normalised=FALSE) bins = binPrep(sites=region1$SiteName, ages=region1$CRA,h=50) region1.spd = spd(x=x,bins=bins,timeRange=c(16000,8000)) region1.spd.corrected = transformSPD(region1.spd) ## End(Not run)## Not run: data(emedyd) region1 = subset(emedyd,Region==1) x = calibrate(x=region1$CRA, errors=region1$Error,normalised=FALSE) bins = binPrep(sites=region1$SiteName, ages=region1$CRA,h=50) region1.spd = spd(x=x,bins=bins,timeRange=c(16000,8000)) region1.spd.corrected = transformSPD(region1.spd) ## End(Not run)
Function for uncalibrating one or more radiocarbon dates.
uncalibrate(x, ...) ## Default S3 method: uncalibrate(x, CRAerrors = 0, roundyear = TRUE, calCurves = "intcal20", ...) ## S3 method for class 'CalGrid' uncalibrate( x, calCurves = "intcal20", eps = 1e-05, compact = TRUE, verbose = TRUE, ... )uncalibrate(x, ...) ## Default S3 method: uncalibrate(x, CRAerrors = 0, roundyear = TRUE, calCurves = "intcal20", ...) ## S3 method for class 'CalGrid' uncalibrate( x, calCurves = "intcal20", eps = 1e-05, compact = TRUE, verbose = TRUE, ... )
x |
Either a vector of calibrated radiocarbon ages or an object of class CalGrid. |
... |
ignored |
CRAerrors |
A vector of standard deviations corresponding to each estimated radiocarbon age (ignored if x is a CalGrid object). |
roundyear |
Whether the randomised estimate is rounded or not. Default is TRUE. |
calCurves |
A string naming a calibration curve already provided with the rcarbon package (currently 'intcal20','intcal13','intcal13nhpine16','shcal20','shcal13','shcal13shkauri16',”marine13','marine20' and 'normal') or a custom curve provided as matrix/data.frame in three columns ("CALBP","C14BP","Error"). The default is the 'intcal20' curve and only one curve can currently be specified for all dates. |
eps |
Cut-off value for density calculation (for CalGrid objects only). |
compact |
A logical variable indicating whether only uncalibrated ages with non-zero probabilities should be returned (for CalGrid objects only). |
verbose |
A logical variable indicating whether extra information on progress should be reported (for CalGrid objects only). |
This function takes one or more calibrated calendars and looks-up the corresponding uncalibrated age, error of the stated calibration curve at that point. It also provides a randomised estimate of the uncalibrate age based on the curve error (and optionally also a hypothetical measurement error.
A data.frame with specifying the original data, the uncalibrated age without the calibration curve error (ccCRA), the calibration curve error at this point in the curve (ccError), a randomised uncalibrated age (rCRA) given both the stated ccError and any further hypothesised instrumental error provided by the CRAerrors argument (rError).
# Uncalibrate two calendar dates uncalibrate(c(3050,2950))# Uncalibrate two calendar dates uncalibrate(c(3050,2950))
Gives the TRUE indices of calibrated dates (CalDates class object) based on Logical expressions of time intervals.
which.CalDates(x, s, p)which.CalDates(x, s, p)
x |
A CalDates class object |
s |
Logical expression indicating dates to keep. The expression should include the term |
p |
Probability mass meeting the condition defined by |
The function subsets CalDates class objects by identifying all dates that have a probability mass larger than p for a user defined logical expression of temporal interval containing the term BP, where BP refers to radiocarbon date. See examples for further detailes
A CalDates class object.
## Generate some calibrated dates x = calibrate(c(12100,5410,5320,3320),errors=c(20,20,30,30)) ## Subsets all dates that have a probability mass above 0.8 before 10000 BP x2 = which.CalDates(x,BP>10000,p=0.8) ## Subsets all dates that have a probability mass above 0.5 between 6000 and 6300 BP x3 = which.CalDates(x,BP>6000&BP<6300,p=0.5)## Generate some calibrated dates x = calibrate(c(12100,5410,5320,3320),errors=c(20,20,30,30)) ## Subsets all dates that have a probability mass above 0.8 before 10000 BP x2 = which.CalDates(x,BP>10000,p=0.8) ## Subsets all dates that have a probability mass above 0.5 between 6000 and 6300 BP x3 = which.CalDates(x,BP>6000&BP<6300,p=0.5)