Package {leafwax}

Title:	Bayesian Inversion of Leaf Wax Hydrogen Isotopes to Precipitation
Version:	0.2.7
Description:	Bayesian inversion of leaf wax hydrogen isotopes to reconstruct precipitation isotopes using hierarchical spatial models. Provides fourteen Bayesian models that vary in their use of spatial Gaussian processes and ancillary covariates (precipitation amount, plant functional type, C4 fraction). Models are pre-computed using 'Stan' and stored as posterior distributions, so prediction does not require 'Stan' to be installed. A 100-draw fixture ships with the package; full 1000-draw posteriors are downloaded from a versioned 'Zenodo' deposit on first use; see Bradley (2026) <doi:10.5281/zenodo.20085465>.
License:	MIT + file LICENSE
URL:	https://github.com/bradleylab/leafwax, https://bradleylab.github.io/leafwax/
BugReports:	https://github.com/bradleylab/leafwax/issues
Encoding:	UTF-8
Language:	en-US
RoxygenNote:	7.3.3
Imports:	jsonlite, stats, utils
Suggests:	rappdirs, knitr, rmarkdown, testthat (≥ 3.0.0)
VignetteBuilder:	knitr
Depends:	R (≥ 3.5.0)
LazyData:	true
LazyDataCompression:	xz
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2026-06-01 17:49:57 UTC; abradley
Author:	Alex Bradley [aut, cre]
Maintainer:	Alex Bradley <abradley@wustl.edu>
Repository:	CRAN
Date/Publication:	2026-06-15 13:40:02 UTC

leafwax: Bayesian Calibration of Leaf Wax Hydrogen Isotope Reconstructions

Description

The leafwax package provides tools for probabilistic inversion of leaf wax hydrogen isotope measurements (delta-2-H) to reconstruct precipitation isotope values. It implements hierarchical Bayesian models that account for multiple sources of uncertainty including measurement error, biological fractionation, and spatial correlation in isotope patterns.

Main Functions

invert_d2H

Bayesian inversion of leaf wax delta2H to precipitation delta2H

available_models

List all available calibration models

load_posteriors

Load posterior distributions for a specific model

get_model_parameters

Get model capabilities and required parameters

validate_model_inputs

Validate inputs for a specific model

Available Models

The package includes 14 calibration models with different capabilities. The v10 fits include precipitation amount (baseline_env* and full* variants), C4 abundance, and PFT cover; none of the v10 variants carry a fitted elevation coefficient despite the historical "elevation_*" naming. Runtime capability flags in load_posteriors() are derived from each model's posterior columns at load time.

• Basic models: baseline, baseline_sp

• Precipitation models: baseline_env, baseline_env_sp

• Vegetation models: baseline_veg, baseline_veg_sp, c4_only_sp

• Combined spatial models: elevation_only_sp, elevation_c4_sp, elevation_c4_interact_sp

• Full models: full, full_sp, full_interact, full_interact_sp

Models with "_sp" suffix use spatial Gaussian processes with 125 knots on a Fibonacci sphere lattice for improved uncertainty quantification.

Model Selection

Pass model = "auto" to predict_d2h_precip() to let select_best_model_from_flags() choose a model based on which covariates the caller has supplied; otherwise pick a model name from available_models() explicitly.

Key Features

• Hierarchical Bayesian framework for uncertainty propagation

• Support for single and multi-location inversions

• Spatial correlation via Gaussian processes

• Automatic handling of missing covariates

Author(s)

Maintainer: Alex Bradley abradley@wustl.edu (ORCID)

References

Bowen, G. J., Cai, Z., Fiorella, R. P., & Putman, A. L. (2019). Isotopes in the water cycle: Regional-to global-scale patterns and applications. Annual Review of Earth and Planetary Sciences, 47, 453-479. doi:10.1146/annurev-earth-053018-060220

Sachse, D., Billault, I., Bowen, G. J., Chikaraishi, Y., Dawson, T. E., Feakins, S. J., ... & Kahmen, A. (2012). Molecular paleohydrology: Interpreting the hydrogen-isotopic composition of lipid biomarkers from photosynthesizing organisms. Annual Review of Earth and Planetary Sciences, 40, 221-249. doi:10.1146/annurev-earth-042711-105535

Bradley, A. (2026). leafwax v10 model posteriors. Zenodo DOI doi:10.5281/zenodo.20085465.

See Also

Useful links:

• https://github.com/bradleylab/leafwax

• https://bradleylab.github.io/leafwax/

