[Mboost-commits] r777 - in pkg/mboostDevel: . R inst man tests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Jul 3 17:47:27 CEST 2014


Author: hofner
Date: 2014-07-03 17:47:27 +0200 (Thu, 03 Jul 2014)
New Revision: 777

Added:
   pkg/mboostDevel/R/confint.R
Modified:
   pkg/mboostDevel/DESCRIPTION
   pkg/mboostDevel/NAMESPACE
   pkg/mboostDevel/inst/CHANGES
   pkg/mboostDevel/man/mboost_package.Rd
   pkg/mboostDevel/tests/regtest-inference.R
Log:
- first (experimental) version of bootstrap CIs for boosting models


Modified: pkg/mboostDevel/DESCRIPTION
===================================================================
--- pkg/mboostDevel/DESCRIPTION	2014-07-03 15:45:02 UTC (rev 776)
+++ pkg/mboostDevel/DESCRIPTION	2014-07-03 15:47:27 UTC (rev 777)
@@ -1,7 +1,7 @@
 Package: mboostDevel
 Title: Model-Based Boosting
-Version: 2.3-0
-Date: 2014-06-26
+Version: 2.4-0
+Date: 2014-xx-yy
 Authors at R: c(person("Torsten", "Hothorn", role = c("aut", "cre"),
                     email = "Torsten.Hothorn at R-project.org"),
              person("Peter", "Buehlmann", role = "aut"),

Modified: pkg/mboostDevel/NAMESPACE
===================================================================
--- pkg/mboostDevel/NAMESPACE	2014-07-03 15:45:02 UTC (rev 776)
+++ pkg/mboostDevel/NAMESPACE	2014-07-03 15:47:27 UTC (rev 777)
@@ -82,5 +82,10 @@
 S3method(extract, bl_tree)
 S3method(residuals, mboost)
 S3method(risk, mboost)
+S3method(confint, mboost)
+S3method(confint, glmboost)
+S3method(plot, mboost.ci)
+S3method(lines, mboost.ci)
+S3method(print, glmboost.ci)
 
 useDynLib(mboostDevel)

Added: pkg/mboostDevel/R/confint.R
===================================================================
--- pkg/mboostDevel/R/confint.R	                        (rev 0)
+++ pkg/mboostDevel/R/confint.R	2014-07-03 15:47:27 UTC (rev 777)
@@ -0,0 +1,183 @@
+
+confint.mboost <- function(object, B = 1000, newdata = NULL,
+                           B.mstop = 25, which = NULL, ...) {
+
+    which <- object$which(which, usedonly = FALSE)
+
+    ## create new data and/or restructure data
+    newdata <- .create_newdata(object, newdata, which)
+
+    outer.folds <- cv(model.weights(object), B = B)
+    predictions <- vector("list", B)
+
+    for (i in 1:B) {
+        ## update model
+        mod <- update(object, weights = outer.folds[, i],
+                      risk = "inbag")
+        if (B.mstop > 0) {
+            ## <FIXME> are the weights handled correctly?
+            cvr <- cvrisk(mod, folds = cv(model.weights(mod), B = B.mstop))
+            mod[mstop(cvr)]
+        }
+        predictions[[i]] <- .predict_confint(mod, newdata = newdata,
+                                             which = which)
+    }
+    res <- list(boot_pred = predictions, data = newdata, model = object)
+    class(res) <- "mboost.ci"
+    return(res)
+}
+
+confint.glmboost <- function(object, B = 1000, B.mstop = 25,
+                             which = NULL, ...) {
+
+    outer.folds <- cv(model.weights(object), B = B)
+    which <- object$which(which, usedonly = FALSE)
+
+    coefficients <- matrix(NA, ncol = length(which), nrow = B)
+    colnames(coefficients) <- names(coef(object, which = which))
+
+    for (i in 1:B) {
+        ## update model
+        mod <- update(object, weights = outer.folds[, i],
+                      risk = "inbag")
+        if (B.mstop > 0) {
+            ## <FIXME> are the weights handled correctly?
+            cvr <- cvrisk(mod, folds = cv(model.weights(mod), B = B.mstop))
+            mod[mstop(cvr)]
+        }
+        coefficients[i, ] <- unlist(coef(mod, which = which, off2int = TRUE))
+    }
+    res <- list(boot_coefs = coefficients, model = object)
+    class(res) <- "glmboost.ci"
+    return(res)
+}
+
+print.glmboost.ci <- function(x, which = NULL, level = 0.95) {
+    quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
+    which <- x$model$which(which, usedonly = FALSE)
+    tmp <- apply(x$boot_coefs[, which], 2, FUN = quantile, probs = quantiles)
+    print(t(tmp))
+}
+
+## ## check for varing...
+## data <- model.frame(x, which = w)[[1]]
+## get_vary <- x$baselearner[[w]]$get_vary
+## vary <- ""
+## if (!is.null(get_vary)) vary <- get_vary()
+## if (!is.null(newdata)) data <- newdata[, colnames(data), drop = FALSE]
+## if (vary != "") {
+##     v <- data[[vary]]
+##     if (is.factor(v)) v <- factor(levels(v)[-1], levels = levels(v))
+##     if (is.numeric(v)) v <- 1
+## }
+
+## ## Aditionally needed: Check for multivariate base-learners (except bols)
+
+
+## check for by variable and bivariate base-learners which both need a different
+## data set for prediction
+.create_newdata <- function(object, newdata = NULL, which, ...) {
+    if (is.null(newdata)) {
+        data <- newdata <- model.frame(object, which = which)
+        for (w in which) {
+            ## make grid!
+            tmp <- data[[w]][rep(1, 100), , drop = FALSE]
+            grid <- function(x) {
+                if (is.numeric(x)) {
+                    return(seq(min(x), max(x), length = 100))
+                } else {
+                    return(rep(levels(x), length.out = 100))
+                }
+            }
+            for (j in 1:ncol(data[[w]]))
+                tmp[, colnames(data[[w]])[j]] <- grid(data[[w]][,j])
+            rownames(tmp) <- NULL
+            newdata[[w]] <- tmp
+        }
+    } else {
+        ## restructure new data
+        data <- model.frame(object, which = which)
+        nms <- lapply(data, colnames)
+        tmp <- lapply(nms, function(x) newdata[, x, drop = FALSE])
+        newdata <- tmp
+    }
+    return(newdata)
+}
+
+## .create_newdata.glmboost <- function(object, newdata, ...) {
+##     if (is.null(newdata)) {
+##         data <- model.frame(object)
+##         ## make grid!
+##         tmp <- data[rep(1, 100), ]
+##         grid <- function(x) {
+##             if (is.numeric(x)) {
+##                 return(seq(min(x), max(x), length = 100))
+##             } else {
+##                 return(rep(levels(x), length.out = 100))
+##             }
+##         }
+##         for (j in 1:ncol(data))
+##             tmp[, colnames(data)[j]] <- grid(data[,j])
+##         newdata <- tmp
+##     }
+##     return(newdata)
+## }
+
+## special prediction function for the construction of confidence intervals:
+.predict_confint <- function(object, newdata = NULL, which, ...) {
+    predictions <- matrix(NA, ncol = length(which), nrow = nrow(newdata[[1]]))
+    for (w in which) {
+        predictions[, w] <- predict(object, newdata[[w]], which = w)
+    }
+    return(predictions)
+}
+
+# .predict_confint.glmboost <- function(object, newdata, which, ...) {
+#     warning("shouldn't we return confints for coef?")
+#     predict(object, newdata = newdata, which = which)
+# }
+
+plot.mboost.ci <- function(x, which, level = 0.95,
+                           ylim = NULL, type = "l", col = "black",
+                           ci.col = "grey",  raw = FALSE, ...) {
+
+    which <- x$model$which(which, usedonly = FALSE)
+
+    if (is.null(ylim)) {
+        preds <- sapply(x$boot_pred, function(p) p[, which])
+        if (!raw) {
+            quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
+            tmp <- apply(preds, 1, FUN = quantile, probs = quantiles)
+            if (is.null(ylim))
+                ylim <- range(tmp)
+        } else {
+            if (is.null(ylim))
+                ylim <- range(preds)
+        }
+    }
+
+    plot(x$model, which = which, type = "n", ylim = ylim,
+         col = col, ...)
+    lines(x, which, level, col = ci.col, raw = raw, ...)
+    lines(x$model, which = which, type = "l", col = col, ...)
+}
+
+lines.mboost.ci <- function(x, which, level = 0.95, col = "grey",
+                            raw = FALSE, ...) {
+    preds <- sapply(x$boot_pred, function(p) p[, which])
+    x.data <- x$data[[which]]
+    if (ncol(x.data) > 1) {
+        stop("Cannot plot lines for more than 1 dimenstion")
+    } else {
+        x.data <- x.data[, 1]
+    }
+    if (!raw) {
+        quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
+        tmp <- apply(preds, 1, FUN = quantile, probs = quantiles)
+        polygon(c(x.data, rev(x.data)),
+                c(tmp[1, ], rev(tmp[2,])),
+                col = col, border = col)
+    } else {
+        matlines(x$data[[which]], preds, col = col, lty = "solid", ...)
+    }
+}

Modified: pkg/mboostDevel/inst/CHANGES
===================================================================
--- pkg/mboostDevel/inst/CHANGES	2014-07-03 15:45:02 UTC (rev 776)
+++ pkg/mboostDevel/inst/CHANGES	2014-07-03 15:47:27 UTC (rev 777)
@@ -1,4 +1,15 @@
+                CHANGES in `mboost' VERSION 2.4-0 (2014-xx-yy, rZZZ)
 
+  o  added functions to compute (bootstrap) confidence intervals
+
+                CHANGES in `mboost' VERSION 2.3-1 (2014-xx-yy, rXYZ)
+
+  o  changed vignette mboost_tutorial to reflect latest changes in mboost.
+
+  o  Bugfixes: 
+     - glmboost()$model.frame() was broken
+     - glmboost()$update() was broken
+
                 CHANGES in `mboost' VERSION 2.3-0 (2014-06-26, r771)
 
   o  stabsel was recoded and now uses different terminology, much more options
@@ -22,6 +33,8 @@
   o  boost_control: added new argument stopintern for internal stopping
      (based on oobag data) during fitting
 
+  o  All data sets have been moved to the new package set TH.data
+
   o  Misc:
      - added new argmument which to variable.names()
      - added new method risk to extract risks

Modified: pkg/mboostDevel/man/mboost_package.Rd
===================================================================
--- pkg/mboostDevel/man/mboost_package.Rd	2014-07-03 15:45:02 UTC (rev 776)
+++ pkg/mboostDevel/man/mboost_package.Rd	2014-07-03 15:47:27 UTC (rev 777)
@@ -15,8 +15,8 @@
 \tabular{ll}{
 Package: \tab mboostDevel\cr
 Type: \tab Package\cr
-Version: \tab 2.3-0\cr
-Date: \tab 2014-06-26\cr
+Version: \tab 2.4-0\cr
+Date: \tab 2014-xx-yy\cr
 License: \tab GPL-2\cr
 LazyLoad: \tab yes\cr
 LazyData: \tab yes\cr

Modified: pkg/mboostDevel/tests/regtest-inference.R
===================================================================
--- pkg/mboostDevel/tests/regtest-inference.R	2014-07-03 15:45:02 UTC (rev 776)
+++ pkg/mboostDevel/tests/regtest-inference.R	2014-07-03 15:47:27 UTC (rev 777)
@@ -220,3 +220,21 @@
                    sampling.type = "SS", assumption = "r-concave")
 stabsel_parameters(p = p, cutoff = res$cutoff - 0.01, q = q, B = B,
                    sampling.type = "SS", assumption = "r-concave")
+
+
+### check confidence intervals
+data("bodyfat", package = "TH.data")
+bodyfat$ID <- factor(sample(1:5, size = nrow(bodyfat), replace = TRUE))
+glm <- glmboost(DEXfat ~ ., data = bodyfat)
+gam <- gamboost(DEXfat ~ ., data = bodyfat)
+
+refit <- glm$update(weights = model.weights(glm), risk = "inbag")
+stopifnot(all.equal(coef(refit), coef(glm)))
+
+confint.glm <- confint(glm, B = 100, B.mstop = 2)
+confint.glm
+
+confint.gam <- confint(gam, B = 100, B.mstop = 2)
+plot(confint.gam, which = 1)
+plot(confint.gam, which = 2)
+plot(confint.gam, which = 3)



More information about the Mboost-commits mailing list