Introduction
The matrix seriation problem in archaeology is based on three
conditions and two assumptions, which Dunnell
(1970) summarizes as follows.
The homogeneity conditions state that all the groups included in a
seriation must:
- Be of comparable duration,
- Belong to the same cultural tradition,
- Come from the same local area.
The mathematical assumptions state that the distribution of any
historical or temporal class:
- Is continuous through time,
- Exhibits the form of a unimodal curve.
Theses assumptions create a distributional model and ordering is
accomplished by arranging the matrix so that the class distributions
approximate the required pattern. The resulting order is
inferred to be chronological.
Reciprocal ranking
Reciprocal ranking iteratively rearrange rows and/or columns
according to their weighted rank in the data matrix until convergence
(Ihm 2005).
For a given incidence matrix \(C\):
- The rows of \(C\) are rearranged in
increasing order of:
\[ x_{i} = \sum_{j = 1}^{p} j
\frac{c_{ij}}{c_{i \cdot}} \]
- The columns of \(C\) are rearranged
in a similar way:
\[ y_{j} = \sum_{i = 1}^{m} i
\frac{c_{ij}}{c_{\cdot j}} \]
These two steps are repeated until convergence. Note that this
procedure could enter into an infinite loop.
## Build an incidence matrix with random data
set.seed(12345)
incidence1 <- matrix(sample(0:1, 400, TRUE, c(0.6, 0.4)), nrow = 20)
## Get seriation order on rows and columns
## If no convergence is reached before the maximum number of iterations (100),
## it stops with a warning.
(indices <- seriate_rank(incidence1, margin = c(1, 2), stop = 100))
#> <RankPermutationOrder>
#> Permutation order for matrix seriation:
#> - Row order: 1 4 20 3 9 16 19 10 13 2 11 7 17 5 6 18 14 15 8 12...
#> - Column order: 1 16 9 4 8 14 3 20 13 2 6 18 7 17 5 11 19 12 15 10...
## Permute matrix rows and columns
incidence2 <- permute(incidence1, indices)
## Plot matrix
heatmap(incidence1, Rowv = NA, Colv = NA, revC = TRUE,
scale = "none", col = c("white", "black"))
heatmap(incidence2, Rowv = NA, Colv = NA, revC = TRUE,
scale = "none", col = c("white", "black"))
The positive difference from the column mean percentage (in french
“écart positif au pourcentage moyen”, EPPM) represents a deviation from
the situation of statistical independence (Desachy 2004). As independence can be
interpreted as the absence of relationships between types and the
chronological order of the assemblages, EPPM is a useful graphical tool
to explore significance of relationship between rows and columns related
to seriation (Desachy 2004).
## Replicates Desachy 2004 results
data("compiegne", package = "folio")
## Get seriation order for columns on EPPM using the reciprocal averaging method
## Expected column order: N, A, C, K, P, L, B, E, I, M, D, G, O, J, F, H
(indices <- seriate_rank(compiegne, EPPM = TRUE, margin = 2))
#> <RankPermutationOrder>
#> Permutation order for matrix seriation:
#> - Row order: 1 2 3 4 5...
#> - Column order: 14 1 3 11 16 12 2 5 9 13 4 7 15 10 6 8...
## Permute columns
compiegne_permuted <- permute(compiegne, indices)
Correspondence
analysis
Seriation
Correspondence Analysis (CA) is an effective method for the seriation
of archaeological assemblages. The order of the rows and columns is
given by the coordinates along one dimension of the CA space, assumed to
account for temporal variation. The direction of temporal change within
the correspondence analysis space is arbitrary: additional information
is needed to determine the actual order in time.
## Get data
data("zuni", package = "folio")
heatmap(as.matrix(zuni), Rowv = NA, Colv = NA, revC = TRUE)
## Get row permutations from CA coordinates
(zun_indices <- seriate_average(zuni, margin = 1))
## Plot CA results
## Rows
dimensio::plot_rows(zun_indices) +
ggplot2::labs(title = "Rows") +
ggplot2::theme_bw()
## Columns
dimensio::plot_columns(zun_indices) +
ggplot2::labs(title = "Columns") +
ggplot2::theme_bw()
#> <AveragePermutationOrder>
#> Permutation order for matrix seriation:
#> - Row order: 372 387 350 367 110 417 364 407 357 160 373 406 35...
#> - Column order: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18...
## Permute data matrix
zuni_permuted <- permute(zuni, zun_indices)
heatmap(as.matrix(zuni_permuted), Rowv = NA, Colv = NA, revC = TRUE)
Refining
Peeples and Schachner (2012) propose a
procedure to identify samples that are subject to sampling error or
samples that have underlying structural relationships and might be
influencing the ordering along the CA space. This relies on a partial
bootstrap approach to CA-based seriation where each sample is replicated
n times. The maximum dimension length of the convex hull
around the sample point cloud allows to remove samples for a given
cutoff value.
According to Peeples and Schachner
(2012), “[this] point removal procedure [results in] a reduced
dataset where the position of individuals within the CA are highly
stable and which produces an ordering consistend with the assumptions of
frequency seriation.”
## Replicates Peeples and Schachner 2012 results
## Samples with convex hull maximum dimension length greater than the cutoff
## value will be marked for removal.
## Define cutoff as one standard deviation above the mean
fun <- function(x) { mean(x) + sd(x) }
## Partial bootstrap CA
## Warning: this may take a few seconds!
(zuni_boot <- dimensio::bootstrap(zun_indices, n = 1000))
#> --- Correspondence Analysis (CA) ------------------------------------------------
#> * Row variable: 420 categories (+ 420000 supplementary).
#> * Column variable: 18 categories (+ 18000 supplementary).
#>
#> * 17 components.
## Plot bootstrap CA results for the columns
dimensio::plot_columns(zuni_boot) +
ggplot2::labs(title = "Columns") +
ggplot2::theme_bw()
(zuni_keep <- refine(zuni_boot, cutoff = fun, margin = 1))
#> <RefineCA>
#> Partial bootstrap CA seriation refinement:
#> - Cutoff value: 2.55
#> - Rows to keep: 346 of 420 (82%)
## Plot convex hull
## blue: convex hull for samples; red: convex hull for types
### All bootstrap samples
rchull <- as.data.frame(zuni_keep[["hull"]])
ggplot2::ggplot(data = rchull) +
ggplot2::aes(x = x, y = y, group = id) +
ggplot2::geom_polygon(fill = "blue", alpha = 0.1) +
ggplot2::geom_vline(xintercept = 0, size = 0.5, linetype = "dashed") +
ggplot2::geom_hline(yintercept = 0, size = 0.5, linetype = "dashed") +
ggplot2::coord_fixed() +
ggplot2::labs(x = "", y = "") +
ggplot2::theme_bw()
### Only retained samples
rchull_sub <- subset(rchull, id %in% zuni_keep[["keep"]])
ggplot2::ggplot(data = rchull_sub) +
ggplot2::aes(x = x, y = y, group = id) +
ggplot2::geom_polygon(fill = "blue", alpha = 0.1) +
ggplot2::geom_vline(xintercept = 0, size = 0.5, linetype = "dashed") +
ggplot2::geom_hline(yintercept = 0, size = 0.5, linetype = "dashed") +
ggplot2::coord_fixed() +
ggplot2::labs(x = "", y = "") +
ggplot2::theme_bw()
## Histogram of convex hull maximum dimension length
hull_length <- data.frame(length = zuni_keep[["length"]])
ggplot2::ggplot(data = hull_length) +
ggplot2::aes(x = length) +
ggplot2::geom_histogram(breaks = seq(0, 4.5, by = 0.5), fill = "grey70") +
ggplot2::geom_vline(xintercept = fun(hull_length$length), colour = "red") +
ggplot2::labs(
title = "Convex hull max. dim.",
x = "Maximum length",
y = "Count"
) +
ggplot2::theme_bw()
