test_that("Bayesian mixed model fits when brms is available", {
  
  skip_if_not_installed("brms")
  skip_on_cran()
  
  set.seed(123)
  
  dat <- data.frame(
    id = factor(rep(1:10, each = 3)),
    x = rep(1:3, 10)
  )
  
  dat$y <- 5 + 2 * dat$x +
    rnorm(30) +
    rnorm(10)[dat$id]
  
  m <- suppressWarnings(
    fit_bayes_mix(
      y ~ x + (1 | id),
      data = dat,
      chains = 1,
      iter = 100,
      warmup = 50,
      refresh = 0
    )
  )
  
  expect_s3_class(m, "biomix_bayes")
  expect_s3_class(m, "biomix")
  expect_true(inherits(get_fit(m), "brmsfit"))
})


test_that("GAMM fits when mgcv is available", {
  
  skip_if_not_installed("mgcv")
  
  set.seed(123)
  
  dat <- data.frame(
    id = factor(rep(1:10, each = 5)),
    x = rep(seq(0, 1, length.out = 5), 10)
  )
  
  dat$y <- sin(2 * pi * dat$x) +
    rnorm(50, sd = 0.2)
  
  m <- fit_gamm(
    y ~ s(x, k = 4),
    data = dat,
    random = list(id = ~1)
  )
  
  expect_s3_class(m, "biomix_gamm")
  expect_true(is.list(m))
  expect_true(is.finite(AIC(get_fit(m)$lme)))
})


test_that("multivariate mixed model function validates input", {
  
  expect_error(
    fit_multivariate(
      list(y1 = 1:5, y2 = 1:5),
      data.frame(x = 1:5)
    ),
    regexp = "formula|data|model",
    ignore.case = TRUE
  )
})


test_that("spatial mixed model validates input", {
  
  expect_error(
    fit_spatial_mix(
      y ~ x,
      data = data.frame(
        y = rnorm(10),
        x = rnorm(10)
      )
    ),
    regexp = "spatial|coordinate|random",
    ignore.case = TRUE
  )
})


test_that("survival mixed model fits when coxme is available", {
  
  skip_if_not_installed("coxme")
  
  set.seed(123)
  
  dat <- data.frame(
    time = rexp(50),
    status = rbinom(50, 1, 0.7),
    x = rnorm(50),
    id = factor(rep(1:10, each = 5))
  )
  
  m <- fit_survival_mix(
    survival::Surv(time, status) ~ x + (1 | id),
    data = dat
  )
  
  expect_s3_class(m, "biomix_survival")
  expect_s3_class(m, "biomix")
  expect_true(is.list(m))
})


test_that("nonlinear mixed model fits", {
  
  skip_if_not_installed("nlme")
  
  set.seed(123)
  
  dat <- data.frame(
    id = factor(rep(1:10, each = 5)),
    x = rep(0:4, 10)
  )
  
  dat$y <- 20 * (1 - exp(-0.5 * dat$x)) +
    rnorm(50, 0, 0.5)
  
  m <- fit_nlme_biomix(
    formula = y ~ SSasymp(x, Asym, R0, lrc),
    data = dat,
    fixed = Asym + R0 + lrc ~ 1,
    random = Asym ~ 1 | id,
    start = c(
      Asym = 20,
      R0 = 0,
      lrc = log(0.5)
    ),
    control = nlme::nlmeControl(
      maxIter = 100,
      msMaxIter = 200
    )
  )
  
  expect_s3_class(m, "biomix_nlmm")
  expect_s3_class(m, "biomix")
  expect_true(is.finite(AIC(m$fit)))
})


test_that("accessor returns underlying model", {
  
  skip_if_not_installed("nlme")
  
  set.seed(123)
  
  dat <- data.frame(
    id = factor(rep(1:5, each = 4)),
    x = rep(1:4, 5)
  )
  
  dat$y <- 5 + 2 * dat$x +
    rnorm(20) +
    rnorm(5)[dat$id]
  
  m <- fit_biomix(
    y ~ x,
    data = dat,
    random = ~1 | id
  )
  
  fit <- get_fit(m)
  
  expect_true(inherits(fit, "lme"))
})


test_that("nonlinear methods are available", {
  
  skip_if_not_installed("nlme")
  
  set.seed(123)
  
  dat <- data.frame(
    id = factor(rep(1:5, each = 4)),
    x = rep(0:3, 5)
  )
  
  dat$y <- 10 * (1 - exp(-0.5 * dat$x)) +
    rnorm(20, 0, 0.2)
  
  m <- fit_nlme_biomix(
    formula = y ~ SSasymp(x, Asym, R0, lrc),
    data = dat,
    fixed = Asym + R0 + lrc ~ 1,
    random = Asym ~ 1 | id,
    start = c(
      Asym = 10,
      R0 = 0,
      lrc = log(0.5)
    ),
    control = nlme::nlmeControl(
      maxIter = 100,
      msMaxIter = 200
    )
  )
  
  expect_s3_class(m, "biomix_nlmm")
  expect_true(is.finite(AIC(m$fit)))
  expect_true(is.finite(BIC(m$fit)))
})

