Getting Started with LactCurveModels

Overview

LactCurveModels is an R package for fitting nonlinear lactation curve models to daily milk yield records from dairy animals. The package fits up to 20 published models simultaneously, evaluates each model using six goodness-of-fit metrics, and produces publication-ready diagnostic figures and summary tables — all in a single function call.

Key features:

Fits 20 nonlinear lactation curve models spanning four mathematical families
Supports batch processing of multiple animals from a single CSV file
Produces 15 diagnostic figures per animal plus 3 combined cross-animal outputs
Reports R², Adjusted R², RMSE, AIC, BIC, and Durbin-Watson for every model
Flexible selection of any subset of models or animals to analyse

Installation

Install the package from CRAN:

install.packages("LactCurveModels")

Or install the development version from GitHub:

# install.packages("devtools")
devtools::install_github("venkatesanraja/LactCurveModels")

Load the package:

library(LactCurveModels)

Input Data Format

Your input CSV file must contain exactly three columns in any column order:

Column	Type	Description
`animal_id`	character	Unique identifier for each animal
`dim`	integer	Days-in-milk time index (typically 1 to 20)
`dmy`	numeric	Daily milk yield in kg/day

A typical CSV file looks like this:

animal_id,dim,dmy
Animal_A,1,9.3
Animal_A,2,9.9
Animal_A,3,10.2
...
Animal_B,1,8.1
Animal_B,2,8.7
...

Note: The dim column represents a fortnightly time index across the lactation period. Each record corresponds to one fortnightly observation of daily milk yield for a given animal.

Basic Usage

Analyse All Animals with All 20 Models

results <- run_lactation_analysis(
  input_csv_path = "path/to/your/data.csv"
)

By default, all animals in the CSV are analysed using all 20 models, and outputs are saved to a folder named lactation_results in the current working directory.

Specify a Custom Output Folder

results <- run_lactation_analysis(
  input_csv_path = "path/to/your/data.csv",
  out_dir        = "D:/my_results/lactation_analysis"
)

Selecting Specific Models

Use the selected_models argument to fit only a subset of the 20 available models. This is useful when you already know which model families are most appropriate for your data, or when you want faster processing.

results <- run_lactation_analysis(
  input_csv_path  = "path/to/your/data.csv",
  selected_models = c("Wood_1967", "Wilmink_k005", "Ali_Schaeffer",
                      "Pollott_Multiplicative")
)

All 20 Available Model Names

The full list of valid model name strings is given below, grouped by mathematical family:

Family	Model Name String	Reference
Brody / Exponential	`Brody_1923`	Brody, 1923
Brody / Exponential	`Brody_1924`	Brody, 1924
Brody / Exponential	`Parabolic_Exp_Sikka`	Sikka, 1950
Brody / Exponential	`Papajcsik_Bodero1`	Papajcsik & Bodero, 1988
Brody / Exponential	`Papajcsik_Bodero3`	Papajcsik & Bodero, 1988
Wood / Gamma-type	`Wood_1967`	Wood, 1967
Wood / Gamma-type	`Wood_Dhanoa`	Dhanoa, 1981
Wood / Gamma-type	`Cappio_Borlino`	Cappio-Borlino et al., 1995
Wood / Gamma-type	`Papajcsik_Bodero2`	Papajcsik & Bodero, 1988
Wood / Gamma-type	`Pollott_Multiplicative`	Pollott, 2000
Polynomial / Linear	`Cobby_LeDu`	Cobby & Le Du, 1978
Polynomial / Linear	`Quadratic_Dave`	Dave, 1971
Polynomial / Linear	`Mixed_Log_GS`	Guo & Swalve, 1995
Polynomial / Linear	`Khandekar`	Khandekar, 1993
Polynomial / Linear	`Inverse_Poly_Nelder`	Nelder, 1966
Wilmink / Log-exp	`Wilmink_k005`	Wilmink, 1987
Wilmink / Log-exp	`Wilmink_k_estimated`	Wilmink, 1987
Wilmink / Log-exp	`Ali_Schaeffer`	Ali & Schaeffer, 1987
Wilmink / Log-exp	`Log_Quadratic_Adediran`	Adediran et al., 2012
Wilmink / Log-exp	`Morant_Gnanasakthy`	Morant & Gnanasakthy, 1989

Selecting Specific Animals

Use the selected_animals argument to analyse only a subset of the animals present in the CSV file. Animal IDs must match the values in the animal_id column exactly.

results <- run_lactation_analysis(
  input_csv_path   = "path/to/your/data.csv",
  selected_animals = c("Animal_A", "Animal_C")
)

Combined Selection of Models and Animals

All four arguments can be used together:

results <- run_lactation_analysis(
  input_csv_path   = "path/to/your/data.csv",
  selected_models  = c("Wood_1967", "Wilmink_k005", "Ali_Schaeffer"),
  selected_animals = c("Animal_A", "Animal_B"),
  out_dir          = "D:/results"
)

Understanding the Outputs

