Assessing grain production resilience under climate stress: temporal holdouts and cross-regional validation

Abstract

Accurately identifying grain production resilience is a prerequisite for improving agricultural risk governance under climate change. The relationship between predictive accuracy for capacity-based indices and observed performance under climate shocks requires direct testing. Building on existing distinctions between potential and revealed resilience, this study compares capacity-based resilience with outcome-based resistance. We integrate grain production, socioeconomic, and climate data for Chinese prefecture-level cities, construct a provincial capacity-based composite index and a measure of resistance to climate shocks based on expected yield estimated using past information, and conduct a machine learning audit using strict temporal holdouts, transfer to nonoverlapping regions, nationwide replication, information-set ablation, probability calibration, and spatial residual diagnostics. The findings are threefold. First, capacity-based composite indices are readily predicted under random splits, but strict temporal and spatial validation substantially weakens this advantage. Direct tests of construct alignment reveal almost no stable correlation between capacity scores or their components and realized resistance. Second, the full city-level model has limited ability to predict continuous variation in realized resistance. Grain production history provides most of the information for below-trend classification, and probability calibration needs improvement. Third, aggregate discrimination is similar across the major producing areas, strict external sample, and nationwide replication, but annual AUC values differ. Comparisons of feature importance indicate similar patterns of information dependence, whereas spatial clustering of test residuals reveals regional prediction biases. We therefore argue that machine learning studies of agricultural resilience should assess construct validity, simple benchmarks, validation that prevents information leakage, external transfer, and calibration before interpreting features. High accuracy is substantively meaningful only when the measurement target aligns with the intended application.

Introduction

Climate variability is an important source of fluctuations in global crop yields (Ray et al., 2015), and climate extremes contribute to yield anomalies across major crops and regions (Vogel et al., 2019). Compound heat and moisture extremes pose additional risks to crop production under climate change (Lesk et al., 2022). Evidence also indicates that anthropogenic climate change has slowed global agricultural productivity growth, with effects varying across regions (Ortiz-Bobea et al., 2021). These findings motivate closer attention to how grain production systems withstand external shocks and maintain their core functions, an important concern for food security and agricultural risk governance.

Resilience provides a framework for analyzing responses to uncertainty. Development research links resilience to the maintenance of human wellbeing in the face of shocks (Barrett and Constas, 2014), while food system research emphasizes sustaining functions across interconnected social, economic, and biophysical processes (Tendall et al., 2015). Within agriculture, farm management research highlights buffering, adaptation, and transformation (Darnhofer, 2014), and farming system assessments distinguish system functions, resilience capacities, and resilience-enhancing attributes (Meuwissen et al., 2019). These perspectives motivate assessments of both available resources and performance under disturbance. For the resource-based composite index examined here, the question is whether higher capacity scores correspond to better production performance during climate shocks.

An agricultural system’s resource base is not fully equivalent to its performance under actual shocks. Zawalińska et al. (2022) explicitly distinguish potential from revealed resilience and examine the latter through productivity changes during crises in Polish agriculture. A review of farming system assessments notes that resilience is variously treated as an outcome, a capacity, or a process, and that conflating these meanings creates conceptual and methodological ambiguities (van der Lee et al., 2022). Building on this literature, we examine whether regional capacity indicators align with realized grain production resistance under climate stress. Greater mechanization, better infrastructure, and stronger economic capacity may improve the conditions for responding to risk, but these advantages do not establish how much production a system maintains during a particular disturbance. Our measurement question is therefore whether the capacity index used in this study captures performance under climate shocks, beyond the regional resource advantages that it summarizes.

This issue deserves particular attention in studies using composite resilience indices. Because these indices generally combine multiple economic, production, and resource variables through weighting, they contain substantial information about stable regional characteristics. Machine learning models used to predict such indices may achieve high accuracy. However, this accuracy need not indicate that the model has identified the capacity of an agricultural system to respond to external shocks. It may instead arise from information embedded in index construction, persistent regional differences, and inherent relationships among variables. In other words, machine learning can improve the ability to predict an indicator without automatically ensuring that the indicator has good construct validity.

Machine learning offers flexible methods for prediction with high-dimensional data and nonlinear relationships (Mullainathan and Spiess, 2017; Athey and Imbens, 2019). Its applications to crop yield prediction are documented in agricultural reviews (Elavarasan et al., 2018; van Klompenburg et al., 2020). However, predictive performance does not by itself establish explanatory validity (Shmueli, 2010). For regional agricultural data with temporal persistence and spatial dependence, random train–test splits may place observations from adjacent years in the same region into both training and evaluation samples. When the intended application concerns unseen years or regions, this overlap can lead to overly favorable performance estimates (Roberts et al., 2017). A key question in machine learning assessments of agricultural resilience is therefore whether high predictive accuracy transfers to unseen years and regions and corresponds to actual resistance to shocks.

This study therefore asks whether machine learning can identify grain production resilience rather than simply predict resilience scores determined by regional advantages and index design. To answer this question, we use grain production, socioeconomic, and climate data for Chinese prefecture-level cities to construct two representations, a capacity-based composite index and realized resistance under climate shocks. We first test their construct alignment directly at the province–year level. We then compare random-split results with simple trend benchmarks, strict temporal holdouts, leave-one-province-out validation, and transfer to nonoverlapping regions. Nationwide replication, information-set ablation, probability calibration, and spatial residual diagnostics are used to systematically evaluate model validity in unseen periods and regions and under different information conditions.

This study makes three main contributions. First, building on the distinction between potential and revealed resilience, it directly tests the alignment of a regional capacity index with realized grain production resistance under climate stress in China. These tests show that the composite score examined here cannot be interpreted as realized resistance without validation against shock outcomes. Second, it integrates simple benchmarks, temporal and spatial holdouts that prevent information leakage, nonoverlapping external validation, nationwide replication, probability calibration, and spatial residual diagnostics into a machine learning audit. This extends model evaluation beyond within-sample accuracy to construct validity, transferability, and deployment reliability. Third, information-set ablation and cross-sample comparisons of feature importance assess the contribution of city-level data to the evaluation of current production performance. By distinguishing predictive dependence, causal mechanisms, and policy actionability, the study provides reproducible empirical evidence to inform data development for agricultural resilience measurement and production risk assessment.

Theoretical framework

Agricultural production systems face disturbances arising from climate, resource constraints, and socioeconomic change. Farming system resilience concerns the maintenance of system functions through robustness, adaptability, and transformability (Meuwissen et al., 2019), within the wider food system perspective on sustaining functions under disturbance (Tendall et al., 2015). Because resilience is multidimensional and context dependent, translating it into measurable and verifiable indicators remains an important challenge.

We develop an analytical framework covering the resilience concept, indicator measurement, and identification through machine learning. The framework addresses whether predicted grain production resilience actually corresponds to an agricultural system’s ability to withstand climate shocks.

The meaning of agricultural resilience

The ecological foundations of resilience emphasize persistence and the ability to absorb disturbances while maintaining essential system relationships (Holling, 1973). Food system research extends this perspective to the maintenance of functions across social, economic, and biophysical processes (Tendall et al., 2015). Grain production is one component of this broader system, and the present study focuses specifically on production performance under climate stress.

Recent research emphasizes that agricultural resilience encompasses not only resistance to individual shocks but also capacities for adaptation and transformation. Meuwissen et al. (2019) argue that farming system resilience should be understood through robustness, adaptability, and transformability, and that assessments should be tied to specific disturbances and system functions. A crucial distinction remains, however, between the conditions available for responding to shocks and actual performance after shocks occur.

The distinction between potential and revealed resilience provides a basis for separating available capacities from observed performance (Zawalińska et al., 2022). This separation also responds to the need to distinguish capacity, process, and outcome in farming system assessments (van der Lee et al., 2022). In this study, capacity-based resilience emphasizes resources available before risks materialize, including agricultural inputs, infrastructure, technological conditions, and economic support. Outcome-based resistance emphasizes the actual production response to a specific disturbance, namely whether grain production remains stable during a climate shock. It captures a specific response dimension rather than all aspects of agricultural resilience. Higher development levels can improve the conditions for responding to risk, but whether those conditions translate into realized resistance depends on disturbance type, production structure, and adaptive behavior.

Potential bias in composite indicator measurement

Agricultural resilience has been assessed through different conceptual and measurement approaches. Meuwissen et al. (2019) distinguish system functions, resilience capacities, and resilience-enhancing attributes, whereas Zampieri et al. (2020) derive a production-based resilience indicator from crop production time series. These approaches should not be treated as interchangeable. For the resource-based composite index examined in this study, whether higher capacity scores correspond to better production performance during climate shocks remains an empirical question.

Substantial agricultural inputs, well-developed infrastructure, and strong economic capacity can indicate favorable conditions for responding to risk. They do not necessarily imply that a region can maintain stable grain output during extreme climate events. This raises a question of construct validity, for which convergent and discriminant validation provide an established methodological foundation (Campbell and Fiske, 1959). An indicator’s validity depends not only on its predictability but also on whether it corresponds to the concept defined by theory. Here, we test a specific implication by examining whether capacity scores align with production outcomes under climate stress. This comparison is a targeted test of construct alignment, rather than a full multitrait–multimethod validation. If the two measures do not align, a predictable capacity index should not automatically be interpreted as a measure of realized resistance.

Validation logic for identifying agricultural resilience with machine learning

Agricultural reviews document the use of machine learning to predict crop yields from agronomic and environmental information (Elavarasan et al., 2018; van Klompenburg et al., 2020). However, predictive performance does not by itself establish that a theoretical concept has been correctly measured, consistent with the distinction between prediction and explanation (Shmueli, 2010). If a resilience prediction target is itself constructed from model inputs, high accuracy may arise from structural relationships among variables rather than the identification of a system’s capacity to withstand shocks. Regional agricultural data may also exhibit temporal persistence and spatial correlation. Random cross-validation can place closely related observations in both training and test sets, producing overly favorable evaluations for applications to unseen years or regions. Roberts et al. (2017) motivate validation designs that account for such structure and match the intended prediction task. We therefore use temporal and spatial holdouts to assess transfer to the settings of interest.

