## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 7,
  fig.height = 4,
  message   = FALSE,
  warning   = FALSE
)

## ----install, eval = FALSE----------------------------------------------------
# # From CRAN (once published)
# install.packages("NRMstatsML")
# 
# # Development version from GitHub
# # install.packages("remotes")
# remotes::install_github("yourorg/NRMstatsML")

## ----load-data----------------------------------------------------------------
library(NRMstatsML)

data(nrm_example)
head(nrm_example)

## ----data-check---------------------------------------------------------------
nrm_data_check(nrm_example)

## ----trend--------------------------------------------------------------------
trend_res <- nrm_trend(nrm_example,
                       time_var  = "year",
                       value_var = "crop_yield",
                       breaks    = TRUE)
print(trend_res)

## ----trend-plot---------------------------------------------------------------
nrm_plot(trend_res)

## ----mk-only------------------------------------------------------------------
mk <- nrm_mann_kendall(nrm_example, time_var = "year",
                       value_var = "soil_OC")
print(mk)

ss <- nrm_sens_slope(nrm_example, time_var = "year",
                     value_var = "soil_OC")
print(ss)

## ----multivariate-------------------------------------------------------------
mv <- nrm_multivariate(nrm_example,
        formula = crop_yield ~ N + P + K + rainfall,
        scale   = TRUE)
print(mv)

## ----pls----------------------------------------------------------------------
pl <- nrm_pls(nrm_example,
              formula = crop_yield ~ N + P + K + rainfall + soil_OC,
              ncomp   = 3)
print(pl)

## ----response-curve-----------------------------------------------------------
rc <- nrm_response_curve(nrm_example,
        input_var    = "N",
        response_var = "crop_yield",
        type         = "quadratic")
print(rc)
nrm_plot(rc)

## ----optimise-----------------------------------------------------------------
opt <- nrm_optimize_input(rc, price_ratio = 0.6)
print(opt)

## ----arima--------------------------------------------------------------------
ar <- nrm_arima(nrm_example, value_var = "crop_yield")
print(ar)

## ----forecast-----------------------------------------------------------------
fc <- nrm_forecast(nrm_example,
                   value_var = "crop_yield",
                   horizon   = 5)
print(fc)
nrm_plot(fc)

## ----panel, eval = FALSE------------------------------------------------------
# # 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)

## ----bootstrap----------------------------------------------------------------
bs <- nrm_bootstrap(nrm_example,
        stat_fn = function(d) mean(d$crop_yield),
        n_iter  = 500)
print(bs)

## ----monte-carlo--------------------------------------------------------------
mc <- nrm_monte_carlo(nrm_example,
        stat_fn   = function(d) mean(d$crop_yield),
        n_iter    = 500,
        noise_sd  = 0.1)
print(mc)

## ----uncertainty--------------------------------------------------------------
unc <- nrm_uncertainty(nrm_example,
         stat_fn = function(d) mean(d$crop_yield),
         method  = "bootstrap",
         n_iter  = 500)
print(unc)

## ----automl, eval = FALSE-----------------------------------------------------
# # 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)

## ----benchmark, eval = FALSE--------------------------------------------------
# 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)

## ----session-info-------------------------------------------------------------
sessionInfo()