R Under development (unstable) (2026-01-25 r89330 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2026 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.

> ###
> # check families
> 
> require("gamboostLSS")
Loading required package: gamboostLSS
Loading required package: gamlss.dist
Loading required package: mboost
Loading required package: parallel
Loading required package: stabs

Attaching package: 'mboost'

The following object is masked from 'package:gamlss.dist':

    Family


Attaching package: 'gamboostLSS'

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

    model.weights

> 
> ### check families with only one offset specified (other to choose via optim)
> set.seed(1907)
> n <- 5000
> x1  <- runif(n)
> x2 <- runif(n)
> mu <- 2 -1*x1 - 3*x2
> sigma <- exp(-1*x1 + 3*x2)
> df <- exp(1 + 3*x1 + 1*x2)
> y <- rTF(n = n, mu = mu, sigma = sigma, nu = df)
> data <- data.frame(y = y, x1 = x1, x2 = x2)
> rm("y", "x1", "x2")
> 
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(mu = 0.5),
+                      data = data,
+                      control = boost_control(mstop = 10), center = TRUE)
> coef(model)
$mu
(Intercept)          x2 
  0.1075190  -0.2170573 
attr(,"offset")
[1] 0.5

$sigma
(Intercept)          x2 
 -0.7575752   1.5293790 
attr(,"offset")
[1] 1.151955

$df
(Intercept)          x2 
  0.1255079  -0.1746240 
attr(,"offset")
[1] 0.8369583

> 
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(sigma = 1),
+                      data = data,
+                      control = boost_control(mstop = 10), center = TRUE)
> coef(model)
$mu
(Intercept)          x2 
  0.2464260  -0.4974803 
attr(,"offset")
[1] 0.5774545

$sigma
(Intercept)          x2 
 -0.4516042   0.9116904 
attr(,"offset")
[1] 0

$df
(Intercept)          x2 
  0.3149720  -0.6358598 
attr(,"offset")
[1] -0.1776425

> 
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(df = 4),
+                      data = data,
+                      control = boost_control(mstop = 10), center = TRUE)
> coef(model)
$mu
(Intercept)          x2 
  0.1074937  -0.2170063 
attr(,"offset")
[1] 0.5774545

$sigma
(Intercept)          x2 
  -0.889911    1.796536 
attr(,"offset")
[1] 1.152757

$df
(Intercept)          x1          x2 
 0.07962737  0.06912987 -0.23090205 
attr(,"offset")
[1] 1.386294

> 
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(mu = 0.5, df = 1),
+                      data = data,
+                      control = boost_control(mstop = 10), center = TRUE)
> coef(model)
$mu
(Intercept)          x2 
  0.1953042  -0.3942766 
attr(,"offset")
[1] 0.5

$sigma
(Intercept)          x2 
 -0.5914862   1.1940815 
attr(,"offset")
[1] 0.548112

$df
(Intercept)          x2 
  0.1911318  -0.2203124 
attr(,"offset")
[1] 0

> 
> ### multivariate minimum for offset
> loss <- function(df, sigma,y, f){
+     -1 * (lgamma((df+1)/2) - log(sigma) - lgamma(1/2) - lgamma(df/2) - 0.5 *
+           log(df) - (df+1)/2 * log(1 + (y-f)^2 / (df * sigma^2)))
+ }
> riski <- function(x, y, w = rep(1, length(y))){
+     f <- x[[1]]
+     df <- exp(x[[2]])
+     sigma <- exp(x[[3]])
+     sum(w * loss(y = y, f = f, df = df, sigma = sigma))
+ }
> 
> res <- optim(fn = riski, par = c(0, 1, 1), y = data$y, w = rep(1, length(data$y)))$par
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(mu = res[1], sigma =
+                                   exp(res[2]), df = exp(res[3])),
+                      data = data,
+                      control = boost_control(mstop = 10), center = TRUE)
> model[100]

	 LSS Models fitted via Model-based Boosting

Call:
glmboostLSS(formula = y ~ x1 + x2, data = data, families = StudentTLSS(mu = res[1],     sigma = exp(res[2]), df = exp(res[3])), control = boost_control(mstop = 10),     center = TRUE)

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

Families:

	 Student's t-distribution: mu (id link) 