Based on this reasoning, we develop a validation framework for identifying agricultural resilience through machine learning.

The framework comprises three steps.

First, capacity-based resilience indicators are constructed from regional resource conditions to test whether machine learning can predict existing composite resilience assessments. Second, yield deviations from historical benchmarks during climate stress events are used to measure realized resistance and assess the model’s ability to predict this production performance. Third, temporal extrapolation and spatial transfer tests evaluate the external validity of predictive performance.

Research design

Study area, data sources, and empirical tasks

This study examines China’s grain production systems, focusing on representations of resilience and their predictability under climate stress in major grain-producing areas. Chinese regions differ substantially in agricultural resource endowments, production structures, technological conditions, and climate exposure. These differences provide a rich empirical setting for comparing capacity-based composite resilience with actual shock responses and testing temporal and spatial extrapolation by machine learning models.

We define two related empirical tasks at different levels of analysis. The first is a provincial audit of capacity-based composite resilience. Using a panel of 13 major grain-producing provinces for 2011–2024, we calculate composite resilience scores from resistance, adaptation, and recovery dimensions using the entropy-weight method and assess their predictability under different validation designs. This task includes 182 province–year observations. The second task predicts city-level realized resistance to climate shocks. Historical grain yields before each target year provide a benchmark based on past information, and deviations of actual yields from that benchmark measure realized resistance during climate stress events. These tasks, respectively, capture the comprehensive capacity conditions of grain production systems and actual production performance under specific shocks, enabling comparison of their empirical meanings within a common framework.

The core city-level development sample comprises 170 cities in the 13 major grain-producing provinces. To assess spatial transferability, we designate in advance a strict external sample of 127 cities in 18 provincial-level administrative regions outside the major producing areas. These cities do not participate in training, validation, or parameter determination for the model developed in the major producing areas. Finally, we extend the entire analytical workflow to 297 cities nationwide to test reproducibility over a broader geographical scope. The city roster comprises 293 prefecture-level cities and four directly administered municipalities, excluding autonomous prefectures, leagues, and county-level cities administered directly by provinces. Throughout the paper, the nationwide sample refers to these 297 study cities (Figure 1).

Study area and distribution of model development and external validation cities. Blue circles indicate the 170 model development cities in the 13 major grain-producing provinces, and orange-red triangles indicate the 127 external validation study cities in 18 nonoverlapping provincial-level regions. Together, the two groups comprise 297 cities and represent the full study-city roster. Of these, 160 cities in the major producing areas and 117 external cities contribute P80 test events during 2021–2023. The basemap is the English standard map with approval number GS(2022)4308. City points form a separate research data layer, and the map retains its original aspect ratio.

The study draws on five main data sources. First, prefecture-level agricultural statistics cover total grain output, grain-sown area, and grain yield per unit area, primarily from the China City Statistical Yearbook and local statistical yearbooks. Second, a city-level socioeconomic panel provides information on agricultural inputs, production structure, economic development, fiscal expenditure, and science and technology expenditure. Third, daily NASA POWER weather data are used to construct city-level temperature, precipitation, and extreme climate indicators. Data requests use the coordinates of representative city points, rounded to 0.5° according to the established procedure. Fourth, an independent extreme climate database is used to assess sensitivity to the source of climate predictors while holding NASA POWER-identified P80 events fixed. Fifth, composite grain production resilience indices and their resistance, adaptation, and recovery components for the 13 major grain-producing provinces during 2011–2024 support the provincial capacity-based analysis and subsequent tests of construct alignment.

Raw prefecture-level grain yield data have 97.87% coverage in the major producing areas. The city-level socioeconomic database provides a complete city–year panel structure, although individual variables have varying degrees of missingness. Missing predictors are handled within the machine learning pipeline to ensure consistent processing rules across training, validation, and testing. Outcome construction requires an observed current yield and at least five valid historical yield observations.

Two measures of grain production resilience

The capacity-based composite resilience index

Grain production resilience concerns not only whether output remains stable when shocks occur but also the production system’s overall capacity to maintain production, adjust resource allocation, and restore functions under resource constraints, environmental change, and risk disturbances. Following this logic, we construct a provincial assessment framework comprising three dimensions, resistance capacity, adaptation capacity, and recovery capacity, with six criteria and 13 indicators.

Resistance capacity captures the foundations for maintaining basic production scale and stable supply. Its production resource conditions criterion includes grain-sown area and the effective irrigation rate, while production factor support includes cultivated-land area, primary-sector employment, and total agricultural machinery power. Adaptation capacity captures the ability to adjust to external change through improved production efficiency, resource allocation, and farmers’ livelihoods. Production efficiency adjustment includes grain yield per unit area and the grain multiple-cropping index, while farmers’ livelihood capacity includes rural per capita disposable income and agricultural science and technology investment. Recovery capacity captures the conditions for restoring production order and capacity and sustaining development after shocks and production fluctuations. Risk and shock constraints include grain-output volatility and the disaster-affected area rate, while fiscal and capital support includes fiscal expenditure on agriculture, forestry, and water affairs and fixed-asset investment in agriculture, forestry, animal husbandry, and fisheries.

To reduce the influence of subjective weighting, we use the entropy-weight method. Standardization parameters and entropy weights are estimated within the relevant training sample, frozen, and applied to holdout observations, after which the associated provincial results are reevaluated. Indicators measured in different units are first min–max normalized. For positive indicators,

z ijt = x ijt − min ( x j ) max ( x j ) − min ( x j )

For negative indicators,

z ijt = max ( x j ) − x ijt max ( x j ) − min ( x j )

where x ijt is the original value of indicator j for province i in year t , and z ijt is its standardized value. We then calculate each indicator’s information entropy, divergence coefficient, and objective weight.

p ijt = z ijt ∑ i , t z ijt , e j = − k ∑ i , t p ijt ln p ijt ,

d j = 1 − e j , w j = d j ∑ j d j .

where k = 1 / ln N , N is the number of province–year observations, and w j is the entropy weight of indicator j . The final composite resilience index is

CR I it = ∑ j = 1 13 w j z ijt .

Resistance, adaptation, and recovery component scores are obtained by summing the corresponding indicators in Table 1. The composite resilience score equals the sum of the three components. This yields 182 province–year observations for the 13 major grain-producing provinces during 2011–2024, used as the empirical representation of capacity-based composite resilience. We do not reconstruct the composite index at the prefecture-level city scale. City-level capacity and production structure information enters only as predictors of realized resistance.

Composite indicator framework for grain production resilience.

Dimension | Criterion | Indicator | Unit | Sign | Entropy weight
Resistance capacity | Production resource conditions | Grain-sown area (V1) | 1,000 ha | + | 0.108
Resistance capacity | Production resource conditions | Effective irrigation rate (V2) | % | + | 0.043
Resistance capacity | Production factor support | Cultivated-land area (V3) | 1,000 ha | + | 0.108
Resistance capacity | Production factor support | Primary-sector employment (V4) | 10,000 persons | + | 0.095
Resistance capacity | Production factor support | Total agricultural machinery power (V5) | 10,000 kW | + | 0.098
Adaptation capacity | Production efficiency adjustment | Grain yield per unit area (V6) | kg/ha | + | 0.035
Adaptation capacity | Production efficiency adjustment | Grain multiple-cropping index (V7) | % | + | 0.106
Adaptation capacity | Farmers’ livelihood capacity | Rural per capita disposable income (V8) | RMB 10,000 | + | 0.074
Adaptation capacity | Farmers’ livelihood capacity | Agricultural science and technology investment (V9) | RMB 100 million | + | 0.018
Recovery capacity | Risk and shock constraints | Grain-output volatility (V10) | % | − | 0.117
Recovery capacity | Risk and shock constraints | Disaster-affected area rate (V11) | % | − | 0.056
Recovery capacity | Fiscal and capital support | Fiscal expenditure on agriculture, forestry, and water affairs (V12) | RMB 100 million | + | 0.128
Recovery capacity | Fiscal and capital support | Fixed-asset investment in agriculture, forestry, animal husbandry, and fisheries (V13) | RMB 100 million | + | 0.014

“+” denotes a positive indicator, and “−” denotes a negative indicator.

Realized resistance under climate stress

The capacity-based composite index captures the overall foundations of grain production systems and the conditions for adjustment and recovery. To examine actual production performance under specific climate pressures, we construct an outcome-based resistance measure from grain yield deviations.

For city i and target year t , we fit a local linear Huber trend using only the five most recent years with valid grain yield observations before the target year. Let W it denote this set of five valid historical years. Then,

ln ( Y i s ) = α it + β it s + ε is , s ∈ W it ,

The yield benchmark for the target year is

Y ̂ it = exp ( α ̂ it + β ̂ it t ) .

This benchmark is a one-step-ahead trend extrapolation based strictly on past information. The historical information set updates annually. For example, the benchmark for 2022 may use observed yield in 2021, and the benchmark for 2023 may use observed yield in 2022. Information from the target year or any later year never enters estimation of that year’s benchmark.

Continuous realized resistance is defined as

R it = 100 × ( Y it Y ̂ it − 1 ) .

A higher R it indicates better realized yield relative to the benchmark based on past information. The measure converts production performance during climate stress events into comparable yield deviations, while also reflecting crop composition, production inputs, other disturbances, and estimation of the historical benchmark. We further define two binary outcomes.

B it = 1 ( R it < 0 ) ,

S it = 1 ( R it < − 5 ) .