• Report bugs at https://github.com/bradleylab/leafwax/issues

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  # List available models
  models <- available_models()
  n_models <- length(models)

  # Simple single-location inversion
  result <- invert_d2H(
    d2H_wax = -150,
    d2H_wax_sd = 3,
    longitude = -120,
    latitude = 40,
    model_name = "baseline",
    verbose = FALSE
  )
})

Assess a paleoclimate claim against the leaf-wax taxonomy

Description

Walks the four-level taxonomy from manuscript Section 4.5.6 and reports the highest level a claim survives at. The taxonomy is:

• Level 1: a leaf-wax delta-2-H change occurred between two intervals. Defensible when the change exceeds analytical uncertainty.

• Level 2: the wax change is consistent with a directional hydroclimate change. Requires (i) sediment-source change ruled out by independent evidence (sediment_source_ruled_out), AND (ii) depositional artifact ruled out by independent evidence (depositional_artifact_ruled_out), AND EITHER (a) named corroborating evidence against vegetation reorganization via corroborating_proxies (the original path), OR (b) demonstration that the observed wax shift exceeds the vegetation-only envelope computed from a user-supplied PFT-change scenario (level2_vegetation_path; see compute_vegetation_envelope() and manuscript Section 4.5.3).

• Level 3: the wax change implies a quantitative delta-2-H_precip magnitude. Requires a defended local effective slope and explicit uncertainty propagation through the inversion. When reconstruction is NULL the function calls invert_d2H() itself.

• Level 4: the magnitude is uniquely attributable to precipitation isotope change rather than to vegetation, source-water seasonality, or evapotranspirative enrichment. Requires independent stationarity evidence for each non- precipitation control over the interval.

Usage

assess_claim(
  record,
  claim,
  reconstruction = NULL,
  longitude = NULL,
  latitude = NULL,
  model_name = "baseline_sp",
  ...
)

Arguments

Details

Use this when a colleague claims that a downcore record shows a specific d2H_precip shift and you want a structured check of which levels of the taxonomy that claim actually clears.

Value

A list with elements:

• highest_level - integer in 0:4. 0 means even Level 1 did not clear.

• levels - data frame, one row per level, with columns level, passed (logical), summary (one-line reason).

• details - per-level lists of computed quantities (e.g., delta_wax, threshold, p_exceed, missing fields).

• claim - the (validated) claim object.

Get available models

Description

Returns a list of all available calibration models for leaf wax hydrogen isotope inversion. Models vary in their complexity and data requirements.

Usage

available_models()

Value

Character vector of available model names. Models include:

• baseline: Basic OIPC model without spatial effects

• baseline_sp: Basic model with spatial Gaussian process

• baseline_env: Includes precipitation-amount effects

• baseline_env_sp: Precipitation-amount effects with spatial GP

• baseline_veg: Includes vegetation interaction effects

• baseline_veg_sp: Vegetation interactions with spatial GP

• c4_only_sp: C4 vegetation effects only (spatial)

• elevation_only_sp: Historical elevation-context variant (spatial)

• elevation_c4_sp: Historical elevation-context + C4 variant

• elevation_c4_interact_sp: Historical elevation/C4-interaction name; C4 effect only

• full: Precipitation amount + vegetation interactions without spatial component

• full_sp: Precipitation amount + vegetation interactions with spatial component

• full_interact: Precipitation amount + vegetation interactions

• full_interact_sp: Full interaction model with spatial GP

Examples

# List all available models
models <- available_models()
print(models)

# Get details for a specific model
model_info <- get_model_parameters("baseline_sp")
print(model_info$description)

Batch predict precipitation d2H for multiple sites

Description

Processes multiple sites with progress indicators and optional parallelization. Handles large datasets efficiently by processing in chunks.

Usage

batch_predict(
  data,
  model = "auto",
  chunk_size = 100,
  parallel = FALSE,
  n_cores = NULL,
  progress = TRUE,
  return_diagnostics = FALSE,
  ...
)

Arguments

Value

Data frame with predictions for all sites

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  data(example_data)
  large_data <- example_data[rep(seq_len(nrow(example_data)), length.out = 12), ]
  row.names(large_data) <- NULL

  # Process in chunks
  results <- batch_predict(
    large_data,
    chunk_size = 6,
    progress = FALSE,
    verbose = FALSE
  )

  # Process with a specific model
  results <- batch_predict(
    large_data,
    model = "baseline_sp",
    chunk_size = 6,
    progress = FALSE,
    verbose = FALSE
  )
})

Check if model data exists in cache

Description

Checks whether the heavy posterior file for a model is present in the user cache populated by download_model_data().

Usage

check_data_cache(
  model_name,
  data_type = c("minimal", "standard", "full"),
  verbose = FALSE
)

Arguments

Value

