Identifying and predicting heat stress events for grazing dairy cows using rumen temperature boluses

Graphical Abstract Summary We present a new approach to developing predictive models of heat stress risk for grazing dairy herds. This proof-of-concept study demonstrates the potential of using a machine learning algorithm to predict heat stress events defined using rumen temperature. The approach has been enabled by animal sensors collecting large amounts of data on individual cows. The success of this approach is encouraging, as the scale and variable nature of farming outdoors in temperate climates has, to date, limited the development of our understanding and management of thermal stress among dairy cattle in these systems. To develop this approach further will require testing with independent datasets from a wider range of environments alongside model tuning.


JDS Communications® TBC; TBC
Abstract: Heat stress events in dairy cows are associated with behavioral and physiological changes such as seeking shade, increased respiration rate and body temperature, reduced milk production and psychological distress.Knowledge of the relationship between weather and animal responses to heat stress enables automated alerts using forecast weather, aiding early provision of shade or other mitigation practices.While numerous heat stress indices for cattle have been developed, these have limitations for cows exposed to wind and solar radiation, i.e., predominately grazing outdoors or managed on pasture.To develop a predictive model for heat stress events in pasture-based dairy systems, rumen temperature data from smaXtec rumen boluses in 443 cows on 3 dairy farms in Northland, New Zealand were used to identify heat stress events and these were matched with automated weather station data collected on or near the farm.Heat stress rate (HSR) was defined as the percentage of cows within an age-breed group having a rumen temperature greater than 3 standard deviations (SD) above an individual cows' mean and heat stress events were defined as HSR > 25%.Single and multiple linear regression models, including published heat stress indices, were able to predict a high proportion of heat stress events (sensitivity 61-68%), but were insufficiently discriminating, predicting also a high number of false positives (precision only 14-18%).A machine learning algorithm, cubist, was the best performing model, predicting 79% of heat stress events with a precision of 52% for this data set.Our proof-of-concept study demonstrates the potential of this approach, using climate data to predict and forecast heat stress events in pasture-based dairy systems.Further work should test the cubist model using independent data, refine data set construction, investigate the value of including known animal variables such as cow age or breed, and incorporate other measures of heat stress such as respiration rate.
W eather conditions in summer can lead to heat stress in dairy cattle with negative impacts on milk production and welfare, particularly in pasture-based dairy systems where there is less ability to manage the environment than in indoor systems.Heat stress risk is expected to be exacerbated as daily temperatures continue to rise and extreme hot weather events become more frequent (Jago et al., 2023).The ability to identify or predict the onset of herd or group level heat stress in dairy cows could enable operational management decisions (such as timing of milking or choice of paddock) to reduce its impact on-farm as well as to allow strategic management planning (such as changing milking frequency, breeding choices, tree planting and installation of portable shade structures) in response to the potential impact of climate change.
Previous research has investigated the relationship between weather conditions and respiration rate for grazing dairy cows (Bryant et al., 2022).Manual observation of animal respiration rate is the gold standard indicator for heat stress response; however, its collection is labor intensive, prone to operator error, performed only at discrete points in time, and not scalable to large numbers of animals, or in situations of high stocking density (Wijffels et al., 2021).This limits the ability to collect a sufficiently large data set across many environments for training and validation of predictive models.For example, the study of Bryant et al., (2022) used observations from only a single geographical region.Therefore, the use of automation would be beneficial to enhance the understanding of the relationship between animal responses to climatic conditions that could cause heat stress in a grazing environment.
Large data sets collected by sensor technologies provide an alternative approach to train models for predicting heat stress risk and determine the efficacy of management mitigations when heat stress conditions exist.One sensor-based approach is to use accelerometer data to predict heavy breathing, and the proportion of cows breathing heavily in a group has been shown to mirror changes in vaginal temperature (Bar et. al., 2019).Similarly, rumen temperature has been shown to be correlated with rectal temperature (Bewley et al., 2009;Boehmer, 2015) and has been used to evaluate the effects of heat stress (Lees et al., 2018).Once inserted, rumen temperature boluses record rumen temperature continuously, and don't require visual assessment of placement, so require less work for data collection over long periods than rectal or vaginal temperature loggers.However, the applicability of rumen temperature boluses to indicate heat stress in pasturebased dairy systems, and the weather conditions that drive it, are relatively unknown.The objective of this study was to explore the feasibility of identifying and predicting heat stress events in grazing dairy cows from automatically monitored weather conditions.We hypothesized that rumen bolus sensor data could be used to identify heat stress events and automated weather data could be used to develop a model to predict these events.
The study used "smaXtec" (smaXtec animal care GmbH, Graz, Austria) rumen bolus data from 443 cows located on 3 farms near Dargaville (Northland, New Zealand), from the period 1 Jan 2021 to 7 Jun 2023.Farm 1, had a subset of 192 cows that had smaXtec boluses from a herd of ~230 cows.Farm 2, located ~17 km NNE from Farm 1, had a subset of 49 cows with boluses from a herd of ~400 cows.Farm 3, located ~4 km ENE from Farm 2, had rumen temperature data from 259 unique cows from a herd of ~180 cows.On Farms 1 and 2 no new boluses were added to replacement animals entering the herd, which differed to Farm 3 where animals that entered the herd to replace culls, sales and deaths received boluses.At Farms 2 and 3, cows were always milked once per day in the morning.No specific heat stress mitigations (e.g., provision of sprinklers) were used on any of the farms.Institutional animal care and use or equivalent approval was not obtained since data were already recorded by the farms for their own use.
Data obtained from the smaXtec bolus included rumen temperature ("temp_without_drink_cycles"), which is the raw rumen temperature corrected by the manufacturer's proprietary algorithm for the effects of drinking events, as well as additional interpretive variables activity ("act"), rumination ("rum_index") and drinking events ("drink_cycles_v2"), stored at both 10-min and 1-h resolution.The current study focused on prediction using rumen temperature.The correction algorithm removed short-temp drops in rumen temperature presumed to be caused by drinking events, however the details of the proprietary algorithm are unknown.
Weather data for the same period was obtained from 2 Davis Wireless Vantage Pro2 Plus weather stations (Davis Instruments, 2023).One station was located at Farm 1, the other was located at Farm 2 until it was moved, for reasons relating to another project, to Farm 3 on 7 Jan 2023.Due to the proximity (4 km) of Farms 2 and 3 the weather was assumed to be the same at each site.Data was accessed via the WeatherLink v2 API (sensor numbers 43,52,53,56,242,243,504).Weather data was available at 15-min intervals and included air temperature (°C), humidity (%), solar radiation (MJ m −2 h −1 ), rainfall (mm), and wind speed (m s −1 ).Additional metrics, such as daily minimum air temperature, cumulative solar radiation, temperature-humidity index (Thom, 1959), and grazing heat load index (Bryant et al., 2022) were calculated from these.
Both data sets were stored in a Snowflake database (Snowflake Inc., Bozeman, MT, USA).Using the dbplyr package in R (R Core Team, 2023; Wickham et al., 2023), the weather data on the hour was joined to the hourly smaXtec bolus data for the 3 farms.The resulting table of 5.28 million rows (cows × hours) was downloaded and saved as a parquet file (size 113 Mb).Variables of interest (as listed above) were examined, and spurious values were replaced with "Not Available" (NA).Cows with evidently implausible activity (1 cow), drinking (1 cow) or rumen temperature (6 cows) data were assumed to have faulty sensors and their data were removed.
Four main steps were undertaken to analyze the data.1) A suitable indicator of heat stress incidence was proposed and heat stress events were identified.2) The data was balanced between heat stressed and non-heat stressed states (Branco et al., 2015).
3) Several regression models were tested for predicting the heat stress rate from weather data.4) These models were compared with existing heat load indexes for the prediction of heat stress events, particularly the grazing heat load index (Bryant et al., 2022).
Since we lacked respiration rate or panting data, heat stress incidence was defined using rumen temperature, which has been previously linked to thermal stress in dairy cattle (Donkersloot et al., 2017;Levit et al., 2021).The cumulative distribution of rumen temperature (hourly; adjusted for drinking) of individual cows varied in mean (gray vertical lines) and SD (Figure 1a).This highlighted that a fixed temperature threshold (e.g., 39.5 C, Liu et al., 2019;Levit et al., 2021) was unlikely to indicate heat stress across all animals.The smaXtec rumen temperature for each cow was therefore scaled so that all cows had a common mean and SD (Figure 1b).This provided the basis for a proposed heat stress rate (HSR) metric for a group or herd, which was defined as the percentage of cows with scaled rumen temperature > 3 SD above the mean (indicated by the vertical line in Figure 1b).Three standard deviations from the mean is often used to indicate "extreme" events (Grafarend, 2006), and was proposed as a proof of concept here.
In this study we chose to predict heat stress rate (i.e., HSR) from weather variables alone, so that in the future the output could be used to provide a regional forecast of expected heat stress events (proposed to occur when HSR exceeded 25%, see below).The key weather variables used were air temperature (AIR_C; °C), daily minimum air temperature (AIR_MIN; °C), daily cumulative mean solar radiation (SOL_CUM; MJ m −2 h −1 ), relative humidity (HUM_PC; %) and windspeed (WIND_MPS; m s −1 ). Figure 2 shows correlations between these weather variables and HSR.Other cow behavior variables (activity, ACT; rumination, RUM; drinking, DRINK) are also shown in Figure 2, as are the hour of the day (HOUR) and the month of the year (MONTH).Also shown are the temperature humidity index (THI) of Thom (1959) and the more recently developed grazing heat load index (GHLI) of Bryant et al. (2022), which is based on temperature, wind speed and solar radiation.THI is in units of temperature (°F) and GHLI is in units of respiration rate (min −1 ).
Based on this data, a threshold at HSR > 25% was proposed as a suitable indicator of heat stress events, as values of HSR > 0% and < 25% can occur under any conditions (e.g., at nighttime), whereas HSR > 25% observations generally matched expected heat stress conditions (high AIR_C, high AIR_MIN, high HUM_PC, high SOL_CUM, low WIND_MPS, low ACT, high DRINK, low RUM, afternoon HOUR, summer MONTH).On this basis, heat stress events were detected on 27, 23, and 35 d at Farms 1, 2 and 3 (out of 324, 330 and 330 d, respectively).
Models based on severely imbalanced training data tend to perform poorly when predicting the minority cases which are of primary interest (in this case, heat stress events; Branco et al., 2015).To avoid this, hours with air temperature ≤ 20°C were first excluded, as heat stress rarely occurs under these conditions (Bohmanova et al., 2007;Bryant et al., 2007).The remaining data were then resampled (Branco et al., 2015) to achieve similar numbers of positive (defined as HSR > 25%) and negative (defined as HSR ≤ 25%) data rows.The probability of retaining a data row was P = 100% when HSR was > 25% (806/806 rows), P = 5% + (100% -5%) × HSR ÷ 25% when HSR was in the range 0-25% (1,478/4,121 rows), and P = 5% when HSR was 0% (2,433/48,022 A range of models were tested for prediction of heat stress rate (HSR; Figure 3 and Rsq values in Table 1) and heat stress events (HSR > 25%; Table 1).These included simple linear regression (using variables such as AIR_C, THI or GHLI), multiple linear regression (lm), generalized additive models (gam; Wood, 2011) and machine learning models from the "caret" package in R (Kuhn, 2008), specifically random forest (rf), cubist, gradient boosted machine (gbm), support vector machine with radial basis (svmRadial), and k-nearest neighbors (knn), trained using the default settings.Three models were chosen for presentation here, and comparison with the GHLI.The models were trained using the 4,717 training cases, then validated against the entire data set of 52,949 cases (which were predominantly negative cases).
Figure 3a shows the model predictions of HSR over the entire data set, compared with the proposed threshold of 25% (red lines).Each data value (y-axis) is plotted against the corresponding model prediction (x-axis).For each model, the following statistics are reported (Table 1; Branco et al., 2015): the coefficient of determination ("Rsq"; the proportion of variance in the data explained by the model), accuracy ("Acc"; the proportion of correct predictions), sensitivity ("Sen"; the proportion of positive events that the model correctly predicted), precision ("Pre"; the proportion of predicted events that were, in fact, positive events according to the data; also known as the positive predictive value, PPV), and the F1 score, which is a summary of the model's ability to predict positive cases.The formulae for these metrics are: where HSR model,i are the model predictions corresponding to the data HSR data,i and HSR data is the mean value of the HSR calculated from the data.TP, TN, FP, FN are the number of true positives (HSR > 25%), true negatives, false positives and false negatives, respectively.Figure 3b shows the variable importance assessed using the "iml" package in R (Molnar et al., 2018), which is the increase in prediction error that occurs when each variable is removed.
While the data points themselves are not without error, the results in Table 1 clearly illustrate the improvement in training and predictive ability as model sophistication increased.Note that, because the training data were not a random sample of the full data set, the model performance statistics can only be compared between models, not between the training and predictive data sets.
The simple GHLI model (a linear model of cow respiration rate based on air temperature, windspeed and solar radiation) is proposed to predict heat stress events in grazing dairy cows, as indicated by marked changes in panting and drooling, when its value is 70 or more (Bryant et al., 2022).When applied to the current data set, the GHLI had a sensitivity of 72%, which was similar to the other models, but a precision of only 4%, indicating a high proportion of false positives (Table 1).Allowing for the different definitions of heat stress, this result indicates that the GHLI predicted several events where there was limited corresponding increase in rumen temperature (HSR ≤ 25%), indicating a need for further understanding about respiration rate as an indicator of when animals experience heat stress and how heat stress should be defined.Similarly, a multivariate linear model (lm) using AIR_C, AIR_MIN, SOL_CUM, HUM_PC and WIND_MPS performed only a little better.An additive nonlinear model using the same variables (gam) improved on this further.However, only machine learning models (e.g., cubist), that could fully exploit the structure of the data, gave predictions at an acceptable level of sensitivity and precision (i.e., detecting a high proportion of events, but without as many false positives; Figure 3 and Table 1).
Training Rsq and accuracy were greater than 0 and 0.50 respectively, as expected, for all models (Table 1).However, when the models were tested by predicting against the full data set, it became clear that Rsq and accuracy were not useful metrics of predictive performance.The large number of negative cases in the full data set resulted in low or negative Rsq values and misleadingly high values of accuracy.By comparison, sensitivity, precision and F1 gave more useful metrics of predictive performance.Only the performance of the cubist model was high enough that it may be useful for practical application (random forest "rf" performed similarly well; data not presented).Note that these results are robust to the selection of the threshold for HSR (HSR > 25%).Using a lower (for example) HSR threshold would affect both the proportion of positive cases and also the proportion of positive predictions, in contrast to examples in the literature where increases in prevalence are often associated with increases in sensitivity and decreases in precision (e.g., Leeflang et al., 2013).
As well as predictive accuracy, modeling allows us to assess the importance of different predictor variables for making predictions.The importance plot (Figure 3b) estimates the relative importance of variables in the multivariate models.The importance of variables varied slightly among the models, but air temperature, cumulative solar radiation and minimum daily air temperature were important in most models, with wind speed being less useful.Bryant et al., (2022) came to a similar conclusion in developing the GHLI index, albeit wind speed was concluded to be of more importance than humidity.
Animal susceptibility to heat stress is known to vary with breed and age (Bryant et al., 2007) and including these as predictor variables in the current models (results not shown here) indicated a weak association.Compared with the variables reported in Figure 3b they were generally the least important variables, and their inclusion improved model performance only slightly (sensitivity, precision and F1 of cubist each improved by 3-4%).Further work

Figure 1 .
Figure 1.Variation in smaXtec rumen temperature (a) for individual cows, (b) after scaling to a common mean, and (c) between cows in a herd at a point in time.Vertical lines in each facet indicate 3 SD above the mean, note this appears as shading in (a).Traces in (a) and (b) are colored by cow number and in (c) are colored by air temperature ranging from 0°C (blue) to 30°C (yellow).
Woodward et al. | PREDICTING HEAT STRESSrows).Resampling reduced the data set from 52,949 to 4,717 rows.While it would have been desirable to hold back some positive data rows for model testing, this would have resulted in too small a training set.