Here, B it indicates that realized yield is below the benchmark based on past information, and S it indicates a yield decline of more than 5% relative to that benchmark. The Huber loss function reduces the influence of local yield outliers on trend fitting, improving the robustness of the historical benchmark.

Climate stress event identification and model information sets

Climate stress indicators

Daily NASA POWER weather data are used to construct city-level climate stress indicators. For consistent comparisons across cities, all underlying variables used to construct climate stress are measured from April through October, capturing exposure within this common seasonal window. Exposure outside this window, including winter wheat overwintering and early spring conditions, remains an area for extending the indicators. Heat stress is characterized by the number of days with maximum temperatures of at least 30 °C and at least 35 °C, and the number of hot nights with minimum temperatures of at least 20 °C. Cold stress is measured by frost days with minimum temperatures below 0 °C. Drought stress combines cumulative April–October precipitation and the maximum number of consecutive dry days, defined as days with precipitation below 1 mm. Excess precipitation stress combines the number of heavy-rain days with at least 20 mm of precipitation and maximum consecutive five-day precipitation (Rx5day).

To make climate exposure comparable across cities, each climate variable is first standardized within cities. For city i , climate indicator k , and year t , the 2008–2018 baseline period defines

z ikt = X ikt − μ i , k , 2008 − 2018 σ i , k , 2008 − 2018 .

The mean and sample standard deviation are calculated from valid observations for the city during the baseline period. For city–indicator series with a standard deviation no greater than 10 − 12 , the denominator is set to one. Deviations are therefore zero when the series is constant during the baseline period, while subsequent changes enter as differences from the baseline mean. If a current climate value is missing, its standardized deviation is set to zero when constructing climate stress. Four nonnegative stress measures are then defined.

Hea t it = max [ z ( Day s ≥ 30 ∘ C ) , z ( Day s ≥ 35 ∘ C ) , z ( HotNights ) , 0 ] ,

Col d it = max [ z ( FrostDays ) , 0 ] ,

Drough t it = max [ − z ( Precipitation ) , z ( ConsecutiveDryDay s ) , 0 ] ,

We t it = max [ z ( HeavyRainDays ) , z ( Rx 5 day ) , 0 ] .

Composite city–year climate stress is the maximum of these four measures.

C it = max ( Hea t it , Col d it , Drough t it , We t it ) .

The P80 and P90 thresholds are determined from the distribution of city–year climate stress scores among study cities during the 2008–2018 baseline period and held fixed thereafter. An event occurs when the score reaches the corresponding threshold. In the major producing areas, the P80 and P90 thresholds are 1.978655 and 2.284024, respectively. Strict external testing uses these same thresholds. Nationwide replication reestimates thresholds from the baseline distribution for all 297 cities, yielding P80 and P90 values of 1.959824 and 2.279158. Each city’s climate standardization parameters are still estimated from its own 2008–2018 observations. Fixed historical thresholds allow subsequent event frequencies to reflect changes relative to the baseline climate distribution. P80 defines events in the main analysis, while P90 is used for sensitivity analysis. Among records meeting the P80 stress criterion with both an available rolling benchmark and observed current yield, the major producing areas contribute 752 city–year climate stress events during 2013–2023, comprising 228 in 2013–2018, 174 in 2019–2020, and 350 in 2021–2023. Among all city–year observations in each period, the shares flagged as climatic P80 events are 22.9%, 52.9%, and 69.2% in training, validation, and testing, respectively. The corresponding shares entering the final modeling samples are 22.4%, 51.2%, and 68.6%.

Predictors and information sets

To align predictive information with the level of analysis in each task, we construct separate information sets for provincial composite resilience and city-level realized resistance.

The provincial capacity-based composite resilience task uses 13 province–year predictors, including agricultural plastic film intensity, fertilizer intensity, pesticide intensity, the agricultural producer price index, rural per capita disposable income, GDP per capita, grain yield per unit area, mechanization, agricultural research expenditure, the labor input ratio, disaster intensity, the multiple-cropping index, and rural fixed-asset investment. This task examines whether models can learn the stable structural information embedded in the composite resilience index under different validation designs.

The full city-level realized resistance model uses 26 predictors from three information categories. The first comprises 13 contemporaneous climate variables describing climate exposure and extreme weather, with definitions, units, and transformations reported in Table A6. The second comprises three yield-history variables generated solely from past information, namely lagged yield, the lagged yield anomaly, and expected yield. The third comprises 10 lagged capacity and production structure variables, including the share of grain-sown area, agricultural machinery power per unit of grain-sown area, fertilizer intensity, rural electricity-use intensity, the share of agriculture in the gross output value of agriculture, forestry, animal husbandry, and fisheries, GDP per capita, the primary-sector share of GDP, the share of science and technology expenditure in fiscal expenditure, primary-sector employment, and fiscal expenditure.

Because the city model uses April–October climate information from the target year, it is intended for assessing current production risk after that information becomes available. To identify the incremental predictive value of different sources, we conduct information-set ablations that exclude the lagged yield anomaly, exclude all yield-history variables, retain climate variables only, or retain capacity and production structure variables only. Comparing performance over the identical test period reveals whether useful signals primarily come from climate exposure, historical production trajectories, or regional capacity and production structure.

Machine learning models and validation strategies

The provincial composite resilience task

The provincial capacity-based composite resilience index is predicted using several continuous-outcome machine learning algorithms, including Ridge, Random Forest, Extra Trees, Histogram Gradient Boosting, and Gradient Boosting with Huber loss. Evaluation extends beyond fit among observations drawn from similar distributions to temporal extrapolation and spatial transfer.

We use four validation and benchmark comparisons. First, random row splits measure interpolation across similar provinces and adjacent years. Second, two simple benchmarks, province-specific linear trends and the last observed value, assess the incremental predictive value of complex models over trends and persistence. Third, strict temporal holdouts assign 2011–2018 to training, 2019–2020 to validation, and 2021–2024 to testing. Fourth, leave-one-province-out validation evaluates generalization to provinces excluded from training. Appendix Table A1 reports the years, sample sizes, and model selection arrangements for each design.

Beyond predictive performance, we test construct alignment between capacity-based composite resilience and city-level realized resistance. Composite scores and the three component scores for the 13 major grain-producing provinces are matched to mean realized resistance among cities experiencing P80 climate stress events in the same province and year. Spearman rank correlations measure alignment in their rankings. Confidence intervals are obtained from 2,000 bootstrap replications clustered by province, and association models control for province and year fixed effects. This analysis assesses the empirical correspondence between capacity-based composite resilience and actual production performance during shocks.

Two additional association tests with predefined specifications assess sensitivity to aggregation and temporal alignment. First, within the set of P80 event cities E pt in each province–year, current grain-sown area A it provides weights for R ¯ pt A = ∑ i ∈ E pt A it R it / ∑ i ∈ E pt A it , which is related to the contemporaneous provincial capacity score. Second, the equally weighted city mean of realized resistance is retained, while the capacity score is replaced by its value in the previous calendar year, C p , t − 1 . Lags are constructed in the complete 2011–2024 provincial panel before matching to event outcomes during 2013–2023, so the capacity scores refer to 2012–2022. Both tests report the composite score and all three components, and matched province–year observations retain equal weights. Resistance is not assigned a value of zero in province–years without events. For each additional correlation, a 95% confidence interval is obtained from 2,000 bootstrap draws of provinces with replacement, retaining all matched years of each selected province and recalculating the rank correlation in every replication. The 2.5th and 97.5th percentiles define the interval. Both tests use the same bootstrap draws with random seed 20260920. For the composite score, we also estimate association models with province and year fixed effects. The additional tests use standard errors clustered by province and inference based on a t distribution with 12 degrees of freedom. These are sensitivity analyses of aggregation and temporal alignment within the event sample, not new prediction tasks. The area-weighted mean is not equivalent to resistance in total provincial grain output, and lagging the indicator does not itself change how the existing capacity scores were constructed.

The city-level realized resistance task

Candidates for continuous realized resistance include Ridge, Random Forest, Extra Trees, Histogram Gradient Boosting, and Gradient Boosting with Huber loss. Candidates for the binary below-trend and severe-loss outcomes include weighted Logit, Random Forest, Extra Trees, Histogram Gradient Boosting, and Gradient Boosting classifiers. Weighted Logit uses L2 regularization with C = 0.5 and balanced class weights (class_weight = “balanced”). Within each fitting sample, an observation in class c receives weight n / ( 2 n c ) , where n is the fitting sample size and n c is the number of observations in that class. Class weights are recalculated using the development sample when the model is refitted.

The city-level task uses a predefined temporal split, with 2013–2018 for training, 2019–2020 for validation, and 2021–2023 for a single final test. Continuous models are selected by validation RMSE, and classifiers by validation ROC AUC. For classification, probabilities predicted for the 2019–2020 validation sample by models trained on 2013–2018 are also used to select a classification threshold. Thresholds from 0.10 to 0.90 in increments of 0.005 are evaluated, and the threshold maximizing validation balanced accuracy is chosen, with ties resolved in favor of the smaller threshold. Once the model and threshold are selected, the model is fitted once more on the 2013–2020 development sample. The numerical threshold fixed during validation is applied directly to the refitted model’s 2021–2023 test probabilities. Models are selected by AUC, while probability error and calibration are evaluated separately. During testing, historical observed yield information may update as years advance, but all other model parameters, preprocessing parameters, and classification thresholds remain fixed.

The strict temporal splits for both tasks are summarized in Table 2.

Temporal splits for the two core tasks.

Empirical task | Unit of analysis | Training period | Validation period | Test period
Provincial capacity-based composite resilience | Province–year | 2011–2018 (104 observations) | 2019–2020 (26 observations) | 2021–2024 (52 observations)
City-level P80 realized resistance | City–year climate stress event | 2013–2018 (228 events) | 2019–2020 (174 events) | 2021–2023 (350 events)

