> install.packages("gamboostLSS")
Installing package into **********************
(as ‘lib’ is unspecified)
Warning message:
package ‘gamboostLSS’ is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

txt <- "axis box coef dlogis dnorm fitted grey lines median optimize plogis
     plot pnorm points predict qlnorm qlogis qnorm qt qweibull rgb
     rmultinom tail update var weighted.mean"

imports_for_undefined_globals <- function(txt, lst, selective = TRUE) {
    if(!missing(txt))
        lst <- scan(what = character(), text = txt, quiet = TRUE)
    nms <- lapply(lst, find)
    ind <- sapply(nms, length) > 0L
    imp <- split(lst[ind], substring(unlist(nms[ind]), 9L))
    if(selective) {
        sprintf("importFrom(%s)",
                vapply(Map(c, names(imp), imp),
                       function(e)
                           paste0("\"", e, "\"", collapse = ", "),
                       ""))
    } else {
        sprintf("import(\"%s\")", names(imp))
    }
}

writeLines(imports_for_undefined_globals(txt))

## families object for new distribution
newStudentT <- Families(mu= newStudentTMu(mu=mu, df=df),
                        df=newStudentTDf(mu=mu, df=df))

#  Error in body(x@response) : 
#   trying to get slot "response" from an object of a basic class ("function") with no slots

> options(error=traceback, deparse.max.lines = 10)
> library("gamboostLSS")
> library("gamlss")
> 
> # y ~ Binomial(p = invlogit(x + z), n= 60)
> 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))
> (m_gamlss <- gamlss(ymat ~ x + z, data = data, family = BB))
GAMLSS-RS iteration 1: Global Deviance = 700.266 
GAMLSS-RS iteration 2: Global Deviance = 553.8168 
GAMLSS-RS iteration 3: Global Deviance = 532.8859 
GAMLSS-RS iteration 4: Global Deviance = 532.8656 
GAMLSS-RS iteration 5: Global Deviance = 532.8656 

Family:  c("BB", "Beta Binomial") 
Fitting method: RS() 

Call:  gamlss(formula = ymat ~ x + z, family = BB, data = data) 

Mu Coefficients:
(Intercept)            x            z  
  -0.005859     1.028064     0.994731  
Sigma Coefficients:
(Intercept)  
     -5.265  

 Degrees of Freedom for the fit: 4 Residual Deg. of Freedom   96 
Global Deviance:     532.866 
            AIC:     540.866 
            SBC:     551.286 
> # gamlss(ymat ~ x + z, data = data, family = BI) # samesame
> 
> Betabin <- as.families(fname="BB")
> gamboostLSS(ymat ~ x + z, data = data, families = Betabin)
Error in eval(expr, envir, enclos) : object 'bd' not found
No traceback available 
> # bd is the "binomial denominator" (i.e., # of trials), not being handed over!
> 
> #-------------------------------------------------------------------------------
> # ugly, hacky fix: set bd-defaults to 60 wherever it's needed
> Betabin60 <- Betabin
> add_bd_60 <- function(f){
+   frmls <- formals(f)
+   frmls$bd <- 60
+   formals(f) <- frmls
+   f
+ }
> with(environment(Betabin60[[1]]@ngradient), {
+   #all slots share same env, so enough to do this for one slot
+   FAM$dldm <- add_bd_60(FAM$dldm)
+   FAM$d2ldm2  <- add_bd_60(FAM$d2ldm2)
+   FAM$d2ldm2  <- add_bd_60(FAM$d2ldm2)
+   FAM$dldd  <- add_bd_60(FAM$dldd)
+   FAM$d2ldd2  <- add_bd_60(FAM$d2ldd2)
+   FAM$G.dev.incr <- add_bd_60(FAM$G.dev.incr)
+   FAM$mu.initial <- quote(mu <- (y + 0.5)/(60 + 1))
+   pdf <- add_bd_60(pdf)
+ })
> with(environment(Betabin60[[2]]@ngradient), {
+   FAM$dldm <- add_bd_60(FAM$dldm)
+   FAM$d2ldm2  <- add_bd_60(FAM$d2ldm2)
+   FAM$d2ldm2  <- add_bd_60(FAM$d2ldm2)
+   FAM$dldd  <- add_bd_60(FAM$dldd)
+   FAM$d2ldd2  <- add_bd_60(FAM$d2ldd2)
+   FAM$G.dev.incr <- add_bd_60(FAM$G.dev.incr)
+   FAM$mu.initial <- quote(mu <- (y + 0.5)/(60 + 1))
+   pdf <- add_bd_60(pdf)
+ })
> #-------------------------------------------------------------------------------
> data$ones <- rep(1, n)
> # matrix-response like gamlss produces multivariate model for success and failure
> # probabilities (I guess?!?), not model for sucess probabilities: 
> coef(gamboostLSS(list(
+   mu = ymat ~ bols(x) + bols(z), 
+   sigma = ymat ~ bols(ones, intercept=FALSE)),
+   data = data, families = Betabin60))
$mu
$mu$`bols(x)`
(Intercept)           x        <NA>        <NA> 
-0.20896247  0.97583503  0.02081331 -0.97472539 

