R/DryingKineticModels.R

In dryingkineticmodels: Drying Kinetic Models Comparison and Analysis

#' @importFrom grDevices dev.off png
#' @importFrom graphics abline hist legend lines par plot
#' @importFrom stats AIC BIC coef deviance df.residual fitted lm median pf predict qt residuals shapiro.test qqline qqnorm
#' @importFrom utils browseURL
"_PACKAGE"

# ============================================================
# TWO-STAGE GRID SEARCH HELPER FUNCTION
# ============================================================
#' Two-Stage Grid Search for NLS Fitting
#'
#' Performs a coarse grid search over a parameter grid and returns the
#' best-fitting nonlinear least squares model using \code{nlsLM}.
#'
#' @param formula A nonlinear model formula.
#' @param data A data frame containing the variables in \code{formula}.
#' @param param_grid A named list of numeric vectors defining the grid of
#'   starting values for each parameter.
#' @param control Control parameters passed to
#'   \code{\link[minpack.lm]{nls.lm.control}}.
#'
#' @return The best-fitting \code{nls} object, or \code{NULL} if all
#'   starting value combinations failed to converge.
#'
#' @keywords internal
grid_search_fit_2stage <- function(formula, data, param_grid,
                                   control = minpack.lm::nls.lm.control(maxiter = 200)) {
  coarse_grid <- expand.grid(param_grid)
  best_fit    <- NULL
  best_rss    <- Inf

  for (i in seq_len(nrow(coarse_grid))) {
    start_vals <- as.list(coarse_grid[i, , drop = FALSE])
    tryCatch({
      fit <- minpack.lm::nlsLM(formula, data = data, start = start_vals, control = control)
      rss <- sum(residuals(fit)^2)
      if (rss < best_rss) {
        best_rss <- rss
        best_fit <- fit
      }
    }, error = function(e) NULL)
  }
  best_fit
}

# ============================================================
# MODEL FITTING FUNCTIONS
# ============================================================

