Modelling the association of heat stress with productive performance in Swiss dairy cows using large but sparse data

Although the link between heat stress and reduced performance is well documented, breed-specific differences in heat tolerance have not been extensively studied in grassland-based systems, and the computational challenges of modelling large-scale, sparse dairy data remain underexplored. Our objective was to quantify the association between temperature–humidity index (THI) and production performance in Swiss dairy systems, while demonstrating how real-world, high-dimensional, sparse dairy data can serve as a testbed to advance generalised additive mixed modelling (GAMM) tools. We compiled a dataset spanning 42 years, comprising over 130M monthly test-day records of milk and energy-corrected milk (ECM) yields from 46,975 farms and more than 4.2M cow-lactations across six breeds. Data were analysed using GAMM fitted separately for each breed to model non-linear effects of THI and days in milk, including fixed effects of parity and year, and random effects of cow, farm, and zip code nested within test-day records. The low annual sampling frequency per cow resulted in largely sparse, high-dimensional random effect matrices, presenting a substantial computational challenge for standard routines such as gam and bam in R. We developed an optimised hybrid solution based on gamm4 in R and MixedModels.jl in Julia, adapted to the long, wide, sparse data structures. Across breeds and parity, peak milk and ECM occurred at THI 47–58, with the magnitude of change in response variables from these peaks to THI 65–68 remaining small ( ∼ 0.12 – 1.50 % ), likely reflecting subtle physiological adjustments rather than immediate production losses. The optimised GAMM framework enabled robust estimation of these non-linear, biologically relevant responses in large-scale, sparse dairy data, providing a foundation for future studies integrating statistical modelling and animal physiology.