Preprocessing for the city task is embedded in the modeling pipeline. During model selection, median imputation, missingness indicators, and any required robust scaling parameters are estimated from the training sample and applied unchanged to validation observations. After the final model is selected, preprocessing parameters are reestimated once on the 2013–2020 development sample and then held fixed for both the major-producing-area test sample and the strict external sample. This separates model development from final testing. The selected continuous model for the major producing areas, Extra Trees, uses 500 trees, a maximum depth of six, a minimum of eight observations per leaf, and a feature sampling fraction of 0.75. The nationwide selected Ridge model uses a regularization parameter of α = 10. Weighted Logit permits up to 3,000 iterations, and all candidate and final models use random seed 20260730.

External validation and model diagnostics

Transferability should be assessed against the spatial or temporal setting in which a model will be used, a general validation principle discussed for ecological models by Wenger and Olden (2012). To assess spatial transferability, the model trained and finalized in the major producing areas is applied directly to 127 external cities. These cities are located in 18 provincial-level administrative regions that do not overlap with the 13 major grain-producing provinces. Of these, 117 cities experience at least one P80 climate stress event during 2021–2023, yielding 237 external test events. None of these cities’ information is used in model development.

We then replicate the full workflow for all 297 cities nationwide, repeating training, validation, and final fitting. The nationwide training and validation periods contain 367 and 259 P80 events, respectively. The 2021–2023 test sample covers 278 cities and 595 P80 climate stress events. The nationwide city roster includes both the major-producing-area and external cities, but event membership is redetermined using nationwide thresholds. Strict external testing evaluates geographical transfer of a fixed model, whereas nationwide replication evaluates the complete workflow after expanding city coverage. The two provide complementary assessments.

Discrimination and calibration describe distinct aspects of predictive performance, a general methodological distinction emphasized in clinical risk prediction by Van Calster et al. (2019). For continuous outcomes, we primarily report R 2 , RMSE, MAE, and Spearman correlations. For binary outcomes, we report ROC AUC, average precision (AP), the Brier score, expected calibration error (ECE), and calibration intercepts and slopes, together with balanced accuracy, sensitivity, and specificity at the fixed classification threshold. No post hoc probability recalibration is performed on the final test sets.

The 95% confidence intervals for classification performance use a percentile bootstrap clustered by city. Each replication samples cities with replacement, drawing the same number of cities as in the original test sample. All test events and frozen predictions for each selected city enter the resample. The procedure is repeated 2,000 times, and the 2.5th and 97.5th percentiles define the intervals. The below-trend analysis uses random seed 20260920. Severe-loss intervals cover three R < −5% endpoints, the major-producing-area P80 sample, the major-producing-area P90 sample, and the nationwide P80 sample. These are resampled consecutively in that order using random seed 20260927, with 2,000 valid replications for each endpoint (Table 3; Appendix Table A4). Each interval characterizes city-sampling variation conditional on a fixed model. We do not conduct inference on the difference between the P80 and P90 AUCs.

Sensitivity to climate stress event thresholds.

Panel A: Continuous deviations and below-trend classification
Event threshold | Cutoff τ | Training N | Validation N | Test N | Continuous R2 | Below-trend AUC | AP | Brier
P80 | 1.9787 | 228 | 174 | 350 | −0.081 | 0.621 | 0.693 | 0.261
P90 | 2.2840 | 113 | 143 | 290 | 0.024 | 0.609 | 0.719 | 0.282

Panel B: Severe-loss classification and estimation precision
Event threshold | Test N | Cities | Severe-loss events | Severe-loss AUC | 95% CI
P80 | 350 | 160 | 32 | 0.597 | [0.492, 0.702]
P90 | 290 | 148 | 28 | 0.640 | [0.526, 0.764]

τ is determined from the distribution of climate stress scores for all city–year observations in the major producing areas during 2008–2018 and held fixed thereafter. The P90 event sample undergoes the same complete training, validation, and testing procedure as P80. Both rows use the full 26-variable information set. AP and Brier in Panel A refer to below-trend classification, and severe loss in Panel B is defined as R < −5%. Validation-selected P90 classification thresholds are 0.870 for below-trend outcomes and 0.520 for severe losses. Test prevalence for below-trend outcomes is 0.583 under P80 and 0.597 under P90. Severe-loss intervals use saved test predictions, 2,000 city-clustered bootstrap resamples, the percentile method, and random seed 20260927. All test events of each selected city are retained. Each of the three existing prediction endpoints yields 2,000 valid replications; the nationwide endpoint interval appears in Appendix Table A4. P80 is the main event definition, while P90 describes predictive performance under a higher climate stress threshold.

For model interpretation, permutation importance is calculated on the fixed test samples to characterize each fixed model’s dependence on individual predictors and compare importance rankings between the major producing areas and the nationwide replication. Single-variable permutation shuffles one variable while holding the others unchanged. Results are therefore interpreted alongside predictor correlations, information-set ablation, and historical-information benchmarks. Because the two samples overlap in city coverage, their ranking comparison describes dependence patterns after coverage is expanded.

Finally, a five-nearest-neighbor spatial weight matrix based on city coordinates is used to examine the spatial structure of prediction errors. Global Moran’s I is calculated separately for realized resistance, prediction residuals, and predicted probabilities, with statistical significance assessed using 499 random permutations. Residual spatial correlation diagnoses regional error structure, while spatial correlation in predicted probabilities describes the geographical distribution of risk scores. Together, these analyses form a linked validation framework covering composite capacity measurement, shock outcome measurement, temporal extrapolation, spatial transfer, probability calibration, and spatial residual diagnostics.

Empirical results

Capacity-based composite resilience, predictability, and construct alignment

Predicting composite resilience under different validation designs

Predictive performance for provincial capacity-based composite resilience varies substantially with the validation design. Under an 80%/20% random row split with random seed 42, Gradient Boosting with Huber loss achieves an R2 of 0.957 on 37 held-out province–year observations. With a strict temporal split using 2011–2018 for training, 2019–2020 for validation, and 2021–2024 for testing, the validation-selected Ridge model yields an R2 of 0.805 on 52 test observations. On the same temporal test sample, the province-specific linear trend benchmark achieves an R2 of 0.869, while the persistence benchmark that carries each province’s 2020 value forward throughout the test period achieves 0.809. Neither performs worse than the selected machine learning model. For prediction in unseen provinces, the pooled R2 from leave-one-province-out validation is 0.531 (Figure 2a).

Predictability and construct alignment of capacity-based composite resilience. Panel a reports R2 under random row splitting, two simple temporal benchmarks, strict temporal holdout, and leave-one-province-out validation. Green indicates the two simple temporal benchmarks, and blue indicates the machine learning audit results. The random split uses Gradient Boosting with Huber loss and seed 42, while temporal and spatial extrapolation use the validation-selected Ridge model. Panel b plots composite scores against mean realized resistance among event cities in the same province and year. Each point represents an equally weighted province–year observation. Panel c reports Spearman correlations between realized resistance and the composite score and its three components. Error bars are 95% confidence intervals from 2,000 bootstrap replications clustered by province. Appendix Table A1 provides comparisons using the same algorithm and the years, sample sizes, and model selection arrangements for each design.

These findings indicate that composite resilience scores contain substantial learnable structure, but high performance under random splitting does not directly extend to temporal or spatial extrapolation. Over the same future test period, province-specific linear trends outperform the selected machine learning model, while persistence performs similarly. Temporal extrapolation and spatial transfer characterize performance in unseen years and provinces, respectively. Appendix Table A1 provides comparisons holding the Ridge algorithm fixed, results across random seeds, and the complete validation arrangements.

Construct alignment between capacity-based composite resilience and realized resistance

Matching composite scores to mean realized resistance among P80 event cities in the same province and year yields 120 province–year observations from 13 provinces. In the baseline specification, the Spearman correlation between composite resilience and realized resistance is 0.022, with a 95% confidence interval of [−0.102, 0.162] from 2,000 bootstrap replications clustered by province. Correlations for the resistance, adaptation, and recovery components are 0.038, −0.042, and 0.095, respectively, and all confidence intervals include zero (Figures 2b,c). After controlling for province and year fixed effects, the baseline composite-score coefficient is 11.051 (province-clustered SE = 15.754; 95% CI [−23.275, 45.377]; p = 0.496), providing no statistically detectable association with realized resistance (Table 4). These results show that the resource, adjustment, and recovery conditions captured by capacity-based composite resilience contain different empirical information from actual production performance during specific climate stress events.

Two-way fixed effects associations under the baseline, area weighting, and lagged capacity specifications.

Specification | Coef. | Province-clustered SE | 95% CI | p (t₁₂) | Province–year N | Province N
Baseline, equally weighted city resistance × current composite score | 11.051 | 15.754 | [−23.275, 45.377] | 0.496 | 120 | 13
Area-weighted realized resistance × current composite score | 5.332 | 13.239 | [−23.513, 34.176] | 0.694 | 120 | 13
Equally weighted city resistance × one-year-lagged composite score | 12.316 | 23.951 | [−39.869, 64.501] | 0.616 | 120 | 13

All three models control for province and year fixed effects. Standard errors are clustered by province, and intervals and p-values use the t(12) distribution. This analysis assesses the robustness of construct alignment to aggregation and temporal alignment.

Two additional tests assess the sensitivity of this construct alignment result to within-province aggregation and temporal alignment. First, realized resistance is reaggregated using grain-sown area as weights among event cities in each province–year. Second, provincial capacity scores are replaced with their values in the previous calendar year, while current realized resistance remains an equally weighted city mean. Both tests retain the same 120 matched province–year observations to ensure direct comparability across specifications.

