[Mboost-commits] r803 - in pkg: mboostDevel mboostDevel/inst mboostDevel/man mboostPatch mboostPatch/R mboostPatch/inst mboostPatch/man mboostPatch/tests mboostPatch/tests/Examples mboostPatch/vignettes

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Oct 2 16:02:14 CEST 2014


Author: hofner
Date: 2014-10-02 16:02:14 +0200 (Thu, 02 Oct 2014)
New Revision: 803

Added:
   pkg/mboostPatch/R/confint.R
   pkg/mboostPatch/R/stabsel.R
   pkg/mboostPatch/inst/NEWS.Rd
   pkg/mboostPatch/man/confint.Rd
   pkg/mboostPatch/man/stabsel.Rd
   pkg/mboostPatch/tests/regtest-inference.R
   pkg/mboostPatch/tests/regtest-inference.Rout.save
Removed:
   pkg/mboostPatch/NEWS
   pkg/mboostPatch/inst/CHANGES
Modified:
   pkg/mboostDevel/DESCRIPTION
   pkg/mboostDevel/inst/NEWS.Rd
   pkg/mboostDevel/man/mboost_package.Rd
   pkg/mboostPatch/.Rbuildignore
   pkg/mboostPatch/.RbuildignoreCRAN
   pkg/mboostPatch/DESCRIPTION
   pkg/mboostPatch/NAMESPACE
   pkg/mboostPatch/R/bmrf.R
   pkg/mboostPatch/R/brad.R
   pkg/mboostPatch/R/btree.R
   pkg/mboostPatch/R/helpers.R
   pkg/mboostPatch/R/mboost.R
   pkg/mboostPatch/R/methods.R
   pkg/mboostPatch/R/plot.R
   pkg/mboostPatch/inst/CITATION
   pkg/mboostPatch/man/baselearners.Rd
   pkg/mboostPatch/man/blackboost.Rd
   pkg/mboostPatch/man/mboost_package.Rd
   pkg/mboostPatch/man/methods.Rd
   pkg/mboostPatch/tests/Examples/mboost-Ex.Rout.save
   pkg/mboostPatch/tests/birds_Biometrics.Rout.save
   pkg/mboostPatch/tests/bugfixes.Rout.save
   pkg/mboostPatch/tests/regtest-baselearner.Rout.save
   pkg/mboostPatch/tests/regtest-blackboost.Rout.save
   pkg/mboostPatch/tests/regtest-family.Rout.save
   pkg/mboostPatch/tests/regtest-gamboost.Rout.save
   pkg/mboostPatch/tests/regtest-glmboost.Rout.save
   pkg/mboostPatch/tests/regtest-hatmatrix.Rout.save
   pkg/mboostPatch/vignettes/SurvivalEnsembles.Rout.save
   pkg/mboostPatch/vignettes/mboost.Rout.save
   pkg/mboostPatch/vignettes/mboost_illustrations.Rout.save
   pkg/mboostPatch/vignettes/mboost_tutorial.Rout.save
Log:
- merged changes of minor release (mboostDevel) to new patch version (mboostPatch)
- updated version numbers in both versions


Modified: pkg/mboostDevel/DESCRIPTION
===================================================================
--- pkg/mboostDevel/DESCRIPTION	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostDevel/DESCRIPTION	2014-10-02 14:02:14 UTC (rev 803)
@@ -1,7 +1,7 @@
 Package: mboostDevel
 Title: Model-Based Boosting
-Version: 2.4-0
-Date: 2014-10-02
+Version: 2.5-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/inst/NEWS.Rd
===================================================================
--- pkg/mboostDevel/inst/NEWS.Rd	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostDevel/inst/NEWS.Rd	2014-10-02 14:02:14 UTC (rev 803)
@@ -1,6 +1,25 @@
 \name{NEWS}
 \title{News for Package 'mboost'}
 
