API Reference

Core Infrastructure

DigiMuh: a toolkit for ingesting, storing, and querying dairy-cow environmental and physiological sensor data.

The package consolidates heterogeneous CSV exports from on-farm monitoring systems into a single normalised SQLite database.

Central constants for all DigiMuh analyses and visualisations.

digimuh.constants.THI_REFERENCE = 68.8

THI mild-stress threshold from Hoffmann et al. (2020).

digimuh.constants.DELTA_STABLE = 1.0

Max |Δ breakpoint| classified as ‘stable’.

digimuh.constants.DELTA_MODERATE = 3.0

Max |Δ breakpoint| classified as ‘decreased’ / ‘increased’.

digimuh.constants.MIN_COHORT_SIZE = 10

Minimum animals for longitudinal Sankey plots.

digimuh.constants.THI_RANGE = (45, 80)

Default x_range for THI → body temp fits.

digimuh.constants.TEMP_RANGE = (5, 35)

Default x_range for barn temp → body temp fits.

digimuh.constants.MIN_READINGS = 50

Minimum readings per animal-year to attempt a fit.

digimuh.constants.GRID_STEPS = 200

Number of breakpoint candidates in grid search.

Hierarchical configuration for DigiMuh analysis pipelines.

Priority (highest wins):

1. CLI arguments (always override everything)

2. .env in the project directory (quick per-project overrides)

3. ~/.config/digimuh/config.yaml (machine-specific, never in repo)

4. Built-in defaults

Usage in an entry point:

from digimuh.config import load_config

cfg = load_config()
# cfg.database   → Path to cow.db
# cfg.output     → Path to results directory
# cfg.tierauswahl → Path to Tierauswahl.xlsx
# cfg.n_jobs     → Number of parallel workers

CLI arguments still work and override everything:

digimuh-extract --db /other/cow.db --out /tmp/test

Setup a new machine:

digimuh-config

This creates ~/.config/digimuh/config.yaml interactively.

class digimuh.config.DigiMuhConfig(database=None, output=<factory>, tierauswahl=None, n_jobs=20, smaxtec_drink_correction=False, _sources=<factory>)

Resolved configuration for a DigiMuh pipeline run.

Parameters:

• database (Path | None)

• output (Path)

• tierauswahl (Path | None)

• n_jobs (int)

• smaxtec_drink_correction (bool)

• _sources (dict)

database: Path | None = None

output: Path

tierauswahl: Path | None = None

n_jobs: int = 20

smaxtec_drink_correction: bool = False

digimuh.config.load_config(cli_args=None, project_root=None)

Load configuration with full priority chain.

Parameters:

• cli_args (Namespace | None) – Parsed argparse namespace (from entry point).

• project_root (Path | None) – Project root for .env lookup. Defaults to cwd.

Returns:

Resolved DigiMuhConfig.

Return type:

DigiMuhConfig

digimuh.config.print_config(cfg)

Log the resolved configuration with sources.

Parameters:

cfg (DigiMuhConfig)

Return type:

None

digimuh.config.setup_interactive()

Interactive setup for a new machine. Creates config.yaml.

Return type:

None

digimuh.config.main()

Entry point for digimuh-config.

Return type:

None

Pretty-printed console output for DigiMuh analysis pipelines.

Uses rich for coloured panels, tables, progress bars, and tree views. Falls back to plain logging if rich is not installed.

Usage:

from digimuh.console import console, section, result_table, progress

console.print("[bold blue]Starting analysis …[/]")
with progress("Fitting animals") as pb:
    task = pb.add_task("Broken-stick", total=220)
    for animal in animals:
        fit(animal)
        pb.advance(task)
section("Results", "Breakpoint analysis complete")
result_table("Convergence", headers, rows)

digimuh.console.setup_logging(level=20)

Configure logging with rich handler if available.

Parameters:

level (int)

Return type:

None

digimuh.console.section(title, subtitle='')

Print a section header.

Parameters:

• title (str)

• subtitle (str)

Return type:

None

digimuh.console.reset_steps()

Reset the step counter (for testing).

Return type:

None

digimuh.console.result_table(title, headers, rows, highlight_col=None)

Print a formatted results table.

Parameters:

• title (str)

• headers (list[str])

• rows (list[list[Any]])

• highlight_col (int | None)

Return type:

None

digimuh.console.stars_styled(stars)

Return rich-styled significance stars.

Parameters:

stars (str)

Return type:

str

digimuh.console.kv(key, value, indent=2)

Print a key-value pair.

Parameters:

• key (str)

• value (Any)

• indent (int)

Return type:

None

digimuh.console.kv_pair(key, val1, val2, sep=' / ')

Print a key with two values (e.g. converged / total).

Parameters:

• key (str)

• val1 (Any)

• val2 (Any)

• sep (str)

Return type:

None

digimuh.console.progress(description='Processing')

Context manager for a rich progress bar.

Parameters:

description (str)

digimuh.console.banner(title, version='0.1.0')

Print the DigiMuh startup banner.

Parameters:

• title (str)

• version (str)

Return type:

None

digimuh.console.done(message='All analyses complete.')

Print completion message.

Parameters:

message (str)

Return type:

None

Central routing of result artefacts into subject-specific subfolders.

The six subjects correspond to the §3.x sections in docs/analysis_00_methods.md:

01_extract — raw per-animal / per-cow-day CSVs pulled from the DB 02_breakpoints — broken-stick fits, diagnostics, Spearman 03_temporal — circadian, crossings, cross-correlation, ETA 04_production — TNF × yield, Wood residuals, milk-yield class 05_composition — MLP thin-milk + dilution partition 06_longitudinal — year-on-year stability, Sankey, raincloud

Call sites use resolve_output(data_dir, filename) on both write and read paths. When a file’s subject cannot be determined the helper falls back to the flat root so unrecognised outputs still land somewhere predictable.

digimuh.paths.resolve_output(data_dir, filename, *, create=True)

Return the canonical path for an output file inside data_dir.

When the filename’s subject is in the routing table, the file lands in (and is read from) the corresponding subject subfolder. When it isn’t, the file stays at the flat root. The target directory is created on demand unless create=False.

Parameters:

• data_dir (Path) – The analysis output directory (e.g. results/broken_stick).

• filename (str) – A basename (no path components) like "rumen_barn.csv" or "diagnostic_thi_top10.svg".

• create (bool) – When True (default), mkdir(parents=True, exist_ok=True) is called on the resolved folder before it is returned, so callers never have to.

Returns:

data_dir / SUBDIRS[subject] / filename when a subject is known, otherwise data_dir / filename.

Return type:

Path

digimuh.paths.resolve_input(data_dir, filename)

Return the path to read filename from.

Preference order:

1. Subject subfolder (if known from the routing table).

2. Flat root (legacy layout, for back-compat with older runs).

Does not create directories. Returns whichever path exists; if neither exists, returns the subject-subfolder path so callers get a consistent “canonical” error via FileNotFoundError.

Parameters:

• data_dir (Path)

• filename (str)

Return type:

Path

Shared database connection, view initialisation, and plotting defaults used across all DigiMuh analysis modules.

Usage:

from digimuh.analysis_utils import connect_db, setup_plotting
con = connect_db(Path("cow.db"))
setup_plotting()

digimuh.analysis_utils.connect_db(db_path, create_views=True)

Open the cow database and optionally create analysis views.

Parameters:

• db_path (Path) – Path to the SQLite database file.

• create_views (bool) – If True, execute create_views.sql to ensure all analysis views exist.

Returns:

An open sqlite3.Connection with Row factory enabled.

Return type:

Connection

digimuh.analysis_utils.query_df(con, sql, params=())

Execute a query and return a pandas DataFrame.

Parameters:

• con (sqlite3.Connection) – Active database connection.

• sql (str) – SQL query string.

• params (tuple) – Bind parameters for the query.

Returns:

A pandas.DataFrame with the query results.

Return type:

pd.DataFrame

digimuh.analysis_utils.setup_plotting()

Configure matplotlib for publication-quality figures.

Sets SVG text as editable (svg.fonttype = "none"), uses a clean style, and configures reasonable defaults for axis labels and tick sizes.

Return type:

None

digimuh.analysis_utils.save_fig(fig, name, out_dir)

Save a figure as SVG, PNG, and close it.

Parameters:

• fig (Figure) – The matplotlib figure to save.

• name (str) – Base filename (without extension).

• out_dir (Path) – Output directory.

Return type:

None

Data Ingestion & Extraction

DigiMuh ingestion engine — CSV to SQLite star-schema builder.

Reads ~8.9 GB of heterogeneous dairy-cow CSV data exported from multiple on-farm sensor systems (smaXtec, HerdePlus, Gouna, HOBO, LoRaWAN, DWD weather) and consolidates them into a single normalised SQLite database with dimension tables for animals, sensors, barns, and source-file provenance.

The resulting database uses a star schema: small dimension tables hold entity metadata; large fact tables hold time-series measurements with foreign keys pointing back to the dimensions. Timestamps are stored as ISO-8601 TEXT (SQLite’s date/time functions work natively on this representation). Composite indexes on (entity_id, timestamp) are created for every fact table to accelerate the most common query pattern: “give me all X for animal Y between dates A and B”.

Usage:

# Smoke test — first 5 files per folder
python -m digimuh.ingest_cow_db /path/to/csv/root --db cow.db --test-n 5

# Full ingestion (~2–3 h depending on disk I/O)
python -m digimuh.ingest_cow_db /path/to/csv/root --db cow.db

# Verbose mode — prints CREATE TABLE SQL for each table
python -m digimuh.ingest_cow_db /path/to/csv/root --db cow.db -v

Dependencies:

• Python ≥ 3.10 (for X | Y union type hints)

• tqdm (optional, for progress bars)

See also

• docs/database_structure.md — full schema documentation

• docs/column_dictionary.md — column-level data dictionary

Authors:

Bart R. H. Geurten, Claude (Anthropic)

digimuh.ingest.FOLDER_CONFIG: dict[str, dict] = {'output_allocations': {'entity_col': 'animal_id', 'entity_source': 'column_animal_id', 'glob': '*.csv', 'table': 'allocations'}, 'outputs_bcs': {'entity_source': 'filename_animal', 'glob': '*.csv', 'table': 'bcs'}, 'outputs_gouna': {'entity_source': 'filename_animal', 'glob': '*.csv', 'table': 'gouna'}, 'outputs_herdeplus_mlp_gemelk_kalbung': {'entity_source': 'filename_animal', 'glob': '*.csv', 'table': 'herdeplus'}, 'outputs_hobo': {'entity_source': None, 'glob': '*.csv', 'table': 'hobo_weather'}, 'outputs_lorawan': {'entity_source': 'filename_sensor', 'glob': '*.csv', 'table': 'lorawan'}, 'outputs_smaxtec_barns': {'entity_source': 'filename_barn', 'glob': '*.csv', 'table': 'smaxtec_barns'}, 'outputs_smaxtec_derived': {'entity_source': 'filename_animal', 'glob': '*.csv', 'table': 'smaxtec_derived'}, 'outputs_smaxtec_events': {'entity_source': 'filename_animal', 'glob': '*.csv', 'table': 'smaxtec_events'}, 'outputs_smaxtec_water_intake': {'entity_source': 'filename_animal', 'glob': '*.csv', 'table': 'smaxtec_water_intake'}}

Mapping of folder names to ingestion configuration.

Each value dict contains:

• table — target SQLite table name

• entity_source — how the entity ID is derived

• entity_col — (optional) CSV column holding the entity ID

• glob — file-glob pattern for discovering CSVs

Type:

dict

digimuh.ingest.STANDALONE_FILES: dict[str, dict] = {'herdeplus_diseases.csv': {'entity_col': 'cow', 'entity_source': 'column_cow', 'table': 'diseases'}, 'outputs_dwd.csv': {'entity_source': None, 'table': 'dwd_weather'}}

Single-file sources that live in the root directory.

Type:

dict

