[Mboost-commits] r795 - in pkg/mboostDevel: R man tests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Sep 25 17:21:42 CEST 2014


Author: hofner
Date: 2014-09-25 17:21:41 +0200 (Thu, 25 Sep 2014)
New Revision: 795

Modified:
   pkg/mboostDevel/R/mboost.R
   pkg/mboostDevel/man/mboost.Rd
   pkg/mboostDevel/tests/regtest-glmboost.R
Log:
- merged changes from mboostPatch to mboostDevel


Modified: pkg/mboostDevel/R/mboost.R
===================================================================
--- pkg/mboostDevel/R/mboost.R	2014-09-25 14:47:41 UTC (rev 794)
+++ pkg/mboostDevel/R/mboost.R	2014-09-25 15:21:41 UTC (rev 795)
@@ -264,7 +264,13 @@
             } else {
                 ## only if no selection of baselearners
                 ## was made via the `which' argument
-                return(offset + matrix(rowSums(pr), ncol = 1))
+                ret <- matrix(rowSums(pr), ncol = 1)
+                if (length(offset) != 1 && !is.null(newdata)) {
+                    warning("Offset not used for prediction when ", sQuote("newdata"), " is specified")
+                } else {
+                    ret <- ret + offset
+                }
+                return(ret)
             }
         }, "cumsum" = {
             if (!nw) {
@@ -276,10 +282,15 @@
                 attr(pr, "offset") <- offset
                 return(pr)
             } else {
-                ret <- 0
-                for (i in 1:max(xselect)) ret <- ret + pfun(i, agg = "none")
-                return(.Call("R_mcumsum", as(ret, "matrix"), PACKAGE = "mboostDevel")
-                       + offset)
+                pr <- 0
+                for (i in 1:max(xselect)) pr <- pr + pfun(i, agg = "none")
+                pr <- .Call("R_mcumsum", as(pr, "matrix"), PACKAGE = "mboostDevel")
+                if (length(offset) != 1 && !is.null(newdata)) {
+                    warning("Offset not used for prediction when ", sQuote("newdata"), " is specified")
+                } else {
+                    pr <- pr + offset
+                }
+                return(pr)
             }
          }, "none" = {
             if (!nw) {
@@ -289,11 +300,11 @@
                 attr(pr, "offset") <- offset
                 return(pr)
             } else {
-                ret <- 0
-                for (i in 1:max(xselect)) ret <- ret + pfun(i, agg = "none")
-                ret <- as(ret, "matrix")
-                attr(ret, "offset") <- offset
-                return(ret)
+                pr <- 0
+                for (i in 1:max(xselect)) pr <- pr + pfun(i, agg = "none")
+                pr <- as(pr, "matrix")
+                attr(pr, "offset") <- offset
+                return(pr)
             }
          })
          return(pr)

Modified: pkg/mboostDevel/man/mboost.Rd
===================================================================
--- pkg/mboostDevel/man/mboost.Rd	2014-09-25 14:47:41 UTC (rev 794)
+++ pkg/mboostDevel/man/mboost.Rd	2014-09-25 15:21:41 UTC (rev 795)
@@ -137,6 +137,18 @@
        mod <- mboost_fit(list(btree(age), bols(waistcirc), bbs(hipcirc)),
                          response = DEXfat))
   plot(mod, ask = FALSE, main = "base-learner")
+
+  ### use stability selection for base-learner selection
+  library("stabs")
+  ## fit a linear model
+  # (as this is quicker; also possible with non-linear model from above)
+  mod <- glmboost(DEXfat ~ ., data = bodyfat)
+  ## now run stability selection
+  (sbody <- stabsel(mod, q = 3, PFER = 1, sampling.type = "MB"))
+  opar <- par(mai = par("mai") * c(1, 1, 1, 2.7))
+  plot(sbody)
+  par(opar)
+  plot(sbody, type = "maxsel", ymargin = 6)
 }
 \keyword{models}
 \keyword{nonlinear}

Modified: pkg/mboostDevel/tests/regtest-glmboost.R
===================================================================
--- pkg/mboostDevel/tests/regtest-glmboost.R	2014-09-25 14:47:41 UTC (rev 794)
+++ pkg/mboostDevel/tests/regtest-glmboost.R	2014-09-25 15:21:41 UTC (rev 795)
@@ -184,7 +184,10 @@
 
 stopifnot(max(abs(predict(gmod) -  predict(gbmod))) < 1e-4)
 
-### predictions:
+################################################################################
+### Check predictions
+################################################################################
+
 set.seed(1907)
 x1 <- rnorm(100)
 x2 <- rnorm(100)