Logical indicating whether the cached posterior file exists.

Examples

local({
  old <- options(leafwax.cache_dir = file.path(tempdir(), "leafwax_cache"))
  on.exit(options(old))

  exists <- check_data_cache("baseline_sp", verbose = FALSE)
})

Clear download cache

Description

Removes downloaded model data from the local cache.

Usage

clear_download_cache(
  model_name = NULL,
  type = c("all", "posteriors"),
  confirm = TRUE
)

Arguments

Value

Invisible NULL

Examples

local({
  old <- options(leafwax.cache_dir = file.path(tempdir(), "leafwax_cache"))
  on.exit({
    unlink(getOption("leafwax.cache_dir"), recursive = TRUE, force = TRUE)
    options(old)
  })

  post_dir <- file.path(getOption("leafwax.cache_dir"), "posteriors")
  dir.create(post_dir, recursive = TRUE, showWarnings = FALSE)
  file.create(file.path(post_dir, "baseline_sp_posterior.rds"))

  # Clear cache for a specific model without prompting
  suppressMessages(clear_download_cache("baseline_sp", confirm = FALSE))
})

Compare predictions across multiple models

Description

Runs predictions using multiple models and compares results.

Usage

compare_models(
  data,
  models = NULL,
  summary_fun = mean,
  return_all = FALSE,
  progress = TRUE,
  ...
)

Arguments

Value

Data frame with ensemble predictions or list of all results

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  data(example_data)

  # Compare multiple models
  comparison <- compare_models(
    example_data,
    models = c("baseline", "baseline_env", "baseline_sp"),
    progress = FALSE
  )

  # Get all individual model results
  all_results <- compare_models(
    example_data,
    models = c("baseline", "baseline_env"),
    return_all = TRUE,
    progress = FALSE
  )
})

Vegetation-only envelope for a paleo wax-isotope record

Description

Computes the posterior-propagated wax-shift envelope expected under a user-supplied PFT-change scenario, holding d2H_precip constant. This is the magnitude path of the Level 2 claim taxonomy described in the accompanying manuscript Section 4.5.3: the calibration's PFT main-effect and PFT-by-\delta^2 H_p interaction coefficients are combined across all posterior draws to bound how much wax change vegetation reorganization alone can produce at the site, with no contribution from precipitation-isotope change.

Usage

compute_vegetation_envelope(
  oipc_ref,
  from,
  to,
  model_name = "full_interact_sp",
  n_draws = NULL,
  verbose = TRUE
)

Arguments

Details

For each posterior draw the per-draw envelope is ⁠envelope = sum_k beta_k * delta_pft_k + sum_k beta_d2Hp_x_k * oipc_ref * delta_pft_k⁠ where k indexes the four PFT classes (tree, shrub, grass, C4) and delta_pft_k = to[k] - from[k]. The \beta_{\delta^2 H_p} \times PFT interaction terms are interactions with the precipitation-isotope calibration slope. The returned envelope_p975_abs = quantile(abs(envelope_draws), 0.975) is the absolute upper bound used by assess_claim() path (b) to test whether an observed ⁠|delta_wax|⁠ exceeds what the supplied PFT scenario alone can produce.

Passing path (b) rejects the vegetation-only null for the supplied scenario at the site. It does not identify the hydroclimate mechanism, quantify the precipitation-isotope change, or address sediment-source change, depositional artifact, compound-source mixing, age-model errors, evapotranspirative regime change, or seasonality shifts (which assess_claim() gates with separate fields). The calibration coefficients are derived from spatial variation across sites; applying them to within-record temporal vegetation change assumes the same response holds through time at one location.

Value

List with

• envelope_draws - numeric vector of per-draw signed envelopes (per mil).

• envelope_median - posterior median of envelope_draws.

• envelope_p975_abs - quantile(abs(envelope_draws), 0.975), the absolute upper bound used by assess_claim() path (b).

• oipc_ref - the oipc_ref value used (echoed for traceability).

• delta_pft - named numeric vector to - from.

• n_draws_used - integer, number of draws contributing to envelope_draws.

• model_name - the model name used.

• details - list with coefs_summary (posterior median of each of the eight coefficients).

Manuscript reference

Section 4.5.3 of the accompanying manuscript defines the vegetation-only envelope and the constant-precipitation framing. Section 4.5.6 places it in the Level 2 claim taxonomy.

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  # Hypothetical: a 30 percentage-point woody-to-grass transition at a
  # site where the OIPC raster lookup gave d2H_precip approximately -60 per mil.
  env <- compute_vegetation_envelope(
    oipc_ref = -60,
    from = c(tree = 0.4, shrub = 0.3, grass = 0.2, C4 = 0.05),
    to   = c(tree = 0.1, shrub = 0.2, grass = 0.5, C4 = 0.20),
    model_name = "full_interact_sp",
    n_draws = 100,
    verbose = FALSE
  )
  vegetation_bound <- env$envelope_p975_abs
})

