[Mboost-commits] r755 - in pkg/mboostDevel: . R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 6 15:31:10 CET 2013
Author: hofner
Date: 2013-11-06 15:31:09 +0100 (Wed, 06 Nov 2013)
New Revision: 755
Modified:
pkg/mboostDevel/NAMESPACE
pkg/mboostDevel/R/inference.R
pkg/mboostDevel/R/methods.R
pkg/mboostDevel/tests/regtest-inference.R
Log:
- improved / corrected computation of error bounds / parameters of stabsel
- code of stabsel refactored: parameters are now computed in a specific,
user-visible function stabsel_parameters() which could also be used in
conjunction with other machine learning algorithms
- added selected.stabsel method
Modified: pkg/mboostDevel/NAMESPACE
===================================================================
--- pkg/mboostDevel/NAMESPACE 2013-10-30 17:56:05 UTC (rev 754)
+++ pkg/mboostDevel/NAMESPACE 2013-11-06 14:31:09 UTC (rev 755)
@@ -16,7 +16,8 @@
boost_control, mstop, Family,
GaussReg, Gaussian, GaussClass, Laplace, Binomial, Poisson, GammaReg, QuantReg,
ExpectReg, NBinomial, PropOdds, Weibull, Loglog, Lognormal, AUC, mboost_fit,
- Huber, AdaExp, Gehan, CoxPH, Hurdle, Multinomial, FP, IPCweights, cvrisk, cv, bbs, stabsel,
+ Huber, AdaExp, Gehan, CoxPH, Hurdle, Multinomial, FP, IPCweights,
+ cvrisk, cv, bbs, stabsel, stabsel_parameters,
bols, bspatial, brandom, btree, bss, bns, brad, bmono, bmrf, buser, survFit, selected,
nuisance, "%+%", "%X%", "%O%", extract, "mstop<-")
###, basesel, fitsel)
@@ -69,7 +70,9 @@
# S3method(selected, glmboost)
S3method(update, mboost)
S3method(print, stabsel)
+S3method(print, stabsel_parameters)
S3method(plot, stabsel)
+S3method(selected, stabsel)
S3method(extract, mboost)
S3method(extract, glmboost)
S3method(extract, blackboost)
Modified: pkg/mboostDevel/R/inference.R
===================================================================
--- pkg/mboostDevel/R/inference.R 2013-10-30 17:56:05 UTC (rev 754)
+++ pkg/mboostDevel/R/inference.R 2013-11-06 14:31:09 UTC (rev 755)
@@ -11,6 +11,69 @@
error.bound <- match.arg(error.bound)
B <- ncol(folds)
+ pars <- stabsel_parameters(cutoff = cutoff, q = q,
+ PFER = PFER, p = p, B = B,
+ verbose = verbose, error.bound = error.bound)
+ 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 (error.bound == "SS") {
+ ## use complementary pairs
+ folds <- cbind(folds, model.weights(object) - folds)
+ }
+ ss <- cvrisk(object, 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(object)
+ }
+ 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(object))
+ if (extends(class(object), "glmboost"))
+ rownames(phat) <- variable.names(object)
+ ret <- list(phat = phat, selected = which((mm <- apply(phat, 1, max)) >= cutoff),
+ max = mm, cutoff = cutoff, q = q, PFER = PFER, error.bound = error.bound)
+ class(ret) <- "stabsel"
+ ret
+}
+
+stabsel_parameters <- function(cutoff, q, PFER, p,
+ B = ifelse(error.bound == "MB", 100, 50),
+ verbose = FALSE, error.bound = c("MB", "SS"),
+ FWER) {
+
+ error.bound <- match.arg(error.bound)
+
## only two of the four arguments can be specified
if ((nmiss <- sum(missing(PFER), missing(cutoff),
missing(q), missing(FWER))) != 2) {
@@ -51,14 +114,10 @@
if (missing(cutoff)) {
if (error.bound == "MB") {
- cutoff <- min(0.9, tmp <- (q^2 / (PFER * p) + 1) / 2)
+ cutoff <- min(1, tmp <- (q^2 / (PFER * p) + 1) / 2)
upperbound <- q^2 / p / (2 * cutoff - 1)
} else {
- objective_cf <- function(cutoff) {
- PFER / p - minD(q, p, cutoff, B)
- }
- root <- uniroot(objective_cf, lower = 0.5, upper = 0.9)$root
- cutoff <- min(0.9, tmp <- floor(root * 2 * B) / (2* B))
+ cutoff <- optimal_cutoff(p, q, PFER, B)
upperbound <- minD(q, p, cutoff, B) * p
}
upperbound <- signif(upperbound, 3)
@@ -71,16 +130,10 @@
if (missing(q)) {
if (error.bound == "MB") {
- q <- ceiling(sqrt(PFER * (2 * cutoff - 1) * p))
+ q <- floor(sqrt(PFER * (2 * cutoff - 1) * p))
upperbound <- q^2 / p / (2 * cutoff - 1)
} else {
- objective_q <- function(q) {
- PFER / p - minD(q, p, cutoff, B)
- }
- root <- uniroot(objective_q, lower = 1,
- upper = min(sqrt((B - 1) / (2 * B) * p^2),
- (B - 1) / (2 * B) * p))$root
- q <- ceiling(root)
+ q <- optimal_q(p, cutoff, PFER, B)
upperbound <- minD(q, p, cutoff, B) * p
}
upperbound <- signif(upperbound, 3)
@@ -98,56 +151,14 @@
}
upperbound <- signif(upperbound, 3)
}
+
if (verbose && PFER >= p)
warning("Upper bound for PFER larger than the number of base-learners.")
- 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 (error.bound == "SS") {
- ## use complementary pairs
- folds <- cbind(folds, model.weights(object) - folds)
- }
- ss <- cvrisk(object, 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(object)
- }
- 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(object))
- if (extends(class(object), "glmboost"))
- rownames(phat) <- variable.names(object)
- ret <- list(phat = phat, selected = which((mm <- apply(phat, 1, max)) >= cutoff),
- max = mm, cutoff = cutoff, q = q, PFER = upperbound)
- class(ret) <- "stabsel"
- ret
+ res <- list(cutoff = cutoff, q = q, PFER = upperbound,
+ error.bound = error.bound)
+ class(res) <- "stabsel_parameters"
+ res
}
print.stabsel <- function(x, decreasing = FALSE, ...) {
@@ -161,9 +172,20 @@
}
cat("\nSelection probabilities:\n")
print(sort(x$max[x$max > 0], decreasing = decreasing))
- cat("\nCutoff: ", x$cutoff, "; ", sep = "")
+ cat("\n")
+ print.stabsel_parameters(x)
+ cat("\n")
+ invisible(x)
+}
+
+print.stabsel_parameters <- function(x, ...) {
+ cat("Cutoff: ", x$cutoff, "; ", sep = "")
cat("q: ", x$q, "; ", sep = "")
- cat("PFER: ", x$PFER, "\n\n")
+ if (x$error.bound == "MB")
+ cat("PFER: ", x$PFER, "\n")
+ else
+ cat("PFER(*): ", x$PFER,
+ "\n (*) or expected number of low selection probability variables\n")
invisible(x)
}
@@ -208,7 +230,6 @@
### http://www.statslab.cam.ac.uk/~rds37/papers/r_concave_tail.R
### or
### http://www.statslab.cam.ac.uk/~rjs57/r_concave_tail.R
-
D <- function(theta, which, B, r) {
## If pi = ceil{ B * 2 * eta} / B + 1/B,..., 1 return the tail probability.
## If pi < ceil{ B * 2 * eta} / B return 1
@@ -220,7 +241,7 @@
s <- 1/r
thetaB <- theta * B
k_start <- (ceiling(2 * thetaB) + 1)
- if(k_start > B)
+ if(k_start >= B)
stop("theta to large")
Find.a <- function(prev_a)
@@ -254,14 +275,52 @@
return(max(cur_optim))
}
+## minD function for error bound in case of r-concavity
minD <- function(q, p, pi, B, r = c(-1/2, -1/4)) {
which <- ceiling(signif(pi / (1/(2* B)), 10))
- #if ((which) %% 1 != 0)
- # stop(sQuote("pi"), " must be a multiple of 1/(2 * B)")
-
- maxQ <- min(sqrt((B - 1) / (2 * B) * p^2),
- (B - 1) / (2 * B) * p)
+ maxQ <- maxQ(p, B)
if (q > maxQ)
stop(sQuote("q"), " must be <= ", maxQ)
min(c(1, D(q^2 / p^2, which - B, B, r[1]), D(q / p, which , 2*B, r[2])))
}
+
+## function to find optimal cutoff in stabsel (when error.bound = "SS")
+optimal_cutoff <- function(p, q, PFER, B) {
+ ## cutoff values can only be multiples of 1/(2B)
+ cutoff <- (2*B):1/(2*B)
+ cutoff <- cfs[cfs >= 0.5]
+ for (i in 1:length(cutoff)) {
+ if (minD(q, p, cutoff[i], B) * p > PFER) {
+ if (i == 1)
+ cutoff <- cutoff[i]
+ else
+ cutoff <- cutoff[i - 1]
+ break
+ }
+ }
+ cutoff[length(cutoff)]
+}
+
+## function to find optimal q in stabsel (when error.bound = "SS")
+optimal_q <- function(p, cutoff, PFER, B) {
+ for (q in 1:maxQ(p, B)) {
+ if (minD(q, p, cutoff, B) * p > PFER) {
+ q <- q - 1
+ break
+ }
+ }
+ max(1, q)
+}
+
+## obtain maximal value possible for q
+maxQ <- function(p, B) {
+ if(B <= 1)
+ stop("B must be at least 2")
+
+ fact_1 <- 4 * B / p
+ tmpfct <- function(q)
+ ceiling(q * fact_1) + 1 - 2 * B
+
+ res <- tmpfct(1:p)
+ length(res[res < 0])
+}
Modified: pkg/mboostDevel/R/methods.R
===================================================================
--- pkg/mboostDevel/R/methods.R 2013-10-30 17:56:05 UTC (rev 754)
+++ pkg/mboostDevel/R/methods.R 2013-11-06 14:31:09 UTC (rev 755)
@@ -427,6 +427,9 @@
selected.mboost <- function(object, ...)
object$xselect()
+selected.stabsel <- function(object, ...)
+ object$selected
+
summary.mboost <- function(object, ...) {
ret <- list(object = object, selprob = NULL)
Modified: pkg/mboostDevel/tests/regtest-inference.R
===================================================================
--- pkg/mboostDevel/tests/regtest-inference.R 2013-10-30 17:56:05 UTC (rev 754)
+++ pkg/mboostDevel/tests/regtest-inference.R 2013-11-06 14:31:09 UTC (rev 755)
@@ -131,3 +131,47 @@
dim(sbody$phat)
(sbody <- stabsel(mod, q = 3, PFER = 0.2, error.bound = "SS"))
dim(sbody$phat)
+
+
+## check stabsel_parameters and (theoretical) error control
+cutoff <- 0.6
+for (i in 1:10) {
+ print(stabsel_parameters(cutoff = cutoff, q = i, p = 100, error.bound = "MB"))
+}
+for (i in 1:10) {
+ print(stabsel_parameters(cutoff = cutoff, q = i, p = 100, error.bound = "SS"))
+}
+
+## check if missing values are determined correctly (especially at the extreme values)
+p <- 100
+B <- 50
+cutoff <- 0.6
+# low PFER
+PFER <- 0.001
+(res <- stabsel_parameters(p = p, cutoff = cutoff, PFER = PFER, B = B, error.bound = "SS"))
+stabsel_parameters(p = p, cutoff = cutoff, q = res$q, B = B, error.bound = "SS")
+# high PFER
+PFER <- 50
+(res <- stabsel_parameters(p = p, cutoff = cutoff, PFER = PFER, B = B, error.bound = "SS"))
+stabsel_parameters(p = p, cutoff = cutoff, q = res$q, B = B, error.bound = "SS")
+# medium PFER
+PFER <- 1
+(res <- stabsel_parameters(p = p, cutoff = cutoff, PFER = PFER, B = B, error.bound = "SS"))
+stabsel_parameters(p = p, cutoff = cutoff, q = res$q, B = B, error.bound = "SS")
+stabsel_parameters(p = p, cutoff = cutoff, q = res$q + 1, B = B, error.bound = "SS")
+
+
+q <- 10
+# high PFER
+PFER <- 5
+(res <- stabsel_parameters(p = p, q = q, PFER = PFER, B = B, error.bound = "SS"))
+stabsel_parameters(p = p, cutoff = res$cutoff, q = q, B = B, error.bound = "SS")
+# low PFER
+PFER <- 0.001
+(res <- stabsel_parameters(p = p, q = q, PFER = PFER, B = B, error.bound = "SS"))
+stabsel_parameters(p = p, cutoff = res$cutoff, q = q, B = B, error.bound = "SS")
+# medium PFER
+PFER <- 1
+(res <- stabsel_parameters(p = p, q = q, PFER = PFER, B = B, error.bound = "SS"))
+stabsel_parameters(p = p, cutoff = res$cutoff, q = q, B = B, error.bound = "SS")
+stabsel_parameters(p = p, cutoff = res$cutoff - 0.01, q = q, B = B, error.bound = "SS")
More information about the Mboost-commits
mailing list