Introduction
NRMstatsML is an R package for statistical and
machine learning-based analysis of long-term Natural Resource Management
(NRM) datasets. It is designed for datasets spanning 10–20 years
covering soil health, water, crop yield, and climate variables.
The package is organised into seven core modules:
trendML |
nrm_trend, nrm_mann_kendall,
nrm_sens_slope, nrm_structural_break |
Trend detection and structural change |
multiSysML |
nrm_multivariate, nrm_pls,
nrm_sem |
Multivariate system modelling |
responseML |
nrm_response_curve,
nrm_optimize_input |
Yield-response and optimisation |
tsML |
nrm_forecast, nrm_arima |
Time series forecasting |
panelML |
nrm_panel, nrm_did |
Panel data and treatment effects |
uncertaintyML |
nrm_uncertainty, nrm_bootstrap,
nrm_monte_carlo |
Uncertainty quantification |
autoML |
nrm_automl, nrm_benchmark |
Automated model selection |
Example Data
The package ships with nrm_example, a synthetic 20-year
dataset that mimics a long-term fertiliser experiment.
library(NRMstatsML)
data(nrm_example)
head(nrm_example)
#> year treatment crop_yield soil_OC N P K rainfall runoff biomass
#> 1 2000 control 3.37 0.867 55.6 23.2 26.3 782.5 76.0 5.32
#> 2 2001 NPK 3.30 0.836 63.1 22.8 34.6 676.4 75.1 4.70
#> 3 2002 control 3.21 0.882 68.8 23.2 35.7 712.6 60.4 4.98
#> 4 2003 NPK 3.26 0.829 79.3 23.5 29.5 657.5 77.0 6.10
#> 5 2004 control 3.75 0.870 66.6 27.5 36.8 574.1 69.3 5.85
#> 6 2005 NPK 3.87 0.783 76.4 25.0 41.2 672.5 80.4 5.89
Always start by validating the data:
nrm_data_check(nrm_example)
#> === NRMstatsML: Data Check ===
#> Rows : 20
#> Columns : 10
#> Numeric : year, crop_yield, soil_OC, N, P, K, rainfall, runoff, biomass
#> Status : OK
Module 3 — Response Curve & Optimisation (responseML)
rc <- nrm_response_curve(nrm_example,
input_var = "N",
response_var = "crop_yield",
type = "quadratic")
print(rc)
#> === NRMstatsML: Response Curve ===
#> Type : quadratic
#> Input : N
#> Response : crop_yield
#> R-squared : 0.3446
nrm_plot(rc)
Economic optimum
Find the nitrogen rate that maximises profit given a cost-to-price
ratio:
opt <- nrm_optimize_input(rc, price_ratio = 0.6)
print(opt)
#> === NRMstatsML: Input Optimisation ===
#> Price ratio : 0.600
#> Biophysical optimum : NA
#> Economic optimum : NA
Module 4 — Time Series Forecasting (tsML)
ar <- nrm_arima(nrm_example, value_var = "crop_yield")
print(ar)
#> === NRMstatsML: ARIMA Model ===
#> Variable : crop_yield
#> Order : ARIMA(0,1,0)
#> AIC = 23.98 BIC = 24.93
fc <- nrm_forecast(nrm_example,
value_var = "crop_yield",
horizon = 5)
print(fc)
#> === NRMstatsML: Time Series Forecast ===
#> Variable : crop_yield
#> Method : auto_arima
#> Horizon : 5 steps
#> Accuracy (in-sample):
#> ME RMSE MAE MPE MAPE MASE ACF1
#> Training set 0.0317 0.4206 0.2887 0.2445 7.5499 0.9506 -0.3701
nrm_plot(fc)
Module 5 — Panel Data & Treatment Effects (panelML)
Panel and DiD analyses require a multi-site, multi-year dataset. See
?nrm_panel and ?nrm_did for worked examples
with appropriate data.
# Requires a panel dataset with site and year identifiers
pm <- nrm_panel(panel_data,
formula = crop_yield ~ N + P + K + rainfall,
index = c("site", "year"),
model = "within")
nrm_summary(pm)
did <- nrm_did(panel_data,
outcome = "crop_yield",
treat_var = "treatment_binary",
time_var = "post_period")
print(did)
Module 6 — Uncertainty & Sensitivity Analysis
(uncertaintyML)
bs <- nrm_bootstrap(nrm_example,
stat_fn = function(d) mean(d$crop_yield),
n_iter = 500)
print(bs)
#> -- Bootstrap --
#> Mean : 3.7513 SD : 0.0916 CI : [3.5735, 3.9325]
mc <- nrm_monte_carlo(nrm_example,
stat_fn = function(d) mean(d$crop_yield),
n_iter = 500,
noise_sd = 0.1)
print(mc)
#> -- Monte Carlo --
#> Mean : 3.7508 SD : 0.0094 CI : [3.7325, 3.7685]
Use nrm_uncertainty() as a unified entry point:
unc <- nrm_uncertainty(nrm_example,
stat_fn = function(d) mean(d$crop_yield),
method = "bootstrap",
n_iter = 500)
print(unc)
#> === NRMstatsML: Uncertainty Analysis ===
#> Method : bootstrap
#> Mean : 3.7558
#> SD : 0.0904
#> 95% CI : [3.5866, 3.9416]
Module 7 — AutoML & Model Benchmarking (autoML)
# Requires caret + method-specific packages (e.g. randomForest, gbm)
am <- nrm_automl(nrm_example,
formula = crop_yield ~ N + P + K + rainfall + soil_OC,
methods = c("lm", "rf", "gbm"),
cv_folds = 5,
seed = 42)
nrm_summary(am)
Benchmarking on a hold-out set
n <- nrow(nrm_example)
train <- nrm_example[seq_len(floor(0.8 * n)), ]
test <- nrm_example[seq(floor(0.8 * n) + 1, n), ]
m_ols <- lm(crop_yield ~ N + P + K, data = train)
bm <- nrm_benchmark(
models = list(ols = m_ols),
test_data = test,
response_var = "crop_yield"
)
print(bm)
Recommended Workflow
nrm_data_check() # 1. Validate data
↓
nrm_trend() # 2. Detect trends and breaks
↓
nrm_multivariate() / nrm_pls() # 3. Multivariate modelling
↓
nrm_response_curve() # 4. Fit yield-response
nrm_optimize_input() # → Identify optimal inputs
↓
nrm_forecast() # 5. Forecast future values
↓
nrm_uncertainty() # 6. Quantify uncertainty
↓
nrm_automl() # 7. Compare candidate models
nrm_benchmark() # → Select best for deployment