@@ -195,22 +198,36 @@
 amod <- glmboost(y ~ -1 + x1 + x2, data = DF, center = FALSE)
 agg <- c("none", "sum", "cumsum")
 whi <- list(NULL, 1, 2, c(1,2))
+pred <- vector("list", length = 4)
 for (i in 1:4){
-    pred <- vector("list", length=3)
+    pred[[i]][[i]] <- vector("list", length=3)
     for (j in 1:3){
-        pred[[j]] <- predict(amod, aggregate=agg[j], which = whi[[i]])
+        pred[[i]][[j]] <- predict(amod, aggregate=agg[j], which = whi[[i]])
     }
     if (i == 1){
-        stopifnot(max(abs(pred[[2]] - pred[[3]][,ncol(pred[[3]])]))  < sqrt(.Machine$double.eps))
-        if ((pred[[2]] - rowSums(pred[[1]]))[1] - attr(coef(amod), "offset") < sqrt(.Machine$double.eps))
+        stopifnot(max(abs(pred[[i]][[2]] - pred[[i]][[3]][,ncol(pred[[i]][[3]])]))  < sqrt(.Machine$double.eps))
+        if ((pred[[i]][[2]] - rowSums(pred[[i]][[1]]))[1] - attr(coef(amod), "offset") < sqrt(.Machine$double.eps))
             warning(sQuote("aggregate = sum"), " adds the offset, ", sQuote("aggregate = none"), " doesn't.")
-        stopifnot(max(abs(pred[[2]] - rowSums(pred[[1]]) - attr(coef(amod), "offset")))   < sqrt(.Machine$double.eps))
+        stopifnot(max(abs(pred[[i]][[2]] - rowSums(pred[[i]][[1]]) - attr(coef(amod), "offset")))   < sqrt(.Machine$double.eps))
     } else {
-        stopifnot(max(abs(pred[[2]] - sapply(pred[[3]], function(obj) obj[,ncol(obj)])))  < sqrt(.Machine$double.eps))
-        stopifnot(max(abs(pred[[2]] - sapply(pred[[1]], function(obj) rowSums(obj))))  < sqrt(.Machine$double.eps))
+        stopifnot(max(abs(pred[[i]][[2]] - sapply(pred[[i]][[3]], function(obj) obj[,ncol(obj)])))  < sqrt(.Machine$double.eps))
+        stopifnot(max(abs(pred[[i]][[2]] - sapply(pred[[i]][[1]], function(obj) rowSums(obj))))  < sqrt(.Machine$double.eps))
     }
 }
+## compare which = NULL and which == c(1, 2)
+# type = "none"
+stopifnot(all.equal(pred[[1]][[1]], pred[[4]][[1]]$x1 + pred[[4]][[1]]$x2, check.attributes = FALSE))
+# type = "sum"
+predictions <- as.matrix(DF[, c("x1", "x2")]) %*% matrix(coef(amod), ncol = 1) +
+    attr(coef(amod), "offset")
+stopifnot(all.equal(pred[[1]][[2]], predictions))
+stopifnot(all.equal(c(pred[[1]][[2]]),
+                    rowSums(pred[[4]][[2]]) + attr(coef(amod), "offset"),
+                    check.attributes = FALSE))
+# type = "cumsum"
+stopifnot(all.equal(pred[[1]][[3]], pred[[4]][[3]]$x1 + pred[[4]][[3]]$x2 + attr(coef(amod), "offset")))
 
