# tests/testthat/test-linear_models.R
# ---- shared fixture -------------------------------------------------------
make_wheat <- function(n = 60, seed = 42) {
  set.seed(seed)
  data.frame(
    yield      = 3 + 0.03 * rep(c(0,60,120,180), each = n/4) +
      rnorm(n, 0, 0.3),
    nitrogen   = rep(c(0, 60, 120, 180), each = n/4),
    phosphorus = rep(c("low","high"), times = n/2),
    block      = rep(c("A","B","C"), length.out = n),
    stringsAsFactors = FALSE
  )
}
dat <- make_wheat()
# ---- fit_linear -----------------------------------------------------------
test_that("fit_linear returns agriLM with lm engine", {
  m <- fit_linear(dat, "yield ~ nitrogen")
  expect_s3_class(m, "agriLM")
  expect_equal(m$engine, "lm")
})
test_that("fit_linear accepts formula object", {
  m <- fit_linear(dat, "yield ~ nitrogen + phosphorus")
  expect_s3_class(m, "agriLM")
})
test_that("fit_linear lm coefficients are named", {
  m <- fit_linear(dat, "yield ~ nitrogen")
  cf <- coef(m)
  expect_true("nitrogen" %in% names(cf))
  expect_true("(Intercept)" %in% names(cf))
})
test_that("fit_linear fitted values have correct length", {
  m <- fit_linear(dat, "yield ~ nitrogen")
  expect_length(fitted(m), nrow(dat))
})
test_that("fit_linear residuals sum to ~0", {
  m <- fit_linear(dat, "yield ~ nitrogen")
  expect_lt(abs(sum(residuals(m))), 1e-8)
})
test_that("fit_linear predict works on new data", {
  m  <- fit_linear(dat, "yield ~ nitrogen")
  nd <- data.frame(nitrogen = c(90, 150))
  pr <- predict(m, newdata = nd)
  expect_length(pr, 2)
})
test_that("fit_linear with random effects uses lmer engine", {
  skip_if_not_installed("lme4")
  m <- fit_linear(dat, "yield ~ nitrogen", random = "(1|block)")
  expect_equal(m$engine, "lmer")
  expect_s3_class(m, "agriLM")
})
test_that("coef.agriLM returns fixef for lmer", {
  skip_if_not_installed("lme4")
  m  <- fit_linear(dat, "yield ~ nitrogen", random = "(1|block)")
  cf <- coef(m)
  expect_named(cf)
  expect_true("nitrogen" %in% names(cf))
})
# ---- print / summary ------------------------------------------------------
test_that("print.agriLM outputs agriReg header", {
  m <- fit_linear(dat, "yield ~ nitrogen")
  expect_output(print(m), "agriReg")
})
test_that("summary.agriLM shows goodness-of-fit section", {
  m <- fit_linear(dat, "yield ~ nitrogen")
  expect_output(summary(m), "Goodness-of-fit")
})
# ---- plot -----------------------------------------------------------------
test_that("plot.agriLM type='fit' returns ggplot", {
  skip_if_not_installed("ggplot2")
  m <- fit_linear(dat, "yield ~ nitrogen")
  p <- plot(m, type = "fit")
  expect_s3_class(p, "gg")
})
test_that("plot.agriLM type='residuals' returns ggplot", {
  skip_if_not_installed("ggplot2")
  m <- fit_linear(dat, "yield ~ nitrogen")
  p <- plot(m, type = "residuals")
  expect_s3_class(p, "gg")
})
test_that("plot.agriLM type='qq' returns ggplot", {
  skip_if_not_installed("ggplot2")
  m <- fit_linear(dat, "yield ~ nitrogen")
  p <- plot(m, type = "qq")
  expect_s3_class(p, "gg")
})
test_that("plot.agriLM type='scale' returns ggplot", {
  skip_if_not_installed("ggplot2")
  m <- fit_linear(dat, "yield ~ nitrogen")
  p <- plot(m, type = "scale")
  expect_s3_class(p, "gg")
})
# ---- fit_polynomial -------------------------------------------------------
test_that("fit_polynomial returns agriLM for best degree", {
  m <- fit_polynomial(dat, x_col = "nitrogen", y_col = "yield", max_degree = 2)
  expect_s3_class(m, "agriLM")
})
test_that("fit_polynomial prints comparison table", {
  expect_output(
    fit_polynomial(dat, "nitrogen", "yield", max_degree = 2),
    "degree"
  )
})