sptest: Spatial Permutation Test of summed probability distributions.
In rcarbon: Calibration and Analysis of Radiocarbon Dates

sptest

R Documentation

Spatial Permutation Test of summed probability distributions.

Description

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).

Usage

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
)

Arguments

`calDates`	A `CalDates` class object.
`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 `binPrep`) function. Bin names should follow the format "x_y", where x refers to a unique location (e.g. a site) and y is a integer value (e.g. "S023_1", "S023_2","S034_1", etc.).
`locations`	A `sf` class object of the site locations.
`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 `h` km from each focal sie, and 0 outside this range. When set "gaussian" the weight declines with an exponential rate equal to exp(-d^2/h^2), where the d is the distance to the focal site. Default is "gaussian".
`rate`	An expression defining how the rate of change is calculated, where `t1` is the summed probability for a focal block, `t2` is the summed probability for next block, and `d` is the duration of the blocks. Default is a geometric growth rate (i.e `expression((t2/t1)^(1/d)-1)`).
`nsim`	The total number of simulations. Default is 1000.
`runm`	The window size of the moving window average. Must be set to `NA` if the rates of change are calculated from the raw SPDs.
`permute`	Indicates whether the permutations should be based on the `"bins"` or the `"locations"`. Default is `"locations"`.
`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 `timeRange` is equal to 1. Default is FALSE.
`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.

Details

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: (SPD_{t}/SPD_{t+1}^{1/\Delta t}-1); 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.

Value

A spatialTest class object

References

Crema, E.R., Bevan, A., Shennan, S. (2017). Spatio-temporal approaches to archaeological radiocarbon dates. Journal of Archaeological Science, 87, 1-9.

Examples

## 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)

rcarbon documentation built on Aug. 24, 2023, 5:11 p.m.