digimuh.ingest.DIMENSION_TABLES_SQL: str = '\n-- ┌──────────────────────────────────────────────────────────────┐\n-- │  Dimension tables                                           │\n-- └──────────────────────────────────────────────────────────────┘\n\nCREATE TABLE IF NOT EXISTS animals (\n    animal_id  INTEGER PRIMARY KEY   -- EU ear tag, IS the rowid\n);\n\nCREATE TABLE IF NOT EXISTS sensors (\n    sensor_id    INTEGER PRIMARY KEY AUTOINCREMENT,\n    sensor_name  TEXT    NOT NULL UNIQUE\n);\n\nCREATE TABLE IF NOT EXISTS barns (\n    barn_id    INTEGER PRIMARY KEY AUTOINCREMENT,\n    barn_name  TEXT    NOT NULL UNIQUE\n);\n\nCREATE TABLE IF NOT EXISTS source_files (\n    file_id   INTEGER PRIMARY KEY AUTOINCREMENT,\n    filename  TEXT    NOT NULL,\n    folder    TEXT    NOT NULL,\n    UNIQUE(filename, folder)\n);\n'

SQL script to create the four dimension tables.

animals uses the 15-digit EU ear tag as its INTEGER PRIMARY KEY which makes it an alias for SQLite’s internal rowid — the fastest possible lookup. The other dimensions use AUTOINCREMENT surrogates because their natural keys are short strings.

Type:

str

digimuh.ingest.extract_entity_id_from_filename(filename)

Extract the first underscore-delimited segment from a filename.

The DigiMuh CSV export encodes the entity identifier (animal ear tag, sensor name, or barn name) as the first segment of the filename, separated from the rest by an underscore.

Parameters:

filename (str) – Basename of the CSV file, e.g. "276001260919234_smaxtec_derived_2021-04-01_2024-09-30.csv".

Returns:

The first segment, e.g. "276001260919234".

Return type:

str

Example

>>> extract_entity_id_from_filename("CU-1_LoRaWAN_raw_2021.csv")
'CU-1'

digimuh.ingest.guess_column_type(values)

Infer an SQLite column type from a sample of string values.

Tries to parse every non-empty value as int, then float. Falls back to TEXT if either parse fails. An all-empty sample also returns TEXT.

Parameters:

values (list[str]) – Sample of raw string values from a CSV column.

Returns:

One of "INTEGER", "REAL", or "TEXT".

Return type:

str

Note

SQLite uses dynamic typing so the column affinity is advisory, but correct affinity improves storage efficiency and query planner decisions.

digimuh.ingest.read_csv_sample(filepath, n_rows=100)

Read headers and up to n_rows data rows from a CSV file.

Used during schema inference to peek at column headers and a small sample of data for type guessing.

Parameters:

• filepath (Path) – Path to the CSV file.

• n_rows (int) – Maximum number of data rows to read. The default of 100 is enough for most narrow tables; wide sparse tables like smaxtec_derived benefit from larger samples (see build_create_table_sql()).

Returns:

A (headers, rows) tuple where headers is a list of raw column-name strings and rows is a list of lists of string values.

Return type:

tuple[list[str], list[list[str]]]

digimuh.ingest.sanitise_column_name(name)

Clean a raw CSV header into a safe SQL column name.

Strips whitespace, lowercases, replaces any character that is not [a-z0-9_] with an underscore, collapses consecutive underscores, and strips leading/trailing underscores.

Parameters:

name (str) – Raw column header string from the CSV.

Returns:

A sanitised, SQL-safe column name.

Return type:

str

Example

>>> sanitise_column_name("BCS-Wert")
'bcs_wert'
>>> sanitise_column_name("21141733_1__Temperature")
'21141733_1_temperature'

digimuh.ingest.build_create_table_sql(table_name, csv_headers, csv_rows, entity_source, entity_col=None)

Build a CREATE TABLE statement by inspecting CSV structure.

Examines the raw headers and a sample of data rows to infer column types, then constructs the DDL statement with appropriate foreign-key references to dimension tables. Column names are sanitised via sanitise_column_name().

Parameters:

• table_name (str) – Target SQLite table name.

• csv_headers (list[str]) – Raw CSV column headers (unsanitised).

• csv_rows (list[list[str]]) – Sample data rows for type inference.

• entity_source (str | None) – How the entity ID is obtained — one of the values documented in FOLDER_CONFIG, or None.

• entity_col (str | None) – If the entity is identified by a column in the CSV itself ("column_animal_id" or "column_cow"), this is the name of that column.

Returns:

A (sql, col_names) tuple. sql is the full CREATE TABLE statement. col_names is an ordered list of the sanitised CSV column names (matching CSV column order), used to build the corresponding INSERT statement.

Return type:

tuple[str, list[str]]

class digimuh.ingest.CowDBIngester(root_dir, db_path, chunk_size=50000, verbose=False, test_n=None)

Orchestrates reading CSV files and inserting into SQLite.

Walks the DigiMuh export directory structure, creates dimension and fact tables on the fly by inspecting the first CSV in each folder, then streams rows in configurable batches using executemany for throughput.

The class caches all dimension-table lookups in memory (animals, sensors, barns, source files) so repeated INSERT-OR-IGNORE round trips are avoided for entities already seen.

Parameters:

• root_dir (Path) – Root directory containing all CSV folders and standalone files.

• db_path (Path) – Path where the SQLite database will be created.

• chunk_size (int) – Number of rows per INSERT batch. Larger values use more memory but reduce transaction overhead.

• verbose (bool) – If True, emit DEBUG-level log messages including full CREATE TABLE SQL.

• test_n (int | None) – If set, only ingest the first test_n files per folder (standalone files are always fully ingested). Useful for rapid schema and pipeline validation.

Example

>>> ingester = CowDBIngester(
...     root_dir=Path("/data/DigiMuh-Export"),
...     db_path=Path("cow.db"),
...     chunk_size=50_000,
...     test_n=5,
... )
>>> ingester.run()

run()

Execute the full ingestion pipeline.

This is the main entry point. It:

1. Creates dimension tables (animals, sensors, barns, source_files).

2. Iterates over FOLDER_CONFIG, creating each fact table from the first CSV’s schema and then streaming all CSV files in that folder.

3. Processes STANDALONE_FILES the same way.

4. Creates composite indexes on all fact tables.

5. Prints a summary of rows inserted, database size, and elapsed time.

Raises:

• sqlite3.Error – On unrecoverable database errors.

• Individual CSV failures are logged and skipped. –

Return type:

None

digimuh.ingest.main()

Parse command-line arguments and run the ingestion pipeline.

This function is the console-script entry point registered in pyproject.toml as digimuh-ingest.

Return type:

None

Quick validation of the ingested cow database.

Checks row counts, null rates, value ranges, referential integrity, and temporal coverage. Run immediately after ingestion to catch problems before launching analysis.

Usage:

python -m digimuh.validate_db --db cow.db

digimuh.validate_db.check_table_counts(con)

Verify all expected tables exist and report row counts.

Returns:

List of warning/error messages (empty = all good).

Parameters:

con (Connection)

Return type:

list[str]

digimuh.validate_db.check_null_rates(con)

Report null rates for key columns.

Returns:

List of warning messages for suspiciously high null rates.

Parameters:

con (Connection)

Return type:

list[str]

digimuh.validate_db.check_value_ranges(con)

Check that numeric values fall within plausible ranges.

Returns:

List of warnings for out-of-range values.

Parameters:

con (Connection)

Return type:

list[str]

digimuh.validate_db.check_temporal_coverage(con)

Report the date range of each timestamped table.

Returns:

List of warnings for unexpected gaps or ranges.

Parameters:

con (Connection)

Return type:

list[str]

digimuh.validate_db.check_referential_integrity(con)

Check foreign key relationships between fact and dimension tables.

Returns:

List of warnings for orphaned references.

Parameters:

con (Connection)

Return type:

list[str]

digimuh.validate_db.main()

Run all validation checks and print a summary.

Return type:

None

Extract rumen temperature, respiration, barn climate, and production data from the DigiMuh database into analysis-ready CSVs.

Run this once. Downstream scripts (stats_runner and viz_runner) read the CSVs and never touch the database.

Outputs (in --out directory):

tierauswahl.csv          Animal selection (cleaned from xlsx)
rumen_barn.csv           Rumen temp + barn THI/temp per 10-min tick
respiration_barn.csv     Respiration + barn THI/temp per reading
production.csv           Mean milk yield + lactation nr per animal
climate_daily.csv        Daily barn climate summary (Jun–Sep)

Usage:

digimuh-extract --db cow.db --tierauswahl Tierauswahl.xlsx \
    --out results/broken_stick

¹ Analysis led by Dr. med. vet. Gundula Hoffmann, ATB Potsdam.

digimuh.extract.load_tierauswahl(path)

Load and clean the collaborator-provided animal selection list.

Parameters:

path (Path)

Return type:

DataFrame

digimuh.extract.extract_rumen_barn(con, tierauswahl, exclude_drinking=True, climate_source='smaxtec')

Extract rumen temp + barn climate for all selected animals.

Parameters:

• con – Database connection.

• tierauswahl (DataFrame) – Animal selection DataFrame.

• exclude_drinking (bool) – If True, exclude drinking events + 15 min padding. If False, rely solely on smaXtec’s built-in temp_without_drink_cycles correction.

• climate_source (str) – "smaxtec" (barn sensors, default) or "hobo" (weather loggers, THI derived via NRC 1971).

Return type:

DataFrame

digimuh.extract.extract_respiration_barn(con, tierauswahl, climate_source='smaxtec')

Extract respiration + barn climate for all selected animals.

Parameters:

• tierauswahl (DataFrame)

• climate_source (str)

Return type:

DataFrame

digimuh.extract.extract_production(con, tierauswahl)

Extract milk yield and lactation number per animal.

Parameters:

tierauswahl (DataFrame)

Return type:

DataFrame

digimuh.extract.extract_daily_milk_yield(con, tierauswahl)

Extract daily milk yield (sum of milkings per day) per animal.

Each row is one animal-day with the total yield across all milkings that day. Typically 2-3 milkings per day.

Parameters:

• con – Database connection.

• tierauswahl (DataFrame) – Cleaned Tierauswahl DataFrame.

Returns:

DataFrame with animal_id, year, date, daily_yield_kg, n_milkings.

Return type:

DataFrame

digimuh.extract.extract_daily_milk_yield_full(con, tierauswahl)

Full-history daily yield per Tierauswahl animal.

Unlike extract_daily_milk_yield() this is not restricted to each animal’s datetime_enter → datetime_exit observation window — we pull every day on which the animal has a valid herdeplus_milked_mkg entry, across her entire HerdePlus history in the DB.

This is the data source for the full-history Wood (1967) fit used by digimuh.stats_lactation_curve.compute_wood_residuals(): a per-lactation curve converges far more often when it sees the full pre-peak / post-peak shape instead of a summer slice.

Parameters:

• con – DB connection.

• tierauswahl (DataFrame) – Cleaned Tierauswahl DataFrame (provides the animal_id set to query).

Returns:

DataFrame with animal_id, date (datetime64[ns]), daily_yield_kg, n_milkings.

Return type:

DataFrame

digimuh.extract.extract_mlp_test_days(con, tierauswahl)

Per-(animal, test-day) MLP composition records for Tierauswahl animals.

HerdePlus stores two interleaved channels in the same table:

• high-frequency per-milking events (herdeplus_milked_* populated, herdeplus_mlp_* null)

• monthly MLP (Milchleistungsprüfung) test-day analytics (herdeplus_mlp_* populated, herdeplus_milked_* typically null)

This helper extracts the MLP channel: fat %, protein %, fat-kg, lactose, somatic cell count, urea, fat-to-protein ratio, lactose-kg, and energy-corrected milk (ECM) at one-month resolution. These are the standard dairy health / composition analytics that would otherwise be missing from the analysis pipeline.

Parameters:

• con – DB connection.

