---
title: "palimpsestr: Probabilistic Decomposition of Archaeological Palimpsests Using Stratigraphic Entanglement Fields"
author: "Enzo Cocca"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{palimpsestr: Probabilistic Decomposition of Archaeological Palimpsests}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5.5, dpi = 150,
                      warning = FALSE, message = FALSE)
```

*Archaeological palimpsests --- deposits in which material from multiple occupation phases is superimposed and partially intermixed --- represent one of the most persistent analytical challenges in field archaeology. The **palimpsestr** R package implements the Stratigraphic Entanglement Field (SEF) framework, a probabilistic method for decomposing palimpsests into latent depositional phases by jointly modelling four evidence domains: spatial proximity, vertical distribution, chronological overlap, and cultural similarity. The package provides a mixed-type finite mixture model --- a diagonal-covariance Gaussian over the numeric features combined with a per-phase categorical distribution over the cultural class --- fitted via Expectation-Maximisation, augmented with taphonomic weighting, a uniform noise/outlier component, and stratigraphic entanglement constraints. Three interpretable diagnostics --- the Stratigraphic Entanglement Index (SEI), Excavation Stratigraphic Energy (ESE), and Palimpsest Dissolution Index (PDI) --- allow practitioners to assess deposit coherence, detect intrusive finds, and evaluate overall phase separability. This article describes the methodology and demonstrates a complete analytical workflow on a realistic multi-period Roman villa dataset.*

## Introduction

Archaeological palimpsests arise when successive episodes of human occupation leave material traces within the same sedimentary matrix, producing deposits in which artefacts, ecofacts, and features from chronologically distinct phases co-occur in close spatial and vertical proximity. The phenomenon is ubiquitous: from deeply stratified Near Eastern tells to shallow open-air Palaeolithic sites where millennia of activity are compressed into a few centimetres of deposit. Disentangling the depositional history of such contexts is essential for any behavioural, technological, or settlement-pattern inference, yet the analytical toolkit available to excavators has remained surprisingly limited.

The Harris Matrix (Harris, 1989) provides a rigorous formalism for recording observed stratigraphic relationships, but its relationships are deterministic and binary --- a unit either overlies another or it does not. It cannot express the gradual, probabilistic mixing that characterises many palimpsests. Conversely, spatial clustering methods such as k-means (Kintigh and Ammerman, 1982) operate on coordinate data alone, ignoring vertical stratigraphy, chronological evidence from radiocarbon or ceramic seriation, and cultural-typological information. Neither approach offers a unified probabilistic framework, and neither produces quantitative diagnostics for deposit integrity or mixing intensity.

**palimpsestr** addresses these limitations by integrating four evidence domains --- horizontal coordinates, vertical elevation, chronological range, and cultural class --- into a single model that produces soft (probabilistic) phase assignments together with three diagnostic statistics. The Stratigraphic Entanglement Index (SEI) quantifies pairwise depositional coherence; Excavation Stratigraphic Energy (ESE) captures local disruption; and the Palimpsest Dissolution Index (PDI) summarises global phase separability. All three are interpretable by field practitioners and can be exported to GIS for spatial interrogation.

## Methodology

### The Stratigraphic Entanglement Index (SEI)

The SEI quantifies how strongly two finds $i$ and $j$ are "entangled" --- that is, how much evidence supports their co-occurrence within the same depositional event. It is defined as:

$$SEI_{ij} = \frac{w_s}{d_{xy}(i,j) + \epsilon} + \frac{w_z}{\max(|z_i - z_j|, z_{floor})} + w_t \cdot O_t(i,j) + w_c \cdot \mathbb{1}[c_i = c_j]$$

where:

- $d_{xy}(i,j)$ is the Euclidean distance in the horizontal plane, and $\epsilon$ is a small constant preventing division by zero. The first term captures spatial proximity: finds close together in plan are more likely to derive from the same depositional event.
- $|z_i - z_j|$ is the absolute vertical separation, with $z_{floor}$ a minimum threshold. The second term captures stratigraphic proximity: finds at the same depth are more likely contemporaneous.
- $O_t(i,j)$ is the chronological overlap ratio, computed from the intersection of the temporal ranges $[t_{min,i}, t_{max,i}]$ and $[t_{min,j}, t_{max,j}]$ relative to their union. Finds whose date ranges overlap strongly receive higher entanglement scores.
- $\mathbb{1}[c_i = c_j]$ is an indicator function equal to 1 when both finds belong to the same cultural or typological class, and 0 otherwise.
- $w_s, w_z, w_t, w_c$ are domain weights (defaulting to equal weighting) that allow the analyst to encode prior knowledge about which evidence domains are most reliable for a given site.

Each SEI component is normalised to [0, 1] by dividing by its within-dataset maximum, making weights directly interpretable as relative importance. The SEI matrix is row-normalised and used both as a regularisation penalty during model fitting and as the basis for the ESE diagnostic.

### The Diagonal Gaussian Mixture EM

Phase decomposition proceeds via a Gaussian mixture model fitted by Expectation-Maximisation (EM). The feature matrix for each find is constructed by concatenating: scaled horizontal coordinates $(x, y)$, vertical coordinate $(z)$, temporal midpoint and temporal span of the date range, and a one-hot encoding of the cultural class variable. This multi-domain feature representation ensures that the clustering draws on all available evidence simultaneously.

The algorithm is initialised with K-means, and the resulting hard assignments are converted to softmax probabilities. The EM procedure then iterates between an E-step (computing posterior phase membership probabilities given current parameters) and an M-step (updating mixture weights, means, and diagonal covariance matrices given current memberships). Two modifications distinguish SEF from a standard Gaussian mixture:

1. **Taphonomic weighting**: if a taphonomic integrity score is available for each find, it is used to down-weight poorly preserved or disturbed items during parameter estimation, reducing their influence on cluster centroids.
2. **Stratigraphic entanglement penalties**: the SEI matrix is incorporated as a soft constraint, encouraging finds with high mutual entanglement to receive similar phase assignments.

Convergence is assessed by monitoring the log-likelihood; the algorithm terminates when the relative change falls below a tolerance threshold or a maximum number of iterations is reached.

### Diagnostic Statistics

Three statistics summarise different aspects of the decomposition:

- **Excavation Stratigraphic Energy (ESE)**: for each find, ESE measures its dissimilarity from its spatial neighbours in terms of phase assignment. A high ESE value indicates that a find is surrounded by items assigned to different phases --- a signature of depositional disruption, bioturbation, or post-depositional displacement.

- **Palimpsest Dissolution Index (PDI)**: a global statistic defined as $PDI = 1 - \bar{H} / \log(K)$, where $\bar{H}$ is the mean Shannon entropy of the phase-probability vectors across all finds and $K$ is the number of phases. PDI ranges from 0 (complete mixing; every find is equally likely to belong to any phase) to 1 (perfect separation; every find is assigned to exactly one phase with certainty). Values above 0.7 generally indicate well-resolved phases.

- **Entropy**: the Shannon entropy of each find's phase-probability vector, $H_i = -\sum_k p_{ik} \log p_{ik}$. Low entropy means confident assignment; high entropy signals ambiguity, typically at phase boundaries or within heavily mixed zones.

- **Intrusion score**: a composite metric combining high entropy, high ESE, and low SEI connectivity. Finds exceeding a threshold on this composite are flagged as candidate intrusions --- items that are spatially or stratigraphically anomalous relative to their assigned phase.

## The Dataset: Villa Romana

We demonstrate the complete analytical workflow using `villa_romana`, a real archaeological dataset from Poggio Gramignano (Lugnano in Teverina, TR), a multi-period Roman villa with an annexed Late Antique infant cemetery. The dataset contains 615 inventoried finds from 54 stratigraphic units, with taphonomic scores assigned by the excavator.

```{r load-data}
library(palimpsestr)
data(villa_romana)
cat("Finds:", nrow(villa_romana), "\n")
cat("Stratigraphic units:", length(unique(villa_romana$context)), "\n")
cat("Material classes:", length(unique(villa_romana$class)), "\n")
cat("Date range:", min(villa_romana$date_min), "to", max(villa_romana$date_max), "\n")
cat("Mean taf_score:", round(mean(villa_romana$taf_score), 3), "\n")
```

The site spans eight archaeological periods:

| Period | Chronology | Description |
|:-------|:-----------|:------------|
| 10 | 650--350 BCE | Pre-Roman (foundation pits) |
| 9 | 30 BCE--80 CE | Augustan--Flavian (villa construction) |
| 8 | 80--137 CE | Flavian |
| 7 | 138--179 CE | Antonine--Severan |
| 6 | 180--324 CE | Severan--Late Antique |
| 5 | 324--425 CE | Late Antique |
| 4 | 425--450 CE | Late Antique (cemetery) |
| 3 | 450--550 CE | Late Antique--Early Medieval |

The deposit includes levelling fills (US 76, 117, 118), backfills, and redeposited material, reflected in the taphonomic scores (taf_score: 0.5 for disturbed, 0.7 for redeposited, 1.0 for in situ).

```{r data-summary}
cat("Material classes:\n")
sort(table(villa_romana$class), decreasing = TRUE)
```

```{r data-structure}
str(villa_romana)
```

## Phase Count Selection

When the number of phases is unknown, `compare_k()` fits models across a range of $K$ values and reports BIC, ICL, and PDI for each. The optimal $K$ minimises BIC while maximising PDI.

```{r compare-k}
ck <- compare_k(
  villa_romana,
  k_values = 2:7,
  tafonomy = "taf_score",
  context  = "context",
  seed     = 42
)
print(ck)
```

```{r gg-compare-k, fig.cap="Model selection diagnostics: BIC, PDI, and mean entropy across candidate phase counts. The dashed line marks the optimal K."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_compare_k(ck)
}
```

## Fitting the Model (K = 4)

Based on the selection diagnostics and archaeological knowledge of the site (4 occupation periods), we fit with K = 4 and multiple random initialisations.

```{r fit-model}
fit <- fit_sef(
  villa_romana,
  k        = 4,
  tafonomy = "taf_score",
  context  = "context",
  n_init   = 10,
  seed     = 42
)
print(fit)
```

```{r model-summary}
summary(fit)
```

### Model Overview

```{r model-stats}
cat("PDI:", pdi(fit), "\n")
sef_summary(fit)
```

### EM Convergence

```{r gg-convergence, fig.cap="EM convergence trace showing log-likelihood stabilisation across iterations."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_convergence(fit)
}
```

## Spatial Diagnostics

### Phase-field Map

```{r gg-phasefield, fig.cap="Dominant phase assignment. Point size reflects assignment confidence (inverse entropy). Larger points indicate higher confidence."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_phasefield(fit)
}
```

The phase-field map reveals the spatial structure of the decomposition. The model assigns each find to one of four chronologically coherent phases despite the presence of post-depositional disturbance.

### Entropy Map

```{r gg-entropy, fig.cap="Spatial distribution of Shannon entropy. High-entropy zones (bright) mark areas of uncertain phase assignment, typically at phase boundaries or in heavily mixed deposits."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_entropy(fit)
}
```

Entropy mapping is particularly useful for identifying the spatial extent of mixing. Zones of elevated entropy correspond to areas where multiple depositional phases are interleaved and cannot be cleanly separated.

### Energy Map (ESE)

```{r gg-energy, fig.cap="Excavation Stratigraphic Energy field. Elevated ESE values (bright) highlight zones of depositional disruption where neighbouring finds belong to different phases."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_energy(fit)
}
```

The energy map highlights local disruptions. Clusters of high-ESE finds may indicate bioturbation channels, pit cuts, or other post-depositional disturbances.

### Intrusion Detection Map

```{r gg-intrusions, fig.cap="Candidate intrusions overlaid on the phase-field map. Circled finds exhibit high entropy, high stratigraphic energy, and low entanglement with their neighbours. Top 10 suspects are labelled."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_intrusions(fit, top_n = 10)
}
```

### Vertical Phase Profile

```{r gg-phase-profile, fig.cap="Vertical phase profile: depth vs horizontal position. Phase 1 (deepest) corresponds to the earliest occupation; Phase 4 (shallowest) to the most recent."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_phase_profile(fit)
}
```

The vertical profile confirms the expected stratigraphic ordering: earlier phases occupy deeper strata, later phases are shallower.

## Model-output diagnostic plots

```{r model-plots, eval=FALSE}
fit <- fit_sef(demo_moderate, k = 3, context = "context", noise = TRUE)
gg_phase_composition(fit)        # per-phase class profile (heatmap)
gg_unit_coherence(fit)           # within-unit coherence
gg_outliers(fit)                 # intrusion ranking (noise posterior)
gg_direction(fit)                # residual vs latent finds
```

## Intrusion Analysis

The `detect_intrusions()` function computes a composite intrusion probability for each find, combining entropy, ESE, and SEI connectivity into a single score.

```{r intrusions}
di <- detect_intrusions(fit)
suspects <- di[di$intrusion_prob > 0.5, ]
cat("Suspected intrusions:", nrow(suspects), "out of", nrow(villa_romana), "finds",
    sprintf("(%.1f%%)\n", 100 * nrow(suspects) / nrow(villa_romana)))