Table 5 reports rank associations under the three specifications. With area weighting, the Spearman correlation between the composite score and realized resistance is 0.036, with a 95% confidence interval of [−0.080, 0.159]. With capacity lagged by 1 year, the correlation is 0.008, with a 95% confidence interval of [−0.124, 0.153]. Point estimates for the three components also remain low, as in the baseline. Alternative aggregation and temporal alignment thus produce consistent rank patterns, supporting the treatment of capacity-based composite resilience and realized resistance during shocks as related but distinct empirical dimensions.

Construct alignment under alternative aggregation and temporal alignment.

Capacity measure | Baseline, equal city weights × current capacity | Area weights × current capacity | Equal city weights × one-year-lagged capacity
Composite resilience score | 0.022 [−0.102, 0.162] | 0.036 [−0.080, 0.159] | 0.008 [−0.124, 0.153]
Resistance capacity score | 0.038 [−0.111, 0.182] | 0.025 [−0.119, 0.152] | 0.065 [−0.091, 0.228]
Adaptation capacity score | −0.042 [−0.195, 0.082] | 0.001 [−0.133, 0.111] | −0.082 [−0.216, 0.040]
Recovery capacity score | 0.095 [−0.051, 0.217] | 0.100 [−0.044, 0.226] | 0.079 [−0.051, 0.199]

Entries report Spearman ρ and 95% confidence intervals from 2,000 bootstrap replications clustered by province. All three specifications use the same 120 province–year observations. Area weights use event cities’ current grain-sown area. The one-year-lagged specification matches provincial capacity scores from the previous calendar year to current event outcomes.

Models with province and year fixed effects provide evidence consistent with the rank associations in both sensitivity specifications (Table 4). The composite-score coefficient is 5.332 when realized resistance is area weighted and 12.316 when capacity is lagged by 1 year. Both 95% confidence intervals cross zero. Together with the rankings in Table 5, these results show weak construct convergence under alternative aggregation and temporal alignment, reinforcing the robustness of the baseline conclusion.

Event counts and city coverage vary across the 120 matched province–years, but both additional tests retain the entire matched sample. Appendix Table A2 reports event counts, the number of cities with observed outcomes, and coverage proportions for each province–year.

Realized resistance under climate stress, strict future prediction, and information sources

Diagnosing the realized resistance benchmark

When a fixed Huber trend estimated over 2008–2018 is extrapolated forward, the Pearson correlation between current and lagged yield anomalies in the major-producing-area event sample is 0.653 during validation and 0.906 during testing. It reaches 0.920 in the nationwide test sample. Using a rolling Huber benchmark based on the five most recent valid yield observations before each target year reduces these correlations to 0.130, 0.074, and 0.051, respectively (Figure 3a). Thus, strong residual persistence under fixed-trend extrapolation is sensitive to the benchmark specification and weakens substantially when the benchmark updates using past information. These diagnostics use observations with both current and lagged outcomes available, so their sample sizes are smaller than those used for full prediction, as detailed in Appendix Table A3.

Realized resistance benchmarks and sources of predictive information. Panel a reports Pearson correlations between current and lagged anomalies under the fixed training-period Huber benchmark and the rolling five-observation Huber benchmark. The three comparisons are, in order, major-producing-area validation, major-producing-area testing, and nationwide testing. Panel b reports below-trend test AUCs for five information sets, with variable counts in parentheses. Models are selected separately for each information set using the same training and validation procedure and identical test periods. The dashed line marks AUC = 0.5. Blue denotes the major producing areas, and green denotes nationwide replication.

Model selection during validation

Table 6 reports candidate performance after training on 2013–2018 and application to the 174 validation events in 2019–2020. For the continuous outcome, Extra Trees has the lowest RMSE at 7.520, closely followed by Ridge at 7.523. Both have validation R2 values slightly below zero. For below-trend classification, class-balanced weighted Logit has the highest validation AUC at 0.667 and is selected under the predefined AUC rule. Its Brier score ranking differs from its AUC ranking, illustrating the distinction between risk ranking and probability accuracy.

Comparison of candidate models during validation.

Panel A: Continuous realized resistance
Model | Validation RMSE | Validation MAE | Validation R2
Ridge | 7.523 | 5.055 | −0.030
Random Forest | 7.904 | 5.236 | −0.138
Extra Trees (selected) | 7.520 | 5.043 | −0.030
Histogram Gradient Boosting | 8.259 | 5.735 | −0.242
Gradient Boosting with Huber loss | 7.822 | 5.307 | −0.114
Mean benchmark (excluded from model selection) | 7.638 | 4.696 | −0.062

Panel B: Below-trend classification
Model | Validation AUC | Validation AP | Validation Brier
Weighted Logit (selected) | 0.667 | 0.572 | 0.263
Random Forest | 0.604 | 0.575 | 0.260
Extra Trees | 0.614 | 0.558 | 0.243
Histogram Gradient Boosting | 0.543 | 0.466 | 0.365
Gradient Boosting | 0.577 | 0.515 | 0.325
Prevalence benchmark (excluded from model selection) | 0.500 | 0.379 | 0.291

Training covers 2013–2018 (228 events), and validation covers 2019–2020 (174 events). Continuous models are selected by the lowest validation RMSE and classifiers by the highest validation AUC. Bold type denotes selected models. The mean and prevalence benchmarks are estimated from the training sample and used only for comparison, not candidate selection. RMSE and MAE are measured in percentage points of realized resistance. Observation-level validation predictions and complete performance metrics for all models have been saved separately. This table does not report test performance.

The shares of below-trend outcomes in the training, validation, and test periods are 61.4%, 37.9%, and 58.3%, respectively. The classification threshold of 0.585 is selected from the training model’s validation predictions. The model is then refitted once on the 2013–2020 development sample, and the threshold is held fixed for testing. Period-specific outcome frequencies and final probability diagnostics jointly inform the assessment of transfer performance under this development procedure.

Strict temporal test performance in the major producing areas

Across 350 P80 events in 160 cities in the major producing areas during 2021–2023, the continuous model achieves an R2 of −0.081, an RMSE of 4.618 percentage points, and an MAE of 2.965 percentage points. The Spearman correlation between predicted and observed values is 0.172. The scatterplot shows predictions concentrated within a narrow range, with large positive and negative deviations inadequately captured (Figure 4a). Relative to the test-sample mean used as the reference for R2, continuous prediction provides no positive explanatory gain.

Continuous prediction, risk discrimination, and probability calibration under strict testing. Panel a plots observed and predicted resistance for the 350 test events in the major producing areas, with the dashed line indicating perfect prediction. Panel b shows below-trend ROC curves for the major-producing-area, strict external, and nationwide samples, with random ranking as the dashed reference. Panel c plots mean predicted probability against observed prevalence for 10 equal-frequency probability groups in the same samples, with the dashed line indicating ideal calibration. Blue circles, orange-red squares, and green diamonds denote the three samples, respectively. Curves are used solely for diagnosis and do not inform model updates.

For below-trend classification, AUC is 0.621 and AP is 0.693, compared with an observed below-trend prevalence of 0.583. A city-clustered bootstrap of the fixed test predictions with 2,000 replications gives an AUC 95% confidence interval of [0.561, 0.678]. At the fixed threshold of 0.585, balanced accuracy is 0.570, sensitivity is 0.235, and specificity is 0.904. Only 48 of the 204 below-trend events are identified. The full model retains some risk-ranking information, but high specificity and low sensitivity at the fixed threshold indicate that it is better at ruling out below-trend outcomes than identifying them. Table 7 consolidates the complete cross-sample metrics and probability diagnostics.

Predictive performance and probability calibration in the three fixed test samples.

Panel A: Continuous realized resistance
Sample | Event N | Cities | R2 | RMSE | MAE | Spearman ρ
Major producing areas, strict test | 350 | 160 | −0.081 | 4.618 | 2.965 | 0.172
Strict external test | 237 | 117 | −0.152 | 7.983 | 4.671 | 0.266
Nationwide replication | 595 | 278 | −0.120 | 6.194 | 4.061 | 0.368

Panel B: Below-trend classification and probability calibration
Sample | N | Prev. | AUC | AP | Brier | ECE | Cal. intercept1 | Cal. slope2 | Bal. acc. | Sens. | Spec.
Major producing areas, strict test | 350 | 0.583 | 0.621 | 0.693 | 0.261 | 0.159 | 0.723 | 0.592 | 0.570 | 0.235 | 0.904
Strict external test | 237 | 0.616 | 0.613 | 0.701 | 0.305 | 0.240 | 1.240 | 0.165 | 0.563 | 0.247 | 0.879
Nationwide replication | 595 | 0.595 | 0.618 | 0.700 | 0.244 | 0.106 | 0.397 | 0.457 | 0.575 | 0.407 | 0.743

The test period is 2021–2023 for all three samples. Strict external testing directly uses the model, preprocessing, P80 event threshold, and classification threshold fixed in the major producing areas. Nationwide analysis repeats the same workflow on the expanded study sample and reestimates the historical P80 threshold. Classification thresholds are 0.585 for the major-producing-area and external samples and 0.570 nationwide, each determined during validation in the corresponding development workflow. The nationwide test sample includes the major-producing-area and external events plus eight new boundary events. AP is average precision. ECE is calculated using 10 equal-frequency groups ordered by predicted probability. RMSE and MAE are measured in percentage points. 1The calibration-in-the-large intercept uses logit(predicted probability) as an offset, with the slope fixed at one; its ideal value is zero. 2The calibration slope comes from a Logit diagnostic estimating both intercept and slope; its ideal value is one. Bounds of [10−6, 1–10−6] provide numerical protection when calculating logits; all probabilities in the three samples lie within these bounds. City-clustered confidence intervals for AUC and R2 appear in Appendix Table A5, and simple probability benchmarks appear in Appendix Table A3.