• tierauswahl (DataFrame) – Cleaned Tierauswahl DataFrame.

Returns:

DataFrame with one row per (animal_id, test-day timestamp) and the eleven MLP columns plus herdeplus_calving_lactation.

Return type:

DataFrame

digimuh.extract.extract_calvings(con)

Extract all calving-confirmation events herd-wide.

Pulls one row per calving from smaxtec_events across all animals — not just the Tierauswahl — so that downstream Wood (1967) lactation-curve fits can determine DIM for any cow-day including lactations that started before the observation window.

Returns:

DataFrame with animal_id and calving_date (pandas datetime64[ns], day resolution).

Return type:

DataFrame

digimuh.extract.extract_climate(con, tierauswahl, climate_source='smaxtec')

Extract daily barn climate for each summer in the dataset.

Parameters:

• tierauswahl (DataFrame)

• climate_source (str)

Return type:

DataFrame

digimuh.extract.main()

Return type:

None

Statistics

Core statistical functions for the broken-stick analysis pipeline.

digimuh.stats_core.p_to_stars(p)

Convert p-value to significance stars.

Returns:

'***' if p < 0.001, '**' if p < 0.01, '*' if p < 0.05, 'n.s.' otherwise.

Parameters:

p (float)

Return type:

str

digimuh.stats_core.run_broken_stick_fits(rumen, resp, frontiers_only=False)

Fit broken-stick models per animal-year.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame.

• resp (DataFrame) – respiration_barn.csv DataFrame.

• frontiers_only (bool) – If True, skip Davies/pscore/Hill fits (Frontiers paper: broken-stick only). These methods are reserved for the COMPAG companion papers.

Returns:

One row per animal-year with breakpoint results.

Return type:

DataFrame

digimuh.stats_core.compute_spearman(rumen, resp)

Per-animal Spearman correlations.

Parameters:

• rumen (DataFrame)

• resp (DataFrame)

Return type:

DataFrame

digimuh.stats_core.compute_below_above(rumen, resp, bs_results)

Per-animal means below/above their individual THI breakpoint.

Parameters:

• rumen (DataFrame)

• resp (DataFrame)

• bs_results (DataFrame)

Return type:

DataFrame

digimuh.stats_core.run_statistical_tests(beh)

Run Fisher resampling tests within each year, BH-FDR corrected.

Uses reRandomStats FisherResamplingTest with medianDiff as the test statistic and 20,000 permutations.

Tests per year: - Body temp below vs above breakpoint

Parameters:

beh (DataFrame) – behavioural_response DataFrame.

Returns:

one row per test, with raw p, adjusted p, stars.

Return type:

DataFrame

Temporal analysis functions for the broken-stick pipeline.

digimuh.stats_temporal.compute_cross_correlation(rumen, bs_results, max_lag=24)

Cross-correlation and cross-covariance of climate vs rumen temp, computed separately below and above each animal’s breakpoint.

For each animal with a converged breakpoint (THI and barn temp), the rumen temperature time series is split at the breakpoint. The normalised cross-correlation function (CCF) and raw cross-covariance are computed for lags -max_lag to +max_lag (in units of the 10-min sampling interval, so max_lag=24 = 4 hours).

Positive lags mean the climate signal leads the rumen temperature response.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame (must have timestamp column).

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

• max_lag (int) – Maximum lag in samples (default 24 = 4 hours).

Returns:

animal_id, year, predictor (thi/temp), region (below/above), lag, xcorr, xcov.

Return type:

DataFrame with columns

digimuh.stats_temporal.compute_circadian_null_model(rumen, bs_results)

Rumen temperature circadian profile on non-stress days.

For each animal with a converged THI breakpoint, identifies days where barn THI stayed below the breakpoint for the entire day (no heat stress). Computes the mean rumen temperature at each clock hour across these cool days.

This gives the circadian null model: what rumen temperature looks like when the cow is entirely in the thermoneutral zone.

Also computes the profile on stress days (THI exceeded breakpoint at some point during the day) for comparison.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame.

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

Returns:

animal_id, year, hour, day_type (cool/stress), body_temp_mean, body_temp_std, n_readings.

Return type:

DataFrame with columns

digimuh.stats_temporal.compute_thi_daily_profile(rumen, bs_results)

Barn THI profile across 24h, by month.

For each month of the observation period, computes the mean barn THI at each clock hour. Overlaid with the herd median breakpoint, this shows when heat stress typically begins and ends each month.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame.

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

Returns:

year, month, month_label, hour, thi_mean, thi_std, thi_q25, thi_q75, n_readings. Also includes herd_median_bp.

Return type:

DataFrame with columns

digimuh.stats_temporal.compute_derivative_ccf(rumen, bs_results, max_lag=24)

Cross-correlation of rate-of-change signals: dTHI/dt vs dTbody/dt.

Instead of correlating the raw levels (contaminated by shared diurnal cycle), we correlate the temporal derivatives. This asks: “when the barn heats up, how long until the cow heats up?”

The derivative removes DC offsets and slow trends while preserving the temporal coupling of changes. The rumen’s thermal inertia (~100 L) means the derivative response is delayed by 30-90 min.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame.

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

• max_lag (int) – Maximum lag in 10-min samples (default 24 = 4 hours).

Returns:

animal_id, year, predictor, region, lag, lag_minutes, dxcorr (normalised CCF of derivatives).

Return type:

DataFrame with columns

digimuh.stats_temporal.compute_event_triggered_average(rumen, bs_results, window=36, min_gap=6, crossing_hour_range=None)

Peri-event average of rumen temperature around THI breakpoint crossings.

Finds moments when barn THI crosses the animal’s breakpoint upward (heat stress onset). Extracts a window of rumen temperature centred on each crossing event and averages across events.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame.

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

• window (int) – Half-window in 10-min samples (default 36 = 6 hours).

• min_gap (int) – Minimum samples between events to avoid overlap (default 6 = 1 hour).

• crossing_hour_range (tuple[int, int] | None) – If set, only include crossing events whose clock hour falls within [start, end). E.g. (8, 11) keeps crossings at 8:00-10:59. None = all hours.

Returns:

animal_id, year, predictor, event_id, relative_lag (samples from crossing), relative_minutes, body_temp, climate_val. Also a summary DataFrame.

Return type:

DataFrame with columns

digimuh.stats_temporal.compute_crossing_times(rumen, bs_results, min_gap=6)

Extract clock times of all breakpoint crossing events.

For each animal, finds moments when barn THI (or barn temp) crosses the animal’s individual breakpoint upward, and records the clock time. Used for the activation raster plot.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame.

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

• min_gap (int) – Minimum samples between events (default 6 = 1 hour).

Returns:

animal_id, year, predictor, breakpoint, crossing_timestamp, clock_hour, clock_minute, day_fraction.

Return type:

DataFrame with columns

digimuh.stats_temporal.compute_climate_eta(rumen, bs_results, window=36, min_gap=6, crossing_hour_range=(8, 11))

Climate signal around breakpoint crossings, normalised to the breakpoint.

For each animal and each crossing event (upward), extracts both barn THI and barn temperature in a ±window around the crossing. The trigger predictor is normalised by subtracting the animal’s breakpoint value, so y=0 corresponds to the threshold.

Two modes: - ‘thi’ crossings: THI normalised, barn temp as companion - ‘temp’ crossings: barn temp normalised, THI as companion

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame.

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

• window (int) – Half-window in 10-min samples (default 36 = 6 hours).

• min_gap (int) – Minimum samples between events (default 6 = 1 hour).

• crossing_hour_range (tuple[int, int] | None) – Restrict to crossings in this clock-hour range. Default (8, 11) = 8:00–10:59.

Returns:

animal_id, year, trigger (thi/temp), event_id, relative_minutes, thi_raw, temp_raw, thi_norm, temp_norm, breakpoint_thi, breakpoint_temp.

Return type:

DataFrame with columns

Production impact analysis for the broken-stick pipeline.

digimuh.stats_production.compute_thermoneutral_fraction(rumen, bs_results)

Compute daily thermoneutral fraction (TNF) per animal.

For each animal with a converged THI breakpoint, count the fraction of 10-min readings per calendar day where barn THI is at or below the individual breakpoint. This is the fraction of the day the animal spends within its thermoneutral zone.

Also computes analogous fraction for barn temperature breakpoints.

Parameters:

• rumen (DataFrame) – rumen_barn.csv DataFrame (needs timestamp, barn_thi, barn_temp).

• bs_results (DataFrame) – broken_stick_results.csv DataFrame.

Returns:

animal_id, year, date, thi_tnf, temp_tnf, n_readings, mean_thi, mean_body_temp.

Return type:

DataFrame with columns

digimuh.stats_production.compute_tnf_yield_analysis(tnf, daily_yield)

Merge daily TNF with daily milk yield for per-cow, per-day pairs.

For each cow, each day produces two floats: - thi_tnf: fraction of the day below the cow’s THI breakpoint (0-1) - relative_yield: that day’s milk yield / cow-specific P95 (0-1ish)

The P95 is the 95th percentile of each individual cow’s daily yields across her entire dataset (all years). This is a robust, cow-specific reference maximum that avoids outlier sensitivity of the absolute max.

Parameters:

• tnf (DataFrame) – Daily TNF DataFrame from compute_thermoneutral_fraction.

• daily_yield (DataFrame) – Daily milk yield DataFrame (animal_id, date, daily_yield_kg, year).

Returns:

animal_id, year, date, thi_tnf, temp_tnf, daily_yield_kg, yield_p95, relative_yield.

Return type:

DataFrame with one row per cow-day

digimuh.stats_production.classify_cow_years_by_wood_residual(wood_residuals, terciles=None, min_days=10)

Assign each (animal_id, year) a low / middle / high class.

The class is based on the cow-year’s mean Wood residual so that a single label follows each animal-lactation across all her cow-days. Tercile boundaries default to Q33/Q67 of the cow-year means themselves (not the cow-day level) so the three classes each hold roughly a third of cow-years.

Parameters:

• wood_residuals (DataFrame) – Output of compute_wood_residuals — one row per cow-day with animal_id, year, yield_residual.

• terciles (tuple[float, float] | None) – Optional (Q33, Q67) boundary override. If None the boundaries are computed from the cow-year means.

• min_days (int) – Require this many cow-days per (animal_id, year) for the mean to be trusted; cow-years below this are dropped.

Returns:

(class_table, (q33, q67)) — class_table has columns animal_id, year, mean_residual, n_days, yield_class (categorical, ordered low < middle < high).

Return type:

tuple[DataFrame, tuple[float, float]]

digimuh.stats_production.compute_tnf_yield_by_class(tnf_yield, class_table, wood=None)

Attach yield_class (and optional DIM-adjusted residual) to a TNF × yield table.

Parameters:

• tnf_yield (DataFrame) – tnf_yield.csv DataFrame (one row per cow-day with thi_tnf, temp_tnf, daily_yield_kg, relative_yield, year).

• class_table (DataFrame) – Output of classify_cow_years_by_wood_residual().

• wood (DataFrame | None) – Optional daily_milk_yield_wood.csv DataFrame. When given, yield_residual, yield_expected and dim are merged onto each cow-day on (animal_id, date) so downstream correlations can use the DIM-adjusted residual as the response and remove the within-cow-year lactation decline from the heat-stress signal.

Returns:

Filtered cow-day DataFrame restricted to cow-years that have a class; gains yield_class, mean_residual and (when wood is supplied) yield_residual, yield_expected, dim.

Return type:

DataFrame

digimuh.stats_production.tnf_yield_correlations_by_class(tnf_by_class, predictors=('thi_tnf', 'temp_tnf'), responses=('daily_yield_kg', 'relative_yield', 'yield_residual'))

Per-class Spearman correlations of TNF vs yield.

Parameters:

• tnf_by_class (DataFrame) – Output of compute_tnf_yield_by_class().

• predictors (tuple[str, ...]) – TNF columns to test.

• responses (tuple[str, ...]) – Yield columns to test.

