## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5.5, dpi = 150, warning = FALSE, message = FALSE) ## ----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") ## ----data-summary------------------------------------------------------------- cat("Material classes:\n") sort(table(villa_romana$class), decreasing = TRUE) ## ----data-structure----------------------------------------------------------- str(villa_romana) ## ----compare-k---------------------------------------------------------------- ck <- compare_k( villa_romana, k_values = 2:7, tafonomy = "taf_score", context = "context", seed = 42 ) print(ck) ## ----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) } ## ----fit-model---------------------------------------------------------------- fit <- fit_sef( villa_romana, k = 4, tafonomy = "taf_score", context = "context", n_init = 10, seed = 42 ) print(fit) ## ----model-summary------------------------------------------------------------ summary(fit) ## ----model-stats-------------------------------------------------------------- cat("PDI:", pdi(fit), "\n") sef_summary(fit) ## ----gg-convergence, fig.cap="EM convergence trace showing log-likelihood stabilisation across iterations."---- if (requireNamespace("ggplot2", quietly = TRUE)) { gg_convergence(fit) } ## ----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) } ## ----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) } ## ----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) } ## ----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) } ## ----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) } ## ----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 ## ----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))) ## ----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.") } ## ----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.") ## ----purity-stats------------------------------------------------------------- n_pure <- sum(us_tab$purity >= 0.9) cat(n_pure, "out of", nrow(us_tab), "units have purity >= 90%\n") ## ----transition-matrix-------------------------------------------------------- tm <- phase_transition_matrix(fit) print(tm) ## ----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) ## ----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) } ## ----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.") ## ----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.") ## ----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------------------------------------------------------------------- H <- harris_from_contexts(villa_romana, z_col = "z", context_col = "context") cat("Harris penalty matrix:", nrow(H), "x", ncol(H), "\n") ## ----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") ## ----harris-validate---------------------------------------------------------- validate_phases_harris(fit) ## ----report-en---------------------------------------------------------------- report_sef(fit, lang = "en") ## ----report-it---------------------------------------------------------------- report_sef(fit, lang = "it") ## ----export, eval=FALSE------------------------------------------------------- # export_results(fit, dir = "results/", format = "csv", prefix = "villa_romana") ## ----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, eval=FALSE-------------------------------------------------- # as_plotly(gg_phasefield(fit)) # as_plotly(gg_intrusions(fit)) ## ----recommend-setup---------------------------------------------------------- rec <- recommend_setup(villa_romana, context = "context", tafonomy = "taf_score") print(rec) ## ----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", ]) ## ----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) ## ----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, ...) ## ----oxcal-adapter------------------------------------------------------------ ranges <- data.frame(id = c("a", "b"), start = c(-700, -50), end = c(-600, 100)) chronology_from_oxcal(ranges) ## ----shiny-app, eval=FALSE---------------------------------------------------- # launch_app() ## ----eval=FALSE--------------------------------------------------------------- # # Finds dated "BM2b" (50-year range) weight more than "BM" (300-year range) # fit <- fit_sef(data, k = 3, chrono_precision = TRUE) ## ----eval=FALSE--------------------------------------------------------------- # # Adds abs(tmid_find - tmid_context_mean) / max_range as a feature # fit <- fit_sef(data, k = 4, context = "context", residuality = TRUE) ## ----eval=FALSE--------------------------------------------------------------- # fit <- fit_sef(data, k = 4, class_scale = TRUE) ## ----eval=FALSE--------------------------------------------------------------- # fit <- fit_sef(data, k = 4, tafonomy = "taf_score", taf_as_feature = TRUE) ## ----eval=FALSE--------------------------------------------------------------- # fit <- fit_sef(data, k = 4, subclass = "subtype") ## ----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") ## ----eval=FALSE--------------------------------------------------------------- # # install.packages("remotes") # remotes::install_github("enzococca/palimpsestr")