require("gamboostLSS")
###negbin dist, linear###
set.seed(2611)
x1 <- rnorm(1000)
x2 <- rnorm(1000)
x3 <- rnorm(1000)
x4 <- rnorm(1000)
x5 <- rnorm(1000)
x6 <- rnorm(1000)
mu <- exp(1.5 + x1^2 +0.5 * x2 - 3 * sin(x3) -1 * x4)
sigma <- exp(-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)
#fit models at number of params + 1
#glmboost
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 3), method = "noncyclic")
#linear baselearner with bols
model <- gamboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 3), method = "noncyclic",
baselearner = "bols")
#nonlinear bbs baselearner
model <- gamboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 3), method = "noncyclic",
baselearner = "bbs")
#reducing model and increasing it afterwards should yield the same fit
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 50), method = "noncyclic")
m_co <- coef(model)
mstop(model) <- 5
mstop(model) <- 50
stopifnot(all.equal(m_co, coef(model)))
model <- gamboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 50), method = "noncyclic",
baselearner = "bols")
m_co <- coef(model)
mstop(model) <- 5
mstop(model) <- 50
stopifnot(all.equal(m_co, coef(model)))
model <- gamboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 50), method = "noncyclic",
baselearner = "bbs")
m_co <- coef(model)
mstop(model) <- 5
mstop(model) <- 50
stopifnot(all.equal(m_co, coef(model)))
model <- gamboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 50), method = "noncyclic",
baselearner = "bbs")
m_co <- coef(model)
mstop(model) <- 5
mstop(model) <- 50
stopifnot(all.equal(m_co, coef(model)))
## check cvrisk for noncyclic models
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 3), method = "noncyclic")
cvr1 <- cvrisk(model, grid = 1:50, cv(model.weights(model), B = 5))
cvr1
plot(cvr1)
risk(model, merge = TRUE)
risk(model, merge = FALSE)
## test that mstop = 0 is possible
compare_models <- function (m1, m2) {
stopifnot(all.equal(coef(m1), coef(m2), check.environment=FALSE))
stopifnot(all.equal(predict(m1), predict(m2), check.environment=FALSE))
stopifnot(all.equal(fitted(m1), fitted(m2), check.environment=FALSE))
stopifnot(all.equal(selected(m1), selected(m2), check.environment=FALSE))
stopifnot(all.equal(risk(m1), risk(m2), check.environment=FALSE))
## remove obvious differences from objects
m1$control <- m2$control <- NULL
m1$call <- m2$call <- NULL
if (!all.equal(m1, m2, check.environment=FALSE))
stop("Objects of offset model + 1 step and model with 1 step not identical")
invisible(NULL)
}
# set up models
mod <- glmboostLSS(y ~ ., data = dat, method = "noncyclic", control = boost_control(mstop = 0))
mod2 <- glmboostLSS(y ~ ., data = dat, method = "noncyclic", control = boost_control(mstop = 1))
mod3 <- glmboostLSS(y ~ ., data = dat, method = "noncyclic", control = boost_control(mstop = 1))
lapply(coef(mod), function(x) stopifnot(is.null(x)))
mstop(mod3) <- 0
mapply(compare_models, m1 = mod, m2 = mod3)
mstop(mod) <- 1
mapply(compare_models, m1 = mod, m2 = mod2)
mstop(mod3) <- 1
mapply(compare_models, m1 = mod, m2 = mod3)
## check selected
set.seed(1907)
x1 <- rnorm(500)
x2 <- rnorm(500)
x3 <- rnorm(500)
x4 <- rnorm(500)
x5 <- rnorm(500)
x6 <- rnorm(500)
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(500)
for( i in 1:500)
y[i] <- rnbinom(1, size = sigma[i], mu = mu[i])
dat <- data.frame(x1, x2, x3, x4, x5, x6, y)
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 10),
center = TRUE, method = "cyclic")
selected(model) # ok (at least in principle)
selected(model, merge = TRUE) # ok
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 10),
center = TRUE, method = "noncyclic")
selected(model) # ok (at least in principle)
selected(model, merge = TRUE) ## BROKEN
## with informative sigma:
sigma <- exp(-0.4 * x3 -0.2 * x4 +0.2 * x5 + 1 * x6)
y <- numeric(500)
for( i in 1:500)
y[i] <- rnbinom(1, size = sigma[i], mu = mu[i])
dat <- data.frame(x1, x2, x3, x4, x5, x6, y)
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 20),
center = TRUE, method = "cyclic")
selected(model) # ok (at least in principle)
selected(model, merge = TRUE) # ok
model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat,
control = boost_control(mstop = 20),
center = TRUE, method = "noncyclic")
selected(model) # ok (at least in principle)
selected(model, merge = TRUE) ## BROKEN
## Check merged risk for reducing mstop to 0, and increasing it again does not contain an NA
stopifnot(all(!is.na(risk(model, merge = TRUE))))
mstop(model) <- 0
stopifnot(all(!is.na(risk(model, merge = TRUE))))
mstop(model) <- 10
stopifnot(all(!is.na(risk(model, merge = TRUE))))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.