lactationcurve

  1from importlib.metadata import PackageNotFoundError, version
  2
  3try:
  4    __version__ = version("lactationcurve")
  5except PackageNotFoundError:
  6    __version__ = "unknown"
  7
  8__doc__ = f"""
  9A package for fitting **dairy animal lactation curves**, evaluating
 10**lactation curve characteristics (LCCs)** (time to peak, peak yield,
 11cumulative yield, persistency), and computing **305-day milk yield**
 12using the **ICAR guideline**.
 13
 14> **Contact:** Meike van Leerdam, mbv32@cornell.edu
 15>
 16> **Authors:** Meike van Leerdam, Douwe de Kok, Judith Osei-Tete, Lucia Trapanese
 17
 18> **Initial authored:** 2025‑08‑12
 19
 20> **Updated:** 2026-10-15
 21
 22---
 23
 24  The 305 day yield for milk, fat, and protein is a widely used metric
 25  in dairy production, and the International Committee for Animal
 26  Recording (ICAR) provides guidelines outlining approved methods for
 27  its calculation. However, a global survey of milk recording
 28  organizations revealed substantial variation in how these methods are
 29  implemented. The Test Interval Method is used by 74% of the
 30  organizations, reflecting a preference for methodological simplicity,
 31  but it comes with trade-offs in estimation accuracy. The use of the
 32  other approved methods showed wide variation in correction factors,
 33  standard lactation curves, test day definitions, minimum sample
 34  requirements, and exclusion criteria. Such inconsistencies can
 35  introduce yield variability that complicates comparisons, for
 36  example in international breeding value evaluation, and limit the
 37  metric’s usefulness in universal models, such as decision support
 38  tools. Thus, the objective of this work was to reformulate the ICAR
 39  guideline section 2, procedure 2, into a unified, transparent, and
 40  accessible software implementation to improve standardization,
 41  enhance documentation, support continuous development, and increase
 42  the accuracy of 305 day yield estimation.
 43
 44  To achieve this, the ICAR guideline was translated into an
 45  open-source Python package that serves as a reference implementation
 46  for 305-day yield calculation, with lactation curve modelling at its
 47  core. In addition to the methods described in the original guideline,
 48  the package incorporates 14 lactation curve models, including both
 49  traditional and Bayesian fitting approaches as well as AI-based
 50  models. It also provides tools to derive biologically relevant
 51  characteristics such as time to peak, peak yield, cumulative yield,
 52  and persistency.
 53
 54  The framework can be directly integrated into analytical workflows,
 55  allowing users to calculate 305-day yields, fit and compare
 56  lactation curves, and derive key lactation characteristics through a
 57  single function call. These functionalities are further supported
 58  through an interactive website and an openly available [GitHub
 59  repository](https://github.com/Bovi-analytics/bovi).
 60  In addition, the project includes an online validation
 61  platform that enables users to use standardized lactation data to
 62  compare self-estimated 305-day yields against both reference
 63  calculations (from the package) and observed daily milk yields.
 64
 65  We encourage everyone to use, test, and contribute to the package,
 66  which is available under the AGPL-3.0-or-later license. We welcome feedback and
 67  suggestions for improvement, and we are committed to maintaining and
 68  updating the package to ensure it remains a valuable resource
 69  for the dairy industry and research community.
 70  For bug reports or feature requests, please submit them on the
 71  [GitHub issues page](https://github.com/Bovi-analytics/bovi/issues)
 72  or contact us directly via [email](mbv32@cornell.edu).
 73
 74
 75## Main Lactation curve models implemented:
 76
 77**Fischer**:
 78A combination of an exponential and a linear component, characterized by an
 79early exponential phase followed by an approximately linear decline.
 80
 81**Wood**:
 82An incomplete gamma-type function combining a power-law growth term t^b with
 83an exponential decay term e^(-ct).
 84
 85**Ali & Schaeffer**:
 86A linear regression model based on quadratic polynomials in standardized time
 87(t/305) and in the log-transformed inverse time log⁡(305/t), where the log
 88term increases flexibility in modelling the ascending phase and peak of
 89lactation.
 90
 91**Wilmink**:
 92A combination of an exponential and a linear component. In contrast to the
 93Fischer model, the additional parameter scaling the exponential term increases
 94flexibility in describing early-lactation dynamics and peak formation.
 95
 96**MilkBot**:
 97An empirical, mechanistically motivated four-parameter model describing
 98lactation as the development and decay of udder capacity. It consists of a
 99ramp-up phase and exponential decline, with parameters controlling scale,
100onset, growth rate, and decay.
101Both Bayesian and frequentist fitting approaches are implemented for the
102MilkBot model, allowing users to choose between traditional optimization
103methods and a probabilistic framework that incorporates prior knowledge.
104
105Additional models available for a.o. symbolic lactation curve
106characteristics (LCC) derivations:
107**Brody**, **Sikka**, **Nelder**, **Dhanoa**, **Emmans**,
108**Hayashi**, **Rook**, **Dijkstra**, **Prasad**.
109
110---
111
112## Model Formulas
113
114* **Wood** : `y(t) = a * t^b * exp(-c * t)`
115* **Wilmink** : `y(t) = a + b * t + c * exp(k * t)` with default `k = -0.05`
116* **Ali & Schaeffer** :  `t_scaled = t / 305`, `L = ln(305 / t)`
117  `y(t) = a + b*t_scaled + c*t_scaled^2 + d*L + k*L^2`
118* **Fischer** : `y(t) = a - b*t - a*exp(-c*t)`
119* **MilkBot** : `y(t) = a * (1 - exp((c - t)/b) / 2) * exp(-d*t)`
120
121---
122
123## Features
124
125- **Frequentist fitting** (numeric optimization):
126  - Wood, Wilmink, Ali & Schaeffer, Fischer, MilkBot
127- **Frequentist fitting** (algebraic least squares):
128  - MilkBot
129- **Bayesian fitting via MilkBot API**:
130  - MilkBot
131- **Lactation Curve Characteristics** — symbolic + numeric:
132  - time_to_peak, peak_yield, cumulative_milk_yield, persistency
133- **ICAR procedures cumulative milk yield:**
134  - Test Interval Method (also called the Fleischmann method)
135  - Centering Date Method
136  - Interpolation Standard Lactation Curve (ISLC) Method
137  - Best Predict Method
138- Input validation/normalization via `validate_and_prepare_inputs`
139
140---
141
142The package is organized into three main modules:
143
1441. `lactationcurve.fitting`
1452. `lactationcurve.characteristics`
1463. `lactationcurve.preprocessing`
147
148---
149
150## Output Types Summary Of Most Important Functions
151
152| Function | Output |
153
154|---------|--------|
155
156| `fit_lactation_curve` | Predicted yields (np.ndarray) |
157
158| `get_lc_parameters` | Tuple of numerical parameters |
159
160| `bayesian_fit_milkbot_single_lactation` | Dict of MilkBot parameters |
161
162| `lactation_curve_characteristic_function` | (expr, params, func) |
163
164| `calculate_characteristic` | float (LCC value) |
165
166| `test_interval_method` | DataFrame with lactation cumulative milk yield per TestId |
167
168| `centering_date_method` | DataFrame with lactation cumulative milk yield totals per TestId |
169
170| `ISLC_method` | DataFrame with lactation cumulative milk yield totals per TestId  |
171
172| `best_predict_method` | DataFrame with lactation cumulative milk yield totals per TestId |
173
174---
175
176## The meaning of a TestId
177
178
179The `TestId` is an identifier for a lactation,
180which can be used to group records belonging to the same lactation together.
181It is not the same as a cow ID, as a cow can have multiple lactations
182(e.g., across different calvings).
183If a `TestId` column is not provided,
184the package will assume all records belong to a single lactation
185and will create a `TestId` column with all values set to 0.
186
187---
188
189## Minimum number of records per lactation needed per method
190
191
192There are many different testing schemes and data collection protocols in the dairy industry,
193which can result in varying numbers of milk yield records per lactation. Furthermore,
194somethimes records are missing or excluded due to quality control,
195which can lead to incomplete lactation data.
196Or the other option is that the lactation is still ongoing and not all records are available yet.
197The minimum number of records needed to technically be able to calculate the 305-day yield
198depends on the method used:
199- Test Interval Method: At least 2 records
200- Centering Date Method: At least 2 records
201- ISLC Method: The method can be applied to lactations without any measurements,
202    in which case the result will be the population mean from the standard curve.
203- Best Predict Method: The method can be applied to lactations without any measurements,
204    in which case the result will be the population mean from the standard curve.
205- Lactation curve models:
206In general you need at least one record more than the number of parameters in the model
207 to be able to fit the model. So Wood and Fischer (3 parameters) need at least 4 records,
208 Wilmink and MilkBot (4 parameters) need at least 5 records, Ali & Schaeffer (5 parameters)
209 needs at least 6 records.
210
211
212However, in practice it is recommended to have more records than the minimum
213to ensure a reliable fit and avoid overfitting.
214Ideally, the records should be well-distributed across the lactation
215period to capture the curve's shape accurately.
216
217---
218
219## Often used abreviations
220
221- **DIM**: Days in Milk (the days since calving in the current lactation)
222- **LC**: Lactation Curve
223- **LCC**: Lactation Curve Characteristic
224- **ISLC**: Interpolation using the Standard Lactation Curve
225- **ICAR**: International Committee for Animal Recording
226- **API**: Application Programming Interface
227
228
229---
230
231## Bayesian Fitting (MilkBot API)
232
233* Set `fitting="bayesian"` and `model="milkbot"` in
234  `fit_lactation_curve` or `calculate_characteristic`.
235* Provide an **API key** via .env
236* Choose priors via custom_priors:
237    - "[CHEN](https://github.com/Bovi-analytics/Chen-et-al-2023b)" → Chen et al. 2023
238      published priors
239    - dict    → Custom priors in MilkBot format (overrides `continent`)
240      Custom priors have a specific format, to help you build them,
241      use the `build_prior` helper function. To make your own prior you need
242      for each MilkBot parameter
243      (scale,ramp,decay,onset) to specify a mean and a standard deviation (std).
244      Also provide a standard deviation for milk yield through seMilk to
245      specify the expected noise in the data.
246      This seMilk is default set to 4 kg to reflect typical day-to-day variation in milk yield,
247      but you can adjust it based on the expected variability in your data.
248
249* if no custom priors are provided, the default MilkBot priors will be used:
250    In that case set cutsom_priors to None and specify the desired continent and breed.
251    continent options:
252    - "USA"   → MilkBot USA priors (default)
253    - "EU"    → MilkBot EU priors > mainly estimates lower milk production
254
255    breed options:
256    - H → Holstein (default)
257    - J → Jersey
258
259    If also parity is provided, the continent-specific priors will be
260    further refined by parity-specific priors.
261    The priors are sensitive to the used metric. The default is kg, but if you use lb,
262    specify milk_unit="lb" to use the appropriate priors.
263
264* The helper `bayesian_fit_milkbot_single_lactation(...)`
265  normalizes differing API responses.
266* The key can be requested by sending an email to Jim Ehrlich
267  [jehrlich@MilkBot.com](mailto:jehrlich@MilkBot.com).
268* More information about the API can be found in the
269  [API documentation](https://api.milkbot.com/), or in the
270  corresponding
271  [paper](https://peerj.com/articles/54/#MainContent).
272
273---
274
275## Tutorials
276
277Hands-on programming example notebooks are available through
278*[GitHub tutorials and example notebooks](https://github.com/Bovi-analytics/bovi/tree/main/packages/models/lactationcurve/notebooks)*.
279The tutorials will guide you through how to work with the package,
280including how to import the package, fit lactation curves,
281calculate 305-day yields, and derive lactation curve characteristics.
282
283---
284## Citing the lactationcurve package
285
286If you use the `lactationcurve` package in your research, please consider citing it as follows:
287
288*van Leerdam, M. B., de Kok, D., Osei-Tete, J. A., &
289Hostens, M. (2026). Bovi-analytics/bovi:
290v.1.1.6. (v.1.1.6). Zenodo.
291https://doi.org/10.5281/zenodo.18715145*
292
293
294If you also use the Bayesian fitting functionality that relies
295on the MilkBot API, please also cite the following paper:
296
297*Ehrlich, J.L., 2013. Quantifying inter-group variability
298in lactation curve shape and magnitude with the MilkBot
299lactation model. PeerJ 1, e54.
300https://doi.org/10.7717/peerj.54*
301
302If you use the 305-day yield calculation methods based on the ICAR guideline,
303please also cite the following paper:
304Best Predict method:
305*VanRaden, P. M. (1997). Lactation yields and accuracies computed from test
306day yields and (co) variances by best prediction.
307Journal of dairy science, 80(11), 3015-3022.*
308
309ISLC:
310*Wilmink, J. B. M. (1987).
311Comparison of different methods of predicting 305-day milk yield using means
312calculated from within-herd lactation curves. Livestock Production Science, 17, 1-17.*
313
314Test Interval Method:
315*Sargent, F. D., V. H. Lyton, and 0. G. Wall, J r . 1968.
316Test interval method of calculating Dairy Herd Improvement Association records.
317Journal of dairy science, 51-170.*
318
319Centering Date Method:
320*U.S. Department of Agriculture, Bureau of Dairying. (1925).
321Cow Testing Association Letter No. 1. Washington, DC.
322
323
324---
325
326## License
327
328[GNU AGPLv3 License](https://github.com/Bovi-analytics/bovi/blob/main/LICENSE)
329
330---
331
332## Current version
333
334{__version__}
335
336"""
337
338
339# import submodules to make them available at the package level
340
341from . import characteristics, fitting, preprocessing
342
343__all__ = ["fitting", "characteristics", "preprocessing"]
344# from .characteristics import (
345#     calculate_characteristic,
346#     lactation_curve_characteristic_function,
347#     numeric_cumulative_yield,
348#     numeric_peak_yield,
349#     numeric_time_to_peak,
350#     persistency_fitted_curve,
351#     persistency_milkbot,
352#     persistency_wood,
353#     test_interval_method,
354# )
355# from .fitting import (
356#     ali_schaeffer_model,
357#     bayesian_fit_milkbot_single_lactation,
358#     brody_model,
359#     dhanoa_model,
360#     dijkstra_model,
361#     emmans_model,
362#     fischer_model,
363#     fit_lactation_curve,
364#     get_chen_priors,
365#     get_lc_parameters,
366#     get_lc_parameters_least_squares,
367#     hayashi_model,
368#     milkbot_model,
369#     nelder_model,
370#     prasad_model,
371#     rook_model,
372#     sikka_model,
373#     wilmink_model,
374#     wood_model,
375#     build_prior,
376# )
377# from .preprocessing import (
378#     PreparedInputs,
379#     standardize_lactation_columns,
380#     validate_and_prepare_inputs,
381# )
382
383# __all__ = [
384#     # Preprocessing
385#     "PreparedInputs",
386#     "standardize_lactation_columns",
387#     "validate_and_prepare_inputs",
388#     # Fitting
389#     "ali_schaeffer_model",
390#     "bayesian_fit_milkbot_single_lactation",
391#     "brody_model",
392#     "dhanoa_model",
393#     "dijkstra_model",
394#     "emmans_model",
395#     "fischer_model",
396#     "fit_lactation_curve",
397#     "get_chen_priors",
398#     "get_lc_parameters",
399#     "get_lc_parameters_least_squares",
400#     "hayashi_model",
401#     "milkbot_model",
402#     "nelder_model",
403#     "prasad_model",
404#     "rook_model",
405#     "sikka_model",
406#     "wilmink_model",
407#     "wood_model",
408#     # Characteristics
409#     "calculate_characteristic",
410#     "lactation_curve_characteristic_function",
411#     "numeric_cumulative_yield",
412#     "numeric_peak_yield",
413#     "numeric_time_to_peak",
414#     "persistency_fitted_curve",
415#     "persistency_milkbot",
416#     "persistency_wood",
417#     "test_interval_method",
418#     "centering_date_method",
419#     "build_prior",
420# ]