rcarbon (version 1.5.1)

sptest: 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
)

Value

A spatialTest class object

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.

References

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

See Also

permTest for a non-spatial permutation test; plot.spatialTest for plotting; spd2rc for computing geometric growth rates.

Examples

## Reproduce Crema et al 2017 ##
if (FALSE) {
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")  
}