Within-record d2H_precip change detection

Description

Given a downcore invert_d2H() reconstruction posterior with full draws (return_full = TRUE), report (a) the posterior probability that the difference in mean d2H_precip between two stratigraphic intervals exceeds user-supplied magnitudes, and (b) the within- record 95\ 4.5.3:

Usage

detect_change(
  reconstruction,
  age,
  baseline_interval,
  test_intervals = NULL,
  sigma_residual,
  sigma_analytical = 3,
  rho_t = NULL,
  beta_eff,
  confidence = 0.95,
  magnitudes = NULL
)

Arguments

Details

\mathrm{threshold}_{precip} = \frac{z_{\alpha/2}\, \sqrt{2 \sigma_{residual}^2 (1 - \rho_t) + 2 \sigma_{analytical}^2}} {\beta_{\mathrm{eff}}}

The threshold is the smallest difference in d2H_precip between two independent samples that can be distinguished from within-record noise at the chosen confidence level. The lag-1 autocorrelation rho_t enters only the residual term, because analytical measurement error is independent between samples by construction. The spatial GP intercept contributes a constant to every sample in the record and cancels in the contrast; the same sigma_residual from the spatial calibration applies (manuscript Section 4.5.3).

Value

A list with elements:

• threshold - the detection threshold on d2H_precip at the requested confidence level.

• formula - the components used: z, rho_t, sigma_residual, sigma_analytical, beta_eff.

• intervals - a data frame with one row per test interval reporting the posterior median and CI of delta and (if magnitudes supplied) the posterior probability of exceeding each magnitude.

Detect model capabilities from static v10 metadata

Description

Detect model capabilities from static v10 metadata

Usage

detect_model_capabilities(model_name)

Arguments

Value

List of capability flags

Download model data from GitHub releases

Description

Downloads model posterior draws and lookup tables from GitHub releases with progress tracking and integrity verification.

Usage

download_model_data(
  model_name,
  version = "latest",
  data_type = c("posteriors"),
  cache_dir = NULL,
  overwrite = FALSE,
  verify = TRUE,
  verbose = TRUE
)

Arguments

Value

Logical indicating success

Examples

cache_dir <- file.path(tempdir(), "leafwax_download_example")
ok <- download_model_data(
  "baseline",
  version = "v1.0.1",
  cache_dir = cache_dir,
  verify = FALSE,
  verbose = FALSE
)

Download file with progress bar

Description

Downloads a file with a text progress bar showing download progress.

Usage

download_with_progress(url, destfile, verbose = TRUE)

Arguments

Value

Logical indicating success

Estimate lag-1 temporal autocorrelation

Description

Estimate the lag-1 autocorrelation rho_t of a leaf-wax record's residuals after a flat-mean detrend, ordering by age. This is the quantity that enters the within-record detection threshold from manuscript Section 4.5.3 (⁠Var(X1 - X2) = 2 sigma^2 (1 - rho_t)⁠).

Usage

estimate_temporal_autocorrelation(
  d2h_wax,
  age,
  method = c("ar1", "lomb_scargle")
)

Arguments

Details

Two methods are supported:

• "ar1" (default): Pearson correlation of resid[-n] with resid[-1] after age-ordering. For irregularly sampled series this is an approximation; see "lomb_scargle" for an alternative.

• "lomb_scargle": not yet implemented. Returns an error pointing the user at "ar1" until the spectral implementation lands. The plan is to estimate rho_t from the dominant timescale of a Lomb-Scargle periodogram on the irregularly sampled series.

Value

Numeric scalar in ⁠[-1, 1]⁠, or NA_real_ when the residuals are constant (e.g., n < 3 finite samples).

Examples

set.seed(1)
n <- 200
rho <- 0.7
e  <- numeric(n); e[1] <- rnorm(1, 0, 5)
for (k in 2:n) e[k] <- rho * e[k-1] + rnorm(1, 0, 5 * sqrt(1 - rho^2))
ag <- seq(0, 10000, length.out = n)
estimate_temporal_autocorrelation(-150 + e, ag)

Example leaf wax hydrogen isotope data

Description

A 10-row data frame of synthetic leaf wax hydrogen isotope measurements bundled for demonstration and testing.

Usage

example_data

Format

A data frame with 10 rows and 10 variables:

site_id

Character, site identifier

longitude