Loss function: { 
     -1 * (lgamma((df + 1)/2) - log(sigma) - lgamma(1/2) - lgamma(df/2) -  
         0.5 * log(df) - (df + 1)/2 * log(1 + (y - f)^2/(df *  
         sigma^2))) 
 } 
 

	 Student's t-distribution: sigma (log link) 

Loss function: { 
     -1 * (lgamma((df + 1)/2) - f - lgamma(1/2) - lgamma(df/2) -  
         0.5 * log(df) - (df + 1)/2 * log(1 + (y - mu)^2/(df *  
         exp(2 * f)))) 
 } 
 

	 Student's t-distribution: df (log link) 

Loss function: { 
     -1 * (lgamma((exp(f) + 1)/2) - log(sigma) - lgamma(1/2) -  
         lgamma(exp(f)/2) - 0.5 * f - (exp(f) + 1)/2 * log(1 +  
         (y - mu)^2/(exp(f) * sigma^2))) 
 } 
 
> coef(model)
$mu
(Intercept)          x1          x2 
  1.2446880  -0.8978331  -1.6342267 
attr(,"offset")
[1] 0.5228725

$sigma
(Intercept)          x1          x2 
 -0.5893432  -0.8476283   2.7681255 
attr(,"offset")
[1] 0.3968482

$df
(Intercept)          x1          x2 
  0.3707909   0.6081607  -0.6771414 
attr(,"offset")
[1] 0.7357846

> 
> ### check as families
> 
> ### Gaussian: two different ways
> model <- glmboostLSS(y ~ x1 + x2, families = as.families("NO"),
+                      data = data,
+                      control = boost_control(mstop = 10), center = TRUE)
> 
> model2 <- glmboostLSS(y ~ x1 + x2, families = GaussianLSS(),
+                       data = data,
+                       control = boost_control(mstop = 10), center = TRUE)
> 
> coef(model, off2int = TRUE)  # as.families("NO")
$mu
(Intercept)          x2 
 0.04778256 -0.13022745 

$sigma
(Intercept)          x1          x2 
  1.0141654  -0.1381848   1.3859472 

> coef(model2, off2int = TRUE) # GaussianLSS()
$mu
(Intercept)          x2 
 0.04778256 -0.13022745 

$sigma
(Intercept)          x1          x2 
  1.0141654  -0.1381848   1.3859472 

> 
> ### change link function inside as.families()
> model2 <- glmboostLSS(abs(y) ~ x1 + x2, families = as.families("NO", mu.link = "log"),
+                       data = data,
+                       control = boost_control(mstop = 10), center = TRUE)
> coef(model2)
$mu
(Intercept)          x2 
 -0.5210051   1.0517957 
attr(,"offset")
[1] 1.248357

$sigma
(Intercept)          x2 
 -0.7673921   1.2107247 
attr(,"offset")
[1] 1.539722

> 
> 
> model3 <- glmboostLSS(abs(y)/(max(y)+.01) ~ x1 + x2, families = as.families("BE", mu.link = "logit",
+                                                                             sigma.link = "log"),
+                       data = data,
+                       control = boost_control(mstop = 10), center = TRUE)
> coef(model3)
$mu
(Intercept)          x2 
 -0.2137522   0.7177991 
attr(,"offset")
[1] -2.931116

$sigma
(Intercept)          x2 
-0.68253369  0.09271978 
attr(,"offset")
[1] -0.6931472

> 
> 
> model4 <- glmboostLSS(y ~ x1 + x2, families = as.families("TF", mu.link = "identity",
+                                                           sigma.link = "log", 
+                                                           nu.link = "log"),
+                       data = data,
+                       control = boost_control(mstop = 10), center = TRUE)
> coef(model4)
$mu
(Intercept)          x2 
 0.05795847 -0.11700550 
attr(,"offset")
[1] -0.01672538

$sigma
(Intercept)          x2 
 -0.6968961   1.0374093 
attr(,"offset")
[1] 1.761574

$nu
(Intercept)          x2 
  0.0572190  -0.1155127 
attr(,"offset")
[1] 2.302585

> 
> ### Additionally use stabilization
> 
> model4 <- glmboostLSS(y ~ x1 + x2, families = as.families("TF", mu.link = "identity",
+                                                           sigma.link = "log", 
+                                                           nu.link = "log", 
+                                                           stabilization = "L2"),
+                       data = data,
+                       control = boost_control(mstop = 10), center = TRUE)
> coef(model4)
$mu
(Intercept)          x2 
  0.4233792  -0.8547101 