Information-set ablation and simple benchmarks

Information-set ablation shows that retaining yield history matters more for classification than removing the lagged yield anomaly alone. Excluding the lagged anomaly changes below-trend AUC in the major producing areas from 0.621 to 0.617. Removing all yield-history information, including expected yield, reduces AUC to 0.495. Using climate variables alone or capacity and production structure variables alone produces AUCs of approximately 0.477. Nationwide replication shows the same pattern. AUC is 0.618 with the full information set, 0.483 without yield history, 0.476 with climate alone, and 0.485 with capacity and structure alone (Figure 3b). Across these ablations in the major producing areas, continuous-outcome R2 ranges from −0.212 to −0.034, and severe-loss AUC ranges from 0.487 to 0.602 (Appendix Table A3).

Simple benchmarks further constrain the incremental value of the full model. A Logit model using only lagged yield and the lagged yield anomaly achieves a test AUC of 0.629 in the major producing areas, similar to 0.621 for the full 26-variable model. A constant prediction equal to development-sample prevalence has an AUC of 0.500. The full model’s Brier score of 0.261 also fails to improve on either 0.243 for the historical-information Logit or 0.248 for the development-prevalence benchmark (Appendix Table A3). Within the candidate models and samples examined here, production history provides the main identifiable predictive signal. Adding climate and capacity variables yields no consistent test improvement.

Adding expected yield to the historical-information Logit raises AUC in the major producing areas to 0.768, above 0.621 for the full 26-variable model. This comparison further highlights the information that the historical production benchmark provides for classifying deviations. Expected yield is generated from past information and also enters the definition of realized resistance. Both the three-variable model and the full model containing this variable therefore reflect dependence on the historical benchmark and outcome construction. We report the three-variable result as an audit of dependence on outcome construction, alongside the other simple benchmarks (Appendix Table A3).

Spatial transfer, nationwide replication, and model diagnostics

Strict external validation, nationwide replication, and probability calibration

Directly applying the model finalized in the major producing areas to external cities excluded from development produces a below-trend AUC of 0.613 across 237 events in 117 cities. Repeating the training and validation procedure for all 297 study cities nationwide yields an AUC of 0.618 across 595 test events in 278 cities. City-clustered bootstrap 95% confidence intervals are [0.529, 0.691] for strict external testing and [0.565, 0.668] for nationwide analysis. Together with the major-producing-area AUC of 0.621, these results indicate similar ranking performance over the full test period. Predicting the magnitude of continuous deviations remains more challenging (Table 7).

The major-producing-area and strict external samples use the same P80 threshold of 1.978655, while the nationwide workflow reestimates the threshold from its historical sample at 1.959824. The lower nationwide threshold adds eight boundary events, producing 350 + 237 + 8 = 595 test events and increasing city coverage from 277 in the union of the two groups to 278. Strict external results provide evidence of transfer to nonoverlapping regions, while nationwide results characterize overall performance after expanding the study area.

Annual results reveal temporal differences beneath this aggregate performance. AUCs for 2021–2023 are 0.775, 0.643, and 0.371 in the major producing areas, 0.714, 0.657, and 0.573 in the strict external sample, and 0.747, 0.627, and 0.469 nationwide. Appendix Table A5 reports annual sample sizes, prevalence, and AUCs, providing a more detailed basis for assessing how model performance varies over time.

Probability calibration further differentiates transfer performance across samples. The strict external sample has a Brier score of 0.305 and an ECE of 0.240, compared with 0.261 and 0.159 in the major producing areas. Its calibration slope is 0.165, compared with 0.592 in the major producing areas and 0.457 nationwide. A constant prediction fixed at the major-producing-area development prevalence of 0.5124 yields an external Brier score of 0.247, lower than the model’s 0.305 (Appendix Table A3). Transfer of ranking information and accurate estimation of probability levels therefore require separate evaluation.

Calibration-in-the-large diagnostics with the slope fixed at one yield intercepts of 0.723, 1.240, and 0.397 for the three samples, respectively. Corresponding mean predicted probabilities are 0.424, 0.383, and 0.509, whereas observed prevalence is 0.583, 0.616, and 0.595. Both sets of evidence indicate systematic underprediction of probabilities, with the largest discrepancy in the external sample. Figures 4b,c show the corresponding discrimination and calibration curves. These diagnostics use fixed predictions and provide a reference for subsequent probability calibration and threshold application.

Cross-sample comparisons of predictive information

Single-variable permutations of the below-trend classifiers’ fixed test predictions rank expected yield first in both the major-producing-area and nationwide analyses. Its permutation reduces AUC by an average of 0.185 and 0.175, respectively. Drought stress, high-temperature days, growing-season precipitation, rural electricity-use intensity, and the share of grain-sown area provide secondary information. The Spearman correlation between the importance rankings of all 26 variables in the two samples is 0.733 (p < 0.001; Figure 5), indicating similar patterns of predictive dependence.

Predictive information and cross-sample importance rankings. Panel a shows the mean AUC decreases for the top 10 variables, ordered by importance in the major-producing-area below-trend classifier, in the major-producing-area and nationwide tests. Error bars represent one standard deviation across repeated permutations, not parameter confidence intervals. Panel b compares importance ranks for all 26 variables across the two test samples. Rank one denotes the most important variable, and the dashed line indicates identical rankings. Dark blue identifies the top five variables in the major producing areas. Importance characterizes predictive dependence in fixed models, and the two test samples overlap.

This result complements information-set ablation and simple benchmarks by highlighting the role of the historical production benchmark in the current models. Because the nationwide test includes events from the major producing areas, ranking agreement here is a descriptive comparison after sample expansion. Lagged and expected yields are highly correlated, so shared information among variables also influences single-variable permutation results. Figure 5 is therefore interpreted as a diagnostic of predictive dependence in fixed models and read alongside the removal of yield history as a group.

Spatial residual structure and differences in data domains

After aggregating 2021–2023 test events by city, Global Moran’s I for prediction residuals is 0.266 in the major producing areas, 0.247 in the strict external sample, and 0.176 nationwide. All p-values from 499 permutations are 0.002, the smallest reportable value at this number of permutations. The residual tests identify regional prediction biases in all three samples. Meanwhile, Moran’s I for predicted below-trend probabilities is 0.239, 0.408, and 0.521, respectively, describing the spatial distribution of model-assessed risk (Figure 6). These diagnostics provide distinct spatial information about prediction errors and risk distributions.

Nationwide test residuals and cross-sample spatial correlation. Panel a maps mean residuals (observed minus predicted values) for the 278 cities in the nationwide test sample, in percentage points. Blue indicates realized resistance above predictions, and orange-red indicates realized resistance below predictions. The color scale is symmetrically capped at the 95th percentile of absolute residuals. Points beyond the limits remain visible in the endpoint colors. The complete standard map with approval number GS(2022)4308 is used as the basemap, with city residuals in a separate point layer. Panel b reports city-level Global Moran’s I for realized resistance, prediction residuals, and below-trend probabilities. Weights are constructed from the five nearest neighbors using longitude and latitude, excluding the city itself, symmetrizing adjacency, and standardizing rows. Significance is assessed with 499 permutations. Section 4.4.2 reports sensitivity to weights based on great-circle distances.

External testing also involves changes in data coverage and predictor ranges. For example, lagged rural electricity-use intensity is missing for 36.7% of the external sample, and agricultural machinery intensity for 32.9%. Grouping events by the number of missing variables gives an AUC of 0.638 for 179 events with no more than two missing variables and 0.600 for 58 events with more than two. Grouping by the number of variables outside the development-sample range gives an AUC of 0.640 for 137 events with no more than one such variable and 0.569 for 100 events with more than one. All groups use the same fixed model, and these descriptive comparisons characterize data-domain differences in the external sample. Together with the spatial diagnostics, they show that external transfer faces regional error structure, missing predictors, and changing predictor ranges. Aggregate AUC should therefore be assessed alongside data coverage and spatial diagnostics.

Robustness and sensitivity analyses

Sensitivity to P80 and P90 event thresholds

Raising the climate stress threshold from P80 to P90 reduces the training, validation, and test samples from 228, 174, and 350 events to 113, 143, and 290 events, respectively. For P90, we repeat candidate model comparisons, validation-based selection, classification-threshold fixing, and final fitting on the development sample. Random Forest is selected for the continuous task, Histogram Gradient Boosting for below-trend classification, and Random Forest for severe-loss classification. Table 3 reports the resulting test performance.

Below-trend AUC under P90 is 0.609, similar in magnitude to 0.621 under P80. AP values are 0.719 and 0.693, and Brier scores are 0.282 and 0.261, respectively. Continuous R2 is 0.024 under P90 and −0.081 under P80, both near zero. For severe-loss classification, the 32 positive P80 events yield an AUC of 0.597, with a city-clustered 95% confidence interval of [0.492, 0.702]. The 28 positive P90 events yield an AUC of 0.640, with an interval of [0.526, 0.764] (Table 3). Both event definitions produce below-trend discrimination of approximately 0.61. The severe-loss point estimates and intervals provide a basis for assessing the precision of ranking rare events.

Sensitivity to climate inputs, data processing, and loss definitions

Holding NASA POWER-identified P80 events and temporal splits unchanged, replacing the original climate variables with quantile features from the independent extreme climate dataset produces a continuous R2 of −0.048 and a below-trend AUC of 0.612. Combining both climate information sources yields −0.079 and 0.646, respectively. Severe-loss AUCs are 0.648 and 0.620. Climate input sources therefore affect some metrics, but the overall assessment of continuous predictive ability remains unchanged. This analysis compares alternative climate feature sources within the same P80 event sample (Appendix Table A4).