```

```{r intrusion-table}
if (nrow(suspects) > 0) {
  # Merge with source data and phase assignments
  phase_tab <- as_phase_table(fit)
  susp_data <- merge(suspects, villa_romana[, c("id", "context", "class")], by = "id")
  susp_data <- merge(susp_data, phase_tab[, c("id", "dominant_phase")], by = "id")
  susp_data <- susp_data[order(susp_data$intrusion_prob, decreasing = TRUE), ]
  susp_show <- susp_data[, c("id", "context", "class", "dominant_phase", "intrusion_prob")]
  names(susp_show) <- c("Find", "US", "Class", "Phase", "Intrusion Prob.")
  knitr::kable(head(susp_show, 15), row.names = FALSE, digits = 3,
               caption = "Top suspected intrusions ranked by intrusion probability.")
}
```

Flagged finds warrant closer examination: they may represent genuinely intrusive material (e.g., items displaced by animal burrowing or construction cuts), or they may sit at genuine phase boundaries where assignment is inherently ambiguous.

## Stratigraphic Unit Purity

The `us_summary_table()` function aggregates diagnostics per stratigraphic unit, reporting dominant phase, purity, mean entropy, energy, and intrusion count.

```{r us-purity}
us_tab <- us_summary_table(fit)
knitr::kable(us_tab, row.names = FALSE, digits = 3,
             caption = "Per-US diagnostics: dominant phase, purity, mean entropy, mean ESE, and number of intrusions.")