$mu$`bols(z)`
  (Intercept)             z          <NA>          <NA> 
 0.0001546734  0.9402122029 -0.2000348943 -0.9390756321 

attr(,"offset")
[1] 0.1940493

$sigma
$sigma$`bols(ones, intercept = FALSE)`
     ones      <NA> 
-3.221820 -3.210209 

attr(,"offset")
[1] 0
>
> coef(gamboostLSS(list(
+   mu = y ~ bols(x) + bols(z),
+   sigma = y ~ bols(ones, intercept=FALSE)),
+   data = data, families = Betabin60))
$mu
$mu$`bols(x)`
(Intercept)           x 
 -0.2089625   0.9758350 

$mu$`bols(z)`
 (Intercept)            z 
0.0001546734 0.9402122029 

attr(,"offset")
[1] 0.1940493

$sigma
$sigma$`bols(ones, intercept = FALSE)`
    ones 
-3.22182 

attr(,"offset")
[1] 0

library(gamboostLSS)

# g_muvsi is a function that allows to see which distribution
# parameter is updated in each iteration.
# If mu is updated, then returns 1, if sigma is updated, then it returns 2.
# The first value in the risk vector is removed as it is the initial value.
g_muvsi <- function(gamboostLSS_object) {
  df1 <- data.frame(risk = risk(gamboostLSS_object)$mu[-1], muvsi = 1)
  df2 <- data.frame(risk = risk(gamboostLSS_object)$sigma[-1], muvsi = 2)
  df <- rbind(df1, df2)
  df <- df[order(df$risk, decreasing = TRUE), ]
  return(df$muvsi)
}

set.seed(1)
n <- 3
x1 <- c(1,2,3)

mu <- x1
sigma <- exp(.1 * x1)

data <- cbind(x1)
y <- rnorm(3, mu, sigma)

# center x1
cx1 <- scale(x1, center = T, scale = F)

b_g <- gamboostLSS(y ~ bols(cx1),
                   control = boost_control(mstop = 10),
                   method = "noncyclic")

> g_muvsi(b_g)
 [1] 1 2 2 1 2 2 1 2 2 1

mstop(b_g) = 2
ngr_mu <- (1/exp(predict(b_g)$sigma)^2) * (y - predict(b_g)$mu)
ngr_si <- -1 + exp(-2 * predict(b_g)$sigma) * (y - predict(b_g)$mu)^2

lm_mu <- lm(ngr_mu ~ cx1)
lm_si <- lm(ngr_si ~ cx1)

risk_mu <- -sum(dnorm(y, predict(b_g)$mu + .1 * lm_mu$fitted.values, exp(predict(b_g)$sigma), log = T))
risk_si <- -sum(dnorm(y, predict(b_g)$mu, exp(predict(b_g)$sigma + .1 * lm_si$fitted.values), log = T))

> risk_mu
[1] 3.515798
> risk_si
[1] 3.517315

> but_risk_mu <- -sum(dnorm(y, predict(b_g)$mu, exp(predict(b_g)$sigma), log = T))
> but_risk_mu
[1] 3.609269

but_risk_si <- -sum(dnorm(y, predict(b_g)$mu, exp(predict(b_g)$sigma), log = T))