Returns:

DataFrame of yield_class × predictor × response with n, n_animals, rs, p, slope, intercept, median_y. The class "pooled" is included as a reference row.

Return type:

DataFrame

digimuh.stats_production.compute_daily_crossing_flags(crossing_times)

Per (animal_id, year, date): did any THI / barn-temp crossing occur?

Consumes crossing_times.csv — one row per upward crossing of the cow’s individual breakpoint — and collapses to a cow-day table with two boolean flags so downstream plots can split cow-days into “days with a crossing” vs “days without”.

Parameters:

crossing_times (DataFrame) – crossing_times.csv DataFrame with animal_id, year, date, predictor ("thi" / "temp").

Returns:

DataFrame with animal_id, year, date, thi_crossed, temp_crossed (booleans), n_thi_crossings, n_temp_crossings (ints).

Return type:

DataFrame

digimuh.stats_production.attach_daily_climate_means(tnf_by_class, rumen)

Ensure cow-day rows carry mean_thi and mean_barn_temp.

mean_thi is usually already present from compute_thermoneutral_fraction; mean_barn_temp was added later so legacy CSVs may lack it. Missing columns are re-derived from rumen_barn.csv via a per-cow-day groupby.

Parameters:

• tnf_by_class (DataFrame)

• rumen (DataFrame)

Return type:

DataFrame

digimuh.stats_production.attach_crossing_flags(tnf_by_class, crossing_flags)

Merge the cow-day crossing flag table onto a TNF × class table.

Parameters:

• tnf_by_class (DataFrame)

• crossing_flags (DataFrame)

Return type:

DataFrame

digimuh.stats_production.crossing_day_comparison(df, response='yield_residual')

Per-predictor effect of “crossed that day” on a yield response.

For each climate predictor (THI, barn-temp) and for the pooled cow-days plus each yield_class (if present), reports:

• group sizes,

• medians,

• Mann-Whitney U statistic and p,

• median-difference effect size (crossed − not-crossed).

Parameters:

• df (DataFrame) – Cow-day DataFrame that has thi_crossed, temp_crossed, the response column, and optionally yield_class.

• response (str) – Response column name.

Returns:

Long DataFrame, one row per (group × predictor).

Return type:

DataFrame

digimuh.stats_production.daily_climate_vs_yield_correlations(df, response='yield_residual')

Spearman rs of daily mean THI / barn-temp vs a yield response.

Parameters:

• df (DataFrame)

• response (str)

Return type:

DataFrame

DIM-adjusted milk yield via per-lactation Wood (1967) curves.

Wood model:

y(t) = a · t^b · exp(-c · t)

fitted in log space by OLS (ln y = ln a + b ln t − c t). Each cow-lactation gets its own (a, b, c); lactations with too few points fall back to a parity-pooled Wood curve. The resulting residuals y − ŷ are approximately symmetric, freeing downstream stratifications (e.g. terciles) from DIM × parity confounding.

References

Wood, P.D.P. (1967) Algebraic model of the lactation curve in

cattle. Nature 216: 164–165. doi:10.1038/216164a0

Wilmink, J.B.M. (1987) Adjustment of test-day milk, fat and

protein yield for age, season and stage of lactation. Livest Prod Sci 16: 335–348.

Macciotta, N.P.P. et al. (2011) Mathematical description of

lactation curves in dairy cattle. Ital J Anim Sci 10: e51.

digimuh.stats_lactation_curve.DIM_MIN: int = 5

Minimum DIM to include (fresh-cow noise dominates 0–4 d).

digimuh.stats_lactation_curve.DIM_MAX: int = 305

Maximum DIM to include (standard 305-day lactation).

digimuh.stats_lactation_curve.MIN_POINTS_PER_LACTATION: int = 30

Fewer than this and we fall back to the parity-pooled curve.

digimuh.stats_lactation_curve.load_calvings(data_dir, db_path=None)

Load one calving-confirmation event per row.

Resolution order:

1. <data_dir>/calvings.csv if present.

2. Query smaxtec_events on db_path and cache to <data_dir>/calvings.csv.

Parameters:

• data_dir (Path) – Directory containing extract-stage CSVs.

• db_path (Path | None) – SQLite DigiMuh database (only used if the CSV is missing).

Returns:

DataFrame with animal_id, calving_date (pandas datetime64).

Return type:

DataFrame

digimuh.stats_lactation_curve.attach_dim(daily_yield, calvings, dim_min=5, dim_max=305)

Attach dim and lactation_nr to each cow-day.

For each cow-day, DIM is the number of days between the most recent prior calving and the observation date. Cow-days with no prior calving on record, or whose DIM is outside [dim_min, dim_max], are dropped.

Parameters:

• daily_yield (DataFrame) – DataFrame with animal_id, date, daily_yield_kg.

• calvings (DataFrame) – Output of load_calvings().

• dim_min (int) – Lower DIM bound (inclusive). Default 5.

• dim_max (int) – Upper DIM bound (inclusive). Default 305.

Returns:

Filtered DataFrame with dim (int) and lactation_nr (1, 2, 3, …) added. Cow-days are dropped if no prior calving exists.

Return type:

DataFrame

digimuh.stats_lactation_curve.fit_wood(dim, yld)

Fit Wood’s incomplete-gamma lactation model via log-linear OLS.

Model:

y = a · t^b · exp(-c · t)
ln y = ln a + b · ln t − c · t

Parameters:

• dim (ndarray) – Days in milk (must be > 0).

• yld (ndarray) – Daily yield (kg/d, must be > 0).

Returns:

Dict with a, b, c, r_squared, n, peak_dim, peak_yield and converged (True when the fit is biologically plausible: b>0, c>0, peak within [5, 120] days).

Return type:

dict

digimuh.stats_lactation_curve.predict_wood(dim, a, b, c)

Evaluate the Wood curve at the given DIMs.

All four arguments broadcast together — scalars or arrays of matching length are accepted. The result is NaN anywhere dim <= 0 or any of the parameters is NaN, evaluated per-element rather than per-array, so a single NaN in the middle of a / b / c does not contaminate the rest of the output.

Return type:

ndarray

digimuh.stats_lactation_curve.fit_wood_per_lactation(df_with_dim, min_points=30)

Fit one Wood curve per (animal_id, calving_date).

Lactations with fewer than min_points valid observations are marked method = "parity_fallback" and will use the parity-pooled curve at prediction time.

Parameters:

• df_with_dim (DataFrame) – Output of attach_dim().

• min_points (int) – Per-lactation point threshold for an individual fit. Default 30.

Returns:

DataFrame of fit results, one row per (animal_id, calving_date), with columns a, b, c, r_squared, n, peak_dim, peak_yield, converged, method, parity.

Return type:

DataFrame

digimuh.stats_lactation_curve.load_daily_yields_for_fitting(data_dir)

Load the full-history daily yield CSV if present.

daily_milk_yield_full.csv is the per-animal full HerdePlus history (no date filter) written by digimuh.extract.extract_daily_milk_yield_full(). It is the recommended fit frame for compute_wood_residuals(): per-lactation Wood fits converge far more often when they see the full pre-peak / post-peak shape instead of a summer slice.

Returns:

DataFrame or None when the file does not exist.

Parameters:

data_dir (Path)

Return type:

DataFrame | None

digimuh.stats_lactation_curve.compute_wood_residuals(daily_yield, calvings, min_points=30, dim_min=5, dim_max=305, fit_yields=None)

End-to-end: attach DIM, fit per-lactation Wood curves, add residuals.

The default behaviour (fit_yields=None) fits Wood curves on the same yield frame used to compute residuals — the original single-frame mode. Passing fit_yields decouples the two: Wood parameters are estimated on fit_yields (typically the full HerdePlus history so each lactation has its pre-peak and post-peak tail visible), then residuals are evaluated on daily_yield (typically the analysis window — for the Frontiers broken-stick pipeline, the Tierauswahl summer slice).

Parameters:

• daily_yield (DataFrame) – The cow-days the caller wants residuals for. Must have animal_id, date, daily_yield_kg, year.

• calvings (DataFrame) – From load_calvings().

• min_points (int) – Per-lactation point threshold for an individual Wood fit (lactations with fewer points use the parity pool).

• dim_min (int) – DIM clipping range applied to both frames.

• dim_max (int) – DIM clipping range applied to both frames.

• fit_yields (DataFrame | None) – Optional frame to estimate Wood parameters on. When given, the returned fits table reports the full-history lactation coverage (typically tens of thousands more points per animal), and yield_expected on each daily_yield cow-day comes from the curve fitted on that richer data. When None, falls back to the legacy single-frame behaviour.

Returns:

(yields, fits) where yields is the daily_yield cow-day DataFrame enriched with dim, lactation_nr, parity, yield_expected, yield_residual, yield_residual_rel (residual / expected), and method ("per_lactation" / "parity_fallback"). fits is the per-lactation fit table.

Return type:

tuple[DataFrame, DataFrame]

digimuh.stats_lactation_curve.residual_terciles(residuals)

Return the (Q33, Q67) boundaries of a residual series.

Parameters:

residuals (Series)

Return type:

tuple[float, float]

Milk composition vs climate analysis on MLP test-days.

The MLP channel (§3.20) supplies once-a-month per-cow test-day rows with fat %, protein %, lactose, fat-kg, ECM, F/E ratio, SCC and urea. Matching these to the cow-day climate from tnf_yield.csv lets us ask whether heat exposure shifts milk composition in the direction the “thin milk” hypothesis predicts:

This module exposes reusable helpers so the milk-yield classifier can orchestrate the join, correlation and plot without re-implementing the logic.

digimuh.stats_milk_composition.load_mlp_test_days(data_dir)

Load mlp_test_days.csv and add a day-resolution date column.

Parameters:

data_dir (Path)

Return type:

DataFrame

digimuh.stats_milk_composition.merge_mlp_with_cowday(mlp, wood, tnf_yield, class_table=None, rumen=None)

Join an MLP test-day table to the cow-day climate + residual frames.

The (animal_id, date) join is exact: MLP timestamps land on the test-day itself. Rows with no matching cow-day (e.g. MLP test-days outside the Tierauswahl summer window) are dropped — there is no climate to correlate against for those.

Parameters:

• mlp (DataFrame) – mlp_test_days.csv DataFrame.

• wood (DataFrame) – daily_milk_yield_wood.csv DataFrame. Provides yield_residual, dim, year.

• tnf_yield (DataFrame) – tnf_yield.csv DataFrame. Provides thi_tnf, temp_tnf, mean_thi (and mean_barn_temp when available).

• class_table (DataFrame | None) – Optional yield_class_per_cow_year.csv frame (from stats_production.classify_cow_years_by_wood_residual()); when supplied, the output carries a yield_class column.

• rumen (DataFrame | None)

Returns:

DataFrame with one row per matched MLP test-day and all MLP, climate, residual and (optionally) class columns.

Return type:

DataFrame

digimuh.stats_milk_composition.mlp_climate_correlations(df, responses=(('herdeplus_mlp_mkg', 'Test-day milk', 'kg/d'), ('herdeplus_mlp_fat_percent', 'Fat %', '%'), ('herdeplus_mlp_protein_percent', 'Protein %', '%'), ('herdeplus_mlp_lactose', 'Lactose %', '%'), ('herdeplus_mlp_fkg', 'Fat', 'kg/d'), ('herdeplus_mlp_ekg_percent', 'Protein (Eiweiß)', 'kg/d'), ('herdeplus_mlp_lkg', 'Lactose', 'kg/d'), ('herdeplus_mlp_f_e', 'F/E ratio', '—'), ('herdeplus_mlp_ecm', 'ECM (energy-corrected milk)', 'kg/d'), ('herdeplus_mlp_cell_count', 'Somatic cell count', '10³/mL'), ('herdeplus_mlp_urea', 'Urea', 'mg/dL')), predictors=(('mean_thi', 'Daily mean barn THI'), ('mean_barn_temp', 'Daily mean barn temperature (°C)'), ('thi_tnf', 'Daily time below THI breakpoint (TNF)')))

