[Mboost-commits] r787 - in pkg/mboostDevel: R inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Jul 28 12:03:01 CEST 2014


Author: hofner
Date: 2014-07-28 12:02:59 +0200 (Mon, 28 Jul 2014)
New Revision: 787

Modified:
   pkg/mboostDevel/R/confint.R
   pkg/mboostDevel/R/mboost.R
   pkg/mboostDevel/R/plot.R
   pkg/mboostDevel/inst/CHANGES
   pkg/mboostDevel/man/confint.Rd
Log:
- improved plot method for varying coefficients (ylim now suitable) and
  base-learners of factor variables.
- confint now works with interaction effects (by) and spatial base-learners
- Fixed newly introduced bug in update()


Modified: pkg/mboostDevel/R/confint.R
===================================================================
--- pkg/mboostDevel/R/confint.R	2014-07-25 14:42:58 UTC (rev 786)
+++ pkg/mboostDevel/R/confint.R	2014-07-28 10:02:59 UTC (rev 787)
@@ -1,7 +1,10 @@
 
 confint.mboost <- function(object, parm = NULL, level = 0.95,
                            B = 1000, B.mstop = 25, newdata = NULL,
-                           which = parm, ...) {
+                           which = parm,
+                           papply = ifelse(B.mstop == 0, mclapply, lapply),
+                           papply.mstop = mclapply,
+                           ...) {
 
     which <- object$which(which, usedonly = FALSE)
     if (!all(which %in% object$which(NULL, usedonly = FALSE)))
@@ -11,22 +14,24 @@
     newdata <- .create_newdata(object, newdata, which)
 
     outer.folds <- cv(model.weights(object), B = B)
-    predictions <- vector("list", B)
 
     cat("Start computing bootstrap confidence intervals... \n")
-    for (i in 1:B) {
+
+    do_update <- function(i) {
+        #for (i in 1:B) {
         cat("\rB =", i)
         ## update model
         mod <- update(object, weights = outer.folds[, i],
                       risk = "inbag", trace = FALSE)
         if (B.mstop > 0) {
             ## <FIXME> are the weights handled correctly?
-            cvr <- cvrisk(mod, folds = cv(model.weights(mod), B = B.mstop))
+            cvr <- cvrisk(mod, folds = cv(model.weights(mod), B = B.mstop),
+                          papply = papply.mstop, ...)
             mod[mstop(cvr)]
         }
-        predictions[[i]] <- .predict_confint(mod, newdata = newdata,
-                                             which = which)
+        .predict_confint(mod, newdata = newdata, which = which)
     }
+    predictions <- papply(1:B, do_update)
     cat("\n")
 
     ## prepare returned object
@@ -55,7 +60,7 @@
                       risk = "inbag")
         if (B.mstop > 0) {
             ## <FIXME> are the weights handled correctly?
-            cvr <- cvrisk(mod, folds = cv(model.weights(mod), B = B.mstop))
+            cvr <- cvrisk(mod, folds = cv(model.weights(mod), B = B.mstop), ...)
             mod[mstop(cvr)]
         }
         coefficients[i, ] <- unlist(coef(mod, which = which, off2int = TRUE))
@@ -108,11 +113,7 @@
 }
 
 
