| Type: | Package |
| Title: | Ensemble Models for Lactation Curves |
| Version: | 1.1 |
| Description: | Lactation curves describe temporal changes in milk yield and are key to breeding and managing dairy animals more efficiently. The use of ensemble modeling, which consists of combining predictions from multiple models, has the potential to yields more accurate and robust estimates of lactation patterns than relying solely on single model estimates. The package EMOTIONS fits 47 models for lactation curves and creates ensemble models using model averaging based on Akaike information criterion (AIC), Bayesian information criterion (BIC), root mean square percentage error (RMSPE) and mean squared error (MAE), variance of the predictions, cosine similarity for each model's predictions, and Bayesian Model Average (BMA). The daily production values predicted through the ensemble models can be used to estimate resilience indicators in the package. The package allows the graphical visualization of the model ranks and the predicted lactation curves. Additionally, the packages allows the user to detect milk loss events and estimate residual-based resilience indicators. |
| Depends: | R (≥ 4.2) |
| Imports: | dplyr, orthopolynom, quantreg, minpack.lm, tidyr, ggplot2, ggridges, parameters, rlang, tidyselect, splines, tibble |
| License: | GPL-3 |
| Encoding: | UTF-8 |
| LazyData: | true |
| RoxygenNote: | 7.3.3 |
| Suggests: | knitr, rmarkdown |
| VignetteBuilder: | knitr |
| NeedsCompilation: | no |
| Packaged: | 2026-01-28 08:30:57 UTC; pablofonseca |
| Author: | Pablo Fonseca [aut, cre], Marcos Prates [aut], Aroa Suarez-Vega [aut], Ruth Arribas Gonzalo [aut], Beatriz Gutierrez-Gil [aut], Juan José Arranz [aut] |
| Maintainer: | Pablo Fonseca <p.fonseca@csic.es> |
| Repository: | CRAN |
| Date/Publication: | 2026-01-28 10:30:11 UTC |
Estimate normalized model's weights based on a Expectation–Maximization (EM) algorithm using a gamma distribution
Description
Estimate normalized model's weights based on a Expectation–Maximization (EM) algorithm using a gamma distribution
Usage
BMAweight_gamma(converged_models, x, trait)
Arguments
converged_models |
A list containing the fitted models for the individual |
x |
data frame containing the daily milking records for the individual |
trait |
The name of the column containing daily milking records |
Value
A vector containing the model normalized weight based on the Bayeasian Model Average (BMA) appaoch proposed by Duan et al. (2006) with an adaptation for to use a gamma distribution
Estimate normalized model‘s weights based on the cosine similarity for each model’s predictions
Description
Estimate normalized model‘s weights based on the cosine similarity for each model’s predictions
Usage
CosSquaredWeight(converged_models, x)
Arguments
converged_models |
A list containing the fitted models for the individual |
x |
data frame containing the daily milking records for the individual |
Value
A vector containing the model normalized weight
Estimate the Akaike information criterion (AIC), Bayeasian information criterion (BIC), root mean square percentage error (RMSPE) and mean squared error (MAE) for the fitted models
Description
Estimate the Akaike information criterion (AIC), Bayeasian information criterion (BIC), root mean square percentage error (RMSPE) and mean squared error (MAE) for the fitted models
Usage
GetLacModelsMetrics(converged_models, x, trait)
Arguments
converged_models |
A list containing the fitted models for the individual |
x |
data frame containing the daily milking records for the individual |
trait |
The name of the column containing daily milking records |
Value
A data frame containing the AIC, BIC, RMSPE and MAE for each fitted model for the individual
A wrap function to the ModelsLac function that allows the fit of lactation curve models based on daily production and days in milk records simultaneously for a list of animals
Description
The function uses a data frame containing the daily milking records as input
Usage
LacCurveFit(
data,
ID,
trait,
dim,
alpha = 0.5,
models = "All",
param_list = NULL,
silent = T
)
Arguments
data |
A data frame containing the daily milking records |
ID |
The name of the column containing the unique IDs of the individuals |
trait |
The name of the column containing daily milking records |
dim |
The name of the column containing days in milk records |
alpha |
A penalization factor, ranging from 0 to 1, for the estimation of the model's weight |
models |
A vector describing the models to be included in the analysis. In total, 47 models are included in EMOTIONS. The default option is "All", which results in the inclusion of the 47 models. Alternatively, a vector containing any subset of the following models can be provided: "MMR","MME","brody23","brody24", "SCH","SCHL","PBE","wood","DHA", "CB","QP","CLD","PapBo1","PapBo2", "PapBo3", "PapBo4", "PapBo6", "GS1", "GS2","LQ", "wil", "wilk", "wilycsml", "BC", "DJK","MG2", "MG4", "MG", "KHN", "AS", "FRP","PTmult","PTmod", "MonoG", "MonoGpw", "DiG", "DiGpw","legpol3", "legpol4", "legpolWil", "cubsplin3", "cubsplin4", "cubsplin5", "cubsplindef", "wilminkPop", "qntReg". |
param_list |
A list composed by the models, named as in the models parameter, and the repective parameters included in the models. |
silent |
A logical string defining if warning should be printed or not during the model fitting. The default value is TRUE (not printing warnings). |
Value
A list containing the fitted models, the model's weigths and ranks for each weighting strategy, and the predicted daily production obtained through the model ensemble for each weighting strategy
A data frame containing the daily milk yield for 100 individuals up to 210 days in milk
Description
A data frame containing the daily milk yield for 100 individuals up to 210 days in milk
Usage
data(LacData)
Format
A data frame with 21821 rows and 3 variables:
ID: Unique individual code
DIM: Days in milk
DMY: Daily milk yield
Examples
data(LacData)
Create a line plot that shows the range of the ranks obtained for each model across the individuals
Description
Create a line plot that shows the range of the ranks obtained for each model across the individuals
Usage
ModelRankRange(LacCurveFit, metric = "AIC_rank")
Arguments
LacCurveFit |
The object obtained from the LacCurveFit function |
metric |
The name of the metric to be use to plot the model´s ranks |
Value
A line plot that shows the range of the ranks obtained for each model across the individuals
Performs the model fitting and the weight assignment based on different strategies for each individual ID
Description
Performs the model fitting and the weight assignment based on different strategies for each individual ID
Usage
ModelsLac(
data,
ID_col,
ID,
trait,
dim,
alpha,
models,
param_list = NULL,
silent = silent
)
Arguments
data |
A data frame containing the daily milking records |
ID_col |
The name of the column containing the unique IDs of the individuals |
ID |
The individual ID that is being analyzed |
trait |
The name of the column containing daily milking records |
dim |
The name of the column containing days in milk records |
alpha |
A penalization factor, ranging from 0 to 1, for the estimation of the model's weight |
models |
A vector describing the models to be included in the analysis. In total, 47 models are included in EMOTIONS. The default option is "All", which results in the inclusion of the 47 models. Alternatively, a vector containing any subset of the following models can be provided: "MMR","MME","brody23","brody24", "SCH","SCHL","PBE","wood","DHA", "CB","QP","CLD","PapBo1","PapBo2", "PapBo3", "PapBo4", "PapBo6", "GS1", "GS2","LQ", "wil", "wilk", "wilycsml", "BC", "DJK","MG2", "MG4", "MG", "KHN", "AS", "FRP","PTmult","PTmod", "MonoG", "MonoGpw", "DiG", "DiGpw","legpol3", "legpol4", "legpolWil", "cubsplin3", "cubsplin4", "cubsplin5", "cubsplindef", "wilminkPop", "qntReg" |
param_list |
A list composed by the models, named as in the models parameter, and the repective parameters included in the models. |
silent |
A logical string defining if warning should be printed or not during the model fitting. The defaul is TRUE (not printing warnings). |
Value
A list containing the fitted models, the model's weigths and ranks, and the predicted daily production obtained through the model ensemble
Define the parameters for the lactation curve models to be fitted
Description
Define the parameters for the lactation curve models to be fitted
Usage
ParDef(x, trait, dim)
Arguments
x |
data frame containing the daily milking records for the individual |
trait |
The name of the column containing daily milking records |
dim |
The name of the column containing days in milk records |
Value
A list containing the parameters to be included in the lactation curve models
Plot the actual daily milk daily production and the predicted values highlighting the detected milk loss events
Description
Plot the actual daily milk daily production and the predicted values highlighting the detected milk loss events
Usage
PlotMilkLoss(
data,
ID,
res.milkloss,
MY_col,
MY_pred,
col = c("red", "blue", "darkgreen"),
id_col = "ID"
)
Arguments
data |
A data frame containing the observed and predicted daily milking records |
ID |
The ID of the individual that will have the daily milking records plotted |
res.milkloss |
The object with the output of milkloss_detect function |
MY_col |
The name of the column containing the observed milk yield |
MY_pred |
The name of the column containing the predicted milk yield |
col |
The colors of the actual, predicted values, and milk loss events. In this order |
id_col |
The name of the column containing the individual IDs |
Value
A plot with the actual daily milk daily production and the predicted values highlighting the detected milk loss events
Plot the actual daily milk daily production and the predicted values obtained by the ensemble model
Description
Plot the actual daily milk daily production and the predicted values obtained by the ensemble model
Usage
PlotWeightLac(
data,
ID,
trait,
metric,
dim,
col = c("red", "blue"),
point_size = 2,
line_size = 1,
axis_text_size = 15,
axis_title_size = 15
)
Arguments
data |
The object generated by the LacCurveFit. |
ID |
The ID of the individual that will have the daily milking records plotted |
trait |
The name of the column containing daily milking records. |
metric |
The name of the strategy used obtained the predicted values through the ensemble model. |
dim |
The name of the column containing days in milk records |
col |
The colors of the actual and predicted values. |
point_size |
Numeric value indicating the size of the observed data points in the plot. |
line_size |
Numeric value indicating the thickness of the regression line. |
axis_text_size |
Numeric value defining the font size of the axis tick labels. |
axis_title_size |
Numeric value defining the font size of the axis titles. |
Value
A plot with the actual and predicted daily milk production across the days in milk
A function to estimate resilience estimators (logarithm of variance, lag1 autocorrelation and skewness) based on daily milk production records
Description
A function to estimate resilience estimators (logarithm of variance, lag1 autocorrelation and skewness) based on daily milk production records
Usage
ResInd(
production_df,
dim_filter_range = c(1, 7, 203, 210),
outlier_sd_threshold = 4,
weight = "weight_AIC",
trait,
DIM,
ID_col
)
Arguments
production_df |
The list containing the data frames with the daily production records (actual or predicted) obtained from the LacCurveFit function |
dim_filter_range |
A vector containing the lower and upper limits to remove lactation records from the begin and end of the lactation, if needed. If it is not necessary to remove daily records, the first two values can be set as the minimum days in milk value and the last two as the maximum days in milk values |
outlier_sd_threshold |
A threshold defining the maximum standard deviations to consider an individual resilience indicator value |
weight |
The name of the column containing the selected ensemble prediction. The default is weight_AIC |
trait |
The name of the column containing daily milking records |
DIM |
The name of the column containing days in milk records |
ID_col |
The name of the column containing the unique IDs of the individuals |
Value
A list containing the daily milk production values after filtering, the list of removed animals, and a data frame with the resilience indicators
The function RidgeModels allows the visualization of the distribution of model's ranks across individuals using ridge density plots
Description
The function RidgeModels allows the visualization of the distribution of model's ranks across individuals using ridge density plots
Usage
RidgeModels(LacCurveFit, metric = "AIC_rank")
Arguments
LacCurveFit |
The object obtained from the LacCurveFit function |
metric |
The name of the metric to be use to plot the model´s ranks |
Value
A ridge density plots for the models included in the ensemble
Estimate normalized model's weights based on the variance of the predictions
Description
Estimate normalized model's weights based on the variance of the predictions
Usage
VarWeight(converged_models, x)
Arguments
converged_models |
A list containing the fitted models for the individual |
x |
data frame containing the daily milking records for the individual |
Value
A vector containing the model normalized weight
Impute missing daily milk yields using the ensemble created
Description
Impute missing daily milk yields using the ensemble created
Usage
imp_my(out.ind = NULL, dim = NULL)
Arguments
out.ind |
The list containing the data frames with the daily production records obtained from the LacCurveFit function |
dim |
A vector with the days in milk where the milk yield will be imputed. It can contain observed and missing DIM |
Value
A data frame containing the imputed milk yields for the days in milk informed.
Identify milk loss events and resilience indicators from daily milk yields
Description
Identify milk loss events and resilience indicators from daily milk yields
Usage
milkloss_detect(
data,
id_col,
dim_col,
MY_col = "MY_real",
MY_pred,
dim_start = 1L,
dim_end = 305L,
rec_mode = c("pctbase", "band", "resid"),
drop_pct = 0.1,
min_len = 1L,
tol = 0.05,
stick = 3L,
rec = 1
)
Arguments
data |
A data frame containing the observed and predicted daily milking records. |
id_col |
The name of the column containing the individual IDs. |
dim_col |
The name of the column containing the days in milk. |
MY_col |
The name of the column containing the observed milk yield. |
MY_pred |
The name of the column containing the predicted milk yield (baseline). |
dim_start |
The first day in milk to consider when identifying milk loss events and resilience indicators. |
dim_end |
The last day in milk to consider when identifying milk loss events and resilience indicators. |
rec_mode |
How "recovery" is defined. One of:
"pctbase": recovery when the observed value reaches a given fraction of the baseline ( |
drop_pct |
Minimum relative drop from the anchor (baseline reference) to accept an episode. |
min_len |
Minimum number of consecutive days with negative residuals required to define an episode. |
tol |
Used when the "band" mode is selected. Half-width of the tolerance band around baseline in relative terms. |
stick |
Minimum number of consecutive days in recovery to consider an episode finished. |
rec |
Minimum relative recovery from the nadir to finish an episode (used in "pctbase" and "resid" modes). |
Details
The function computes several descriptors of milk-yield perturbation episodes.
1) Nadir (day of minimum)
The worst day inside the episode (deepest point of the perturbation).
t_hat = argmin_{t in [t_start, t_end]} obs(t)
Nadir = obs(t_hat)
where t_start and t_end are the episode boundaries.
2) Amplitude (drop)
Depth of the dip relative to the baseline at the episode start.
A = baseline(t_start) - obs(t_hat)
Some variants use baseline(t_hat) instead of baseline(t_start); here the
start of the episode is used as the reference.
3) ML_per_event (AUD)
Total milk lost (in baseline units) over the episode, i.e., the integrated milk deficit.
ML_per_event = AUD = sum_{t=t_start..t_end} [baseline(t) - obs(t)]
In discrete data, AUD is computed with day-weighting: each observation contributes
(baseline(t) - obs(t)) * delta_days
where delta_days is the gap to the next observed DIM (last day weight = 1).
4) Time-to-baseline (TTB)
Time after the nadir until the profile returns to (and stays near) the baseline.
Recovery is declared when obs(t) re-enters a tolerance band around the
baseline and stays there for stick consecutive days (controlled by tol and stick).
We find the smallest tau >= 0 such that for all u in the interval from
t_hat + tau to t_hat + tau + stick - 1:
abs(obs(u) - baseline(u)) <= tol * baseline(u)
Then: TTB = tau.
If this condition is never satisfied before DIM 305, TTB is set to NA
(right-censored).
5) Recovery half-life (t_1_2)
Earliest time after nadir when half of the drop has been recovered.
With amplitude A as above, define the half-recovery level:
L_half = baseline(t_start) - A / 2
Then:
t_1_2 = min{tau >= 0 : obs(t_hat + tau) >= L_half}
6) Slopes (decline and recovery)
Average daily change during the decline into the nadir and during early recovery, summarizing the episode shape.
For a K-day local window:
DeclineSlope = (obs(min(t_hat, t_start + K)) - obs(t_start)) / (min(t_hat, t_start + K) - t_start)
RecoverySlope = (obs(min(t_end, t_hat + K)) - obs(t_hat)) / (min(t_end, t_hat + K) - t_hat)
7) AUC_deviation
Trapezoidal area under the curve of the milk deficit baseline(t) - obs(t)
across the whole episode. It summarizes how much milk was lost and for how long.
Conceptually:
AUC_deviation = integral_{t_start..t_end} [baseline(t) - obs(t)] dt
In practice this is approximated via the trapezoidal rule on discrete DIMs.
8) prod_decline_slope_amp
Product of the decline slope (anchor -> nadir) and the amplitude (anchor - nadir). It combines speed and depth of the decline into a single indicator of how "aggressive" the drop is.
prod_decline_slope_amp = DeclineSlope * A
9) prod_recovery_slope_TTB
Product of the recovery slope (nadir -> recovery) and time-to-baseline (TTB). It combines how fast the animal recovers with how long recovery takes, summarizing recovery efficiency.
prod_recovery_slope_TTB = RecoverySlope * TTB
Value
A list with two data frames:
-
episodes: individual milk loss events and their resilience indicators; -
aggregates: milk loss events aggregated per individual.
The resilience indicators identified are described in the Details section.
A data frame containing the models included in the EMOTIONS package that can have the parameters edited
Description
A data frame containing the models included in the EMOTIONS package that can have the parameters edited
Usage
data(model_pars)
Format
A data frame with 42 rows and 2 variables:
Model_acronym: The acronym of the models that can have the parameters edited
Parameters: The list of parameters included in each model
Examples
data(model_pars)
A data frame containing the daily milk yield for 100 individuals up to 210 days in milk
Description
A data frame containing the daily milk yield for 100 individuals up to 210 days in milk
Usage
data(models_EMOTIONS)
Format
A data frame with 47 rows and 5 variables:
Model: Model description
Model_acronym: Model acronym used in the EMOTIONS functions
Authors: List of authors of each model available in EMOTIONS
Year: Year of publication of each model available in EMOTIONS
Reference: Link of the original publication of each model available in EMOTIONS
Examples
data(models_EMOTIONS)