Spearman rs + p + OLS slope per (group × predictor × response).

group is "pooled" by default; when the input has a yield_class column each of its categories is added as a separate group row.

Returns:

Long DataFrame with columns group, predictor, response, response_label, unit, n, n_animals, rs, p, slope, intercept.

Parameters:

• df (DataFrame)

• responses (tuple[tuple[str, str, str], ...])

• predictors (tuple[tuple[str, str], ...])

Return type:

DataFrame

digimuh.stats_milk_composition.compute_dilution_partition(merged)

Partition each composition response into dilution + rumen effects.

For every cow-day row we compute what the fat % and protein % would be if the cow produced her personal reference amount of fat and protein (in kg) but diluted them into the actually- observed volume. Deviations from that dilution prediction tell us how much of the observed composition drop on hot days is explained by added water alone versus how much extra suppression comes from reduced rumen output.

Per-cow references are the animal’s own mean test-day fat_kg / protein_kg / volume across all of her MLP rows in the merged frame — so “dilution-predicted” is self-referential to the cow and does not borrow composition across animals.

Added columns:

• volume_ratio = observed_volume / ref_volume

• fat_percent_diluted = fat_kg_ref / observed_volume × 100

• protein_percent_diluted = protein_kg_ref / observed_volume × 100

• fat_percent_rumen = observed − dilution-predicted fat %

• ``protein_percent_rumen``= observed − dilution-predicted protein %

• fat_kg_ref, protein_kg_ref, volume_ref (the per-cow references, repeated on every row)

Parameters:

merged (DataFrame)

Return type:

DataFrame

digimuh.stats_milk_composition.dilution_partition_summary(df, predictor='mean_thi')

Compare observed vs dilution-predicted composition slopes.

For each of fat % and protein %, computes three correlations (rs + OLS slope) against the climate predictor:

1. observed — the actual composition drop

2. dilution-predicted — the drop implied by the cow’s own

   reference fat/protein kg diluted into the observed volume

3. rumen residual — observed − dilution-predicted;

   positive means rumen over-produced relative to the reference, negative means rumen was suppressed below the reference

A cleanly dilutive explanation has a near-zero rumen-residual slope and a dilution-predicted slope that matches the observed one. A suppression-dominant story has an observed slope noticeably steeper than the dilution-only one.

Returns a long DataFrame with columns nutrient ("fat" or "protein"), component ("observed" / "dilution_predicted" / "rumen_residual"), rs, p, slope, n.

Parameters:

• df (DataFrame)

• predictor (str)

Return type:

DataFrame

digimuh.stats_milk_composition.thin_milk_verdict(correlations, predictor='mean_thi')

Distil the MLP × climate table into a one-line verdict.

The “thin milk” hypothesis predicts, under a hotter cow-day:

• milk volume (mkg) rises with climate — positive rs

• fat % and protein % drop — negative rs

• F/E ratio drops — negative rs

Returns a dict with the pooled rs for each component plus a single "verdict" string ("supported" / "partial" / "refuted").

Parameters:

• correlations (DataFrame)

• predictor (str)

Return type:

dict

Longitudinal breakpoint stability analysis.

digimuh.stats_longitudinal.compute_breakpoint_icc(bs_results, wood_yield=None, production=None)

One-way random-effects ICC(1,1) on multi-year breakpoints.

For each predictor in {THI, barn temperature} we run two flavours of the ICC:

• raw — directly on the converged breakpoint values.

• residual — on OLS residuals from regressing the breakpoint on parity bucket (1/2/3+) plus mean DIM across that cow-year’s summer window. Requires daily_milk_yield_wood.csv-style frame; silently skipped when not available.

The cohort for each row is the subset of animals with at least two converged breakpoints for that predictor.

Returns one row per (predictor, mode) with point estimate, 95% CI, F, df, p, n_animals, n_obs, mean k (the harmonic- mean k0 used by Shrout & Fleiss for unbalanced designs).

Parameters:

• bs_results (DataFrame)

• wood_yield (DataFrame | None)

• production (DataFrame | None)

Return type:

DataFrame

digimuh.stats_longitudinal.compute_stability(bs_results)

ICC and paired data for repeat animals.

Returns:

(pairs DataFrame, ICC value).

Parameters:

bs_results (DataFrame)

Return type:

tuple[DataFrame, float]

digimuh.stats_longitudinal.make_summary_table(bs)

Generate Table 1 for the manuscript.

Parameters:

bs (DataFrame)

Return type:

DataFrame

Year-round milk-yield statistics for the large Neubau cohort.

digimuh.stats_annual_yield.MIN_DAYS_PER_COW_YEAR: int = 30

Fewer milking days than this in a calendar year and the per-cow-year z-score (mean/SD) is too unstable to trust — that cow-year is dropped.

digimuh.stats_annual_yield.DOY_BIN_DAYS: int = 7

Day-of-year bin width for the median surfaces (weekly).

digimuh.stats_annual_yield.DIM_BIN_DAYS: int = 15

DIM bin width for the (day-of-year × DIM) surface.

digimuh.stats_annual_yield.LACTATION_CAP: int = 8

Lactation numbers ≥ this are pooled into one top bin (sparse tail).

digimuh.stats_annual_yield.MIN_COWS_PER_CELL: int = 3

Surface cells backed by fewer cows than this are left empty (NaN).

digimuh.stats_annual_yield.BASELINE_LOOKBACK_DAYS: int = 3

How far back to search for a valid pre-run baseline yield when the calendar day immediately before a heat-stress run has no milking record.

digimuh.stats_annual_yield.compute_annual_zscores(daily_full, calvings, dim_min=5, dim_max=305, min_days=30)

Z-score each cow’s daily yield within each calendar year.

For every (animal, calendar-year) the daily yield is standardised to that cow-year’s own mean and SD, so cows of different production level and lactation stage become comparable before pooling. DIM and lactation number are attached for the downstream surfaces.

Parameters:

• daily_full (DataFrame) – Full-history daily yield (animal_id, date, daily_yield_kg).

• calvings (DataFrame) – From stats_lactation_curve.load_calvings().

• dim_min (int) – DIM clipping (defaults 5–305).

• dim_max (int) – DIM clipping (defaults 5–305).

• min_days (int) – Minimum milking days in a cow-year to keep it.

Returns:

DataFrame with animal_id, year, date, doy, dim, lactation_nr, daily_yield_kg, z_yield.

Return type:

DataFrame

digimuh.stats_annual_yield.build_median_surface(z_df, y_col, value_col, doy_bin=7, dim_bin=15, min_cows=3, n_boot=2000, seed=12345)

Median-across-cows yield per (day-of-year bin × y bin) cell.

Each cow contributes once per cell (her mean value_col in that cell), so cows with many days do not dominate; the cell value is the median over those per-cow means — “the median cow’s yield”. A bootstrap 95% CI of that median (resampling over cows) is attached so the same surface can be drawn as a heatmap or as CI-banded line plots.

Parameters:

• z_df (DataFrame) – Output of compute_annual_zscores().

• y_col (str) – "lactation_nr" or "dim".

• value_col (str) – "z_yield" or "daily_yield_kg".

• doy_bin (int) – Day-of-year bin width.

• dim_bin (int) – DIM bin width (only used when y_col == "dim").

• min_cows (int) – Cells with fewer backing cows are set NaN.

• n_boot (int) – Bootstrap resamples for the per-cell median CI.

• seed (int) – RNG seed (reproducible).

Returns:

doy_bin (left edge), y_bin, median_value, ci_lo, ci_hi, n_cows, n_obs.

Return type:

Long-form DataFrame

digimuh.stats_annual_yield.compute_heatstress_duration_deltas(yield_wood, breakpoints, climate_daily, lookback=3)

Yield change across runs of consecutive heat-stress days.

A day is heat-stress for a cow when the barn’s daily-maximum THI reaches or exceeds that cow’s individual (converged) THI breakpoint. Within each cow-year we find maximal runs of consecutive calendar days that are all heat-stress, take the last non-stress day before the run as the baseline, and record the yield change on the 1st, 2nd, 3rd … day of the run.

The primary response is the Wood residual (DIM-adjusted), so the lactation decline over the spell is not mistaken for a heat effect; raw kg is carried alongside.

Parameters:

• yield_wood (DataFrame) – Summer Wood frame (animal_id, date, year, dim, daily_yield_kg, yield_residual).

• breakpoints (DataFrame) – broken_stick_results with animal_id, year, thi_breakpoint, thi_converged.

• climate_daily (DataFrame) – day, barn_thi_max (herd-level barn climate).

• lookback (int) – Max days to search back for a valid baseline yield.

Returns:

animal_id, year, run_id, streak_day, date, dim, daily_yield_kg, yield_residual, baseline_residual, baseline_kg, delta_residual, delta_kg, run_length.

Return type:

Long DataFrame, one row per (cow-year run-day with valid yield)

digimuh.stats_annual_yield.aggregate_per_cow(deltas, max_streak=7)

Collapse run-day deltas to one median per cow per streak-day.

Stage one of the per-cow → inter-cow estimator: each cow contributes her own median delta at each streak position, so a cow with many heat-stress runs does not dominate the herd-level summary. Feed the result to summarise_duration_deltas() to obtain the inter-cow median (median across cows) with a cow-level bootstrap CI.

Parameters:

• deltas (DataFrame) – Output of compute_heatstress_duration_deltas().

• max_streak (int) – Keep streak days 1…max_streak.

Returns:

animal_id, streak_day, delta_residual, delta_kg (each the cow’s median at that streak day), n_obs.

Return type:

DataFrame

digimuh.stats_annual_yield.summarise_duration_deltas(deltas, max_streak=7, n_boot=20000, seed=12345)

Median yield delta per streak-day with bootstrap 95% CI.

Operates on whatever rows it is given: pass the raw run-day deltas for the pooled median, or the aggregate_per_cow() table for the inter-cow median (median of per-cow medians, bootstrap over cows).

Parameters:

• deltas (DataFrame) – Run-day deltas, or per-cow medians from aggregate_per_cow().

• max_streak (int) – Report streak days 1…max_streak (deeper runs pooled out — the sample thins quickly).

• n_boot (int) – Bootstrap resamples for the median CI.

• seed (int) – RNG seed (reproducible).

Returns:

streak_day, n, n_cows, median_delta_residual, ci_lo_residual, ci_hi_residual, median_delta_kg, ci_lo_kg, ci_hi_kg.

Return type:

DataFrame

digimuh.stats_annual_yield.test_duration_medians(per_cow, max_streak=7)

Wilcoxon signed-rank test that each streak-day median ≠ 0, BH-FDR.

One test per streak day, on the per-cow median deltas (one value per cow — independent observations), asking whether the centre of that per-cow distribution differs from zero. Raw p-values are Benjamini-Hochberg FDR corrected within each response metric (the seven streak-day tests form the family).

Parameters:

• per_cow (DataFrame) – Output of aggregate_per_cow() (per-cow medians).

• max_streak (int) – Number of streak days tested (the test family size).

Returns:

metric (“residual”/”kg”), streak_day, n_cows, median, wilcoxon_stat, p_raw, p_fdr, stars.

Return type:

DataFrame

CLI orchestration for the broken-stick statistical pipeline.

Usage:

digimuh-stats --data results/broken_stick

digimuh.stats_runner.main()

Return type:

None

Visualisation

Shared matplotlib configuration and figure I/O for all plots.

digimuh.viz_base.setup_figure()

Configure matplotlib for publication-quality figures.

Return type:

None

digimuh.viz_base.save_figure(fig, name, out_dir)

Save figure as SVG + PNG and close.

The target folder is resolved through digimuh.paths.resolve_output(), which routes the file into a subject-specific subfolder (e.g. 03_temporal) when the stem is known. Callers stay unaware of the subfolder mapping — they pass the same name / out_dir as before.

Parameters:

• name (str)

• out_dir (Path)

