modelTest: Monte-Carlo simulation test for SPDs

Description Usage Arguments Details Value Note References Examples

View source: R/tests.R


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,
  raw = FALSE, model = c("exponential"), method = c("uncalsample"),
  predgrid = NA, datenormalised = FALSE, spdnormalised = FALSE,
  ncores = 1, fitonly = FALSE, a = 0, b = 0, verbose = TRUE)



A CalDates object containing calibrated radiocarbon ages


A vector of errors corresponding to each radiocarbon age


Number of simulations


A vector indicating which bin each radiocarbon date is assigned to.


A number indicating the window size of the moving average to smooth both observed and simulated SPDs. If set to NA no moving average is applied.Default is NA.


A vector of length 2 indicating the start and end date of the analysis in cal BP.


A logical variable indicating whether all permuted SPDs should be returned or not. Default is FALSE.


A vector indicating the model to be fitted. Currently the acceptable options are 'uniform', 'linear', 'exponential' and 'custom'. Default is 'exponential'.


Method for the creation of random dates from the fitted model. Either 'uncalsample' or 'calsample'. Default is 'uncalsample'. See below for details.


A data.frame containing calendar years (column calBP) and associated summed probabilities (column PrDens). Required when model is set to 'custom'.


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 x were calibrated without normalisation (via normalised=FALSE in the calibrate function), in which case setting datenormalised=TRUE here will rescale each dates probability mass to sum to unity before aggregating the dates, while setting datenormalised=FALSE will ensure unnormalised dates are used for both observed and simulated SPDs. Default is FALSE.


A logical variable indicating whether the total probability mass of the SPD is normalised to sum to unity for both observed and simulated data.


Number of cores used for for parallel execution. Default is 1.


A logical variable. If set to TRUE, only the the model fitting is executed and returned. Default is FALSE.


Starter value for the exponential fit with the nls function using the formula y ~ exp(a + b * x) where y is the summed probability and x is the date. Default is 0.


Starter value for the exponential fit with the nls function using the formula y ~ exp(a + b * x) where y is the summed probability and x is the date. Default is 0.


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). 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). Alternatively, it is 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 n samples, with n 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 'intcal13' the probability of that particular bin being assigned to 'intcal13' is 0.75.


An object of class SpdModelTest with the following elements


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.


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: 
caldates <- calibrate(x=emedyd$CRA, errors=emedyd$Error, normalised=FALSE, calMatrix=TRUE)
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

## End(Not run)

rcarbon documentation built on Oct. 1, 2018, 5:03 p.m.