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.

> ###
> # check glmboostLSS()
> 
> 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.

> 
> 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)
> 
> ### check subset method
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(),
+                      control = boost_control(mstop = 10),
+                      center = TRUE)
> model2 <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(),
+                           control = boost_control(mstop = 20),
+                           center = TRUE)
> model[20]

	 LSS Models fitted via Model-based Boosting

Call:
glmboostLSS(formula = y ~ x1 + x2, families = StudentTLSS(),     control = boost_control(mstop = 10), center = TRUE)

Number of boosting iterations (mstop):  mu = 20, sigma = 20, df = 20 
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))) 
 } 
 
> stopifnot(all.equal(coef(model),coef(model2)))
> stopifnot(length(coef(model2, aggregate = "none")[[1]][[1]]) ==
+           length(coef(model, aggregate = "none")[[1]][[1]]))
> stopifnot(length(coef(model2, aggregate = "none")[[1]][[1]]) == 20)
> 
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(),
+                      control = boost_control(mstop = 10),
+                      center = TRUE)
> model2[10]

	 LSS Models fitted via Model-based Boosting

Call:
glmboostLSS(formula = y ~ x1 + x2, families = StudentTLSS(),     control = boost_control(mstop = 20), center = TRUE)

Number of boosting iterations (mstop):  mu = 10, sigma = 10, df = 10 
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))) 
 } 
 
> stopifnot(all.equal(coef(model),coef(model2)))
> stopifnot(length(coef(model2, aggregate = "none")[[1]][[1]]) ==
+           length(coef(model, aggregate = "none")[[1]][[1]]))
> stopifnot(length(coef(model2, aggregate = "none")[[1]][[1]]) == 10)
> 
> ### check trace
> model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(),
+                      control = boost_control(mstop = 10, trace =TRUE),
+                      center = TRUE)
[  1] .........
Final risk: 13591.79 
> model[100]
[ 11] ........................................ -- risk: 12681.81 
[ 51] ........................................ -- risk: 12508.02 
[ 91] .........
Final risk: 12491.57 

	 LSS Models fitted via Model-based Boosting

Call:
glmboostLSS(formula = y ~ x1 + x2, families = StudentTLSS(),     control = boost_control(mstop = 10, trace = TRUE), 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))) 
 } 
 
> 
> ### check formula-interface with lists
> set.seed(1907)
> n <- 5000
> x1  <- runif(n)
> x2 <- runif(n)
> mu <- 2 - 3*x2
> sigma <- exp(-1*x1 + 3*x2)
> df <- exp(1 + 3*x1)
> y <- rTF(n = n, mu = mu, sigma = sigma, nu = df)
> model <- glmboostLSS(list(mu = y ~ x2,
+                           sigma = y ~ x1 + x2,
+                           df = y ~ x1),
+                      families = StudentTLSS(),
+                      control = boost_control(mstop = 10, trace =TRUE),
+                      center = TRUE)
[  1] .........
Final risk: 13810.72 
> 
> stopifnot(all.equal(lapply(coef(model, which = ""), function(x) names(x)[-1]),
+                     list(mu = "x2", sigma = c("x1", "x2"), df = "x1")))
> 
> model <- glmboostLSS(list(mu = y ~ x2,
+                           df = y ~ x1,
+                           sigma = y ~ x1 + x2),
+                      families = StudentTLSS(),
+                      control = boost_control(mstop = 10, trace =TRUE),
+                      center = TRUE)
[  1] .........
Final risk: 13810.72 
> 
> stopifnot(all.equal(lapply(coef(model, which = ""), function(x) names(x)[-1]),
+                     list(mu = "x2", sigma = c("x1", "x2"), df = "x1")))
> 
> ### check formula-interface with lists and different responses
> ### (not really sensible with the current families)
> y2 <- y + 1
> model2 <- glmboostLSS(list(mu = y ~ x2,
+                            sigma = y ~ x1 + x2,
+                            df = y2 ~ x1),
+                       families = StudentTLSS(),
+                       control = boost_control(mstop = 10, trace =TRUE),
+                       center = TRUE)
[  1] .........
Final risk: 14135.24 
Warning message:
In mboostLSS_fit(formula = formula, data = data, families = families,  :
  responses differ for the components
> sapply(model, function(comps) comps$offset)
      mu    sigma       df 
1.160244 1.177889 0.765566 
> sapply(model2, function(comps) comps$offset)
       mu     sigma        df 
1.1602444 1.1778894 0.8176417 
> stopifnot((model2[[1]]$offset - model[[1]]$offset) < sqrt(.Machine$double.eps))
> stopifnot((model2[[2]]$offset - model[[2]]$offset) < sqrt(.Machine$double.eps))
> stopifnot(model2[[3]]$offset != model[[3]]$offset)
> 
> ### even better check for offset-issue when different responses are used
> 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=1, n=length(mu))
> y2 <- rnorm(mean=0, sd=sigma, n=length(sigma))
> dat <- data.frame(x1, x2, x3, x4, y, y2)
> model <- glmboostLSS(formula=list(mu=y~x1+x2, sigma=y2~x3+x4),
+                      families=GaussianLSS(), data=dat)
Warning message:
In mboostLSS_fit(formula = formula, data = data, families = families,  :
  responses differ for the components
> ## offset for mu must be equal to the mean of y per default
> stopifnot((model$mu$offset - mean(y)) < sqrt(.Machine$double.eps))
> ## offset for sigma must be equal to the log of the standard deviation of y2 per
> ## default
> stopifnot((model$sigma$offset - log(sd(y2))) < sqrt(.Machine$double.eps))
> 
> 
> ### check weights interface
> set.seed(1907)
> n <- 2500
> x1  <- runif(n)
> x2 <- runif(n)
> mu <- 2 - 3*x2
> sigma <- exp(-1*x1 + 3*x2)
> df <- exp(1 + 3*x1)
> y <- rTF(n = n, mu = mu, sigma = sigma, nu = df)
> dat <- data.frame(x1, x2, y)
> dat2 <- rbind(dat, dat) # data frame with duplicate entries
> 
> model <- glmboostLSS(list(mu = y ~ x2,
+                           sigma = y ~ x1 + x2,
+                           df = y ~ x1),
+                      data = dat, weights = rep(2, nrow(dat)),
+                      families = StudentTLSS(),
+                      control = boost_control(mstop = 10, trace =TRUE),
+                      center = TRUE)
[  1] .........
Final risk: 13810.99 
> 
> model2 <- glmboostLSS(list(mu = y ~ x2,
+                            sigma = y ~ x1 + x2,
+                            df = y ~ x1),
+                       data = dat2, families = StudentTLSS(),
+                       control = boost_control(mstop = 10, trace =TRUE),
+                       center = TRUE)
[  1] .........
Final risk: 13810.99 
> 
> stopifnot(all.equal(coef(model), coef(model2)))
> 
> proc.time()
   user  system elapsed 
   6.46    0.29    6.76