attr(,"offset")
[1] -0.01672538

$sigma
(Intercept)          x2 
 -0.6004438   0.8826834 
attr(,"offset")
[1] 1.761574

$nu
(Intercept)          x2 
  0.2748352  -0.5548323 
attr(,"offset")
[1] 2.302585

> 
> 
> model4 <- glmboostLSS(y ~ x1 + x2, families = as.families("TF", mu.link = "identity",
+                                                           sigma.link = "log", 
+                                                           nu.link = "log", 
+                                                           stabilization = "MAD"),
+                       data = data,
+                       control = boost_control(mstop = 10), center = TRUE)
> coef(model4)
$mu
(Intercept)          x2 
   1.146445   -2.314422 
attr(,"offset")
[1] -0.01672538

$sigma
(Intercept)          x1          x2 
  -1.586238   -0.701680    2.484648 
attr(,"offset")
[1] 1.761574

$nu
(Intercept)          x1          x2 
 -0.5822664   3.5450017  -1.0655973 
attr(,"offset")
[1] 2.302585

> 
> 
> ### check survival families
> if (require("survival")) {
+   #set.seed(22256)
+ 
+   x1 <- runif(1000)
+   x2 <- runif(1000)
+   x3 <- runif(1000)
+   w <- rnorm(1000)
+   
+   time <-  exp(3 + 1*x1 +2*x2  + exp(0.2 * x3) * w)
+   status <- rep(1, 1000)
+   
+   ## check if error messages are correct
+   try(glmboost(time ~ x1 + x2 + x3, family = Lognormal(),
+                control = boost_control(trace = TRUE), center = TRUE))
+   try(glmboostLSS(time ~ x1 + x2 + x3, families = LogNormalLSS(),
+                   control = boost_control(trace = TRUE), center = TRUE))
+   try(glmboostLSS(list(mu = Surv(time, status) ~ x1 + x2 + x3,
+                        sigma = time ~ x1 + x2 + x3), families = LogNormalLSS(),
+                   control = boost_control(trace = TRUE), center = TRUE))
+   
+   ## check results
+   # LogNormalLSS()
+   (m1 <- survreg(Surv(time, status) ~ x1 + x2 + x3, dist="lognormal"))
+   model <- glmboostLSS(Surv(time, status) ~ x1 + x2 + x3, families = LogNormalLSS(),
+                        control = boost_control(trace = TRUE), center = TRUE)
+   stopifnot(sum(abs(coef(model, off2int = TRUE)[[1]] - c(3, 1, 2, 0)))
+             < sum(abs(coef(m1) - c(3, 1, 2, 0))))
+   stopifnot(sum(abs(coef(model, off2int = TRUE)[[2]] - c(0, 0, 0, 0.2))) < 0.25)
+   
+   # LogLogLSS()
+   etamu <- 3 + 1*x1 +2*x2
+   etasi <- exp(rep(0.2, 1000))
+   for (i in 1:1000)
+     time[i] <- exp(rlogis(1, location = etamu[i], scale = etasi[i]))
+   status <- rep(1, 1000)
+   (m1 <- survreg(Surv(time, status) ~ x1 + x2 + x3, dist="loglogistic"))
+   model <- glmboostLSS(Surv(time, status) ~ x1 + x2 + x3, families = LogLogLSS(),
+                        control = boost_control(trace = TRUE, mstop = 350), center = TRUE)
+   stopifnot(sum(abs(coef(model, off2int = TRUE, which ="")[[1]] - c(3, 1, 2, 0)))
+             < sum(abs(coef(m1) - c(3, 1, 2, 0))))
+   stopifnot(sum(abs(coef(model, off2int = TRUE)[[2]] - c(0.2, 0, 0, 0))) < 0.25)
+   
+   # WeibullLSS()
+   etamu <- 3 + 1*x1 +2*x2
+   etasi <- exp(rep(0.2, 1000))
+   status <- rep(1, 1000)
+   time <- rep(NA, 1000)
+   for (i in 1:1000)
+     time[i] <- rweibull(1, shape = exp(- 0.2), scale = exp(etamu[i]))
+   (m1 <- survreg(Surv(time, status) ~ x1 + x2 + x3, dist="weibull"))
+   model <- glmboostLSS(Surv(time, status) ~ x1 + x2 + x3,
+                        families = WeibullLSS(),
+                        control = boost_control(trace = TRUE, mstop = 300), center = TRUE)
+   stopifnot(sum(abs(coef(model, off2int = TRUE, which ="")[[1]] - c(3, 1, 2, 0)))
+             < sum(abs(coef(m1) - c(3, 1, 2, 0))))
+   stopifnot(sum(abs(coef(model, off2int = TRUE)[[2]] - c(0.2, 0, 0, 0))) < 0.4)
+ }
Loading required package: survival
Error in family@check_y(y) : 
  response is not an object of class 'Surv' but 'family = Lognormal()'