test_that("model_summary works for nonlinear mixed-effects model", {
  
  skip_if_not_installed("nlme")
  
  set.seed(123)
  
  dat <- data.frame(
    id = factor(rep(1:5, each = 5)),
    x = rep(0:4, 5)
  )
  
  dat$y <- 20 * (1 - exp(-0.5 * dat$x)) +
    rnorm(25, 0, 0.3)
  
  m <- fit_nlme_biomix(
    formula = y ~ SSasymp(x, Asym, R0, lrc),
    data = dat,
    fixed = Asym + R0 + lrc ~ 1,
    random = Asym ~ 1 | id,
    start = c(
      Asym = 20,
      R0 = 0,
      lrc = log(0.5)
    ),
    control = nlme::nlmeControl(
      maxIter = 100,
      msMaxIter = 200
    )
  )
  
  s <- model_summary(m)
  
  expect_type(s, "list")
  expect_true(all(c(
    "model_type",
    "n",
    "AIC",
    "BIC",
    "logLik",
    "fixed_effects",
    "random_effects"
  ) %in% names(s)))
  
  expect_true(is.character(s$model_type))
  expect_true(is.numeric(s$n))
  expect_true(is.finite(s$AIC))
  expect_true(is.finite(s$BIC))
  expect_true(is.finite(s$logLik))
  expect_true(is.numeric(s$fixed_effects))
  expect_true(is.data.frame(s$random_effects))
})

test_that("model_summary.biomix_nlmm rejects non-nlmm model", {
  
  m <- fit_biomix(
    y ~ x,
    data = data.frame(
      y = rnorm(10),
      x = rnorm(10)
    )
  )
  
  expect_error(
    model_summary.biomix_nlmm(m)
  )
})

# ---------------------------------------------------------
# Bayesian model validation
# ---------------------------------------------------------

test_that("Bayesian model validates formula", {
  
  skip_if_not_installed("brms")
  dat <- data.frame(
    y = rnorm(10),
    x = rnorm(10)
  )
  
  expect_error(
    fit_bayes_mix(
      formula = "y ~ x",
      data = dat
    ),
    "formula must be a formula"
  )
})


test_that("Bayesian model validates data", {
  
  skip_if_not_installed("brms")
  expect_error(
    fit_bayes_mix(
      formula = y ~ x,
      data = list(
        y = rnorm(10),
        x = rnorm(10)
      )
    ),
    "data must be a data frame"
  )
})


# ---------------------------------------------------------
# Bayesian family validation
# ---------------------------------------------------------

test_that("Bayesian model rejects unsupported family", {
  
  skip_if_not_installed("brms")
  dat <- data.frame(
    y = rnorm(10),
    x = rnorm(10)
  )
  
  expect_error(
    fit_bayes_mix(
      y ~ x,
      data = dat,
      family = "unsupported_family"
    ),
    "Unsupported Bayesian family"
  )
})


test_that("Bayesian model rejects multiple character families", {
  
  skip_if_not_installed("brms")
  dat <- data.frame(
    y = rnorm(10),
    x = rnorm(10)
  )
  
  expect_error(
    fit_bayes_mix(
      y ~ x,
      data = dat,
      family = c("gaussian", "poisson")
    ),
    "Unsupported Bayesian family"
  )
})


test_that("Bayesian model rejects invalid family object", {
  
  skip_if_not_installed("brms")
  dat <- data.frame(
    y = rnorm(10),
    x = rnorm(10)
  )
  
  expect_error(
    fit_bayes_mix(
      y ~ x,
      data = dat,
      family = list(name = "gaussian")
    ),
    "'family' must be a supported character name or a brms family object"
  )
})


# ---------------------------------------------------------
# Bayesian brms family object
# ---------------------------------------------------------

test_that("Bayesian model accepts a brms family object", {
  
  skip_if_not_installed("brms")
  skip_on_cran()
  set.seed(456)
  
  dat <- data.frame(
    y = rnorm(20),
    x = rnorm(20)
  )
  
  m <- suppressWarnings(
    fit_bayes_mix(
      y ~ x,
      data = dat,
      family = brms::brmsfamily("gaussian"),
      chains = 1,
      iter = 100,
      warmup = 50,
      refresh = 0
    )
  )
  
  expect_s3_class(m, "biomix_bayes")
  expect_s3_class(m, "biomix")
  expect_true(inherits(m$fit, "brmsfit"))
  expect_true(inherits(m$family, "brmsfamily"))
})


# ---------------------------------------------------------
# Bayesian model stores supplied components
# ---------------------------------------------------------

test_that("Bayesian model stores supplied components", {
  
  skip_if_not_installed("brms")
  skip_on_cran()
  
  set.seed(789)
  
  dat <- data.frame(
    y = rnorm(20),
    x = rnorm(20)
  )
  
  m <- suppressWarnings(
    fit_bayes_mix(
      y ~ x,
      data = dat,
      family = "gaussian",
      chains = 1,
      iter = 100,
      warmup = 50,
      refresh = 0
    )
  )
  
  expect_equal(
    m$formula,
    y ~ x
  )
  
  expect_identical(
    m$data,
    dat
  )
  
  expect_equal(
    m$family,
    "gaussian"
  )
  
  expect_equal(
    m$model_type,
    "Bayesian mixed model"
  )
})


# ---------------------------------------------------------
# Supported Bayesian families
# ---------------------------------------------------------

test_that("Bayesian model accepts supported character families", {
  
  skip_if_not_installed("brms")
  
  # Test the validation branch without fitting a Stan model.
  # brms::brm() is mocked so this test remains fast.
  
  local_mocked_bindings(
    brm = function(...) {
      structure(
        list(),
        class = "brmsfit"
      )
    },
    .package = "brms"
  )
  
  dat <- data.frame(
    y = rnorm(10),
    x = rnorm(10)
  )
  
  families <- c(
    "gaussian",
    "bernoulli",
    "binomial",
    "poisson",
    "negbinomial",
    "Gamma",
    "student",
    "lognormal"
  )
  
  for (fam in families) {
    
    m <- fit_bayes_mix(
      y ~ x,
      data = dat,
      family = fam
    )
    
    expect_s3_class(m, "biomix_bayes")
    expect_equal(m$family, fam)
  }
})