```

Units with purity > 90% indicate consistent phase assignment within the context. Lower purity and higher mean entropy correspond to disturbed contexts where bioturbation or construction cuts have introduced material from adjacent phases.

```{r purity-stats}
n_pure <- sum(us_tab$purity >= 0.9)
cat(n_pure, "out of", nrow(us_tab), "units have purity >= 90%\n")
```

## Phase Transition Matrix

Shows how often finds from phase $i$ sit directly above finds from phase $j$, revealing spatial adjacency between phases:

```{r transition-matrix}
tm <- phase_transition_matrix(fit)
print(tm)
```

High off-diagonal values indicate zones of contact between phases. This is expected at the boundary between successive occupation periods.

## Model Validation

When true phase assignments are known (e.g., from simulation), palimpsestr provides validation metrics. We demonstrate with the simulated `demo_moderate` dataset:

```{r validation}
data(demo_moderate)
fit_val <- fit_sef(demo_moderate, k = 3, tafonomy = "taf_score", context = "context", seed = 42)
ari <- adjusted_rand_index(fit_val, demo_moderate$true_phase)
cat("Adjusted Rand Index:", round(ari, 3), "\n")
confusion_matrix(fit_val, demo_moderate$true_phase)
```

```{r gg-confusion, fig.cap="Confusion matrix heatmap on simulated data. Strong diagonal = correct phase recovery."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_confusion(fit_val, demo_moderate$true_phase)
}
```

For the real Poggio Gramignano data, no ground truth is available. Validation relies on archaeological coherence: the model should recover phases consistent with the site periodisation, and flag known disturbed contexts (levelling fills, backfills) as mixed.

## Sensitivity Analysis of Domain Weights

The model combines four evidence domains (spatial, vertical, temporal, cultural) with configurable weights `c(ws, wz, wt, wc)`. Since version 0.14.0 these weights scale the corresponding feature dimensions and therefore enter the mixture likelihood directly (in earlier versions they affected only the SEI matrix). To assess how weight choices affect the results, we compare three configurations:

```{r sensitivity}
configs <- list(
  equal    = c(ws = 1, wz = 1, wt = 1, wc = 1),
  spatial  = c(ws = 2, wz = 1, wt = 0.5, wc = 0.5),
  temporal = c(ws = 0.5, wz = 0.5, wt = 2, wc = 1)
)

