Ethics statement

The present study only used data obtained from pre-existing databases; thus, approval by an institutional animal care and use committee was not necessary.

Data source

The data were obtained from a commercial dairy farm located in northern México, which is a Holstein-intensive production system with two milkings per day, and an annual average of 1,878 milking cows, excluding dry cows. The dairy farm averaged 31.6 L/cow per day (9,638 L/cow per 305 days), had a 21-day pregnancy rate of 25% and an annual culling rate of 42%. The data were extracted from the dairy management software Dairy Comp 305 (Valley Agricultural Software, Tulare, CA, USA).

Data arrangement and binary classification design

To construct an ML binary classification (Li and Tong, Reference Li and Tong2020), the common classification of culling into voluntary or involuntary is not useful (Fetrow et al., Reference Fetrow, Nordlund and Norman2006). So, the culling classification was adapted from Fetrow (Reference Fetrow1987). The first category is forced culling (non-data-based), which represents cows that have to leave the farm immediately or soon and for which there is no reasonable solution other than culling. This category includes diseases (i.e., contagious mastitis, and lameness), trauma and death. These cows were not considered in this study.

The second category is economic culling (data-based), which is the objective of this study. This category represents those cows selected by the dairy manager at the end of each month during 2020. The applied criteria for each cow considered the features of milk yield (L/d), days in milk (d), lactation number (n), pregnancy status (Pregnant/Non-pregnant), days open (d) and days pregnant (d). These cows were labelled as (i) Cull (n = 408) and represented the first label of the binary classification.

The second label represents those cows that were not selected to be culled and were labelled as (ii) Stay (n = 20,800) these enable a ratio to be determined (Cull/Stay) 1:50. As described by Krawczyk (Reference Krawczyk2016), an imbalanced dataset represents a common challenge in ML, where the optimal scenario in a binary classification will be a ratio of 1:1, but this is uncommon in real-life scenarios. To reduce this imbalance without losing power in the Cull label, an undersampling technique was applied to the Stay label to achieve a ratio of (Cull/Stay) 1:10, by the ‘RandomUnderSampler’ procedure of the library imbalanced-learn (Lemaître et al., Reference Lemaître, Nogueira and Aridas2017), this represents a ratio closer to that which the algorithms would normally have to deal with, on a typical dairy farm (Liu et al., Reference Liu, Wu and Zhou2008). Therefore, the final database for the ML binary classification resulted in 408 cows for the label Cull and 3993 for the label Stay.

Data preprocessing and statical analysis

