Install
The model engine uses R's standard libraries. R 4.2 or later is required.
install.packages("remotes")
remotes::install_github("MR-Eini/yieldlag", ref = "v0.3.0",
build_vignettes = FALSE)
library(yieldlag)
Alternatively, download yieldlag_0.3.0.tar.gz from the release and use install.packages(path, repos = NULL, type = "source"). The release archive includes rendered tutorials. The package name is yieldlag.
Prepare data
The directory interface reads the bundled study layout:
barley <- read_crop_data(poland_example_path(), "barley",
harvest_month = 6)
wheat <- read_crop_data(poland_example_path(), "wheat",
harvest_month = 7)
The harvest month is the endpoint of a 12-month agricultural window, not a physiological phenology model. Barley uses July–June; wheat uses August–July. Seven variables × twelve months produce 84 terms.
For your own data, use prepare_crop_data() with one row per region/year, finite weather terms, and yield values or NA for prediction years. Map names explicitly rather than adopting the Poland identifiers.
panel <- data.frame(
province = c("A", "A", "B", "B"),
year = c(2020, 2021, 2020, 2021),
yield = c(40, 42, 38, 41),
temperature_june = c(18, 20, 17, 19),
rainfall_june = c(65, 40, 72, 48)
)
manifest <- data.frame(
Term = c("temperature_june", "rainfall_june"),
Variable = c("TMP", "PCP"),
CalendarMonth = c(6, 6), AgriculturalPosition = c(1, 1)
)
data <- prepare_crop_data(panel, manifest,
region = "province", year = "year", yield = "yield")
This small table illustrates the schema only; a chronological comparison requires at least twelve observed years. Optional positive area weights and independent national observations can also be supplied. Missing area weights produce equal-weight aggregation and this is recorded.
Compare methods
config <- crop_model_config(
methods = crop_model_methods(),
stress_temperature = "MAX",
stress_precipitation = "PCP",
robust = TRUE
)
comparison <- compare_crop_models(barley, config,
last_yield_year = 2018, predict_through = 2019)
summary(comparison)
plot(comparison)
The method is selected by minimum province-year tuning RMSE. Later evaluation scores are reported separately. In the website runs, initial history is 1999–2005, tuning is 2006–2011, and evaluation is 2012–2018. At each origin, coefficients and preprocessing are refitted using earlier yields; evaluation-period hyperparameters remain fixed.
last_yield_year controls the final fitting cutoff. Existing 2019 yields remain available as observations for plots, but do not enter fitting or tuning. The 2019 predictions use complete agricultural-year weather, so they are retrospective.
Choose a subset of methods for a faster comparison. Supply a named grids list to change the tuning candidates; document those choices before interpreting results.
Predict and diagnose
selected <- get_crop_model(comparison)
future <- subset(barley$panel, Year == 2019)
predict(selected, future, interval = "prediction")
prediction_support(selected, future)
stress <- get_crop_model(comparison, "stress_lag")
explain_crop_prediction(stress, future)
coef(stress)
Prediction intervals are empirical 80% bands calibrated from tuning errors. They are available after comparison; direct fit_crop_model() does not independently calibrate them. National bands use national tuning errors rather than averaging provincial interval widths.
prediction_support() counts weather terms outside pooled training ranges and identifies unseen regions. explain_crop_prediction() supports stress-lag models and returns components in the response unit that sum to the prediction. Components are centred relative to training features; they are not causal attribution.
Fit directly
Use fixed parameters when tuning is not required:
fictional <- synthetic_stress_example()
model <- fit_crop_model(fictional, "stress_lag",
parameters = list(lambda = 1, smooth_ratio = 0,
threshold = 0.75, temperature_variable = "TMP",
precipitation_variable = "PCP"), robust = FALSE)
Direct fitting uses all finite responses in the supplied data. Mask or remove future yields first when a forecast-origin cutoff is required. The deterministic example has a declared nonlinear equation and provides a recovery test, not real-crop validation.
Export and reproduce
write_crop_results(comparison, "results/barley")
saveRDS(comparison, "results/barley/comparison.rds")
Exports include model metrics, row-level rolling and final predictions, candidate tuning scores, selected settings, coefficients, ensemble weights, prediction components, support diagnostics, configuration, fitted objects, plots, session information, and checksums. Existing artifacts require explicit overwrite = TRUE.
The website calculations can be reproduced with run_poland_website.R: install YieldLag 0.2.0 and jsonlite, then run the script from a working directory with Rscript run_poland_website.R barley or Rscript run_poland_website.R wheat. The script uses the published package defaults and enables MAX/PCP compound terms. Download the result tables to inspect the exact configuration and settings.
Interpret results
RMSE, MAE and bias retain the yield unit; bias is prediction minus observation. Predictive R² is 1 − SSE/SST and can be negative. Squared correlation measures association, not agreement. National scores and provincial scores describe different targets.
A flexible estimator need not beat a simpler model. The Poland runs show crop-dependent rankings, and tuning selection can differ from the best later-period result. Provincial errors are correlated by region and year; 112 province-years do not represent 112 independent experiments. See validation assumptions and data terms.