[Mboost-commits] r793 - in pkg/mboostPatch: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 25 16:37:22 CEST 2014
Author: hofner
Date: 2014-09-25 16:37:22 +0200 (Thu, 25 Sep 2014)
New Revision: 793
Modified:
pkg/mboostPatch/R/mboost.R
pkg/mboostPatch/man/mboost.Rd
pkg/mboostPatch/tests/regtest-glmboost.R
Log:
- fixed behavior of predict for models with non-scalar offset
- added an example using "stabs" to manual
Modified: pkg/mboostPatch/R/mboost.R
===================================================================
--- pkg/mboostPatch/R/mboost.R 2014-09-22 17:13:50 UTC (rev 792)
+++ pkg/mboostPatch/R/mboost.R 2014-09-25 14:37:22 UTC (rev 793)
@@ -259,7 +259,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) {
@@ -271,10 +277,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 = "mboost")
- + offset)
+ pr <- 0
+ for (i in 1:max(xselect)) pr <- pr + pfun(i, agg = "none")
+ pr <- .Call("R_mcumsum", as(pr, "matrix"), PACKAGE = "mboost")
+ 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) {
@@ -284,11 +295,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/mboostPatch/man/mboost.Rd
===================================================================
--- pkg/mboostPatch/man/mboost.Rd 2014-09-22 17:13:50 UTC (rev 792)
+++ pkg/mboostPatch/man/mboost.Rd 2014-09-25 14:37:22 UTC (rev 793)
@@ -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/mboostPatch/tests/regtest-glmboost.R
===================================================================
--- pkg/mboostPatch/tests/regtest-glmboost.R 2014-09-22 17:13:50 UTC (rev 792)
+++ pkg/mboostPatch/tests/regtest-glmboost.R 2014-09-25 14:37:22 UTC (rev 793)
@@ -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