sens <- do.call(rbind, lapply(names(configs), function(nm) {
  w <- configs[[nm]]
  f <- fit_sef(villa_romana, k = 4, weights = w, seed = 42,
               tafonomy = "taf_score", context = "context")
  di <- detect_intrusions(f)
  data.frame(
    config = nm,
    pdi = pdi(f),
    mean_entropy = mean(f$entropy, na.rm = TRUE),
    n_intrusions = sum(di$intrusion_prob > 0.5),
    stringsAsFactors = FALSE
  )
}))
knitr::kable(sens, digits = 4,
             caption = "Sensitivity of PDI, mean entropy, and intrusion count to weight configuration.")
```

Because the weights now reshape the feature space, they genuinely change the recovered phases — emphasising the most informative domain for a given site improves separability, while down-weighting a noisy domain reduces spurious structure. The `optimize_weights()` function provides a data-driven approach to weight selection via cross-validated held-out log-likelihood; it fits under `var_structure = "spherical"` (the structure under which the weights are statistically identifiable) and corrects for the Jacobian of the weighting so that configurations are directly comparable. The returned weights are intended to be used with `fit_sef(..., var_structure = "spherical")`.

## Bootstrap Confidence Intervals

Uncertainty on key statistics can be quantified via bootstrap resampling:

```{r bootstrap}
bs <- bootstrap_sef(fit, n_boot = 50, verbose = FALSE)
knitr::kable(bs, digits = 4,
             caption = "Bootstrap confidence intervals (50 replicates) for key model statistics.")
