## ----setup, include = FALSE---------------------------------------------------
#file.edit(normalizePath("~/.Renviron"))
LOCAL <- identical(Sys.getenv("LOCAL"), "TRUE")
#LOCAL=TRUE
knitr::opts_chunk$set(purl = LOCAL)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
library(SelectBoost.gamlss)
library(gamlss)

## ----cache=TRUE, eval=LOCAL---------------------------------------------------
bc <- NULL
if (requireNamespace("pscl", quietly = TRUE)) {
  data("bioChemists", package = "pscl")
  bc <- pscl::bioChemists
  needed <- c("art", "kid5", "phd", "fem", "mar", "ment")
  missing_cols <- setdiff(needed, names(bc))
  if (length(missing_cols)) {
    cat(
      "bioChemists dataset is missing columns:",
      paste(missing_cols, collapse = ", "),
      "— skipping ZIP example.\n"
    )
    bc <- NULL
  } else {
    keep <- stats::complete.cases(bc[, needed])
    bc <- bc[keep, , drop = FALSE]
    bc$fem <- factor(bc$fem)
    bc$mar <- factor(bc$mar)
    bc$kid5 <- factor(bc$kid5)

    # ZIP engine via stepwise (mu) + grouped/elastic features
    set.seed(10)
    fit_zip <- sb_gamlss(
      art ~ 1, data = bc, family = gamlss.dist::ZIP(),
      mu_scope = ~ kid5 + fem + mar + phd + ment,
      engine = "grpreg", grpreg_penalty = "grLasso",
      B = 80, pi_thr = 0.6, pre_standardize = FALSE, parallel = "auto", trace = FALSE
    )
    sel_zip <- selection_table(fit_zip)
    sel_zip <- sel_zip[order(sel_zip$parameter, -sel_zip$prop), , drop = FALSE]
    head(sel_zip, n = min(8L, nrow(sel_zip)))
    stable_zip <- sel_zip[sel_zip$prop >= fit_zip$pi_thr, , drop = FALSE]
    if (nrow(stable_zip)) {
      stable_zip
    } else {
      cat("No terms reached the stability threshold of", fit_zip$pi_thr, "for this run.\n")
    }
  }
} else {
  cat("Package 'pscl' not installed; skipping bioChemists example.\n")
}

## ----cache=TRUE, eval=LOCAL---------------------------------------------------
if (requireNamespace("pscl", quietly = TRUE) && !is.null(bc)) {
  set.seed(11)
  fit_zinb <- sb_gamlss(
    art ~ 1, data = bc, family = gamlss.dist::ZINBI(),
    mu_scope = ~ kid5 + fem + mar + phd + ment,
    engine = "grpreg", grpreg_penalty = "grLasso",
    B = 80, pi_thr = 0.6, pre_standardize = FALSE, parallel = "auto", trace = FALSE
  )
  sel_zinb <- selection_table(fit_zinb)
  sel_zinb <- sel_zinb[order(sel_zinb$parameter, -sel_zinb$prop), , drop = FALSE]
  head(sel_zinb, n = min(8L, nrow(sel_zinb)))
  stable_zinb <- sel_zinb[sel_zinb$prop >= fit_zinb$pi_thr, , drop = FALSE]
  if (nrow(stable_zinb)) {
    stable_zinb
  } else {
    cat("No terms reached the stability threshold of", fit_zinb$pi_thr, "for this run.\n")
  }
} else if (!requireNamespace("pscl", quietly = TRUE)) {
  cat("Package 'pscl' not installed; skipping bioChemists example.\n")
} else {
  cat("bioChemists dataset unavailable; skipping ZINBI example.\n")
}