+## same with names
 agg <- c("none", "sum", "cumsum")
 whi <- list(NULL, "x1", "x2", c("x1","x2"))
 for (i in 1:4){
@@ -250,7 +267,65 @@
 stopifnot(ncol(pr) == 3 || all(pr[,c(1,ncol)] == 0))
 amod[100]
 
-# compare predictions with gamboost
+## check predictions with offset
+amod <- glmboost(y ~ -1 + x1 + x2, offset = 10, data = DF, center=FALSE)
+bmod <- glmboost(y ~ -1 + x1 + x2, offset = fitted(amod), data = DF, center=FALSE)
+## specify newdata
+newdata <- data.frame(x1 = rnorm(10), x2 = rnorm(10), x3 = rnorm(10))
+agg <- c("none", "sum", "cumsum")
+whi <- list(NULL, "x1", "x2", c("x1","x2"))
+
+preda <- predb <- vector("list", length = 4)
+for (i in 1:4){
+    preda[[i]][[i]] <- vector("list", length=3)
+    for (j in 1:3){
+        preda[[i]][[j]] <- predict(amod, aggregate=agg[j], which = whi[[i]],
+                                   newdata = newdata)
+        ## a warning is issued for "sum" and "cumsum"
+        predb[[i]][[j]] <- predict(bmod, aggregate=agg[j], which = whi[[i]],
+                                   newdata = newdata)
+    }
+    if (i == 1){
+        ## checks for amod
+        stopifnot(max(abs(preda[[i]][[2]] - preda[[i]][[3]][,ncol(preda[[i]][[3]])]))  < sqrt(.Machine$double.eps))
+        if ((preda[[i]][[2]] - rowSums(preda[[i]][[1]]))[1] - attr(coef(amod), "offset") < sqrt(.Machine$double.eps))
+            warning(sQuote("aggregate = sum"), " adds the offset, ", sQuote("aggregate = none"), " doesn't.")
+        stopifnot(max(abs(preda[[i]][[2]] - rowSums(preda[[i]][[1]]) - attr(coef(amod), "offset")))   < sqrt(.Machine$double.eps))
+        ## same for bmod
+        stopifnot(max(abs(predb[[i]][[2]] - predb[[i]][[3]][,ncol(predb[[i]][[3]])]))  < sqrt(.Machine$double.eps))
+        if ((predb[[i]][[2]] - rowSums(predb[[i]][[1]]))[1] < sqrt(.Machine$double.eps))
+            warning(sQuote("aggregate = sum"), " adds the offset, ", sQuote("aggregate = none"), " doesn't.")
+        ## here, the offset is not used in both cases!
+        stopifnot(max(abs(predb[[i]][[2]] - rowSums(predb[[i]][[1]])))   < sqrt(.Machine$double.eps))
+    } else {
+        ## checks for amod
+        stopifnot(max(abs(preda[[i]][[2]] - sapply(preda[[i]][[3]], function(obj) obj[,ncol(obj)])))  < sqrt(.Machine$double.eps))
+        stopifnot(max(abs(preda[[i]][[2]] - sapply(preda[[i]][[1]], function(obj) rowSums(obj))))  < sqrt(.Machine$double.eps))
+        ## same for bmod
+        stopifnot(max(abs(predb[[i]][[2]] - sapply(predb[[i]][[3]], function(obj) obj[,ncol(obj)])))  < sqrt(.Machine$double.eps))
+        stopifnot(max(abs(predb[[i]][[2]] - sapply(predb[[i]][[1]], function(obj) rowSums(obj))))  < sqrt(.Machine$double.eps))
+    }
+}
+## compare which = NULL and which == c(1, 2);
+## The offset should be always dropped for bmod, but kept for amod
+# type = "none"
+stopifnot(all.equal(preda[[1]][[1]], preda[[4]][[1]]$x1 + preda[[4]][[1]]$x2, check.attributes = FALSE))
+stopifnot(all.equal(predb[[1]][[1]], predb[[4]][[1]]$x1 + predb[[4]][[1]]$x2, check.attributes = FALSE))
+# type = "sum"
+predictionsA <- as.matrix(newdata[, c("x1", "x2")]) %*% matrix(coef(amod), ncol = 1) +
+    attr(coef(amod), "offset")
+predictionsB <- as.matrix(newdata[, c("x1", "x2")]) %*% matrix(coef(bmod), ncol = 1)
+stopifnot(all.equal(preda[[1]][[2]], predictionsA))
+stopifnot(all.equal(predb[[1]][[2]], predictionsB))
+stopifnot(all.equal(c(preda[[1]][[2]]), rowSums(preda[[4]][[2]]) + attr(coef(amod), "offset"),
+                    check.attributes = FALSE))
+stopifnot(all.equal(c(predb[[1]][[2]]), rowSums(predb[[4]][[2]]),
+                    check.attributes = FALSE))
+# type = "cumsum"
+stopifnot(all.equal(predb[[1]][[3]], predb[[4]][[3]]$x1 + predb[[4]][[3]]$x2))
+stopifnot(all.equal(predb[[1]][[3]], predb[[4]][[3]]$x1 + predb[[4]][[3]]$x2))
+
+### compare predictions with gamboost
 mod1 <- glmboost(y ~ -1 + x1 + x2 + x3, data = DF, center=FALSE)
 mod2 <- gamboost(y ~ x1 + x2 + x3, data = DF, baselearner= function(x) bols(x, intercept=FALSE))
 pr1_2 <- predict(mod2, aggre = "cumsum")
@@ -286,7 +361,10 @@
 pr2 <- predict(logit, type="response")
 stopifnot(pr - pr2 < sqrt(.Machine$double.eps))
 
-### coefficients:
+################################################################################
+### Check coefficients
+################################################################################
+
 set.seed(1907)
 x1 <- rnorm(100)
 x2 <- rnorm(100)



More information about the Mboost-commits mailing list