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.

> ###
> # Bugfixes
> 
> 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

> require("gamlss")
Loading required package: 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

Attaching package: 'gamlss.dist'

The following object is masked from 'package:mboost':

    Family

Loading required package: nlme
 **********   GAMLSS Version 5.4-22  ********** 
For more on GAMLSS look at https://www.gamlss.com/
Type gamlssNews() to see new features/changes/bug fixes.

> 
> ## subset method was missing if initial mstop = 1
> set.seed(1907)
> x1 <- rnorm(1000)
> x2 <- rnorm(1000)
> x3 <- rnorm(1000)
> x4 <- rnorm(1000)
> x5 <- rnorm(1000)
> x6 <- rnorm(1000)
> mu    <- exp(1.5 + x1^2 +0.5 * x2 - 3 * sin(x3) -1 * x4)
> sigma <- exp(-0.2 * x4 +0.2 * x5 +0.4 * x6)
> y <- numeric(1000)
> for( i in 1:1000)
+     y[i] <- rnbinom(1, size = sigma[i], mu = mu[i])
> dat <- data.frame(x1, x2, x3, x4, x5, x6, y)
> 
> model <- gamboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
+                      control = boost_control(mstop = 1))
> model[10]

	 LSS Models fitted via Model-based Boosting

Call:
gamboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(),     control = boost_control(mstop = 1))

Number of boosting iterations (mstop):  mu = 10, sigma = 10 
Step size:  mu = 0.1, sigma = 0.1 

Families:

	 Negative Negative-Binomial Likelihood: mu (log link) 

Loss function: -(lgamma(y + sigma) - lgamma(sigma) - lgamma(y + 1) + sigma *  
     log(sigma) + y * f - (y + sigma) * log(exp(f) + sigma)) 
 

	 Negative Negative-Binomial Likelihood: sigma (log link) 

Loss function: -(lgamma(y + exp(f)) - lgamma(exp(f)) - lgamma(y + 1) + exp(f) *  
     f + y * log(mu) - (y + exp(f)) * log(mu + exp(f))) 
 
> 
> 
> ## selected() was broken (didn't call mboost-function)
> stopifnot(all.equal(selected(model),
+                     list(mu = mboost::selected(model[[1]]),
+                          sigma = mboost::selected(model[[2]]))))
> 
> # If the families argument is not specified explicitly in mboostLSS one gets an
> # error in cvrisk.mboostLSS() (spotted by Almond Stöcker).
> # (https://github.com/boost-R/gamboostLSS/issues/9)
> set.seed(1907)
> x1 <- rnorm(1000)
> mu    <- 2*x1
> sigma <- rep(1, 1000)
> y <- numeric(1000)
> for( i in 1:1000)
+        y[i] <- rnorm(1, mean = mu[i], sd=sigma[i])
> dat <- data.frame(x1=x1, y=y)
> ## model with default families
> model <- mboostLSS(y ~ bbs(x1), data = dat, control = boost_control(mstop = 2))
> ## cvrisk.mboostLSS() does not work as families was not specified in model call
> cvr <- try(cvrisk(model, folds = cv(model.weights(model), B=3), trace=FALSE),
+            silent = TRUE)
Warning message:
In make.grid(mstop(object)) :
  Duplicates produced; Only unique values are returned
> if (inherits(cvr, "try-error"))
+     stop("cvrisk does not work if no family was (explicitly) chosen")
> 
> 
> # Make sure that gamlss.dist::BB and gamlss.dist::BI work (spotted by F. Scheipl)
> # (https://github.com/boost-R/gamboostLSS/issues/12)
> set.seed(123)
> n <- 100
> x <- rnorm(n)
> z <- rnorm(n)
> data <- data.frame(y = rbinom(n, p = plogis(x + z), size = 60), x = x, z= z)
> data$ymat <- with(data, cbind(success = data$y, fail = 60 - data$y))
> 
> m1 <- gamlss(ymat ~ x + z, data = data, family = BB)
GAMLSS-RS iteration 1: Global Deviance = 699.795 
GAMLSS-RS iteration 2: Global Deviance = 534.6991 
GAMLSS-RS iteration 3: Global Deviance = 517.8675 
GAMLSS-RS iteration 4: Global Deviance = 517.8664 
GAMLSS-RS iteration 5: Global Deviance = 517.8664 
> m2 <- gamlss(ymat ~ x + z, data = data, family = BI)
GAMLSS-RS iteration 1: Global Deviance = 517.8684 
GAMLSS-RS iteration 2: Global Deviance = 517.8684 
> # same with boosting
> m3 <- glmboostLSS(ymat ~ x + z, data = data, families = as.families("BB"))
> m4 <- glmboost(ymat ~ x + z, data = data, family = as.families("BI"))
Warning message:
In as.families("BI") :
  For boosting one-parametric families, please use the mboost package.
