## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.align = "centre",
  warning   = FALSE,
  message   = FALSE
)
library(agriReg)

## ----load-data----------------------------------------------------------------
wheat <- load_example_data("wheat_trial")
head(wheat)

## ----clean--------------------------------------------------------------------
wheat_clean <- clean_agri_data(wheat,
                                yield_col     = "yield",
                                flag_outliers = TRUE)
summary(wheat_clean)

## ----normalise----------------------------------------------------------------
wheat_norm <- yield_normalize(wheat_clean,
                               yield_col = "yield",
                               method    = "zscore",
                               group_by  = "block")
head(wheat_norm[, c("block", "yield", "yield_norm")])

## ----outlier------------------------------------------------------------------
wheat_clean <- outlier_flag(wheat_clean, col = "yield", method = "modz")
table(wheat_clean$yield_outlier)

## ----lm-simple----------------------------------------------------------------
m_lm <- fit_linear(wheat_clean, "yield ~ nitrogen + phosphorus")
print(m_lm)

## ----lm-plots, fig.height=4---------------------------------------------------
plot(m_lm, type = "fit")
plot(m_lm, type = "residuals")
plot(m_lm, type = "qq")

## ----residual-panel, fig.height=6---------------------------------------------
residual_check(m_lm)

## ----lmer---------------------------------------------------------------------
m_lmer <- fit_linear(wheat_clean,
                     formula = "yield ~ nitrogen + phosphorus + variety",
                     random  = "(1|block)")
summary(m_lmer)

## ----poly---------------------------------------------------------------------
m_poly <- fit_polynomial(wheat_clean,
                          x_col      = "nitrogen",
                          y_col      = "yield",
                          max_degree = 3)
coef(m_poly)

## ----maize-load---------------------------------------------------------------
maize <- load_example_data("maize_growth")
# Use well-watered plants only
maize_ww <- maize[maize$treatment == "WW", ]
head(maize_ww)

## ----logistic-----------------------------------------------------------------
m_log <- fit_nonlinear(maize_ww,
                        x_col = "days",
                        y_col = "biomass_g",
                        model = "logistic")
summary(m_log)

## ----logistic-plot, fig.height=4----------------------------------------------
plot(m_log)

## ----gompertz, eval = FALSE---------------------------------------------------
# m_gomp <- fit_nonlinear(maize_ww, "days", "biomass_g", "gompertz")
# coef(m_gomp)

## ----asymptotic---------------------------------------------------------------
m_asym <- fit_nonlinear(wheat_clean,
                         x_col = "nitrogen",
                         y_col = "yield",
                         model = "asymptotic")
plot(m_asym)

## ----quadratic----------------------------------------------------------------
m_quad <- fit_nonlinear(wheat_clean,
                         x_col = "nitrogen",
                         y_col = "yield",
                         model = "quadratic")

opt <- optimum_dose(m_quad)
cat(sprintf("Optimum nitrogen rate : %.1f kg/ha\n", opt["x_opt"]))
cat(sprintf("Predicted max yield   : %.2f t/ha\n",  opt["y_max"]))

## ----lp-----------------------------------------------------------------------
m_lp <- fit_nonlinear(wheat_clean,
                       x_col = "nitrogen",
                       y_col = "yield",
                       model = "linear_plateau")
cat("Critical point (plateau onset):", optimum_dose(m_lp)["x_opt"], "kg/ha\n")

## ----herb-load----------------------------------------------------------------
herb <- load_example_data("herbicide_trial")
# Subset to one species
amaranth <- herb[herb$species == "Amaranth", ]

## ----drc----------------------------------------------------------------------
m_dr <- fit_dose_response(amaranth,
                           dose_col = "dose_g_ha",
                           resp_col = "weed_control_pct",
                           fct      = "LL.4")
summary(m_dr)

## ----drc-plot, fig.height=4---------------------------------------------------
plot(m_dr, log_dose = TRUE)

## ----ed-----------------------------------------------------------------------
ed_estimates(m_dr, respLev = c(10, 50, 90))

## ----compare------------------------------------------------------------------
cmp <- compare_models(
  ols        = m_lm,
  mixed      = m_lmer,
  asymptotic = m_asym,
  quadratic  = m_quad,
  lp         = m_lp
)
print(cmp)

## ----broom--------------------------------------------------------------------
library(broom)
tidy(m_lm$fit)
glance(m_lm$fit)

## ----save-plot, eval=FALSE----------------------------------------------------
# p <- plot(m_log)
# ggplot2::ggsave("logistic_growth.png", p, width = 7, height = 4.5, dpi = 300)

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