Seriation

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(ca)
library(plotrix)
library(kableExtra)

Attaching package: 'kableExtra'

The following object is masked from 'package:dplyr':

    group_rows

# read data
data <- read.csv("seriate.allen.phase.csv")

# Row sums
(RS <- rowSums(data[, -(1)]))

[1]  1  7 44 13 27

# Column sums
CS <- colSums(data[, -(1)])
as.matrix(CS)

          [,1]
Poynor      21
Patton      45
Hume.Eng    10
H.L.Eff      5
Kill.Pin     1
Bull.Brus    7
Plain        3

burial_ct <- data[, -(1)]
labels <- paste0(data$Burial)
rownames(burial_ct) <- labels
burial_pct <- burial_ct / RS  * 100

burial_ca <- ca(burial_ct)
plot.ca(burial_ca, labels=c(0, 2), cex = .75)

# Correspondence Analysis Ordering on Dimension 1
(ca_ord <- order(burial_ca$rowcoord[, 1]))

[1] 3 2 1 4 5

library(vegan)

Loading required package: permute

Loading required package: lattice

This is vegan 2.6-4

burial_dec <- decorana(burial_ct)
# Decorana Ordering on Projection 1
(dec_ord <- order(burial_dec$rproj[, 1]))

[1] 4 5 1 2 3

(cor.test(ca_ord, dec_ord, method = "kendall"))

    Kendall's rank correlation tau

data:  ca_ord and dec_ord
T = 5, p-value = 1
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau 
  0 

detach("package:vegan")
detach("package:permute")

library(seriation)

Registered S3 method overwritten by 'seriation':
  method         from 
  reorder.hclust vegan

Attaching package: 'seriation'

The following object is masked from 'package:lattice':

    panel.lines

burial_dist <- dist(burial_pct, method="manhattan")
# Distance matrix:
round(burial_dist, 1)

                 AN21-E-Owens-1Bur CE6-Hackney-2Bur CE12-Allen-18Bur
CE6-Hackney-2Bur              85.7                                  
CE12-Allen-18Bur             118.2             53.2                 
AN13-Jowell-Pur              123.1             65.9             73.1
AN26-Patton-Pur               74.1             57.1             74.4
                 AN13-Jowell-Pur
CE6-Hackney-2Bur                
CE12-Allen-18Bur                
AN13-Jowell-Pur                 
AN26-Patton-Pur             71.2

set.seed(42)
burial_ser1 <- seriate(burial_dist, method="ARSA")
burial_ser2 <- seriate(burial_dist, method="BBURCG")
burial_ser3 <- seriate(burial_dist, method="BBWRCG")
ord1 <- get_order(burial_ser1)
ord2 <- get_order(burial_ser2)
ord3 <- get_order(burial_ser3)
# Suggested orderings
rbind(ord1, ord2, ord3)

     AN13-Jowell-Pur CE12-Allen-18Bur CE6-Hackney-2Bur AN26-Patton-Pur
ord1               4                3                2               5
ord2               1                5                2               3
ord3               1                5                2               3
     AN21-E-Owens-1Bur
ord1                 1
ord2                 4
ord3                 4

Crit <- c("AR_events", "AR_deviations", "Gradient_raw", "Gradient_weighted")
burial_ord1 <- permute(burial_dist, ord1)
burial_ord2 <- permute(burial_dist, ord2)
burial_ord3 <- permute(burial_dist, ord3)
gahagan_dec_ord <- permute(burial_dist, dec_ord)
gahagan_ca_ord <- permute(burial_dist, ca_ord)

# Comparison between original and ordered:
Results <- list(Original=burial_dist, Ordered1=burial_ord1, Ordered2=burial_ord2, Ordered3=burial_ord3, CA=gahagan_ca_ord, DEC=gahagan_dec_ord)
round(sapply(Results, criterion, method=Crit), 1)

                  Original Ordered1 Ordered2 Ordered3     CA   DEC
AR_events              8.0      3.0      3.0      3.0   10.0   9.0
AR_deviations        152.4     12.2     12.2     12.2  333.2 267.9
Gradient_raw           4.0     14.0     14.0     14.0    0.0   2.0
Gradient_weighted    135.3    534.1    534.1    534.1 -162.7   6.7

battleship.plot(burial_pct[ord2, ], 
                col="#798E87", 
                cex.labels = .75)

Recalibrate legacy date

# check for update
library(oxcAAR)
library(rcarbon)
quickSetupOxcal()

NULL

# an21
an21 <- data.frame(bp = c(300),
                   std = c(40),
                   names = c("41AN21")
)
an21.cal <- oxcalCalibrate(an21$bp,
                           an21$std,
                           an21$names)
an21.cal

=============================
    R_Date: 41AN21
=============================


BP = 300, std = 40


unmodelled:                    posterior:

one sigma                      
1516 AD - 1590 AD (50.37%)     
1620 AD - 1648 AD (17.9%)      

two sigma                      
1480 AD - 1660 AD (95.45%)     

three sigma                    
1450 AD - 1678 AD (97.92%)     
1764 AD - 1800 AD (1.81%)      

Calibrated with:
     Atmospheric data from Reimer et al (2020) 

## plot
plot(an21.cal)