bug_BB_BI.R

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(ymat ~ x + z, data = data, family = BI) # samesame

Betabin <- as.families(fname="BB")
gamboostLSS(ymat ~ x + z, data = data, families = Betabin)
# 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))

#This works:
coef(gamboostLSS(list(
  mu = y ~ bols(x) + bols(z),
  sigma = y ~ bols(ones, intercept=FALSE)),
  data = data, families = Betabin60))