R Under development (unstable) (2025-02-23 r87804 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ### > # test functionality of stabilization > > require("gamboostLSS") Loading required package: gamboostLSS Loading required package: mboost Loading required package: parallel Loading required package: stabs Attaching package: 'gamboostLSS' The following object is masked from 'package:stats': model.weights > > ## simulate Gaussian data > set.seed(0804) > x1 <- runif(1000) > x2 <- runif(1000) > x3 <- runif(1000) > x4 <- runif(1000) > mu <- 1.5 +1 * x1 +4 * x2 > sigma <- exp(1 - 0.2 * x3 -0.4 * x4) > y <- rnorm(mean = mu, sd = sigma, n = length(mu)) > dat <- data.frame(x1, x2, x3, x4, y) > > ## do not use stabilization > m1 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families=GaussianLSS(), + data=dat) > coef(m1) $mu (Intercept) x1 x2 -1.9916597 0.2601562 3.6573948 attr(,"offset") [1] 4.075487 $sigma (Intercept) x1 x2 x3 x4 0.07184243 -0.08466051 0.14507711 -0.20507878 -0.34966371 attr(,"offset") [1] 0.8534752 > > ## use stabilization via options (for backwards compatibility) > options(gamboostLSS_stab_ngrad = TRUE) > m2 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families=GaussianLSS(), + data=dat) Warning messages: 1: In check_stabilization(stabilization) : Usage of 'options(gamboostLSS_stab_ngrad = TRUE)' is deprecated. Use argument 'stabilization' in the fitting family. See ?Families for details. 2: In check_stabilization(stabilization) : 'stabilization' is set to "MAD" > coef(m2) $mu (Intercept) x1 x2 -2.5558491 0.8170992 4.2184347 attr(,"offset") [1] 4.075487 $sigma (Intercept) x1 x2 x3 x4 0.1071746 -0.1188396 0.1373902 -0.2372693 -0.3712830 attr(,"offset") [1] 0.8534752 > options(gamboostLSS_stab_ngrad = FALSE) > > ## now use novel interface via families: > m3 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = GaussianLSS(stabilization = "MAD"), + data=dat) > stopifnot(all.equal(coef(m3), coef(m2))) > > ## check if everything is handled correctly > GaussianLSS(stabilization = "MAD") $mu Normal distribution: mu(id link) Loss function: -dnorm(x = y, mean = f, sd = sigma, log = TRUE) $sigma Normal distribution: sigma (log link) Loss function: -dnorm(x = y, mean = mu, sd = exp(f), log = TRUE) attr(,"class") [1] "families" attr(,"qfun") function (p, mu = 0, sigma = 1, lower.tail = TRUE, log.p = FALSE) { qnorm(p = p, mean = mu, sd = sigma, lower.tail = lower.tail, log.p = log.p) } attr(,"name") [1] "Gaussian" > GaussianLSS(stabilization = "L2") $mu Normal distribution: mu(id link) Loss function: -dnorm(x = y, mean = f, sd = sigma, log = TRUE) $sigma Normal distribution: sigma (log link) Loss function: -dnorm(x = y, mean = mu, sd = exp(f), log = TRUE) attr(,"class") [1] "families" attr(,"qfun") function (p, mu = 0, sigma = 1, lower.tail = TRUE, log.p = FALSE) { qnorm(p = p, mean = mu, sd = sigma, lower.tail = lower.tail, log.p = log.p) } attr(,"name") [1] "Gaussian" > GaussianLSS(stabilization = "none") $mu Normal distribution: mu(id link) Loss function: -dnorm(x = y, mean = f, sd = sigma, log = TRUE) $sigma Normal distribution: sigma (log link) Loss function: -dnorm(x = y, mean = mu, sd = exp(f), log = TRUE) attr(,"class") [1] "families" attr(,"qfun") function (p, mu = 0, sigma = 1, lower.tail = TRUE, log.p = FALSE) { qnorm(p = p, mean = mu, sd = sigma, lower.tail = lower.tail, log.p = log.p) } attr(,"name") [1] "Gaussian" > res <- try(GaussianLSS(stabilization = "test"), silent = TRUE) > res [1] "Error in match.arg(stabilization) : \n 'arg' should be one of \"none\", \"MAD\", \"L2\"\n" attr(,"class") [1] "try-error" attr(,"condition") > > > ############################################################ > ## continue these checks for other families > dat$y <- runif(1000, min = 0.01, max = 0.99) > FAMILIES <- list( + GaussianLSS, + GammaLSS, + BetaLSS, + StudentTLSS) > for (i in 1:length(FAMILIES)) { + m_none <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "none"), + data=dat) + m_MAD <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "MAD"), + data=dat) + m_L2 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "L2"), + data=dat) + + stopifnot(tail(risk(m_none, merge = TRUE), 1) != tail(risk(m_MAD, merge = TRUE), 1)) + cat('Risks:\n stabilization = "none":', + tail(risk(m_none, merge = TRUE), 1), + '\n stabilization = "MAD":', + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") + } Risks: stabilization = "none": 126.3867 stabilization = "MAD": 126.3865 stabilization = "L2": 126.3866 Risks: stabilization = "none": 152.9363 stabilization = "MAD": 151.7991 stabilization = "L2": 152.3966 Risks: stabilization = "none": -16.74875 stabilization = "MAD": -16.85658 stabilization = "L2": -16.81103 Risks: stabilization = "none": 126.3869 stabilization = "MAD": 126.3865 stabilization = "L2": 126.3868 > > ## check as.families interface for 2:4 parametric families > dat$y <- rnorm(1000, mean = 10, sd = 1) > FAMILIES <- list( + "NO", + "TF") > for (i in 1:length(FAMILIES)) { + m_none <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "none"), + data=dat) + m_MAD <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "MAD"), + data=dat) + m_L2 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "L2"), + data=dat) + cat('Risks:\n stabilization = "none":', + tail(risk(m_none, merge = TRUE), 1), + '\n stabilization = "MAD":', + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") + } Risks: stabilization = "none": 1417.42 stabilization = "MAD": 1417.416 stabilization = "L2": 1417.42 Risks: stabilization = "none": 1424.665 stabilization = "MAD": 1417.082 stabilization = "L2": 1419.756 > > FAMILIES <- list("BCT") > require("gamlss.dist") Loading required package: gamlss.dist Attaching package: 'gamlss.dist' The following object is masked from 'package:mboost': Family > dat$y <- rBCT(1000, mu = 100, sigma = 0.1, nu = 0, tau = 2) > for (i in 1:length(FAMILIES)) { + m_none <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "none"), + data=dat) + m_MAD <- try(glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "MAD"), + data=dat), silent = TRUE) + if (inherits(m_MAD, "try-error")) { + warning("BCT cannot be fitted with stabilization = 'MAD'", immediate. = TRUE) + break + } + + m_L2 <- try(glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "L2"), + data=dat), silent = TRUE) + if (inherits(m_L2, "try-error")) { + warning("BCT cannot be fitted with stabilization = 'L2'", immediate. = TRUE) + break + } + cat('Risks:\n stabilization = "none":', + tail(risk(m_none, merge = TRUE), 1), + '\n stabilization = "MAD":', + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") + } Risks: stabilization = "none": 4336.164 stabilization = "MAD": 4289.131 stabilization = "L2": 4303.955 > > ## check survival families > dat$zens <- sample(c(0, 1), 1000, replace = TRUE) > FAMILIES <- list( + LogNormalLSS, + WeibullLSS, + LogLogLSS) > families <- c("LogNormalLSS", "WeibullLSS", "LogLogLSS") > require(survival) Loading required package: survival > for (i in 1:length(FAMILIES)) { + cat(families[i], "\n\n") + m_none <- glmboostLSS(Surv(y, zens) ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "none"), + data=dat) + m_MAD <- try(glmboostLSS(Surv(y, zens) ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "MAD"), + data=dat), silent = TRUE) + + if (inherits(m_MAD, "try-error")) { + warning(families[i], " cannot be fitted with stabilization = 'MAD'", immediate. = TRUE) + break + } + m_L2 <- try(glmboostLSS(Surv(y, zens) ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "L2"), + data=dat), silent = TRUE) + + if (inherits(m_L2, "try-error")) { + warning(families[i], " cannot be fitted with stabilization = 'L2'", immediate. = TRUE) + break + } + cat('Risks:\n stabilization = "none":', + tail(risk(m_none, merge = TRUE), 1), + '\n stabilization = "MAD":', + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") + } LogNormalLSS Risks: stabilization = "none": 341.8219 stabilization = "MAD": 341.7978 stabilization = "L2": 342.3934 WeibullLSS Risks: stabilization = "none": 3610.826 stabilization = "MAD": 705.8757 stabilization = "L2": 599.4504 LogLogLSS Risks: stabilization = "none": 522.9556 stabilization = "MAD": 538.4426 stabilization = "L2": 518.0424 > > ## check count data > dat$y <- rnbinom(1000, size=10, mu=5) > FAMILIES <- list( + NBinomialLSS, + ZIPoLSS, + ZINBLSS) > for (i in 1:length(FAMILIES)) { + m_none <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "none"), + data=dat) + m_MAD <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "MAD"), + data=dat) + m_L2 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "L2"), + data=dat) + + cat('Risks:\n stabilization = "none":', + tail(risk(m_none, merge = TRUE), 1), + '\n stabilization = "MAD":', + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") + } Risks: stabilization = "none": 2364.811 stabilization = "MAD": 2363.998 stabilization = "L2": 2364.279 Risks: stabilization = "none": 2414.483 stabilization = "MAD": 2401.002 stabilization = "L2": 2393.06 Risks: stabilization = "none": 2366.663 stabilization = "MAD": 2362.5 stabilization = "L2": 2364.284 > > proc.time() user system elapsed 26.60 0.70 27.29