+\section{Changes in mboost version 2.5-0 (2014-xx-yy)}{
+  \subsection{User-visible changes}{
+    \itemize{
+      \item 
+    }
+  }
+  \subsection{Miscellaneous}{
+    \itemize{
+      \item 
+    }
+  }
+  \subsection{Bug-fixes}{
+    \itemize{
+      \item 
+    }
+  }
+}
+	
+
 \section{Changes in mboost version 2.4-0 (2014-10-02)}{
   \subsection{User-visible changes}{
     \itemize{

Modified: pkg/mboostDevel/man/mboost_package.Rd
===================================================================
--- pkg/mboostDevel/man/mboost_package.Rd	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostDevel/man/mboost_package.Rd	2014-10-02 14:02:14 UTC (rev 803)
@@ -18,8 +18,8 @@
 \tabular{ll}{
 Package: \tab mboost\cr
 Type: \tab Package\cr
-Version: \tab 2.4-0\cr
-Date: \tab 2014-10-02\cr
+Version: \tab 2.5-0\cr
+Date: \tab 2014-xx-yy\cr
 License: \tab GPL-2\cr
 LazyLoad: \tab yes\cr
 LazyData: \tab yes\cr

Modified: pkg/mboostPatch/.Rbuildignore
===================================================================
--- pkg/mboostPatch/.Rbuildignore	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/.Rbuildignore	2014-10-02 14:02:14 UTC (rev 803)
@@ -2,3 +2,5 @@
 to_do_list.txt
 ^\..*
 .*/auto
+inst/NEWS.html
+inst/NEWS.pdf

Modified: pkg/mboostPatch/.RbuildignoreCRAN
===================================================================
--- pkg/mboostPatch/.RbuildignoreCRAN	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/.RbuildignoreCRAN	2014-10-02 14:02:14 UTC (rev 803)
@@ -4,3 +4,5 @@
 ^\..*
 .*/auto
 vignettes/.*\.Rout\.save$
+inst/NEWS.html
+inst/NEWS.pdf
\ No newline at end of file

Modified: pkg/mboostPatch/DESCRIPTION
===================================================================
--- pkg/mboostPatch/DESCRIPTION	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/DESCRIPTION	2014-10-02 14:02:14 UTC (rev 803)
@@ -1,6 +1,6 @@
 Package: mboost
 Title: Model-Based Boosting
-Version: 2.3-1
+Version: 2.4-1
 Date: 2014-xx-yy
 Authors at R: c(person("Torsten", "Hothorn", role = c("aut", "cre"),
                     email = "Torsten.Hothorn at R-project.org"),
@@ -15,10 +15,10 @@
   component-wise (penalised) least squares estimates or regression
   trees as base-learners for fitting generalized linear, additive
   and interaction models to potentially high-dimensional data.
-Depends: R (>= 2.14.0), methods, stats, parallel
+Depends: R (>= 2.14.0), methods, stats, parallel, stabs
 Imports: Matrix, survival, splines, lattice, nnls, quadprog, utils
 Suggests: party (>= 1.0-3), TH.data, MASS, fields, BayesX, gbm, mlbench,
-        RColorBrewer, rpart (>= 4.0-3), stabs
+        RColorBrewer, rpart (>= 4.0-3)
 LazyData: yes
 License: GPL-2
 URL: http://r-forge.r-project.org/projects/mboost/

Modified: pkg/mboostPatch/NAMESPACE
===================================================================
--- pkg/mboostPatch/NAMESPACE	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/NAMESPACE	2014-10-02 14:02:14 UTC (rev 803)
@@ -1,8 +1,10 @@
 
 import(methods)
 import(stats)
+import(stabs)
 
 import(Matrix)
+import(parallel)
 importFrom(survival, Surv, survfit)
 importFrom(splines, bs, splineDesign)
 importFrom(lattice, levelplot)
@@ -20,7 +22,8 @@
        Huber, AdaExp, Gehan, CoxPH, Hurdle, Multinomial, FP, IPCweights,
        cvrisk, cv, bbs,
        bols, bspatial, brandom, btree, bss, bns, brad, bmono, bmrf, buser, survFit, selected,
-       nuisance, "%+%", "%X%", "%O%", extract, risk, "mstop<-")
+       nuisance, "%+%", "%X%", "%O%", extract, risk, "mstop<-",
+       stabsel.mboost, stabsel_parameters.mboost, confint.mboost, confint.glmboost)
        ###, basesel, fitsel)
 exportClasses("boost_family")
 exportMethods("show")
@@ -78,5 +81,12 @@
 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)
+S3method(stabsel, mboost)
+S3method(stabsel_parameters, mboost)
 
 useDynLib(mboost)

Deleted: pkg/mboostPatch/NEWS
===================================================================
--- pkg/mboostPatch/NEWS	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/NEWS	2014-10-02 14:02:14 UTC (rev 803)
@@ -1 +0,0 @@
-link inst/CHANGES
\ No newline at end of file

Modified: pkg/mboostPatch/R/bmrf.R
===================================================================
--- pkg/mboostPatch/R/bmrf.R	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/R/bmrf.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -53,13 +53,12 @@
     return(ret)
 }
 
-hyper_bmrf <-
-function (mf, vary, bnd = NULL, df = 4, lambda = NULL, center = FALSE)
-{
+hyper_bmrf <- function (mf, vary, bnd = NULL, df = 4, lambda = NULL,
+                        center = FALSE) {
     if (is.null(bnd))
         stop("Neighbourhood relations must be given in matrix or boundary format.")
     else if (inherits(bnd, "bnd"))
-        K <- bnd2gra(bnd)
+        K <- BayesX::bnd2gra(bnd)
     else if (isMATRIX(bnd) &&
              nrow(bnd) == ncol(bnd) &&
              nlevels(mf[[1]]) <= nrow(bnd) &&
@@ -75,9 +74,7 @@
         center = center)
 }
 
-X_bmrf <-
-function (mf, vary, args)
-{
+X_bmrf <- function (mf, vary, args) {
     K <- args$K
     districts <- rownames(K)
     X <- Diagonal(nrow(K))

Modified: pkg/mboostPatch/R/brad.R
===================================================================
--- pkg/mboostPatch/R/brad.R	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/R/brad.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -1,9 +1,9 @@
 brad <- function(..., by = NULL, index = NULL, knots = 100, df = 4, lambda = NULL,
-                 covFun = stationary.cov,
+                 covFun = fields::stationary.cov,
                  args = list(Covariance = "Matern", smoothness = 1.5, theta = NULL)) {
 
-    if (!require("fields"))
-        stop("cannot load ", sQuote("fields"))
+    # if (!require("fields"))
+    #     stop("cannot load ", sQuote("fields"))
 
     cll <- match.call()
     cll[[1]] <- as.name("brad")
@@ -116,7 +116,10 @@
     ## first we need to build a correct matrix of mf
     x <- as.matrix(mf[which(colnames(mf) != vary)])
     if (length(knots) == 1) {
-        knots <- cover.design(R = unique(x), nd = knots)$design
+        if (!require("fields"))
+            stop("Cannot load package", sQuote("fields"),
+                 ", which is needed for the automatic knot placement")
+        knots <- fields::cover.design(R = unique(x), nd = knots)$design
     }
     if ("theta" %in% names(args) && is.null(args$theta)){
         ## (try to) compute effective range
@@ -135,16 +138,16 @@
 #       rho( max(x_(i) - x_(j)), smoothness, theta = max(x_(i) - x_(j))/c ) = 0.001
 #  <==> rho(c, smoothness, theta = 1) = 0.001
 effective_range <- function(x, eps = 0.001, interval = c(0.1, 100),
-                            covFun = stationary.cov, args = list()){
+                            covFun = fields::stationary.cov, args = list()){
 
-    if ( !( length(deparse(covFun)) == length(deparse(stationary.cov))
-           && all(deparse(covFun) == deparse(stationary.cov)) ) &&
-         !( length(deparse(covFun)) == length(deparse(Exp.cov))
-           && all(deparse(covFun) == deparse(Exp.cov)) ) ){
+    if ( !( length(deparse(covFun)) == length(deparse(fields::stationary.cov))
+           && all(deparse(covFun) == deparse(fields::stationary.cov)) ) &&
+         !( length(deparse(covFun)) == length(deparse(fields::Exp.cov))
+           && all(deparse(covFun) == deparse(fields::Exp.cov)) ) ){
         ## if cov.funcion is not one of stationary.cov and Exp.cov
         warning(sQuote("effective_range()"), " is only implemented for ",
                 sQuote("stationary.cov"), " and ", sQuote("Exp.cov"),
-                " from package:fields.")
+                " from package ", sQuote("fields"))
         return(NULL)
     }
 

Modified: pkg/mboostPatch/R/btree.R
===================================================================
--- pkg/mboostPatch/R/btree.R	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/R/btree.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -1,8 +1,8 @@
 
 ### the classical tree-based baselearner; stumps by default
 ### (also fits an additive model)
-btree <- function(..., 
-    tree_controls = ctree_control(stump = TRUE, 
+btree <- function(...,
+    tree_controls = party::ctree_control(stump = TRUE,
                                   mincriterion = 0,
                                   savesplitstats = FALSE)) {
 
@@ -52,7 +52,7 @@
         df <- mf
         df[[rname]] <- y
         object <- party:::ctreedpp(fm, data = df)
-        fitmem <- ctree_memory(object, TRUE)
+        fitmem <- party::ctree_memory(object, TRUE)
         where <- rep.int(0, nrow(mf))
         storage.mode(where) <- "integer"
         storage.mode(weights) <- "double"

Added: pkg/mboostPatch/R/confint.R
===================================================================
--- pkg/mboostPatch/R/confint.R	                        (rev 0)
+++ pkg/mboostPatch/R/confint.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -0,0 +1,326 @@
+
+confint.mboost <- function(object, parm = NULL, level = 0.95,
+                           B = 1000, B.mstop = 25, newdata = NULL,
+                           which = parm,
+                           papply = ifelse(B.mstop == 0, mclapply, lapply),
+                           cvrisk_options = list(), ...) {
+
+    which <- object$which(which, usedonly = FALSE)
+    if (!all(which %in% object$which(NULL, usedonly = FALSE)))
+        stop(sQuote("which"), " is wrongly specified")
+
+    if (!is.list(cvrisk_options))
+        stop(sQuote("cvrisk_options"), " must be a named list")
+    if (length(cvrisk_options) > 0 && is.null(names(cvrisk_options)))
+        stop(sQuote("cvrisk_options"), " must be a named list")
+    if ("folds" %in% names(cvrisk_options))
+        stop("One cannot modify the folds of the inner bootstrap")
+    if ("object" %in% names(cvrisk_options))
+        stop("One cannot specify the model (object) of the inner bootstrap")
+
+    ## create new data and/or restructure data
+    newdata <- .create_newdata(object, newdata, which)
+
+    outer.folds <- cv(model.weights(object), B = B)
+
+    cat("Start computing bootstrap confidence intervals... \n")
+
+    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 <- do.call("cvrisk",
+                           args = c(list(object = mod,
+                                       folds = cv(model.weights(mod), B = B.mstop)),
+                                       cvrisk_options))
+            mod[mstop(cvr)]
+        }
+        .predict_confint(mod, newdata = newdata, which = which)
+    }
+    predictions <- papply(1:B, do_update, ...)
+    cat("\n")
+
+    ## prepare returned object
+    res <- list(level = level, boot_pred = predictions, data = newdata,
+                model = object)
+    attr(res, "which") <- which
+    class(res) <- "mboost.ci"
+    return(res)
+}
+
+confint.glmboost <- function(object, parm = NULL, level = 0.95,
+                             B = 1000, B.mstop = 25,
+                             which = parm, ...) {
+
+    outer.folds <- cv(model.weights(object), B = B)
+    which <- object$which(which, usedonly = FALSE)
+    if (!all(which %in% object$which(NULL, usedonly = FALSE)))
+        stop(sQuote("which"), " is wrongly specified")
+
+    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))
+    }
+
+    ## prepare returned object
+    res <- list(confint = .ci_glmboost(coefficients, level = level, which = which),
+                level = level, boot_coefs = coefficients, model = object)
+    attr(res, "which") <- which
+    class(res) <- "glmboost.ci"
+    return(res)
+}
+
+.ci_glmboost <- function(coefficients, level, which = NULL) {
+    quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
+
+    tmp <- apply(coefficients, 2, FUN = quantile, probs = quantiles)
+    CI <- as.data.frame(t(tmp))[which, ]
+    return(CI)
+}
+
+## pe = add point estimte
+print.glmboost.ci <- function(x, which = NULL, level = x$level, pe = FALSE, ...) {
+
+    if (is.null(which)) {
+        which <- attr(x, "which")
+    } else {
+        which <- x$model$which(which, usedonly = FALSE)
+        if (!all(which %in% attr(x, "which")))
+            stop(sQuote("which"), " is wrongly specified")
+    }
+
+    if (!is.null(level) && level != x$level) {
+        CI <- .ci_glmboost(x$boot_coefs,  level = level, which = which)
+    } else {
+        CI <- x$confint[which, ]
+    }
+
+    if (pe) {
+        tmp <- data.frame(beta = coef(x$model, which))
+        CI <- cbind(tmp, CI)
+    }
+    if (length(which) > 1) {
+        cat("\tBootstrap Confidence Intervals\n")
+    } else {
+        cat("\tBootstrap Confidence Interval\n")
+    }
+    print(CI, ...)
+    return(invisible(x))
+}
+
+
+## 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)
+        nms <- names(object$baselearner)[which]
+
+        for (w in which) {
+            ## get data from w-th base-learner
+            tmp <- data[[w]][rep(1, 100), , drop = FALSE]
+
+            ## are there varying coefficients (i.e. argument by)
+            get_vary <- object$baselearner[[w]]$get_vary
+            vary <- ""
+            if (!is.null(get_vary)) vary <- get_vary()
+            if (vary != "") {
+                if (is.factor(tmp[[vary]])) {
+                    if (nlevels(tmp[[vary]]) > 2)
+                        stop("Atomatic data creation for ", sQuote("by"),
+                             " variables with more than two levels is",
+                             " currently not supported;",
+                             " Specify ", sQuote("newdata"), " instead.")
+                    data[[w]][[vary]] <- factor(levels(data[[w]][[vary]])[-1],
+                                                levels = levels(data[[w]][[vary]]))
+                }
+                if (is.numeric(tmp[[vary]]))
+                    data[[w]][[vary]] <- 1
+            }
+
+            ## now make grid
+            grid <- function(x) {
+                if (is.numeric(x)) {
+                    return(seq(min(x), max(x), length = 100))
+                } else {
+                    return(rep(unique(x), length.out = 100))
+                }
+            }
+
+            for (j in 1:ncol(data[[w]]))
+                tmp[, colnames(data[[w]])[j]] <- grid(data[[w]][,j])
+
+            ## FIXME: what about btree and bmrf?
+
+            ## check if any base-learner is a bivariate smooth effect, i.e. if
+            ## base-learner is multivariate and bbs, bspatial or brad
+            which.vary <- colnames(tmp) == vary
+            multivar <- grepl("bbs|bspatial|brad", nms[w]) &
+                ncol(tmp[!which.vary]) >= 2
+            if (multivar) {
+                ## make grid
+                egrid <- expand.grid(tmp[!which.vary])
+                if (vary != "") {
+                    x.vary <- tmp[which.vary]
+                    rownames(x.vary) <- NULL
+                    tmp <- cbind(egrid, x.vary)
+                } else {
+                    tmp <- egrid
+                }
+            }
+            ## reset rownames
+            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)
+}
+
+## special prediction function for the construction of confidence intervals:
+.predict_confint <- function(object, newdata = NULL, which, ...) {
+    nrows <- sapply(newdata, nrow)
+    predictions <- sapply(which, function(w)
+                          matrix(predict(object, newdata[[w]], which = w),
+                                               ncol = 1, nrow = nrows[w]))
+    if (is.matrix(predictions))
+        predictions <- as.data.frame(predictions)
+    names(predictions) <- names(newdata[which])
+    return(predictions)
+}
+
+
+### plot functions
+plot.mboost.ci <- function(x, which, level = x$level,
+                           ylim = NULL, type = "l", col = "black",
+                           ci.col = rgb(170, 170, 170, alpha = 85,
+                                        maxColorValue = 255),
+                           raw = FALSE, print_levelplot = TRUE, ...) {
+
+    if (missing(which)) {
+        which <- attr(x, "which")
+    } else {
+        which <- x$model$which(which, usedonly = FALSE)
+        if (!all(which %in% attr(x, "which")))
+            stop(sQuote("which"), " is wrongly specified")
+    }
+    if (length(which) > 1)
+        stop("Specify a single base-learner using ", sQuote("which"))
+
+    CI <- .ci_mboost(x$boot_pred, level = level, which = which, raw = raw)
+
+    ## 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, ...)
+    }
+}
+
+lines.mboost.ci <- function(x, which, level = x$level,
+                            col = rgb(170, 170, 170, alpha = 85,
+                                      maxColorValue = 255),
+                            raw = FALSE, ...) {
+
+    if (missing(which)) {
+        which <- attr(x, "which")
+    } else {
+        which <- x$model$which(which, usedonly = FALSE)
+        if (!all(which %in% attr(x, "which")))
+            stop(sQuote("which"), " is wrongly specified")
+    }
+    if (length(which) > 1)
+        stop("Specify a single base-learner using ", sQuote("which"))
+
+
+    CI <- .ci_mboost(x$boot_pred, level = level, which = which, raw = raw)
+
+    x.data <- x$data[[which]]
+    ## 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 (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 {
+        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", ...)
+        }
+    }
+}
+
+.ci_mboost <- function(predictions, level, which = NULL, raw = FALSE) {
+
+    preds <- sapply(predictions, function(p) p[[which]])
+    if (!raw) {
+        quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
+        preds <- apply(preds, 1, FUN = quantile, probs = quantiles)
+    }
+
+    return(preds)
+}

Modified: pkg/mboostPatch/R/helpers.R
===================================================================
--- pkg/mboostPatch/R/helpers.R	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/R/helpers.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -233,4 +233,4 @@
     cf <- solve.QP(Dmat = XtX, dvec = as.vector(Xty), Amat = t(D),
                    bvec = rep(0, nrow(D)))$solution
     cf
-}
\ No newline at end of file
+}