Numeric, longitude in decimal degrees (-180 to 180)

latitude

Numeric, latitude in decimal degrees (-90 to 90)

elevation

Numeric, elevation in meters above sea level

d2h_wax

Numeric, leaf wax hydrogen isotope value in per mil VSMOW

d2h_wax_sd

Numeric, analytical uncertainty in per mil

c4_fraction

Numeric, C4 vegetation fraction (0-1)

pft_tree

Numeric, tree plant functional type fraction (0-1)

pft_shrub

Numeric, shrub plant functional type fraction (0-1)

pft_grass

Numeric, grass plant functional type fraction (0-1)

Source

Synthetic values designed to span the calibration range.

Examples

data(example_data)
head(example_data)

Generate Fibonacci sphere points

Description

Generates evenly distributed points on a sphere using the Fibonacci spiral method.

Usage

generate_fibonacci_sphere(n_points = N_SPATIAL_KNOTS)

Arguments

Value

Matrix with columns "lon" and "lat" in degrees

Generate human-readable model description

Description

Generate human-readable model description

Usage

generate_model_description(model_name, capabilities)

Arguments

Value

Character description

Get all model metadata

Description

Returns a list of all 14 models with their descriptions and properties. The models are based on the spatial leafwax hierarchical Bayesian models.

Usage

get_all_model_metadata()

Value

Named list of available models with metadata

Get leafwax data cache directory

Description

Returns the path to the local cache directory for leafwax model data. Uses rappdirs for platform-specific paths or a user-specified directory.

Usage

get_cache_dir(create = TRUE)

Arguments

Value

Character string with the cache directory path

Examples

local({
  old <- options(leafwax.cache_dir = file.path(tempdir(), "leafwax_cache"))
  on.exit(options(old))

  cache_dir <- get_cache_dir(create = FALSE)
  dir.exists(cache_dir)
})

Get cache files for a model

Description

Internal helper that returns the cached file paths for a model. The v0.2 download layout ships a single posterior per model at ⁠posteriors/<model>_posterior.rds⁠, so the returned vector has at most one element. The data_type argument is accepted for API compatibility but does not affect the result.

Usage

get_cache_files(model_name, data_type, cache_dir)

Arguments

Value

Character vector of cached file paths that exist on disk.

Get cache size information

Description

Reports the disk space used by cached model data.

Usage

get_cache_info(by_model = FALSE, by_type = FALSE)

Arguments

Value

Data frame with cache size information

Examples

local({
  old <- options(leafwax.cache_dir = file.path(tempdir(), "leafwax_cache"))
  on.exit(options(old))

  # Get total cache size
  cache_info <- get_cache_info()

  # Get size by model and type
  cache_info <- get_cache_info(by_model = TRUE, by_type = TRUE)
})

Get data manifest

Description

Loads or downloads the data manifest with file checksums. Returns NULL (with a warning()) when the manifest is unreachable and there is no cached copy on disk; callers must treat that as "checksum verification skipped" rather than "no checksums found".

Usage

get_data_manifest()

Value

Parsed manifest list, or NULL if no manifest is available locally and the download failed.

Get path to data file

Description

Returns the full path to a data file based on the data source.

Usage

get_data_path(filename, data_source = "auto")

Arguments

Value

Character string with the file path

Get data download URLs

Description

Constructs download URLs for model data from GitHub releases.

Usage

get_data_url(model_name, version = "latest", data_type = c("posteriors"))

Arguments

Value

List of download URLs and filenames

Examples

# Get URLs for latest version
urls <- get_data_url("baseline_sp", "latest")

# Get URLs for specific version
urls <- get_data_url("baseline_sp", "v1.0.1")

Get model info

Description

Get model info

Usage

get_model_info(model_name)

Arguments

Value

List with model metadata

Get model parameters

Description

Returns which parameters each model expects and provides metadata about model capabilities for validation.

Usage

get_model_parameters(model_name)

Arguments

Value

List with model parameters and capabilities

Get URL configuration

Description

Loads the URL configuration from the package data.

Usage

get_url_config()

Value

List with URL configuration

Ensemble predictions across multiple models

Description

Runs invert_d2H against each model in models and combines the per-draw reconstructions per site, preserving the per-site dimension. Useful for downstream uncertainty estimates that should span structural model uncertainty rather than condition on one calibration variant.

Usage

invert_d2H_ensemble(
  ...,
  models = c("full_sp", "full_interact_sp", "elevation_c4_interact_sp"),
  ensemble_method = c("equal", "all")
)

Arguments

Value