Return type:

None

digimuh.viz_base.add_significance_bracket(ax, x1, x2, y, stars, h=0.02, lw=1.2)

Draw a significance bracket with stars between two x positions.

Parameters:

• ax – Matplotlib axes.

• x1 (float) – Left and right x positions.

• x2 (float) – Left and right x positions.

• y (float) – Y position of the bracket bottom (data coords).

• stars (str) – Text to display (e.g. ‘***’, ‘n.s.’).

• h (float) – Bracket height as fraction of y-range.

• lw (float) – Line width.

Return type:

None

Breakpoint analysis figures for the Frontiers manuscript.

digimuh.viz_breakpoints.plot_grouped_boxplots(bs, out_dir)

Side-by-side boxplots: rumen vs respiration breakpoints per year, with within-year Wilcoxon signed-rank tests (BH-FDR corrected).

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_paired_below_above(beh, tests, out_dir)

Paired boxplots: below vs above breakpoint with significance brackets.

Y-axes are unified across all years (and across panels of the same response variable) so per-year figures are directly comparable when laid out side by side: the body-temperature panel uses the same [min, max] in 2021, 2022, 2023, 2024, and the respiration panel uses its own shared range. The shared limits also re-anchor the significance brackets so they sit just below the top of every panel rather than tracking each year’s local data range.

Parameters:

• beh (DataFrame)

• tests (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_paired_rumen_vs_resp(bs, out_dir)

Paired boxplot: rumen temp breakpoint vs respiration breakpoint.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_spearman(spearman, out_dir)

Spearman correlation distributions.

Parameters:

• spearman (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_climate(climate, out_dir)

Daily barn THI and temperature time series per summer.

Parameters:

• climate (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_predictors(bs, out_dir)

Scatter plots: breakpoint vs production parameters.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_stability(pairs, icc, out_dir)

Year-to-year breakpoint stability scatter.

Parameters:

• pairs (DataFrame)

• icc (float)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_examples(rumen, resp, bs, out_dir, *, select_top=None, show_hill=True, show_davies=True, predictors=None)

Example broken-stick panels per signal/predictor combination.

Two selection modes:

• Scenario mode (select_top=None, default): 3 diagnostic panels per predictor illustrating (A) BS+Hill converged, (B) BS failed / Hill rescues, (C) Davies n.s. — requires Davies/Hill columns.

• Top-N mode (select_top=N): renders the N BS-converged animal-years with highest R² per predictor in a grid. Useful for picking publication example fits.

Parameters:

• rumen (DataFrame) – Rumen + barn climate data from extract stage.

• resp (DataFrame) – Respiration + barn climate data (pass an empty DataFrame to skip respiration predictors, e.g. in --no-resp mode).

• bs (DataFrame) – Broken-stick results DataFrame.

• out_dir (Path) – Output directory for figures.

• select_top (int | None) – If given, switch to top-N mode and render N panels of highest-R² BS-converged animal-years per predictor.

• show_hill (bool) – If False, skip Hill fit overlay and onset line. Set False for Frontiers figures (BS-only).

• show_davies (bool) – If False, skip Davies/Pscore annotation.

• predictors (tuple[str, ...] | None) – Optional subset of prefixes to include ("thi", "temp", "resp_thi", "resp_temp"). Defaults to all four.

Return type:

None

digimuh.viz_breakpoints.plot_bodytemp_vs_resp_scatter(bs, out_dir)

Scatter: rumen temp THI breakpoint vs resp THI breakpoint.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_breakpoints.plot_thi_vs_temp_scatter(bs, out_dir)

Scatter: THI breakpoint vs barn temp breakpoint.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

Correlation and event-triggered figures for the Frontiers manuscript.

digimuh.viz_correlation.plot_cross_correlation(out_dir)

Plot cross-correlation AND cross-covariance curves below vs above bp.

Error bands are standard error of the mean (SEM). Peak lag annotated with vertical coloured line and horizontal arrow from lag=0 to peak, labelled ‘rumen temperature N min later’.

Parameters:

out_dir (Path)

Return type:

None

digimuh.viz_correlation.plot_derivative_ccf(out_dir)

Plot derivative CCF: d(climate)/dt vs d(body_temp)/dt.

Parameters:

out_dir (Path)

Return type:

None

digimuh.viz_correlation.plot_event_triggered_average(out_dir, traces_file='event_triggered_traces.csv', suffix='', title_extra='')

Plot peri-event average of rumen temp around breakpoint crossings.

Three panels per predictor: A) Climate signal (THI or barn temp) aligned to crossing B) Raw rumen temperature C) Rumen temperature baseline-subtracted (acute onset)

Parameters:

• out_dir (Path) – Output directory.

• traces_file (str) – Name of the traces CSV file.

• suffix (str) – Appended to output filename (e.g. ‘_filtered’).

• title_extra (str) – Appended to figure title (e.g. ‘ (8-11h crossings)’).

Return type:

None

digimuh.viz_correlation.plot_climate_eta(out_dir)

Climate signal around breakpoint crossings, normalised to breakpoint.

Two figures: 1. THI-triggered crossings: left y = THI − THI_breakpoint,

right y = barn temperature (raw °C)

2. Barn temp-triggered crossings: left y = barn_temp − temp_breakpoint, right y = THI (raw)

In both cases y=0 on the left axis is the breakpoint threshold.

Parameters:

out_dir (Path)

Return type:

None

Production impact figures for the Frontiers manuscript.

digimuh.viz_production.plot_tnf_yield(out_dir)

Scatter plots: daily thermoneutral fraction vs daily milk yield.

Each dot is one cow-day. Panel A shows absolute yield, Panel B shows relative yield (daily yield / cow-specific P95).

Parameters:

out_dir (Path)

Return type:

None

digimuh.viz_production.plot_tnf_yield_by_class(tnf_by_class, correlations, out_dir, response='daily_yield_kg')

Scatter grid — 3 yield classes × 2 TNF predictors.

Each panel shows cow-day points coloured by class, an OLS reference line, and a Spearman annotation taken directly from the correlations table so the figure and the console summary never disagree.

Parameters:

• tnf_by_class (DataFrame) – Output of stats_production.compute_tnf_yield_by_class().

• correlations (DataFrame) – Output of stats_production.tnf_yield_correlations_by_class().

• out_dir (Path) – Figure output directory.

• response (str) – Column name for the y-axis — usually "daily_yield_kg" (default) or "relative_yield".

Return type:

None

digimuh.viz_production.plot_crossing_day_raincloud(df, comparison, out_dir, response='yield_residual', group='pooled')

Raincloud pair per climate predictor — crossed-day vs not.

One figure, two panels: THI crossings (left), barn-temp crossings (right). For each, two rainclouds compare the response on cow-days with ≥1 breakpoint crossing of that predictor vs cow-days with none.

Parameters:

• df (DataFrame) – Cow-day DataFrame with thi_crossed, temp_crossed, the response column, and (if filtering) yield_class.

• comparison (DataFrame) – Output of stats_production.crossing_day_comparison() with the same response. Used to pull the Mann-Whitney annotation.

• out_dir (Path) – Figure output directory.

• response (str) – Response column, default yield_residual.

• group (str) – Which class row to plot; "pooled" or one of "low", "middle", "high".

Return type:

None

digimuh.viz_production.plot_daily_climate_vs_yield(df, correlations, out_dir, response='yield_residual')

Scatter of daily-mean THI / barn-temp vs a yield response.

One figure, two panels — mean THI left, mean barn-temp right. A single OLS fit across ALL cow-days per panel; the Spearman rs, p and slope come from correlations so figure and console agree.

Parameters:

• df (DataFrame)

• correlations (DataFrame)

• out_dir (Path)

• response (str)

Return type:

None

Figures for the MLP composition × climate analysis.

Two figures:

• mlp_thin_milk_hypothesis.{svg,png} — a compact 2×2 scatter directly testing the hypothesis: milk volume and ECM on the top row (expected to rise with climate if dilution is real), fat % and protein % on the bottom row (expected to fall). Each panel has an OLS reference line and a Spearman annotation.

• mlp_composition_heatmap.{svg,png} — a signed-rs heatmap across every MLP response × every climate predictor × each Wood-residual yield class, so patterns in SCC / urea / lactose and class-specific differences come out at a glance.

digimuh.viz_milk_composition.plot_thin_milk_hypothesis(merged, correlations, out_dir, predictor='mean_thi', predictor_label='Daily mean barn THI')

2×2 scatter panel directly testing the dilution / thin-milk story.

Parameters:

• merged (DataFrame)

• correlations (DataFrame)

• out_dir (Path)

• predictor (str)

• predictor_label (str)

Return type:

None

digimuh.viz_milk_composition.plot_dilution_partition(df, summary, out_dir, predictor='mean_thi', predictor_label='Daily mean barn THI')

Overlay observed vs dilution-predicted fat % / protein %.

Two panels (fat % top, protein % bottom) showing three OLS reference lines per panel:

• observed (vermillion) — the actual drop

• dilution predicted (blue, dashed) — the drop the

  volume increase alone can explain if the cow’s reference fat/protein kg were held constant

• rumen residual (grey, dotted) — observed − dilution;

  how much extra drop comes from reduced absolute fat/protein output

A zero-rs rumen residual means pure dilution; a negative residual means rumen suppression on top of dilution.

Parameters:

• df (DataFrame)

• summary (DataFrame)

• out_dir (Path)

• predictor (str)

• predictor_label (str)

Return type:

None

digimuh.viz_milk_composition.plot_composition_heatmap(correlations, out_dir, predictor='mean_thi')

Signed-rs heatmap: rows = MLP metrics, columns = yield class.

Parameters:

• correlations (DataFrame)

• out_dir (Path)

• predictor (str)

Return type:

None

Longitudinal and pipeline-summary figures.

digimuh.viz_longitudinal.plot_longitudinal_breakpoints(bs, out_dir)

Track how individual breakpoints change across years.

Only animals present in 2+ years are included. Shows both absolute breakpoints and change relative to the animal’s first year.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_longitudinal.plot_breakpoint_raincloud(out_dir)

Raincloud plots of annual breakpoint crossing counts per animal.

For each animal-year, counts how many times the barn THI (or barn temp) crossed that animal’s individual breakpoint upward. Reads crossing events from crossing_times.csv (produced by stats).

Each year gets a horizontal raincloud: half-violin (KDE), jittered scatter of integer counts, and boxplot. Shows whether cows experience more threshold exceedances over time (climate signal).

Parameters:

out_dir (Path)

Return type:

None

digimuh.viz_longitudinal.plot_breakpoint_value_raincloud(bs, out_dir)

Raincloud plots of the individual breakpoint value per year.

Companion to plot_breakpoint_raincloud(), which shows crossing counts. Here each converged animal contributes its fitted THI (or barn-temp) breakpoint for the year, so the rows show how the herd’s threshold distribution shifts across 2021-2024. Same per-year colour code and Kruskal-Wallis test as the crossing-count raincloud.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_longitudinal.plot_breakpoint_retest(bs, out_dir)

This-summer vs last-summer breakpoint scatter (predictability).

Third companion to the per-year value and crossing-count rainclouds: each point is one cow’s consecutive-summer pair (x = last summer, y = this summer), coloured by the later year. The identity line is perfect stability; the flat fitted line and near-zero r show the breakpoint does not carry over between years.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_longitudinal.plot_breakpoint_icc(out_dir, csv_stem='breakpoint_icc', figure_stem='breakpoint_icc_forest', title_suffix='')

Forest plot of one-way ICC(1,1) for multi-year breakpoints.

Reads the cohort-specific ICC CSV (default breakpoint_icc.csv) and renders a forest with point estimates, 95% CIs, the n / k label next to each row, and three reference lines: ICC = 0 (no repeatability — within = between), ICC = 0.5 (“moderate” repeatability, Koo & Li 2016), and ICC = 0.75 (“good”). Pass csv_stem='breakpoint_icc_strict' and figure_stem='breakpoint_icc_strict_forest' to render a cohort-suffixed variant.