```

```{r gg-bootstrap, fig.cap="Bootstrap confidence intervals for PDI, mean entropy, mean energy, log-likelihood, and ARI."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  gg_bootstrap(bs)
}
```

## Harris Matrix Integration

palimpsestr can incorporate stratigraphic information from a Harris Matrix as a penalty during model fitting. The `harris_from_contexts()` function auto-generates a penalty matrix from the mean depth of finds in each context. Alternatively, `read_harris()` imports an external Harris Matrix from a CSV edge list.

```{r harris}
H <- harris_from_contexts(villa_romana, z_col = "z", context_col = "context")
cat("Harris penalty matrix:", nrow(H), "x", ncol(H), "\n")
```

```{r harris-fit}
fit_h <- fit_sef(villa_romana, k = 4, harris = H,
                 tafonomy = "taf_score", context = "context",
                 n_init = 5, seed = 42)
cat("PDI without Harris:", round(pdi(fit), 4), "\n")
cat("PDI with Harris:   ", round(pdi(fit_h), 4), "\n")
di_no <- detect_intrusions(fit)
di_h  <- detect_intrusions(fit_h)
cat("Intrusions without Harris:", sum(di_no$intrusion_prob > 0.5), "\n")
cat("Intrusions with Harris:   ", sum(di_h$intrusion_prob > 0.5), "\n")
```

### Validate phases against stratigraphy

```{r harris-validate}
validate_phases_harris(fit)
```

The `consistent` column flags whether the dominant phase ordering respects stratigraphic depth. `FALSE` indicates an inversion that may warrant investigation.

## Interpretive Report

The `report_sef()` function generates a plain-language summary of the analysis, suitable for inclusion in excavation reports. It is available in English and Italian.

```{r report-en}
report_sef(fit, lang = "en")
```

```{r report-it}
report_sef(fit, lang = "it")
```

## Structured Exports

All results can be exported to CSV for integration with external databases and GIS:

```{r export, eval=FALSE}
export_results(fit, dir = "results/", format = "csv", prefix = "villa_romana")
```

This creates 4 files:

- `villa_romana_phases.csv` --- phase assignments and diagnostics per find
- `villa_romana_intrusions.csv` --- intrusion probabilities
- `villa_romana_us_summary.csv` --- per-US aggregated statistics
- `villa_romana_model_summary.csv` --- model-level metrics (n, k, PDI, BIC, etc.)

## GIS Export

For integration with GIS workflows, palimpsestr can export results as `sf` spatial objects:

```{r gis-export, eval=FALSE}
if (requireNamespace("sf", quietly = TRUE)) {
  sf_pts   <- as_sf_phase(fit, crs = 32632L, dims = "XYZ")
  sf_links <- as_sf_links(fit, quantile_threshold = 0.90, crs = 32632L)
  sf::st_write(sf_pts,   "villa_romana.gpkg", layer = "phases")
  sf::st_write(sf_links, "villa_romana.gpkg", layer = "sei_links")
}
```

## Interactive Visualization

When **plotly** is available, any `gg_*` plot can be converted to an interactive visualization using `as_plotly()`:

```{r interactive, eval=FALSE}
as_plotly(gg_phasefield(fit))
as_plotly(gg_intrusions(fit))
```

## Comparison with Traditional Approaches

| Aspect | Harris Matrix | k-means | palimpsestr (SEF) |
|:-------|:-------------|:--------|:------------------|
| Evidence domains | Stratigraphy only | Spatial only | Spatial + vertical + temporal + cultural |
| Assignment type | Deterministic | Hard | Probabilistic (soft) |
| Mixing detection | No | No | Yes (entropy, ESE) |
| Intrusion detection | Manual | No | Automatic |
| Scalability | Manual | Good | Good |
| Interpretability | High | Low | High (PDI, reports) |

The SEF framework does not replace the Harris Matrix. Rather, it extends it. The Harris Matrix remains indispensable for documenting observed physical relationships between depositional units. What palimpsestr adds is a quantitative, probabilistic layer of analysis that operates on find-level data within those units --- precisely the scale at which palimpsest formation is most problematic.

## Setup recommendations (v0.13.0)

Real-world archaeological datasets often have subtle structural properties
that strongly affect which `fit_sef()` flags should be enabled --- for
example, coordinates assigned at the unit centroid rather than per find,
chronology inherited from the unit's archaeological period, highly
imbalanced class distributions, or informative taphonomic scores.
The `recommend_setup()` helper inspects a dataset and returns the
suggested flags plus diagnostic warnings:

```{r recommend-setup}
rec <- recommend_setup(villa_romana, context = "context", tafonomy = "taf_score")
print(rec)
```

For `villa_romana` the helper detects that (i) coordinates are at unit
centroid, (ii) chronology is US-tied (every find inherits its unit's
period), (iii) taphonomic scores are informative (range 0.30--1.00),
and that the empirically optimal setup is `class_scale = TRUE` plus
`taf_as_feature = TRUE`. The same call also issues an explicit warning
that `detect_intrusions()` direction will be `in_context` for every
row on this dataset --- directional diagnostics require per-find
typological dating, not unit-level chronology.

## Directional intrusion diagnostics (v0.13.0)

```{r directional-intrusions, fig.cap = "Distribution of intrusion directions across stratigraphic units."}
fit_vr <- fit_sef(villa_romana, k = 4, context = "context",
                  tafonomy = "taf_score",
                  class_scale = TRUE, taf_as_feature = TRUE)
