---
title: "EMOTIONS: Ensemble Models fOr lacTatION curveS"
author: 
- name: "Pablo A. S. Fonseca"
  affiliation: 
  - "Instituto de Ganadería de Montaña (CSIC-Univ. de León), 24346 Grulleros, León"
  email: "p.fonseca@csic.es"
- name: "Marcos Prates"
  affiliation: 
  - "Department of statistics, Universidade Federal de Minas Gerai, Brasil"
- name: "Aroa Suárez-Vega"
  affiliation: "Dpto. Producción Animal, Facultad de Veterinaria, Universidad de León (24007)"
- name: "Ruth Arribas Gonzalo"
  affiliation: "Dpto. Producción Animal, Facultad de Veterinaria, Universidad de León (24007)"
- name: "Beatriz Gutierrez-Gil"
  affiliation: "Dpto. Producción Animal, Facultad de Veterinaria, Universidad de León (24007)"
- name: "Juan José Arranz"
  affiliation: "Dpto. Producción Animal, Facultad de Veterinaria, Universidad de León (24007)" 

date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{EMOTIONS}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(EMOTIONS)
```

# **Introduction**

The mathematical representation of lactation curves has significant applications in various areas of animal science. Models that simulate milk production under different conditions are valuable for physiologists, nutritionists, and geneticists, allowing them to study mammary gland function and test hypotheses. These models also support management decisions related to timing and efficiency. Over the past decades, numerous models have been proposed for representing lactation curves in dairy species. These models differ mainly in their regression type (linear or nonlinear), number of parameters, relationships among parameters, and ability to represent lactation patterns such as peak yield, time at peak, and persistency.

Despite the advantages of having a diverse range of models, selecting a single model to represent the lactation curve can present limitations. Typically, model selection is based on comparing metrics such as the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). However, this metric-based approach assumes that the chosen metric is an optimal criterion for model selection. In practice, these metrics may fail to correctly identify the best model, especially when no single model is clearly superior. Furthermore, using a single model selected based on these metrics can lead to overfitting and introduce bias, particularly when the dataset is noisy or contains numerous variables.

Ensemble modeling and model averaging are powerful techniques that enhance robustness, accuracy, and generalization in predictive modeling. Instead of relying on a single model, ensemble methods combine multiple models to reduce variance, mitigate overfitting, and improve predictive performance. Model averaging incorporates model uncertainty by weighting predictions according to their posterior probabilities, leading to more reliable estimates. These approaches are particularly valuable when individual models exhibit varying performance across different datasets or conditions. By integrating multiple models, ensemble methods improve stability and resilience, making them especially useful in complex biological and ecological systems where uncertainty quantification is crucial.

The **EMOTIONS** package provides a set of tools for fitting 47 different lactation curve models previously reported in the literature. Some of these models and the pre-defined starting parameters were obtained from the lactcurves R package (https://CRAN.R-project.org/package=lactcurves). Once the data is fitted to each model, ensemble predictions are generated using bagging based on AIC, BIC, root mean square percentage error (RMSPE), mean absolute error (MAE), and variance. Additionally, the package provides predictions for daily milk records using Bayesian Model Averaging (BMA) and calculates cosine similarity for each model's predictions. The ranking of models across individual predictions can be visualized using the **RidgeModels** and **ModelRankRange** functions, which help users better understand the weight assigned to each model. Furthermore, the **PlotWeightLac** function allows users to compare predicted and actual daily milk records. Lastly, **EMOTIONS** enables the estimation of resilience indicators based on lag-1 autocorrelation, logarithm of residual variance, and residual skewness using predicted daily milking records.

## **Installing the Package**

Use the following code to install **EMOTIONS**:

```{r, eval=FALSE}
install.packages("EMOTIONS")
```

## **Analysis**

### **Loading the Package**

```{r}
library(EMOTIONS)
```

### **Input Data**

**EMOTIONS** includes a dummy dataset containing daily milk records (up to 210 days) for 100 unique individuals. This dataset can be accessed as follows:

```{r}
# Load the dummy dataset
data("LacData")

