rcarbon (version 1.4.0)

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.

Value

A spatialTest class object

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/\Delta 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.

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; spweights for computing spatial weights; spd2rc for computing geometric growth rates.

Examples

# NOT RUN {
## 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")  
# }