-## ## Aditionally needed: Check for multivariate base-learners (except bols)
-
-## FIXME: check for by variable and bivariate base-learners which both need a
-## different data set for prediction
-## FIXME: what about factor variables? do we get the correct levels?
+## FIXME: Aditionally needed: Does multivariate bols base-learners work correctly?
 .create_newdata <- function(object, newdata = NULL, which, ...) {
     if (is.null(newdata)) {
         data <- newdata <- model.frame(object, which = which)
@@ -193,11 +194,6 @@
     if (is.matrix(predictions))
         predictions <- as.data.frame(predictions)
     names(predictions) <- names(newdata[which])
-
-    ## ###### FIXME FIXME FIXME
-    ## WAS ist mit predict == 0?
-    ## Wie geht es in .ci_mboost weiter?
-    ## ###### FIXME FIXME FIXME
     return(predictions)
 }
 
@@ -207,7 +203,7 @@
                            ylim = NULL, type = "l", col = "black",
                            ci.col = rgb(170, 170, 170, alpha = 85,
                                         maxColorValue = 255),
-                           raw = FALSE, ...) {
+                           raw = FALSE, print_levelplot = TRUE, ...) {
 
     if (missing(which)) {
         which <- attr(x, "which")
@@ -221,14 +217,41 @@
 
     CI <- .ci_mboost(x$boot_pred, level = level, which = which, raw = raw)
 
-    if (is.null(ylim)) {
-        ylim <- range(CI)
+    ## check if data (without by variable, which is not varying in the plot
+    ## data) has more than one column
+    varying <- which(sapply(x$data[[which]], function(x) length(unique(x))) > 1)
+    if (ncol(x$data[[which]]) > 1 && length(varying) > 1) {
+
+        if (length(varying) > 2)
+            stop("Plots only implemented for more than 2 variables.")
+
+        p1 <- plot(x$model, which = which, newdata = x$data[[which]],
+                   main = "Mean surface", ...)
+        ## make level plots for upper and lower CI
+        fm <- as.formula(paste("pr ~ ", paste(names(varying), collapse = "*"), sep = ""))
+        pr <- CI[1, ]  ## lower CI
+        p2 <- levelplot(fm, data = x$data[[which]], main = paste(rownames(CI)[1], "CI surface"), ...)
+        pr <- CI[2, ]  ## upper CI
+        p3 <- levelplot(fm, data = x$data[[which]], main = paste(rownames(CI)[2], "CI surface"), ...)
+        if (print_levelplot) {
+                   ## position = left, bottom, right, top
+            print(p1, position=c(0, 0, 0.33, 1), more=TRUE)
+            print(p2, position=c(0.33, 0, 0.66, 1), more=TRUE)
+            print(p3, position=c(0.66, 0, 1, 1))
+            warning("The scale is not the same")
+        } else {
+            return(list(mean = p1, lowerCI = p2, upperCI = p3))
+        }
+    } else {
+        if (is.null(ylim)) {
+            ylim <- range(CI)
+        }
+        plot(x$model, which = which, newdata = x$data[[which]], rug = FALSE,
+             type = "n", ylim = ylim, col = col, ...)
+        lines(x, which, level, col = ci.col, raw = raw, ...)
+        lines(x$model, which = which, newdata = x$data[[which]], rug = FALSE,
+              type = "l", col = col, ...)
     }
-
-    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 = x$level,
@@ -250,18 +273,34 @@
     CI <- .ci_mboost(x$boot_pred, level = level, which = which, raw = raw)
 
     x.data <- x$data[[which]]
-    if (ncol(x.data) > 1) {
+    ## check if data (without by variable, which is not varying in the plot
+    ## data) has more than one column
+    if (ncol(x.data) > 1 &&
+        sum(sapply(x.data, function(x) length(unique(x))) > 1) > 1) {
         stop("Cannot plot lines for more than 1 dimension")
     } else {
         x.data <- x.data[, 1]
     }
 
-    if (!raw) {
-        polygon(c(x.data, rev(x.data)),
-                c(CI[1, ], rev(CI[2,])),
-                col = col, border = col)
+    if (is.factor(x.data)) {
+        if (raw)
+            warning("plotting raw values is currently not implemented",
+                    " for factors")
+        pData <- cbind(x.data, t(CI))
+        pData <- unique(pData)
+        for (i in 1:nrow(pData)) {
+            polygon(x = pData[i, 1] + c(-0.35, 0.35, 0.35, -0.35),
+                    y = rep(pData[i, 2:3], each = 2),
+                    col = col, border = col, ...)
+        }
     } else {
-        matlines(x$data[[which]], CI, col = col, lty = "solid", ...)
+        if (!raw) {
+            polygon(c(x.data, rev(x.data)),
+                    c(CI[1, ], rev(CI[2,])),
+                    col = col, border = col)
+        } else {
+            matlines(x$data[[which]], CI, col = col, lty = "solid", ...)
+        }
     }
 }
 

Modified: pkg/mboostDevel/R/mboost.R
===================================================================
--- pkg/mboostDevel/R/mboost.R	2014-07-25 14:42:58 UTC (rev 786)
+++ pkg/mboostDevel/R/mboost.R	2014-07-28 10:02:59 UTC (rev 787)
@@ -178,7 +178,8 @@
                            trace = NULL) {
 
         control$mstop <- mstop
-        control$risk <- risk
+        if (!is.null(risk))
+            control$risk <- risk
         if (!is.null(trace))
             control$trace <- trace
         ### use user specified offset only (since it depends on weights otherwise)

Modified: pkg/mboostDevel/R/plot.R
===================================================================
--- pkg/mboostDevel/R/plot.R	2014-07-25 14:42:58 UTC (rev 786)
+++ pkg/mboostDevel/R/plot.R	2014-07-28 10:02:59 UTC (rev 787)
@@ -10,12 +10,6 @@
 
     which <- x$which(which, usedonly = is.null(which))
 
-    pr <- predict(x, which = which, newdata = newdata)
-    if (is.null(ylim)) ylim <- range(pr, na.rm = TRUE)
-    ## <FIXME> default ylim not suitable for plotting varying coefficient
-    ##         base-learners; Users need to specify suitable values themselves
-
-    ## FIXED?
     if (is.null(xlab)){
         userspec <- FALSE
         xlab <- variable.names(x)
@@ -57,20 +51,48 @@
 
         plot_helper <- function(xl, yl){
             pr <- predict(x, newdata = data, which = w)
+            if (is.null(ylim)) ylim <- range(pr, na.rm = TRUE)
+
             if (vary != "") {
                 datavary <- data[, colnames(data) == vary, drop = FALSE]
                 data <- data[, colnames(data) != vary, drop = FALSE]
             }
-
             if (ncol(data) == 1) {
-                if (!add){
-                    plot(sort(data[[1]]), pr[order(data[[1]], na.last = NA)], type = type,
-                         xlab = xl, ylab = yl, ylim = ylim, ...)
+                if (!add) {
+                    if (is.factor(data[[1]])) {
+                        xVals <- unique(sort(data[[1]]))
+                        xValsN <- as.numeric(xVals)
+                        yVals <- unique(pr[order(data[[1]], na.last = NA)])
+                        if (length(pr) == 1 && pr == 0) {
+                            yVals <- rep(0, length(xVals))
+                        }
+                        plot(xValsN, yVals,
+                             type = "n", xaxt = "n",
+                             xlim = range(as.numeric(xVals)) + c(-0.5, 0.5),
+                             xlab = xl, ylab = yl, ylim = ylim)
+                        axis(1, at = xValsN, labels = levels(xVals))
+                        for (i in 1:length(xVals)) {
+                            lines(x = rep(xValsN[i], 2) + c(-0.35, 0.35),
+                                  y = rep(yVals[i], 2), ...)
+                        }
+                    } else {
+                        plot(sort(data[[1]]), pr[order(data[[1]], na.last = NA)], type = type,
+                             xlab = xl, ylab = yl, ylim = ylim, ...)
+                    }
                     if (rug) rug(data[[1]], col = rugcol)
                 } else {
                     if (is.factor(data[[1]])){
-                        boxplot(pr[order(data[[1]], na.last = NA)] ~ sort(data[[1]]),
-                                add = TRUE, ...)
+                        xVals <- unique(sort(data[[1]]))
+                        xValsN <- as.numeric(xVals)
+                        yVals <- unique(pr[order(data[[1]], na.last = NA)])
+                        if (length(pr) == 1 && pr == 0) {
+                            yVals <- rep(0, length(xVals))
+                        }
+                        axis(1, at = xValsN, labels = levels(xVals))
+                        for (i in 1:length(xVals)) {
+                            lines(x = rep(xValsN[i], 2) + c(-0.35, 0.35),
+                                  y = rep(yVals[i], 2), ...)
+                        }
                     } else {
                         lines(sort(data[[1]]), pr[order(data[[1]], na.last = NA)], type =
                               type, ...)

Modified: pkg/mboostDevel/inst/CHANGES
===================================================================
--- pkg/mboostDevel/inst/CHANGES	2014-07-25 14:42:58 UTC (rev 786)
+++ pkg/mboostDevel/inst/CHANGES	2014-07-28 10:02:59 UTC (rev 787)
@@ -1,12 +1,20 @@
                 CHANGES in `mboost' VERSION 2.4-0 (2014-xx-yy, rZZZ)
 
-  o  added functions to compute (bootstrap) confidence intervals
+  o  added confint function to compute (bootstrap) confidence intervals
+     together with plot and print methods
 
+  o  improved plot method for varying coefficients (ylim now suitable) and
+     base-learners of factor variables.
+
+  o  tweaked update function: we now can turn the trace off and specify
+     the type of risk as well as the oobweight to update()
+
+
                 CHANGES in `mboost' VERSION 2.3-1 (2014-xx-yy, rXYZ)
 
   o  changed vignette mboost_tutorial to reflect latest changes in mboost.
 
-  o  Bugfixes: 
+  o  Bugfixes:
      - glmboost()$model.frame() was broken
      - glmboost()$update() was broken
 

Modified: pkg/mboostDevel/man/confint.Rd
===================================================================
--- pkg/mboostDevel/man/confint.Rd	2014-07-25 14:42:58 UTC (rev 786)
+++ pkg/mboostDevel/man/confint.Rd	2014-07-28 10:02:59 UTC (rev 787)
@@ -15,10 +15,12 @@
 }
 \usage{
 \method{confint}{mboost}(object, parm = NULL, level = 0.95, B = 1000,
-        B.mstop = 25, newdata = NULL, which = parm, ...)
+        B.mstop = 25, newdata = NULL, which = parm,
+        papply = ifelse(B.mstop == 0, mclapply, lapply),
+        papply.mstop = mclapply, ...)
 \method{plot}{mboost.ci}(x, which, level = x$level, ylim = NULL, type = "l", col = "black",
      ci.col = rgb(170, 170, 170, alpha = 85, maxColorValue = 255),
-     raw = FALSE, ...)
+     raw = FALSE, print_levelplot = TRUE,...)
 \method{lines}{mboost.ci}(x, which, level = x$level,
      col = rgb(170, 170, 170, alpha = 85, maxColorValue = 255),
      raw = FALSE, ...)
@@ -47,14 +49,30 @@
     number of outer bootstrap replicates used to compute the empirical
     bootstrap confidence intervals.
   }
+  \item{B.mstop}{
+    number of inner bootstrap replicates used to determine the optimal
+    mstop on each of the \code{B} bootstrap samples.
+  }
   \item{newdata}{
     optionally, a data frame on which to compute the predictions for the
     confidence intervals.
   }
-  \item{B.mstop}{
-    number of inner bootstrap replicates used to determine the optimal
-    mstop on each of the \code{B} bootstrap samples.
+  \item{papply}{
+    (parallel) apply function for the outer bootstrap, defaults to
+    \code{\link[parallel]{mclapply}} if no inner bootstrap is used to
+    determine the optimal stopping iteration. For details see
+    argument \code{papply} in \code{\link{cvrisk}}. Be careful with your
+    computing resources if you use parallel computing for both, the
+    inner and the outer bootstrap.
   }
+  \item{papply.mstop}{
+    (parallel) apply function for the inner bootstrap, defaults to
+    \code{\link[parallel]{mclapply}} if we use an inner bootstrap is to
+    determine the optimal stopping iteration. For details see
+    argument \code{papply} in \code{\link{cvrisk}}. Be careful with your
+    computing resources if you use parallel computing for both, the
+    inner and the outer bootstrap.
+  }
   \item{x}{
     a confidence interval object.
   }
@@ -62,10 +80,11 @@
     limits of the y scale. Per default computed from the data to plot.
   }
   \item{type}{
-    type of graphic for the point estimate, i.e., the predicted function.
-methods  }
+    type of graphic for the point estimate, i.e., for the predicted
+    function. Per default a line is plotted.
+  }
   \item{col}{
-    color of the point estimate, i.e., the predicted function.
+    color of the point estimate, i.e., for the predicted function.
   }
   \item{ci.col}{
     color of the confidence interval.
@@ -74,6 +93,14 @@
     logical, should the raw function estimates or the derived confidence
     estimates be plotted?
   }
+  \item{print_levelplot}{
+    logical, should the \pkg{lattice} \code{\link{levelplot}} be printed
+    or simply returned for further modifications. This argument is only
+    considered if bivariate effect estimates are plotted. If
+    \code{print_levelplot} is set to \code{FALSE}, a list with objects
+    \code{mean}, \code{lowerPI} and \code{upperPI} is returned
+    containing the three \code{\link{levelplot}} objects.
+  }
   \item{pe}{
     logical, should the point estimtate (PE) be also returned?
   }
@@ -90,9 +117,10 @@
   An object of class \code{glmboost.ci} or \code{mboost.ci} with special
   \code{print} and/or \code{plot} functions.
 }
-\references{
-  %% ~put references to the literature/web site here ~
-}
+%% \references{
+%% Benjamin Hofner , Thomas Kneib and Torsten Hothorn (2014),
+%% A Unified Framework of Constrained Regression.
+%% }
 \author{
   Benjamin Hofner <benjamin.hofner at fau.de>
 }



More information about the Mboost-commits mailing list