bic + 2 * sum(entropy) (was
bic - 2 * sum(entropy)). With the lower-is-better BIC
convention the entropy term must be added so that
overlapping/fuzzy partitions are penalised (Biernacki, Celeux &
Govaert 2000); the previous sign rewarded them. compare_k()
ICL values change accordingly (BIC, PDI and the case study are
unaffected; ICL is not used in any manuscript figure).archaeo_sim(): the simulated date bounds
date_min and date_max are now returned as
integer-valued years (via floor()), consistent with
archaeological dating conventions (reviewer suggestion, PCI Archaeo
#1019).harris_from_contexts() gains
exclude_contexts: a character vector of context labels
exempt from the verticality (depth-rank) penalty, for the fills of cuts
and multi-period contexts where depth does not track chronology. Default
NULL leaves the behaviour unchanged. The help page now
documents the verticality assumption and its failure modes (reviewer
suggestion, PCI Archaeo #1019).compare_k() and pdi() help
pages now explain the BIC and PDI model-selection metrics and reference
Schwarz (1978); the manuscript adds the same definitions and stresses
combining statistical criteria with archaeological judgement.read_pyarchinit() gains reperti_table
(default "pyarchinit_reperti") and
reperti_geometry: when a material find is plotted as a
point in pyarchinit_reperti (matched by site +
numero_inventario = id_rep, the
pyarchinit_reperti_view join), its own point
x, y (and z from
quota or the geometry Z) are used instead of the US-polygon
centroid. Finds without a point fall back to the US centroid as before.
Read via sf/GDAL so SpatiaLite and PostGIS geometries are parsed
correctly. This enables find-level spatial resolution where
piece-plotting is available.read_pyarchinit(): when taf is
NULL (the default) and the chronology_table
(palimpsest_chronology) carries a per-US taf
column, those values now populate taf_score (units without
a value default to 0.5). This lets the pyArchInit QGIS algorithms (fit /
intrusions / report) honour a per-US taphonomic score stored alongside
the absolute chronology, without any change to the calling
.rsx.read_pyarchinit() gains chronology_table
(default "palimpsest_chronology"): an optional per-US table
of absolute dating with columns sito, area,
us, start, end (calendar years,
BCE negative) — e.g. OxCal-calibrated ranges from
chronology_from_oxcal(). When present it overrides the
free-text datazione for the matching units (envelope of
multiple dates per US), joined on (sito, area, us) with a
(sito, us) fallback.r:palimpsestrfit,
r:palimpsestrintrusions, r:palimpsestrreport)
accept an optional PG_connection parameter (a libpq DSN)
and read from PostgreSQL/PostGIS when it is set, falling back to the
SQLite/Spatialite Database_file otherwise.
read_pyarchinit() itself was already database-agnostic
(plain DBI). The pyArchInit tab passes its active connection string, so
PostgreSQL pyArchInit projects are supported.export_sef_report(): assembles a complete, narrated
report from a sef_fit — the report_sef()
interpretive text, all applicable gg_* diagnostic plots,
and diagnostic tables (US summary, phase-transition matrix, per-US phase
assignments, model statistics) — rendered to PDF and/or DOCX from a
shipped RMarkdown template (inst/rmarkdown/sef_report.Rmd),
in Italian or English. A markdown narrative sidecar
(<file>.md) is always written; when pandoc/LaTeX are
unavailable the function degrades gracefully to that markdown plus PNG
figures. tinytex added to Suggests.read_pyarchinit() now sources elevation per
archaeological practice: US elevation from pyarchinit_quote
(mean quota_q per unit, with a (sito, us) join
fallback), material elevation from quota_usm, and pottery
inheriting its US elevation. The deprecated
us_table.quota_abs form field is no longer used.source = "both"|"materials"|"pottery" argument:
reads inventario_materiali_table and/or
pottery_table (pottery class taken from the
first non-empty of
ware/material/form); a
find_source column is added when
source = "both".r:palimpsestrreport algorithm
(palimpsestr_report_db.rsx): database → fit → narrated
PDF/DOCX report, with Source / Language / Format parameters.r:palimpsestrfit and
r:palimpsestrintrusions gain a Source
parameter.villa_romana example database builder now populates
quota_usm, pottery_table, and
pyarchinit_quote to exercise the corrected sourcing.read_pyarchinit(): reads the
inventario_materiali_table (finds) and
us_table (stratigraphic units) of a pyArchInit database
(SQLite/Spatialite or PostgreSQL/PostGIS) and assembles a
fit_sef()-ready data.frame. Each find inherits its
stratigraphic unit’s coordinates (centroid of the US polygon supplied
via us_geometry, e.g. the
pyunitastratigrafiche layer), elevation, chronology, and
context. Free-text period strings (e.g. “II sec. a.C.”, “I/II
sec. d.C.”, “età romana”) are resolved to numeric date_min/
date_max by a built-in archaeological-date parser
(extendable via date_labels); units without a digitised
polygon fall back to synthetic per-unit coordinates so the analysis can
still run at unit resolution. Validated against real pyArchInit
Spatialite databases.chronology_from_oxcal(): OxCal counterpart of
chronology_from_rcarbon(). Accepts either an
oxcAAR::oxcalCalibrate() result
(oxcAARCalibratedDatesList, reduced by HPD / median+IQR /
weighted mean) or a generic data.frame of
start/end calibrated ranges, and returns the
date_min/date_max/date_mid
columns expected by fit_sef(). oxcAAR is in
Suggests.detect_intrusions() gains an
intrusion_threshold argument and an
intrusion_type factor column that couples the intrusion
magnitude (the noise-component posterior) with the chronological
direction: not_flagged, residual (flagged and
older-than-context), latent_feature (flagged and
younger-than-context), or outlier_in_context. Existing
columns are unchanged.gg_outliers() now colours each find by
intrusion_type.gg_phase_composition(): per-phase class profile from
the multinomial model (cat_prob), as a heatmap (default) or
stacked bars.gg_direction(): diverging lollipop of the directional
intrusion diagnostic (residual / older-than-context vs latent-feature /
younger-than-context), with a graceful placeholder when every find is in
context.gg_outliers(): non-spatial ranking of the intrusion
probability, using the noise-component posterior when available,
otherwise the heuristic composite.gg_unit_coherence(): within-stratigraphic-unit phase
coherence, with an optional comparison fit (e.g. the legacy
Gaussian-class model).All four are exposed in a new “Model profile” tab of the Shiny app
and are compatible with as_plotly().
Bug fixes for noise-component model selection (addresses code review on the 0.17.0 pull request).
cv_sef() and optimize_weights() scored only
the Gaussian phases of a noise = TRUE fit, so the uniform
component was omitted and outlier held-out points were mis-scored (their
mass was missing and the Gaussian weights no longer summed to one). The
fit now stores the noise mixing weight (noise_weight), its
log-density (noise_logdens), and the global class
frequencies (cat_global), and both evaluators add the noise
component to the held-out likelihood.compare_k(..., noise = TRUE) no longer
gives noise models a free parameter relative to non-noise fits.New
fit_sef(noise = FALSE, noise_prior = 0.05). When
TRUE, a uniform background component is added to the
mixture (Fraley & Raftery, 1998). Finds that fit no Gaussian phase
accumulate posterior probability on it, so the fit gains a
noise_prob vector that is a genuine posterior
probability of being an outlier/intrusion — unlike the previous
composite score, it is not forced onto [0, 1] by rescaling,
so a clean assemblage need not contain a probability-1 intrusion. As a
side effect the noise component absorbs extreme finds, keeping them out
of the phase centroid and variance estimates (robust phase estimation).
The reported phase_prob is then the distribution over
phases conditional on a find not being noise.
detect_intrusions() now returns the
model-based posterior in its intrusion_prob column
whenever the fit was made with noise = TRUE, falling back
to the heuristic composite (rescaled entropy + energy + inverse local
SEI) otherwise. The directional columns are unchanged.
The option is propagated through bootstrap_sef() and
exposed in the Shiny app.
This release adds two model-based treatments of evidence that the EM engine previously simplified away: per-find dating uncertainty, and the stratigraphic constraint. Both are opt-in and off by default, so existing fits are unchanged.
fit_sef(chrono_uncertainty = FALSE). When
TRUE, a find’s dating interval width is propagated into the
likelihood instead of collapsing the interval to its mid-point. The
mid-date is modelled as observed with a measurement variance
(date_max - date_min)^2 / 12 (a uniform prior over the
interval), added to the temporal dimension’s component variance in the
E-step and removed from the component variance estimate by
method-of-moments deconvolution in the M-step. Finds with wide dating
ranges then rely less on chronology and carry their dating uncertainty
into the phase probabilities. Unlike chrono_precision
(which adds 1/tspan as a feature), this is a model-based
treatment of the uncertainty. The internal
diag_log_density() gains an extra_var
argument.New
fit_sef(strat_dynamic = FALSE, strat_beta = 1).
The stratigraphic context constraint was previously a static penalty
frozen from the k-means initialisation: if the EM moved assignments, the
penalty still reflected the starting clusters. With
strat_dynamic = TRUE it becomes a Neighborhood-EM /
hidden-Markov-random-field term recomputed from the current phase
posteriors at every E-step (Ambroise & Govaert, 1997): each find is
rewarded for sharing the phase of its stratigraphic-unit neighbours, so
units stay coherent as the fit evolves. strat_beta sets the
field strength; any harris penalty continues to be applied
statically.
Both options are propagated through bootstrap_sef()
and exposed in the Shiny app.
The class label is now modelled as a per-phase
categorical distribution (a Gaussian-times-multinomial
mixed-type mixture), controlled by the new class_model
argument to fit_sef(). "multinomial" is the
default; "gaussian" reproduces the previous one-hot
behaviour and is kept for backward compatibility.
Up to v0.14.0 the class label was one-hot encoded into the Gaussian feature block. Because zero/one dummies make finds of different classes almost perfectly separable, this had two undesirable effects: per-component variances on the dummy columns collapsed to the floor (over-confident assignments), and — more importantly — the class label could split a single stratigraphic unit across several phases, even though a unit is one depositional event. Under the multinomial model the numeric evidence (space, depth, chronology) drives the phase assignment and the class label only tilts it softly, so units with finds of mixed classes stay coherent.
New class_smoothing argument
(default 1): a Dirichlet pseudo-count, shrunk towards the global class
frequencies, that keeps every per-phase class probability strictly
positive. 0 removes the shrinkage.
New phase_composition() accessor:
returns the estimated per-phase class profile (row j =
probability that a find in phase j belongs to each class).
The model-based, directly interpretable counterpart of a phase-by-class
cross-tabulation, stored on the fit as $cat_prob.
cv_sef(), optimize_weights(), and
bootstrap_sef() propagate the class model and include the
categorical term in the held-out likelihood;
reorder_phases() reorders the per-phase class profiles;
BIC/ICL use the correct mixed-type parameter count (k(L-1)
categorical parameters instead of treating the dummies as
Gaussian).
Note: on datasets recorded at stratigraphic-unit resolution (finds
inheriting US-centroid coordinates and US-tied chronology, e.g. the
bundled villa_romana), phase assignments remain
near-certain under either model because there are only as many distinct
numeric points as there are units; that certainty reflects the data
resolution, not the model. The multinomial model’s benefit there is
keeping units coherent and yielding interpretable per-phase class
profiles.
This release is a statistical-correctness pass on the core engine.
Several parameters that were silently inert are now active, and three
diagnostics that depended on the measurement scale or were distorted by
single outliers have been made scale-invariant and robust. Fits produced
with non-default weights, and the output of
optimize_weights(), will differ from 0.13.0; default-weight
fits are affected only by the more careful BIC/ICL accounting and the
new multi-start default.
Domain weights now enter the model.
weights = c(ws, wz, wt, wc) previously affected only the
SEI matrix, leaving fit_sef() — and hence the
cross-validated objective of optimize_weights() —
completely invariant to them. The weights now scale the corresponding
standardised feature dimensions (ws: planar coordinates,
wz: depth, wt: chronology-derived features,
wc: the one-hot class block), so they genuinely influence
the mixture likelihood. Weights must be strictly positive; they are
stored on the fit and propagated by bootstrap_sef(),
cv_sef(), and optimize_weights().
New var_structure argument to
fit_sef() ("diagonal", default, or
"spherical"). Free diagonal variances absorb any
per-dimension rescaling, so domain weights are not identifiable under a
diagonal fit; the "spherical" structure (one shared
variance per component) makes them identifiable.
optimize_weights() now selects weights under
var_structure = "spherical" and corrects each held-out
log-likelihood by the Jacobian of the weighting
(n_test * sum(log(w_d))) so configurations with different
weights are directly comparable.
Taphonomic down-weighting is no longer applied
twice. fit_sef() passed
weights_obs = 1 - 0.5 * taf and the EM multiplied
by (1 - 0.5 * taf) again, so a find with
taf = 1 contributed at weight 0.25 instead of the
documented 0.5. It is now applied exactly once.
Removed the inert taf_weight_e
term. It subtracted a per-row constant from every mixture
component, which cancels exactly in the softmax and so had no effect on
the posteriors (only on the reported log-likelihood). The dead
taf, taf_weight_m, and
taf_weight_e arguments were removed from the internal
EM.
BIC and ICL are now built from the true unpenalized
mixture likelihood. With a stratigraphic penalty the EM
objective is a penalized criterion, not a log-likelihood; using it
biased comparisons across K and across context settings.
model_stats$loglik is the unpenalized likelihood; the
penalized EM objective is available as
model_stats$loglik_penalized. The "spherical"
structure also gets its own (smaller) parameter count in the
BIC.
SEI and ESE spatial/vertical components are now
scale-invariant and robust to duplicate coordinates. The former
1/d kernel normalised by its maximum was dominated by the
single closest pair, so one near-duplicate pair (e.g. two finds sharing
a grid square) collapsed every other spatial affinity towards zero. Both
components now use a bounded exponential kernel exp(-d / h)
with a data-driven bandwidth h (median pairwise
separation). ese() consequently no longer depends on the
measurement unit of the site grid. The eps and
z_floor arguments of sei_matrix() are
deprecated and ignored.
detect_intrusions() directional envelope is
now robust. The leave-one-out unit envelope used a strict
min/max over the other finds, so a single chronological outlier in a
unit could mask genuinely residual/intrusive neighbours. A new
envelope argument (default c(0.05, 0.95)) uses
robust quantiles; pass c(0, 1) for the previous min/max
behaviour.
n_init now defaults to 5 (was 1)
and k-means initialisation uses nstart = 5. The EM
objective is multimodal, so a single start risked a poor local
optimum.
detect_intrusions() now returns two new columns:
direction: factor classifying each find as
older_than_context, in_context, or
younger_than_context based on interval overlap between the
find’s chronology and the leave-one-out chronological envelope of its
stratigraphic unit.chrono_gap: signed offset in years (negative =
residual, positive = latent intrusion, 0 = in-context). Backward
compatible: the existing intrusion_prob column is
unchanged.New type_longevity() function: per-class temporal
envelope of the phases with mean posterior weight above a configurable
threshold (default 0.1). Returns longevity_min,
longevity_max, longevity_span,
dominant_phase, n_finds, and a list-column
weight_matrix with the per-phase posterior
weights.
New gg_longevity() Gantt-style plot for the
type_longevity() output.
New chronology_from_rcarbon() adapter: converts an
rcarbon::CalDates object into the
date_min/date_max/date_mid
columns expected by fit_sef(). Three reduction methods: HPD
(default), median_iqr, weighted_mean. BCE/CE
sign convention by default. rcarbon is in
Suggests.
New recommend_setup() helper: inspects a dataset and
returns recommended fit_sef() flags
(class_scale, taf_as_feature,
chrono_precision, residuality) based on
diagnostic heuristics (number of classes, singleton classes, US-centroid
coordinates, US-tied chronology, taf_score informativeness). Reports
caveats (e.g. directional intrusions will all be in_context
when the chronology is US-tied) and prints a ready-to-paste
fit_sef() call. Closes #2.
direction and
chrono_gap.No breaking changes. All v0.12.0 user code continues to work as before.
gg_cv() (cross-validation diagnostics),
gg_bootstrap() (confidence interval forest plot),
gg_weights() (weight sensitivity heatmap).villa_romana dataset: realistic 300-find Roman villa
with 4 phases, bioturbation, construction cuts, and residual
pottery.cv_sef(),
optimize_weights()).reorder_phases()
applied to each replicate before ARI computation, solving label
switching in confidence intervals.neighbourhood is set, energy is divided by actual neighbour
count instead of n-1, removing peripheral bias.feature_matrix() now
accepts external center/scale parameters and
preserves scaling attributes for proper cross-validation.adjusted_rand_index() and
confusion_matrix() for comparing estimated vs true
phases.us_summary_table() aggregates
diagnostics per stratigraphic unit;
phase_transition_matrix() reveals vertical phase ordering;
export_results() writes all results to CSV.harris_from_contexts()
auto-generates penalties from depth ordering; read_harris()
imports CSV edge lists; validate_phases_harris() checks
phase-stratigraphy consistency.gg_convergence() (EM trace),
gg_phase_profile() (depth vs phase),
gg_confusion() (confusion matrix heatmap with ARI).bootstrap_sef() for bootstrap confidence intervals on
PDI, entropy, energy, loglik, and ARI.cv_sef() for k-fold cross-validation of phase count
selection.optimize_weights() for data-driven SEI weight
estimation via grid search cross-validation.reorder_phases() relabels phases by mean depth (phase 1
= deepest/oldest) to solve label switching.sei_sparse() wrapper for sparse SEI computation on
large datasets (n > 1000).fit_sef() gains n_init parameter for
multiple random initialisations (default: 1). Best run by log-likelihood
is retained.fit_sef() default em_iter increased from
25 to 100 for better convergence.$converged flag in
sef_fit objects; warning issued on non-convergence.fit_sef() and compare_k() model selection
metrics.sei_matrix() gains max_dist parameter for
sparse computation (zeroes beyond threshold).sei_matrix() and ese() fully vectorized
(50-100x faster on large datasets).as_plotly() with
enriched archaeological tooltips (ID, context, phase, dating, class,
entropy, energy, intrusion probability).gg_* functions now embed tooltip data for seamless
plotly conversion.@examples,
@family, and @seealso for all exported
functions._pkgdown.yml for documentation website.cran-comments.md for CRAN submission.demo_easy,
demo_moderate, demo_compressed.summary.sef_fit() and
print.sef_fit() (NAMESPACE now uses S3method()
instead of export()).as_sf_phase() and as_sf_links() CRS
default (NA → NA_integer_) for compatibility
with sf >= 1.0.Authors@R in DESCRIPTION.LazyData field from DESCRIPTION.fit_sef().as_sf_links() and enhanced
as_sf_phase() with XY/XYZ support.
Need a high-speed mirror for your open-source project?
Contact our mirror admin team at info@clientvps.com.
This archive is provided as a free public service to the community.
Proudly supported by infrastructure from VPSPulse , RxServers , BuyNumber , UnitVPS , OffshoreName and secure payment technology by ArionPay.