Bayesian hierarchical modelling for regional crop yield prediction under climate uncertainty

Introduction

Yield prediction models are only as useful as their honesty about what they don't know. A model that forecasts 6.0 t ha⁻¹ without stating whether it could equally be 5.0 or 7.0 is of limited value to insurers, commodity traders, or policy planners [1] . Yet most operational crop forecasting systems whether regression-based, machine-learning, or process-based report point estimates without formal uncertainty bounds [2] . Bayesian hierarchical models (BHMs) address this gap by treating parameters as probability distributions rather than fixed values [3] . The hierarchical structure allows information to be shared across spatial units (municipalities, counties) while respecting local differences a property called "partial pooling" that is well suited to agricultural data where nearby fields share climate but differ in soil and management [4] . Despite these advantages, BHMs have seen limited application in operational crop forecasting, largely because they are computationally expensive and require statistical expertise [5] . This research built and validated a BHM for winter wheat yield across southern Sweden and benchmarked it against three widely used alternatives.

Theoretical Background

In a Bayesian framework, unknown parameters θ are assigned prior distributions p (θ) reflecting existing knowledge. Observed data D update these priors via Bayes' theorem to produce posterior distributions p (θ|D) ∝ p (D|θ) p (θ). For hierarchical models, parameters at one level (e.g. municipal intercepts) are drawn from a shared higher-level distribution (e.g. regional mean), enabling partial pooling [3] . The posterior predictive distribution p(y*|D) integrates over all parameter uncertainty, producing probabilistic forecasts with credible intervals rather than single values. Inference was performed by Markov chain Monte Carlo (MCMC) sampling using the No-U-Turn Sampler (NUTS) in Stan [6] .

Simulation parameters

The BHM was specified as a two-level Gaussian model: y_ij = α_j + X_ijβ + ε_ij, where y_ij is yield in year i and municipality j, α_j ~ N(μ_α, σ²_α) are municipality random intercepts, X_ij is a vector of six covariates (April-June precipitation, July mean temperature, cumulative GDD, soil clay fraction, nitrogen application rate, and sowing date index), β are fixed-effect slopes, and ε_ij ~ N(0, σ²). Weakly informative priors were used: N (0, 10) for β, Half-Cauchy (0, 5) for variance parameters. Four MCMC chains of 4,000 iterations (2,000 warmup) were run; convergence was confirmed by R-hat < 1.01 and effective sample size > 400 for all parameters.

Material and Methods

Material

Municipal-level winter wheat yield data (2001-2020) for 48 municipalities in Götaland were obtained from Statistics Sweden (SCB). Climate data came from the Swedish Meteorological and Hydrological Institute (SMHI) gridded product (4 km resolution). Soil texture was extracted from the SGU national soil map. Nitrogen application rates were estimated from SCB fertiliser surveys.

Methods

Model performance was assessed by leave-one-year-out crossvalidation (LOYO-CV): each year in turn was held out for prediction while the remaining 19 years served as training data. Metrics included RMSE, MAPE, and coverage of 95% credible intervals. Three benchmark models were fitted on identical data: random forest (500 trees, mtry = 2), ridge regression (λ chosen by 10-fold CV), and the APSIM-Wheat process model calibrated to local cultivars. All statistical analyses were performed in R v4.2.1 with the brms package (Stan backend). The BHM achieved the lowest RMSE (0.31 t ha⁻¹) and highest R² (0.87) among all models (Table 1). Its 95% credible intervals captured observed yields in 93.5% of held-out cases, close to the nominal 95%. July temperature had a negative effect on yield (posterior mean -0.187 t ha⁻¹ per °C), consistent with heat stress during grain filling (Table 2). Nitrogen rate and spring precipitation were the strongest positive predictors.

Results

Discussion

The BHM's 30% RMSE improvement over random forest mirrors findings from Bayesian yield models in the U.S. Corn Belt [7] and France [8] , confirming that partial pooling across spatial units is particularly effective for agricultural data with strong spatial autocorrelation. The well-calibrated credible intervals are arguably more valuable than the point-estimate improvement, because they allow downstream users to quantify economic risk directly from the predictive distribution [9] . The negative July temperature coefficient reflects the sensitivity of Scandinavian winter wheat to heat stress during anthesis and early grain filling, a vulnerability projected to intensify under climate change [10] . The non-significant soil clay effect suggests that within Götaland's relatively narrow textural range, climate and management dominate yield variation.

Conclusion

A Bayesian hierarchical model for winter wheat yield in southern Sweden outperformed random forest, ridge regression, and APSIM in cross-validation accuracy and provided wellcalibrated probabilistic forecasts. The 95% credible intervals captured observed yields in 93.5% of held-out years. July temperature and nitrogen rate were the strongest yield predictors. BHMs offer a practical path toward probabilistic crop forecasting that honestly communicates prediction uncertainty to insurers, traders, and policy planners.

References

1. Asseng S, Ewert F, Martre P, et al. Applications, and accuracies. Adv Agron. 2019; 154:201-255.
2. Van Klompenburg T, Kassahun A, Bijman C. Crop yield prediction using machine learning: a systematic literature review. Comput Electron Agric. 2020; 177:105709.
3. Gelman A, Hill J. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge: Cambridge University Press; 2006.
4. Ogle K, Barber JJ, and Barron-Gafford GA, et al. Quantifying ecological memory in plant and ecosystem processes. Ecol Lett. 2015; 18(3):221-235.
5. Crespo Cuaresma J, Hlouskova J, Obersteiner M. Forecasting with Bayesian methods: an application to crop production. J Appl Stat. 2008; 35(5):571-582.
6. Carpenter B, Gelman A, Hoffman MD, et al. Stan: a probabilistic programming language. J Stat Softw. 2017; 76(1):1-32.
7. Peng B, Guan K, Tang J, et al. Towards a multiscale crop modelling framework for climate adaptation assessment. Nat Plants. 2020; 6(4):338-348.
8. Martre P, Wallach D, Asseng S, et al. Multimodel ensembles of wheat growth: many models are better than one. Glob Change Biol. 2015; 21(2):911-925.
9. Lobell DB, Burke MB. On the use of statistical models to predict crop yield responses to climate change. Agric For Meteorol. 2010; 150(11):1443-1452.
10. Asseng S, Ewert F, Martre P, et al. rising temperatures reduce global wheat production. Nat Clim Change. 2015; 5(2):143-147.
11. Holzworth DP, Huth NI, DeVoil PG, et al. APSIM evolution towards a new generation of agricultural systems simulation. Environ Model Softw. 2014; 62:327-350.
12. Schwalm CR, Anderegg WRL, Michalak AM, et al. Global patterns of drought recovery. Nature. 2017; 548(7666):202205.
13. Hoffman MD, Gelman A. The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J Mach Learn Res. 2014; 15(1):1593-1623.
14. Iizumi T, Yokozawa M, Sakurai G, et al. Historical changes in global yields: major cereal and legume crops from 1982 to 2006. Global Ecol Biogeogr. 2014; 23(3):346-357.