Per-Animal Outputs

For each successfully processed animal, a sub-folder named after the animal ID is created inside out_dir. It contains the following files:

CSV Tables

File	Contents
`AnimalID_parameter_estimates.csv`	Fitted parameters and standard errors
`AnimalID_summary_metrics.csv`	R², Adj. R², AIC, BIC, RMSE, Durbin-Watson
`AnimalID_actual_and_predicted_values.csv`	Observed vs predicted milk yield at each DIM
`AnimalID_residuals.csv`	Residuals at each time point for all models

Diagnostic Figures (15 PNG files)

Figure	Description
Fig01	All 20 fitted curves overlaid on observed data
Fig02	Four-panel ranked model fits (best to worst R²)
Fig03	Best 3 vs Worst 3 models side by side
Fig04	Models grouped and coloured by mathematical family
Fig05	Residual diagnostics for top 4 models (4 sub-panels each)
Fig06	Residual bubble chart across all models and time points
Fig07	R² and Adjusted R² bar chart for all models
Fig08	RMSE bar chart for all models
Fig09	AIC and BIC bar chart for all models
Fig10	Durbin-Watson statistic bar chart for all models
Fig11	R² vs RMSE performance bubble chart coloured by model family
Fig12	Multi-metric ranking dot-plot (6 metrics simultaneously)
Fig13	Predicted vs actual scatter plots (one panel per model)
Fig14	Pearson correlation heatmap of model predictions
Fig15	Predicted peak yield and time-to-peak bar charts

Combined Cross-Animal Outputs

When more than one animal is analysed, the following combined outputs are saved directly in out_dir:

File / Figure	Contents
`COMBINED_summary_metrics_all_animals.csv`	All metrics for all animals and models
`COMBINED_parameter_estimates_all_animals.csv`	All parameter estimates across all animals
`COMBINED_best_model_per_animal.csv`	Best fitting model per animal with key metrics
`COMBINED_Fig_R2_CrossAnimal.png`	Grouped bar chart comparing R² across animals
`COMBINED_Fig_RMSE_CrossAnimal.png`	Grouped bar chart comparing RMSE across animals
`COMBINED_Fig_BestModel_Overlay.png`	Best-fit curve overlay for all animals

Accessing Results in R

The function returns a named list (one element per animal) which can be accessed directly in R:

# View goodness-of-fit metrics for Animal_A
results[["Animal_A"]]$metrics_df

# View parameter estimates for Animal_A
results[["Animal_A"]]$param_table

# Access the fitted nls object for Wood (1967) model
results[["Animal_A"]]$model_fits[["Wood_1967"]]$model

# View model coefficients
coef(results[["Animal_A"]]$model_fits[["Wood_1967"]]$model)

# View individual metrics
results[["Animal_A"]]$model_fits[["Wood_1967"]]$metrics

Using fit_lactation_models() Directly

For advanced users, the lower-level function fit_lactation_models() fits models to a single animal’s pre-formatted data frame:

# Build a data frame for one animal
animal_df <- data.frame(
  x = 1:20,
  y = c(9.3, 9.9, 10.2, 10.4, 10.3, 10.1, 9.9, 9.6, 9.3, 9.0,
        8.7, 8.4, 8.1, 7.8, 7.5, 7.2, 6.9, 6.6, 6.3, 6.0),
  z = (1:20) / 365
)

# Fit all 20 models
fits <- fit_lactation_models(animal_df)

# Fit selected models only
fits <- fit_lactation_models(
  animal_df,
  selected_models = c("Wood_1967", "Wilmink_k005", "Ali_Schaeffer")
)

# Access results
fits[["Wood_1967"]]$metrics
fits[["Wood_1967"]]$predictions
coef(fits[["Wood_1967"]]$model)

The data frame passed to fit_lactation_models() must contain:

Column	Description
`x`	Integer time index (days-in-milk, e.g. 1 to 20)
`y`	Numeric daily milk yield in kg/day
`z`	Numeric, equal to `x / 365` (required by Morant-Gnanasakthy model)

Goodness-of-Fit Metrics

The following metrics are computed for every successfully fitted model:

Metric	Description
R²	Coefficient of determination
Adjusted R²	R² penalised for the number of parameters
RMSE	Root mean square error (kg/day)
AIC	Akaike information criterion — lower is better
BIC	Bayesian information criterion — lower is better
Durbin-Watson	Test statistic for autocorrelation in residuals (ideal value = 2)

Citation

If you use LactCurveModels in your research, please cite it as:

citation("LactCurveModels")

Raja TV, Lalremruati PC, Priyadharshini P, Nidhishree NS, Gurjar D, Alex R, Vohra V (2025). LactCurveModels: Lactation Curve Model Fitting for Dairy Animals. R package version 0.1.0. https://cran.r-project.org/package=LactCurveModels

Getting Help

Use ?run_lactation_analysis and ?fit_lactation_models for function-level help
Report bugs at: https://github.com/venkatesanraja/LactCurveModels/issues