## ----cache=TRUE, eval=LOCAL---------------------------------------------------
if (requireNamespace("pscl", quietly = TRUE)) {
  if (exists("fit_zip", inherits = FALSE) && exists("fit_zinb", inherits = FALSE)) {
    zip_g <- fit_zip$final_fit
    zinb_g <- fit_zinb$final_fit
    if (inherits(zip_g, "gamlss") && inherits(zinb_g, "gamlss")) {
      safe_scalar <- function(expr) {
        out <- tryCatch(expr, error = function(e) NA_real_)
        if (length(out) == 0L) NA_real_ else as.numeric(out[1L])
      }
      met <- data.frame(
        model = c("ZIP", "ZINBI"),
        AIC = c(safe_scalar(AIC(zip_g)), safe_scalar(AIC(zinb_g))),
        GAIC_k2 = c(safe_scalar(gamlss::GAIC(zip_g, k = 2)), safe_scalar(gamlss::GAIC(zinb_g, k = 2))),
        deviance = c(safe_scalar(zip_g$deviance), safe_scalar(zinb_g$deviance))
      )
      if (any(!is.na(met$AIC)) || any(!is.na(met$GAIC_k2)) || any(!is.na(met$deviance))) {
        print(met)
      } else {
        cat("Model fit statistics unavailable; skipping metric table.\n")
      }
    } else {
      cat("Final fits not found; metrics skipped.\n")
    }
  } else {
    cat("Models not available; run the fitting chunks first.\n")
  }
} else {
  cat("Package 'pscl' not installed; skipping bioChemists example.\n")
}

## ----cache=TRUE, eval=LOCAL---------------------------------------------------
aq <- subset(airquality, !is.na(Ozone) & !is.na(Wind) & !is.na(Temp) & !is.na(Solar.R))
aq$Month <- factor(aq$Month)

set.seed(20)
fit_zaga <- sb_gamlss(
  Ozone ~ 1, data = aq, family = gamlss.dist::ZAGA(),
  mu_scope    = ~ pb(Temp) + pb(Wind) + pb(Solar.R) + Month,
  sigma_scope = ~ pb(Temp),
  engine = "grpreg", grpreg_penalty = "grLasso",
  B = 80, pi_thr = 0.6, pre_standardize = TRUE, parallel = "auto", trace = FALSE
)
sel_zaga <- selection_table(fit_zaga)
sel_zaga <- sel_zaga[order(sel_zaga$parameter, -sel_zaga$prop), , drop = FALSE]
head(sel_zaga, n = min(8L, nrow(sel_zaga)))
stable_zaga <- sel_zaga[sel_zaga$prop >= fit_zaga$pi_thr, , drop = FALSE]
if (nrow(stable_zaga)) {
  stable_zaga
} else {
  cat("No terms reached the stability threshold of", fit_zaga$pi_thr, "for this run.\n")
}

## ----cache=TRUE, eval=LOCAL---------------------------------------------------
if (requireNamespace('ggplot2', quietly = TRUE)) {
  print(effect_plot(fit_zaga, 'Temp', aq, what = 'mu'))
}

## ----cache=TRUE, eval=LOCAL---------------------------------------------------
if (requireNamespace("nlme", quietly = TRUE)) {
  data_ok <- tryCatch({
    utils::data("Orthodont", package = "nlme", envir = environment())
    TRUE
  }, error = function(e) FALSE)
  if (!data_ok || !exists("Orthodont", inherits = FALSE)) {
    cat("Dataset 'Orthodont' unavailable; skipping longitudinal example.\n")
  } else {
    od <- get("Orthodont", inherits = FALSE)
    rm(Orthodont)
    od$Subject <- factor(od$Subject)
    od$Sex <- factor(od$Sex)
    set.seed(30)
    fit_long <- sb_gamlss(
      distance ~ gamlss::random(Subject),     # random intercept
      data = od, family = gamlss.dist::NO(),
      mu_scope = ~ pb(age) + Sex + pb(age):Sex,      # select fixed effects
      engine = "grpreg", grpreg_penalty = "grLasso",
      B = 80, pi_thr = 0.6, pre_standardize = TRUE, parallel = "auto", trace = FALSE
    )
  sel_long <- selection_table(fit_long)
  sel_long <- sel_long[order(sel_long$parameter, -sel_long$prop), , drop = FALSE]
  head(sel_long, n = min(8L, nrow(sel_long)))
  stable_long <- sel_long[sel_long$prop >= fit_long$pi_thr, , drop = FALSE]
  if (nrow(stable_long)) {
    stable_long
  } else {
    cat("No terms reached the stability threshold of", fit_long$pi_thr, "for this run.\n")
  }
    if (requireNamespace("ggplot2", quietly = TRUE)) {
      gg <- tryCatch(
        effect_plot(fit_long, "age", od, what = "mu"),
        error = function(e) {
          message("effect_plot() skipped for random-effect example: ", conditionMessage(e))
          NULL
        }
      )
      if (!is.null(gg)) {
        print(gg)
      }
    }
  }
} else {
  cat("Package 'nlme' not installed; skipping longitudinal example.\n")
}