intr <- detect_intrusions(fit_vr)
table(intr$direction, useNA = "ifany")
head(intr[intr$direction == "younger_than_context", ])
```

The `direction` column distinguishes **older-than-context** finds (residuals
redeposited from earlier deposits) from **younger-than-context** finds
(latent intrusions, often signalling unrecognised cuts and refills).
The `chrono_gap` column gives the signed mismatch in years.

## Type longevity (v0.13.0)

```{r type-longevity, fig.cap = "Posterior-weighted longevity envelope per cultural class."}
tl <- type_longevity(fit_vr)
head(tl[order(-tl$longevity_span), c("class", "longevity_min", "longevity_max", "longevity_span", "n_finds")])
if (requireNamespace("ggplot2", quietly = TRUE)) gg_longevity(tl, fit_vr)
```

Long-lived classes (e.g. `tegula`, `dolium`) span multiple phases
with non-negligible posterior weight; short-lived classes are confined
to a single phase.

## Calibrated radiocarbon and OxCal import (v0.13.0, v0.19.0)

```{r rcarbon-adapter, eval = requireNamespace("rcarbon", quietly = TRUE)}
cal <- rcarbon::calibrate(x = c(2500, 2400, 2350), errors = c(30, 30, 30), verbose = FALSE)
chrono_df <- chronology_from_rcarbon(cal, method = "hpd")
chrono_df
# Merge into your main data:
# my_data <- merge(my_data, chrono_df, by = "id")
# fit_sef(my_data, ...)
```

The adapter takes calibrated probability densities from
`rcarbon::calibrate()` and returns the chronological-interval columns
(`date_min`, `date_max`, `date_mid`) used by `fit_sef()`. Three
reduction methods are available: HPD (default), median+IQR, and
weighted mean +/- 2 SD.

Since v0.19.0, `chronology_from_oxcal()` provides the same service for
OxCal-modelled dates. It accepts either an `oxcAAR::oxcalCalibrate()` result
(reduced by the same three methods) or a generic data.frame of calibrated
`start`/`end` ranges exported from OxCal:

```{r oxcal-adapter}
ranges <- data.frame(id = c("a", "b"), start = c(-700, -50), end = c(-600, 100))
chronology_from_oxcal(ranges)
```

## Model Limitations and Methodological Notes

### Diagonal Covariance Assumption

The EM uses diagonal Gaussians, assuming conditional independence between features within each phase. This is a deliberate parsimony choice: for $p$ features and $k$ phases, diagonal covariance requires $k(2p+1)-1$ parameters versus $k(p(p+3)/2+1)-1$ for full covariance. With typical archaeological datasets (n = 50--500, p = 7--10), full covariance risks overfitting.

### SEI Component Normalisation

Each SEI component is normalised to [0, 1] by dividing by its within-dataset maximum. This makes weights interpretable but means **absolute SEI values are not comparable across datasets**. Cross-site comparisons should use derived statistics (PDI, mean entropy) which are scale-invariant.

### PDI Interpretation

The PDI is a descriptive measure, not a formal test statistic. The interpretive guideline (PDI > 0.7 = well-resolved) is empirically derived from simulation studies and should not be treated as a formal decision boundary. Bootstrap confidence intervals provide uncertainty quantification.

### Intrusion Detection

The intrusion probability score is a heuristic combining rescaled entropy, ESE, and inverse local SEI. The 0.5 threshold is an exploratory guideline. Users should examine flagged observations individually using domain knowledge.

### Phase Labelling

Phase numbers are arbitrary in mixture models (label switching). Use `reorder_phases()` to ensure phase 1 = deepest (oldest) stratum. The bootstrap procedure applies phase reordering to each replicate.

### Cross-Validation

The `cv_sef()` and `optimize_weights()` functions standardise test data using the training set's center and scale, ensuring consistent feature scales across folds.

### Scalability

The SEI and ESE matrices are $O(n^2)$ in memory. For datasets exceeding ~2000 finds, use `sei_sparse()` or the `max_dist` parameter in `sei_matrix()` to limit computation to spatial neighbours.

## Interactive Shiny Dashboard

palimpsestr includes a built-in Shiny application for interactive analysis. The dashboard supports CSV upload, SQLite and PostgreSQL connections, and provides all diagnostic plots, tables, and maps in a point-and-click interface.

```{r shiny-app, eval=FALSE}
launch_app()
```

The app has six tabs: Data Input (CSV/SQLite/PostgreSQL with column mapping), Analysis (parameter selection and model fitting), Results (interactive plots via plotly), Tables (phase assignments, intrusions, US purity), Maps (overlay on excavation plan geometries), and Report (interpretive report generation and CSV/ZIP download).

## Advanced Feature Space Options

Version 0.12.0 introduces five optional improvements to the feature space that address limitations observed in real-data testing (Poggio Gramignano villa, Solarolo terramare).

### Chronological Precision Weighting

Finds with precise dating (narrow date range) should contribute more to temporal clustering than finds with generic dating. When `chrono_precision = TRUE`, the feature `1/tspan` is added, giving higher weight to precisely dated finds.

```{r eval=FALSE}
# Finds dated "BM2b" (50-year range) weight more than "BM" (300-year range)
fit <- fit_sef(data, k = 3, chrono_precision = TRUE)
```

### Residuality Detection

When a find's dating is inconsistent with its stratigraphic context (e.g., Late Antique ceramic in a Pre-Roman pit fill), the `residuality` feature captures this mismatch. Requires `context`.

```{r eval=FALSE}
# Adds abs(tmid_find - tmid_context_mean) / max_range as a feature
fit <- fit_sef(data, k = 4, context = "context", residuality = TRUE)
```

### Class Scaling

By default, one-hot encoded class columns (0/1) are not scaled, causing them to dominate the clustering when many classes are present. `class_scale = TRUE` multiplies one-hot columns by $1/\sqrt{n_{classes}}$, putting their total energy on par with a single numeric feature.

```{r eval=FALSE}
fit <- fit_sef(data, k = 4, class_scale = TRUE)
```

### Taphonomic Score as Feature

In addition to using `taf_score` as a weight in the EM algorithm, it can be added as an explicit feature dimension, so that finds with similar taphonomic conditions cluster together.

```{r eval=FALSE}
fit <- fit_sef(data, k = 4, tafonomy = "taf_score", taf_as_feature = TRUE)
```

### Sub-class Encoding

When finer typological distinctions are available (e.g., "ceramic_impasto" vs "ceramic_sigillata" instead of just "ceramic"), the `subclass` parameter uses them for one-hot encoding without changing the main `class` column.

```{r eval=FALSE}
fit <- fit_sef(data, k = 4, subclass = "subtype")
```

### Using All Five Together

```{r eval=FALSE}
fit <- fit_sef(data, k = 4,
               tafonomy = "taf_score", context = "context",
               chrono_precision = TRUE,
               taf_as_feature = TRUE,
               residuality = TRUE,
               class_scale = TRUE,
               subclass = "subtype")
```

## Conclusion

The **palimpsestr** package offers a unified probabilistic framework for decomposing archaeological palimpsests, addressing a long-standing gap between qualitative stratigraphic reasoning and quantitative spatial analysis. By jointly modelling spatial proximity, vertical distribution, chronological overlap, and cultural similarity, the Stratigraphic Entanglement Field approach produces soft phase assignments that acknowledge --- rather than suppress --- the inherent uncertainty of mixed deposits. The three diagnostic statistics (SEI, ESE, PDI) provide interpretable, actionable summaries at the find, deposit, and assemblage levels respectively.

The package is available on GitHub at <https://github.com/enzococca/palimpsestr>.

```{r eval=FALSE}
# install.packages("remotes")
remotes::install_github("enzococca/palimpsestr")
```