Descriptive statistics were performed, and a scatterplot matrix was created as described by Emerson et al. (Reference Emerson, Green, Schloerke, Crowley, Cook, Hofmann and Wickham2013). This was done to visualize the multivariate data that included categorical and quantitative features from both labels using the function ‘ggpairs’ from the package ‘GGally 2.2.1’, in R v4.4.1 statistical software (https://www.R-project.org/). To avoid highly correlated features that can affect the performance of the learning process, a correlation analysis was performed among the six features before the training and testing phase using the library ‘seaborn’ (Waskom, Reference Waskom2021). A heat-map plot and the corresponding values are shown in Figure S1, where correlations >0.9 were not found, and thus all the features were included in the learning models.

Feature scaling was considered because the features have different ranges, to which some ML algorithms, including k-nearest neighbors (k-NN), support vector machine (SVM) and multilayer perceptron (MLP) are sensitive. Two methods were tested, and the one with the best performance was selected for each algorithm (Ali and Faraj, Reference Ali and Faraj2014). The first method is standardization, which removes the mean, and scales to one unit of standard deviation, and is computed as follows:

\begin{equation*}x_{std}^{\left( i \right)} = \,\frac{{{x^{\left( i \right)}}{ -\mu_x}}}{{_x}}\end{equation*}

The second method was max–min normalization, which scales the features in the range 0–1, following the formula:

\begin{equation*}x_{norm}^{\left( i \right)} = \frac{{{x^{\left( i \right)}} - {x_{min}}}}{{{x_{max}} - {x_{min}}}}\end{equation*}

ML algorithms and implementation

In the present study, six ML classifier algorithms were compared to represent a range of algorithmic approaches, from linear algorithms to those including ensembles to neural network-based methods. Choosing the correct ones depends on several factors, including dataset size, feature types, data dimensionality, interpretability, linearity and computational constraints. It is a common practice to compare several algorithms that can theoretically perform the desired task (Shahinfar et al., Reference Shahinfar, Page, Guenther, Cabrera, Fricke and Weigel2014), and from the performance comparison select the best one. As described by Wolpert (Reference Wolpert1996): ‘No single classifier works best across all possible scenarios’, according to the ‘No Free Lunch’ theorem.

The preprocessed database (scaled for k-NN, SVM and MLP) was randomly split into training and testing sets using ‘train_test_split’ from ‘sklearn.model_selection’ (Pedregosa et al., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot and Duchesnay2011), with 70% allocated to train and 30% to test, with a fixed random seed to ensure reproducibility.

To optimize the performance, hyperparameter tuning was conducted using grid search with cross-validation. The GridSearchCV procedure was coupled with a tenfold stratified k-fold cross-validation to ensure robust performance estimates across class-balanced folds, accuracy was defined as the scoring metric, and parallel processing enabled. The process involved fitting the base model with all combinations of the parameter grid across the training set defined ahead for each classifier, yielding the best-performing model based on the highest mean cross-validated accuracy. The optimal hyperparameters and corresponding cross-validated score were retrieved, and the best estimator was used for the subsequent evaluation.

To run the supervised ML models, the Google Colaboratory (Bisong, Reference Bisong2019) was used. This platform is a hosted Jupyter Notebook, running Python 3.10.12 (Python Software Foundation). The libraries ‘pandas’, ‘matplotlib’ and ‘NumPy’, were used to handle databases and plot results. The ML classifier algorithms were obtained from the library ‘scikit-learn’ (Pedregosa et al., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot and Duchesnay2011) and a brief description and implementation of each one are described below.

Model performance evaluation

To evaluate the abilities of the ML algorithms employed on the supervised culling models task, several quantitative measures of performance were computed on the outputs of the testing phase after the models were trained (Luque et al., Reference Luque, Carrasco, Martín and Heras A2019; Chicco and Jurman, Reference Chicco and Jurman2020). The first metric is a 2 × 2 confusion matrix, which describes the outcome of the binary classification task, as follows:

\begin{equation*}{\text{Confusion Matrix}}\,{\mathbf{M}} = \,\left( {\begin{array}{*{20}{c}} {{\text{TP}}}&{{\text{FN}}} \\ {{\text{FP}}}&{{\text{TN}}} \end{array}} \right)\end{equation*}

where TP is the true positive (Cull cows correctly classified), TN is the true negative (Stay cows correctly classified), FP is the false positive (actual Stay cows wrongly classified as Cull) and FN is the false negative (actual Cull cows wrongly classified as Stay). A normalization was applied to scale the values, typically so that each row sums to 1 (or 100%), making it easier to interpret the model's performance across classes, especially when class imbalances exist.

From this confusion matrix several metrics were computed (Luque et al., Reference Luque, Carrasco, Martín and Heras A2019). Accuracy is the most commonly used among researchers, represents the fraction of correct classifications of the algorithm, is the ratio of correctly classified positive and negative, for all instances of the dataset, where 0 is the worst value and 1 the best value. It is computed as follows:

\begin{equation*}{\text{Accuracy}} = { }\frac{{{\text{TP}} + {\text{TN}}}}{{{{\text{n}}^ + } + {{\text{n}}^ - }}} = \frac{{{\text{TP}} + {\text{TN}}}}{{{\text{TP}} + {\text{TN}} + {\text{FP}} + {\text{FN}}}}\end{equation*}

Precision or positive predictive value is the ratio of true positive to all positive (true and false), so it determines the ability of a classifier to correctly predict positive instances, and ranges from a worst value of 0 to a best value of 1, calculated as follows:

\begin{equation*}{\text{Precision}} = { }\frac{{{\text{TP}}}}{{{\text{TP}} + {\text{FP}}}}\end{equation*}

Recall or sensitivity is the ratio of true positive to true positive and false negative. It is the ability of a classifier to find the positive predictions over all the results that are positive, and was estimated as follows:

\begin{equation*}{\text{Recall}} = \frac{{{\text{TP}}}}{{{\text{TP}} + {\text{FN}}}}\end{equation*}

The F₁-score is one of the most used tests to evaluate the performance of ML predictions, and is defined as the harmonic mean of precision and recall. F₁-score ranges from 0 to 1, and has the following shape:

\begin{equation*}{{\text{F}}_1}{\text{ score}} = 2 \cdot { }\frac{{{\text{precision }} \cdot {\text{ recall}}}}{{{\text{precision}} + {\text{recall}}}}\end{equation*}

The alternative measure to the F₁-score is the Matthews correlation coefficient (MCC), which has the advantage of being unaffected by imbalanced datasets (Chicco and Jurman, Reference Chicco and Jurman2020). It ranges in the interval −1 for perfect misclassification to +1 for perfect classification, and is computed as follows:

\begin{equation*}{\text{MCC}} = \frac{{{\text{TP}} \cdot {\text{TN}} - {\text{FP}} \cdot {\text{FN}}}}{{\sqrt {\left( {{\text{TP}} + {\text{FP}}} \right) \cdot \left( {{\text{TP}} + {\text{FN}}} \right) \cdot \left( {{\text{TN}} + {\text{FP}}} \right) \cdot \left( {{\text{TN}} + {\text{FN}}} \right)} }}\end{equation*}

A model inspection technique (permutation feature importance) was performed to estimate the contribution of each feature to the performance of the ML algorithms by using the procedure of the library scikit-learn permutation importance with the following hyperparameters: the number of repeats = 30, and random state = 42 (Fisher et al., Reference Fisher, Rudin and Dominici2019).

Labels description

Table 1 shows the summary and descriptive statistics of both labels (Cull and Stay) for quantitative features, while Table 2 shows the summary of the categorical features of the labels (Cull and Stay). Figure 1 shows generalized pairs plots for exploratory data analysis. In this learning approach, six features were employed per cow, three for productive characteristics and another three for reproductive performance.

Figure 1.

Generalized pairs plot for the six features. Notes: The quantitative–quantitative (milk yield, days in milk open days and days of pregnancy) pairs are shown in scatterplots and density plots. The categorical–categorical (pregnancy status and lactation) pairs are shown in mosaic plots. The quantitative-categorical pairs are shown in boxplots and barcode plots. The diagonal (from the left upper to right lower) reflects the marginal distributions of the features. Green = data for the label Cull and blue = data for the label Stay.

Table 1.

Summary and descriptive statics of both labels (Cull and Stay), for quantitative features

SD, standard deviation; CI, confidence interval at 95%.

Table 2.

Summary of the categorical features of the binary labels (Cull and Stay)

The mean of the feature milk yield was lower in the Cull label (22.0 ± 5.7 L/d) than in the Stay label (31.7 ± 6.3 L/d). This can be explained by the fact that milk yield is one of the most important metrics related to voluntary culling. For days in milk, the mean was higher in the Cull label than in the Stay label (249.7 ± 93.2 vs 219.5 ± 79.0 d) in the Stay label; this trend in the Cull label is usually accompanied by a low milk yield and a failure to conceive. The feature lactation number had a similar distribution in cows with both labels. In both labels, most of the cows were in first lactation, 50% in the Cull label and 41% in the Stay label. The Stay label shows the distribution of lactations in the herd; this pattern is commonly observed in herds with high reproductive performance (Hare et al., Reference Hare, Norman and Wright2006; Fricke et al., Reference Fricke, Wiltbank and Pursley2023). Even then, it can be supposed that in the Cull label, most of the cows would have three lactations or more, and this can be observed in involuntary culling, where the ages play a key role in the incidence of diseases (e.g., Rilanto et al., Reference Rilanto, Reimus, Orro, Emanuelson, Viltrop and Mõtus2020). It is important to consider that a first lactating cow does not guarantee a continued presence in the herd until the productive and reproductive parameters are met.

The reproductive parameters that comprise the secondary reason of importance are related to voluntary culling according to productive parameters such as milk yield (Weigel et al., Reference Weigel, Palmer and Caraviello2003). The first feature is pregnancy status (pregnant or non-pregnant), and the pregnant cows were higher in the Stay label 71% vs 41% in the Cull label. These indicate that a confirmed pregnancy does not guarantee permanence in the herd, basically due to other influencing factors. The mean of days of pregnancy were lower in the Cull label than in the Stay label (48.4 ± 61.8 vs 98.5 ± 75.6 days). This indicates that if a pregnant cow is culled, it would have an early pregnancy, which can be explained by those cows that have too many days in milk and low milk yield. Also, this tendency can be seen in the feature days open which represents the period since calving to conception, where the Cull label has a mean of 193.4 ± 100.3 vs Stay label with 115.5 ± 55.5 (days). This indicates that for some reason the cows in the Cull label experienced more problems with being pregnant, so even if a cow was pregnant, they were too advanced in their lactation curve.

Classification performance

Figure 2 shows the confusion matrix of the best cross-validation model of all the algorithms employed. In Table 3, several metrics were calculated from the confusion matrix and presented as described by Chicco and Jurman (Reference Chicco and Jurman2020). The performance metrics revealed significant variability across the six classifiers, reflecting differences in their underlying assumptions, capacity to model complex decision boundaries and robustness to data characteristics such as feature interactions, noise and class imbalance.

Figure 2.

Normalized confusion matrix of the classifier algorithms employed. Notes: Each row corresponds to an instance of the true class (Stay or Cull), while each column corresponds to an instance of the predicted class. As a result, the diagonal elements indicate the proportion of correctly identified classes. Misclassification is reflected in the off-diagonal elements, which represent instances incorrectly assigned to a different class due to confusion.

Table 3.

Results of the metrics for the evaluation of the prediction performance of the ML classifiers

^a These metrics have a scale of 0–1.

^b This metric has a scale of −1 to +1.

Several factors likely contributed to the variation in performance across classifiers. The first is model complexity and expressiveness of a low complexity (LR, GNB): linear (LR) or independence-based (GNB) models underperform due to limited capacity to capture non-linearities or feature dependencies, yielding lower MCC (0.725, 0.503) and high complexity (k-NN, SVM, RF, MLP): non-parametric (k-NN, RF) or non-linear (SVM with RF, MLP) models excel by adapting to intricate patterns, achieving MCC > 0.91. RF and MLP had high F₁ scores (0.992, 0.983) reflecting their ability to fit the data latent structure. Another factor is the decision boundary flexibility; LR linear boundary and GNB Gaussian assumption constrain adaptability, and missing edge cases (higher FN). k-NN local, RF piecewise, SVM’s kernel-based and MLP’s layered boundaries align better with the data separability, minimizing errors.

Regarding robustness of assumptions: GNB independence violation degrades performance (lowest recall = 0.459). LR assumes linearity, underfitting non-linear effects. k-NN, RF and MLP impose fewer assumptions, thriving on the data empirical structure (precision ≥ 0.992). Also, the dataset likely has strong clusterability (k-NN success), non-linear interactions (RF, MLP) and moderate noise (SVM robustness). The variation arises from each algorithm’s capacity to model the data underlying distribution and decision boundary. RF and k-NN outperform due to their flexibility and data-driven nature, while NB and LR lag due to restrictive assumptions. SVM and MLP offer a middle ground, excelling with tuning. These results suggest the culling task benefits from non-linear, high-capacity models, with RF as the top performer due to its ensemble robustness and feature interaction handling. In Figure S2, the learning curves are presented and reveal that RF, k-NN and MLP are the best fits for this dataset, with RF leading due to its high accuracy and stability, with a cross-validation accuracy of ∼0.995 and low variability, making it the most reliable choice for this dataset. The learning curve for the RF model shows that the training samples (ranging from 0 to 2,500) represent the individual data points used to train the model. As the number of training samples increases, the training score remains consistently high at around 0.995, while the cross-validation score improves and stabilizes at approximately 0.985. By the time the model reaches 2,500 samples, it demonstrates robustness to data-related issues, as the gap between the training and cross-validation scores narrows – indicating minimal overfitting. This suggests that the model generalizes well when provided with a sufficient amount of data, with performance plateauing beyond this point. LR and GNB struggle, likely due to their inability to capture the data complexity.

The feature importance results (Figures S3–S8) highlight that milk yield was the primary driver of culling decisions, followed by reproductive features like days open and days pregnant. Lactation number has minimal impact, and days in milk and pregnancy status vary in importance depending on the model. Overall, RF, k-NN and MLP provided the most reliable insights due to their high performance, while GNB struggles due to its restrictive assumptions.

To summarize, we aimed to transfer the knowledge of culling cows of an expert (farmer or manager) entity (regardless of whether the expert entity was correct or not) to an AI, to prove that specific culling criteria can be learned and, in the future, used as a tool for culling dairy cows in individual dairy farm. It is noteworthy to say that this methodology is for individual herds since each one has its culling criteria. For example, milk yield is one of the most important metrics related to culling, where the level of this production at which the cows are candidates for culling, changes between herds, usually, this is determined by economic estimations like for example income over cost. Also, these criteria are affected by and differs among or according to geographical regions, milk price, feeding costs, replacement costs, climate factors (heat stress) and herd feeding management (Buza et al., Reference Buza, Holden, White and Ishler2014).

To our knowledge, a methodology like the one here proposed has not previously been reported in the scientific literature, but other approaches using AI have been used. For example, Lopez-Suarez et al. (Reference Lopez-Suarez, Armengol, Calsamiglia, Castillejos, Iliadis, Maglogiannis and Plagianakos2018) used DT to predict the milk yield in the first lactation to identify those cows poorly productive in the early stages with an accuracy of 83%. Ouseleye et al. (Reference Ouseleye, Rodriguez and Araujo2023) used RF and they were able to predict the longevity of dairy cows using productive and genetic features with a coefficient of determination of 86%. Adamczyk et al. (Reference Adamczyk, Zaborski, Grzesiak, Makulska and Jagusiak2016) tested several ML algorithms to identify the culling reasons but with unfavourable results.

The number of features that are incorporated into the training and testing process of AI modelling will depend on the data availability. In this study, six features for each label were taken into account because these are the features that the expert entity (farmer or manager) of these farms takes into account for the decision-making in voluntary culling. There may be a case where more or fewer features need to be taken into account, for example, the presence of some disease, where animals that present this disease have more weight in the decision of culling. But a priori a common practice in ML is to make a dimensionality reduction by several algorithms, such as principal component analysis and their variations, where only the features that explain the maximum amount of variance are incorporated in the learning process, these have the advantage of increasing the efficacy and reduce the computational time (Bharadiya, Reference Bharadiya2023).

The feasibility of this methodology applied under on-farm conditions would need an ML-based software system following these conditions as described by Foidl and Felderer (Reference Foidl and Felderer2019): (i) An expert entity (human rationale, i.e., farmer or manager) would be needed to load the information in the interphase as the data source, indicating which cows leave the heard so the system knows that is for voluntary culling. Also, all the information of all the cows would be online, so a stable internet connection would also be needed. (ii) The engineering environment of the pre-training steps includes data integration, data storage, transformation and normalization. (iii) A final step where the selected algorithms are continuously learning from data, the culling criteria of the herd, in the first months would not be precise, but over time as more information is loaded the algorithms would increase their accuracy. (iv) If all the prior steps work correctly, the farmer or the technician would be able to consult in the interface the results of the trained ML algorithms, which would be cows that do not justify their place in the dairy herd based on the criteria learned from data. Even though this approach can look complicated, the more feasible way is that this type of engineering is linked with the daily management software of the farm like Dairy Comp 305 (Valley Agricultural Software, Tulare, CA, USA), AfiFarm (Afimilk, Kibbutz Afikim, Israel) or SmartDairy (Boumatic, Boumatic, Madison, WI, USA) as examples, so the farmer doesn’t have to load the information in several platforms.

In conclusion, the learning approach proposed in this research shows that data-based culling criteria can be learned and predicted by AI, with variable performance among ML algorithms. In this study, with the specific features used by this herd, RF is the optimal candidate to be deployed in on-farm conditions. To increase the reliability of these procedures, these approaches need to be tested in several farms, under different scenarios (i.e., geographic and management variability) and varieties of features.