fit_lewis <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ exp(-k * time), data = dat,
      param_grid = list(k = c(0.001, 0.01, 0.1, 0.5)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_page <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ exp(-k * time^n), data = dat,
      param_grid = list(
        k = c(0.001, 0.01, 0.1),
        n = c(0.5, 1.0, 1.5, 2.0)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_modified_page <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ exp(-(k * time)^n), data = dat,
      param_grid = list(
        k = c(0.001, 0.01, 0.05),
        n = c(0.5, 1.0, 1.5, 2.0)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_hp <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time), data = dat,
      param_grid = list(
        a = c(0.9, 1.0, 1.1),
        k = c(0.001, 0.01, 0.1)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_logarithmic <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time) + c, data = dat,
      param_grid = list(
        a = c(0.9, 1.0, 1.1),
        k = c(0.001, 0.01, 0.1),
        c = c(-0.1, 0.0, 0.1)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_two_term <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k0 * time) + b * exp(-k1 * time), data = dat,
      param_grid = list(
        a  = c(0.3, 0.6, 0.9),
        b  = c(0.1, 0.4, 0.7),
        k0 = c(0.001, 0.01, 0.1),
        k1 = c(0.001, 0.01, 0.1)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_two_term_exp <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time) + (1 - a) * exp(-(k * a * time)),
      data = dat,
      param_grid = list(
        a = c(0.3, 0.5, 0.7),
        k = c(0.001, 0.01, 0.1)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_diffusion_approx <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time) + (1 - a) * exp(-(k * b * time)),
      data = dat,
      param_grid = list(
        a = c(0.3, 0.6, 0.9),
        b = c(0.2, 0.5, 0.8),
        k = c(0.001, 0.01, 0.1)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_wang_singh <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ 1 + a * time + b * time^2, data = dat,
      param_grid = list(
        a = c(-0.1, -0.01, -0.001),
        b = c(0.0001, 0.001, 0.01)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_midilli <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time^n) + b * time, data = dat,
      param_grid = list(
        a = c(0.9, 1.0, 1.1),
        k = c(0.001, 0.01, 0.1),
        n = c(0.8, 1.2, 1.8),
        b = c(-0.001, -0.0001, 0.0001)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_modified_hp <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time) + b * exp(-g * time) + c * exp(-h * time),
      data = dat,
      param_grid = list(
        a = c(0.4, 0.7),
        b = c(0.2, 0.4),
        c = c(0.1, 0.2),
        k = c(0.01, 0.1),
        g = c(0.001, 0.01),
        h = c(0.0001, 0.001)),
      control = minpack.lm::nls.lm.control(maxiter = 200))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_verma <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time) + (1 - a) * exp(-g * time), data = dat,
      param_grid = list(
        a = c(0.3, 0.6, 0.9),
        k = c(0.001, 0.01, 0.1),
        g = c(0.001, 0.005, 0.05)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_weibull <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ exp(-(time / alpha)^beta), data = dat,
      param_grid = list(
        alpha = c(10, 50, 100, 200),
        beta  = c(0.5, 1.0, 1.5, 2.0)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_aghbashlo <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ exp(-k1 * time / (1 + k2 * time)), data = dat,
      param_grid = list(
        k1 = c(0.001, 0.01, 0.1),
        k2 = c(0.0001, 0.001, 0.01)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_jena_das <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time^n) + b * time + c, data = dat,
      param_grid = list(
        a = c(0.9, 1.0, 1.1),
        k = c(0.001, 0.01),
        n = c(0.8, 1.2, 1.8),
        b = c(-0.001, -0.0001),
        c = c(-0.05, 0.0, 0.05)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_hii <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time^n) + b * exp(-g * time^m), data = dat,
      param_grid = list(
        a = c(0.6, 0.9),
        k = c(0.001, 0.01),
        n = c(0.8, 1.2),
        b = c(0.1, 0.3),
        g = c(0.001, 0.005),
        m = c(0.5, 1.0)),
      control = minpack.lm::nls.lm.control(maxiter = 200))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_parabolic <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a + b * time + c * time^2, data = dat,
      param_grid = list(
        a = c(0.9, 1.0, 1.1),
        b = c(-0.1, -0.01, -0.001),
        c = c(0.00001, 0.0001, 0.001)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_thompson <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ exp((-a + sqrt(a^2 + 4 * b * time)) / (2 * b)), data = dat,
      param_grid = list(
        a = c(-50, -20, -10, -5),
        b = c(0.5, 1.0, 2.0, 5.0)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_demir <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time^n) + b, data = dat,
      param_grid = list(
        a = c(0.9, 1.0, 1.1),
        k = c(0.001, 0.01, 0.1),
        n = c(0.8, 1.2, 1.8),
        b = c(-0.05, 0.0, 0.05)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

fit_thin_layer_explin <- function(dat) {
  tryCatch({
    fit <- grid_search_fit_2stage(
      MR ~ a * exp(-k * time) + b * time + c, data = dat,
      param_grid = list(
        a = c(0.7, 0.9, 1.0),
        k = c(0.001, 0.01, 0.1),
        b = c(-0.001, -0.0001),
        c = c(-0.05, 0.0, 0.05)))
    if (is.null(fit)) stop("Grid search failed to converge")
    list(success = TRUE, fit = fit, summary = summary(fit),
         coefficients = coef(fit), predicted = predict(fit),
         residuals = residuals(fit))
  }, error = function(e) list(success = FALSE, message = e$message))
}

# ============================================================
# RUN ALL MODELS
# ============================================================
#' Fit All Twenty Drying Kinetic Models
#'
#' Calls all individual model-fitting functions and returns results as a
#' named list.
#'
#' @param dat A data frame with columns \code{time} and \code{MR}.
#' @param models Optional character vector of model names to fit. If
#'   \code{NULL} (default), all 20 models are fitted, preserving the
#'   original behavior. Valid names are \code{"Lewis"}, \code{"Page"},
#'   \code{"Modified_Page"}, \code{"Henderson_Pabis"},
#'   \code{"Logarithmic"}, \code{"Two_Term"},
#'   \code{"Two_Term_Exponential"}, \code{"Diffusion_Approximation"},
#'   \code{"Wang_Singh"}, \code{"Midilli_Kucuk"},
#'   \code{"Modified_Henderson_Pabis"}, \code{"Verma"}, \code{"Weibull"},
#'   \code{"Aghbashlo"}, \code{"Jena_Das"}, \code{"Hii_et_al"},
#'   \code{"Parabolic"}, \code{"Thompson"}, \code{"Demir_et_al"}, and
#'   \code{"Thin_Layer_ExpLin"}.
#'
#' @return A named list. Each element is the list returned by the
#'   corresponding \code{fit_*()} function. Length 20 if \code{models}
#'   is \code{NULL}, otherwise \code{length(models)}.
#'
#' @keywords internal
fit_all_models <- function(dat, models = NULL) {
  fit_functions <- list(
    Lewis                    = fit_lewis,
    Page                     = fit_page,
    Modified_Page            = fit_modified_page,
    Henderson_Pabis          = fit_hp,
    Logarithmic              = fit_logarithmic,
    Two_Term                 = fit_two_term,
    Two_Term_Exponential     = fit_two_term_exp,
    Diffusion_Approximation  = fit_diffusion_approx,
    Wang_Singh               = fit_wang_singh,
    Midilli_Kucuk            = fit_midilli,
    Modified_Henderson_Pabis = fit_modified_hp,
    Verma                    = fit_verma,
    Weibull                  = fit_weibull,
    Aghbashlo                = fit_aghbashlo,
    Jena_Das                 = fit_jena_das,
    Hii_et_al                = fit_hii,
    Parabolic                = fit_parabolic,
    Thompson                 = fit_thompson,
    Demir_et_al              = fit_demir,
    Thin_Layer_ExpLin        = fit_thin_layer_explin
  )

  if (!is.null(models)) {
    unknown <- setdiff(models, names(fit_functions))
    if (length(unknown) > 0) {
      stop("Unknown model name(s): ", paste(unknown, collapse = ", "),
           "\nValid names are: ", paste(names(fit_functions), collapse = ", "))
    }
    fit_functions <- fit_functions[models]
  }

  lapply(fit_functions, function(f) f(dat))
}

# ============================================================
# COMPARISON TABLE
# ============================================================
#' Compare Fitted Drying Kinetic Models
#'
#' Computes goodness-of-fit statistics for all successfully fitted models
#' and returns a ranked comparison table.
#'
#' @param results A named list of model results as returned by
#'   \code{fit_all_models()}.
#'
#' @return A data frame with one row per successfully fitted model
#'   (models with \eqn{R^2 \le 0} are excluded), sorted by a composite
#'   rank score based on \eqn{R^2}, RMSE, and \eqn{\chi^2}. Columns are:
#'   \describe{
#'     \item{Kinetic Drying Models}{Model name.}
#'     \item{Functional Form}{Equation of the model.}
#'     \item{R2}{Coefficient of determination.}
#'     \item{RMSE}{Root mean square error.}
#'     \item{MAE}{Mean absolute error.}
#'     \item{chi2}{Reduced chi-squared statistic.}
#'     \item{RSS}{Residual sum of squares.}
#'     \item{Estimated Coefficients}{Parameter estimates as a character string.}
#'   }
#'
#' @keywords internal
compare_models <- function(results) {

  functional_forms <- c(
    Lewis                    = "MR = exp(-k*t)",
    Page                     = "MR = exp(-k*(t^n))",
    Modified_Page            = "MR = exp(-(k*t)^n)",
    Henderson_Pabis          = "MR = a*exp(-k*t)",
    Logarithmic              = "MR = a*exp(-k*t) + c",
    Two_Term                 = "MR = a*exp(-k0*t) + b*exp(-k1*t)",
    Two_Term_Exponential     = "MR = a*exp(-k*t) + (1-a)*exp(-(k*a*t))",
    Diffusion_Approximation  = "MR = a*exp(-k*t) + (1-a)*exp(-(k*b*t))",
    Wang_Singh               = "MR = 1 + a*t + b*(t^2)",
    Midilli_Kucuk            = "MR = a*exp(-k*(t^n)) + b*t",
    Modified_Henderson_Pabis = "MR = a*exp(-k*t) + b*exp(-g*t) + c*exp(-h*t)",
    Verma                    = "MR = a*exp(-k*t) + (1-a)*exp(-g*t)",
    Weibull                  = "MR = exp(-(t/alpha)^beta)",
    Aghbashlo                = "MR = exp((-k1*t) / (1 + k2*t))",
    Jena_Das                 = "MR = a*exp(-k*(t^n)) + b*t + c",
    Hii_et_al                = "MR = a*exp(-k*(t^n)) + b*exp(-g*(t^m))",
    Parabolic                = "MR = a + b*t + c*(t^2)",
    Thompson                 = "MR = exp((-a + sqrt(a^2 + 4b*t)) / 2b)",
    Demir_et_al              = "MR = a*exp(-k*(t^n)) + b",
    Thin_Layer_ExpLin        = "MR = a*exp(-k*t) + b*t + c"
  )

  rows <- lapply(names(results), function(model_name) {
    m <- results[[model_name]]
    if (!m$success) {
      return(data.frame(
        `Kinetic Drying Models`  = model_name,
        `Functional Form`        = functional_forms[model_name],
        R2 = NA, RMSE = NA, MAE = NA, chi2 = NA, RSS = NA,
        `Estimated Coefficients` = NA,
        stringsAsFactors = FALSE, check.names = FALSE))
    }
    obs    <- m$predicted + m$residuals
    res    <- m$residuals
    n      <- length(obs)
    p      <- length(m$coefficients)
    ss_res <- deviance(m$fit)
    ss_tot <- sum((obs - mean(obs))^2)
    r2_val <- 1 - ss_res / ss_tot

    coef_str <- paste(
      names(m$coefficients),
      round(m$coefficients, 4),
      sep = " = ", collapse = ", ")

    data.frame(
      `Kinetic Drying Models`  = model_name,
      `Functional Form`        = functional_forms[model_name],
      R2   = round(r2_val,                                    6),
      RMSE = round(sqrt(ss_res / (n - p)),                   6),
      MAE  = round(mean(abs(res)),                            6),
      chi2 = round((ss_res / (n - p)) / (ss_tot / (n - 1)), 6),
      RSS  = round(ss_res,                                    6),
      `Estimated Coefficients` = coef_str,
      stringsAsFactors = FALSE, check.names = FALSE)
  })

  tbl <- do.call(rbind, rows)
  tbl <- tbl[!is.na(tbl$R2) & tbl$R2 > 0, ]
  tbl$rank_r2   <- rank(-tbl$R2,   ties.method = "min")
  tbl$rank_rmse <- rank( tbl$RMSE, ties.method = "min")
  tbl$rank_chi2 <- rank( tbl$chi2, ties.method = "min")
  tbl$score     <- tbl$rank_r2 + tbl$rank_rmse + tbl$rank_chi2
  tbl <- tbl[order(tbl$score), ]
  tbl$rank_r2 <- tbl$rank_rmse <- tbl$rank_chi2 <- tbl$score <- NULL
  rownames(tbl) <- NULL
  tbl
}
#' ANOVA Table for the Best NLS Model
#'
#' Computes and prints a regression ANOVA table for a fitted nonlinear
#' least squares model.
#'
#' @param best_model A model result list (the \code{success = TRUE} list
#'   returned by any \code{fit_*()} function).
#' @param dat A data frame with columns \code{time} and \code{MR}.
#'
#' @return A data frame containing the ANOVA table with columns
#'   \code{Source}, \code{Sum of Squares}, \code{df},
#'   \code{Mean Square}, \code{F Value}, and \code{Pr(>F)}.
#'   The table is also printed to the console.
#'
#' @keywords internal
anova_nls <- function(best_model, dat) {
  fit    <- best_model$fit
  obs    <- dat$MR
  pred   <- fitted(fit)
  res    <- obs - pred
  n      <- length(obs)
  p      <- length(coef(fit))
  ss_res <- deviance(fit)
  df_res <- df.residual(fit)
  ss_tot <- sum((obs - mean(obs))^2)
  ss_reg <- ss_tot - ss_res
  df_reg <- p
  df_tot <- n - 1
  ms_reg <- ss_reg / df_reg
  ms_res <- ss_res / df_res
  f_val  <- ms_reg / ms_res
  p_val  <- pf(f_val, df_reg, df_res, lower.tail = FALSE)

  anova_tbl <- data.frame(
    Source           = c("Regression", "Residual", "Total"),
    `Sum of Squares` = round(c(ss_reg, ss_res, ss_tot), 6),
    df               = c(df_reg, df_res, df_tot),
    `Mean Square`    = c(round(ms_reg, 6), round(ms_res, 6), NA),
    `F Value`        = c(round(f_val, 4), NA, NA),
    `Pr(>F)`         = c(format.pval(p_val, digits = 4, eps = 0.0001), "", ""),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  return(anova_tbl)
}

# ============================================================
# DIAGNOSTIC INTERPRETATION
# ============================================================
#' Interpret Residual Diagnostic Tests
#'
#' Runs four residual diagnostic tests on the best model and prints
#' a plain-language interpretation of each result.
#'
#' @param best_model A model result list (the \code{success = TRUE} list
#'   returned by any \code{fit_*()} function).
#' @param verbose Logical. If \code{TRUE}, prints the interpretation to
#'   the console via \code{message()}. Default \code{FALSE}.
#'
#' @return A list (returned invisibly) with elements \code{shapiro},
#'   \code{dw}, \code{bp}, and \code{runs} — the raw test objects from
#'   \code{\link[stats]{shapiro.test}},
#'   \code{\link[lmtest]{dwtest}},
#'   \code{\link[lmtest]{bptest}}, and
#'   \code{\link[tseries]{runs.test}} respectively.
#'
#' @keywords internal
diagnostic_interpretation <- function(best_model, verbose = FALSE) {
  res    <- best_model$residuals
  fitted <- best_model$predicted
  shapiro_res <- shapiro.test(res)
  dw_res      <- lmtest::dwtest(lm(res ~ fitted))
  bp_res      <- lmtest::bptest(lm(res ~ fitted))
  runs_res    <- tseries::runs.test(as.factor(ifelse(res >= median(res), "+", "-")))
  if (verbose) {
    message(paste(
      "====================================",
      "MODEL DIAGNOSTIC INTERPRETATION",
      "====================================\n",
      "1. NORMALITY CHECK",
      paste("Shapiro-Wilk p-value :", round(shapiro_res$p.value, 4)),
      paste("Interpretation :", ifelse(shapiro_res$p.value > 0.05,
                                       "Residuals are approximately normal.",
                                       "Residuals deviate from normality."), "\n"),
      "2. AUTOCORRELATION CHECK",
      paste("Durbin-Watson statistic :", round(dw_res$statistic, 4)),
      paste("Durbin-Watson p-value :",   round(dw_res$p.value,   4)),
      paste("Interpretation :", ifelse(dw_res$p.value > 0.05,
                                       "Residuals appear independent.",
                                       "Residuals show autocorrelation (dependency)."), "\n"),
      "3. HOMOSCEDASTICITY CHECK",
      paste("Breusch-Pagan p-value :", round(bp_res$p.value, 4)),
      paste("Interpretation :", ifelse(bp_res$p.value > 0.05,
                                       "Homoscedasticity assumption satisfied (constant variance).",
                                       "Heteroscedasticity detected (non-constant variance)."), "\n"),
      "4. RUNS TEST",
      paste("Runs Test statistic :", round(runs_res$statistic, 4)),
      paste("Runs Test p-value   :", round(runs_res$p.value,   4)),
      paste("Interpretation :", ifelse(runs_res$p.value > 0.05,
                                       "Residuals are randomly distributed (no pattern).",
                                       "Non-random pattern detected in residuals.")),
      sep = "\n"))
  }
  invisible(list(shapiro = shapiro_res, dw = dw_res,
                 bp = bp_res, runs = runs_res))
}

# ============================================================
# RESIDUAL DIAGNOSTICS PLOTS
# ============================================================
#' Plot Residual Diagnostics
#'
#' Produces a 2x2 panel of residual diagnostic plots: residuals vs
#' fitted values, histogram of residuals, normal Q-Q plot, and a
#' residual sequence plot.
#'
#' @param best_model A model result list (the \code{success = TRUE} list
#'   returned by any \code{fit_*()} function).
#'
#' @return Called for its side effect (plots). Returns \code{NULL}
#'   invisibly.
#'
#' @keywords internal
residual_diagnostics <- function(best_model) {
  res    <- best_model$residuals
  fitted <- best_model$predicted
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  par(mfrow = c(2, 2))
  plot(fitted, res, main = "Residuals vs Fitted",
       xlab = "Fitted Values", ylab = "Residuals")
  abline(h = 0, lty = 2)
  hist(res, main = "Histogram of Residuals", xlab = "Residuals")
  qqnorm(res); qqline(res)
  plot(res, type = "b", main = "Residual Sequence Plot",
       xlab = "Observation", ylab = "Residuals")
}

# ============================================================
# BEST MODEL
# ============================================================
#' Extract and Summarise the Best Model
#'
#' Identifies the top-ranked model from the comparison table, prints its
#' coefficients, fit statistics, hypothesis test results, and runs
#' residual diagnostics.
#'
#' @param results A named list of model results as returned by
#'   \code{fit_all_models()}.
#' @param comparison A data frame as returned by \code{compare_models()}.
#' @param dat A data frame with columns \code{time} and \code{MR}.
#' @param verbose Logical. If \code{TRUE}, prints coefficients, fit
#'   statistics, and diagnostic interpretation to the console via
#'   \code{message()}. Default \code{FALSE}.
#'
#' @return The model result list for the best model (the
#'   \code{success = TRUE} list returned by the corresponding
#'   \code{fit_*()} function).
#'
#' @keywords internal
get_best_model <- function(results, comparison, dat, verbose = FALSE) {
  best_name  <- comparison$`Kinetic Drying Models`[1]
  best_model <- results[[best_name]]
  obs <- dat$MR
  res    <- best_model$residuals
  n      <- length(res)
  p      <- length(best_model$coefficients)
  ss_res <- sum(res^2)
  ss_tot <- sum((obs - mean(obs))^2)
  r2   <- round(1 - ss_res / ss_tot, 6)
  rmse <- round(sqrt(ss_res / (n - p)), 6)
  if (verbose) {
    message(paste(
      "====================================",
      paste("BEST MODEL :", best_name),
      "====================================\n",
      "Coefficients:",
      paste("R2  :", r2),
      paste("RMSE:", rmse),
      paste("AIC :", round(AIC(best_model$fit), 4)),
      paste("BIC :", round(BIC(best_model$fit), 4)),
      "====================================",
      "HYPOTHESIS TESTING OF PARAMETERS",
      "====================================",
      sep = "\n"))
  }
  hyp_raw <- coef(summary(best_model$fit))
  hyp_df  <- data.frame(
    Parameter      = rownames(hyp_raw),
    Estimate       = round(hyp_raw[, "Estimate"],   6),
    Std_Error      = round(hyp_raw[, "Std. Error"], 6),
    t_value        = round(hyp_raw[, "t value"],    4),
    p_value        = round(hyp_raw[, "Pr(>|t|)"],   4),
    Interpretation = ifelse(hyp_raw[, "Pr(>|t|)"] < 0.05,
                            "Significant", "Not Significant"),
    stringsAsFactors = FALSE)
  diagnostic_interpretation(best_model, verbose = verbose)
  return(list(
    best_model = best_model,
    coefficients = best_model$coefficients,
    model_statistics = list(
      R2 = r2,
      RMSE = rmse,
      AIC = round(AIC(best_model$fit), 4),
      BIC = round(BIC(best_model$fit), 4)
    ),
    hypothesis_test = hyp_df
  ))
}

# ============================================================
# WORD EXPORT
# ============================================================
#' Export Analysis Results to a Word Document
#'
#' Assembles all analysis outputs into a formatted Word document and
#' opens it for saving.
#'
#' @param comparison_table A data frame as returned by
#'   \code{compare_models()}.
#' @param best_model A model result list for the best model.
#' @param best_name Character. Name of the best model.
#' @param anova_tbl A data frame as returned by \code{anova_nls()}.
#' @param diag_results A list of diagnostic test objects as returned by
#'   \code{diagnostic_interpretation()}.
#' @param plot_path Character. File path to the saved residual
#'   diagnostics PNG.
#' @param obs_pred_plot_path Character. File path to the saved predicted
#'   vs observed PNG.
#' @param merged_df A data frame of observed, predicted, residual, and
#'   95\% CI values.
#' @param drying_curve_path Character. File path to the saved drying
#'   curve PNG.
#' @param dat A data frame with columns \code{time} and \code{MR}.
#' @param verbose Logical. If \code{TRUE}, prints a confirmation message
#'   to the console. Default \code{FALSE}.
#'
#' @return Called for its side effect. Returns \code{NULL} invisibly.
#'
#' @keywords internal
#' @noRd
export_to_word <- function(comparison_table, best_model, best_name,
                           anova_tbl, diag_results,
                           plot_path, obs_pred_plot_path, merged_df,
                           drying_curve_path, dat, verbose = FALSE) {
  doc <- officer::read_docx()

  doc <- doc |>
    officer::body_add_par("Drying Kinetics Model Analysis", style = "heading 1") |>
    officer::body_add_par("", style = "Normal")

  doc <- doc |>
    officer::body_add_par("1. Model Comparison", style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    officer::body_add_par("1a. Statistical Comparison", style = "heading 3") |>
    officer::body_add_par("", style = "Normal") |>
    flextable::body_add_flextable(
      flextable::flextable(comparison_table[, c("Kinetic Drying Models",
                                                "R2", "RMSE", "MAE", "chi2", "RSS")]) |>
        flextable::set_header_labels(
          `Kinetic Drying Models` = "Kinetic Drying Models",
          R2   = "R\u00b2", RMSE = "RMSE", MAE = "MAE",
          chi2 = "\u03c7\u00b2", RSS  = "RSS") |>
        flextable::bold(part = "header") |>
        flextable::bg(i = 1, bg = "#D9E1F2", part = "body") |>
        flextable::autofit() |> flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal") |>
    officer::body_add_par("1b. Model Details", style = "heading 3") |>
    officer::body_add_par("", style = "Normal") |>
    flextable::body_add_flextable(
      flextable::flextable(comparison_table[, c("Kinetic Drying Models",
                                                "Functional Form",
                                                "Estimated Coefficients")]) |>
        flextable::set_header_labels(
          `Kinetic Drying Models`  = "Kinetic Drying Models",
          `Functional Form`        = "Functional Form",
          `Estimated Coefficients` = "Estimated Coefficients") |>
        flextable::bold(part = "header") |>
        flextable::bg(i = 1, bg = "#D9E1F2", part = "body") |>
        flextable::width(j = "Kinetic Drying Models",  width = 1.5) |>
        flextable::width(j = "Functional Form",        width = 2.0) |>
        flextable::width(j = "Estimated Coefficients", width = 2.5) |>
        flextable::fontsize(size = 9, part = "all") |>
        flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal")

  obs    <- dat$MR
  res    <- best_model$residuals
  n      <- length(res)
  p      <- length(best_model$coefficients)
  ss_res <- sum(res^2)
  ss_tot <- sum((obs - mean(obs))^2)
  r2     <- round(1 - ss_res / ss_tot, 6)
  rmse   <- round(sqrt(ss_res / (n - p)), 6)
  mae    <- round(mean(abs(res)), 6)

  coef_df <- data.frame(
    Parameter = names(best_model$coefficients),
    Estimate  = round(as.numeric(best_model$coefficients), 6),
    stringsAsFactors = FALSE)

  fit_stats_df <- data.frame(
    Metric = c("R\u00b2", "RMSE", "MAE"),
    Value  = c(r2, rmse, mae),
    stringsAsFactors = FALSE)

  doc <- doc |>
    officer::body_add_par(paste("2. Best Model:", best_name), style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    officer::body_add_par("Model Coefficients", style = "heading 3") |>
    flextable::body_add_flextable(
      flextable::flextable(coef_df) |>
        flextable::autofit() |> flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal") |>
    officer::body_add_par("Goodness-of-Fit Statistics", style = "heading 3") |>
    flextable::body_add_flextable(
      flextable::flextable(fit_stats_df) |>
        flextable::autofit() |> flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal")

  hyp_raw <- coef(summary(best_model$fit))
  hyp_df  <- data.frame(
    Parameter      = rownames(hyp_raw),
    Estimate       = round(hyp_raw[, "Estimate"],   6),
    Std_Error      = round(hyp_raw[, "Std. Error"], 6),
    t_value        = round(hyp_raw[, "t value"],    4),
    p_value        = round(hyp_raw[, "Pr(>|t|)"],   4),
    Interpretation = ifelse(hyp_raw[, "Pr(>|t|)"] < 0.05,
                            "Significant", "Not Significant"),
    stringsAsFactors = FALSE)

  doc <- doc |>
    officer::body_add_par("Hypothesis Testing of Parameters", style = "heading 3") |>
    officer::body_add_par("", style = "Normal") |>
    flextable::body_add_flextable(
      flextable::flextable(hyp_df) |>
        flextable::set_header_labels(
          Parameter = "Parameter", Estimate = "Estimate",
          Std_Error = "Std. Error", t_value = "t-value",
          p_value   = "p-value", Interpretation = "Interpretation") |>
        flextable::bold(part = "header") |>
        flextable::autofit() |> flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal")

  anova_tbl_clean <- anova_tbl
  anova_tbl_clean[is.na(anova_tbl_clean)] <- ""
  doc <- doc |>
    officer::body_add_par("3. ANOVA Table", style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    flextable::body_add_flextable(
      flextable::flextable(anova_tbl_clean) |>
        flextable::autofit() |> flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal")

  diag_df <- data.frame(
    Test = c("Shapiro-Wilk (Normality)",
             "Durbin-Watson (Autocorrelation)",
             "Breusch-Pagan (Homoscedasticity)",
             "Runs Test (Randomness)"),
    Statistic = c(round(diag_results$shapiro$statistic, 4),
                  round(diag_results$dw$statistic,      4),
                  round(diag_results$bp$statistic,      4),
                  round(diag_results$runs$statistic,    4)),
    `p-value` = c(round(diag_results$shapiro$p.value, 4),
                  round(diag_results$dw$p.value,      4),
                  round(diag_results$bp$p.value,      4),
                  round(diag_results$runs$p.value,    4)),
    Interpretation = c(
      ifelse(diag_results$shapiro$p.value > 0.05, "Normal",        "Non-normal"),
      ifelse(diag_results$dw$p.value      > 0.05, "Independent",   "Autocorrelated"),
      ifelse(diag_results$bp$p.value      > 0.05, "Homoscedastic", "Heteroscedastic"),
      ifelse(diag_results$runs$p.value    > 0.05, "Random",        "Non-random")),
    stringsAsFactors = FALSE, check.names = FALSE)

  doc <- doc |>
    officer::body_add_par("4. Residual Diagnostic Tests", style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    flextable::body_add_flextable(
      flextable::flextable(diag_df) |>
        flextable::autofit() |> flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal")

  doc <- doc |>
    officer::body_add_par("5. Residual Diagnostic Plots", style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    officer::body_add_img(src = plot_path, width = 6, height = 4.5) |>
    officer::body_add_par("", style = "Normal")

  doc <- doc |>
    officer::body_add_par("6. Predicted vs Observed Values with 95% CI", style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    flextable::body_add_flextable(
      flextable::flextable(merged_df) |>
        flextable::set_header_labels(
          Time         = "Time",
          Observed_MR  = "Observed MR",
          Predicted_MR = "Predicted MR",
          Residual     = "Residual",
          Lower_95CI   = "Lower 95% CI",
          Upper_95CI   = "Upper 95% CI") |>
        flextable::bold(part = "header") |>
        flextable::autofit() |> flextable::theme_vanilla()) |>
    officer::body_add_par("", style = "Normal")

  doc <- doc |>
    officer::body_add_par("7. Drying Curve", style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    officer::body_add_img(src = drying_curve_path, width = 6, height = 4.5) |>
    officer::body_add_par("", style = "Normal")

  doc <- doc |>
    officer::body_add_par("8. Predicted vs Observed Plot", style = "heading 2") |>
    officer::body_add_par("", style = "Normal") |>
    officer::body_add_img(src = obs_pred_plot_path, width = 5, height = 4) |>
    officer::body_add_par("", style = "Normal")

  tmp <- tempfile(fileext = ".docx")
  print(doc, target = tmp)
  if (interactive()) {
    browseURL(tmp)
  }
  if (verbose) {
    message("Word document saved to: ", tmp)
  }
}

# ============================================================
# MASTER FUNCTION
# ============================================================
#' Drying Kinetic Model Comparison and Analysis
#'
#' Fits 20 thin-layer drying kinetic models to experimental moisture
#' ratio data, ranks them by a composite goodness-of-fit criterion
#' (R\eqn{^2}, RMSE, and \eqn{\chi^2}), identifies the best model,
#' predicts moisture ratio (MR) using the best model,
#' performs residual diagnostic tests, and exports a fully formatted
#' report to a Word document.
#'
#' The function expects the input to have at least two numeric columns.
#' The **first** numeric column is treated as drying time as 'time' and the
#' **second** as moisture ratio as 'MR'. Reorder your columns in Excel if
#' needed before calling the function.
#'
#' @param file_path Either a character string giving the path to an
#'   \code{.xlsx} file, or a data frame. In both cases the first two
#'   numeric columns are used as \code{time} and \code{MR}.
#' @param verbose Logical. If \code{TRUE}, also prints detailed
#'   diagnostic and progress information to the console via
#'   \code{message()}. Default \code{FALSE}. Regardless of this
#'   setting, a short note confirming which columns were used as
#'   TIME and MR is always printed, since choosing the wrong columns
#'   silently would give incorrect results.
#' @param models Optional character vector naming which of the 20
#'   models to fit (see Details for valid names). If \code{NULL}
#'   (default), all 20 models are fitted, preserving prior behavior.
#'   Fitting fewer models is faster and is mainly useful for quick
#'   checks or small examples.
#' @param export_word Logical. If \code{TRUE}, builds and exports a
#'   Word document summarizing the model comparison, diagnostics, and
#'   plots. If \code{FALSE} (default), skips Word export and only
#'   returns the results list.
#'
#' @return A list (returned invisibly) with five elements:
#'   \describe{
#'     \item{comparison_table}{Data frame of all model fit statistics,
#'       sorted by composite rank.}
#'     \item{best_model}{The result list for the best-fitting model.}
#'     \item{anova_table}{ANOVA table for the best model.}
#'     \item{predicted_MR}{Data frame containing observed time, observed
#'       MR, predicted MR, residuals, and 95\% confidence intervals from
#'       the best-fitting model.}
#'     \item{columns_used}{Data frame recording which input columns were
#'       used as TIME and MR. Always check this to confirm correct
#'       column selection, especially when \code{verbose = FALSE}.}
#'   }
#'
#' @details
#' The 20 models fitted are: Lewis, Page, Modified Page,
#' Henderson & Pabis, Logarithmic, Two-Term, Two-Term Exponential,
#' Diffusion Approximation, Wang & Singh, Midilli-Kucuk,
#' Modified Henderson & Pabis, Verma, Weibull, Aghbashlo et al.,
#' Jena & Das, Hii et al., Parabolic, Thompson, Demir et al., and
#' Thin-Layer Exponential-Linear. The corresponding names to use with
#' the \code{models} argument are \code{"Lewis"}, \code{"Page"},
#' \code{"Modified_Page"}, \code{"Henderson_Pabis"},
#' \code{"Logarithmic"}, \code{"Two_Term"},
#' \code{"Two_Term_Exponential"}, \code{"Diffusion_Approximation"},
#' \code{"Wang_Singh"}, \code{"Midilli_Kucuk"},
#' \code{"Modified_Henderson_Pabis"}, \code{"Verma"}, \code{"Weibull"},
#' \code{"Aghbashlo"}, \code{"Jena_Das"}, \code{"Hii_et_al"},
#' \code{"Parabolic"}, \code{"Thompson"}, \code{"Demir_et_al"}, and
#' \code{"Thin_Layer_ExpLin"}.
#'
#' Starting values are obtained via a coarse grid search using
#' \code{\link[minpack.lm]{nlsLM}}, making the fitting robust across
#' a wide range of datasets without requiring manual starting value
#' specification.
#'
#' Goodness-of-fit statistics reported are R\eqn{^2}, RMSE, MAE,
#' \eqn{\chi^2}, and RSS. Model ranking uses the combined rank of
#' R\eqn{^2}, RMSE, and \eqn{\chi^2} following Goyal et al. (2007).
#'
#' Residual diagnostics include the Shapiro-Wilk test (normality),
#' Durbin-Watson test (autocorrelation), Breusch-Pagan test
#' (homoscedasticity), and Runs test (randomness).
#'
#' @references
#' Goyal, R. K., Kingsly, A. R. P., Manikantan, M. R., & Ilyas, S. M.
#' (2007). Mathematical modelling of thin layer drying kinetics of plum
#' in a tunnel dryer. \emph{Journal of Food Engineering}, 79(1), 176--180.
#' \doi{10.1016/j.jfoodeng.2006.01.041}
#'
#' @examples
#' # Small, fast toy example (runs automatically during R CMD check)
#' toy_df <- data.frame(
#'   time = c(0, 30, 60, 90, 120, 150),
#'   MR   = c(1.00, 0.72, 0.51, 0.36, 0.26, 0.18)
#' )
#' toy_result <- dryingkineticmodels(toy_df, models = c("Lewis", "Page"),
#'                                    export_word = FALSE)
#' toy_result$comparison_table
#'
#' \donttest{
#' # Fuller example using a bundled sample dataset (slower, includes Word export)
#' sample_file <- system.file("extdata", "sample_drying_data.csv",
#'                            package = "dryingkineticmodels")
#' df <- read.csv(sample_file)
#' result <- dryingkineticmodels(df, export_word = TRUE)
#' }
#'
#' @export
dryingkineticmodels <- function(file_path, verbose = FALSE, models = NULL,
                                export_word = FALSE) {
  if (is.character(file_path)) {
    dat <- readxl::read_excel(file_path)
  } else {
    dat <- file_path
  }

  num_cols <- names(dat)[sapply(dat, is.numeric)]
  if (length(num_cols) < 2) {
    stop("At least 2 numeric columns required. Please check your Excel file.")
  }
  time_col <- num_cols[1]
  mr_col   <- num_cols[2]

  if (verbose) {
    message(paste(
      "NOTE: This function assumes:",
      "  - 1st numeric column = TIME (e.g., time)",
      "  - 2nd numeric column = MR   (e.g., moisture ratio)",
      "  Please ensure your Excel file follows this column order.\n",
      paste0("Using column '", time_col, "' as TIME"),
      paste0("Using column '", mr_col,   "' as MR"),
      "If this is wrong, reorder your Excel columns and try again.",
      sep = "\n"))
  }

  dat        <- dat[, c(time_col, mr_col)]
  names(dat) <- c("time", "MR")

  results          <- fit_all_models(dat, models = models)
  comparison_table <- compare_models(results)

  if (verbose) {
    message("Model comparison table computed (see returned $comparison_table).")
  }

  best_name  <- comparison_table$`Kinetic Drying Models`[1]
  best_model_out <- get_best_model(results, comparison_table, dat, verbose = verbose)
  best_model     <- best_model_out$best_model
  anova_tbl      <- anova_nls(best_model, dat)

  tmp_plot <- tempfile(fileext = ".png")
  png(tmp_plot, width = 800, height = 600, res = 120)
  residual_diagnostics(best_model)
  dev.off()
  residual_diagnostics(best_model)

  obs  <- dat$MR
  pred <- best_model$predicted
  r2   <- round(1 - sum(best_model$residuals^2) /
                  sum((obs - mean(obs))^2), 4)

  tmp_obs_pred <- tempfile(fileext = ".png")
  png(tmp_obs_pred, width = 600, height = 500, res = 120)
  plot(obs, pred,
       xlab = "Observed MR", ylab = "Predicted MR",
       main = "Predicted vs Observed MR",
       pch = 16, col = "steelblue")
  abline(0, 1, lty = 2, col = "red")
  legend("topleft", legend = paste0("R\u00b2 = ", r2), bty = "n")
  dev.off()

  plot(obs, pred,
       xlab = "Observed MR", ylab = "Predicted MR",
       main = "Predicted vs Observed MR",
       pch = 16, col = "steelblue")
  abline(0, 1, lty = 2, col = "red")
  legend("topleft", legend = paste0("R\u00b2 = ", r2), bty = "n")

  pred_vals <- predict(best_model$fit)
  se_pred   <- sqrt(sum(best_model$residuals^2) /
                      (length(best_model$residuals) -
                         length(best_model$coefficients)))
  t_val     <- qt(0.975, df = length(best_model$residuals) -
                    length(best_model$coefficients))
  upper_ci  <- pred_vals + t_val * se_pred
  lower_ci  <- pred_vals - t_val * se_pred
  obs_vals  <- best_model$predicted + best_model$residuals

  tmp_drying_curve <- tempfile(fileext = ".png")
  png(tmp_drying_curve, width = 700, height = 500, res = 120)
  plot(dat$time, obs_vals,
       xlab = "Time", ylab = "Moisture Ratio (MR)",
       main = "Drying Curve",
       pch  = 16, col = "darkgreen",
       ylim = range(c(obs_vals, upper_ci, lower_ci)))
  lines(dat$time, pred_vals, col = "steelblue",  lwd = 2)
  lines(dat$time, upper_ci,  col = "tomato",     lwd = 1.5, lty = 2)
  lines(dat$time, lower_ci,  col = "darkorange", lwd = 1.5, lty = 2)
  legend("topright",
         legend = c("Observed", "Predicted", "Upper 95% CI", "Lower 95% CI"),
         col    = c("darkgreen", "steelblue", "tomato", "darkorange"),
         pch    = c(16, NA, NA, NA),
         lty    = c(NA, 1, 2, 2),
         lwd    = c(NA, 2, 1.5, 1.5),
         bty    = "n")
  dev.off()

  plot(dat$time, obs_vals,
       xlab = "Time", ylab = "Moisture Ratio (MR)",
       main = "Drying Curve",
       pch  = 16, col = "darkgreen",
       ylim = range(c(obs_vals, upper_ci, lower_ci)))
  lines(dat$time, pred_vals, col = "steelblue",  lwd = 2)
  lines(dat$time, upper_ci,  col = "tomato",     lwd = 1.5, lty = 2)
  lines(dat$time, lower_ci,  col = "darkorange", lwd = 1.5, lty = 2)
  legend("topright",
         legend = c("Observed", "Predicted", "Upper 95% CI", "Lower 95% CI"),
         col    = c("darkgreen", "steelblue", "tomato", "darkorange"),
         pch    = c(16, NA, NA, NA),
         lty    = c(NA, 1, 2, 2),
         lwd    = c(NA, 2, 1.5, 1.5),
         bty    = "n")

  res    <- best_model$residuals
  fitted <- best_model$predicted
  diag_results <- list(
    shapiro = shapiro.test(res),
    dw      = lmtest::dwtest(lm(res ~ fitted)),
    bp      = lmtest::bptest(lm(res ~ fitted)),
    runs    = tseries::runs.test(as.factor(ifelse(res >= median(res), "+", "-"))))

  merged_df <- data.frame(
    Time         = dat$time,
    Observed_MR  = round(best_model$predicted + best_model$residuals, 6),
    Predicted_MR = round(pred_vals,                                   6),
    Residual     = round(best_model$residuals,                        6),
    Lower_95CI   = round(pred_vals - t_val * se_pred,                 6),
    Upper_95CI   = round(pred_vals + t_val * se_pred,                 6),
    stringsAsFactors = FALSE)

  if (verbose) {
    cat("\nModel Comparison Table:\n")
    print(comparison_table)

    cat("\nANOVA Table:\n")
    print(anova_tbl)

    cat("\nHypothesis Testing of Parameters:\n")
    print(best_model_out$hypothesis_test)

    cat("\nPredicted vs Observed Values with 95% CI:\n")
    print(merged_df)
  }

  if (export_word) {
    tryCatch({
      export_to_word(comparison_table, best_model, best_name,
                     anova_tbl, diag_results,
                     plot_path          = tmp_plot,
                     obs_pred_plot_path = tmp_obs_pred,
                     merged_df          = merged_df,
                     drying_curve_path  = tmp_drying_curve,
                     dat                = dat,
                     verbose            = verbose)
    }, error = function(e) {
      warning(paste("Word export failed:", e$message), call. = FALSE)
    })
  }

  return(invisible(list(
    comparison_table = comparison_table,
    best_model       = best_model,
    anova_table      = anova_tbl,
    predicted_MR     = merged_df,
    columns_used     = data.frame(
      role   = c("TIME", "MR"),
      column = c(time_col, mr_col),
      stringsAsFactors = FALSE))))
}