Error in family@check_y(y) : 
  response is not an object of class 'Surv' but 'families = LogNormalLSS()'
Error in family@check_y(y) : 
  response is not an object of class 'Surv' but 'family = LogNormalLSS()'
In addition: Warning message:
In mboostLSS_fit(formula = formula, data = data, families = families,  :
  responses differ for the components
[   1] ........................................ -- risk: 1552.695 
[  41] ........................................ -- risk: 1543.948 
[  81] ...................
Final risk: 1543.569 
[   1] ........................................ -- risk: 2246.758 
[  41] ........................................ -- risk: 2240.168 
[  81] ........................................ -- risk: 2237.48 
[ 121] ........................................ -- risk: 2236.279 
[ 161] ........................................ -- risk: 2235.629 
[ 201] ........................................ -- risk: 2235.268 
[ 241] ........................................ -- risk: 2235.067 
[ 281] ........................................ -- risk: 2234.955 
[ 321] .............................
Final risk: 2234.906 
[   1] ........................................ -- risk: 2134.384 
[  41] ........................................ -- risk: 1904.389 
[  81] ........................................ -- risk: 1826.087 
[ 121] ........................................ -- risk: 1790.255 
[ 161] ........................................ -- risk: 1777.757 
[ 201] ........................................ -- risk: 1774.481 
[ 241] ........................................ -- risk: 1773.758 
[ 281] ...................
Final risk: 1773.656 
>   
> if (require("DirechletReg")) {
+   set.seed(1907)
+   # Check Dirichlet family
+   n <- 150
+   p <- 10
+   
+   x <- matrix(runif(p * n, 0,1), n)
+   
+   x <- data.frame(x)
+   
+   a1 <- exp(2.5*x[,1] - x[,2] + 3*x[,3]) 
+   a2 <- exp(2*x[,4] + 2*x[,5] - x[,6])
+   a3 <- exp(1.5*x[,7] -  1.5*x[,8] + x[,9])
+   A <- cbind(a1,a2,a3)
+   
+   y <- DirichletReg::rdirichlet(nrow(A),A)
+   
+   colnames(y) <- c("y1","y2","y3")
+   
+   
+   # Check cyclical 
+   model <- glmboostLSS(y ~ ., data = x,
+                        families = DirichletLSS(K=3, stabilization = "none"),
+                        control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1))
+   model2 <- glmboostLSS(y ~ X1 + X2 + X3, data = x,
+                         families = DirichletLSS(K=3, stabilization = "none"),
+                         control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1))
+   # Check noncyclical
+   model3 <- glmboostLSS(y ~ ., data = x,
+                         families = DirichletLSS(K=3, stabilization = "none"),
+                         control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1), method = "noncyclic")
+   model4 <- glmboostLSS(y ~ X1 + X2 + X3, data = x,
+                         families = DirichletLSS(K=3, stabilization = "none"),
+                         control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1), method = "noncyclic")
+   
+   model[300]
+   model2[300]
+   model3[300]
+   model4[300]
+   coef(model, off2int = TRUE)
+   
+   # Check stabilization for Dirichlet family
+   model1.5 <- glmboostLSS(y ~ ., data = x,
+                           families = DirichletLSS(K=3, stabilization = "MAD"),
+                           control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1))
+   model2.5 <- glmboostLSS(y ~ X1 + X2 + X3, data = x,
+                           families = DirichletLSS(K=3, stabilization = "L2"),
+                           control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1))
+   model1.5[300]
+   model2.5[300]
+   coef(model1.5, off2int = TRUE)
+   
+   # Check gamboostLSS for Dirichlet family
+   x.train <- matrix(rnorm(p * n, 0,1), n)
+   x.train <- data.frame(x.train)
+   
+   a1.train <- exp(x.train[,1]**2) 
+   a2.train <- exp(2*tanh(3*x.train[,2]))
+   a3.train <- exp(cos(x.train[,3]))
+   
+   A <- cbind(a1.train,a2.train,a3.train)
+   
+   y.train <- DirichletReg::rdirichlet(nrow(A),A)
+   
+   x <- paste(c(paste("bbs(X", 1:p, ")", sep = "")), collapse = "+")
+   form <- as.formula((paste("y.train ~",  x)))
+   
+   model5 <- gamboostLSS(form, data = x.train,
+                         families = DirichletLSS(K = 3, stabilization = "none"),
+                         control = boost_control(trace = TRUE, mstop = 500, nu = 0.1), method = 'noncyclic')
+   model5[300]
+   coef(model5[200])
+   
+   # Check stabsel for Dirichlet family
+   modstabs <- stabsel(model, cutoff = 0.9, PFER = 5)
+   modstabs
+   
+   # Check for correct error message for Dirichlet family
+   try(glmboostLSS(y ~ ., data = x,
+                   families = DirichletLSS(),
+                   control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1)))
+ }
Loading required package: DirechletReg
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'DirechletReg'
> 
> ### Check that "families"-object contains a response function
> NBinomialMu2 <- function(...){
+     RET <- NBinomialMu(...)
+     RET@response <- function(f) NA
+     RET
+ }
> 
> NBinomialSigma2 <- function(...){
+     RET <- NBinomialSigma(...)
+     RET@response <- function(f) NA
+     RET
+ }
> 
> NBinomialLSS2 <- function(mu = NULL, sigma = NULL){
+     if ((!is.null(sigma) && sigma <= 0) || (!is.null(mu) && mu <= 0))
+         stop(sQuote("sigma"), " and ", sQuote("mu"),
+              " must be greater than zero")
+     RET <- Families(mu = NBinomialMu2(mu = mu, sigma = sigma),
+                         sigma = NBinomialSigma2(mu = mu, sigma = sigma))
+     return(RET)
+ }
> 
> try(NBinomialLSS2())
Error in Families(mu = NBinomialMu2(mu = mu, sigma = sigma), sigma = NBinomialSigma2(mu = mu,  : 
  response function not specified in families for:
	mu, sigma
> 
> #detach("package:gamboostLSS", unload = TRUE)
> 
> 
> ### Check that weighed median works well, which is needed if the negative
> ### gradient is stabilized using MAD
> set.seed(122)
> 
> ## some tests
> x <- c(1, 2, 3)
> stopifnot(weighted.median(x) == median(x))
> 
> ## this doesn't work a.t.m.
> w <- c(0, 1, 2)
> stopifnot(weighted.median(x, w) == median(rep(x, w)))
> 
> x <- c(1, 2, 3, 4)
> stopifnot(weighted.median(x) == median(x))
> 
> w <- c(0, 1, 2, 3)
> stopifnot(weighted.median(x, w) == median(rep(x, w)))
> 
> w <- rep(0:1, each = 50)
> x <- rnorm(100)
> stopifnot(weighted.median(x, w) == median(rep(x, w)))
> 
> w <- rep(0:1, each = 51)
> x <- rnorm(102)
> stopifnot(weighted.median(x, w) == median(rep(x, w)))
> 
> ## hope this is correct:
> w <- runif(100, 0, 1)
> x <- rnorm(100)
> (wm <- weighted.median(x, w))
[1] 0.0965101
> 
> ## check weighted.median with NAs.
> 
> # Missing in weights
> tmp <- w[2]
> w[2] <- NA
> stopifnot(is.na(weighted.median(x, w)))
> stopifnot(weighted.median(x, w, na.rm = TRUE) == wm)
> ## not always true but here it is
> 
> # Missing in x
> w[2] <- tmp
> x[5] <- NA
> stopifnot(is.na(weighted.median(x, w)))
> stopifnot(weighted.median(x, w, na.rm = TRUE) == wm)
> ## not always true but here it is
> 
> proc.time()
   user  system elapsed 
   9.45    0.53    9.96