Thomas, J., Mayr, A., Bischl, B., Schmid, M., Smith, A., and Hofner, B. (2016), 
Stability selection for component-wise gradient boosting in multiple dimensions. 
Preliminary version: \url{http://arxiv.org/abs/1611.10171}.

Thomas, J., Mayr, A., Bischl, B., Schmid, M., Smith, A., and Hofner, B. (2017), 
Gradient boosting for distributional regression - faster tuning and improved 
variable selection via noncyclical updates. 
Accepted for publication in \emph{Statistics and Computing}.\cr
Preliminary version: \url{http://arxiv.org/abs/1611.10171}.

library("gamboostLSS")
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")

## error:
predint(gb, which = "x1")

## update value of nuisance parameters
## use response(fitted()) as this is much quicker than fitted(, type = response)
 for (k in mods[-j])
     assign(names(fit)[k], families[[k]]@response(fitted(fit[[k]])),
                environment(get("ngradient", environment(fit[[j]]$subset))))

for (k in mods[-j]){
    if (!is.null(fit[[k]]))
          assign(names(fit)[k], fitted(fit[[k]], type = "response"),
               environment(families[[j]]@ngradient))
}

library(gamboostLSS)

### Data generating process:
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)

### linear model with y ~ . for both components: 400 boosting iterations
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
                     control = boost_control(mstop = 400),
                     center = TRUE, method = "noncyclic")

## works fine
summary(model)

## summary() breaks if for one component no base-learner was selected
summary(model[10])

## breaks as well 
summary(model[0])

### linear model with y ~ . for both components: 400 boosting iterations
model2 <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
                     control = boost_control(mstop = 400),
                     center = TRUE, method = "cyclic")

## works fine
summary(model2[0])

data$ones <- rep(1, n)
gamboostLSS(list( .... ~ ...,  bla = response ~ bols(ones, intercept =FALSE)), 
  data=data, families=SomethingCraycray())

 bla = response ~ 1

Error: processing vignette 'gamboostLSS_Tutorial.Rnw' failed with diagnostics:
 chunk 23 

Error in rep(NA, mmo - length(rsk)) : invalid 'times' argument
Execution halted

library(gamboostLSS)

### Data generating process
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)

### linear model with y ~ . for both components: 400 boosting iterations
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
                     control = boost_control(mstop = 400),
                     center = TRUE)


### summary just tells you that the model is a list
summary(model)


### summary function based on print.mboostLSS() and summary.mboost()
summary.mboostLSS <- function(object, ...) {

  ####cl <- match.call()
  cat("\n")
  cat("\t LSS Models fitted via Model-based Boosting\n")
  cat("\n")
  if (!is.null(attr(object, "call"))) 
    cat("Call:\n", deparse(attr(object, "call")), "\n\n", sep = "")
  cat("Number of boosting iterations: mstop =", mstop(object), "\n")
  cat("Step size: ", sapply(1:length(object), function(i) object[[i]]$control$nu), "\n\n")

  cat("Families:\n")
  lapply(object, function(xi){

    show(xi$family)

    nm <- variable.names(xi)
    selprob <- tabulate(selected(xi), nbins = length(nm)) / length(selected(xi))
    names(selprob) <- names(nm)
    selprob <- sort(selprob, decreasing = TRUE)
    selprob <- selprob[selprob > 0] 
    show(selprob) 

  } )
  invisible(object)
}


### summary of the model 
summary(model)

library(gamboostLSS)
library(partykit)

# Init boost control params
boost_ctrl <- boost_control(trace = TRUE)

# Init tree control params
tree_ctrl <- ctree_control()

# Fit blackboost model
bb <- blackboost(dist ~ speed, 
                      data = cars,
                      control = boost_ctrl,
                      tree_controls = tree_ctrl)

### plot fit
plot(dist ~ speed, data = cars, main = "Blackboost Fit")
lines(cars$speed, predict(bb), col = "red")

# Fit blackboostLSS model
bb_lss <- blackboostLSS(formula = list(mu = dist ~ speed,
                                           sigma = dist ~ speed),
                            data = cars,
                            method = "cyclic",
                            control = boost_ctrl,
                            tree_controls = tree_ctrl)

# plot fit
plot(dist ~ speed, data = cars, main = "BlackboostLSS fit")
lines(cars$speed, predict(bb_lss)$mu, col = "red")