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 families > > 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. > > > ### 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 > > 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)) Error in family@check_y(y) : response is not an object of class 'Surv' but 'family = Lognormal()' > try(glmboostLSS(time ~ x1 + x2 + x3, families = LogNormalLSS(), + control = boost_control(trace = TRUE), center = TRUE)) Error in family@check_y(y) : response is not an object of class 'Surv' but 'families = LogNormalLSS()' > require("survival") Loading required package: survival > try(glmboostLSS(list(mu = Surv(time, status) ~ x1 + x2 + x3, + sigma = time ~ x1 + x2 + x3), families = LogNormalLSS(), + control = boost_control(trace = TRUE), center = TRUE)) 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 > > ## check results > > # LogNormalLSS() > (m1 <- survreg(Surv(time, status) ~ x1 + x2 + x3, dist="lognormal")) Call: survreg(formula = Surv(time, status) ~ x1 + x2 + x3, dist = "lognormal") Coefficients: (Intercept) x1 x2 x3 2.88588933 1.15548098 2.03642866 0.09869623 Scale= 1.135837 Loglik(model)= -6051 Loglik(intercept only)= -6199.6 Chisq= 297.2 on 3 degrees of freedom, p= <2e-16 n= 1000 > model <- glmboostLSS(Surv(time, status) ~ x1 + x2 + x3, families = LogNormalLSS(), + control = boost_control(trace = TRUE), center = TRUE) [ 1] ........................................ -- risk: 1552.695 [ 41] ........................................ -- risk: 1543.948 [ 81] ................... Final risk: 1543.569 > 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")) Call: survreg(formula = Surv(time, status) ~ x1 + x2 + x3, dist = "loglogistic") Coefficients: (Intercept) x1 x2 x3 3.2109423 0.6543336 2.0799870 -0.2505811 Scale= 1.257854 Loglik(model)= -6681.6 Loglik(intercept only)= -6721.9 Chisq= 80.57 on 3 degrees of freedom, p= <2e-16 n= 1000 > model <- glmboostLSS(Surv(time, status) ~ x1 + x2 + x3, families = LogLogLSS(), + control = boost_control(trace = TRUE), center = TRUE) [ 1] ........................................ -- risk: 2246.758 [ 41] ........................................ -- risk: 2240.168 [ 81] ................... Final risk: 2238.548 > model[350] [ 101] ........................................ -- risk: 2236.777 [ 141] ........................................ -- risk: 2235.907 [ 181] ........................................ -- risk: 2235.422 [ 221] ........................................ -- risk: 2235.153 [ 261] ........................................ -- risk: 2235.003 [ 301] ........................................ -- risk: 2234.92 [ 341] ......... Final risk: 2234.906 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = Surv(time, status) ~ x1 + x2 + x3, families = LogLogLSS(), control = boost_control(trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 350, sigma = 350 Step size: mu = 0.1, sigma = 0.1 Families: Log-Logistic AFT Model: mu (id link) Loss function: { logfw <- function(pred) dlogis(pred, log = TRUE) logSw <- function(pred) plogis(pred, lower.tail = FALSE, log.p = TRUE) eta <- (log(y[, 1]) - f)/sigma -y[, 2] * (logfw(eta) - log(sigma)) - (1 - y[, 2]) * logSw(eta) } Log-Logistic AFT Model: sigma (log link) Loss function: { logfw <- function(pred) dlogis(pred, log = TRUE) logSw <- function(pred) pnorm(pred, lower.tail = FALSE, log.p = TRUE) eta <- (log(y[, 1]) - mu)/exp(f) -y[, 2] * (logfw(eta) - f) - (1 - y[, 2]) * logSw(eta) } > 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")) Call: survreg(formula = Surv(time, status) ~ x1 + x2 + x3, dist = "weibull") Coefficients: (Intercept) x1 x2 x3 3.047471670 1.008282745 1.934326789 -0.007107055 Scale= 1.230652 Loglik(model)= -5567.4 Loglik(intercept only)= -5684.9 Chisq= 235.11 on 3 degrees of freedom, p= <2e-16 n= 1000 > model <- glmboostLSS(Surv(time, status) ~ x1 + x2 + x3, + families = WeibullLSS(), + control = boost_control(trace = TRUE), center = TRUE) [ 1] ........................................ -- risk: 2134.384 [ 41] ........................................ -- risk: 1904.389 [ 81] ................... Final risk: 1858.757 > model[300] [ 101] ........................................ -- risk: 1804.158 [ 141] ........................................ -- risk: 1782.137 [ 181] ........................................ -- risk: 1775.546 [ 221] ........................................ -- risk: 1773.984 [ 261] ....................................... Final risk: 1773.656 LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = Surv(time, status) ~ x1 + x2 + x3, families = WeibullLSS(), control = boost_control(trace = TRUE), center = TRUE) Number of boosting iterations (mstop): mu = 300, sigma = 300 Step size: mu = 0.1, sigma = 0.1 Families: Weibull AFT Model: mu (id link) Loss function: { logfw <- function(pred) pred - exp(pred) logSw <- function(pred) -exp(pred) eta <- (log(y[, 1]) - f)/sigma -y[, 2] * (logfw(eta) - log(sigma)) - (1 - y[, 2]) * logSw(eta) } Weibull AFT Model: sigma (log link) Loss function: { logfw <- function(pred) pred - exp(pred) logSw <- function(pred) -exp(pred) eta <- (log(y[, 1]) - mu)/exp(f) -y[, 2] * (logfw(eta) - f) - (1 - y[, 2]) * logSw(eta) } > 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) > > # 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)) [ 1] ........................................ -- risk: -304.9326 [ 41] ........................................ -- risk: -370.0999 [ 81] ........................................ -- risk: -412.0619 [ 121] ........................................ -- risk: -429.0099 [ 161] ........................................ -- risk: -433.5379 [ 201] ........................................ -- risk: -434.8448 [ 241] ........................................ -- risk: -435.4784 [ 281] ........................................ -- risk: -435.8342 [ 321] ........................................ -- risk: -436.0625 [ 361] ........................................ -- risk: -436.2113 [ 401] ........................................ -- risk: -436.3145 [ 441] ........................................ -- risk: -436.3887 [ 481] ........................................ -- risk: -436.4417 [ 521] ........................................ -- risk: -436.48 [ 561] ........................................ -- risk: -436.5073 [ 601] ........................................ -- risk: -436.5273 [ 641] ........................................ -- risk: -436.5419 [ 681] ........................................ -- risk: -436.5527 [ 721] ........................................ -- risk: -436.5607 [ 761] ........................................ -- risk: -436.5666 [ 801] ........................................ -- risk: -436.571 [ 841] ........................................ -- risk: -436.5742 [ 881] ........................................ -- risk: -436.5767 [ 921] ........................................ -- risk: -436.5785 [ 961] ....................................... Final risk: -436.5798 > model2 <- glmboostLSS(y ~ X1 + X2 + X3, data = x, + families = DirichletLSS(K=3, stabilization = "none"), + control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1)) [ 1] ........................................ -- risk: -249.5704 [ 41] ........................................ -- risk: -263.2822 [ 81] ........................................ -- risk: -267.3187 [ 121] ........................................ -- risk: -268.709 [ 161] ........................................ -- risk: -269.2659 [ 201] ........................................ -- risk: -269.5121 [ 241] ........................................ -- risk: -269.6215 [ 281] ........................................ -- risk: -269.6699 [ 321] ........................................ -- risk: -269.691 [ 361] ........................................ -- risk: -269.7003 [ 401] ........................................ -- risk: -269.7043 [ 441] ........................................ -- risk: -269.706 [ 481] ........................................ -- risk: -269.7067 [ 521] ........................................ -- risk: -269.7071 [ 561] ........................................ -- risk: -269.7072 [ 601] ........................................ -- risk: -269.7073 [ 641] ........................................ -- risk: -269.7073 [ 681] ........................................ -- risk: -269.7073 [ 721] ........................................ -- risk: -269.7073 [ 761] ........................................ -- risk: -269.7073 [ 801] ........................................ -- risk: -269.7073 [ 841] ........................................ -- risk: -269.7073 [ 881] ........................................ -- risk: -269.7073 [ 921] ........................................ -- risk: -269.7073 [ 961] ....................................... Final risk: -269.7073 > # 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") [ 1] ........................................ -- risk: -244.9846 [ 41] ........................................ -- risk: -290.5706 [ 81] ........................................ -- risk: -322.7685 [ 121] ........................................ -- risk: -346.0673 [ 161] ........................................ -- risk: -366.1363 [ 201] ........................................ -- risk: -383.5644 [ 241] ........................................ -- risk: -398.0555 [ 281] ........................................ -- risk: -409.496 [ 321] ........................................ -- risk: -418.0534 [ 361] ........................................ -- risk: -424.0393 [ 401] ........................................ -- risk: -427.9449 [ 441] ........................................ -- risk: -430.5161 [ 481] ........................................ -- risk: -432.1258 [ 521] ........................................ -- risk: -433.1229 [ 561] ........................................ -- risk: -433.7447 [ 601] ........................................ -- risk: -434.1523 [ 641] ........................................ -- risk: -434.4898 [ 681] ........................................ -- risk: -434.7637 [ 721] ........................................ -- risk: -434.9875 [ 761] ........................................ -- risk: -435.1693 [ 801] ........................................ -- risk: -435.325 [ 841] ........................................ -- risk: -435.4596 [ 881] ........................................ -- risk: -435.5783 [ 921] ........................................ -- risk: -435.681 [ 961] ....................................... Final risk: -435.7718 > 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") [ 1] ........................................ -- risk: -232.4491 [ 41] ........................................ -- risk: -253.7881 [ 81] ........................................ -- risk: -262.7404 [ 121] ........................................ -- risk: -266.1183 [ 161] ........................................ -- risk: -267.6713 [ 201] ........................................ -- risk: -268.546 [ 241] ........................................ -- risk: -269.0508 [ 281] ........................................ -- risk: -269.341 [ 321] ........................................ -- risk: -269.5042 [ 361] ........................................ -- risk: -269.5947 [ 401] ........................................ -- risk: -269.645 [ 441] ........................................ -- risk: -269.6722 [ 481] ........................................ -- risk: -269.6872 [ 521] ........................................ -- risk: -269.6957 [ 561] ........................................ -- risk: -269.7006 [ 601] ........................................ -- risk: -269.7035 [ 641] ........................................ -- risk: -269.7051 [ 681] ........................................ -- risk: -269.706 [ 721] ........................................ -- risk: -269.7066 [ 761] ........................................ -- risk: -269.7069 [ 801] ........................................ -- risk: -269.7071 [ 841] ........................................ -- risk: -269.7072 [ 881] ........................................ -- risk: -269.7072 [ 921] ........................................ -- risk: -269.7073 [ 961] ....................................... Final risk: -269.7073 > > model[300] Model first reduced to mstop = 300. Now continue ... LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = x, families = DirichletLSS(K = 3, stabilization = "none"), control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1)) Number of boosting iterations (mstop): a1 = 300, a2 = 300, a3 = 300 Step size: a1 = 0.1, a2 = 0.1, a3 = 0.1 Families: Dirichlet Distribution: a1 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a2 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a3 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } > model2[300] Model first reduced to mstop = 300. Now continue ... LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ X1 + X2 + X3, data = x, families = DirichletLSS(K = 3, stabilization = "none"), control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1)) Number of boosting iterations (mstop): a1 = 300, a2 = 300, a3 = 300 Step size: a1 = 0.1, a2 = 0.1, a3 = 0.1 Families: Dirichlet Distribution: a1 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a2 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a3 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } > model3[300] LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = x, families = DirichletLSS(K = 3, stabilization = "none"), control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1), method = "noncyclic") Number of boosting iterations (mstop): a1 = 128, a2 = 102, a3 = 70 Step size: a1 = 0.1, a2 = 0.1, a3 = 0.1 Families: Dirichlet Distribution: a1 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a2 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a3 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } > model4[300] LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ X1 + X2 + X3, data = x, families = DirichletLSS(K = 3, stabilization = "none"), control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1), method = "noncyclic") Number of boosting iterations (mstop): a1 = 147, a2 = 80, a3 = 73 Step size: a1 = 0.1, a2 = 0.1, a3 = 0.1 Families: Dirichlet Distribution: a1 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a2 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a3 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } > coef(model, off2int = TRUE) $a1 (Intercept) X1 X2 X3 X4 X5 0.12366754 2.20742465 -0.84604037 2.72297063 0.11312073 -0.18841669 X6 X7 X9 X10 0.04669496 0.01709184 0.33719507 -0.03613969 $a2 (Intercept) X1 X2 X3 X4 X5 0.032057413 -0.205850458 0.077374223 -0.257496144 2.072967631 1.838452216 X6 X7 X8 X9 X10 -0.700765330 0.103210232 0.076433826 -0.001009495 0.139590009 $a3 (Intercept) X1 X3 X4 X5 X6 -0.20514066 -0.39745165 -0.25825827 0.26291390 0.12488929 -0.03068018 X7 X8 X9 X10 1.67682595 -1.23806523 1.19132323 0.09625697 > > # 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)) [ 1] ........................................ -- risk: -377.8566 [ 41] ........................................ -- risk: -421.7341 [ 81] ........................................ -- risk: -431.8549 [ 121] ........................................ -- risk: -434.2193 [ 161] ........................................ -- risk: -434.9821 [ 201] ........................................ -- risk: -435.4337 [ 241] ........................................ -- risk: -435.7481 [ 281] ........................................ -- risk: -435.9701 [ 321] ........................................ -- risk: -436.1271 [ 361] ........................................ -- risk: -436.2396 [ 401] ........................................ -- risk: -436.3203 [ 441] ........................................ -- risk: -436.3817 [ 481] ........................................ -- risk: -436.4283 [ 521] ........................................ -- risk: -436.4642 [ 561] ........................................ -- risk: -436.4915 [ 601] ........................................ -- risk: -436.5123 [ 641] ........................................ -- risk: -436.5282 [ 681] ........................................ -- risk: -436.5405 [ 721] ........................................ -- risk: -436.5499 [ 761] ........................................ -- risk: -436.5571 [ 801] ........................................ -- risk: -436.5627 [ 841] ........................................ -- risk: -436.5671 [ 881] ........................................ -- risk: -436.5706 [ 921] ........................................ -- risk: -436.5733 [ 961] ....................................... Final risk: -436.5755 > 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)) [ 1] ........................................ -- risk: -252.5769 [ 41] ........................................ -- risk: -264.885 [ 81] ........................................ -- risk: -267.8088 [ 121] ........................................ -- risk: -268.816 [ 161] ........................................ -- risk: -269.2599 [ 201] ........................................ -- risk: -269.4723 [ 241] ........................................ -- risk: -269.5808 [ 281] ........................................ -- risk: -269.6392 [ 321] ........................................ -- risk: -269.6706 [ 361] ........................................ -- risk: -269.6875 [ 401] ........................................ -- risk: -269.6967 [ 441] ........................................ -- risk: -269.7016 [ 481] ........................................ -- risk: -269.7042 [ 521] ........................................ -- risk: -269.7056 [ 561] ........................................ -- risk: -269.7064 [ 601] ........................................ -- risk: -269.7068 [ 641] ........................................ -- risk: -269.707 [ 681] ........................................ -- risk: -269.7072 [ 721] ........................................ -- risk: -269.7072 [ 761] ........................................ -- risk: -269.7073 [ 801] ........................................ -- risk: -269.7073 [ 841] ........................................ -- risk: -269.7073 [ 881] ........................................ -- risk: -269.7073 [ 921] ........................................ -- risk: -269.7073 [ 961] ....................................... Final risk: -269.7073 > model1.5[300] Model first reduced to mstop = 300. Now continue ... LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ ., data = x, families = DirichletLSS(K = 3, stabilization = "MAD"), control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1)) Number of boosting iterations (mstop): a1 = 300, a2 = 300, a3 = 300 Step size: a1 = 0.1, a2 = 0.1, a3 = 0.1 Families: Dirichlet Distribution: a1 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a2 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a3 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } > model2.5[300] Model first reduced to mstop = 300. Now continue ... LSS Models fitted via Model-based Boosting Call: glmboostLSS(formula = y ~ X1 + X2 + X3, data = x, families = DirichletLSS(K = 3, stabilization = "L2"), control = boost_control(trace = TRUE, mstop = 1000, nu = 0.1)) Number of boosting iterations (mstop): a1 = 300, a2 = 300, a3 = 300 Step size: a1 = 0.1, a2 = 0.1, a3 = 0.1 Families: Dirichlet Distribution: a1 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a2 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a3 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } > coef(model1.5, off2int = TRUE) $a1 (Intercept) X1 X2 X3 X4 X5 0.09152632 2.22354811 -0.85645460 2.73403444 0.14324511 -0.20357164 X6 X7 X8 X9 X10 0.03376635 0.04069848 -0.00956021 0.36748792 -0.02100004 $a2 (Intercept) X1 X2 X3 X4 X5 0.001937808 -0.188186611 0.068589947 -0.245830853 2.096956582 1.831306384 X6 X7 X8 X9 X10 -0.711626649 0.125916400 0.064535474 0.025587299 0.146601675 $a3 (Intercept) X1 X2 X3 X4 X5 -0.25192208 -0.38563083 -0.01309575 -0.24450920 0.29728795 0.13271531 X6 X7 X8 X9 X10 -0.04671728 1.71079410 -1.25154703 1.22421343 0.10922756 > > # 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') [ 1] ........................................ -- risk: -451.5869 [ 41] ........................................ -- risk: -487.0215 [ 81] ........................................ -- risk: -508.6955 [ 121] ........................................ -- risk: -524.4345 [ 161] ........................................ -- risk: -535.9484 [ 201] ........................................ -- risk: -544.85 [ 241] ........................................ -- risk: -551.7029 [ 281] ........................................ -- risk: -557.2145 [ 321] ........................................ -- risk: -561.7546 [ 361] ........................................ -- risk: -565.6992 [ 401] ........................................ -- risk: -569.1306 [ 441] ........................................ -- risk: -572.2762 [ 481] ................... Final risk: -573.7317 > model5[300] LSS Models fitted via Model-based Boosting Call: gamboostLSS(formula = form, data = x.train, families = DirichletLSS(K = 3, stabilization = "none"), control = boost_control(trace = TRUE, mstop = 500, nu = 0.1), method = "noncyclic") Number of boosting iterations (mstop): a1 = 107, a2 = 109, a3 = 84 Step size: a1 = 0.1, a2 = 0.1, a3 = 0.1 Families: Dirichlet Distribution: a1 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a2 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } Dirichlet Distribution: a3 (log link) Loss function: { A <- NULL eval(parse(text = paste0("A <- cbind(", gsub("f", "exp(f)", paste0(grep("a|f", arg, value = TRUE), collapse = ",")), ")"))) result <- -(lgamma(rowSums(A)) - rowSums(lgamma(A)) + rowSums((A - 1) * log(y))) return(result) } > coef(model5[200]) $a1 $a1$`bbs(X1)` 1 2 3 4 5 6 6.93942505 6.11814255 5.29559087 4.47299308 3.66291785 2.88205900 7 8 9 10 11 12 2.14932657 1.47700377 0.87500009 0.35711061 -0.05563997 -0.34517924 13 14 15 16 17 18 -0.48292012 -0.45667847 -0.27337305 0.05451148 0.50660324 1.06489009 19 20 21 22 23 24 1.70100473 2.38841065 3.09784156 3.81452808 4.53426114 5.25427648 $a1$`bbs(X3)` 1 2 3 4 5 6 0.042945893 0.038631150 0.034343361 0.030200817 0.026199147 0.022227895 7 8 9 10 11 12 0.018278344 0.014496224 0.011162264 0.008258213 0.005625062 0.003438298 13 14 15 16 17 18 0.001964153 0.001306919 0.001311758 0.002165746 0.004372004 0.007837411 19 20 21 22 23 24 0.012471688 0.018037023 0.024330616 0.031158822 0.038177189 0.045201186 $a1$`bbs(X6)` 1 2 3 4 5 6 -0.149174083 -0.127251084 -0.105308832 -0.083251081 -0.061071888 -0.039014734 7 8 9 10 11 12 -0.017445944 0.002639152 0.019154586 0.029992805 0.035630245 0.038108769 13 14 15 16 17 18 0.039580688 0.041937241 0.046429713 0.053331193 0.061722547 0.070400713 19 20 21 22 23 24 0.078412123 0.084864236 0.089534904 0.092694860 0.094893431 0.096926042 $a1$`bbs(X7)` 1 2 3 4 5 6 0.058384166 0.058211929 0.058056235 0.057626149 0.056284125 0.053767664 7 8 9 10 11 12 0.049868376 0.044478303 0.037970446 0.030290013 0.021238137 0.012990714 13 14 15 16 17 18 0.007706306 0.005292116 0.005136440 0.008384840 0.016773072 0.030294523 19 20 21 22 23 24 0.047770829 0.067964270 0.089768356 0.112511002 0.135714617 0.158992534 $a1$`bbs(X9)` 1 2 3 4 5 6 0.518567886 0.459813238 0.401134871 0.342914194 0.285608896 0.229694681 7 8 9 10 11 12 0.177567198 0.135092321 0.106667825 0.092252780 0.087453851 0.088509408 13 14 15 16 17 18 0.089554243 0.084870810 0.069631811 0.041959771 0.006294796 -0.030757083 19 20 21 22 23 24 -0.063934056 -0.090915353 -0.112113788 -0.129643969 -0.145664008 -0.161441873 attr(,"offset") [1] 0.0005403703 $a2 $a2$`bbs(X1)` 1 2 3 4 5 6 -0.247540007 -0.216563416 -0.185540315 -0.154604262 -0.124454540 -0.095959389 7 8 9 10 11 12 -0.069946548 -0.046835682 -0.026929024 -0.010727132 0.001914549 0.011716275 13 14 15 16 17 18 0.018970988 0.022411843 0.020236798 0.011519636 -0.003535872 -0.024483907 19 20 21 22 23 24 -0.050753452 -0.081281092 -0.114781952 -0.150067645 -0.186235629 -0.222549535 $a2$`bbs(X2)` 1 2 3 4 5 -2.2310762727 -2.1833091302 -2.1357986081 -2.0936488324 -2.0637928121 6 7 8 9 10 -2.0411512613 -2.0042055499 -1.9226613006 -1.7681339475 -1.5183620815 11 12 13 14 15 -1.1416398606 -0.6146186284 0.0004611164 0.5891757944 1.0608648820 16 17 18 19 20 1.3795665611 1.5592983435 1.6419584366 1.6612543615 1.6396666852 21 22 23 24 1.5993473828 1.5503670472 1.4961588703 1.4413656067 $a2$`bbs(X4)` 1 2 3 4 5 6 -0.091009761 -0.084073574 -0.077040006 -0.069422156 -0.060635744 -0.050106553 7 8 9 10 11 12 -0.037367575 -0.022881923 -0.008827343 0.003232943 0.012314355 0.017532799 13 14 15 16 17 18 0.018660205 0.017978364 0.018318129 0.020635252 0.024735711 0.030436399 19 20 21 22 23 24 0.037062470 0.044000491 0.051469635 0.059488174 0.067611508 0.075720867 $a2$`bbs(X6)` 1 2 3 4 5 6 0.35271896 0.28706907 0.22138881 0.15622556 0.09253621 0.03212959 7 8 9 10 11 12 -0.01891941 -0.05171425 -0.05448282 -0.01541501 0.06118288 0.15270441 13 14 15 16 17 18 0.23882289 0.30520827 0.34218560 0.34318334 0.30670168 0.23179666 19 20 21 22 23 24 0.11798839 -0.02879612 -0.20134338 -0.39318467 -0.59593384 -0.80053215 $a2$`bbs(X7)` 1 2 3 4 5 -0.0476902009 -0.0408250749 -0.0339226712 -0.0269449618 -0.0201117418 6 7 8 9 10 -0.0139177500 -0.0087093001 -0.0045191777 -0.0010152837 0.0024458706 11 12 13 14 15 0.0058620667 0.0089273463 0.0113087490 0.0132752646 0.0148290054 16 17 18 19 20 0.0156804303 0.0154513957 0.0137769521 0.0105774902 0.0061127269 21 22 23 24 0.0007367916 -0.0052937303 -0.0116014843 -0.0179302923 $a2$`bbs(X10)` 1 2 3 4 5 6 -0.064884028 -0.056175411 -0.047457541 -0.039020042 -0.031196953 -0.023852029 7 8 9 10 11 12 -0.017060714 -0.010222937 -0.002483168 0.006508571 0.016315613 0.025923445 13 14 15 16 17 18 0.034754196 0.042494467 0.048484179 0.052395984 0.054844050 0.057017275 19 20 21 22 23 24 0.060316776 0.064900652 0.070374450 0.076391409 0.082677085 0.089007547 attr(,"offset") [1] 6.203405e-11 $a3 $a3$`bbs(X1)` 1 2 3 4 5 6 -0.845818595 -0.728883180 -0.612208524 -0.497486145 -0.387319588 -0.284539512 7 8 9 10 11 12 -0.192020121 -0.112170958 -0.046649657 0.003812507 0.038235056 0.055273703 13 14 15 16 17 18 0.054436652 0.037577822 0.007590207 -0.032983809 -0.082081311 -0.137273684 19 20 21 22 23 24 -0.195738635 -0.256248694 -0.317819909 -0.379773018 -0.441888435 -0.504016663 $a3$`bbs(X2)` 1 2 3 4 5 6 7 0.45261891 0.41905018 0.38573436 0.35344145 0.32267265 0.29211431 0.26089187 8 9 10 11 12 13 14 0.23040676 0.20296272 0.17638594 0.14706962 0.11144120 0.07049690 0.03302735 15 16 17 18 19 20 21 0.01264193 0.02022673 0.05704520 0.11491116 0.18056150 0.24559726 0.30713910 22 23 24 0.36488260 0.42071138 0.47635826 $a3$`bbs(X3)` 1 2 3 4 5 6 -1.003306565 -0.791765931 -0.580094379 -0.367341476 -0.157769501 0.041520671 7 8 9 10 11 12 0.223859000 0.377876516 0.491481069 0.561039454 0.592714071 0.586029049 13 14 15 16 17 18 0.536879864 0.447360410 0.328743122 0.184581002 0.009973509 -0.191347044 19 20 21 22 23 24 -0.410938657 -0.637362516 -0.863596844 -1.089375902 -1.314666306 -1.539845562 $a3$`bbs(X4)` 1 2 3 4 5 6 0.041621958 0.034572483 0.027621427 0.021260884 0.016081365 0.012639737 7 8 9 10 11 12 0.011169542 0.011584552 0.013703797 0.017452851 0.023641353 0.032647600 13 14 15 16 17 18 0.043337892 0.052258542 0.056325702 0.055511964 0.051464490 0.043814266 19 20 21 22 23 24 0.031432081 0.014173866 -0.007005386 -0.030959572 -0.056159689 -0.081548266 $a3$`bbs(X6)` 1 2 3 4 5 6 0.222170865 0.197564465 0.173022203 0.148885306 0.125895803 0.105387903 7 8 9 10 11 12 0.087898773 0.073712534 0.064380628 0.060950706 0.061518836 0.064371123 13 14 15 16 17 18 0.066750838 0.064714500 0.055810210 0.040188207 0.019097967 -0.006063992 19 20 21 22 23 24 -0.033246327 -0.062115055 -0.092400173 -0.122308526 -0.151023369 -0.179512766 $a3$`bbs(X8)` 1 2 3 4 5 -0.0669629004 -0.0590703224 -0.0510803053 -0.0425133154 -0.0328377025 6 7 8 9 10 -0.0218467017 -0.0108004350 -0.0009695191 0.0095490748 0.0226456458 11 12 13 14 15 0.0374657861 0.0514899332 0.0623849633 0.0691159767 0.0712861154 16 17 18 19 20 0.0693439069 0.0639149377 0.0555246589 0.0453104579 0.0345529862 21 22 23 24 0.0244792102 0.0153864610 0.0069557999 -0.0013592738 $a3$`bbs(X9)` 1 2 3 4 5 6 -0.128996637 -0.109181692 -0.089432189 -0.070075336 -0.051503781 -0.034112388 7 8 9 10 11 12 -0.018531836 -0.005811518 0.003289367 0.009201205 0.013475504 0.017211712 13 14 15 16 17 18 0.020736849 0.023785441 0.025762492 0.028176007 0.032420424 0.039228050 19 20 21 22 23 24 0.049097180 0.061347426 0.074675655 0.088194626 0.101624484 0.115034395 attr(,"offset") [1] 4.223816e-05 > > # Check stabsel for Dirichlet family > modstabs <- stabsel(model, cutoff = 0.9, PFER = 5) Run stabsel .................................................................................................... > modstabs Stability Selection with unimodality assumption Selected variables: (Intercept).a1 X1.a1 X2.a1 X3.a1 (Intercept).a2 1 2 3 4 12 X4.a2 X5.a2 X6.a2 (Intercept).a3 X1.a3 16 17 18 23 24 X3.a3 X7.a3 X8.a3 X9.a3 26 30 31 32 Selection probabilities: X9.a1 X8.a1 X6.a1 X7.a1 X6.a3 0.02 0.04 0.06 0.06 0.06 X4.a1 X10.a1 X7.a2 X4.a3 X5.a3 0.08 0.09 0.11 0.12 0.21 X8.a2 X10.a3 X10.a2 X5.a1 X2.a2 0.25 0.27 0.33 0.38 0.47 X9.a2 X2.a3 X1.a2 X3.a2 X3.a3 0.51 0.55 0.77 0.83 0.91 X6.a2 X2.a1 X1.a3 X9.a3 (Intercept).a1 0.96 0.98 0.99 0.99 1.00 X1.a1 X3.a1 (Intercept).a2 X4.a2 X5.a2 1.00 1.00 1.00 1.00 1.00 (Intercept).a3 X7.a3 X8.a3 1.00 1.00 1.00 --- Cutoff: 0.9; q: 19; PFER (*): 4.72 (*) or expected number of low selection probability variables PFER (specified upper bound): 5 PFER corresponds to signif. level 0.143 (without multiplicity adjustment) > > # 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))) Error in DirichletLSS() : Number of categories 'K' must be specified. > > ### 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 67.56 2.43 69.98