Parameters:

• out_dir (Path)

• csv_stem (str)

• figure_stem (str)

• title_suffix (str)

Return type:

None

digimuh.viz_longitudinal.plot_longitudinal_sankey(bs, out_dir)

Alluvial plot tracking breakpoint adaptation across years.

Closed cohort: only animals with a converged breakpoint in every study year are included. Column 0 is the baseline year (all animals start as “stable”). Subsequent columns show the year- to-year Δ category (strongly decreased / decreased / stable / increased / strongly increased). Bezier bands track individual animals between columns.

Duplicate entries per animal-year (multiple date_enter) are resolved by averaging breakpoints before computing deltas.

Parameters:

• bs (DataFrame)

• out_dir (Path)

Return type:

None

digimuh.viz_longitudinal.plot_threshold_sankey(bs, out_dir)

Plotly Sankey diagrams showing how animals flow through the threshold detection pipeline: Davies/pscore → broken-stick → Hill.

Produces one figure for rumen temperature vs THI and one for respiration rate vs THI.

Parameters:

• bs (DataFrame) – broken_stick_results.csv DataFrame.

• out_dir (Path) – Output directory.

Return type:

None

Temporal analysis figures: circadian and crossing activation.

digimuh.viz_temporal.plot_circadian_null_model(out_dir)

Plot the rumen-temperature circadian null model (2×2 grid).

Four panels, all plotted over clock hour on cool vs stress days (stratified by each animal’s own THI breakpoint):

• A — Rumen temperature, with THI-crossing density overlay.

• B — Stress − cool Δ rumen temperature (heat-stress dose), with THI-crossing density overlay.

• C — Barn THI on cool vs stress days + THI-crossing density.

• D — Barn temperature on cool vs stress days + barn-temp-crossing density.

Shaded regions:

• Coloured ribbon around each mean line = ±1 SEM across animals (cool=blue, stress=vermillion, Δ=pink).

• Grey vertical bands at 04:00–07:00 and 16:00–19:00 = excluded milking windows.

• Green ribbon on the right y-axis = kernel density estimate of breakpoint-crossing clock times (smoothed histogram of when, during the 24 h cycle, each cow’s individual breakpoint is crossed upward).

Grid lines are disabled to keep the overlays readable. A companion figure circadian_null_model_stacked is written with the same variables on rows sharing a single hour-of-day x-axis.

Parameters:

out_dir (Path)

Return type:

None

digimuh.viz_temporal.plot_thi_daily_profile(out_dir)

Plot barn THI and barn temperature across 24h by month.

Two stacked panels per year (top: barn temperature, bottom: THI), so the reader can compare the two heat-load axes side by side. Months coloured consistently across panels.

Parameters:

out_dir (Path)

Return type:

None

digimuh.viz_temporal.plot_crossing_raster(out_dir)

Raster plot of breakpoint crossing events across the 24h cycle.

Left panel: activation raster — each row is one cow (sorted by breakpoint value), each dot is a crossing event at that clock time. Dot colour = breakpoint value.

Right panel: raincloud — half-violin (KDE) + jittered scatter + boxplot of crossing clock times.

Parameters:

out_dir (Path)

Return type:

None

Year-round yield line plots and the heat-stress duration-response figure.

digimuh.viz_annual_yield.plot_yield_lines(grid, value_kind, y_kind, name, out_dir)

Median-cow yield over the year as CI-banded lines per band.

One line per y-band (median yield vs day-of-year) with a bootstrap 95% CI ribbon. Lactation bands get the discrete Wong palette and a legend; DIM bands get a viridis gradient and a colourbar (too many bands for a legend).

Parameters:

• grid (DataFrame) – Long-form grid from stats_annual_yield.build_median_surface() (doy_bin, y_bin, median_value, ci_lo, ci_hi, n_cows).

• value_kind (str) – "z" (z-scored, zero reference line) or "kg" (raw kg).

• y_kind (str) – "lactation" or "dim" — sets band colouring.

• name (str) – Output stem (also the CSV companion stem).

• out_dir (Path) – Results directory (routed to 07_annual_yield).

Return type:

None

digimuh.viz_annual_yield.plot_heatstress_duration(summary, deltas, value_tag, name, out_dir, tests=None)

Yield delta vs consecutive heat-stress days: median+CI over rainclouds.

Parameters:

• summary (DataFrame) – From stats_annual_yield.summarise_duration_deltas().

• deltas (DataFrame) – Per-cow medians from stats_annual_yield.aggregate_per_cow().

• value_tag (str) – "residual" (DIM-adjusted) or "kg".

• name (str) – Output stem.

• out_dir (Path) – Results directory.

• tests (DataFrame | None) – Optional per-streak-day significance table from stats_annual_yield.test_duration_medians() (already filtered to this metric); annotates BH-FDR stars on the median.

Return type:

None

CLI orchestration for figure generation.

Usage:

digimuh-plots --data results/broken_stick

digimuh.viz_runner.main()

Return type:

None

Pipelines & Orchestrators

End-to-end runner for the Frontiers in Animal Science paper.

Calls extract, stats, and plots in sequence. Equivalent to the old scripts/run_00_broken_stick_ana.sh but as a Python entry point with proper argument forwarding.

digimuh.run_frontiers_2026.main()

Return type:

None

CLI: year-round milk-yield analysis (z-scores, surfaces, duration).

digimuh.run_annual_yield.main()

Return type:

None

Milk yield classification analysis for the Frontiers manuscript.

Reads daily_milk_yield.csv from the extract stage and compares its pooled distribution against two published tercile schemes and the herd’s own 33rd / 67th percentiles. Writes two figures (pooled and per-year) plus a console summary.

Usage:

digimuh-milk-yield --data results/broken_stick

References

Müschner-Siemens et al. (2020) J Thermal Biol 88:102484. Yan et al. (2021) J Thermal Biol 100:103041.

digimuh.milk_yield_classification.MUESCHNER_SIEMENS: tuple[float, float] = (28.8, 38.4)

terciles of their German Holstein cohort at ATB Potsdam. Same region / breed as our herd.

Type:

Müschner-Siemens et al. (2020)

digimuh.milk_yield_classification.YAN: tuple[float, float] = (26.0, 39.0)

826-cow Chinese Holstein cohort.

Type:

Yan et al. (2021)

digimuh.milk_yield_classification.classify(values, low_bd, middle_bd)

Bucket a yield series into low / middle / high.

Parameters:

• values (Series) – Daily milk yield (kg/d) as a pandas Series.

• low_bd (float) – Upper bound of the “low” bucket (inclusive).

• middle_bd (float) – Upper bound of the “middle” bucket (inclusive).

Returns:

Categorical series of labels "low", "middle", "high".

Return type:

Series

digimuh.milk_yield_classification.load_daily_yields(data_dir)

Load and clean daily_milk_yield.csv.

Drops rows with non-positive or NaN daily_yield_kg.

Parameters:

data_dir (Path) – Directory containing the extract-stage CSVs.

Returns:

DataFrame with at least animal_id, year, date, daily_yield_kg.

Return type:

DataFrame

digimuh.milk_yield_classification.print_summary(df, our_terciles)

Print per-year, pooled, and class-count tables to the console.

Uses the rich console helpers if available, otherwise plain text.

Parameters:

• df (DataFrame)

• our_terciles (tuple[float, float])

Return type:

None

digimuh.milk_yield_classification.plot_pooled_histogram(df, our_terciles, out_dir)

Single-panel histogram of all animal-days with all three schemes.

Parameters:

• df (DataFrame)

• our_terciles (tuple[float, float])

• out_dir (Path)

Return type:

None

digimuh.milk_yield_classification.plot_per_year_histogram(df, our_terciles, out_dir)

Per-year histograms stacked vertically with a shared x-axis.

Parameters:

• df (DataFrame)

• our_terciles (tuple[float, float])

• out_dir (Path)

Return type:

None

digimuh.milk_yield_classification.print_wood_summary(wood, fits)

Pretty-print Wood-curve coverage, fit stats, and residual terciles.

Parameters:

• wood (DataFrame)

• fits (DataFrame)

Return type:

None

digimuh.milk_yield_classification.plot_residual_histogram(wood, out_dir)

Histogram of Wood residuals with tercile and zero lines.

Returns the (Q33, Q67) residual tercile boundaries so the caller can persist them or compare against the raw-yield terciles.

Parameters:

• wood (DataFrame)

• out_dir (Path)

Return type:

tuple[float, float]

digimuh.milk_yield_classification.plot_wood_example_fits(wood, fits, out_dir, n_examples=9, fit_points=None)

Grid of example per-lactation Wood fits with raw data overlay.

Parameters:

• wood (DataFrame) – daily_milk_yield_wood.csv (summer-window residual frame).

• fits (DataFrame) – wood_curve_fits.csv (per-lactation parameters).

• out_dir (Path) – Figure output directory.

• n_examples (int) – Number of panels to draw.

• fit_points (DataFrame | None) – Optional full-history cow-days with dim attached (e.g. the output of attach_dim on daily_milk_yield_full.csv). When supplied, the scatter uses these rows so the Wood curve sits on top of all the cow-days that actually informed the fit rather than the much smaller summer-window slice.

Return type:

None

digimuh.milk_yield_classification.main()

Return type:

None

Standalone Analyses

Subclinical ketosis risk scoring from multi-sensor fusion.

Uses the v_analysis_ketosis view which joins daily milking data (HerdePlus MLP test-day results), smaXtec rumen metrics, water intake, and disease ground truth.

The analysis:

1. Extracts days with MLP test-day data (FPR available).

2. Computes a composite ketosis risk score from: - Fat-to-protein ratio (FPR > 1.4 = energy deficit) - Rumination index (lower = reduced feed intake) - Milk yield deviation from rolling cow mean - Rumen pH (low pH + high FPR = metabolic confusion)

3. Validates against disease records (ground truth).

4. Trains a Random Forest classifier and reports feature importance and cross-validated performance.

Usage:

python -m digimuh.analysis_01_ketosis --db cow.db --out results/ketosis

References

Oetzel (2013) — FPR thresholds for subclinical ketosis. Kaufman et al. (2016) J Dairy Sci 99:5604–18 — rumination

time association with subclinical ketosis.

digimuh.analysis_01_ketosis.load_ketosis_data(con)

Load ketosis analysis view and add derived features.

Parameters:

con – Active database connection with views created.

Returns:

DataFrame with one row per animal per MLP test day, enriched with rolling-mean deviations and risk scores.

Return type:

DataFrame

digimuh.analysis_01_ketosis.train_ketosis_classifier(df, out_dir)

Train a Random Forest to classify ketosis risk.

Uses FPR > 1.4 as the positive label (subclinical ketosis indicator) and evaluates against disease records where available.

Parameters:

• df (DataFrame) – DataFrame from load_ketosis_data().

• out_dir (Path) – Directory for saving results.

Returns:

Dict with performance metrics and feature importances.

Return type:

dict

digimuh.analysis_01_ketosis.plot_ketosis_overview(df, out_dir)

Generate overview plots for ketosis analysis.

Parameters:

• df (DataFrame) – DataFrame from load_ketosis_data().

• out_dir (Path) – Output directory for figures.

Return type:

None

digimuh.analysis_01_ketosis.main()

Entry point for ketosis analysis.

Return type:

None

Per-animal heat stress response modelling.

Uses the v_analysis_heat_stress view which combines daily smaXtec rumen data, DWD weather, gouna respiration, and milking production.

The analysis:

1. Builds per-animal Z-scored rumen temperature following the approach of the NZ smaXtec study (JDS Communications, 2024): each cow’s temperature distribution is scaled to a common mean and SD before thresholding.

2. Fits per-animal thermoregulatory dose-response curves: rumen_temp_z = f(THI) using sigmoidal regression. The inflection point and slope characterise each animal’s heat tolerance.

3. Computes a daily heat load index fusing rumen temp, respiration rate, activity suppression, and water intake.