Modified: pkg/mboostPatch/R/mboost.R
===================================================================
--- pkg/mboostPatch/R/mboost.R	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/R/mboost.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -174,9 +174,14 @@
     )
 
     ### update to new weights; just a fresh start
-    RET$update <- function(weights = NULL, oobweights = NULL, risk = "oobag") {
+    RET$update <- function(weights = NULL, oobweights = NULL, risk = "oobag",
+                           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)
         if (!is.null(offsetarg)) offsetarg <- offset
         mboost_fit(blg = blg, response = response, weights = weights,
@@ -523,7 +528,7 @@
 
 ### just one single tree-based baselearner
 blackboost <- function(formula, data = list(),
-    tree_controls = ctree_control(teststat = "max",
+    tree_controls = party::ctree_control(teststat = "max",
                                testtype = "Teststatistic",
                                mincriterion = 0,
                                maxdepth = 2, savesplitstats = FALSE),
@@ -618,9 +623,11 @@
     ### save standard update function for re-use
     update <- ret$update
     ### needs a specialized update function as well
-    ret$update <- function(weights = NULL, oobweights = NULL, risk = "oobag") {
+    ret$update <- function(weights = NULL, oobweights = NULL, risk = "oobag",
+                           trace = NULL) {
         ## call standard update function
-        res <- update(weights = weights, oobweights = oobweights, risk = risk)
+        res <- update(weights = weights, oobweights = oobweights, risk = risk,
+                      trace = trace)
         ## now re-set all special arguments
         res$newX <- newX
         res$assign <- assign
@@ -715,9 +722,11 @@
     ### save standard update function for re-use
     update <- ret$update
     ### needs a specialized update function as well
-    ret$update <- function(weights = NULL, oobweights = NULL, risk = "oobag") {
+    ret$update <- function(weights = NULL, oobweights = NULL, risk = "oobag",
+                           trace = NULL) {
         ## call standard update function
-        res <- update(weights = weights, oobweights = oobweights, risk = risk)
+        res <- update(weights = weights, oobweights = oobweights, risk = risk,
+                      trace = trace)
         ## now re-set all special arguments
         ret$newX <- newX
         res$assign <- assign

Modified: pkg/mboostPatch/R/methods.R
===================================================================
--- pkg/mboostPatch/R/methods.R	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/R/methods.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -233,8 +233,10 @@
 
 mstop.mboost <- function(object, ...) object$mstop()
 
-update.mboost <- function(object, weights, ...)
-    object$update(weights)
+update.mboost <- function(object, weights, oobweights = NULL,
+                          risk = NULL, trace = NULL, ...)
+    object$update(weights = weights, oobweights = oobweights,
+                  risk = risk, trace = trace)
 
 model.frame.mboost <- function(formula, ...)
     formula$model.frame(...)
@@ -449,9 +451,6 @@
 selected.mboost <- function(object, ...)
     object$xselect()
 
-selected.stabsel <- function(object, ...)
-    object$selected
-
 summary.mboost <- function(object, ...) {
 
     ret <- list(object = object, selprob = NULL)

Modified: pkg/mboostPatch/R/plot.R
===================================================================
--- pkg/mboostPatch/R/plot.R	2014-10-02 09:05:36 UTC (rev 802)
+++ pkg/mboostPatch/R/plot.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -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, ...)
@@ -128,7 +150,7 @@
                         ## xlab not user specified
                         plot_helper(paste(xl, "=", v[i]), yl)
                     } else {
-                        plot_helper(paste(xl, "(", vary, "=", v[i], ")"), yl)
+                        plot_helper(xl, yl)
                     }
                 }
             }

Added: pkg/mboostPatch/R/stabsel.R
===================================================================
--- pkg/mboostPatch/R/stabsel.R	                        (rev 0)
+++ pkg/mboostPatch/R/stabsel.R	2014-10-02 14:02:14 UTC (rev 803)
@@ -0,0 +1,87 @@
+## stabsel method for mboost; requires stabs
+stabsel.mboost <- function(x, cutoff, q, PFER,
+                    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 <- length(variable.names(x))
+    ibase <- 1:p
+
+    sampling.type <- match.arg(sampling.type)
+    if (sampling.type == "MB")
+        assumption <- "none"
+    else
+        assumption <- match.arg(assumption)
+
+    B <- ncol(folds)
+
+    pars <- stabsel_parameters(p = p, cutoff = cutoff, q = q,
+                               PFER = PFER, B = B,
+                               verbose = verbose, sampling.type = sampling.type,
+                               assumption = assumption, FWER = FWER)
+    ## return parameter combination only if eval == FALSE
+    if (!eval)
+        return(pars)
+
+    cutoff <- pars$cutoff
+    q <- pars$q
+    PFER <- pars$PFER
+
+    fun <- function(model) {
+        xs <- selected(model)
+        qq <- sapply(1:length(xs), function(x) length(unique(xs[1:x])))
+        xs[qq > q] <- xs[1]
+        xs
+    }
+    if (sampling.type == "SS") {
+        ## use complementary pairs
+        folds <- cbind(folds, model.weights(x) - folds)
+    }
+    ss <- cvrisk(x, fun = fun,
+                 folds = folds,
+                 papply = papply, ...)
+
+    if (verbose){
+        qq <- sapply(ss, function(x) length(unique(x)))
+        sum_of_violations <- sum(qq < q)
+        if (sum_of_violations > 0)
+            warning(sQuote("mstop"), " too small in ",
+                    sum_of_violations, " of the ", ncol(folds),
+                    " subsampling replicates to select ", sQuote("q"),
+                    " base-learners; Increase ", sQuote("mstop"),
+                    " bevor applying ", sQuote("stabsel"))
+    }
+
+
+    ## if grid specified in '...'
+    if (length(list(...)) >= 1 && "grid" %in% names(list(...))) {
+        m <- max(list(...)$grid)
+    } else {
+        m <- mstop(x)
+    }
+    ret <- matrix(0, nrow = length(ibase), ncol = m)
+    for (i in 1:length(ss)) {
+        tmp <- sapply(ibase, function(x)
+            ifelse(x %in% ss[[i]], which(ss[[i]] == x)[1], m + 1))
+        ret <- ret + t(sapply(tmp, function(x) c(rep(0, x - 1), rep(1, m - x + 1))))
+    }
+
+    phat <- ret / length(ss)
+    rownames(phat) <- names(variable.names(x))
+    if (extends(class(x), "glmboost"))
+        rownames(phat) <- variable.names(x)
+    ret <- list(phat = phat, selected = which((mm <- apply(phat, 1, max)) >= cutoff),
+                max = mm, cutoff = cutoff, q = q, PFER = PFER,
+                sampling.type = sampling.type, assumption = assumption,
+                call = cll)
+    ret$call[[1]] <- as.name("stabsel")
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/mboost -r 803


More information about the Mboost-commits mailing list