If ensemble_method = "equal", a list with: posterior_draws (an ⁠n_draws x n_sites⁠ matrix of pooled per-site, per-draw reconstructions), ensemble_summary (a data frame with one row per site: mean, median, sd, ci_90_lower/ci_90_upper, ci_95_lower/ci_95_upper, n_models_used), model_results (the per-model output from invert_d2H() with return_full = TRUE), models_used (the models actually pooled), and ensemble_method. If ensemble_method = "all", only model_results and ensemble_method are returned.

Invert leaf wax d2H to precipitation d2H

Description

Uses Bayesian posterior draws to invert leaf wax hydrogen isotope values to precipitation isotope values, accounting for all fitted model components including vegetation effects and spatial correlations where applicable.

Usage

invert_d2h(
  d2h_wax,
  d2h_wax_err = NULL,
  longitude,
  latitude,
  elevation = NULL,
  c4_percent = NULL,
  pft_tree = NULL,
  pft_shrub = NULL,
  pft_grass = NULL,
  model_name = "baseline",
  n_draws = NULL,
  return_full = FALSE,
  credible_level = 0.9,
  verbose = TRUE,
  record_id = NULL,
  slope = NULL
)

invert_d2H(
  d2H_wax,
  d2H_wax_sd = NULL,
  longitude,
  latitude,
  elevation = NULL,
  elevation_sd = 100,
  c4_fraction = NULL,
  c4_fraction_sd = 10,
  pft_tree = NULL,
  pft_shrub = NULL,
  pft_grass = NULL,
  model_name = "baseline",
  n_posterior_draws = NULL,
  return_full = FALSE,
  credible_level = 0.9,
  verbose = TRUE,
  record_id = NULL,
  slope = NULL
)

Arguments

Value

If return_full is FALSE, a data frame with columns:

The interval is the posterior predictive specified in manuscript supplement Section S4.1, Eq. 7: the wax-error draw combines analytical uncertainty with the model's posterior residual SD sigma. For within-record change detection, the spatial GP intercept's contribution cancels in any contrast computed from the returned posterior_draws (manuscript Section 4.5.3); the same sigma applies in both regimes.

If return_full is TRUE, a list with:

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  # Simple inversion with base model
  results <- invert_d2h(
    d2h_wax = c(-150, -140, -130),
    d2h_wax_err = c(3, 3, 3),
    longitude = c(-120, -110, -100),
    latitude = c(40, 35, 30),
    elevation = c(1000, 1500, 500),
    model_name = "baseline",
    verbose = FALSE
  )

  # Inversion with spatial model
  results <- invert_d2h(
    d2h_wax = c(-150, -140, -130),
    d2h_wax_err = c(3, 3, 3),
    longitude = c(-120, -110, -100),
    latitude = c(40, 35, 30),
    elevation = c(1000, 1500, 500),
    model_name = "baseline_sp",
    return_full = TRUE,
    verbose = FALSE
  )
})

Get leafwax configuration

Description

Returns current configuration options for the leafwax package.

Usage

leafwax_config(option = NULL)

Arguments

Value

List of options or single option value

Examples

# Get all configuration options
cfg <- leafwax_config()

# Get specific option
auto_download <- leafwax_config("auto_download")

Set leafwax configuration

Description

Sets configuration options for the leafwax package.

Usage

leafwax_set_config(..., persist = TRUE)

Arguments

Value

Invisible NULL

Examples

# Enable auto-download
old <- options()
leafwax_set_config(auto_download = TRUE, persist = FALSE)
options(old)

# Set multiple options
old <- options()
leafwax_set_config(
  auto_download = TRUE,
  cache_dir = "~/my_leafwax_cache",
  verbose = FALSE,
  persist = FALSE
)
options(old)

List available models in cache

Description

Lists all models that have been downloaded to the local cache.

Usage

list_cached_models(data_type = NULL, verbose = TRUE)

Arguments

Value

Character vector of available model names

Examples

local({
  old <- options(leafwax.cache_dir = file.path(tempdir(), "leafwax_cache"))
  on.exit(options(old))

  # List all cached models
  models <- list_cached_models(verbose = FALSE)

  # List models with full data
  models_full <- list_cached_models(data_type = "full", verbose = FALSE)
})

List model names

Description

Returns the names of every model the package can resolve. Prefers the heavy posteriors directory when present (development install), otherwise falls back to the lightweight posteriors directory that ships with every install. The user cache is intentionally not enumerated here so the answer is stable regardless of what has been downloaded.

Usage

list_model_names()

Value

Character vector of model names. Empty if neither directory contains posterior files.

List available models with details

Description

Returns information about the v10 models available in the leafwax package, including which covariates each model uses.

Usage

list_models(check_data = TRUE, verbose = TRUE)

