Seasonal effects.
Explicit assumptions.

Eight complementary estimators

Shared estimation

YieldLag models annual regional yield y_it using a standardized time trend, monthly weather terms, and regional effects. Preprocessing and coefficients are estimated within each training window. Ridge penalties control coefficient size; second-difference penalties smooth consecutive seasonal-month coefficients. Robust fitting uses iteratively reweighted Student-t-style residual weights. Optional pooled AR(1) correction uses residuals preceding a prediction year.

Baselines are independent regional trends and last-observation persistence. Panel ridge pools weather slopes; local ridge estimates regional slopes; PCR estimates training-window principal components; the hierarchical estimator partially pools regional deviations around a global seasonal response.

Seasonal stress-lag basis

For variable v, seasonal month m, region i, and training set T:

mu_i,v,m = mean_T(x_i,t,v,m)
s_v,m^2 = sum_T((x_i,t,v,m - mu_i,v,m)^2) / (n_T - number_of_regions)
z_i,t,v,m = (x_i,t,v,m - mu_i,v,m) / s_v,m
L(z) = z
H(z) = max(z - tau, 0)
D(z) = max(-z - tau, 0)
C(z_T, z_P) = H(z_T) * D(z_P)

Near-zero scales are replaced by one. For an unseen region, pooled training means replace regional means; regional effects are zero. This fallback enables computation but requires separate validation of spatial transfer.

Every declared weather variable receives linear, upper-tail, and lower-tail terms. Compound terms are optional and require explicit temperature and precipitation variable mappings. They match calendar month and seasonal position. A temperature/precipitation pair contributes one product per matched month. Tau is positive and expressed in units of within-region weather variability. It is not an absolute temperature or physiological threshold.

Expanded basis columns are centered and scaled within the training set. The penalized model is

y_it = intercept + b * standardized_year
       + regional_intercept_i + regional_trend_i * standardized_year
       + sum_j(theta_j * standardized_basis_it,j) + error_it

The panel-ridge solver estimates the coefficients. Each variable/basis component has its own seasonal second-difference penalty; compound products are smoothed as a separate component. A common ridge strength and smoothing ratio control all weather components. Regional intercept and trend penalties are 0.05 and 0.25 times the weather ridge penalty; the global trend penalty is 0.05 times that penalty. Lambda, smoothing ratio, and tau are tuned on earlier rolling forecasts. Climatology and basis scaling are refitted inside every forecast-origin training window.

The basis captures nonlinear associations in monthly climate. It does not reconstruct daily degree-days, physiological stress, phenology, or irrigation. The sign of a coefficient is estimated; damage signs are not imposed.

Prediction interpretation

explain_crop_prediction() returns contributions for intercept, trend, regional effects, linear anomalies, upper/lower tails, hot-dry interaction, and AR correction. Their sum reproduces predict() numerically. Centering means a tail contribution can be nonzero even when its raw hinge is zero; the reference is the training basis mean. These are additive statistical components, not causal attribution or measured crop losses.

prediction_support() reports the number of weather terms outside pooled training minima/maxima, the maximum absolute standardized anomaly, and unseen region status. Marginal ranges do not diagnose every unsupported multivariate combination and do not provide a coverage guarantee.

Ensemble and uncertainty

The ensemble fits non-negative weights summing to one. Shrinkage toward equal weighting is tuned by leaving out complete tuning years. Base hyperparameters are themselves chosen using the tuning block, so this is not a fully nested estimate of model-selection uncertainty. Weights are frozen for later evaluation.

Prediction interval half-widths are empirical tuning-error quantiles. National half-widths are separately calibrated from national tuning errors. The nominal 80% label describes the calibration target rather than a formal coverage guarantee.

References

• Conradt (2022), ABSOLUT regional weather-aggregate selection and regression: https://doi.org/10.1007/s00484-022-02356-5. YieldLag provides a different, regularized panel framework; the evaluation adapter runs upstream GPL code separately and does not incorporate that implementation into the MIT package.
• Schlenker and Roberts (2009), nonlinear weather-yield associations: Read the study.
• Gasparrini, Armstrong, and Kenward (2010), distributed lag nonlinear models: Read the study.

YieldLag implements penalized monthly hinge/product features motivated by these established nonlinear-response and distributed-lag ideas.

Validation and interpretation

Chronological design

Observed years are divided into initial history, tuning, and later evaluation. A prediction for year t uses model coefficients and preprocessing fitted to observations with years less than t. Hyperparameters are selected from rolling tuning scores and remain fixed in later evaluation; coefficients and preprocessing are refitted at each later forecast origin. Automatic selection minimizes province-year tuning RMSE and does not use evaluation rankings.

The stress-lag climatology, anomaly scale, expanded feature scaling, and regional effects obey the same training-window rule. Compound mappings are declared by the user. The ensemble combines requested base estimators using tuning data.

By default, tuning errors are reused for hyperparameter/method selection and interval calibration. Configure calibration_years > 0 to reserve years after tuning and before evaluation for separate interval calibration. Calibration responses never select method parameters or the automatic method. These temporally dependent empirical bands carry no distribution-free guarantee. The later block estimates the selected procedure on this dataset and period. Once used to guide development, it should be described as a development benchmark.

Diagnostics

RMSE, MAE, and bias retain the response unit. Bias is prediction minus observation. Predictive R-squared can be negative; squared correlation does not assess calibration. Coverage and width describe empirical prediction bands. Province-year observations share years and regions and are correlated. Year/block resampling is appropriate when estimating uncertainty for comparisons. compare_crop_errors() resamples paired consecutive circular year blocks, keeping all regions in each sampled year. It preserves the caller's random state. Its interval describes RMSE differences on the supplied period; seven annual blocks remain a small basis for inference.

Prediction support flags marginal weather extrapolation and unseen regions. Its absence does not establish support of every weather combination. Decompositions reproduce the fitted equation, not causal or physiological effects.

Poland example

With a 2018 training cutoff, the default periods are 1999--2005 initial history, 2006--2011 tuning, and 2012--2018 evaluation. Weather through 2019 permits a retrospective 2019 forecast. National evaluation has seven annual observations. The example uses fixed area weights and area-weighted provincial observed yields, rather than independently observed national yields or origin-specific annual areas.

Synthetic nonlinear example

synthetic_stress_example() supplies a fully declared deterministic equation with asymmetric temperature/precipitation responses and a compound hot-dry term. A later-period comparison checks that the estimator can recover this kind of nonlinear signal. This controlled result does not demonstrate superiority on observed crops or independently validate a new scientific method.

External evaluation

A scientific comparison should include untouched later data, spatial holdouts, persistence and trend baselines, uncertainty for paired metric differences, sensitivity to seasons/grids, and information available at the intended issue date. Failures and weak crop targets should be reported alongside gains. Operational forecasting requires weather observations or forecasts actually available at that date and appropriately time-varying aggregation weights.