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 multivariate mstop > > #detach("package:gamboostLSS", unload = TRUE) > 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 > > ### create some data first > 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) > > > ## check if model with different mstops is what we expect > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = 10), + center = TRUE) > mstop(model) mu sigma 10 10 > model2 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = list(mu = 10, sigma = 20)), + center = TRUE) > mstop(model2) mu sigma 10 20 > > f1 <- fitted(model, parameter = "mu", type = "response") > f2 <- fitted(model, parameter = "sigma", type = "response") > model3 <- glmboost(y ~ ., family = NBinomialSigma(mu = f1, sigma = f2, + stabilization = "none"), + data = dat, + control = boost_control(mstop = 10), + center = TRUE) > > tmp <- coef(model3) + coef(model)$sigma > stopifnot(max(abs(tmp - coef(model2)$sigma)) < sqrt(.Machine$double.eps)) > > # Plot > layout(matrix(c(1:4, 6, 5), byrow = TRUE, ncol = 2)) > plot(model, xlim = c(0,20), ylim = range(sapply(coef(model2), range))) > plot(model2, xlim = c(0, 20), ylim = range(sapply(coef(model2), range))) > plot(model, xlim = c(0,20), ylim = range(sapply(coef(model2), range)), + parameter = "sigma") > cp <- coef(model3, aggregate = "cumsum") > cp <- matrix(unlist(cp), nrow = length(cp), byrow = TRUE) > cp <- cp + coef(model)$sigma > cp <- cbind(coef(model)$sigma, cp) > cf <- cp[, ncol(cp)] > col <- hcl(h = 40, l = 50, c = abs(cf)/max(abs(cf)) * 490) > matlines(10:20, t(cp), type = "l", xlab = "Number of boosting iterations", + ylab = "Coefficients", col = col) > > > > ### check subset method > ms <- list(mu = 10, sigma = 20) > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ................... Final risk: 3197.832 > > model[c(20, 30)] # check if two values can be specified Model first reduced to mstop = 10. Now continue ... [ 11] ................... Final risk: 3150.898 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(), control = boost_control(mstop = ms, trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 20, sigma = 30 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))) > ms <- list(mu = 20, sigma = 30) > modela <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ............................. Final risk: 3150.898 > stopifnot(max(abs(coef(model)[[1]] - coef(modela)[[1]])) + < sqrt(.Machine$double.eps)) > stopifnot(max(abs(coef(model)[[2]] - coef(modela)[[2]])) + < sqrt(.Machine$double.eps)) > > model[40] # check if one value can be specified Model first reduced to mstop = 20. Now continue ... [ 21] ................... Final risk: 3098.51 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(), control = boost_control(mstop = ms, trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 40, sigma = 40 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))) > mstop(model) mu sigma 40 40 > modelb <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = 40, trace = TRUE), + center = TRUE) [ 1] ....................................... Final risk: 3098.51 > stopifnot(all.equal(risk(model), risk(modelb))) > > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = 10, trace = TRUE), + center = TRUE) [ 1] ......... Final risk: 3209.164 > model[20] [ 11] ......... Final risk: 3162.029 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(), control = boost_control(mstop = 10, trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 20, sigma = 20 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))) > model2 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = 20, trace = TRUE), + center = TRUE) [ 1] ................... Final risk: 3162.029 > stopifnot(all.equal(risk(model), risk(model2))) > > ms <- list(mu = 10, sigma = 20) > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ................... Final risk: 3197.832 > model[c(5,10)] Model first reduced to mstop = 5. Now continue ... [ 6] .... Final risk: 3232.999 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(), control = boost_control(mstop = ms, trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 5, 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))) > ms <- list(mu = 5, sigma = 10) > model2 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ......... Final risk: 3232.999 > stopifnot(all.equal(risk(model), risk(model2))) > > > ### check subset method where only bigger mstop-value is touched > # increase model > ms <- list(mu = 10, sigma = 20) > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ................... Final risk: 3197.832 > model[c(10,25)] [21] .... Final risk: 3193.834 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(), control = boost_control(mstop = ms, trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 10, sigma = 25 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))) > mstop(model) mu sigma 10 25 > ms <- list(mu = 10, sigma = 25) > model2 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ........................ Final risk: 3193.834 > stopifnot(all.equal(risk(model), risk(model2))) > > # reduce model > ms <- list(mu = 10, sigma = 20) > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ................... Final risk: 3197.832 > model[c(10,15)] Model first reduced to mstop = 15. Now continue ... LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(), control = boost_control(mstop = ms, trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 10, sigma = 15 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))) > mstop(model) mu sigma 10 15 > ms <- list(mu = 10, sigma = 15) > model2 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] .............. Final risk: 3202.876 > stopifnot(all.equal(risk(model), risk(model2))) > > # reduce such that mu needs to be touced again > ms <- list(mu = 10, sigma = 20) > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ................... Final risk: 3197.832 > model[c(10,9)] Model first reduced to mstop = 9. Now continue ... Final risk: 3210.562 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = dat, families = NBinomialLSS(), control = boost_control(mstop = ms, trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 10, sigma = 9 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))) > mstop(model) mu sigma 10 9 > ms <- list(mu = 10, sigma = 9) > model2 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = ms, trace = TRUE), + center = TRUE) [ 1] ......... Final risk: 3210.562 > stopifnot(all.equal(risk(model), risk(model2))) > > ### check multiple values of nu > > nus <- list(mu = 0, sigma = 0.2) > model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, + control = boost_control(mstop = 10, nu = nus, trace = TRUE), + center = TRUE) [ 1] ......... Final risk: 3261.092 > stopifnot(all(coef(model)[[1]] == 0)) > stopifnot(any(coef(model)[[2]] != 0)) > > proc.time() user system elapsed 2.90 0.18 3.07