4. Quantifies the production impact: milk yield loss per unit of heat load.

Usage:

python -m digimuh.analysis_03_heat_stress --db cow.db --out results/heat

References

Identifying and predicting heat stress events for grazing dairy cows using rumen temperature boluses. JDS Comm. 2024.

digimuh.analysis_03_heat_stress.load_heat_data(con)

Load heat stress view and add per-animal Z-scored temperature.

Parameters:

con – Database connection with views.

Returns:

DataFrame with per-animal Z-scored rumen temperature and a composite heat load index.

Return type:

DataFrame

digimuh.analysis_03_heat_stress.sigmoid(x, L, k, x0, b)

Four-parameter sigmoid: L / (1 + exp(-k*(x - x0))) + b.

Parameters:

• x (ndarray)

• L (float)

• k (float)

• x0 (float)

• b (float)

Return type:

ndarray

digimuh.analysis_03_heat_stress.fit_dose_response(df, min_days=30)

Fit per-animal sigmoid dose-response: rumen_temp_z = f(THI).

The inflection point (x0) represents the THI at which the cow’s temperature begins to rise sharply — its personal heat tolerance threshold.

Parameters:

• df (DataFrame) – DataFrame from load_heat_data().

• min_days (int) – Minimum number of observation days per animal.

Returns:

animal_id, thi_threshold (x0), slope (k), n_days.

Return type:

DataFrame with one row per animal, columns

digimuh.analysis_03_heat_stress.compute_production_impact(df)

Estimate milk yield loss attributable to heat load.

Bins the heat_load_index into quartiles and computes mean milk yield per bin.

Parameters:

df (DataFrame) – DataFrame from load_heat_data().

Returns:

Summary DataFrame with heat load bins and mean production.

Return type:

DataFrame

digimuh.analysis_03_heat_stress.plot_heat_overview(df, dose, out_dir)

Generate heat stress analysis plots.

Parameters:

• df (DataFrame) – Full analysis DataFrame.

• dose (DataFrame) – Per-animal dose-response fit results.

• out_dir (Path) – Output directory.

Return type:

None

digimuh.analysis_03_heat_stress.main()

Entry point for heat stress analysis.

Return type:

None

Rumen mechanical–chemical–production coupling analysis.

Uses v_analysis_digestive which joins daily smaXtec motility/pH with HerdePlus MLP milk composition.

The key insight: reticulorumen contraction patterns (motility) drive mixing, mixing drives fermentation rate, fermentation determines the volatile fatty acid profile, and VFA ratios directly shape milk fat and protein. This pipeline has a multi-day lag.

The analysis:

1. Computes time-lagged cross-correlations between daily motility metrics and the next available MLP test-day values.

2. Builds a digestive efficiency score from the motility–pH coupling: animals where motility and pH co-vary tightly have well-functioning rumens.

3. Tests whether digestive efficiency predicts milk composition at the next MLP test day.

Usage:

python -m digimuh.analysis_06_digestive --db cow.db --out results/digestive

digimuh.analysis_06_digestive.load_digestive_data(con)

Load digestive analysis view.

Parameters:

con – Database connection with views.

Returns:

DataFrame with daily motility/pH and sparse MLP composition.

Return type:

DataFrame

digimuh.analysis_06_digestive.compute_lagged_correlations(df, predictor_cols, target_cols, max_lag_days=14)

Compute cross-correlations between daily rumen metrics and MLP test-day values at various time lags.

For each animal, MLP test days are identified (non-null target), and the mean of each predictor over the preceding N days is correlated with the target value.

Parameters:

• df (DataFrame) – Full digestive DataFrame, sorted by animal_id + day.

• predictor_cols (list[str]) – Daily rumen metrics to use as predictors.

• target_cols (list[str]) – MLP test-day columns (sparse).

• max_lag_days (int) – Maximum look-back window in days.

Returns:

lag × predictor × target → correlation.

Return type:

DataFrame

digimuh.analysis_06_digestive.compute_digestive_efficiency(df, window=7)

Compute a per-animal rolling digestive efficiency score.

Efficiency is defined as the strength of coupling between motility (contraction interval) and rumen pH over a rolling window. In a well-functioning rumen, shorter contraction intervals (faster mixing) correspond to lower pH (more active fermentation) — a negative correlation.

Parameters:

• df (DataFrame) – Full digestive DataFrame.

• window (int) – Rolling window size in days.

Returns:

DataFrame with animal_id, day, digest_eff, and digest_eff_rank (percentile within herd-day).

Return type:

DataFrame

digimuh.analysis_06_digestive.plot_digestive_results(lagged, efficiency, out_dir)

Generate digestive analysis plots.

Parameters:

• lagged (DataFrame) – Lagged correlation results.

• efficiency (DataFrame) – Digestive efficiency scores.

• out_dir (Path) – Output directory.

Return type:

None

digimuh.analysis_06_digestive.main()

Entry point for digestive efficiency analysis.

Return type:

None

Circadian rhythm analysis as a general-purpose welfare biomarker.

Uses v_analysis_circadian which provides hourly aggregates of rumen temperature, activity, and rumination per animal per day, joined with disease ground truth.

Biological rationale: healthy ruminants exhibit strong ~24h rhythms in core body temperature (nadir early morning, peak late afternoon), activity (bimodal: dawn and dusk feeding bouts), and rumination (complementary to activity — peaks at rest). Circadian amplitude collapse or phase shift is a well-established early marker of sickness in human chronobiology but barely explored in cattle.

The analysis:

1. For each animal-day, fits a single-harmonic Fourier model (24h period) to the hourly profile of each signal.

2. Extracts amplitude (strength of rhythm) and acrophase (time of peak) as daily biomarkers.

3. Computes a Circadian Disruption Index (CDI) = deviation from the animal’s own healthy-period baseline.

4. Tests whether CDI elevation precedes clinical disease onset.

Usage:

python -m digimuh.analysis_11_circadian --db cow.db --out results/circadian

digimuh.analysis_11_circadian.load_circadian_data(con)

Load hourly circadian data.

Parameters:

con – Database connection with views.

Returns:

DataFrame with hourly temp/activity/rumination per animal-day.

Return type:

DataFrame

digimuh.analysis_11_circadian.fit_circadian_harmonic(hours, values)

Fit a single 24h-period harmonic to hourly data.

Model: y(t) = A * cos(2π/24 * t - φ) + M

Uses the closed-form DFT approach (no iterative fitting): compute the first Fourier coefficient at frequency 1/24h.

Parameters:

• hours (ndarray) – Array of hour-of-day values (0–23).

• values (ndarray) – Corresponding signal values.

Returns:

Dict with amplitude, acrophase_h (hour of peak, 0–24), mesor (24h mean), n_hours (data points), and relative_amplitude (amplitude / mesor).

Return type:

dict

digimuh.analysis_11_circadian.extract_circadian_features(df)

Extract circadian features for each animal-day.

For each of temperature, activity, and rumination, computes the 24h Fourier amplitude, acrophase, and relative amplitude.

Parameters:

df (DataFrame) – Hourly circadian DataFrame.

Returns:

DataFrame with one row per animal-day and columns for each signal’s circadian parameters.

Return type:

DataFrame

digimuh.analysis_11_circadian.compute_disruption_index(features, baseline_days=30)

Compute Circadian Disruption Index (CDI) per animal-day.

CDI measures how far each day’s circadian parameters deviate from the animal’s own baseline (first baseline_days of healthy-period data). A Mahalanobis-like distance across amplitude, phase, and mesor of all three signals.

Parameters:

• features (DataFrame) – DataFrame from extract_circadian_features().

• baseline_days (int) – Number of initial healthy days to define the per-animal baseline.

Returns:

DataFrame with animal_id, day, cdi.

Return type:

DataFrame

digimuh.analysis_11_circadian.plot_circadian_results(features, cdi, out_dir)

Generate circadian analysis plots.

Parameters:

• features (DataFrame) – Circadian feature DataFrame.

• cdi (DataFrame) – Circadian Disruption Index DataFrame.

• out_dir (Path) – Output directory.

Return type:

None

digimuh.analysis_11_circadian.main()

Entry point for circadian analysis.

Return type:

None

Reticulorumen contraction entropy as a novel welfare biomarker.

Uses v_analysis_motility which extracts raw motility time series (contraction interval and pulse width) from smaXtec derived data.

Biological rationale: in a healthy rumen, reticulorumen contractions are quasi-periodic with modest beat-to-beat variability (analogous to healthy heart rate variability). Pathological states — acidosis, inflammation, impaction — disrupt this regularity. Very low entropy (rigid, uncoupled contractions) and very high entropy (chaotic, disorganised contractions) both indicate dysfunction.

This is directly analogous to heart rate variability (HRV) analysis in cardiology, applied to the rumen motor complex. As far as we know, this approach has not been published.

The analysis:

1. Computes sample entropy and permutation entropy of the contraction interval (mot_period) series in sliding windows.

2. Derives daily summary statistics: mean entropy, entropy SD, and entropy trend (slope).

3. Correlates entropy features with concurrent rumen pH, rumination index, and disease status.

4. Tests whether entropy changes precede clinical diagnosis.

Usage:

python -m digimuh.analysis_12_motility_entropy --db cow.db --out results/entropy

digimuh.analysis_12_motility_entropy.sample_entropy(x, m=2, r=None)

Compute sample entropy of a time series.

Sample entropy (SampEn) quantifies the regularity of a signal. Lower values indicate more self-similarity (regularity); higher values indicate more complexity/randomness.

Parameters:

• x (ndarray) – 1-D time series (must have len > m+1).

• m (int) – Embedding dimension (template length). Default 2, standard for physiological signals.

• r (float | None) – Tolerance radius. Default is 0.2 * std(x), the standard choice from Richman & Moorman (2000).

Returns:

Sample entropy value. Returns np.nan if the series is too short or constant.

Return type:

float

References

Richman JS, Moorman JR. Am J Physiol Heart Circ Physiol. 2000;278:H2039–49.

digimuh.analysis_12_motility_entropy.permutation_entropy(x, order=3, delay=1, normalize=True)

Compute permutation entropy of a time series.

Permutation entropy (PE) captures the complexity of a signal based on the ordinal patterns of consecutive values. It is robust to noise and monotonic transformations.

Parameters:

• order (int) – Embedding order (permutation length). Default 3.

• delay (int) – Embedding delay. Default 1.

• normalize (bool) – If True, normalise by log(order!) to [0, 1].

• x (ndarray)

Returns:

Permutation entropy value.

Return type:

float

References

Bandt C, Pompe B. Phys Rev Lett. 2002;88:174102.

digimuh.analysis_12_motility_entropy.compute_daily_entropy(con, window_size=50)

Compute daily motility entropy features per animal.

For each animal-day, collects the mot_period readings, computes sample entropy and permutation entropy, and derives summary statistics.

Parameters:

• con – Database connection with views.

• window_size (int) – Minimum number of motility readings per day to compute entropy (default: 50).

Returns:

DataFrame with daily entropy features per animal.

Return type:

DataFrame

digimuh.analysis_12_motility_entropy.compute_predisease_trends(entropy_df, lookback_days=7)

Test whether entropy changes before disease onset.

For each disease event, computes the mean entropy in the lookback_days before diagnosis and compares to the animal’s healthy-period baseline.

Parameters:

• entropy_df (DataFrame) – Daily entropy features.

• lookback_days (int) – Days before diagnosis to examine.

Returns:

DataFrame with pre-disease vs. baseline entropy comparison.

Return type:

DataFrame

digimuh.analysis_12_motility_entropy.plot_entropy_results(entropy_df, trends, out_dir)

Generate entropy analysis plots.

Parameters:

• entropy_df (DataFrame) – Daily entropy features.

• trends (DataFrame) – Pre-disease trend results.

• out_dir (Path) – Output directory.

Return type:

None

digimuh.analysis_12_motility_entropy.main()

Entry point for motility entropy analysis.

Return type:

None