# Display the first rows
head(LacData)
```

The dataset consists of three columns:

- **ID:** Unique individual identifier
- **DIM:** Days in milk
- **DMY:** Daily milk yield

### **Fitting Lactation Curve Models and Generating Ensembles**

The core function of **EMOTIONS** is **LacCurveFit**. This wrapper function integrates multiple supplementary functions that fit 47 lactation curve models. It takes daily milk production and days in milk from **LacData** as input. The following arguments must be provided:

- **data:** Data frame containing daily milking records
- **ID:** Column name containing unique individual IDs
- **trait:** Column name containing daily milking records
- **dim:** Column name containing days in milk records
- **alpha:** Penalization factor (ranging from 0 to 1) for model weight estimation (this parameter should be fine-tuned based on predictions)
- **models:** Vector specifying models to include. The default is **"All"**, which includes all 47 models. Alternatively, users can specify a subset of models.
- **param_list:** List specifying model parameters (optional)

Example usage:

```{r, warning=FALSE}
# Running model fitting and ensemble modeling
out.ensemble <- LacCurveFit(
  data = LacData, ID = "ID", trait = "DMY",
  dim = "DIM", alpha = 0.1,
  models = "All", param_list = NULL, silent=TRUE
)
```

The output of **LacCurveFit** consists of three main lists:

1. **converged_models**: Stores the converged models for each individual.
2. **models_weight**: Stores model weights and ranks for all ensemble methods.
3. **production**: Stores predicted daily milk records for all ensemble methods.

Checking the first six fitted models for individual **ID2**:

```{r}
head(out.ensemble$converged_models$ID2)
```

Checking the first six rows of the **models_weight** list:

```{r}
head(out.ensemble$models_weight$ID2)
```

Checking the first six rows of the **production** list:

```{r}
head(out.ensemble$production$ID2)
```

### **Available Weighting Methods**

The ensemble predictions are obtained using different weighting methods:

- **weight_AIC:** Weighted by AIC
- **weight_BIC:** Weighted by BIC
- **weight_RMSPE:** Weighted by RMSPE
- **weight_MAE:** Weighted by MAE
- **weight_BMA:** Bayesian Model Averaging
- **weight_Var:** Weighted by residual variance
- **weight_CosSquared:** Weighted by cosine similarity
- **SMA:** Simple Model Average (equal weights for all models)

### **Visualizing Model Ranking**

Model ranking across individuals can be visualized using **RidgeModels** and **ModelRankRange**.

```{r, fig.width=7, fig.height=9}
RidgeModels(out.ensemble, metric = "AIC_rank")
```

```{r, fig.width=9, fig.height=9}
ModelRankRange(out.ensemble, metric = "AIC_rank")
```

Note that the function ModelRankRange displays the number of individuals which the model was converged in front of each line.

### **Plotting Actual vs. Predicted Daily Milk Yield**

Another visualization option provided by EMOTIONS is the **PlotWeightLac** function. This function plots the actual daily milk daily production and the predicted values obtained by the ensemble model.The arguments that must be provided to this function are:
  
  - *data:* The list generated by the LacCurveFit with the predicted daily milking records
- *ID:* The ID of the individual that will have the daily milking records plotted
- *trait:* The name of the column containing actual 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

```{r, fig.width=7, fig.height=8}
PlotWeightLac(
  data = out.ensemble, ID = "ID2",
  trait = "DMY", metric = "weight_AIC",
  dim = "DIM", col = c("red", "blue")
)
```

### **Customizing Model Selection**

The function LacCurveFit allows the customization of the models to be included in the ensemble as well as the parameters of the models.

The EMOTIONS package has a dataset with the list of models, and its respective acronyms, available in the package. This dataset can be accessed using the following command.

```{r}
data("models_EMOTIONS")
head(models_EMOTIONS)
```

Users can select specific models to include in the ensemble. For example, the following models will be used to generate the ensemble: `wil`, `wilk`, `wilycsml`, `DiG`, `DiGpw`, `legpol3`, `legpol4`, `legpolWil`, `cubsplin3`, `cubsplin4`, `cubsplin5`, `cubsplindef`, `wilminkPop`, and `qntReg`.

```{r, warning=FALSE}
out.ensemble.sub <- LacCurveFit(
  data = LacData,
  ID = "ID",
  trait = "DMY",
  dim = "DIM",
  alpha = 0.1,
  models = c("wil", "wilk", "wilycsml", "DiG", "DiGpw", "legpol3", 
             "legpol4", "legpolWil", "cubsplin3", "cubsplin4", "cubsplin5", 
             "cubsplindef", "wilminkPop", "qntReg"),
  param_list = NULL
)
```

### Evaluating Model Ranks

The ranking of the selected models can be evaluated using a ridge density plot based on AIC scores:

```{r, fig.width=7, fig.height=8}
RidgeModels(out.ensemble.sub, metric = "AIC_rank")
```

### Customizing Model Parameters

Users can access the list of models that support parameter customization and their respective parameters using the `model_pars` dataset:

```{r}
data(model_pars)
head(model_pars)
```

For instance, the starting values of the parameters for models `MM` and `wil` can be modified as follows:

```{r, warning=FALSE}
edited_list <- list(
  MM = c(a = 20, b = -2),
  wil = c(a = 35, b = -5, c = -0.01, k = 0.2)
)