> 
> round(data.frame(BB_gamlss = coef(m1),
+                  BI_gamlss = coef(m2),
+                  BB_gamboostLSS = coef(m3, off2int = TRUE, parameter = "mu"),
+                  BI_gamboostLSS = coef(m4, off2int = TRUE)), 3)
            BB_gamlss BI_gamlss BB_gamboostLSS BI_gamboostLSS
(Intercept)    -0.014    -0.014         -0.006         -0.014
x               1.009     1.009          0.964          1.009
z               0.965     0.965          0.925          0.965
> 
> 
> ## make sure that combined_risk is not written to the global environment 
> ## (and thus replaced if another model is fitted)
> set.seed(1907)
> x1 <- rnorm(1000)
> x2 <- rnorm(1000)
> x3 <- rnorm(1000)
> x4 <- rnorm(1000)
> x5 <- rnorm(1000)
> x6 <- rnorm(1000)
> mu    <- exp(1.5 +1 * x1 +0.5 * x2 -0.5 * x3 -1 * x4)
> sigma <- exp(-0.4 * x3 -0.2 * x4 +0.2 * x5 +0.4 * x6)
> y <- numeric(1000)
> for( i in 1:1000)
+     y[i] <- rnbinom(1, size = sigma[i], mu = mu[i])
> dat <- data.frame(x1, x2, x3, x4, x5, x6, y)
> 
> model_0 <- mboostLSS(x1 ~ ., families = GaussianLSS(), data = dat, 
+                      control = boost_control(mstop = 30),method = "noncyclic")
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 8 row(s) of zeros
2: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 8 row(s) of zeros
> length_0 <- length(risk(model_0, merge = TRUE))
> if (length_0 != 32)
+     stop("combined risk not correct.")
> 
> model_1 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, 
+                        control = boost_control(mstop = 20), method = "noncyclic")
> length_1 <- length(risk(model_1, merge = TRUE))
> if (length_1 != 22)
+     stop("combined risk not correct.")
> if (length(risk(model_0, merge = TRUE)) != length_0)
+     stop("Combined risk overwritten by new model. Scoping error.")
> 
> model_2 <- mboostLSS(x1 ~ ., families = GaussianLSS(), data = dat, 
+                      control = boost_control(mstop = 11),method = "cyclic")
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 8 row(s) of zeros
2: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 8 row(s) of zeros
> if (length(risk(model_2, merge = TRUE)) != 24)
+     stop("combined risk not correct.")
> if (length(risk(model_0, merge = TRUE)) != length_0)
+     stop("Combined risk overwritten by new model. Scoping error.")
> if (length(risk(model_1, merge = TRUE)) != length_1)
+     stop("Combined risk overwritten by new model. Scoping error.")
> 
> # fix predint when multiple variables are included within a single baselearner #55
> # (https://github.com/boost-R/gamboostLSS/issues/55)
> 
> set.seed(7)
> data <- data.frame(x1 = runif(400, -1, 1),
+                    x2 = runif(400, -1, 1),
+                    x3 = runif(400, -1, 1),
+                    x4 = runif(400, -1, 1))
> data$y <- rnorm(400, 4 * data$x1 + 2 * data$x3)
> 
> gb <- gamboostLSS(y ~ btree(x1, x2, x3, x4), data = data, method = "noncyclic", control = boost_control(mstop = 10))
> p <- predint(gb, which = "x1")
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   2.21    0.21    2.43