[Mboost-commits] r756 - in pkg/mboostDevel: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 8 16:16:37 CET 2013
Author: hofner
Date: 2013-11-08 16:16:37 +0100 (Fri, 08 Nov 2013)
New Revision: 756
Modified:
pkg/mboostDevel/R/inference.R
pkg/mboostDevel/man/stabsel.Rd
pkg/mboostDevel/tests/regtest-inference.R
Log:
- bugfix in minD function needed for improved error bounds
- minor bugfixes in stabsel related functions
- added stabsel_parameters to manual
Modified: pkg/mboostDevel/R/inference.R
===================================================================
--- pkg/mboostDevel/R/inference.R 2013-11-06 14:31:09 UTC (rev 755)
+++ pkg/mboostDevel/R/inference.R 2013-11-08 15:16:37 UTC (rev 756)
@@ -5,6 +5,7 @@
papply = mclapply, verbose = TRUE, FWER,
error.bound = c("MB", "SS"), ...) {
+ call <- match.call()
p <- length(variable.names(object))
ibase <- 1:p
@@ -62,7 +63,8 @@
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)
+ max = mm, cutoff = cutoff, q = q, PFER = PFER, error.bound = error.bound,
+ call = call)
class(ret) <- "stabsel"
ret
}
@@ -118,7 +120,7 @@
upperbound <- q^2 / p / (2 * cutoff - 1)
} else {
cutoff <- optimal_cutoff(p, q, PFER, B)
- upperbound <- minD(q, p, cutoff, B) * p
+ upperbound <- tmp <- minD(q, p, cutoff, B) * p
}
upperbound <- signif(upperbound, 3)
if (verbose && tmp > 0.9 && upperbound - PFER > PFER/2) {
@@ -231,21 +233,22 @@
### 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
+ ## compute upper tail of r-concave distribution function
+ ## If q = ceil{ B * 2 * theta} / B + 1/B,..., 1 return the tail probability.
+ ## If q < ceil{ B * 2 * theta} / B return 1
- if (which <= 0)
- return(1)
- ### pi muss ein vielfaches von 1/(2 * B) sein, oder?
-
s <- 1/r
thetaB <- theta * B
k_start <- (ceiling(2 * thetaB) + 1)
- if(k_start >= B)
+
+ if (which < k_start)
+ return(1)
+
+ if(k_start > B)
stop("theta to large")
Find.a <- function(prev_a)
- uniroot(Calc.a, lower = 0.0001, upper = prev_a,
+ uniroot(Calc.a, lower = 0.00001, upper = prev_a,
tol = .Machine$double.eps^0.75)$root
Calc.a <- function(a) {
@@ -254,7 +257,7 @@
num / denom - thetaB
}
- OptimInt <- function(a) {
+ OptimInt <- function(a, t, k, thetaB, s) {
num <- (k + 1 - thetaB) * sum((a + 0:(t-1))^s)
denom <- sum((k + 1 - (0:k)) * (a + 0:k)^s)
1 - num / denom
@@ -267,16 +270,19 @@
for(k in k_start:B)
a_vec[k] <- Find.a(a_vec[k-1])
- t <- which
cur_optim <- rep(0, B)
for (k in k_start:(B-1))
cur_optim[k] <- optimize(f=OptimInt, lower = a_vec[k+1],
- upper = a_vec[k], maximum = TRUE)$objective
+ upper = a_vec[k],
+ t = which, k = k, thetaB = thetaB, s = s,
+ maximum = TRUE)$objective
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)) {
+ ## get the integer valued multiplier W of
+ ## pi = W * 1/(2 * B)
which <- ceiling(signif(pi / (1/(2* B)), 10))
maxQ <- maxQ(p, B)
if (q > maxQ)
@@ -288,7 +294,7 @@
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]
+ cutoff <- cutoff[cutoff >= 0.5]
for (i in 1:length(cutoff)) {
if (minD(q, p, cutoff[i], B) * p > PFER) {
if (i == 1)
Modified: pkg/mboostDevel/man/stabsel.Rd
===================================================================
--- pkg/mboostDevel/man/stabsel.Rd 2013-11-06 14:31:09 UTC (rev 755)
+++ pkg/mboostDevel/man/stabsel.Rd 2013-11-08 15:16:37 UTC (rev 756)
@@ -1,5 +1,6 @@
\name{stabsel}
\alias{stabsel}
+\alias{stabsel_parameters}
\title{
Stability Selection
}
@@ -12,6 +13,13 @@
B = ifelse(error.bound == "MB", 100, 50)),
papply = mclapply, verbose = TRUE, FWER,
error.bound = c("MB", "SS"), ...)
+
+## function to compute missing parameter from the other two parameters
+## (internally used within stabsel)
+stabsel_parameters(cutoff, q, PFER, p,
+ B = ifelse(error.bound == "MB", 100, 50),
+ verbose = FALSE, error.bound = c("MB", "SS"),
+ FWER)
}
\arguments{
\item{object}{an \code{mboost} object.}
@@ -22,7 +30,14 @@
specifies the amount of falsely selected base-learners, which is
tolerated. See details.}
\item{folds}{ a weight matrix with number of rows equal to the number
- of observations, see \code{\link{cvrisk}}.}
+ of observations, see \code{\link{cvrisk}}.}
+ \item{B}{ number of subsampling replicates. Per default, this is 100
+ for the error bound derived in Meinshausen & Buehlmann (2010) and
+ 50 for the error bound of Shah & Samworth (2013). In the latter
+ case, complementray pairs are used, thus leading to \eqn{2B}
+ subsamples.}
+ \item{p}{ number of possible predictors (including intercept if
+ applicable) }.
\item{papply}{ (parallel) apply function, defaults to \code{\link[parallel]{mclapply}}.
Alternatively, \code{parLapply} can be used. In the
latter case, usually more setup is needed (see example for some
@@ -31,10 +46,8 @@
\code{warnings} should be issued. }
\item{FWER}{ deprecated. Only for compatibility with older versions,
use PFER instead.}
- \item{error.bound}{
- use error bound of Meinshausen & Buehlmann (2010) (\dQuote{"MB"}) or
- of Shah & Samworth (2013) (\dQuote{"SS"}).
- }
+ \item{error.bound}{ use error bound of Meinshausen & Buehlmann (2010)
+ ("MB") or of Shah & Samworth (2013) ("SS"). }
\item{\dots}{additional arguments to \code{\link{cvrisk}}.}
}
\details{
@@ -80,8 +93,15 @@
### low-dimensional example
mod <- glmboost(DEXfat ~ ., data = bodyfat)
- (sbody <- stabsel(mod, q = 3, PFER = 1,
- folds = cv(model.weights(mod), type = "subsampling", B = 100)))
+
+ ## compute cutoff ahead of running stabsel to see if it is a sensible
+ ## parameter choice.
+ ## p = ncol(bodyfat) - 1 (= Outcome) + 1 ( = Intercept)
+ stabsel_parameters(q = 3, PFER = 1, p = ncol(bodyfat) - 1 + 1)
+
+ ## now run stability selection; to make results reproducible
+ set.seed(1234)
+ (sbody <- stabsel(mod, q = 3, PFER = 1))
opar <- par(mai = par("mai") * c(1, 1, 1, 2.7))
plot(sbody)
par(opar)
Modified: pkg/mboostDevel/tests/regtest-inference.R
===================================================================
--- pkg/mboostDevel/tests/regtest-inference.R 2013-11-06 14:31:09 UTC (rev 755)
+++ pkg/mboostDevel/tests/regtest-inference.R 2013-11-08 15:16:37 UTC (rev 756)
@@ -15,13 +15,14 @@
## If pi < ceil{ B * 2 * eta} / B return 1
MAXa <- 100000
- MINa <- 0.0001
+ MINa <- 0.00001
s <- -1/r
etaB <- eta * B
k_start <- (ceiling(2 * etaB) + 1)
- if(k_start > B)
- stop("eta is too large")
+ output <- rep(1, B)
+ if (k_start > B)
+ return(output)
a_vec <- rep(MAXa,B)
@@ -38,20 +39,22 @@
for(k in k_start:B)
a_vec[k] <- Find.a(a_vec[k-1])
+ # NB this function makes use of several gloabl variables
OptimInt <- function(a) {
num <- (k + 1 - etaB) * sum((a + 0:(t-1))^(-s))
denom <- sum((k + 1 - (0:k)) * (a + 0:k)^(-s))
1 - num / denom
}
- output <- rep(1, B)
-
prev_k <- k_start
for(t in k_start:B) {
cur_optim <- rep(0, B)
- for (k in prev_k:(B-1))
- cur_optim[k] <- optimize(f=OptimInt, lower = a_vec[k+1],
- upper = a_vec[k], maximum = TRUE)$objective
+ cur_optim[B] <- OptimInt(a_vec[B])
+ if (prev_k <= (B-1)) {
+ for (k in prev_k:(B-1))
+ cur_optim[k] <- optimize(f=OptimInt, lower = a_vec[k+1],
+ upper = a_vec[k], maximum = TRUE)$objective
+ }
output[t] <- max(cur_optim)
prev_k <- which.max(cur_optim)
}
@@ -92,21 +95,35 @@
points((40:100)/100, bound, col = "green")
stopifnot(all((bound - bound_ss) < sqrt(.Machine$double.eps)))
+## test r-concave bound
+B <- 50
+x <- (1:(2 * B))/(2 * B)
+p <- 1000
+q <- 490
+theta <- q/p
+
+## r-concave bound of Shah & Samworth (2013)
+bound_ss <- (pminD(theta, B) * p)[40:100]
+plot(x[40:100], bound_ss, xlab = "pi")
+## Bound of Meinshausen & Buehlmann (2010)
+points(x[40:100], q^2 / (2 * x[40:100] - 1) / p, col = "red")
+## now our implementation
+bound <- rep(NA, 61)
+for (i in 40:100) {
+ bound[i - 39] <- minD(q, p, i/100, B) * p
+}
+points((40:100)/100, bound, col = "green")
+stopifnot(all((bound - bound_ss) < sqrt(.Machine$double.eps)))
+
### computation of q from other values
cutoff <- 0.6
PFER <- 0.2
B <- 50
p <- 200
-objective <- function(q) {
- PFER / p - minD(q, p, cutoff, B)
-}
-root <- uniroot(objective, lower = 1,
- upper = min(sqrt((B - 1) / (2 * B) * p^2),
- (B - 1) / (2 * B) * p))$root
-(q <- ceiling(root))
+(q <- optimal_q(p = p, cutoff = cutoff, PFER = PFER, B = B))
# check:
-round(minD(q - 1, p, cutoff, B) * p, 3)
round(minD(q, p, cutoff, B) * p, 3)
+round(minD(q + 1, p, cutoff, B) * p, 3)
### computation of cutoff from other values
@@ -114,14 +131,10 @@
B <- 50
p <- 200
q <- 7
-objective <- function(cutoff) {
- PFER / p - minD(q, p, cutoff, B)
-}
-root <- uniroot(objective, lower = 0.5, upper = 0.9)$root
-(cutoff <- floor(root * 2 * B) / (2* B))
+(cutoff <- optimal_cutoff(p = p, q = q, PFER = PFER, B = B))
# check:
round(minD(q, p, cutoff, B) * p, 3)
-round(minD(q, p, cutoff + 1e-5, B) * p, 3)
+round(minD(q, p, cutoff - 1e-2, B) * p, 3)
### check stabsel interface
More information about the Mboost-commits
mailing list