Arguments

Value

Data frame with model information

Examples

# List all models
models <- list_models(verbose = FALSE)
models_head <- head(models)

Load posterior draws for a model

Description

Loads posterior draws for one of the 14 leafwax v10 models. The function searches three tiers in order:

Usage

load_posteriors(model_name, n_draws = NULL, verbose = TRUE)

Arguments

Details

1. Heavy posteriors at ⁠inst/extdata/posteriors/⁠ (1000 draws, development install only; excluded from the CRAN tarball).

2. Cache populated by download_model_data() under get_cache_dir().

3. Preview posteriors at ⁠inst/extdata/posteriors_light/⁠. These are a 100-draw stratified subsample shipped with every install so examples and tests run offline. They are intended as a fixture for code-path verification, not for inference: tail probabilities and 95% credible intervals are noisy at this sample size. The package issues a warning whenever the preview tier is in use; downstream functions (invert_d2H(), assess_claim(), detect_change()) repeat the warning so it is visible at the call that actually matters.

For inference, run download_model_data() once to populate the cache and then call load_posteriors() again – the cache tier wins over the preview tier and no further downloads are needed.

Value

A leafwax_posterior object: a list with draws, metadata (including metadata$tier, one of "heavy", "cache", "light"), optional spatial, and accessor closures.

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  # Load a model (preview tier on a fresh install)
  model <- load_posteriors("baseline", verbose = FALSE)

  # Spatial model with limited draws
  model_fast <- load_posteriors("baseline_sp", n_draws = 50, verbose = FALSE)
})

Local effective slope at a paleo-reconstruction site

Description

Returns a per-draw vector of the local \beta_{\delta^2 H_p} calibration slope at a single site, combining the global precipitation-isotope slope with the spatial slope GP prediction at that site. The returned vector is the raw posterior at the site; every draw the calibration produced is preserved without modification.

Usage

local_effective_slope(
  longitude,
  latitude,
  model_name,
  override = NULL,
  n_draws = NULL,
  verbose = FALSE
)

Arguments

Details

Two modes:

• Default: returns the model's per-draw slope at the site.

• Override (single value or per-draw vector) replaces the model slope with a defended local value (e.g., from independent evidence about source-water seasonality, leaf-water enrichment, or vegetation).

Pass the returned vector to invert_d2H(..., slope = ...) to propagate it through the inversion.

Mechanistic reference values (e.g. the simple two-pool stationarity bound alpha = 1 + epsilon_app/1000 ~ 0.88 under ⁠epsilon_app ~= -120 permil⁠; Sessions 2005) are documented for interpretation but are never applied to the returned draws. The frequency of draws above any chosen reference is computable directly from the returned vector (mean(slope > 0.88)) and carries scientific information about how often the calibration implicates non-stationarity at the site.

Value

Numeric vector of length n_draws, the per-draw effective slope at the site (after override, if any).

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  # St. Louis with the baseline_sp model
  s <- local_effective_slope(
    longitude = -90, latitude = 38,
    model_name = "baseline_sp",
    n_draws = 200
  )
  slope_summary <- summary(s)

  # How often does the calibration imply a slope above the simple-model
  # stationarity bound at this site?
  fraction_above_bound <- mean(s > 0.88)

  # Override with a defended local slope
  s_fixed <- local_effective_slope(
    longitude = -90, latitude = 38,
    model_name = "baseline_sp",
    override = 0.55
  )

  # Pass through to the inversion. The slope vector and the
  # inversion's posterior must use the same n_draws: pair
  # local_effective_slope(..., n_draws = N) with
  # invert_d2H(..., n_posterior_draws = N, slope = s), or pass a
  # single point estimate (e.g., median(s)).
  result <- invert_d2H(d2H_wax = -180, d2H_wax_sd = 3,
                       longitude = -90, latitude = 38,
                       model_name = "baseline_sp",
                       n_posterior_draws = 200,
                       slope = s,
                       verbose = FALSE)
})

Predict precipitation d2H from leaf wax d2H

Description

Main user-facing function for inverting leaf wax hydrogen isotopes to precipitation isotopes. Automatically selects appropriate model based on available data and returns results in a tidy format.

Usage

predict_d2h_precip(
  data = NULL,
  d2h_wax = NULL,
  longitude = NULL,
  latitude = NULL,
  d2h_wax_err = NULL,
  elevation = NULL,
  c4_fraction = NULL,
  pft_tree = NULL,
  pft_shrub = NULL,
  pft_grass = NULL,
  model = "auto",
  n_draws = NULL,
  credible_level = 0.9,
  return_draws = FALSE,
  progress = TRUE,
  verbose = TRUE
)