out.ensemble.edited <- LacCurveFit(
  data = LacData,
  ID = "ID",
  trait = "DMY",
  dim = "DIM",
  alpha = 0.1,
  models = "All",
  param_list = edited_list
)
```

Model convergence can be checked using:

```{r}
out.ensemble.edited$converged_models$ID2[["MM"]]
out.ensemble.edited$converged_models$ID2[["wil"]]
```

### Calculating Resilience Indicators

The `ResInd` function calculates resilience estimators, including logarithm of variance, lag-1 autocorrelation, and skewness, based on the weighted predictions. The required parameters are:

- `production_df`: A list containing data frames with daily production records (either actual or predicted) from `LacCurveFit`.
- `dim_filter_range`: A vector defining the lower and upper limits for filtering lactation records at the beginning and end of lactation. If no filtering is needed, set the first two values as the minimum DIM and the last two as the maximum DIM.
- `outlier_sd_threshold`: A threshold specifying the maximum standard deviations for identifying outlier resilience indicators.
- `weight`: The column name containing the selected ensemble prediction (default: `weight_AIC`).
- `trait`: The column name containing daily milk yield records.
- `DIM`: The column name containing days in milk records.
- `ID_col`: The column name containing unique animal IDs.

```{r}
out.res <- ResInd(
  out.ensemble$production,
  dim_filter_range = c(1, 7, 203, 210),
  outlier_sd_threshold = 4,
  weight = "weight_AIC",
  trait = "DMY",
  DIM = "DIM",
  ID_col = "ID"
)
```

The output of `ResInd` is a list containing:

1. `ri_filtered`: A data frame with daily milk production after filtering and the estimated resilience indicators.
2. `dev_list`: A list of deviations for each animal and DIM.
3. `removed_samples`: A list of animals identified as outliers and removed from analysis.
4. `ri_stats`: A data frame summarizing the resilience indicators.

The first six rows of the filtered dataset can be accessed as follows:

```{r}
head(out.res$ri_filtered)
```

The summary statistics for the resilience indicators can be retrieved using:

```{r}
out.res$ri_stats
```

### Imputation of missing records

The EMOTIONS package also allows the uer to impute missing MY using the ensemble created. For this, the function **imp_my** must be used. The function receive as input the list containing the data frames with the daily production records obtained from the LacCurveFit function and the a vector with the days in milk where the milk yield will be imputed. It can contain observed and missing DIM. Below we will use a example the records from `ID2`. The MY from days 32,34,36,37 and 40 will be masked and the whole 1-305 DIM MY will be predicted using the ensemble.

```{r, warning=FALSE}

#Saving the original values
original.MY<-LacData[which(LacData$ID=="ID2" & LacData$DIM
                           %in%c(32,34,36,37,40)),
                     "DMY"]

#Masking the original values
LacData.imp<-LacData

LacData.imp[which(LacData.imp$ID=="ID2" & LacData.imp$DIM
                           %in%c(32,34,36,37,40)),
                     "DMY"]<-NA

LacData.imp<-LacData.imp[!is.na(LacData.imp$DMY),]

# Running model fitting and ensemble modeling
out.ensemble.imp <- LacCurveFit(
  data = LacData.imp, ID = "ID", trait = "DMY",
  dim = "DIM", alpha = 0.1,
  models = "All", param_list = NULL, silent=TRUE
)

#Imputing missing masked MY
out.imp<-imp_my(out.ind = out.ensemble.imp, dim = 1:305)

#Checking output
out.imp[which(out.imp$ID=="ID2" & out.imp$DIM
                           %in%c(32,34,36,37,40)),]

#Comparing with the original values
original.MY
```

### Identification of milk loss events and estimation of resilience indicators

The function `milkloss_detect` allows the users to identify milk loss events based on different criteria.

- "pctbase": recovery when the observed value reaches a given fraction of the baseline (rec), for a given number of consecutive days (stick); 

- "band": recovery when the observation is inside a tolerance band around the baseline (±tol), for at least stick consecutive days; 

- "resid": recovery when the residual has improved enough from the nadir (by a fraction rec of the nadir’s absolute residual) for stick consecutive days.

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 over t in [t_start, t_end] of 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’s 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 over t from t_start to t_end of [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 'r' consecutive days (controlled by arguments such as '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 + r - 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 = minimum tau >= 0 such that 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 (volume × duration).

Conceptually:

AUC_deviation = integral of [baseline(t) - obs(t)] dt from t_start to t_end

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

The function returns a list with two data frames containing:

*episodes:* individual milk loss events and their resilience indicators;

*aggregates:* milk loss events aggregated per individual.

```{r}
#Creating a input file based on data frame with all ensemble predictions
out.imp.ensemb <- do.call(rbind, out.ensemble$production)

#Estimating the milk loss events
res.ensem <- milkloss_detect(
  data = out.imp.ensemb,
  id_col = "ID",
  dim_col = "DIM",
  MY_col = "DMY",
  MY_pred = "weight_AIC",
  dim_start = 11, dim_end = 294,
  drop_pct = 0.10, min_len = 1,
  rec_mode = "pctbase", rec = 0.99, stick = 2
)

#Checking the output
str(res.ensem)
```

The EMOTIONS package also contains a function, named `PlotMilkLoss`, that allows to plot the milk loss events identified by the `milkloss_detect` function. The arguments of `PlotMilkLoss` are:

*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

```{r, fig.width=17, fig.height=8}
PlotMilkLoss(data=out.imp.ensemb, 
                       ID="ID2", 
                       res.milkloss=res.ensem,
                       MY_col="DMY", 
                       MY_pred="weight_AIC",
                       col = c("red", "blue", "darkgreen"),
                       id_col = "ID")
```