Advanced Real Data Examples: Zero-Inflation, Semicontinuous, and Longitudinal Growth

library(SelectBoost.gamlss)
library(gamlss)
#> Loading required package: splines
#> Loading required package: gamlss.data
#> 
#> Attaching package: 'gamlss.data'
#> The following object is masked from 'package:datasets':
#> 
#>     sleep
#> Loading required package: gamlss.dist
#> Loading required package: nlme
#> Loading required package: parallel
#>  **********   GAMLSS Version 5.5-0  **********
#> For more on GAMLSS look at https://www.gamlss.com/
#> Type gamlssNews() to see new features/changes/bug fixes.

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")
}

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")
}

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")
}