Arguments

Value

A data frame with predictions (or list if return_draws = TRUE):

d2h_precip_mean

Mean predicted precipitation d2H

d2h_precip_median

Median predicted precipitation d2H

d2h_precip_sd

Standard deviation of the posterior predictive interval

d2h_precip_lower

Lower bound of the credible interval

d2h_precip_upper

Upper bound of the credible interval

prediction_interval_width

Width of the credible interval

model_used

Name of model used for prediction

The interval is the posterior predictive specified in manuscript supplement Section S4.1, Eq. 7 (analytical uncertainty plus the model's posterior residual SD).

Examples

local({
  old <- options(leafwax.suppress_preview_warning = TRUE)
  on.exit(options(old))

  # Using data frame input
  data(example_data)
  results <- predict_d2h_precip(example_data, verbose = FALSE)

  # Using individual vectors
  results <- predict_d2h_precip(
    d2h_wax = c(-150, -140, -130),
    longitude = c(-120, -110, -100),
    latitude = c(40, 35, 30),
    elevation = c(1000, 1500, 500),
    verbose = FALSE
  )

  # Specify model explicitly
  results <- predict_d2h_precip(
    example_data,
    model = "baseline_env_sp",
    verbose = FALSE
  )

  # Get full posterior draws
  results <- predict_d2h_precip(
    example_data,
    return_draws = TRUE,
    verbose = FALSE
  )
})

Predict an mPP Gaussian-process random effect at a new location

Description

Single-GP version. Used internally by predict_spatial_dual_gp() for each of the two (intercept, slope) fields. Matches the Matern 3/2 kernel and standardized-coordinate convention from the v10 Stan model.

Usage

predict_one_gp_mpp(
  coords_new,
  knot_coords,
  z_knots,
  sigma_draws,
  ls_km_draws,
  scaling,
  jitter = 1e-04
)

Arguments

Value

matrix(n_draws, n_obs) of predicted GP values at the new sites.

Predict both spatial intercept and spatial slope at new locations

Description

v10 carries two independent GPs. Both share knot coordinates and a single length scale parameter (ls_intercept_km == ls_slope_km in v10's posterior, two names for the same draw), but have distinct sigma_intercept_spatial and sigma_slope_spatial, and distinct ⁠z_intercept_spatial[*]⁠ and ⁠z_slope_spatial[*]⁠ knot effects.

Usage

predict_spatial_dual_gp(coords_new, knot_coords, draws, scaling)

Arguments

Value

list with two matrices, each n_draws x n_obs: intercept (additive contribution to beta_0 in standardized d2H_wax space) and slope (additive contribution to the local \beta_{\delta^2 H_p} slope).

Print method for leafwax_posterior

Description

Print method for leafwax_posterior

Usage

## S3 method for class 'leafwax_posterior'
print(x, ...)

Arguments

Value

The input leafwax_posterior object x, invisibly. Called for its side effect of printing a one-line model summary to the console.

Process chunks in parallel

Description

Process chunks in parallel

Usage

process_parallel(data, chunks, model, n_cores, progress, ...)

Arguments

Value

List of results for each chunk

Process chunks sequentially with progress bar

Description

Process chunks sequentially with progress bar

Usage

process_sequential(data, chunks, model, progress, ...)

Arguments

Value

List of results for each chunk

Select best model based on available data

Description

Automatically selects the most appropriate v10 model name from the 14 shipped variants given which covariates the user has available. Spatial-aware models are preferred when prefer_spatial = TRUE.

Usage

select_best_model_from_flags(
  has_elevation = FALSE,
  has_c4 = FALSE,
  has_pft = FALSE,
  prefer_spatial = TRUE,
  verbose = FALSE
)

Arguments

Value

Character string with selected v10 model name

Validate input data for inversion

Description

Checks that input data meets requirements for the specified model and returns cleaned, validated data.

Usage

validate_inputs(
  d2h_wax,
  longitude,
  latitude,
  d2h_wax_err = NULL,
  elevation = NULL,
  c4_fraction = NULL,
  pft_tree = NULL,
  pft_shrub = NULL,
  pft_grass = NULL,
  model_name = "baseline"
)

Arguments

Value

List of validated inputs

Validate inputs for a specific model

Description

Checks that all required predictors are provided and warns about unused ones.

Usage

validate_model_inputs(
  model_name,
  d2h_wax,
  longitude,
  latitude,
  elevation = NULL,
  c4_percent = NULL,
  pft_tree = NULL,
  pft_shrub = NULL,
  pft_grass = NULL,
  verbose = TRUE
)

Arguments

Value

List with validation results