diff --git a/.travis.yml b/.travis.yml index 013112c..721c4d4 100644 --- a/.travis.yml +++ b/.travis.yml @@ -9,6 +9,7 @@ language: r sudo: required +dist: trusty ## change directory before installation ## as R packages are not available elsewise diff --git a/DESCRIPTION b/DESCRIPTION index 6ab2749..a64d718 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,16 +1,16 @@ Package: gamboostLSS Type: Package Title: Boosting Methods for 'GAMLSS' -Version: 1.2-2 -Date: 2016-10-18 -Author: Benjamin Hofner, Andreas Mayr, Nora Fenske, Matthias Schmid +Version: 2.0-0 +Date: 2016-xx-yy +Author: Benjamin Hofner, Andreas Mayr, Nora Fenske, Janek Thomas, Matthias Schmid Maintainer: Benjamin Hofner Description: Boosting models for fitting generalized additive models for location, shape and scale ('GAMLSS') to potentially high dimensional data. -Depends: R (>= 2.10.0), mboost (>= 2.3-0), parallel +Depends: R (>= 2.10.0), mboost (>= 2.3-0), stabs (>= 0.5-0), parallel Imports: graphics, grDevices, stats, utils -Suggests: gamlss, BayesX, gamlss.dist, survival, R2BayesX +Suggests: gamlss, gamlss.dist, survival, BayesX, R2BayesX LazyLoad: yes LazyData: yes License: GPL-2 diff --git a/NAMESPACE b/NAMESPACE index 004beea..50dee1b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,5 @@ import("mboost") +importFrom("stabs", "stabsel", "subsample", "run_stabsel") import("parallel") importFrom("graphics", "axis", "box", "lines", "plot", "points") @@ -16,6 +17,8 @@ export(glmboostLSS, model.weights, PI, predint, plot.predint, make.grid, cvrisk.mboostLSS, + selected, selected.mboostLSS, selected.stabsel_mboostLSS, + stabsel.mboostLSS, Families, NBinomialLSS, NBinomialMu, NBinomialSigma, StudentTLSS, StudentTMu, StudentTSigma, StudentTDf, @@ -37,6 +40,7 @@ S3method("[", mboostLSS) S3method(print, mboostLSS) S3method(selected, mboostLSS) S3method(risk, mboostLSS) +S3method(risk, nc_mboostLSS) S3method(summary, mboostLSS) S3method(mstop, oobag) @@ -46,10 +50,15 @@ S3method(model.weights, default) S3method(model.weights, mboostLSS) S3method(cvrisk, mboostLSS) +S3method(cvrisk, nc_mboostLSS) S3method(print, cvriskLSS) S3method(plot, cvriskLSS) +S3method(plot, nc_cvriskLSS) S3method(mstop, cvriskLSS) +S3method(stabsel, mboostLSS) +S3method(selected, stabsel_mboostLSS) + S3method(coef, mboostLSS) S3method(coef, glmboostLSS) diff --git a/R/AAA.R b/R/AAA.R index 2692bd9..705b537 100644 --- a/R/AAA.R +++ b/R/AAA.R @@ -3,3 +3,6 @@ ### the variance of the negative gradient options(gamboostLSS_stab_ngrad = FALSE) } + +# get rid of NOTEs in R CMD check for "undefined global functions or variables" +globalVariables(c("ngradient", "y", "combined_risk")) diff --git a/R/as.families.R b/R/as.families.R index f490da1..d8355bd 100755 --- a/R/as.families.R +++ b/R/as.families.R @@ -1,505 +1,649 @@ - -################################################################################ -### Family wrapper for gamlss families - -## to use stability correction for gradients -## please set -## options(gamboostLSS_stab_ngrad = TRUE) - - - -################################################################################ -## constructor - -## a wrapper to as.families (for compatibility reasons -gamlss.Families <- function(...) - as.families(...) - -as.families <- function(fname = "NO", - mu = NULL, sigma = NULL, nu = NULL, tau = NULL, - stabilization = c("none", "MAD")) { - - ## require gamlss.dist - if (!requireNamespace("gamlss.dist", quietly = TRUE)) - stop("Please install package 'gamlss.dist' for using gamlss families.") - - if (is.function(fname)) - fname <- fname() - - if (inherits(fname, "gamlss.family")) - fname <- fname$family[1] - - if (mode(fname) != "character" && mode(fname) != "name") - fname <- as.character(substitute(fname)) - - gamlss.fam <- try(gamlss.dist::gamlss.family(fname), silent = TRUE) - if (inherits(gamlss.fam, "try-error")) - stop(sQuote("fname"), " specifies no valid gamlss family") - - stabilization <- check_stabilization(stabilization) - - npar <- gamlss.fam$nopar - switch(npar, { - ## 1 parameter - fun <- gamlss1parMu(mu = mu, fname = fname) - warning("For boosting one-parametric families,", - " please use the mboost package.") - if (stabilization != "none") - warning("Stabilization is ignored for one-parametric families.") - }, { - ## 2 parameters - fun <- gamlss2parFam(mu = mu, sigma = sigma, - stabilization = stabilization, fname = fname) - }, { - ## 3 parameters - fun <- gamlss3parFam(mu = mu, sigma = sigma, nu = nu, - stabilization = stabilization, fname = fname) - }, { - ## 4 parameters - fun <- gamlss4parFam(mu = mu, sigma = sigma, nu = nu, tau = tau, - stabilization = stabilization, fname = fname) - }) - fun -} - - -################################################################################ -## 1 parameter - -gamlss1parMu <- function(mu = NULL, fname = "EXP") { - - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - dfun <- paste("gamlss.dist::d", fname, sep = "") - pdf <- eval(parse(text = dfun)) - - ## get the loss - loss <- function(y, f, w = 1) { - -pdf(x = y, mu = FAM$mu.linkinv(f), log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f)) - } - ## get the ngradient: mu is linkinv(f) - ## we need dl/deta = dl/dmu*dmu/deta - ngradient <- function(y, f, w = 1) { - FAM$dldm(y = y, mu = FAM$mu.linkinv(f)) * FAM$mu.dr(eta = f) - } - - ## get the offset -> we take the starting values of gamlss - offset <- function(y, w = 1) { - if (!is.null(mu)) { - RET <- FAM$mu.linkfun(mu) - } else { - eval(FAM$mu.initial) - RET <- FAM$mu.linkfun(mean(mu)) - } - return(RET) - } - - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$mu.linkinv(f), offset = offset, - name = paste(FAM$family[2], "(mboost family)")) -} - - -################################################################################ -## 2 parameters - -gamlss2parMu <- function(mu = NULL, sigma = NULL, - stabilization, fname = "NO") { - - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, sigma, w = 1) { - -pdf(x = y, mu = FAM$mu.linkinv(f), sigma = sigma, log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, sigma = sigma)) - } - - ## get the ngradient: mu is linkinv(f) - ## we need dl/deta = dl/dmu*dmu/deta - ngradient <- function(y, f, w = 1) { - ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma) * FAM$mu.dr(eta = f) - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - - ## get the offset -> we take the starting values of gamlss - offset <- function(y, w = 1) { - if (!is.null(mu)) { - RET <- FAM$mu.linkfun(mu) - } else { - eval(FAM$mu.initial) - RET <- FAM$mu.linkfun(mean(mu)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$mu.linkinv(f), offset = offset, - name = paste(FAM$family[2], "1st parameter (mu)")) -} - - -gamlss2parSigma <- function(mu = NULL, sigma = NULL, - stabilization, fname = "NO") { - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, w = 1, mu) { - -pdf(x = y, mu = mu, sigma = FAM$sigma.linkinv(f), log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, mu = mu)) - } - ## get the ngradient: sigma is linkinv(f) - ## we need dl/deta = dl/dsigma*dsigma/deta - ngradient <- function(y, f, w = 1) { - ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f)) * FAM$sigma.dr(eta = f) - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - ## get the offset - offset <- function(y, w = 1) { - if (!is.null(sigma)) { - RET <- FAM$sigma.linkfun(sigma) - } else { - eval(FAM$sigma.initial) - RET <- FAM$sigma.linkfun(mean(sigma)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$sigma.linkinv(f), offset = offset, - name = paste(FAM$family[2], "2nd parameter (sigma)")) -} - -## Build the Families object -gamlss2parFam <- function(mu = NULL, sigma = NULL, stabilization, fname = "NO") { - Families(mu = gamlss2parMu(mu = mu, sigma = sigma, stabilization, fname = fname), - sigma = gamlss2parSigma(mu = mu, sigma = sigma, stabilization, fname = fname), - qfun = get_qfun(fname), - name = fname) -} - -################################################################################ -## 3 parameters - -## sub-family for Mu -gamlss3parMu <- function(mu = NULL, sigma = NULL, nu = NULL, - stabilization, fname = "TF") { - - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, sigma, nu, w = 1) { - -pdf(x = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu, log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, sigma = sigma, nu = nu)) - } - ## get the ngradient: mu is linkinv(f) - ## we need dl/deta = dl/dmu*dmu/deta - ngradient <- function(y, f, w = 1) { - if (FAM$type == "Mixed") { - ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma) * FAM$mu.dr(eta = f) - } else { - ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu) * FAM$mu.dr(eta = f) - } - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - - ## get the offset -> we take the starting values of gamlss - offset <- function(y, w = 1) { - if (!is.null(mu)) { - RET <- FAM$mu.linkfun(mu) - } else { - eval(FAM$mu.initial) - RET <- FAM$mu.linkfun(mean(mu)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$mu.linkinv(f), offset = offset, - name = paste(FAM$family[2], "1st parameter (mu)")) -} - - -gamlss3parSigma <- function(mu = NULL, sigma = NULL, nu = NULL, - stabilization, fname = "TF") { - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, w = 1, mu, nu) { - -pdf(x = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu, log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, mu = mu, nu = nu)) - } - ## get the ngradient: sigma is linkinv(f) - ## we need dl/deta = dl/dsigma*dsigma/deta - ngradient <- function(y, f, w = 1) { - if (FAM$type == "Mixed") { - ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f)) * FAM$sigma.dr(eta = f) - } else { - ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu) * FAM$sigma.dr(eta = f) - } - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - ## get the offset - offset <- function(y, w = 1) { - if (!is.null(sigma)) { - RET <- FAM$sigma.linkfun(sigma) - } else { - eval(FAM$sigma.initial) - RET <- FAM$sigma.linkfun(mean(sigma)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$sigma.linkinv(f), offset = offset, - name = paste(FAM$family[2], "2nd parameter (sigma)")) -} - -gamlss3parNu <- function(mu = NULL, sigma = NULL, nu = NULL, - stabilization, fname = "TF") { - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, w = 1, mu, sigma) { - -pdf(x = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f), log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, mu = mu, sigma = sigma)) - } - ## get the ngradient: sigma is linkinv(f) - ## we need dl/deta = dl/dsigma*dsigma/deta - ngradient <- function(y, f, w = 1) { - if (FAM$type == "Mixed") { - ngr <- FAM$dldv(y = y, nu = FAM$nu.linkinv(f)) * FAM$nu.dr(eta = f) - } else { - ngr <- FAM$dldv(y = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f)) * FAM$nu.dr(eta = f) - } - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - ## get the offset - offset <- function(y, w = 1) { - if (!is.null(nu)) { - RET <- FAM$nu.linkfun(nu) - } else { - eval(FAM$nu.initial) - RET <- FAM$nu.linkfun(mean(nu)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$nu.linkinv(f), - offset = offset, name = paste(FAM$family[2], "3rd parameter (nu)")) -} - -## Build the Families object -gamlss3parFam <- function(mu = NULL, sigma = NULL, nu = NULL, stabilization, fname = "TF") { - Families(mu = gamlss3parMu(mu = mu, sigma = sigma, nu = nu, stabilization, fname = fname), - sigma = gamlss3parSigma(mu = mu, sigma = sigma, nu = nu, stabilization, fname = fname), - nu = gamlss3parNu(mu = mu, sigma = sigma, nu = nu, stabilization, fname = fname), - qfun = get_qfun(fname), - name = fname) -} - - -################################################################################ -## 4 parameters - -gamlss4parMu <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, - stabilization, fname = "BCPE") { - - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, sigma, nu, tau, w = 1) { - -pdf(x = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu, tau = tau, log = TRUE) - } - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, sigma = sigma, nu = nu, tau = tau)) - } - ## get the ngradient: mu is linkinv(f) - ## we need dl/deta = dl/dmu*dmu/deta - ngradient <- function(y, f, w = 1) { - ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu, tau = tau) * - FAM$mu.dr(eta = f) - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - ## get the offset -> we take the starting values of gamlss - offset <- function(y, w = 1) { - if (!is.null(mu)) { - RET <- FAM$mu.linkfun(mu) - } else { - eval(FAM$mu.initial) - RET <- FAM$mu.linkfun(mean(mu)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$mu.linkinv(f), offset = offset, - name = paste(FAM$family[2], "1st parameter (mu)")) -} - - -gamlss4parSigma <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, - stabilization, fname = "BCPE") { - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, w = 1, mu, nu, tau) { - -pdf(x = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu, tau = tau, - log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, mu = mu, nu = nu, tau = tau)) - } - ## get the ngradient: sigma is linkinv(f) - ## we need dl/deta = dl/dsigma*dsigma/deta - ngradient <- function(y, f, w = 1) { - ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu, - tau = tau) * FAM$sigma.dr(eta = f) - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - ## get the offset - offset <- function(y, w = 1) { - if (!is.null(sigma)) { - RET <- FAM$sigma.linkfun(sigma) - } else { - eval(FAM$sigma.initial) - RET <- FAM$sigma.linkfun(mean(sigma)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$sigma.linkinv(f), offset = offset, - name = paste(FAM$family[2], "2nd parameter (sigma)")) -} - -gamlss4parNu <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, - stabilization, fname = "BCPE") { - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, w = 1, mu, sigma, tau) { - -pdf(x = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f), tau = tau, - log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, mu = mu, sigma = sigma, tau = tau)) - } - ## get the ngradient: sigma is linkinv(f) - ## we need dl/deta = dl/dsigma*dsigma/deta - ngradient <- function(y, f, w = 1) { - ngr <- FAM$dldv(y = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f), - tau = tau) * FAM$nu.dr(eta = f) - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - ## get the offset - offset <- function(y, w = 1) { - if (!is.null(nu)) { - RET <- FAM$nu.linkfun(nu) - } else { - eval(FAM$nu.initial) - RET <- FAM$nu.linkfun(mean(nu)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$nu.linkinv(f), offset = offset, - name = paste(FAM$family[2], "3rd parameter (nu)")) -} - - - -gamlss4parTau <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, - stabilization, fname = "BCPE") { - FAM <- gamlss.dist::as.gamlss.family(fname) - NAMEofFAMILY <- FAM$family - pdf <- get_pdf(fname) - - ## get the loss - loss <- function(y, f, w = 1, mu, sigma, nu) { - -pdf(x = y, mu = mu, sigma = sigma, nu = nu, tau = FAM$tau.linkinv(f), - log = TRUE) - } - - ## compute the risk - risk <- function(y, f, w = 1) { - sum(w * loss(y = y, f = f, mu = mu, sigma = sigma, nu = nu)) - } - ## get the ngradient - ngradient <- function(y, f, w = 1) { - ngr <- FAM$dldt(y = y, mu = mu, sigma = sigma, tau = FAM$tau.linkinv(f), - nu = nu) * FAM$tau.dr(eta = f) - ngr <- stabilize_ngradient(ngr, w = w, stabilization) - ngr - } - ## get the offset - offset <- function(y, w = 1) { - if (!is.null(tau)) { - RET <- FAM$tau.linkfun(tau) - } else { - eval(FAM$tau.initial) - RET <- FAM$tau.linkfun(mean(tau)) - } - return(RET) - } - Family(ngradient = ngradient, risk = risk, loss = loss, - response = function(f) FAM$tau.linkinv(f), offset = offset, - name = paste(FAM$family[2], "4th parameter (tau)")) -} - -## Build the Families object -gamlss4parFam <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, stabilization, fname = "BCPE") { - Families(mu = gamlss4parMu(mu = mu, sigma = sigma, nu = nu, tau = tau, stabilization, fname = fname), - sigma = gamlss4parSigma(mu = mu, sigma = sigma, nu = nu, tau = tau, stabilization, fname = fname), - nu = gamlss4parNu(mu = mu, sigma = sigma, nu = nu, tau = tau, stabilization, fname = fname), - tau = gamlss4parTau(mu = mu, sigma = sigma, nu = nu, tau = tau, stabilization, fname = fname), - qfun = get_qfun(fname), - name = fname) -} +################################################################################ +### Family wrapper for gamlss families + + + + +################################################################################ +## constructor + +## a wrapper to as.families (for compatibility reasons +gamlss.Families <- function(...) + as.families(...) + +as.families <- function(fname = "NO", stabilization = c("none", "MAD", "L2"), + mu = NULL, sigma = NULL, nu = NULL, tau = NULL, + mu.link = NULL, sigma.link = NULL, nu.link = NULL, + tau.link = NULL) { + + ## require gamlss.dist + if (!requireNamespace("gamlss.dist", quietly = TRUE)) + stop("Please install package 'gamlss.dist' for using gamlss families.") + + if (is.function(fname)) + fname <- fname() + + if (inherits(fname, "gamlss.family")) + fname <- fname$family[1] + + if (mode(fname) != "character" && mode(fname) != "name") + fname <- as.character(substitute(fname)) + + gamlss.fam <- try(gamlss.dist::gamlss.family(fname), silent = TRUE) + if (inherits(gamlss.fam, "try-error")) + stop(sQuote("fname"), " specifies no valid gamlss family") + + stabilization <- check_stabilization(stabilization) + + npar <- gamlss.fam$nopar + switch(npar, { + ## 1 parameter + fun <- gamlss1parMu(mu = mu, fname = fname, mu.link = mu.link) + warning("For boosting one-parametric families,", + " please use the mboost package.") + if (stabilization != "none") + warning("Stabilization is ignored for one-parametric families.") + }, { + ## 2 parameters + fun <- gamlss2parFam(mu = mu, sigma = sigma, mu.link = mu.link, + sigma.link = sigma.link, + stabilization = stabilization, fname = fname) + }, { + ## 3 parameters + fun <- gamlss3parFam(mu = mu, sigma = sigma, nu = nu, + mu.link = mu.link, sigma.link = sigma.link, + nu.link = nu.link, + stabilization = stabilization, fname = fname) + }, { + ## 4 parameters + fun <- gamlss4parFam(mu = mu, sigma = sigma, nu = nu, tau = tau, + mu.link = mu.link, sigma.link = sigma.link, + nu.link = nu.link, tau.link = tau.link, + stabilization = stabilization, fname = fname) + }) + fun +} + + +################################################################################ +## 1 parameter + +gamlss1parMu <- function(mu = NULL, fname = "EXP", mu.link = mu.link) { + + ## check if a specific link was provided + if (!is.null(mu.link)) fname_link <- eval(parse(text = paste(fname, + "(mu.link = '", mu.link, "')", sep =""))) + else fname_link <- fname + + FAM <- gamlss.dist::as.gamlss.family(fname_link) + NAMEofFAMILY <- FAM$family + dfun <- paste("gamlss.dist::d", fname, sep = "") + pdf <- eval(parse(text = dfun)) + is.bdfamily <- "bd" %in% names(formals(pdf)) + + + ## get the loss + loss <- function(y, f, w = 1) { + if (is.bdfamily) { + bd <- rowSums(y) + y <- y[,1] + return( -pdf(x = y, mu = FAM$mu.linkinv(f), log = TRUE, bd = bd)) + } + -pdf(x = y, mu = FAM$mu.linkinv(f), log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f)) + } + ## get the ngradient: mu is linkinv(f) + ## we need dl/deta = dl/dmu*dmu/deta + ngradient <- function(y, f, w = 1) { + if (is.bdfamily) { + if (!is.matrix(y)) stop("y should be a matrix for this family") + bd <- rowSums(y) + y <- y[,1] + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), bd = bd) * FAM$mu.dr(eta = f) + } else { + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f)) * FAM$mu.dr(eta = f) + } + return(ngr) + } + + ## get the offset -> we take the starting values of gamlss + offset <- function(y, w = 1) { + if (!is.null(mu)) { + RET <- FAM$mu.linkfun(mu) + } else { + if (is.bdfamily) { + if (!is.matrix(y)) stop("y should be a matrix for this family") + bd <- rowSums(y) + y <- y[,1] + } + eval(FAM$mu.initial) + RET <- FAM$mu.linkfun(mean(mu)) + } + return(RET) + } + + mboost::Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$mu.linkinv(f), offset = offset, + name = paste(FAM$family[2], "(mboost family)")) +} + + +################################################################################ +## 2 parameters + +gamlss2parMu <- function(mu = NULL, sigma = NULL, mu.link = NULL, FAM = FAM, + stabilization = stabilization, fname = "NO") { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + is.bdfamily <- "bd" %in% names(formals(pdf)) + + ## get the loss + loss <- function(y, f, sigma, w = 1) { + if (is.bdfamily) { + #if (!is.matrix(y)) stop("y should be a matrix for this family") + bd <- rowSums(y) + y <- y[,1] + return( -pdf(x = y, mu = FAM$mu.linkinv(f), sigma = sigma, log = TRUE, bd = bd)) + } + + -pdf(x = y, mu = FAM$mu.linkinv(f), sigma = sigma, log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, sigma = sigma)) + } + + ## get the ngradient: mu is linkinv(f) + ## we need dl/deta = dl/dmu*dmu/deta + ngradient <- function(y, f, w = 1) { + if (is.bdfamily) { + if (!is.matrix(y)) stop("y should be a matrix for this family") + bd <- rowSums(y) + y <- y[,1] + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma, bd = bd) * FAM$mu.dr(eta = f) + } else { + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma) * FAM$mu.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + + ## get the offset -> we take the starting values of gamlss + offset <- function(y, w = 1) { + if (!is.null(mu)) { + RET <- FAM$mu.linkfun(mu) + } else { + if (is.bdfamily) { + if (!is.matrix(y)) stop("y should be a matrix for this family") + bd <- rowSums(y) + y <- y[,1] + } + eval(FAM$mu.initial) + RET <- FAM$mu.linkfun(mean(mu)) + } + return(RET) + } + + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$mu.linkinv(f), offset = offset, + name = paste(FAM$family[2], "1st parameter (mu)")) +} + + +gamlss2parSigma <- function(mu = NULL, sigma = NULL, + stabilization = stabilization, fname = "NO", FAM = FAM) { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + is.bdfamily <- "bd" %in% names(formals(pdf)) + + ## get the loss + loss <- function(y, f, w = 1, mu ) { + # check if bd exists in this family + if (is.bdfamily) { + bd <- rowSums(y) + y <- y[,1] + return( -pdf(x = y, mu = mu, sigma = FAM$sigma.linkinv(f), log = TRUE, bd = bd)) + } + + -pdf(x = y, mu = mu, sigma = FAM$sigma.linkinv(f), log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, mu = mu)) + } + ## get the ngradient: sigma is linkinv(f) + ## we need dl/deta = dl/dsigma*dsigma/deta + ngradient <- function(y, f, w = 1) { + + # check if bd exists in this family + if (is.bdfamily) { + if (!is.matrix(y)) stop("y should be a matrix for this family") + bd <- rowSums(y) + y <- y[,1] + ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f), bd = bd) * FAM$sigma.dr(eta = f) + + } else { + ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f)) * FAM$sigma.dr(eta = f) + } + + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + ## get the offset + offset <- function(y, w = 1) { + if (!is.null(sigma)) { + RET <- FAM$sigma.linkfun(sigma) + } else { + eval(FAM$sigma.initial) + RET <- FAM$sigma.linkfun(mean(sigma)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$sigma.linkinv(f), offset = offset, + name = paste(FAM$family[2], "2nd parameter (sigma)")) +} + +## Build the Families object +gamlss2parFam <- function(mu = NULL, sigma = NULL, mu.link = mu.link, + sigma.link = sigma.link, stabilization, fname = "NO") { + + FAM <- gamlss.dist::as.gamlss.family(fname) + # check if any link function was set + if (any(!is.null(mu.link), !is.null(sigma.link))) { + # if some are null, set default + if (is.null(mu.link)) mu.link <- FAM$mu.link + if (is.null(sigma.link)) sigma.link <- FAM$sigma.link + fname_link <- paste("gamlss.dist::",fname, "(mu.link = '", mu.link, "', ", + "sigma.link = '", sigma.link, "'", ")", sep ="") + # build family with link + FAM <- gamlss.dist::as.gamlss.family(eval(parse(text = fname_link))) + } + Families(mu = gamlss2parMu(mu = mu, sigma = sigma, FAM = FAM, + stabilization = stabilization, fname = fname), + sigma = gamlss2parSigma(mu = mu, sigma = sigma, FAM = FAM, + stabilization = stabilization, fname = fname), + qfun = get_qfun(fname), + name = fname) +} + +################################################################################ +## 3 parameters + +## sub-family for Mu +gamlss3parMu <- function(mu = NULL, sigma = NULL, nu = NULL, FAM = FAM, + stabilization = stabilization, fname = "TF") { + + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + + ## get the loss + loss <- function(y, f, sigma, nu, w = 1) { + -pdf(x = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu, log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, sigma = sigma, nu = nu)) + } + ## get the ngradient: mu is linkinv(f) + ## we need dl/deta = dl/dmu*dmu/deta + ngradient <- function(y, f, w = 1) { + if (FAM$type == "Mixed") { + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma) * FAM$mu.dr(eta = f) + } else { + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu) * FAM$mu.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + + ## get the offset -> we take the starting values of gamlss + offset <- function(y, w = 1) { + if (!is.null(mu)) { + RET <- FAM$mu.linkfun(mu) + } else { + eval(FAM$mu.initial) + RET <- FAM$mu.linkfun(mean(mu)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$mu.linkinv(f), offset = offset, + name = paste(FAM$family[2], "1st parameter (mu)")) +} + + +gamlss3parSigma <- function(mu = NULL, sigma = NULL, nu = NULL, FAM = FAM, + stabilization = stabilization, fname = "TF") { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + + ## get the loss + loss <- function(y, f, w = 1, mu, nu) { + -pdf(x = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu, log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, mu = mu, nu = nu)) + } + ## get the ngradient: sigma is linkinv(f) + ## we need dl/deta = dl/dsigma*dsigma/deta + ngradient <- function(y, f, w = 1) { + if (FAM$type == "Mixed") { + ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f)) * FAM$sigma.dr(eta = f) + } else { + ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu) * FAM$sigma.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + ## get the offset + offset <- function(y, w = 1) { + if (!is.null(sigma)) { + RET <- FAM$sigma.linkfun(sigma) + } else { + eval(FAM$sigma.initial) + RET <- FAM$sigma.linkfun(mean(sigma)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$sigma.linkinv(f), offset = offset, + name = paste(FAM$family[2], "2nd parameter (sigma)")) +} + +gamlss3parNu <- function(mu = NULL, sigma = NULL, nu = NULL, + stabilization = stabilization, fname = "TF", FAM = FAM) { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + + ## get the loss + loss <- function(y, f, w = 1, mu, sigma) { + -pdf(x = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f), log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, mu = mu, sigma = sigma)) + } + ## get the ngradient: sigma is linkinv(f) + ## we need dl/deta = dl/dsigma*dsigma/deta + ngradient <- function(y, f, w = 1) { + if (FAM$type == "Mixed") { + ngr <- FAM$dldv(y = y, nu = FAM$nu.linkinv(f)) * FAM$nu.dr(eta = f) + } else { + ngr <- FAM$dldv(y = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f)) * FAM$nu.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + ## get the offset + offset <- function(y, w = 1) { + if (!is.null(nu)) { + RET <- FAM$nu.linkfun(nu) + } else { + eval(FAM$nu.initial) + RET <- FAM$nu.linkfun(mean(nu)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$nu.linkinv(f), + offset = offset, name = paste(FAM$family[2], "3rd parameter (nu)")) +} + +## Build the Families object +gamlss3parFam <- function(mu = NULL, sigma = NULL, nu = NULL, + mu.link = NULL, sigma.link = NULL, nu.link = NULL, + stabilization = stabilization, fname = "TF") { + + FAM <- gamlss.dist::as.gamlss.family(fname) + # check if any link function was set + if(any(!is.null(mu.link), !is.null(sigma.link), !is.null(nu.link))){ + # if some are null, set default + if (is.null(mu.link)) mu.link <- FAM$mu.link + if (is.null(sigma.link)) sigma.link <- FAM$sigma.link + if (is.null(nu.link)) nu.link <- FAM$nu.link + + fname_link <- paste("gamlss.dist::",fname, "(mu.link = '", mu.link, "', ", + "sigma.link = '", sigma.link, "',", + "nu.link = '", nu.link, "')", sep ="") + # build family with link + FAM <- gamlss.dist::as.gamlss.family(eval(parse(text = fname_link))) + } + + Families(mu = gamlss3parMu(mu = mu, sigma = sigma, nu = nu, + stabilization = stabilization, fname = fname, FAM = FAM), + sigma = gamlss3parSigma(mu = mu, sigma = sigma, nu = nu, + stabilization = stabilization, fname = fname, FAM = FAM), + nu = gamlss3parNu(mu = mu, sigma = sigma, nu = nu, + stabilization = stabilization, fname = fname, FAM = FAM), + qfun = get_qfun(fname), + name = fname) +} + + +################################################################################ +## 4 parameters + +gamlss4parMu <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, + stabilization = stabilization, fname = "BCPE", FAM = FAM) { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + + ## get the loss + loss <- function(y, f, sigma, nu, tau, w = 1) { + -pdf(x = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu, tau = tau, log = TRUE) + } + ## compute the risk + risk <- function(y, f, w = 1) { + ## get the ngradient: mu is linkinv(f) + sum(w * loss(y = y, f = f, sigma = sigma, nu = nu, tau = tau)) + } + ## we need dl/deta = dl/dmu*dmu/deta + ngradient <- function(y, f, w = 1) { + if (FAM$type == "Mixed") { + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma) * FAM$mu.dr(eta = f) + } else { + ngr <- FAM$dldm(y = y, mu = FAM$mu.linkinv(f), sigma = sigma, nu = nu, tau = tau) * + FAM$mu.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + ## get the offset -> we take the starting values of gamlss + offset <- function(y, w = 1) { + if (!is.null(mu)) { + RET <- FAM$mu.linkfun(mu) + } else { + eval(FAM$mu.initial) + RET <- FAM$mu.linkfun(mean(mu)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$mu.linkinv(f), offset = offset, + name = paste(FAM$family[2], "1st parameter (mu)")) +} + + +gamlss4parSigma <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, FAM = FAM, + stabilization = stabilization, fname = "BCPE") { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + + ## get the loss + loss <- function(y, f, w = 1, mu, nu, tau) { + -pdf(x = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu, tau = tau, + log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, mu = mu, nu = nu, tau = tau)) + } + ## get the ngradient: sigma is linkinv(f) + ## we need dl/deta = dl/dsigma*dsigma/deta + ngradient <- function(y, f, w = 1) { + if (FAM$type == "Mixed") { + ngr <- FAM$dldm(y = y, mu = mu, sigma = FAM$sigma.linkinv(f)) * FAM$sigma.dr(eta = f) + } else { + ngr <- FAM$dldd(y = y, mu = mu, sigma = FAM$sigma.linkinv(f), nu = nu, + tau = tau) * FAM$sigma.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + ## get the offset + offset <- function(y, w = 1) { + if (!is.null(sigma)) { + RET <- FAM$sigma.linkfun(sigma) + } else { + eval(FAM$sigma.initial) + RET <- FAM$sigma.linkfun(mean(sigma)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$sigma.linkinv(f), offset = offset, + name = paste(FAM$family[2], "2nd parameter (sigma)")) +} + +gamlss4parNu <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, FAM = FAM, + stabilization = stabilization, fname = "BCPE") { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + + ## get the loss + loss <- function(y, f, w = 1, mu, sigma, tau) { + -pdf(x = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f), tau = tau, + log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, mu = mu, sigma = sigma, tau = tau)) + } + ## get the ngradient: sigma is linkinv(f) + ## we need dl/deta = dl/dsigma*dsigma/deta + ngradient <- function(y, f, w = 1) { + if (FAM$type == "Mixed") { + ngr <- FAM$dldv(y = y, nu = FAM$nu.linkinv(f), tau = tau) * FAM$nu.dr(eta = f) + } else { + ngr <- FAM$dldv(y = y, mu = mu, sigma = sigma, nu = FAM$nu.linkinv(f), + tau = tau) * FAM$nu.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + ## get the offset + offset <- function(y, w = 1) { + if (!is.null(nu)) { + RET <- FAM$nu.linkfun(nu) + } else { + eval(FAM$nu.initial) + RET <- FAM$nu.linkfun(mean(nu)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$nu.linkinv(f), offset = offset, + name = paste(FAM$family[2], "3rd parameter (nu)")) +} + + + +gamlss4parTau <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, FAM = FAM, + stabilization = stabilization, fname = "BCPE") { + + NAMEofFAMILY <- FAM$family + pdf <- get_pdf(fname) + + ## get the loss + loss <- function(y, f, w = 1, mu, sigma, nu) { + -pdf(x = y, mu = mu, sigma = sigma, nu = nu, tau = FAM$tau.linkinv(f), + log = TRUE) + } + + ## compute the risk + risk <- function(y, f, w = 1) { + sum(w * loss(y = y, f = f, mu = mu, sigma = sigma, nu = nu)) + } + ## get the ngradient + ngradient <- function(y, f, w = 1) { + if (FAM$type == "Mixed") { + ngr <- FAM$dldt(y = y, nu = nu, tau = FAM$tau.linkinv(f)) * FAM$tau.dr(eta = f) + } else { + ngr <- FAM$dldt(y = y, mu = mu, sigma = sigma, tau = FAM$tau.linkinv(f), + nu = nu) * FAM$tau.dr(eta = f) + } + ngr <- stabilize_ngradient(ngr, w = w, stabilization) + ngr + } + ## get the offset + offset <- function(y, w = 1) { + if (!is.null(tau)) { + RET <- FAM$tau.linkfun(tau) + } else { + eval(FAM$tau.initial) + RET <- FAM$tau.linkfun(mean(tau)) + } + return(RET) + } + Family(ngradient = ngradient, risk = risk, loss = loss, + response = function(f) FAM$tau.linkinv(f), offset = offset, + name = paste(FAM$family[2], "4th parameter (tau)")) +} + +## Build the Families object +gamlss4parFam <- function(mu = NULL, sigma = NULL, nu = NULL, tau = NULL, + mu.link = NULL, sigma.link = NULL, nu.link = NULL, tau.link = NULL, + stabilization = stabilization, fname = "BCPE") { + + FAM <- gamlss.dist::as.gamlss.family(fname) + # check if any link function was set + if (any(!is.null(mu.link), !is.null(sigma.link), !is.null(nu.link), + !is.null(tau.link))) { + # if some are null, set default + if (is.null(mu.link)) mu.link <- FAM$mu.link + if (is.null(sigma.link)) sigma.link <- FAM$sigma.link + if (is.null(nu.link)) nu.link <- FAM$nu.link + if (is.null(tau.link)) tau.link <- FAM$tau.link + + + fname_link <- paste("gamlss.dist::",fname, "(mu.link = '", mu.link, "', ", + "sigma.link = '", sigma.link, "',", + "nu.link = '", nu.link, "',", + "tau.link = '", tau.link, "')", sep ="") + # build family with link + FAM <- gamlss.dist::as.gamlss.family(eval(parse(text = fname_link))) + } + + + Families(mu = gamlss4parMu(mu = mu, sigma = sigma, nu = nu, tau = tau, + stabilization = stabilization, fname = fname, FAM = FAM), + sigma = gamlss4parSigma(mu = mu, sigma = sigma, nu = nu, tau = tau, + stabilization = stabilization, fname = fname, FAM = FAM), + nu = gamlss4parNu(mu = mu, sigma = sigma, nu = nu, tau = tau, + stabilization = stabilization, fname = fname, FAM = FAM), + tau = gamlss4parTau(mu = mu, sigma = sigma, nu = nu, tau = tau, + stabilization = stabilization, fname = fname, FAM = FAM), + qfun = get_qfun(fname), + name = fname) +} diff --git a/R/cvrisk.R b/R/cvrisk.R index 5ecf8fd..f85c70d 100644 --- a/R/cvrisk.R +++ b/R/cvrisk.R @@ -4,7 +4,7 @@ make.grid <- function(max, length.out = 10, min = NULL, log = TRUE, dense_mu_grid = TRUE) { - + if (is.null(min)) min <- rep(1, length(max)) if (length(min) == 1) @@ -19,15 +19,15 @@ make.grid <- function(max, length.out = 10, min = NULL, log = TRUE, stop(sQuote("min"), " must be either scalar or a vector of the same length as ", sQuote("max")) - + if (log == TRUE) { min <- log(min) max <- log(max) } - + if (any(sapply(1:length(max), function(i) min[i] >= max[i]))) stop("All min values must be smaller than the respectiv max value.") - + ## single paramter family if (length(max) == 1) { RET <- seq(from = min, to = max, length.out = length.out) @@ -41,20 +41,20 @@ make.grid <- function(max, length.out = 10, min = NULL, log = TRUE, } return(RET) } - + ## ELSE: multiple parameter family - + if (is.null(names(max))) stop(sQuote("max"), " must be a named vector") - + RET <- lapply(1:length(max), function(i) - seq(from = min[i], to = max[i], - length.out = length.out[i])) + seq(from = min[i], to = max[i], + length.out = length.out[i])) if (log == TRUE) RET <- lapply(RET, exp) ## round to get integer values RET <- lapply(RET, round) - + if (any(sapply(RET, function(x) any(duplicated(x))))) { warning("Duplicates produced; Only unique values are returned") RET <- lapply(RET, unique) @@ -63,11 +63,11 @@ make.grid <- function(max, length.out = 10, min = NULL, log = TRUE, RET <- expand.grid(RET) if (!is.null(names(max))) colnames(RET) <- names(max) - + ## no way to produce dense grid if (dense_mu_grid && max(RET[,1]) <= min(RET[, -1])) dense_mu_grid <- FALSE - + ## if dense grid for mu is requested: if (dense_mu_grid) { ## find all grid points where at least one mu grid point is greater than @@ -99,7 +99,7 @@ cvrisk.mboostLSS <- function(object, folds = cv(model.weights(object)), grid = make.grid(mstop(object)), papply = mclapply, trace = TRUE, fun = NULL, ...) { - + weights <- model.weights(object) if (any(weights == 0)) warning("Zero weights in ", sQuote("object")) @@ -112,10 +112,10 @@ cvrisk.mboostLSS <- function(object, folds = cv(model.weights(object)), stop(sQuote("grid"), " must be a matrix with the same number of columns", " as parameters in", sQuote("object")) - + if (!is.null(fun)) stopifnot(is.function(fun)) - + ### WHAT ABOUT: ## fam_name <- object$family@name call <- deparse(attr(object, "call")) @@ -134,10 +134,10 @@ cvrisk.mboostLSS <- function(object, folds = cv(model.weights(object)), mod <- update(object, weights = weights, oobweights = oobweights, risk = "oobag", trace = FALSE, mstop = apply(grid, 2, min)) - + ## now we need to increase mstop (stupid or clever) risks <- vector("numeric", nrow(grid)) - + ## fitting loop j <- 1 while (j <= nrow(grid)) { @@ -178,7 +178,7 @@ cvrisk.mboostLSS <- function(object, folds = cv(model.weights(object)), mod <- update(object, weights = weights, oobweights = oobweights, risk = "oobag", trace = FALSE, mstop = apply(grid, 2, min)) - + res <- vector("list", nrow(grid)) for (i in 1:nrow(grid)) { mod[grid[i, ]] @@ -188,12 +188,12 @@ cvrisk.mboostLSS <- function(object, folds = cv(model.weights(object)), return(res) } } - + OOBweights <- matrix(rep(weights, ncol(folds)), ncol = ncol(folds)) OOBweights[folds > 0] <- 0 oobrisk <- papply(1:ncol(folds), - function(i) dummyfct(i, weights = folds[, i], - oobweights = OOBweights[, i]), ...) + function(i) dummyfct(i, weights = folds[, i], + oobweights = OOBweights[, i]), ...) ## get errors if mclapply is used if (any(idx <- sapply(oobrisk, is.character))) stop(sapply(oobrisk[idx], function(x) x)) @@ -209,7 +209,7 @@ cvrisk.mboostLSS <- function(object, folds = cv(model.weights(object)), attr(oobrisk, "call") <- call attr(oobrisk, "mstop") <- grid attr(oobrisk, "type") <- ifelse(!is.null(attr(folds, "type")), - attr(folds, "type"), "user-defined") + attr(folds, "type"), "user-defined") class(oobrisk) <- "cvriskLSS" oobrisk } @@ -217,7 +217,7 @@ cvrisk.mboostLSS <- function(object, folds = cv(model.weights(object)), print.cvriskLSS <- function(x, ...) { #cat("\n\t Cross-validated", attr(x, "risk"), "\n\t", cat("\n\t Cross-validated risk\n\t", - attr(x, "call"), "\n\n") + attr(x, "call"), "\n\n") print(colMeans(x)) cat("\n\t Optimal number of boosting iterations:", mstop(x), "\n") return(invisible(x)) @@ -227,10 +227,10 @@ plot.cvriskLSS <- function(x, type = c("heatmap", "lines"), xlab = NULL, ylab = NULL, ylim = range(x), main = attr(x, "type"), - ...) { - + ...) { + type <- match.arg(type) - + nms <- names(attr(x, "mstop")) if (type == "lines") { if (is.null(xlab)) @@ -244,7 +244,7 @@ plot.cvriskLSS <- function(x, type = c("heatmap", "lines"), if (is.null(ylab)) ylab <- paste0("Number of boosting iterations (", nms[2], ")") } - + if (type == "lines") { cm <- colMeans(x) plot(1:ncol(x), cm, xlab = xlab, ylab = ylab, @@ -304,7 +304,7 @@ mstop.cvriskLSS <- function(object, parameter = NULL, ...) { if (FALSE) { library(gamboostLSS) - + ## check cvrisk set.seed(1907) x1 <- rnorm(1000) @@ -327,9 +327,9 @@ if (FALSE) { abline(0,1) cvr <- cvrisk(model, folds = cv(model.weights(model), B = 5), grid = grid, papply = lapply) - + ### check timings: - + ### Data generating process: set.seed(1907) x1 <- rnorm(1000) @@ -344,7 +344,7 @@ if (FALSE) { 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) - + system.time({ ## linear model with y ~ . for both components: 1 boosting iterations model <- glmboostLSS(y ~ ., families = NBinomialLSS(), data = dat, @@ -361,9 +361,9 @@ if (FALSE) { control = boost_control(mstop = 1), center = TRUE) model[c(1000, 10)] - + }) - - + + ## what about the warnings in 3d? } diff --git a/R/cvrisk.nc_mboostLSS.R b/R/cvrisk.nc_mboostLSS.R new file mode 100644 index 0000000..0743e34 --- /dev/null +++ b/R/cvrisk.nc_mboostLSS.R @@ -0,0 +1,106 @@ +cvrisk.nc_mboostLSS <- function(object, folds = cv(model.weights(object)), + grid = 1:sum(mstop(object)), + papply = mclapply, trace = TRUE, + mc.preschedule = FALSE, fun = NULL, ...) { + + weights <- model.weights(object) + if (any(weights == 0)) + warning("Zero weights in ", sQuote("object")) + if (is.null(folds)) { + folds <- rmultinom(25, length(weights), weights/sum(weights)) + } else { + stopifnot(is.matrix(folds) && nrow(folds) == length(weights)) + } + + if (!is.null(fun)) + stopifnot(is.function(fun)) + + ### WHAT ABOUT: + ## fam_name <- object$family@name + call <- deparse(attr(object, "call")) + oobrisk <- matrix(0, nrow = ncol(folds), ncol = length(grid)) + if (trace) + cat("Starting cross-validation...") + if (is.null(fun)) { + dummyfct <- function(i, weights, oobweights) { + if (trace) + cat("\n[Fold: ", i, "]\n", sep = "") + ## make model with new weights and max mstop + mod <- update(object, weights = weights, oobweights = oobweights, + risk = "oobag", + mstop = max(grid), trace = trace) + + + risks <- attr(mod, "combined_risk")()[grid] + names(risks) <- grid + risks + + } + } + else { + stop("currently not implemented") + } + + OOBweights <- matrix(rep(weights, ncol(folds)), ncol = ncol(folds)) + OOBweights[folds > 0] <- 0 + if (all.equal(papply, mclapply) == TRUE) { + oobrisk <- papply(1:ncol(folds), + function(i) dummyfct(i = i, + weights = folds[, i], + oobweights = OOBweights[, i]), + mc.preschedule = mc.preschedule, + ...) + } else { + oobrisk <- papply(1:ncol(folds), + function(i) try(dummyfct(i = i, + weights = folds[, i], + oobweights = OOBweights[, i])), + ...) + } + ## if any errors occured remove results and issue a warning + if (any(idx <- sapply(oobrisk, is.character))) { + warning(sum(idx), " fold(s) encountered an error. ", + "Results are based on ", ncol(folds) - sum(idx), + " folds only.\n", + "Original error message(s):\n", + sapply(oobrisk[idx], function(x) x)) + oobrisk[idx] <- NULL + } + if (!is.null(fun)) + return(oobrisk) + oobrisk <- t(as.data.frame(oobrisk)) + oobrisk <- oobrisk / colSums(OOBweights) + colnames(oobrisk) <- grid + rownames(oobrisk) <- 1:nrow(oobrisk) + #attr(oobrisk, "risk") <- fam_name + attr(oobrisk, "call") <- call + attr(oobrisk, "mstop") <- grid + attr(oobrisk, "type") <- ifelse(!is.null(attr(folds, "type")), + attr(folds, "type"), "user-defined") + class(oobrisk) <- c("nc_cvriskLSS", "cvrisk") + oobrisk +} + + +risk.nc_mboostLSS <- function(object, merge = FALSE, + parameter = names(object), ...){ + + if(merge) attr(object, "combined_risk")() + else{ + risk.mboostLSS(object, merge = FALSE, parameter = names(object), ...) + } + +} + + +plot.nc_cvriskLSS <- function(x, type = "lines", + xlab = "Number of boosting iterations", ylab = NULL, + ylim = range(x), + main = attr(x, "type"), + ...) { + if (type != "lines") + warning("Only ", sQuote('type = "lines"'), " supported for noncyclical fitting") + plot.cvriskLSS(x = x, type = "lines", xlab = xlab, ylab = ylab, + ylim = ylim, main = main, ...) +} + diff --git a/R/families.R b/R/families.R index 7a7cdf3..8c5c9a0 100644 --- a/R/families.R +++ b/R/families.R @@ -8,7 +8,7 @@ Families <- function(..., qfun = NULL, name = NULL) { check <- sapply(RET, function(x){ bdy <- body(x@response) (length(x@response(c(-1,0,1))) != 3 && class(bdy) != "call" && - length(bdy) == 1 && is.na(bdy)) + length(bdy) == 1 && is.na(bdy)) }) if (any(check)) stop("response function not specified in families for:\n\t", @@ -48,7 +48,7 @@ NBinomialMu <- function(mu = NULL, sigma = NULL, stabilization) { } return(RET) } - + Family(ngradient = ngradient, offset = offset, risk = risk, loss = loss, response = function(f) exp(f), check_y = function(y) { @@ -64,14 +64,14 @@ NBinomialSigma <- function(mu = NULL, sigma = NULL, stabilization) { # this family boosts sigma therefore f is sigma loss <- function(mu, y, f) -(lgamma(y + exp(f)) - lgamma(exp(f)) - lgamma(y + 1) + exp(f) * f + - y * log(mu) - (y + exp(f)) * log(mu + exp(f))) + y * log(mu) - (y + exp(f)) * log(mu + exp(f))) risk <- function(y, f, w = 1){ RET <- sum(w * loss(y = y, f = f, mu = mu)) return(RET) } ngradient <- function(y, f, w = 1) { # f is sigma ! ngr <- exp(f)*(digamma(y +exp(f)) - digamma(exp(f)) + log(exp(f)) + 1 - - log(mu +exp(f)) - (exp(f) + y)/(mu +exp(f))) + log(mu +exp(f)) - (exp(f) + y)/(mu +exp(f))) ngr <- stabilize_ngradient(ngr, w = w, stabilization) return(ngr) } @@ -87,7 +87,7 @@ NBinomialSigma <- function(mu = NULL, sigma = NULL, stabilization) { } return(RET) } - + Family(ngradient = ngradient, offset = offset, risk = risk, loss = loss, response = function(f) exp(f), check_y = function(y) { @@ -101,7 +101,7 @@ NBinomialSigma <- function(mu = NULL, sigma = NULL, stabilization) { NBinomialLSS <- function(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { if ((!is.null(sigma) && sigma <= 0) || (!is.null(mu) && mu <= 0)) stop(sQuote("sigma"), " and ", sQuote("mu"), " must be greater than zero") @@ -127,7 +127,7 @@ qNBinomial <- function(p, mu = 0, sigma = 1, lower.tail = TRUE, log.p = FALSE) { StudentTMu <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { 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))) + log(df) - (df+1)/2 * log(1 + (y-f)^2 / (df * sigma^2))) } risk <- function(y, f, w = 1){ sum(w * loss(y = y, f = f, df = df, sigma = sigma)) @@ -136,8 +136,8 @@ StudentTMu <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { ngr <- (df+1)*(y-f)/(df*sigma^2 +(y-f)^2) ngr <- stabilize_ngradient(ngr, w = w, stabilization) return(ngr) - } - + } + offset <- function(y, w){ if (!is.null(mu)){ RET <- mu @@ -150,7 +150,7 @@ StudentTMu <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { } return(RET) } - + Family(ngradient = ngradient, risk = risk, loss = loss, offset=offset, response = function(f) f, @@ -164,7 +164,7 @@ StudentTMu <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { StudentTSigma <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { loss <- function(df, mu, y, f){ -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)))) + (df+1)/2 * log(1 + (y-mu)^2 / (df * exp(2 * f)))) } risk <- function(y, f, w = 1){ sum(w * loss(y = y, f = f, df = df, mu = mu)) @@ -188,7 +188,7 @@ StudentTSigma <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { } return(RET) } - + Family(ngradient = ngradient, risk = risk, loss = loss, offset=offset, response = function(f) exp(f), @@ -202,7 +202,7 @@ StudentTSigma <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { StudentTDf <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { loss <- function(sigma, mu,y, f){ -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))) + 0.5*f - (exp(f)+1)/2 * log(1 + (y-mu)^2 / (exp(f)*sigma^2))) } risk <- function(y, f, w = 1){ sum(w * loss(y = y, f = f, sigma = sigma, mu = mu)) @@ -210,7 +210,7 @@ StudentTDf <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { ngradient <- function(y, f, w = 1) { ngr <- exp(f)/2 * (digamma((exp(f)+1)/2)-digamma(exp(f)/2)) - 0.5 - (exp(f)/2 * log(1+ (y-mu)^2/(exp(f)*sigma^2)) - - (y-mu)^2/(sigma^2 + (y-mu)^2/exp(f)) * (exp(-f) +1)/2 ) + (y-mu)^2/(sigma^2 + (y-mu)^2/exp(f)) * (exp(-f) +1)/2 ) ngr <- stabilize_ngradient(ngr, w = w, stabilization) return(ngr) } @@ -228,7 +228,7 @@ StudentTDf <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { } return(RET) } - + Family(ngradient = ngradient, risk = risk, loss = loss, offset=offset, response = function(f) exp(f), @@ -242,7 +242,7 @@ StudentTDf <- function(mu = NULL, sigma = NULL, df = NULL, stabilization) { StudentTLSS <- function(mu = NULL, sigma = NULL, df = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { if ((!is.null(sigma) && sigma <= 0) || (!is.null(df) && df <= 0)) stop(sQuote("sigma"), " and ", sQuote("df"), " must be greater than zero") @@ -298,7 +298,7 @@ LogNormalMu <- function (mu = NULL, sigma = NULL, stabilization){ } return(RET) } - + Family(ngradient = ngradient, risk = risk, offset = offset, loss = loss, response = function(f) f, check_y = function(y) { @@ -339,7 +339,7 @@ LogNormalSigma <- function(mu = NULL, sigma = NULL, stabilization){ } return(RET) } - + Family(ngradient = ngradient, risk = risk, offset = offset, loss = loss, response = function(f) exp(f), check_y = function(y) { @@ -351,7 +351,7 @@ LogNormalSigma <- function(mu = NULL, sigma = NULL, stabilization){ } LogNormalLSS <- function(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { if ((!is.null(sigma) && sigma <= 0)) stop(sQuote("sigma"), " must be greater than zero") stabilization <- check_stabilization(stabilization) @@ -374,17 +374,17 @@ LogLogMu <- function (mu = NULL, sigma = NULL, stabilization){ loss <- function(sigma, y, f) { logfw <- function(pred) dlogis(pred, log = TRUE) - #pred - 2 * (1 + exp(pred)) + #pred - 2 * (1 + exp(pred)) logSw <- function(pred) plogis(pred, lower.tail = FALSE, log.p = TRUE) - #1/(1 + exp(pred)) + #1/(1 + exp(pred)) eta <- (log(y[,1]) - f)/sigma -y[,2] * (logfw(eta) - log(sigma)) - (1 - y[,2]) * logSw(eta) } - + risk <- function(y, f, w = 1) sum(w * loss(y = y, f = f, sigma = sigma)) - + ngradient <- function(y, f, w = 1) { eta <- (log(y[,1]) - f)/sigma nom <- (exp(-eta) + 1) @@ -403,8 +403,8 @@ LogLogMu <- function (mu = NULL, sigma = NULL, stabilization){ } return(RET) } - - + + Family(ngradient = ngradient, risk = risk, offset = offset, loss = loss, response = function(f) f, check_y = function(y) { @@ -419,10 +419,10 @@ LogLogSigma <- function (mu = NULL, sigma = NULL, stabilization){ loss <- function(mu, y, f) { logfw <- function(pred) dlogis(pred, log = TRUE) - #exp(pred)/(1 + exp(pred))^2 + #exp(pred)/(1 + exp(pred))^2 logSw <- function(pred) pnorm(pred, lower.tail = FALSE, log.p = TRUE) - #1/(1 + exp(pred)) + #1/(1 + exp(pred)) eta <- (log(y[,1]) - mu)/exp(f) -y[,2] * (logfw(eta) - f) - (1 - y[,2]) * logSw(eta) } @@ -446,7 +446,7 @@ LogLogSigma <- function (mu = NULL, sigma = NULL, stabilization){ } return(RET) } - + Family(ngradient = ngradient, risk = risk, offset = offset, loss = loss, response = function(f) exp(f), check_y = function(y) { @@ -458,13 +458,13 @@ LogLogSigma <- function (mu = NULL, sigma = NULL, stabilization){ } LogLogLSS <- function(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { if ((!is.null(sigma) && sigma <= 0)) stop(sQuote("sigma"), " must be greater than zero") stabilization <- check_stabilization(stabilization) Families(mu = LogLogMu(mu = mu, sigma = sigma, stabilization = stabilization), sigma = LogLogSigma(mu = mu, sigma = sigma, stabilization = stabilization), -# qfun = qLogLog, + # qfun = qLogLog, name = "Log-Log") } @@ -506,7 +506,7 @@ WeibullMu <- function (mu = NULL, sigma = NULL, stabilization){ } return(RET) } - + Family(ngradient = ngradient, risk = risk, offset = offset, loss = loss, response = function(f) f, check_y = function(y) { @@ -526,7 +526,7 @@ WeibullSigma <- function (mu = NULL, sigma = NULL, stabilization){ -exp(pred) eta <- (log(y[,1]) - mu)/exp(f) -y[,2] * (logfw(eta) - f) - (1 - y[,2]) * logSw(eta) - } + } risk <- function(y, f, w = 1) sum(w * loss(y = y, f = f, mu = mu)) ngradient <- function(y, f, w = 1) { @@ -548,7 +548,7 @@ WeibullSigma <- function (mu = NULL, sigma = NULL, stabilization){ } return(RET) } - + Family(ngradient = ngradient, risk = risk, offset = offset, loss = loss, response = function(f) exp(f), check_y = function(y) { @@ -560,7 +560,7 @@ WeibullSigma <- function (mu = NULL, sigma = NULL, stabilization){ } WeibullLSS <- function(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { if ((!is.null(sigma) && sigma <= 0)) stop(sQuote("sigma"), " must be greater than zero") stabilization <- check_stabilization(stabilization) @@ -577,7 +577,7 @@ qWeibull <- function(p, mu = 1, sigma = 1, lower.tail = TRUE, log.p = FALSE) { GaussianMu <- function(mu = NULL, sigma = NULL, stabilization){ - + loss <- function(sigma, y, f) -dnorm(x=y, mean=f, sd=sigma, log=TRUE) risk <- function(y, f, w = 1) { sum(w * loss(y = y, f = f, sigma = sigma)) @@ -587,7 +587,7 @@ GaussianMu <- function(mu = NULL, sigma = NULL, stabilization){ ngr <- stabilize_ngradient(ngr, w = w, stabilization) return(ngr) } - + offset <- function(y, w){ if (!is.null(mu)){ RET <- mu @@ -596,14 +596,14 @@ GaussianMu <- function(mu = NULL, sigma = NULL, stabilization){ } return(RET) } - + Family(ngradient = ngradient, risk = risk, loss = loss, response = function(f) f, offset=offset, name = "Normal distribution: mu(id link)") } GaussianSigma <- function(mu = NULL, sigma = NULL, stabilization){ - + loss <- function(y, f, mu) - dnorm(x=y, mean=mu, sd=exp(f), log=TRUE) risk <- function(y, f, w = 1) { sum(w * loss(y = y, f = f, mu = mu)) @@ -621,7 +621,7 @@ GaussianSigma <- function(mu = NULL, sigma = NULL, stabilization){ } return(RET) } - + Family(ngradient = ngradient, risk = risk, loss = loss, response = function(f) exp(f), offset=offset, name = "Normal distribution: sigma (log link)") @@ -629,7 +629,7 @@ GaussianSigma <- function(mu = NULL, sigma = NULL, stabilization){ GaussianLSS <- function(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { if ((!is.null(sigma) && sigma <= 0)) stop(sQuote("sigma"), " must be greater than zero") stabilization <- check_stabilization(stabilization) @@ -665,7 +665,7 @@ GammaMu <-function (mu = NULL, sigma = NULL, stabilization) { if (is.null(sigma)) sigma <<- mean(y)^2/(var(y)) RET <- optimize(risk, interval = c(0, max(y^2, na.rm = TRUE)), - y = y, w = w)$minimum + y = y, w = w)$minimum } return(RET) } @@ -703,7 +703,7 @@ GammaSigma <- function(mu = NULL, sigma = NULL, stabilization) { if (is.null(mu)) mu <<- mean(y) RET <- optimize(risk, interval = c(0, max(y, na.rm = TRUE)), - y = y, w = w)$minimum + y = y, w = w)$minimum } return(RET) } @@ -720,7 +720,7 @@ GammaSigma <- function(mu = NULL, sigma = NULL, stabilization) { } GammaLSS <- function (mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { if ((!is.null(sigma) && sigma <= 0)) stop(sQuote("sigma"), " must be greater than zero") if ((!is.null(mu) && mu <= 0)) @@ -743,13 +743,13 @@ qGamma <- function(p, mu = 0, sigma = 1, lower.tail = TRUE, log.p = FALSE) { BetaMu <- function(mu = NULL, phi = NULL, stabilization){ - + # loss is negative log-Likelihood, f is the parameter to be fitted with # logit link -> exp(f) = mu loss <- function(phi, y, f) { - 1 * (lgamma(phi) - lgamma(plogis(f) * phi) - - lgamma((1 - plogis(f)) * phi) + (plogis(f) * phi - 1) * log(y) + - ((1 - plogis(f)) * phi - 1) * log(1 - y)) + lgamma((1 - plogis(f)) * phi) + (plogis(f) * phi - 1) * log(y) + + ((1 - plogis(f)) * phi - 1) * log(1 - y)) } # risk is sum of loss risk <- function(y, f, w = 1) { @@ -769,11 +769,11 @@ BetaMu <- function(mu = NULL, phi = NULL, stabilization){ } else { if (is.null(phi)) - phi <<- mean(y) * (1 - mean(y)) /var(y) - 1 + phi <<- mean(y) * (1 - mean(y)) /var(y) - 1 ### look for starting value of f = qlogis(mu) in "interval" ### i.e. mu ranges from 0.000001 to 0.999999 RET <- optimize(risk, interval = c(qlogis(0.000001), qlogis(0.999999)), y = y, - w = w)$minimum + w = w)$minimum } return(RET) } @@ -792,14 +792,14 @@ BetaMu <- function(mu = NULL, phi = NULL, stabilization){ # Sub-Family for phi BetaPhi <- function(mu = NULL, phi = NULL, stabilization){ - + # loss is negative log-Likelihood, f is the parameter to be fitted with # log-link: exp(f) = phi loss <- function(mu, y, f) { #y <- (y * (length(y) - 1) + 0.5)/length(y) -1 * (lgamma(exp(f)) - lgamma(mu * exp(f)) - - lgamma((1 - mu) * exp(f)) + (mu * exp(f) - 1) * log(y) + - ((1 - mu) * exp(f) - 1) * log(1 - y)) + lgamma((1 - mu) * exp(f)) + (mu * exp(f) - 1) * log(y) + + ((1 - mu) * exp(f) - 1) * log(1 - y)) } # risk is sum of loss risk <- function(y, f, w = 1) { @@ -809,7 +809,7 @@ BetaPhi <- function(mu = NULL, phi = NULL, stabilization){ ngradient <- function(y, f, w = 1) { #y <- (y * (length(y) - 1) + 0.5)/length(y) ngr <- +1 * exp(f) * (mu * (qlogis(y) - (digamma(mu * exp(f)) - digamma((1 - mu) * exp(f)))) + - digamma(exp(f)) - digamma((1 - mu) * exp(f)) + log(1 - y)) + digamma(exp(f)) - digamma((1 - mu) * exp(f)) + log(1 - y)) ngr <- stabilize_ngradient(ngr, w = w, stabilization) return(ngr) } @@ -823,7 +823,7 @@ BetaPhi <- function(mu = NULL, phi = NULL, stabilization){ ### look for starting value of f = log(phi) in "interval" ### i.e. phi possibly ranges from 1e-10 to 1e10 RET <- optimize(risk, interval = c(log(1e-10), log(1e10)), - y = y, w = w)$minimum + y = y, w = w)$minimum } return(RET) } @@ -842,7 +842,7 @@ BetaPhi <- function(mu = NULL, phi = NULL, stabilization){ # families object for new distribution BetaLSS <- function (mu = NULL, phi = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { stabilization <- check_stabilization(stabilization) Families(mu = BetaMu(mu = mu, phi = phi, stabilization = stabilization), phi = BetaPhi(mu = mu, phi = phi, stabilization = stabilization), @@ -860,23 +860,23 @@ qBeta <- function(p, mu = 0, phi = 1, lower.tail = TRUE, log.p = FALSE) { # Zero-inflated Poisson model ZIPoLSS <- function(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { fam <- as.families(fname = "ZIP", mu = mu, sigma = sigma, stabilization = stabilization) - + fam$mu@name <- "Zero-inflated Poisson model: count data component" fam$sigma@name <- "Zero-inflated Poisson model: zero component" - + fam } # Zero-inflated negative binomial model ZINBLSS <- function(mu = NULL, sigma = NULL, nu = NULL, - stabilization = c("none", "MAD")) { + stabilization = c("none", "MAD", "L2")) { fam <- as.families(fname = "ZINBI", mu = mu, sigma = sigma, nu = nu, stabilization = stabilization) - + fam$mu@name <- "Zero-inflated negative binomial model: location parameter for count data component" fam$sigma@name <- "Zero-inflated negative binomial model: scale parameter for count data component" fam$nu@name <- "Zero-inflated negative binomial model: zero component" - + fam } diff --git a/R/helpers.R b/R/helpers.R index 80b1f85..d7a593c 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -1,11 +1,11 @@ ## helper functions check <- function(what, what_char, names) { - + errormsg <- paste0(sQuote(what_char), " can be either a scalar, a (named) vector or a (named) list", - " of ", what_char, " values with same names as ", sQuote("families"), "in ", - sQuote("boost_control")) - + " of ", what_char, " values with same names as ", sQuote("families"), "in ", + sQuote("boost_control")) + if (is.list(what)) { if (is.null(names(what)) && length(what) == length(names)) names(what) <- names @@ -28,7 +28,7 @@ check <- function(what, what_char, names) { what <- what[names] ## sort in order of families } } - + return(what) } @@ -41,7 +41,7 @@ get_data <- function(x, which = NULL) { data <- try(sapply(which, get, env = parent.frame(2)), silent = TRUE) if (inherits(data, "try-error")) - stop("No data set found.") + stop("No data set found.") data <- as.data.frame(data) } return(data) @@ -119,20 +119,20 @@ rescale_weights <- function(w) { ## helper function in a modified version based on mboost_2.2-3 ## print trace of boosting iterations -do_trace <- function(current, mstart, risk, - linebreak = options("width")$width / 2, mstop = 1000) { - current <- current - mstart +do_trace <- function(current, risk, mstart, + linebreak = round(options("width")$width / 2), mstop = 1000) { + + current <- current - mstart + if (current != mstop) { - if ((current - 1) %/% linebreak == (current - 1) / linebreak) { + if (current %% linebreak == 1) { mchr <- formatC(current + mstart, format = "d", width = nchar(mstop) + 1, big.mark = "'") cat(paste("[", mchr, "] ",sep = "")) - } else { - if ((current %/% linebreak != current / linebreak)) { - cat(".") - } else { - cat(" -- risk:", risk[current + mstart], "\n") - } + } + cat(".") + if (current %% linebreak == 0) { + cat(" -- risk:", risk[current + mstart], "\n") } } else { cat("\nFinal risk:", risk[current + mstart], "\n") @@ -166,11 +166,17 @@ stabilize_ngradient <- function(ngr, w = 1, stabilization) { div <- ifelse(div < 0.0001, 0.0001, div) ngr <- ngr / div } + if (stabilization == "L2") { + div <- sqrt(weighted.mean(ngr^2, w =w, na.rm = TRUE)) + div <- ifelse(div < 1e-04, 1e-04, div) + div <- ifelse(div > 1e+04, 1e+04, div) + ngr <- ngr / div + } ngr } -check_stabilization <- function(stabilization = c("none", "MAD")) { +check_stabilization <- function(stabilization = c("none", "MAD", "L2")) { stabilization <- match.arg(stabilization) ## check if old stabilization interface is used and issue a warning if (getOption("gamboostLSS_stab_ngrad")) { @@ -178,7 +184,7 @@ check_stabilization <- function(stabilization = c("none", "MAD")) { " is deprecated.\n", "Use argument ", sQuote("stabilization"), " in the fitting family. See ?Families for details.") if (stabilization == "none") - warning(sQuote("stabilization"), " is set to ", dQuote("MAD")) + warning(sQuote("stabilization"), " is set to ", dQuote("MAD")) } stabilization } diff --git a/R/mboostLSS.R b/R/mboostLSS.R index 833aba1..3e9fd19 100644 --- a/R/mboostLSS.R +++ b/R/mboostLSS.R @@ -4,46 +4,67 @@ ### (glm/gam/m/black)boostLSS functions mboostLSS <- function(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...){ + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...){ + cl <- match.call() if(is.null(cl$families)) cl$families <- families + method <- match.arg(method) + fit <- mboostLSS_fit(formula = formula, data = data, families = families, control = control, weights = weights, ..., - fun = mboost, funchar = "mboost", call = cl) + fun = mboost, funchar = "mboost", call = cl, + method = method) return(fit) } glmboostLSS <- function(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...){ + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...){ + cl <- match.call() if(is.null(cl$families)) cl$families <- families + method <- match.arg(method) + fit <- mboostLSS_fit(formula = formula, data = data, families = families, control = control, weights = weights, ..., - fun = glmboost, funchar = "glmboost", call = cl) + fun = glmboost, funchar = "glmboost", call = cl, + method = method) + return(fit) } gamboostLSS <- function(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...){ + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...){ + cl <- match.call() if(is.null(cl$families)) cl$families <- families + method <- match.arg(method) + fit <- mboostLSS_fit(formula = formula, data = data, families = families, control = control, weights = weights, ..., - fun = gamboost, funchar = "gamboost", call = cl) + fun = gamboost, funchar = "gamboost", call = cl, + method = method) + return(fit) } blackboostLSS <- function(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...){ + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...){ cl <- match.call() if(is.null(cl$families)) cl$families <- families + method <- match.arg(method) + fit <- mboostLSS_fit(formula = formula, data = data, families = families, control = control, weights = weights, ..., - fun = blackboost, funchar = "blackboost", call = cl) + fun = blackboost, funchar = "blackboost", call = cl, + method = method) return(fit) } @@ -51,7 +72,8 @@ blackboostLSS <- function(formula, data = list(), families = GaussianLSS(), ### work horse for fitting (glm/gam/m/black)boostLSS models mboostLSS_fit <- function(formula, data = list(), families = GaussianLSS(), control = boost_control(), weights = NULL, - fun = mboost, funchar = "mboost", call = NULL, ...){ + fun = mboost, funchar = "mboost", call = NULL, + method, ...){ if (length(families) == 0) stop(sQuote("families"), " not specified") @@ -70,8 +92,8 @@ mboostLSS_fit <- function(formula, data = list(), families = GaussianLSS(), if (length(unique(ynames)) > 1) warning("responses differ for the components") response <- lapply(formula, function(fm) - eval(as.expression(fm[[2]]), envir = data, - enclos = environment(fm))) + eval(as.expression(fm[[2]]), envir = data, + enclos = environment(fm))) #response <- eval(as.expression(formula[[1]][[2]]), envir = data, # enclos = environment(formula[[1]])) } else { @@ -85,8 +107,10 @@ mboostLSS_fit <- function(formula, data = list(), families = GaussianLSS(), } mstop <- mstoparg <- control$mstop - control$mstop <- 1 - mstop <- check(mstop, "mstop", names(families)) + control$mstop <- 0 + if (method == "cyclic") + mstop <- check(mstop, "mstop", names(families)) + nu <- control$nu nu <- check(nu, "nu", names(families)) @@ -133,11 +157,9 @@ mboostLSS_fit <- function(formula, data = list(), families = GaussianLSS(), } for (j in mods){ ## update value of nuisance parameters in families - ## use response(fitted()) as this is much quicker than - ## fitted(, type = response) as the latter uses predict() for (k in mods[-j]){ if (!is.null(fit[[k]])) - assign(names(fit)[k], families[[k]]@response(fitted(fit[[k]])), + assign(names(fit)[k], fitted(fit[[k]], type = "response"), environment(families[[j]]@ngradient)) } ## use appropriate nu for the model @@ -148,120 +170,245 @@ mboostLSS_fit <- function(formula, data = list(), families = GaussianLSS(), control=control, weights = w, ...)) } - if (trace) - do_trace(current = 1, mstart = 0, - mstop = max(mstop), - risk = fit[[length(fit)]]$risk()) - ### set up a function for iterating boosting steps - iBoost <- function(niter) { + iBoost <- function(niter, method) { + start <- sapply(fit, mstop) - mvals <- vector("list", length(niter)) - for (j in 1:length(niter)){ - mvals[[j]] <- rep(start[j] + niter[j], max(niter)) - if (niter[j] > 0) - mvals[[j]][1:niter[j]] <- (start[j] + 1):(start[j] + niter[j]) + # initialize combined_risk + combined_risk <- NA + + if (method == "noncyclic") { + ### noncyclical fitting ### + # this is the case for boosting from the beginning + if (is.null(attr(fit, "combined_risk")) | niter == 0) { + combined_risk <- vapply(fit, risk, numeric(1)) + } + + best <- which(names(fit) == tail(names(combined_risk), 1)) + } else { + ### cyclical fitting ### + mvals <- vector("list", length(niter)) + for (j in 1:length(niter)) { + mvals[[j]] <- rep(start[j] + niter[j], max(niter)) + if (niter[j] > 0) + mvals[[j]][1:niter[j]] <- (start[j] + 1):(start[j] + niter[j]) + } } ENV <- lapply(mods, function(j) environment(fit[[j]]$subset)) + ## main loop starts here ## for (i in 1:max(niter)){ - for (j in mods){ - ## update value of nuisance parameters - ## use response(fitted()) as this is much quicker than - ## fitted(, type = response) as the latter uses predict() - for (k in mods[-j]) - assign(names(fit)[k], families[[k]]@response(fitted(fit[[k]])), - environment(get("ngradient", environment(fit[[j]]$subset)))) - ## update value of u, i.e. compute ngradient with new nuisance parameters + if (method == "noncyclic") { + ### noncyclical fitting ### - ENV[[j]][["u"]] <- ENV[[j]][["ngradient"]](ENV[[j]][["y"]], ENV[[j]][["fit"]], - ENV[[j]][["weights"]]) - # same as: - # evalq(u <- ngradient(y, fit, weights), environment(fit[[j]]$subset)) + ## update value of nuisance parameters + ## use response(fitted()) as this is much quicker than fitted(, type = response) + for( k in mods[-best]) { + assign(names(fit)[best], families[[best]]@response(fitted(fit[[best]])), + environment(get("ngradient", environment(fit[[k]]$subset)))) + } + + risks <- numeric(length(fit)) + for(b in 1:length(fit)){ + st <- mstop(fit[[b]]) + mstop(fit[[b]]) = st + 1 + risks[b] <- tail(risk(fit[[b]]), 1) + #evalq({riskfct(y, fit, weights)}, envir = ENV[[b]]) + #risks[b] <- ENV[[b]][["riskfct"]](ENV[[b]][["y"]], ENV[[b]][["fit"]], ENV[[b]][["weights"]]) + fit[[b]][st] + + ## fit[[b]][st] is not enough to reduce the model back to beginning, so + ## so all these values have to be reduced, so that they are calculated + ## correctly the next time + evalq({ + xselect <- xselect[seq_len(mstop)]; + mrisk <- mrisk[seq_len(mstop + 1)]; + ens <- ens[seq_len(mstop)]; + nuisance <- nuisance[seq_len(mstop)] + }, + environment(fit[[b]]$subset)) + } + + best <- which.min(risks) - ## update j-th component to "m-th" boosting step - fit[[j]][mvals[[j]][i]] + ## update value of u, i.e. compute ngradient with new nuisance parameters + evalq({u <- ngradient(y, fit, weights)}, ENV[[best]]) + #ENV[[best]][["u"]] <- ENV[[best]][["ngradient"]](ENV[[best]][["y"]], ENV[[best]][["fit"]], ENV[[best]][["weights"]]) + + ## update selected component by 1 + fit[[best]][mstop(fit[[best]]) + 1] + + ## update risk list + combined_risk[(length(combined_risk) + 1)] <- tail(risk(fit[[best]]), 1) + names(combined_risk)[length(combined_risk)] <- names(fit)[best] + combined_risk <<- combined_risk + + } else { + ### cyclical fitting ### + for (j in mods){ + ## 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)))) + ## update value of u, i.e. compute ngradient with new nuisance parameters + + ENV[[j]][["u"]] <- ENV[[j]][["ngradient"]](ENV[[j]][["y"]], ENV[[j]][["fit"]], + ENV[[j]][["weights"]]) + # same as: + # evalq(u <- ngradient(y, fit, weights), environment(fit[[j]]$subset)) + + ## update j-th component to "m-th" boosting step + fit[[j]][mvals[[j]][i]] + } } + if (trace){ - ## which is the current risk? rev() needed to get the last - ## list element with maximum length - whichRisk <- names(which.max(rev(lapply(lapply(fit, function(x) x$risk()), length)))) - do_trace(current = max(sapply(mvals, function(x) x[i])), - mstart = ifelse(firstRun, 0, max(start)), - mstop = ifelse(firstRun, max(niter) + 1, max(niter)), - risk = fit[[whichRisk]]$risk()) + if (method == "noncyclic") { + do_trace(current = length(combined_risk) - length(fit), mstart = sum(start), + risk = combined_risk[-(1:length(fit))], mstop = niter) + } else { + ## which is the current risk? rev() needed to get the last + ## list element with maximum length + whichRisk <- names(which.max(rev(lapply(fit, function(x) length(risk(x)))))) + do_trace(current = max(sapply(mvals, function(x) x[i])), + mstart = max(start), + mstop = max(niter), + risk = risk(fit[[whichRisk]])[-1]) + } + } } return(TRUE) } - if (any(mstop > 1)){ - ## actually go for initial mstop iterations! - firstRun <- TRUE - tmp <- iBoost(mstop - 1) - } + if (any(mstop > 0)) + tmp <- iBoost(mstop, method = method) - firstRun <- FALSE - - class(fit) <- c(paste(funchar, "LSS", sep=""), "mboostLSS") + class(fit) <- c(paste0(funchar, "LSS"), "mboostLSS") + if(method != "cyclic"){ + class(fit) <- c("nc_mboostLSS", class(fit)) + } ### update to a new number of boosting iterations mstop ### i <= mstop means less iterations than current ### i > mstop needs additional computations ### updates take place in THIS ENVIRONMENT, ### some models are CHANGED! - attr(fit, "subset") <- function(i) { + if(method == "cyclic") { + attr(fit, "subset") <- function(i) { + + i <- check(i, "mstop", names(families)) + + msf <- mstop(fit) + niter <- i - msf + if (all(niter == 0)) { + ## make nothing happen with the model + minStart <- max(msf) + } else { + minStart <- min(msf[niter != 0], i[niter != 0]) + } + + ## check if minStart bigger than mstop of parameters that are not touched + #if (length(msf[niter == 0]) > 0 && minStart < min(msf[niter == 0])) + #minStart <- min(msf[niter == 0]) + + ## reduce models first (when necessary) + if (any(msf > minStart)){ + cf <- class(fit) + class(fit) <- "list" ## needed to use [] operator for lists + + #cat("processed parameters: ", paste(names(fit[msf > minStart]), + # collapse = ", "), "\n") + + lapply(fit[msf > minStart], + function(obj) obj$subset(minStart)) + + ## remove additional boosting iterations from environments + lapply(fit[msf > minStart], function(obj){ + evalq({xselect <- xselect[seq_len(mstop)]; + mrisk <- mrisk[seq_len(mstop + 1)]; + ens <- ens[seq_len(mstop)]; + nuisance <- nuisance[seq_len(mstop)]}, + environment(obj$subset)) + }) + + class(fit) <- cf + if (trace) + cat("Model first reduced to mstop = ", minStart, ".\n", + "Now continue ...\n", sep ="") + } - i <- check(i, "mstop", names(families)) + ## now increase models (when necessary) + if (any(i > minStart)){ + ## set negative values to 0 + ## (only applicable if some parameters do not need to be touched + inc <- ifelse(i - minStart > 0, i - minStart, 0) + tmp <- iBoost(inc, method = method) + } - msf <- mstop(fit) - niter <- i - msf - if (all(niter == 0)) { - ## make nothing happen with the model - minStart <- max(msf) - } else { - minStart <- min(msf[niter != 0], i[niter != 0]) + mstop <<- i } + } + else { + attr(fit, "subset") <- function(i) { + msf <- sum(mstop(fit)) + niter <- i - msf + + ## check if minStart bigger than mstop of parameters that are not touched + #if (length(msf[niter == 0]) > 0 && minStart < min(msf[niter == 0])) + #minStart <- min(msf[niter == 0]) + + ## reduce models first (when necessary) + if (niter < 0 ){ + cf <- class(fit) + d = length(fit) + class(fit) <- "list" ## needed to use [] operator for lists + #cat("processed parameters: ", paste(names(fit[msf > minStart]), + # collapse = ", "), "\n") + #reduce the combined risk values + combined_risk <- attr(fit, "combined_risk")() + combined_risk <<- attr(fit, "combined_risk")()[seq_len(i + d)] + new_stop_value <- table(names(attr(fit, "combined_risk")())) - 1 + + + for(o in names(new_stop_value)){ + fit[[o]]$subset(new_stop_value[o]) + } + + ## remove additional boosting iterations from environments + lapply(fit, function(obj){ + evalq({xselect <- xselect[seq_len(mstop)]; + mrisk <- mrisk[seq_len(mstop + 1)]; + ens <- ens[seq_len(mstop)]; + nuisance <- nuisance[seq_len(mstop)]}, + environment(obj$subset)) + }) + + #update ALL nuisance parameters to last update + ENV <- lapply(mods, function(j) environment(fit[[j]]$subset)) + for(j in names(new_stop_value)){ + for( k in setdiff(names(new_stop_value), j)){ + assign(k, families[[k]]@response(fitted(fit[[k]])), + environment(get("ngradient", environment(fit[[j]]$subset)))) + } + } + + for(k in mods){ + evalq(u <- get("ngradient")(get("y"), fit, weights), ENV[[k]]) + } + + + class(fit) <- cf + } - ## check if minStart bigger than mstop of parameters that are not touched - #if (length(msf[niter == 0]) > 0 && minStart < min(msf[niter == 0])) - #minStart <- min(msf[niter == 0]) - - ## reduce models first (when necessary) - if (any(msf > minStart)){ - cf <- class(fit) - class(fit) <- "list" ## needed to use [] operator for lists - - #cat("processed parameters: ", paste(names(fit[msf > minStart]), - # collapse = ", "), "\n") - - lapply(fit[msf > minStart], - function(obj) obj$subset(minStart)) - - ## remove additional boosting iterations from environments - lapply(fit[msf > minStart], function(obj){ - evalq({xselect <- xselect[1:mstop]; - mrisk <- mrisk[1:mstop]; - ens <- ens[1:mstop]; - nuisance <- nuisance[1:mstop]}, - environment(obj$subset)) - }) - - class(fit) <- cf - if (trace) - cat("Model first reduced to mstop = ", minStart, ".\n", - "Now continue ...\n", sep ="") - } + ## now increase models (when necessary) + else if (niter > 0){ + tmp <- iBoost(niter, method = method) + } - ## now increase models (when necessary) - if (any(i > minStart)){ - ## set negative values to 0 - ## (only applicable if some parameters do not need to be touched - inc <- ifelse(i - minStart > 0, i - minStart, 0) - tmp <- iBoost(inc) + mstop <<- i } - - mstop <<- i } ## make call in submodels nicer: @@ -294,11 +441,15 @@ mboostLSS_fit <- function(formula, data = list(), families = GaussianLSS(), mboostLSS_fit(formula = formula, data = data, families = eval(call[["families"]]), weights = weights, control = control, fun = fun, funchar = funchar, - call = call, oobweights = oobweights) + call = call, oobweights = oobweights, + method = method) } attr(fit, "control") <- control attr(fit, "call") <- call attr(fit, "data") <- data attr(fit, "families") <- families + if(method != "cyclic") + attr(fit, "combined_risk") <- function() combined_risk + return(fit) } diff --git a/R/methods.R b/R/methods.R index 220d492..ec52214 100644 --- a/R/methods.R +++ b/R/methods.R @@ -1,7 +1,7 @@ ## Methods "[.mboostLSS" <- function(x, i, return = TRUE, ...) { - stopifnot((length(i) == 1 | length(i) == length(x)) && i > 0) + stopifnot(length(i) == 1 | length(i) == length(x)) attr(x, "subset")(i) if (return) return(x) invisible(NULL) @@ -13,9 +13,9 @@ coef.mboostLSS <- function(object, which = NULL, if (is.character(parameter)) parameter <- extract_parameter(object, parameter) RET <- lapply(parameter, function(i, object) - coef(object[[i]], which = which, - aggregate = aggregate, ...), - object = object) + coef(object[[i]], which = which, + aggregate = aggregate, ...), + object = object) if (length(RET) == 1) RET <- RET[[1]] return(RET) @@ -31,15 +31,15 @@ coef.glmboostLSS <- function(object, which = NULL, risk.mboostLSS <- function(object, merge = FALSE, parameter = names(object), ...){ if (is.character(parameter)) parameter <- extract_parameter(object, parameter) - + lo <- length(unique(mstop(object))) if (merge) { get_rsk <- function(i, object) { - mmo <- max(mstop(object)) + 1 - rsk <- object[[i]]$risk() + mmo <- max(mstop(object)) + 1 + rsk <- object[[i]]$risk() if (length(rsk) != mmo) { - if (length(rsk) > mmo) - stop("risk cannot contain more entries than mstop") + if (length(rsk) > mmo) + stop("risk cannot contain more entries than mstop") rsk <- c(rsk, rep(NA, mmo - length(rsk))) } rsk @@ -68,6 +68,13 @@ mstop.mboostLSS <- function(object, parameter = names(object), ...){ names(RET) <- names(object)[parameter] if (length(RET) == 1) RET <- RET[[1]] + # change mstop for noncyclical fitting to scalar value with attribute of partitions + if (inherits(object, "nc_mboostLSS")) { + partitions <- RET + RET <- sum(RET) + attr(RET, "partitions") <- partitions + } + return(RET) } @@ -82,16 +89,74 @@ mstop.oobag <- function(object, parameter = names(object), ...){ return(RET) } -selected.mboostLSS <- function(object, parameter = names(object), ...){ +selected.mboostLSS <- function(object, merge = FALSE, parameter = names(object), + ...){ if (is.character(parameter)) parameter <- extract_parameter(object, parameter) - RET <- lapply(parameter, function(i, object) - selected(object[[i]]), - object = object) - names(RET) <- names(object)[parameter] - if (length(RET) == 1) - RET <- RET[[1]] - return(RET) + + #merge is different for noncyclical fitting + if (merge) { + if (inherits(object, "nc_mboostLSS")){ + RET <- names(attr(object, "combined_risk")()) + + for(p in names(parameter)){ + RET[RET == p] <- object[[p]]$xselect() + } + RET <- as.numeric(RET) + names(RET) <- names(attr(object, "combined_risk")()) + return(RET) + } + else { + get_sel <- function(i, object) { + mmo <- max(mstop(object)) + sel <- object[[i]]$xselect() + if (length(sel) != mmo) { + sel <- c(sel, rep(NA, mmo - length(sel))) + } + sel + } + + RET <- sapply(parameter, get_sel, + object = object) + + RET <- as.vector(t(RET)) + names(RET) <- rep(names(parameter), mstop(object)[1]) + lo <- length(unique(mstop(object))) + if (lo != 1) + RET <- RET[!is.na(RET)] + return(RET) + } + + } + else { + RET <- lapply(parameter, function(i, object) + selected(object[[i]]), + object = object) + names(RET) <- names(object)[parameter] + if (length(RET) == 1) + RET <- RET[[1]] + return(RET) + } + +} + +selected.stabsel_mboostLSS <- function(object, parameter = NULL, ...) { + ## extract parameters + param_count <- table(gsub(".*\\.", "", rownames(object$phat))) + ## set up a named list for the results + res <- vector("list", length(param_count)) + names(res) <- names(param_count) + ## subtract the number of base-learners of the first components + param_count <- c(0, cumsum(param_count)) + for (i in 1:length(res)) { + idx <- gsub(".*\\.", "", names(object$selected)) == names(res)[[i]] + res[[i]] <- object$selected[idx] - param_count[i] + } + if (!is.null(parameter)) + res <- res[parameter] + if (length(res) == 1) + res <- res[[1]] + return(res) } @@ -100,8 +165,8 @@ plot.glmboostLSS <- function(x, main = names(x), parameter = names(x), if (is.character(parameter)) parameter <- extract_parameter(x, parameter) lapply(parameter, function(i, x, main, off2int, ...) - plot(x[[i]], main = main[[i]], off2int = off2int, ...), - x = x, main = main, off2int = off2int, ...) + plot(x[[i]], main = main[[i]], off2int = off2int, ...), + x = x, main = main, off2int = off2int, ...) invisible(coef(x, aggregate = "cumsum", off2int = off2int)) } @@ -110,8 +175,8 @@ plot.gamboostLSS <- function(x, main = names(x), parameter = names(x), ...){ if (is.character(parameter)) parameter <- extract_parameter(x, parameter) RET <- lapply(parameter, function(i, x, main, ...) - plot(x[[i]], main = main[[i]], ...), - x = x, main = main, ...) + plot(x[[i]], main = main[[i]], ...), + x = x, main = main, ...) if (any(sapply(RET, class) == "trellis")) { return(RET) } else { @@ -122,29 +187,29 @@ plot.gamboostLSS <- function(x, main = names(x), parameter = names(x), ...){ plot.predint <- function(x, main = "Marginal Prediction Interval(s)", xlab = NULL, ylab = NULL, lty = c("solid", "dashed"), lcol = c("black", "black"), log = "", ...) { - + pi <- attr(x, "pi") which <- attr(x, "which") rawdata <- attr(x, "rawdata") - + if (length(lty) != length(pi) + 1) lty <- c(lty, rep(tail(lty, 1), (length(pi) + 1) - length(lty))) if (length(lcol) != length(pi) + 1) lcol <- c(lcol, rep(tail(lcol, 1), (length(pi) + 1) - length(lcol))) - + if (is.null(xlab)) xlab <- which if (is.null(ylab)) ylab <- "prediction" - + plot(rawdata$x, rawdata$y, pch = 20, col = rgb(0.5, 0.5, 0.5, 0.5), xlab = xlab, ylab = ylab, main = main, log = log, ...) - + lines(x[, which], x$"Prediction (Median)", lty = lty[1], col = lcol[1], ...) - + for (i in seq_along(pi)) { lines(x[, which], x[, paste0(pi[i] * 100, "% PI (lower)")], lty = lty[i + 1], col = lcol[i + 1], ...) @@ -156,20 +221,20 @@ plot.predint <- function(x, main = "Marginal Prediction Interval(s)", PI <- predint <- function(x, which, pi = 0.9, newdata = NULL, ...) { qfun <- get_qfun(x) - + if (length(which) != 1 || !is.character(which)) stop("Please specify the variable for the marginal prediction interval.") - + var <- get_data(x, which = which) if (ncol(var) > 1 || is.factor(var)) stop("Prediction intervals are currently only implemented for ", "base-learners of one numeric variable") - + pred_vars <- lapply(x, extract, what = "variable.names") pred_vars <- unique(unlist(pred_vars)) if ("(Intercept)" %in% pred_vars) pred_vars <- pred_vars[pred_vars != "(Intercept)"] - + if (is.null(newdata)) { tmp <- get_data(x, which = pred_vars) i <- grepl(which, names(tmp)) @@ -187,20 +252,20 @@ PI <- predint <- function(x, which, pi = 0.9, newdata = NULL, ...) { if (nrow(unique(newdata[, !i])) != 1) warning("All variables but", sQuote("which"), "should be constant") } - + newdata[, names(x)] <- predict(x, newdata = newdata, type = "response") newdata$"Prediction (Median)" <- do.call(qfun, args = c(p = 0.5, - newdata[, names(x)])) - + newdata[, names(x)])) + for (i in seq_along(pi)) { newdata[, paste0(pi[i] * 100, "% PI (lower)")] <- do.call(qfun, - args = c(p = (1 - pi[i])/2, newdata[, names(x)])) + args = c(p = (1 - pi[i])/2, newdata[, names(x)])) newdata[, paste0(pi[i] * 100, "% PI (upper)")] <- do.call(qfun, - args = c(p = 1 - (1 - pi[i])/2, newdata[, names(x)])) + args = c(p = 1 - (1 - pi[i])/2, newdata[, names(x)])) } # drop predictions of parameters newdata <- newdata[, !names(newdata) %in% names(x)] - + class(newdata) <- c("predint", "data.frame") attr(newdata, "pi") <- pi attr(newdata, "which") <- which @@ -217,8 +282,11 @@ print.mboostLSS <- function(x, ...){ cat("\n") if (!is.null(attr(x, "call"))) cat("Call:\n", deparse(attr(x, "call")), "\n\n", sep = "") + m <- mstop(x) + if (inherits(x, "nc_mboostLSS")) + m <- attr(m, "partitions") cat("Number of boosting iterations (mstop): ", - paste(names(mstop(x)), mstop(x), sep = " = ", collapse = ", "), "\n") + paste(names(x), m, sep = " = ", collapse = ", "), "\n") nus <- sapply(x, function(xi) xi$control$nu) cat("Step size: ", paste(names(nus), nus, sep = " = ", collapse = ", "), "\n\n") @@ -232,7 +300,7 @@ fitted.mboostLSS <- function(object, parameter = names(object), ...){ if (is.character(parameter)) parameter <- extract_parameter(object, parameter) myApply(parameter, function(i, mod, ...) fitted(mod[[i]], ...), - mod = object, ...) + mod = object, ...) } @@ -244,9 +312,9 @@ predict.mboostLSS <- function(object, newdata = NULL, if (is.character(parameter)) parameter <- extract_parameter(object, parameter) myApply(parameter, function(i, mod, ...) - predict(mod[[i]], newdata = newdata, type = type, which = which, - aggregate = aggregate, ...), - mod = object, ...) + predict(mod[[i]], newdata = newdata, type = type, which = which, + aggregate = aggregate, ...), + mod = object, ...) } update.mboostLSS <- function(object, weights, oobweights = NULL, @@ -272,17 +340,18 @@ summary.mboostLSS <- function(object, ...) { cat("\n") if (!is.null(attr(object, "call"))) cat("Call:\n", deparse(attr(object, "call")), "\n\n", sep = "") + m <- mstop(object) + if (inherits(object, "nc_mboostLSS")) + m <- attr(m, "partitions") cat("Number of boosting iterations (mstop): ", - paste(names(mstop(object)), mstop(object), - sep = " = ", collapse = ", "), - "\n") + paste(names(object), m, sep = " = ", collapse = ", "), "\n") nus <- sapply(object, function(xi) xi$control$nu) cat("Step size: ", paste(names(nus), nus, sep = " = ", collapse = ", "), "\n\n") - + cat("Families:\n") lapply(object, function(xi) show(xi$family)) - + if (inherits(object, "glmboostLSS")) { cat("Coefficients:\n") cf <- coef(object, off2int = TRUE) @@ -292,7 +361,7 @@ summary.mboostLSS <- function(object, ...) { cat("\n") } } - + cat("Selection frequencies:\n") for (i in 1:length(object)) { cat("Parameter ", names(object)[i], ":\n", sep = "") @@ -307,6 +376,202 @@ summary.mboostLSS <- function(object, ...) { invisible(object) } +stabsel.mboostLSS <- function(x, cutoff, q, PFER, + mstop = NULL, + folds = subsample(model.weights(x), B = B), + B = ifelse(sampling.type == "MB", 100, 50), + assumption = c("unimodal", "r-concave", "none"), + sampling.type = c("SS", "MB"), + papply = mclapply, verbose = TRUE, FWER, eval = TRUE, ...) { + + cll <- match.call() + p <- sum(sapply(x, function(obj) length(variable.names(obj)))) + n <- if(inherits(x, "FDboostLSS")) { + x[[1]]$ydim[1] + } else { + nrow(attr(x, "data")) + } + + + ## extract names of base-learners (and add paramter name) + nms <- lapply(x, function(obj) variable.names(obj)) + nms <- lapply(1:length(nms), function(i) + paste(nms[[i]], names(nms)[i], sep = ".")) + + sampling.type <- match.arg(sampling.type) + if (sampling.type == "MB") + assumption <- "none" + else + assumption <- match.arg(assumption) + + + if (inherits(x, "nc_mboostLSS")) { + ## check mstop + if (is.null(mstop)) + mstop <- sum(mstop(x)) + if (length(mstop) != 1 | mstop %% 1 != 0 | mstop < length(x)) { + stop(sQuote("mstop"), " has to be an integer larger than ", + length(x)) + } + mstop_min <- length(x) + } + else { + if (is.null(mstop)) + mstop <- mstop(x) + mstop <- check(mstop, "mstop", names(x)) + mstop_min <- 1 + } + + if (length(unique(mstop)) != 1) + warning("Usually one should use the same ", sQuote("mstop"), + " value for all components for stability selection.") + + if (verbose) + cat("Run stabsel ") + + ## set mstop = 1 to speed things up + x <- update(x, weights = model.weights(x), mstop = mstop_min) + + ## define the fitting function (args.fitfun is not used but needed for + ## compatibility with run_stabsel + fit_model <- function(i, folds, q, args.fitfun) { + if (verbose) + cat(".") + ## start by fitting 1 step in each component + mod <- update(x, weights = folds[, i], mstop = mstop_min) + ## make sure dispatch works correctly + class(mod) <- class(x) + xs <- selected(mod) + nsel <- length(mod) + ## now update model until we obtain q different base-learners altogether + for (m in mstop_min:max(mstop)) { + if (nsel >= q) + break + + mstop(mod) <- m + xs <- selected(mod) + nsel <- sum(sapply(xs, function(selection) + length(unique(selection)))) + if (nsel >= q) + break + } + #this changes nothing for method = "cyclic" but fixes mstop for method = "noncyclic" + mstop <- check(mstop, "mstop", names(x)) + ## complete paths + if (any(sapply(xs, length) < mstop)) { + for (j in 1:length(xs)) { + start <- length(xs[[j]]) + 1 + xs[[j]][start:mstop[j]] <- xs[[j]][1] + } + } + + selected <- lapply(xs, unique) + ret <- lapply(1:length(selected), function(i) { + res <- logical(length(nms[[i]])) + names(res) <- nms[[i]] + res[selected[[i]]] <- TRUE + res + }) + ret <- unlist(ret) + + ## compute selection paths + #merging for method cyclic + if(!inherits(x, "nc_mboostLSS")){ + sequences <- lapply(1:length(xs), function(i) { + res <- matrix(FALSE, nrow = length(nms[[i]]), ncol = mstop[[i]]) + rownames(res) <- nms[[i]] + for (j in 1:mstop[[i]]) + res[xs[[i]][j], j:mstop[[i]]] <- TRUE + res + }) + + if (any(mstop < max(mstop))) + stop("simply add the last column to the smaller matrices") + ## now merge sequences + for (i in 1:ncol(sequences[[1]])) { + for (j in 1:length(sequences)) { + if (i == 1) { + if (j == 1) { + other_params <- rep(FALSE, sum(sapply(sequences, nrow)[-1])) + sequence <- matrix(c(sequences[[i]][, j], + other_params)) + } else { + tmp <- unlist(lapply(sequences[1:j], function(x) x[, + i])) + other_params <- rep(FALSE, sum(sapply(sequences, + nrow)[-(1:j)])) + tmp <- c(tmp, other_params) + sequence <- cbind(sequence, tmp) + } + } else { + if (j < length(sequences)) { + tmp <- unlist(c(lapply(sequences[1:j], function(x) x[, i]), + lapply(sequences[(j+1):length(sequences)], + function(x) x[, i - 1]))) + } else { + tmp <- unlist(lapply(sequences[1:j], function(x) x[, + i])) + } + sequence <- cbind(sequence, tmp) + } + } + } + } + else { + sequence <- matrix(FALSE, nrow = p, ncol = mstop[1]) + rownames(sequence) <- unlist(nms) + + sel <- selected(mod, merge = TRUE) + + for (i in names(mod)) { + varnames <- variable.names(mod[[i]]) + for(j in seq_along(varnames)){ + pos <- which(names(sel) == i & sel == j) + if(length(pos) > 0) + sequence[paste(varnames[j], i, sep = "."), min(pos):mstop[1]] <- TRUE + } + + } + + } + + colnames(sequence) <- 1:ncol(sequence) + ret <- list(selected = ret, path = sequence) + ## was mstop to small? + attr(ret, "violations") <- ifelse(sum(ret$selected) < q, TRUE, FALSE) + return(ret) + } + + ret <- run_stabsel(fitter = fit_model, args.fitter = list(), + n = n, p = p, cutoff = cutoff, q = q, + PFER = PFER, folds = folds, B = B, assumption = assumption, + sampling.type = sampling.type, papply = papply, + verbose = verbose, FWER = FWER, eval = eval, + names = unlist(nms), ...) + + if (verbose) + cat("\n") + + if (!eval) + return(ret) + + violations <- FALSE + if (!is.null(attr(ret, "violations"))) + violations <- attr(ret, "violations") + + if (any(violations)) + warning(sQuote("mstop"), " too small in ", + sum(violations), " of the ", length(violations), + " subsampling replicates to select ", sQuote("q"), + " base-learners; Increase ", sQuote("mstop"), + " bevor applying ", sQuote("stabsel")) + + ret$call <- cll + ret$call[[1]] <- as.name("stabsel") + class(ret) <- c("stabsel_mboostLSS", "stabsel") + ret +} + ################################################################################ ### helpers @@ -335,11 +600,11 @@ weighted.sd <- function(x, w, ...) { weighted.median <- function (x, w = 1, na.rm = FALSE) { if (length(w) == 1) w <- rep(w, length(x)) - + ## remove observations with zero weights x <- x[w != 0] w <- w[w != 0] - + ## remove NAs if na.rm = TRUE if (na.rm) { keep <- !is.na(x) & !is.na(w) @@ -349,22 +614,22 @@ weighted.median <- function (x, w = 1, na.rm = FALSE) { if (any(is.na(x)) | any(is.na(w))) return(NA) } - + ## sort data and weights ind <- order(x) x <- x[ind] w <- w[ind] - + ## first time that fraction of weights is above 0.5 ind1 <- min(which(cumsum(w)/sum(w) > 0.5)) - + ## first time that fraction of weights is below 0.5 ind2 <- ifelse(ind1 == 1, 1, max(which(cumsum(w)/sum(w) <= 0.5))) - + ## if sum of weights is an even integer if(sum(w) %% 1 == 0 && sum(w) %% 2 == 0) return(mean(c(x[ind1], x[ind2]))) - + ## else return return(max(c(x[ind1], x[ind2]))) } diff --git a/README.md b/README.md index 8d9710c..367a9bb 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,10 @@ gamboostLSS =========== -[![Build Status (Linux)](https://travis-ci.org/boost-R/gamboostLSS.svg?branch=master)](https://travis-ci.org/boost-R/gamboostLSS) -[![Build status (Windows)](https://ci.appveyor.com/api/projects/status/373t0tvx5v1i5ooq/branch/master?svg=true)](https://ci.appveyor.com/project/hofnerb/gamboostlss-s2whe/branch/master) +[![Build Status (Linux)](https://travis-ci.org/boost-R/gamboostLSS.svg?branch=devel)](https://travis-ci.org/boost-R/gamboostLSS) +[![Build status (Windows)](https://ci.appveyor.com/api/projects/status/373t0tvx5v1i5ooq/branch/devel?svg=true)](https://ci.appveyor.com/project/hofnerb/gamboostlss-s2whe/branch/devel) [![CRAN Status Badge](http://www.r-pkg.org/badges/version/gamboostLSS)](http://cran.r-project.org/package=gamboostLSS) -[![Coverage Status](https://coveralls.io/repos/github/boost-R/gamboostLSS/badge.svg?branch=master)](https://coveralls.io/github/boost-R/gamboostLSS?branch=master) +[![Coverage Status](https://coveralls.io/repos/github/boost-R/gamboostLSS/badge.svg?branch=devel)](https://coveralls.io/github/boost-R/gamboostLSS?branch=devel) [![](http://cranlogs.r-pkg.org/badges/gamboostLSS)](http://cran.rstudio.com/web/packages/gamboostLSS/index.html) `gamboostLSS` implements boosting algorithms for fitting generalized linear, @@ -22,6 +22,9 @@ For installation instructions see below. Instructions on how to use `gamboostLSS` can be found here: - [gamboostLSS tutorial](https://www.jstatsoft.org/article/view/v074i01). +Details on the noncyclical fitting method can be found here: +- [noncyclical fitting](https://arxiv.org/abs/1611.10171); This is a preliminary version currently under review. + ## Issues & Feature Requests For issues, bugs, feature requests etc. please use the [GitHub Issues](https://github.com/boost-R/gamboostLSS/issues). diff --git a/inst/CITATION b/inst/CITATION index 846b40c..0c4012e 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -68,4 +68,24 @@ citEntry( ) +citEntry( + entry= "TechReport", + title = "Stability selection for component-wise gradient boosting in multiple dimensions", + author = personList(as.person("Janek Thomas"), + as.person("Andreas Mayr"), + as.person("Bernd Bischl"), + as.person("Matthias Schmid"), + as.person("Adam Smith"), + as.person("Benjamin Hofner")), + year = "2016", + institution = "ArXiv", + header = "To cite the noncyclical fitting method of 'gamboostLSS' use:", + url = "https://arxiv.org/abs/1611.10171", + textVersion = + paste("Thomas, J., Mayr, A., Bischl, B., Schmid, M., Smith, A., and Hofner, B. (2016).", + "Stability selection for component-wise gradient boosting in multiple dimensions.", + "arXiv preprint arXiv:1611.10171.") +) + + citFooter('\nUse ', sQuote('toBibtex(citation("gamboostLSS"))'), ' to extract BibTeX references.') diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 6627fa9..a2619cb 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -1,16 +1,41 @@ \name{NEWS} \title{News for Package 'gamboostLSS'} - -\section{Changes in gamboostLSS version 1.2-2 (2016-10-18)}{ +\section{Changes in gamboostLSS version 1.5-0 (2016-yy-zz)}{ + \subsection{User-visible changes}{ + \itemize{ + \item Added new non-cyclic fitting \code{method}s for \code{gamboostLSS} models: + \itemize{ + \item \code{method = "cycling"} is the standard approach which was also + used previously. + \item \code{method = "inner"} and \code{"outer"} provide two new + non-cyclic algorithms. For details see XXX. + \item Non-cyclic methods allow faster cross-validation and better results + for stability selection (see below). + } + \item Stability selection (\code{stabsel}) implemented for \code{gamboostLSS} + models. + \item The (\code{as.families}) interface to incorporate 'gamlss' distributions for \code{gamboostLSS} + models was adapted to allow for other link functions. + \item Added option \code{stabilization = "L2"} to use the mean L2 norm of + the negative gradient to make the updated for the different distribution parameters comparable. + } + } + \subsection{Bug-fixes}{ + \itemize{ + \item Fixed \code{as.families("BB")} and \code{as.families("BI")}, which is actually an mboost family + Closes issue + \href{https://github.com/boost-R/gamboostLSS/issues/12}{#12}. + } + } \subsection{Miscellaneous}{ \itemize{ + \item Added Janek Thomas as new author. \item Updated \file{inst/CITATION} due to release of JSS tutorial paper. \item Updated URL in \file{DESCRIPTION}. } } - } - +} \section{Changes in gamboostLSS version 1.2-1 (2016-03-11)}{ \subsection{User-visible changes}{ diff --git a/man/as.families.Rd b/man/as.families.Rd index e713c12..1ffa83e 100644 --- a/man/as.families.Rd +++ b/man/as.families.Rd @@ -32,8 +32,10 @@ \code{\link{glmboostLSS}}. } \usage{ -as.families(fname = "NO", mu = NULL, sigma = NULL, - nu = NULL, tau = NULL, stabilization = c("none", "MAD")) +as.families(fname = "NO", stabilization = c("none", "MAD", "L2"), + mu = NULL, sigma = NULL, nu = NULL, tau = NULL, + mu.link = NULL, sigma.link = NULL, nu.link = NULL, + tau.link = NULL) ## a wrapper to as.families: gamlss.Families(...) @@ -51,6 +53,10 @@ gamlss.Families(...) \item{sigma}{possible offset value for parameter \code{sigma}.} \item{nu}{possible offset value for parameter \code{nu}.} \item{tau}{possible offset value for parameter \code{tau}.} + \item{mu.link}{different link function for parameter \code{mu}.} + \item{sigma.link}{different link function for parameter \code{sigma}.} + \item{nu.link}{different link function for parameter \code{nu}.} + \item{tau.link}{different link function for parameter \code{tau}.} \item{stabilization}{ governs if the negative gradient should be standardized in each boosting step. It can be either "none" or "MAD". For details see \code{\link{Families}}. @@ -67,6 +73,9 @@ gamlss.Families(...) corresponding \code{mboost}-like sub-families and the final \code{families} object, which can be then used with the fitting functions \code{\link{gamboostLSS}} and \code{\link{glmboostLSS}}. + + If no different link functions are specified, the standard links for the + corresponding family in \code{gamlss.dist} are applied. To extract the necessary information regarding partial derivatives (for the \code{ngradient} - see \code{\link{Family}} for details) and diff --git a/man/cvrisk.Rd b/man/cvrisk.Rd index 3fd58fb..23db946 100644 --- a/man/cvrisk.Rd +++ b/man/cvrisk.Rd @@ -1,8 +1,10 @@ \name{cvrisk.mboostLSS} \alias{cvrisk} \alias{cvrisk.mboostLSS} +\alias{cvrisk.nc_mboostLSS} \alias{make.grid} \alias{plot.cvriskLSS} +\alias{plot.nc_cvriskLSS} \title{ Cross-Validation } \description{ @@ -20,6 +22,10 @@ make.grid(max, length.out = 10, min = NULL, log = TRUE, \method{plot}{cvriskLSS}(x, type = c("heatmap", "lines"), xlab = NULL, ylab = NULL, ylim = range(x), main = attr(x, "type"), ...) + +\method{plot}{nc_cvriskLSS}(x, type = "lines", + xlab = "Number of boosting iterations", ylab = NULL, + ylim = range(x), main = attr(x, "type"), ...) } \arguments{ @@ -34,10 +40,15 @@ make.grid(max, length.out = 10, min = NULL, log = TRUE, bootstrap samples. } \item{grid}{ - a matrix of stopping parameters the empirical risk is to be - evaluated for. Each row represents a parameter combination. The - number of columns must be equal to the number of parameters of the - GAMLSS family. Per default, \code{make.grid(mstop(object))} is used. + If the model was fitted with \code{method = "cyclic"}, grid is + a matrix of stopping parameters the empirical risk is to be evaluated for. + Each row represents a parameter combination. The number of columns must be + equal to the number of parameters of the GAMLSS family. Per default, + \code{make.grid(mstop(object))} is used. + + Otherwise (i.e., for \code{method = "noncyclic"}) grid + is a vector of mstop values. Per default all steps up to the intial stopping + iteration, i.e., \code{1:mstop(object)} are used. } \item{papply}{ (parallel) apply function, defaults to \code{\link[parallel]{mclapply}}. @@ -120,6 +131,8 @@ make.grid(max, length.out = 10, min = NULL, log = TRUE, additionally shows information on the variability of the risk from fold to fold. The heatmap shows only the average risk but in a nicer fashion. + + For the \code{method = "noncyclic"} only the line plot for exists. Hofner et al. (2016) provide a detailed description of cross-validation for \code{\link{gamboostLSS}} models and show a diff --git a/man/families.Rd b/man/families.Rd index de4b72a..465108d 100644 --- a/man/families.Rd +++ b/man/families.Rd @@ -62,40 +62,40 @@ # Gaussian distribution GaussianLSS(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) # Student's t-distribution StudentTLSS(mu = NULL, sigma = NULL, df = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) ############################################################ # Families for continuous non-negative response # Gamma distribution GammaLSS(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) ############################################################ # Families for fractions and bounded continuous response # Beta distribution BetaLSS(mu = NULL, phi = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) ############################################################ # Families for count data # Negative binomial distribution NBinomialLSS(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) # Zero-inflated Poisson distribution ZIPoLSS(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) # Zero-inflated negative binomial distribution ZINBLSS(mu = NULL, sigma = NULL, nu = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) ############################################################ # Families for survival models (accelerated failure time @@ -103,15 +103,15 @@ ZINBLSS(mu = NULL, sigma = NULL, nu = NULL, # Log-normal distribution LogNormalLSS(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) # Log-logistic distribution LogLogLSS(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) # Weibull distribution WeibullLSS(mu = NULL, sigma = NULL, - stabilization = c("none", "MAD")) + stabilization = c("none", "MAD", "L2")) ############################################################ # Constructor function for new GAMLSS distributions @@ -130,8 +130,8 @@ Families(..., qfun = NULL, name = NULL) \item{df}{ offset value for df. } \item{nu}{ offset value for nu. } \item{stabilization}{ governs if the negative gradient should be - standardized in each boosting step. It can be either "none" or - "MAD". See also Details below. + standardized in each boosting step. It can be either "none", + "MAD" or "L2". See also Details below. } } \details{ @@ -183,8 +183,8 @@ Families(..., qfun = NULL, name = NULL) negative gradients one can use the argument \code{stabilization} of the families. If \code{stabilization = "MAD"}, the negative gradient is divided by its (weighted) median absolute deviation \deqn{median_i - (|u_{k,i} - median_j(u_{k,j})|)} in each boosting step. See Hofner et - al. (2016) for details. + (|u_{k,i} - median_j(u_{k,j})|)} in each boosting step. See Hofner et + al. (2016) for details. An alternative is \code{stabilization = "L2"}, where the gradient is divided by its (weighted) mean L2 norm. This results in negative gradient vectors (and hence also updates) of similar size for each distribution parameter, but also for every boosting iteration. } \value{ @@ -246,7 +246,7 @@ newStudentTMu <- function(mu, df){ risk <- function(y, f, w = 1) { sum(w * loss(y = y, f = f, df = df)) } - # ngradient is the negative derivate w.r.t. mu + # ngradient is the negative derivate w.r.t. mu (=f) ngradient <- function(y, f, w = 1) { (df + 1) * (y - f)/(df + (y - f)^2) } @@ -271,7 +271,8 @@ newStudentTDf <- function(mu, df){ risk <- function(y, f, w = 1) { sum(w * loss(y = y, f = f, mu = mu)) } - # ngradient is the negative derivate w.r.t. df + # ngradient is the negative derivate of the loss w.r.t. f + # in this case, just the derivative of the log-likelihood ngradient <- function(y, f, w = 1) { exp(f)/2 * (digamma((exp(f) + 1)/2) - digamma(exp(f)/2)) - 0.5 - (exp(f)/2 * log(1 + (y - mu)^2 / (exp(f) )) - diff --git a/man/gamboostLSS-package.Rd b/man/gamboostLSS-package.Rd index 977c49f..02a437f 100644 --- a/man/gamboostLSS-package.Rd +++ b/man/gamboostLSS-package.Rd @@ -43,7 +43,7 @@ the covariate has on each distribution parameter in the final GAMLSS. } \author{ - Benjamin Hofner, Andreas Mayr, Nora Fenske, Matthias Schmid + Benjamin Hofner, Andreas Mayr, Nora Fenske, Janek Thomas, Matthias Schmid Maintainer: Benjamin Hofner } diff --git a/man/mboostLSS.Rd b/man/mboostLSS.Rd index 3024d41..348e0a3 100644 --- a/man/mboostLSS.Rd +++ b/man/mboostLSS.Rd @@ -17,18 +17,22 @@ } \usage{ mboostLSS(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...) + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...) glmboostLSS(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...) + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...) gamboostLSS(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...) + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...) blackboostLSS(formula, data = list(), families = GaussianLSS(), - control = boost_control(), weights = NULL, ...) + control = boost_control(), weights = NULL, + method = c("cyclic", "noncyclic"), ...) ## fit function: mboostLSS_fit(formula, data = list(), families = GaussianLSS(), control = boost_control(), weights = NULL, - fun = mboost, funchar = "mboost", call = NULL, ...) + fun = mboost, funchar = "mboost", call = NULL, method, ...) } \arguments{ @@ -57,6 +61,9 @@ mboostLSS_fit(formula, data = list(), families = GaussianLSS(), \item{call}{ used to forward the call from \code{mboostLSS}, \code{glmboostLSS}, \code{gamboostLSS} and \code{blackboostLSS}. This argument should not be directly specified by users!} + \item{method}{ fitting method, currently two methods are supported: + \code{"cyclic"} and \code{"noncyclic"}. The latter two requires a one dimensional \code{mstop} + value.} \item{\dots}{Further arguments to be passed to \code{mboostLSS_fit}. In \code{mboostLSS_fit}, \code{\dots} represent further arguments to be passed to \code{\link{mboost}} and \code{\link{mboost_fit}}. So @@ -95,7 +102,7 @@ mboostLSS_fit(formula, data = list(), families = GaussianLSS(), function without \code{LSS}. For further possible arguments see these functions as well as \code{\link{mboost_fit}}. - In all four fitting functions it is possible to specify one or + For \code{method = "cyclic"} it is possible to specify one or multiple \code{mstop} and \code{nu} values via \code{\link{boost_control}}. In the case of one single value, this value is used for all distribution parameters of the GAMLSS model. @@ -110,6 +117,12 @@ mboostLSS_fit(formula, data = list(), families = GaussianLSS(), more-dimensional stopping, one can specify, e.g., \code{mstop = list(mu = 100, sigma = 200)} (see examples). + If \code{method} is set to \code{"noncyclic"}, \code{mstop} has + to be a one dimensional integer. Instead of cycling through all distribution + parameters, in each iteration only the best baselearner is used. One baselearner of every + parameter is selected via RSS, the distribution parameter is then chosen via the loss. + For details on the noncyclic fitting method see Thomas et. al. (2016). + To (potentially) stabilize the model estimation by standardizing the negative gradients one can use the argument \code{stabilization} of the families. See \code{\link{Families}} for details. @@ -142,6 +155,10 @@ Statistical Society, Series C (Applied Statistics), 54, 507-554. Buehlmann, P. and Hothorn, T. (2007), Boosting algorithms: Regularization, prediction and model fitting. Statistical Science, 22(4), 477--505. + +Thomas, J., Mayr, A., Bischl, B., Schmid, M., Smith, A., and Hofner, B. (2016), +Stability selection for component-wise gradient boosting in multiple dimensions. +arXiv preprint arXiv:1611.10171. } \seealso{ diff --git a/man/methods.Rd b/man/methods.Rd index 0b9a0bd..7c63c02 100644 --- a/man/methods.Rd +++ b/man/methods.Rd @@ -79,7 +79,7 @@ PI(x, which, pi = 0.9, newdata = NULL, ...) \method{risk}{mboostLSS}(object, merge = FALSE, parameter = names(object), ...) ### extract selected base-learners -\method{selected}{mboostLSS}(object, parameter = names(object), ...) +\method{selected}{mboostLSS}(object, merge = FALSE, parameter = names(object), ...) ### extract fitted values \method{fitted}{mboostLSS}(object, parameter = names(object), ...) diff --git a/man/stabsel.mboostLSS.Rd b/man/stabsel.mboostLSS.Rd new file mode 100644 index 0000000..84800a0 --- /dev/null +++ b/man/stabsel.mboostLSS.Rd @@ -0,0 +1,145 @@ +\name{stabsel} +\alias{stabsel.mboostLSS} +\alias{selected.stabsel_mboostLSS} + +\title{ + Stability Selection +} +\description{ + Selection of influential variables or model components with error control. +} +\usage{ +## a method to compute stability selection paths for fitted mboostLSS models +\method{stabsel}{mboostLSS}(x, cutoff, q, PFER, mstop = NULL, + folds = subsample(model.weights(x), B = B), + B = ifelse(sampling.type == "MB", 100, 50), + assumption = c("unimodal", "r-concave", "none"), + sampling.type = c("SS", "MB"), + papply = mclapply, verbose = TRUE, FWER, eval = TRUE, ...) +## a method to get the selected parameters +\method{selected}{stabsel_mboostLSS}(object, parameter = NULL, ...) +} +\arguments{ + \item{x}{an fitted model of class \code{"mboostLSS"} or \code{"nc_mboostLSS"}.} + \item{cutoff}{cutoff between 0.5 and 1. Preferably a value between 0.6 + and 0.9 should be used.} + \item{q}{number of (unique) selected variables (or groups of variables + depending on the model) that are selected on each subsample.} + \item{PFER}{upper bound for the per-family error rate. This + specifies the amount of falsely selected base-learners, which is + tolerated. See details.} + \item{mstop}{mstop value to use, if no value is supplied the mstop value of the fitted model is used.} + \item{folds}{ a weight matrix with number of rows equal to the number + of observations, see \code{\link{cvrisk}} and + \code{\link[stabs]{subsample}}. Usually one should not + change the default here as subsampling with a fraction of \eqn{1/2} + is needed for the error bounds to hold. One usage scenario where + specifying the folds by hand might be the case when one has + dependent data (e.g. clusters) and thus wants to draw clusters + (i.e., multiple rows together) not individuals.} + \item{assumption}{ Defines the type of assumptions on the + distributions of the selection probabilities and simultaneous + selection probabilities. Only applicable for + \code{sampling.type = "SS"}. For \code{sampling.type = "MB"} we + always use code{"none"}.} + \item{sampling.type}{ use sampling scheme of of Shah & Samworth + (2013), i.e., with complementarty pairs (\code{sampling.type = "SS"}), + or the original sampling scheme of Meinshausen & Buehlmann (2010).} + \item{B}{ number of subsampling replicates. Per default, we use 50 + complementary pairs for the error bounds of Shah & Samworth (2013) + and 100 for the error bound derived in Meinshausen & Buehlmann + (2010). As we use \eqn{B} complementray pairs in the former case + this leads to \eqn{2B} subsamples.} + \item{papply}{ (parallel) apply function, defaults to + \code{\link[parallel]{mclapply}}. Alternatively, \code{parLapply} + can be used. In the latter case, usually more setup is needed (see + example of \code{\link{cvrisk}} for some details).} + \item{verbose}{ logical (default: \code{TRUE}) that determines wether + \code{warnings} should be issued. } + \item{FWER}{ deprecated. Only for compatibility with older versions, + use PFER instead.} + \item{eval}{ logical. Determines whether stability selection is + evaluated (\code{eval = TRUE}; default) or if only the parameter + combination is returned.} + \item{object}{ a object of class \code{"stabsel_mboostLSS"}.} + \item{parameter}{ select one or multiple effects.} + \item{\dots}{ additional arguments to parallel apply methods such as + \code{\link{mclapply}} and to \code{\link{cvrisk}}.} +} +\details{ + + For details see \code{\link[stabs]{stabsel}} in package \pkg{stabs} + and Hofner et al. (2014). + +} +\value{ + An object of class \code{stabsel} with a special \code{print} method. + The object has the following elements: + \item{phat}{selection probabilities.} + \item{selected}{elements with maximal selection probability greater + \code{cutoff}.} + \item{max}{maximum of selection probabilities.} + \item{cutoff}{cutoff used.} + \item{q}{average number of selected variables used.} + \item{PFER}{per-family error rate.} + \item{sampling.type}{the sampling type used for stability selection.} + \item{assumption}{the assumptions made on the selection + probabilities.} + \item{call}{the call.} +} +\references{ + + B. Hofner, L. Boccuto and M. Goeker (2014), + Controlling false discoveries in high-dimensional situations: Boosting + with stability selection. \emph{Technical Report}, arXiv:1411.1285.\cr + \url{http://arxiv.org/abs/1411.1285}. + + N. Meinshausen and P. Buehlmann (2010), Stability selection. + \emph{Journal of the Royal Statistical Society, Series B}, + \bold{72}, 417--473. + + R.D. Shah and R.J. Samworth (2013), Variable selection with error + control: another look at stability selection. \emph{Journal of the Royal + Statistical Society, Series B}, \bold{75}, 55--80. + +} +\seealso{ + \code{\link[stabs]{stabsel}} and + \code{\link[stabs]{stabsel_parameters}} +} +\examples{ + +### Data generating process: +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) + +### 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") + +\donttest{### Do not test the following code per default on CRAN as it takes some time to run: + +#run stability selection +(s <- stabsel(model, q = 5, PFER = 1)) +#get selected effects +selected(s) + +#visualize selection frequencies +plot(s) + +### END (don't test automatically) +} +} +\keyword{nonparametric} \ No newline at end of file diff --git a/tests/bugfixes.R b/tests/bugfixes.R index a6b5e09..b8eef2e 100644 --- a/tests/bugfixes.R +++ b/tests/bugfixes.R @@ -2,6 +2,7 @@ # Bugfixes require("gamboostLSS") +require("gamlss") ## subset method was missing if initial mstop = 1 set.seed(1907) @@ -28,8 +29,7 @@ stopifnot(all.equal(selected(model), list(mu = mboost::selected(model[[1]]), sigma = mboost::selected(model[[2]])))) - -# If the families argument is not specified explicitly in mboostLSSone gets an +# If the families argument is not specified explicitly in mboostLSS one gets an # error in cvrisk.mboostLSS() (spotted by Almond Stöcker). # (https://github.com/boost-R/gamboostLSS/issues/9) set.seed(1907) @@ -47,3 +47,24 @@ cvr <- try(cvrisk(model, folds = cv(model.weights(model), B=3), trace=FALSE), silent = TRUE) if (inherits(cvr, "try-error")) stop("cvrisk does not work if no family was (explicitly) chosen") + + +# Make sure that gamlss.dist::BB and gamlss.dist::BI work (spotted by F. Scheipl) +# (https://github.com/boost-R/gamboostLSS/issues/12) +set.seed(123) +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)) + +m1 <- gamlss(ymat ~ x + z, data = data, family = BB) +m2 <- gamlss(ymat ~ x + z, data = data, family = BI) +# same with boosting +m3 <- glmboostLSS(ymat ~ x + z, data = data, families = as.families("BB")) +m4 <- glmboost(ymat ~ x + z, data = data, family = as.families("BI")) + +round(data.frame(BB_gamlss = coef(m1), + BI_gamlss = coef(m2), + BB_gamboostLSS = coef(m3, off2int = TRUE, parameter = "mu"), + BI_gamboostLSS = coef(m4, off2int = TRUE)), 3) diff --git a/tests/regtest-families.R b/tests/regtest-families.R index 3cf0482..dd82a72 100644 --- a/tests/regtest-families.R +++ b/tests/regtest-families.R @@ -57,6 +57,61 @@ model <- glmboostLSS(y ~ x1 + x2, families = StudentTLSS(mu = res[1], sigma = model[100] coef(model) +### 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") +coef(model2, off2int = TRUE) # GaussianLSS() + +### 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) + + +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) + + +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) + +### 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) + + +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) + + ### check survival families x1 <- runif(1000) diff --git a/tests/regtest-noncyclic_fitting.R b/tests/regtest-noncyclic_fitting.R new file mode 100644 index 0000000..a312024 --- /dev/null +++ b/tests/regtest-noncyclic_fitting.R @@ -0,0 +1,125 @@ +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))) + stopifnot(all.equal(predict(m1), predict(m2))) + stopifnot(all.equal(fitted(m1), fitted(m2))) + stopifnot(all.equal(selected(m1), selected(m2))) + stopifnot(all.equal(risk(m1), risk(m2))) + ## remove obvious differences from objects + m1$control <- m2$control <- NULL + m1$call <- m2$call <- NULL + if (!all.equal(m1, m2)) + 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) + diff --git a/tests/regtest-stabilization.R b/tests/regtest-stabilization.R index d774d43..9056b48 100644 --- a/tests/regtest-stabilization.R +++ b/tests/regtest-stabilization.R @@ -36,6 +36,7 @@ stopifnot(all.equal(coef(m3), coef(m2))) ## check if everything is handled correctly GaussianLSS(stabilization = "MAD") +GaussianLSS(stabilization = "L2") GaussianLSS(stabilization = "none") res <- try(GaussianLSS(stabilization = "test"), silent = TRUE) res @@ -56,11 +57,17 @@ for (i in 1:length(FAMILIES)) { m_MAD <- glmboostLSS(y ~ x1 + x2 + x3 + x4, families = FAMILIES[[i]](stabilization = "MAD"), data=dat) + m_L2 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "L2"), + data=dat) + stopifnot(tail(risk(m_none, merge = TRUE), 1) != tail(risk(m_MAD, merge = TRUE), 1)) cat('Risks:\n stabilization = "none":', tail(risk(m_none, merge = TRUE), 1), '\n stabilization = "MAD":', - tail(risk(m_MAD, merge = TRUE), 1), "\n") + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") } ## check as.families interface for 2:4 parametric families @@ -75,10 +82,15 @@ for (i in 1:length(FAMILIES)) { m_MAD <- glmboostLSS(y ~ x1 + x2 + x3 + x4, families = as.families(FAMILIES[[i]], stabilization = "MAD"), data=dat) + m_L2 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "L2"), + data=dat) cat('Risks:\n stabilization = "none":', tail(risk(m_none, merge = TRUE), 1), '\n stabilization = "MAD":', - tail(risk(m_MAD, merge = TRUE), 1), "\n") + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") } FAMILIES <- list("BCT") @@ -92,13 +104,23 @@ for (i in 1:length(FAMILIES)) { families = as.families(FAMILIES[[i]], stabilization = "MAD"), data=dat), silent = TRUE) if (inherits(m_MAD, "try-error")) { - warning("BCT cannot be fitted with stabilization", immediate. = TRUE) + warning("BCT cannot be fitted with stabilization = 'MAD'", immediate. = TRUE) break } + + m_L2 <- try(glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = as.families(FAMILIES[[i]], stabilization = "L2"), + data=dat), silent = TRUE) + if (inherits(m_L2, "try-error")) { + warning("BCT cannot be fitted with stabilization = 'L2'", immediate. = TRUE) + break + } cat('Risks:\n stabilization = "none":', tail(risk(m_none, merge = TRUE), 1), '\n stabilization = "MAD":', - tail(risk(m_MAD, merge = TRUE), 1), "\n") + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") } ## check survival families @@ -119,13 +141,23 @@ for (i in 1:length(FAMILIES)) { data=dat), silent = TRUE) if (inherits(m_MAD, "try-error")) { - warning(families[i], "cannot be fitted with stabilization", immediate. = TRUE) + warning(families[i], " cannot be fitted with stabilization = 'MAD'", immediate. = TRUE) break } - cat('Risks:\n stabilization = "none":', + m_L2 <- try(glmboostLSS(Surv(y, zens) ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "L2"), + data=dat), silent = TRUE) + + if (inherits(m_L2, "try-error")) { + warning(families[i], " cannot be fitted with stabilization = 'L2'", immediate. = TRUE) + break + } + cat('Risks:\n stabilization = "none":', tail(risk(m_none, merge = TRUE), 1), '\n stabilization = "MAD":', - tail(risk(m_MAD, merge = TRUE), 1), "\n") + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") } ## check count data @@ -141,8 +173,14 @@ for (i in 1:length(FAMILIES)) { m_MAD <- glmboostLSS(y ~ x1 + x2 + x3 + x4, families = FAMILIES[[i]](stabilization = "MAD"), data=dat) + m_L2 <- glmboostLSS(y ~ x1 + x2 + x3 + x4, + families = FAMILIES[[i]](stabilization = "L2"), + data=dat) + cat('Risks:\n stabilization = "none":', tail(risk(m_none, merge = TRUE), 1), '\n stabilization = "MAD":', - tail(risk(m_MAD, merge = TRUE), 1), "\n") + tail(risk(m_MAD, merge = TRUE), 1), + '\n stabilization = "L2":', + tail(risk(m_L2, merge = TRUE), 1), "\n") }