Removing anomalous yields above twice or below half the previous year’s value produces a continuous R2 of −0.079, a below-trend AUC of 0.627, and a severe-loss AUC that remains 0.597 to three decimal places. In a sensitivity sample using linear interpolation for missing raw agricultural outcomes, the test sample increases to 352 events, continuous R2 is 0.010, and below-trend and severe-loss AUCs are 0.598 and 0.589, respectively (Appendix Table A4). Linear interpolation provides a data-processing sensitivity comparison that complements the main analysis based on raw observations and past information.

The nationwide severe-loss analysis further distinguishes ranking and probability performance across deviation magnitudes. For R below −3%, −5%, and −10%, AUCs are 0.630, 0.747, and 0.573, with 135, 63, and 26 positive events, respectively. The city-clustered 95% confidence interval at −5% is [0.669, 0.817]. Model Brier scores at the three thresholds are 0.195, 0.278, and 0.143, compared with 0.177, 0.100, and 0.043 for constant predictions based on the corresponding development-sample prevalence (Appendix Table A4). The −5% outcome exhibits relatively strong risk ranking, while the simple benchmarks have lower mean squared probability error, again showing the need to evaluate event ranking and probability levels separately.

Sensitivity to spatial weights provides further evidence. Defining the five nearest neighbors by great-circle distance produces residual Moran’s I values of 0.222, 0.294, and 0.181 for the major-producing-area, external, and nationwide samples, respectively. All permutation p-values are 0.002, consistent with residual spatial clustering under the original weights.

Area weighting and one-year-lagged capacity scores both leave the rank association between capacity-based composite resilience and realized resistance weak. Across event thresholds, climate input sources, outlier treatment, missing-data treatment, loss definitions, and spatial weights, the city-level findings on continuous prediction, classification rankings, and spatial errors are broadly qualitatively consistent. Production history supplies most of the ranking information for below-trend classification, and aggregate discrimination remains similar in the strict external sample. Annual evaluations and probability benchmarks reveal temporal differences and calibration needs. We therefore characterize predictive performance through continuous deviation magnitudes, classification rankings, and probability accuracy, providing evidence for improving grain production risk assessment under climate stress.

Discussion

Accuracy depends on the target and validation task

The central finding is not that one algorithm outperforms another, but that the meaning of accuracy depends on what is predicted and where predictions must transfer. The provincial composite index is smooth, persistent, and partly defined by variables subsequently used as model inputs. Random row splitting places every test province in the training sample and positions test observations alongside nearby years from the same province. Thus, R2 = 0.957 accurately answers a narrow interpolation question. Can a flexible model reconstruct familiar capacity scores for familiar provinces and nearby years? It does not answer the policy-relevant question of whether these variables predict actual production failures under future shocks.

This distinction helps reconcile apparently different findings in the literature. Crane-Droesch (2018) demonstrates the value of flexible nonlinear models for crop yield prediction, including evaluation in held-out years. Hu et al. (2023), however, show why predictive accuracy alone does not ensure reliable interpretation or extrapolation when crop yield models are used to assess climate change impacts. More generally, validation should reflect the intended transfer task, and information leakage can distort performance estimates (Roberts et al., 2017; Stock et al., 2023). Evidence from shoreline mapping also indicates that high segmentation accuracy need not imply accurate shoreline positioning (Wang et al., 2026). Our study adds a measurement concern. Even a sophisticated learner cannot establish that a capacity-based target measures realized resistance merely by predicting it accurately.

The province-specific trend benchmark makes this point clearer. Under strict temporal holdout, simple provincial trends outperform the machine learning model. The composite index has structure, but that structure is dominated by persistent interprovincial differences and smooth temporal evolution. Once evaluation concerns future transfer, complexity provides no additional value. This is consistent with broader evidence that spatially random validation can yield optimistic environmental predictions (Meyer et al., 2018; Ploton et al., 2020).

Strict nonoverlapping external testing adds evidence on geographical transfer of the fixed model. Of the 127 external cities, 117 contribute 237 P80 events, yielding a below-trend AUC of 0.613, close to 0.621 in the major-producing-area temporal test. Continuous R2 is −0.152, and sensitivity at the fixed threshold is 0.247. The external model’s Brier score is 0.305, compared with 0.247 for a constant prediction fixed at major-producing-area development prevalence. Risk ranking and probability levels therefore require separate evaluation. Weaker performance among external observations with more missing or out-of-range inputs provides a concrete basis for improving data coverage and geographical calibration (Section 4.3.3 and Appendix Table A3).

Annual evaluation reveals temporal differences within the aggregate results. Below-trend AUCs in 2021–2023 are 0.775, 0.643, and 0.371 in the major producing areas, 0.714, 0.657, and 0.573 in the strict external sample, and 0.747, 0.627, and 0.469 nationwide (Appendix Table A5). Aggregate AUC thus characterizes overall ranking during this test period, while annual metrics provide complementary evidence on variation in risk identification over time.

Why realized resistance is difficult to predict

Outcome-based resistance requires predicting a city’s deviation from its own evolving benchmark under adverse climate exposure. Potential influences include seasonal temperature and precipitation, crop composition, sowing dates, varieties, irrigation, soil moisture, pests and diseases, local disasters, and management responses. Research in China documents adaptation through growing-season adjustments (Cui and Xie, 2022), changes in cropland allocation and crop composition (Cui and Zhong, 2024), adjustments in agricultural inputs such as labor, fertilizer, and machinery (Chen and Gong, 2021), and irrigation that mitigates yield losses from extreme heat (Wang et al., 2024). Field monitoring and modeling in arid northwestern China highlight the importance of soil texture and groundwater conditions for maize irrigation scheduling (Yi et al., 2022). Annual prefecture-level statistics capture these processes only incompletely and often with a one-year lag.

These factors offer possible explanations for continuous predictive performance. Aggregate prefecture-level grain yield combines crops with different phenologies and climate sensitivities. A common April–October window cannot precisely align exposure for rice, wheat, maize, and other crops, while city-centroid climate measures miss cropland-weighted spatial variation. Measurement error in outcomes and predictors may weaken stable signals, and unobserved adaptation may cause similar weather exposure to generate different outcomes across years.

Information-set ablation provides direct evidence about the sources of predictive information. Removing only the previous anomaly barely changes below-trend AUC, whereas removing all yield history reduces ranking to near chance in both the major-producing-area and nationwide samples. In the major producing areas, historical-information Logit achieves an AUC of 0.629, rising to 0.768 when expected yield is added, compared with 0.621 for the full 26-variable model. Both the three-variable and full models contain expected yield, which enters the outcome definition. These comparisons therefore jointly demonstrate the information value of the historical production benchmark for classifying deviations. Adding climate and capacity variables provides no consistent incremental gain in the current test samples (Appendix Table A3).

Construct validity comes before interpretable machine learning

Feature attribution describes how a model uses inputs under a given data distribution. It does not identify the effects of changing those inputs. Under random row splitting, attribution may largely reflect province identity, persistent production structures, and common trends. When strict-test AUC approaches 0.5 after outcome history is removed, even detailed feature-attribution rankings must be interpreted in light of predictive performance.

Machine learning studies of agricultural resilience should proceed in the following order. First, specify the system, function, shock, time horizon, and response dimension. Second, distinguish capacity indicators from realized outcomes and avoid treating components of the target as external explanations. Third, select temporal and spatial holdouts consistent with the intended application. Fourth, compare against simple trend, persistence, mean, or prior-probability benchmarks. Fifth, establish incremental predictive value and calibration. Feature contributions should be interpreted only after these conditions are met.

This sequence does not reject interpretable machine learning. It places interpretation after a successful predictive audit. Predictions should also be combined with knowledge of agronomic processes, crop calendars, and policy variation with independent identification. Machine learning can identify where stable signals exist, while causal or structural designs establish whether an input can change outcomes through policy.

The cross-sample ranking comparison in Figure 5 supports this interpretive sequence. The major-producing-area and nationwide analyses show similar information dependence, but the nationwide sample contains major-producing-area events. The comparison therefore describes agreement after coverage expands. Shared information between lagged and expected yields also affects single-variable permutation results. Considered alongside joint removal of yield history, continuous prediction, and probability calibration, importance rankings can be interpreted as predictive dependence in fixed models and used to identify agricultural processes requiring further study.

Implications for resilience measurement

Capacity indices and realized resistance answer different questions. The former asks whether a region possesses resources believed to support robustness or adaptation. The latter asks how the system performs under a defined disturbance relative to a reasonable counterfactual. Both are useful, but without validation they should not be treated as interchangeable manifestations of the same latent variable.

Monitoring systems can use capacity dashboards to track machinery, infrastructure, public resources, and production structures. Evaluation systems should use actual shock outcomes to test whether high-capacity regions maintain core functions. Repeated divergence between capacity scores and realized outcomes is not simply noise. It may reveal measurement gaps, critical bottlenecks, or context dependence in the conversion of resources into responses. A mature evidence system for resilience should retain separate modules for capacity, exposure, response, and recovery.

The rolling Huber benchmark supports this distinction by forming expectations solely from past information. It is transparent and substantially reduces the residual persistence generated by a fixed early-period trend. Its value lies in aligning the information set with the prediction task and making sensitivity to detrending explicit. The resulting yield deviations jointly reflect production changes and benchmark estimation. Studies interpreting yield anomalies as shocks should conduct comparable tests using benchmarks based on past information.

The province–year matching tests provide direct evidence for this distinction. Rank correlations between realized resistance and the capacity-based composite index and its components are all close to zero, and clustered bootstrap intervals do not exclude zero. This finding does not deny the theoretical importance of capacity reserves. It shows that capacity can be interpreted as a proxy for actual resilience only after external validation against clearly specified shock outcomes. Measurement systems should therefore report both capacity scores and outcome indicators, rather than substitute one index for another form of evidence.

Implications for policy and data systems

