rcarbon: sptest – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

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,
  breaks,
  spatialweights,
  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 `SpatialPoints` or a `SpatialPointsDataFrame` class object. Rownames of each point should much the first part of the bin names supplied (e.g. "S023","S034")
`breaks`	A vector of break points for defining the temporal slices.
`spatialweights`	A `spatialweights` class object defining the spatial weights between the locations (cf. `spweights`)
`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) 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 provided by the spatialweights class object; 2) define temporal slices (using breaks as break values), then compute the total probability mass within each slice; 3) compute the rate of change between abutting temporal slices by using the formula: (SPD_{t}/SPD_{t+1}^{1/Δ t}-1); 4) randomise the location of individual bins or the entire sequence of bins associated with a given location and carry out steps 1 to 3; 5) repeat step 4 nsim times and generate, for each location, a distribution of growth rates under the null hypothesis (i.e. spatial independence); 6) compare, for each location, the observed growth rate with the distribution under the null hypothesis and compute the p-values; and 7) 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 8000 to 5000 Cal BP (c7200-4200 C14BP)
edge=800
timeRange=c(8000,5000)
euroevol2=subset(euroevol,C14Age<=c(timeRange[1]-edge)&C14Age>=c(timeRange[2]-edge))

## define chronological breaks
breaks=seq(8000,5000,-500)

## Create a SpatialPoints class object 
library(sp)
sites = unique(data.frame(SiteID=euroevol2$SiteID,
Longitude=euroevol2$Longitude,Latitude=euroevol2$Latitude))
locations=data.frame(Longitude=sites$Longitude,Latitude=sites$Latitude)
rownames(locations)=sites$SiteID
coordinates(locations)<-c("Longitude","Latitude")
proj4string(locations)<- CRS("+proj=longlat +datum=WGS84")

## Compute Distance and Spatial Weights 
distSamples=spDists(locations,locations,longlat = TRUE)
spatialweights=spweights(distSamples,h=100) #using a kernal bandwidth of 100km

## Calibration and binning
bins=binPrep(sites=euroevol2$SiteID,ages=euroevol2$C14Age,h=200)  
calDates=calibrate(x=euroevol2$C14Age,errors=euroevol2$C14SD,
timeRange=timeRange,normalised=FALSE)

## Main Analysis (over 2 cores; requires doSnow package) 
## NOTE: the number of simulations should be ideally larger 
## to ensure a better resolution of the p/q-values.
res.locations=sptest(calDates,timeRange=timeRange,bins=bins,locations=locations,
spatialweights=spatialweights,breaks=breaks,ncores=2,nsim=100,
permute="locations",datenormalised=FALSE)

## Plot results
library(rworldmap)
base=getMap(resolution="low") #optionally add base map
#retrieve coordinate limits#
xrange=bbox(res.locations$locations)[1,] 
yrange=bbox(res.locations$locations)[2,]

par(mfrow=c(2,2))  
par(mar=c(0.1,0.1,0,0.5))
plot(base,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)) 
par(mar=c(0.1,0.1,0,0.5))
plot(base,col="antiquewhite3",border="antiquewhite3",xlim=xrange,ylim=yrange) 
plot(res.locations,index=4,add=TRUE,legend=TRUE,option="test")  

## End(Not run)

rcarbon

Calibration and Analysis of Radiocarbon Dates

v1.4.2

GPL (>= 2)

Authors

Andrew Bevan [aut] (<https://orcid.org/0000-0001-7967-3117>), Enrico Crema [aut, cre] (<https://orcid.org/0000-0001-6727-5138>), R. Kyle Bocinsky [ctb], Martin Hinz [ctb], Philip Riris [ctb], Fabio Silva [ctb]

Initial release