The results do not support directly ranking mechanization, research expenditure, fertilizer, or fiscal policy priorities by feature importance in models of the composite index. These variables may affect agricultural adaptation, but this study neither identifies their causal effects nor establishes that they independently predict ordinary below-trend outcomes. Policy prioritization requires research designs using credible exogenous variation or directly observed interventions.

For current production risk assessment and the subsequent development of early-warning systems, the priorities for data improvement are clearer. They include crop-specific sown area and yield, planting and harvest dates, irrigation coverage and water availability, soil moisture, remotely sensed vegetation conditions, actual disaster types and affected areas, agricultural insurance payouts, crop varieties, and weather aligned with crop development stages. Matching administrative statistics to cropland grids instead of city centroids would also improve spatial alignment.

Practical application can proceed gradually by region and outcome type. For study cities in the 13 major grain-producing provinces, model rankings can be combined with local agricultural conditions, crop composition, and disaster information to assess risk, with their supporting value evaluated through actual use. For the 127 external cities outside the major producing areas, local probability calibration and validation in later years can connect model outputs to local decision thresholds. For the nationwide severe-loss endpoint, the AUC of 0.747 motivates further research on ranking, while probability benchmarks and calibration results indicate a need to improve absolute risk estimates. In a different institutional setting, Zhou (2026) examines cognitive, institutional, technological, and social conditions associated with AI adoption in environmental governance. This offers a cross-sector perspective on implementation conditions; its applicability to agricultural risk systems requires separate evaluation.

Monitoring systems can also organize data modules by dominant exposure type. High-temperature risk areas could prioritize crop phenology, hot nights, and irrigation timing. Drought-prone areas could add root-zone soil moisture and available irrigation water, while areas exposed to excess precipitation could add inundation duration, drainage capacity, and disaster-affected area. These are priorities for data development corresponding to observed measurement gaps, not inferences about the causal effects of particular inputs.

The nationwide severe-loss results provide a direction for further investigation. At the −5% threshold, AUC is 0.747, with a city-clustered 95% confidence interval of [0.669, 0.817]. The Brier score is 0.278, compared with 0.100 for a constant benchmark using the corresponding development prevalence. The divergence between ranking performance and probability error indicates that further improvements should address both probability levels and decision thresholds, with prospective validation in later years and independent regions (Appendix Table A4).

Broader contributions to the literature

This study connects research on resilience measurement, yield benchmarks, and predictive validation. Farming system research distinguishes functions, capacities, processes, and outcomes (Meuwissen et al., 2019; van der Lee et al., 2022), while Zawalińska et al. (2022) link potential and revealed resilience to productivity in Polish agriculture during crises. Yield time-series studies show that detrending choices affect estimates of production anomalies and disaster losses (Lu et al., 2017; Meng and Qian, 2024). Production-based resilience measures can also be derived from crop production time series (Zampieri et al., 2020), although these measures differ from the event-specific resistance outcome used here. Research on structured data demonstrates the importance of validation design for estimates of performance and transferability (Roberts et al., 2017; Meyer et al., 2018; Ploton et al., 2020). Building on these contributions, we test construct alignment under climate stress in Chinese grain production and integrate benchmarks based on past information, information-set ablation, calibration, and external validation into a reproducible audit.

The study thus provides verifiable evidence on prediction and measurement. When the target is defined as realized resistance under explicit climate stress and evaluated using predefined temporal and spatial designs, existing annual prefecture-level statistics have limited ability to predict continuous loss magnitudes. Meanwhile, parsimonious historical benchmarks and benchmarks including expected yield show that classification AUC can exploit local yield persistence and information embedded in target construction. These findings clarify the information sources and improvement priorities for existing data and support combining capacity-based composite indices and city-level prediction models with direct validation against actual shock outcomes.

Robustness, limitations, and scope

We examine P80/P90 event thresholds, NASA and independent climate sources, major-producing-area and nationwide samples, information-set ablation, strict temporal holdouts, province holdouts, and the removal of anomalous jumps. Only three city–year observations in the raw major-producing-area sample have adjacent-year yield ratios above two or below 0.5, namely Chaoyang in Liaoning in 2009 and 2010, and Huludao in 2014. Removing these three outcome observations from both the targets and subsequent rolling histories changes continuous resistance R2 from −0.081 to −0.079 and below-trend AUC from 0.621 to 0.627, while severe-loss AUC remains 0.597 to three decimal places. The main conclusions are unchanged.

Severe-loss analysis adds information about precision across both event and loss thresholds. AUCs in the major producing areas are 0.597 under P80 (city-clustered 95% CI [0.492, 0.702]) and 0.640 under P90 ([0.526, 0.764]), based on 32 and 28 positive events, respectively. Nationwide AUC at the −5% threshold is 0.747 ([0.669, 0.817]). All three intervals use fixed test predictions and 2,000 city-level bootstrap resamples, with precision reported separately for each event definition. The nationwide analysis also includes existing audits using province and macroregional blocks, calibration diagnostics, and −3%, −5%, and −10% loss thresholds. These results reveal performance differences across outcome definitions and regional compositions.

Appendix Table A4 reports nationwide severe-loss ranking performance and probability errors at the −3%, −5%, and −10% thresholds, together with the confidence interval at −5%. This analysis complements the main outcome by documenting prediction differences across loss-magnitude definitions.

The linear-interpolation sensitivity analysis reruns the complete workflow. Test events increase from 350 to 352, continuous R2 is 0.010, below-trend AUC is 0.598, and severe-loss AUC is 0.589. Relative to the raw-data analysis, continuous prediction remains near zero and classification point estimates are slightly lower. This additional missing-data comparison appears in Appendix Table A4. The main analysis continues to use raw yield observations and a benchmark based on past information.

The study’s scope is limited in three main respects. First, outcome and exposure resolution is limited. Aggregate grain yields combine crops with different phenologies and climate sensitivities. Representative city points and a common April–October window provide a consistent basis for comparison, but finer data are needed on cropland distribution and winter and spring growing processes. Second, process observation and temporal coverage are limited. Annual city data for 2008–2023 support analysis of current production performance, whereas recovery, long-term adaptation, and transformation require continued observation and additional information on irrigation, soils, varieties, pests and diseases, crop calendars, and disasters. Third, transferability and sampling precision remain bounded. Final testing covers 3 years, and annual AUCs reveal temporal variation. City-clustered confidence intervals characterize city-sampling variation conditional on fixed predictions. Shared shocks across cities, changes across years, and uncertainty from model reestimation can be evaluated further with longer time series and independent spatial samples.

Data availability in external validation helps interpret transfer performance. Some input and employment variables have high missingness, while temperature, accumulated temperature, the share of grain-sown area, and extreme climate indicators fall outside development-sample ranges to varying degrees. Median imputation and missingness indicators allow the common pipeline to be applied externally. More complete, crop-specific, and consistently defined data would help distinguish the effects of missing measurements, covariate shifts, and changes in predictive relationships (Section 4.3.3).

Spatial clustering of residuals suggests regional structures that merit further investigation. Production structures, statistical definitions, topography and hydrology, and regional institutions may generate spatial association through observed or omitted variables. Future research could compare regional hierarchical models, spatial effects, and finer cropland-grid exposure measures under independent spatial blocking to identify how different information sources explain regional prediction biases.

Finally, this is a study of predictive and measurement validity. It does not identify causal effects of mechanization, research expenditure, fertilizer, irrigation, or fiscal policy. This boundary is important. Weak prediction does not imply causal irrelevance, and strong prediction does not imply policy actionability. Future work should combine prediction that prevents leakage, crop-specific process data, and causal designs focused on specific adaptation investments.

Conclusion

This study reassesses whether machine learning can identify actual grain production resilience by examining measurement validity and external transferability. The capacity-based composite index achieves R2 = 0.957 under random splitting, but the strict temporal model does not outperform a simple province-specific trend benchmark, and leave-one-province-out R2 is 0.531. More importantly, direct matching of 120 province–year observations from 13 major grain-producing provinces yields a Spearman correlation of only 0.022 between the composite capacity score and realized resistance under climate shocks. Clustered bootstrap intervals for all three components also cross zero. High predictability of the composite index therefore mainly reflects indicator composition, regional persistence, and temporal smoothness, and cannot be taken as evidence that it validly measures actual shock responses.

When the target is changed to realized resistance constructed using only past information, the full 26-variable model provides limited predictive performance. In the major producing areas, continuous R2 is −0.081 and below-trend AUC is 0.621. Applied to the 117 cities with P80 events among the 127 external cities, the fixed model yields an AUC of 0.613 and continuous R2 of −0.152. Nationwide replication across 278 test cities yields an AUC of 0.618 and continuous R2 of −0.120. Similar aggregate AUCs coexist with annual variation, and calibration slopes are below one in all three samples. The external model has a Brier score of 0.305, compared with 0.247 for the corresponding constant development-prevalence benchmark. Comparisons between parsimonious historical Logit and the full model, together with the AUC improvement after adding expected yield, highlight the role of historical benchmarks and outcome-construction information in classification.

Spatial diagnostics reveal regional features of model performance. Global Moran’s I for prediction residuals is 0.266 in the major producing areas, 0.247 in the strict external sample, and 0.176 in the nationwide replication, with p = 0.002 in all 499-permutation tests. The spatial distribution of predicted probabilities and regional residual biases provide different forms of information. Topography and hydrology, crop composition, statistical definitions, and regional conditions are potential explanations for future investigation. Alternative climate inputs, event thresholds, removal of anomalous jumps, and interpolation sensitivity complement the assessment of continuous prediction, risk ranking, and probability accuracy. Severe-loss intervals add evidence about precision for rare outcomes.

By connecting capacity measurement, shock outcomes, and external validation in a predictive audit framework, this study provides direct evidence for assessing the empirical meaning of composite indices and the scope of city-level model applications. Capacity