[Mboost-commits] r791 - in pkg: mboostDevel mboostDevel/R mboostDevel/inst mboostDevel/man mboostDevel/tests mboostPatch mboostPatch/R mboostPatch/inst mboostPatch/man mboostPatch/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 18 17:36:37 CEST 2014
Author: hofner
Date: 2014-09-18 17:36:36 +0200 (Thu, 18 Sep 2014)
New Revision: 791
Removed:
pkg/mboostDevel/R/inference.R
pkg/mboostDevel/man/stabsel.Rd
pkg/mboostPatch/R/inference.R
pkg/mboostPatch/man/stabsel.Rd
pkg/mboostPatch/tests/regtest-inference.R
pkg/mboostPatch/tests/regtest-inference.Rout.save
Modified:
pkg/mboostDevel/DESCRIPTION
pkg/mboostDevel/NAMESPACE
pkg/mboostDevel/R/methods.R
pkg/mboostDevel/inst/CHANGES
pkg/mboostDevel/tests/regtest-inference.R
pkg/mboostDevel/tests/regtest-inference.Rout.save
pkg/mboostPatch/DESCRIPTION
pkg/mboostPatch/NAMESPACE
pkg/mboostPatch/inst/CHANGES
Log:
- stability selection was moved to a separate dedicated package "stabs"
Modified: pkg/mboostDevel/DESCRIPTION
===================================================================
--- pkg/mboostDevel/DESCRIPTION 2014-08-20 14:56:51 UTC (rev 790)
+++ pkg/mboostDevel/DESCRIPTION 2014-09-18 15:36:36 UTC (rev 791)
@@ -18,7 +18,7 @@
Depends: R (>= 2.14.0), methods, stats, parallel
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)
+ RColorBrewer, rpart (>= 4.0-3), stabs
LazyData: yes
License: GPL-2
URL: http://r-forge.r-project.org/projects/mboost/
Modified: pkg/mboostDevel/NAMESPACE
===================================================================
--- pkg/mboostDevel/NAMESPACE 2014-08-20 14:56:51 UTC (rev 790)
+++ pkg/mboostDevel/NAMESPACE 2014-09-18 15:36:36 UTC (rev 791)
@@ -18,7 +18,7 @@
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, stabsel_parameters,
+ cvrisk, cv, bbs,
bols, bspatial, brandom, btree, bss, bns, brad, bmono, bmrf, buser, survFit, selected,
nuisance, "%+%", "%X%", "%O%", extract, risk, "mstop<-")
###, basesel, fitsel)
@@ -70,10 +70,6 @@
S3method(cvrisk, mboost)
# 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)
Deleted: pkg/mboostDevel/R/inference.R
===================================================================
--- pkg/mboostDevel/R/inference.R 2014-08-20 14:56:51 UTC (rev 790)
+++ pkg/mboostDevel/R/inference.R 2014-09-18 15:36:36 UTC (rev 791)
@@ -1,457 +0,0 @@
-
-stabsel <- function(object, cutoff, q, PFER,
- folds = cv(model.weights(object), type = "subsampling",
- 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, ...) {
-
- call <- match.call()
- p <- length(variable.names(object))
- 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)
- ## 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(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,
- sampling.type = sampling.type, assumption = assumption,
- call = call)
- class(ret) <- "stabsel"
- ret
-}
-
-stabsel_parameters <- function(p, cutoff, q, PFER,
- B = ifelse(sampling.type == "MB", 100, 50),
- assumption = c("unimodal", "r-concave", "none"),
- sampling.type = c("SS", "MB"),
- verbose = FALSE, FWER) {
-
- sampling.type <- match.arg(sampling.type)
- if (sampling.type == "MB")
- assumption <- "none"
- else
- assumption <- match.arg(assumption)
-
-
- ## only two of the four arguments can be specified
- if ((nmiss <- sum(missing(PFER), missing(cutoff),
- missing(q), missing(FWER))) != 2) {
- if (nmiss > 2)
- stop("Two of the three argumnets ",
- sQuote("PFER"), ", ", sQuote("cutoff"), " and ", sQuote("q"),
- " must be specifed")
- if (nmiss < 2)
- stop("Only two of the three argumnets ",
- sQuote("PFER"), ", ", sQuote("cutoff"), " and ", sQuote("q"),
- " can be specifed at the same time")
- }
-
- if (!missing(FWER)) {
- if (!missing(PFER))
- stop(sQuote("FWER"), " and ", sQuote("PFER"),
- " cannot be spefified at the same time")
- PFER <- FWER
- warning(sQuote("FWER"), " is deprecated. Use ", sQuote("PFER"),
- " instead.")
- }
-
- if ((!missing(PFER) || !missing(FWER)) && PFER < 0)
- stop(sQuote("PFER"), " must be greater 0")
-
- if (!missing(cutoff) && (cutoff < 0.5 | cutoff > 1))
- stop(sQuote("cutoff"), " must be between 0.5 and 1")
-
- if (!missing(q)) {
- if (p < q)
- stop("Average number of selected base-learners ", sQuote("q"),
- " must be smaller \n than the number of base-learners",
- " specified in the model ", sQuote("object"))
- if (q < 0)
- stop("Average number of selected base-learners ", sQuote("q"),
- " must be greater 0")
- }
-
- if (missing(cutoff)) {
- if (assumption == "none") {
- cutoff <- min(1, tmp <- (q^2 / (PFER * p) + 1) / 2)
- upperbound <- q^2 / p / (2 * cutoff - 1)
- } else {
- if (assumption == "unimodal") {
- cutoff <- tmp <- optimal_cutoff(p, q, PFER, B,
- assumption = assumption)
- upperbound <- q^2 / p / um_const(cutoff, B, theta = q/p)
- } else {
- cutoff <- tmp <- optimal_cutoff(p, q, PFER, B,
- assumption = assumption)
- upperbound <- minD(q, p, cutoff, B) * p
- }
- }
- upperbound <- signif(upperbound, 3)
- if (verbose && tmp > 0.9 && upperbound - PFER > PFER/2) {
- warning("Upper bound for PFER > ", PFER,
- " for the given value of ", sQuote("q"),
- " (true upper bound = ", round(upperbound, 2), ")")
- }
- }
-
- if (missing(q)) {
- if (assumption == "none") {
- q <- floor(sqrt(PFER * (2 * cutoff - 1) * p))
- upperbound <- q^2 / p / (2 * cutoff - 1)
- } else {
- if (assumption == "unimodal") {
- q <- optimal_q(p, cutoff, PFER, B, assumption = assumption)
- upperbound <- q^2 / p / um_const(cutoff, B, theta = q/p)
- } else {
- q <- optimal_q(p, cutoff, PFER, B, assumption = assumption)
- upperbound <- minD(q, p, cutoff, B) * p
- }
- }
- upperbound <- signif(upperbound, 3)
- if (verbose && upperbound - PFER > PFER/2)
- warning("Upper bound for PFER > ", PFER,
- " for the given value of ", sQuote("cutoff"),
- " (true upper bound = ", upperbound, ")")
- }
-
- if (missing(PFER)) {
- if (assumption == "none") {
- upperbound <- PFER <- q^2 / p / (2 * cutoff - 1)
- } else {
- if (assumption == "unimodal") {
- upperbound <- PFER <- q^2 / p / um_const(cutoff, B, theta = q/p)
- } else {
- upperbound <- PFER <- minD(q, p, cutoff, B) * p
- }
- }
- upperbound <- signif(upperbound, 3)
- }
-
- if (verbose && PFER >= p)
- warning("Upper bound for PFER larger than the number of base-learners.")
-
- res <- list(cutoff = cutoff, q = q, PFER = upperbound,
- sampling.type = sampling.type, assumption = assumption)
- class(res) <- "stabsel_parameters"
- res
-}
-
-print.stabsel <- function(x, decreasing = FALSE, ...) {
-
- cat("\tStability Selection")
- if (x$assumption == "none")
- cat(" without further assumptions\n")
- if (x$assumption == "unimodal")
- cat(" with unimodality assumption\n")
- if (x$assumption == "r-concave")
- cat(" with r-concavity assumption\n")
- if (length(x$selected) > 0) {
- cat("\nSelected base-learners:\n")
- print(x$selected)
- } else {
- cat("\nNo base-learner selected\n")
- }
- cat("\nSelection probabilities:\n")
- print(sort(x$max[x$max > 0], decreasing = decreasing))
- cat("\n")
- print.stabsel_parameters(x, heading = FALSE)
- cat("\n")
- invisible(x)
-}
-
-print.stabsel_parameters <- function(x, heading = TRUE, ...) {
- if (heading) {
- cat("Stability Selection")
- if (x$assumption == "none")
- cat(" without further assumptions\n")
- if (x$assumption == "unimodal")
- cat(" with unimodality assumption\n")
- if (x$assumption == "r-concave")
- cat(" with r-concavity assumption\n")
- }
- cat("Cutoff: ", x$cutoff, "; ", sep = "")
- cat("q: ", x$q, "; ", sep = "")
- if (x$sampling.type == "MB")
- cat("PFER: ", x$PFER, "\n")
- else
- cat("PFER(*): ", x$PFER,
- "\n (*) or expected number of low selection probability variables\n")
- invisible(x)
-}
-
-plot.stabsel <- function(x, main = deparse(x$call), type = c("paths", "maxsel"),
- col = NULL, ymargin = 10, np = sum(x$max > 0),
- labels = NULL, ...) {
-
- type <- match.arg(type)
-
- if (is.null(col))
- col <- hcl(h = 40, l = 50, c = x$max / max(x$max) * 490)
-
- if (type == "paths") {
- ## if par(mar) not set by user ahead of plotting
- if (all(par()[["mar"]] == c(5, 4, 4, 2) + 0.1))
- ..old.par <- par(mar = c(5, 4, 4, ymargin) + 0.1)
- h <- x$phat
- h <- h[rowSums(h) > 0, , drop = FALSE]
- matplot(t(h), type = "l", lty = 1,
- xlab = "Number of boosting iterations",
- ylab = "Selection probability",
- main = main, col = col[x$max > 0], ylim = c(0, 1), ...)
- abline(h = x$cutoff, lty = 1, col = "lightgray")
- if (is.null(labels))
- labels <- rownames(x$phat)
- axis(4, at = x$phat[rowSums(x$phat) > 0, ncol(x$phat)],
- labels = labels[rowSums(x$phat) > 0], las = 1)
- } else {
- ## if par(mar) not set by user ahead of plotting
- if (all(par()[["mar"]] == c(5, 4, 4, 2) + 0.1))
- ..old.par <- par(mar = c(5, ymargin, 4, 2) + 0.1)
- if (np > length(x$max))
- stop(sQuote("np"), "is set too large")
- inc_freq <- x$max ## inclusion frequency
- plot(tail(sort(inc_freq), np), 1:np,
- type = "n", yaxt = "n", xlim = c(0, 1),
- ylab = "", xlab = expression(hat(pi)),
- main = main, ...)
- abline(h = 1:np, lty = "dotted", col = "grey")
- points(tail(sort(inc_freq), np), 1:np, pch = 19,
- col = col[tail(order(inc_freq), np)])
- if (is.null(labels))
- labels <- names(x$max)
- axis(2, at = 1:np, labels[tail(order(inc_freq), np)], las = 2)
- ## add cutoff
- abline(v = x$cutoff, col = "grey")
- }
- if (exists("..old.par"))
- par(..old.par) # reset plotting settings
-}
-
-
-
-fitsel <- function(object, newdata = NULL, which = NULL, ...) {
- fun <- function(model) {
- tmp <- predict(model, newdata = newdata,
- which = which, agg = "cumsum")
- ret <- c()
- for (i in 1:length(tmp))
- ret <- rbind(ret, tmp[[i]])
- ret
- }
- ss <- cvrisk(object, fun = fun, ...)
- ret <- matrix(0, nrow = nrow(ss[[1]]), ncol = ncol(ss[[1]]))
- for (i in 1:length(ss))
- ret <- ret + sign(ss[[i]])
- ret <- abs(ret) / length(ss)
- ret
-}
-
-
-### Modified version of the code accompanying the paper:
-### Shah, R. D. and Samworth, R. J. (2013), Variable selection with error
-### control: Another look at Stability Selection, J. Roy. Statist. Soc., Ser.
-### B, 75, 55-80. DOI: 10.1111/j.1467-9868.2011.01034.x
-###
-### Original code available from
-### 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) {
- ## 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
-
- s <- 1/r
- thetaB <- theta * B
- k_start <- (ceiling(2 * thetaB) + 1)
-
- if (which < k_start)
- return(1)
-
- if(k_start > B)
- stop("theta to large")
-
- Find.a <- function(prev_a)
- uniroot(Calc.a, lower = 0.00001, upper = prev_a,
- tol = .Machine$double.eps^0.75)$root
-
- Calc.a <- function(a) {
- denom <- sum((a + 0:k)^s)
- num <- sum((0:k) * (a + 0:k)^s)
- num / denom - thetaB
- }
-
- 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
- }
-
- ## initialize a
- a_vec <- rep(100000, B)
-
- ## compute a values
- for(k in k_start:B)
- a_vec[k] <- Find.a(a_vec[k-1])
-
- 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],
- 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)
- 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 sampling.type = "SS")
-optimal_cutoff <- function(p, q, PFER, B, assumption = "unimodal") {
- if (assumption == "unimodal") {
- ## cutoff values can only be multiples of 1/(2B)
- cutoffgrid <- 1/2 + (2:B)/(2*B)
- c_min <- min(0.5 + (q/p)^2, 0.5 + 1/(2*B) + 0.75 * (q/p)^2)
- cutoffgrid <- cutoffgrid[cutoffgrid > c_min]
- upperbound <- rep(NA, length(cutoffgrid))
- for (i in 1:length(cutoffgrid))
- upperbound[i] <- q^2 / p / um_const(cutoffgrid[i], B, theta = q/p)
- cutoff <- cutoffgrid[upperbound < PFER][1]
- return(cutoff)
- } else {
- ## cutoff values can only be multiples of 1/(2B)
- cutoff <- (2*B):1/(2*B)
- cutoff <- cutoff[cutoff >= 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
- }
- }
- return(tail(cutoff, 1))
- }
-}
-
-## function to find optimal q in stabsel (when sampling.type = "SS")
-optimal_q <- function(p, cutoff, PFER, B, assumption = "unimodal") {
- if (assumption == "unimodal") {
- if (cutoff <= 0.75) {
- upper_q <- max(p * sqrt(cutoff - 0.5),
- p * sqrt(4/3 * (cutoff - 0.5 - 1/(2*B))))
- ## q must be an integer < upper_q
- upper_q <- ceiling(upper_q - 1)
- } else {
- upper_q <- p
- }
- q <- uniroot(function(q)
- q^2 / p / um_const(cutoff, B, theta = q/p) - PFER,
- lower = 1, upper = upper_q)$root
- return(floor(q))
- } else {
- for (q in 1:maxQ(p, B)) {
- if (minD(q, p, cutoff, B) * p > PFER) {
- q <- q - 1
- break
- }
- }
- return(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])
-}
-
-## obtain constant for unimodal bound
-um_const <- function(cutoff, B, theta) {
- if (cutoff <= 3/4) {
- if (cutoff < 1/2 + min(theta^2, 1 / (2*B) + 3/4 * theta^2))
- stop ("cutoff out of bounds")
- return( 2 * (2 * cutoff - 1 - 1/(2*B)) )
- } else {
- if (cutoff > 1)
- stop ("cutoff out of bounds")
- return( (1 + 1/B)/(4 * (1 - cutoff + 1 / (2*B))) )
- }
-}
Modified: pkg/mboostDevel/R/methods.R
===================================================================
--- pkg/mboostDevel/R/methods.R 2014-08-20 14:56:51 UTC (rev 790)
+++ pkg/mboostDevel/R/methods.R 2014-09-18 15:36:36 UTC (rev 791)
@@ -451,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/mboostDevel/inst/CHANGES
===================================================================
--- pkg/mboostDevel/inst/CHANGES 2014-08-20 14:56:51 UTC (rev 790)
+++ pkg/mboostDevel/inst/CHANGES 2014-09-18 15:36:36 UTC (rev 791)
@@ -9,7 +9,9 @@
o tweaked update function: we now can turn the trace off and specify
the type of risk as well as the oobweight to update()
+ o stabsel has been moved to the new package stabs
+
CHANGES in `mboost' VERSION 2.3-1 (2014-xx-yy, rXYZ)
o changed vignette mboost_tutorial to reflect latest changes in mboost.
Deleted: pkg/mboostDevel/man/stabsel.Rd
===================================================================
--- pkg/mboostDevel/man/stabsel.Rd 2014-08-20 14:56:51 UTC (rev 790)
+++ pkg/mboostDevel/man/stabsel.Rd 2014-09-18 15:36:36 UTC (rev 791)
@@ -1,159 +0,0 @@
-\name{stabsel}
-\alias{stabsel}
-\alias{stabsel_parameters}
-\alias{stabsel_parameters.default}
-\alias{stabsel_parameters.mboost}
-\alias{plot.stabsel}
-\title{
- Stability Selection
-}
-\description{
- Selection of influential variables or model components with error control.
-}
-\usage{
-stabsel(object, cutoff, q, PFER,
- folds = cv(model.weights(object), type = "subsampling",
- 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, ...)
-
-## function to compute missing parameter from the other two parameters
-## (internally used within stabsel)
-stabsel_parameters(p, cutoff, q, PFER,
- B = ifelse(sampling.type == "MB", 100, 50),
- assumption = c("unimodal", "r-concave", "none"),
- sampling.type = c("SS", "MB"),
- verbose = FALSE, FWER)
-
-\method{plot}{stabsel}(x, main = deparse(x$call), type = c("paths", "maxsel"),
- col = NULL, ymargin = 10, np = sum(x$max > 0),
- labels = NULL, ...)
-}
-\arguments{
- \item{object}{an \code{mboost} object.}
- \item{cutoff}{cutoff between 0.5 and 1. Preferably a value between 0.6
- and 0.9 should be used.}
- \item{q}{number of (unique) selected base-learners per boosting run.}
- \item{PFER}{upper bound for the per-family error rate. This
- 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}}.}
- \item{assumption}{ Defines the type of assumptions on the
- distributions of the selection probabilities and simultaneous
- selection probabilities. Only applicable for
- \code{sampling.type = "SS"}. For \code{sampling.type = "MB"} we
- always use code{"none"}.}
- \item{sampling.type}{ use sampling scheme of of Shah & Samworth
- (2013), i.e., with complementarty pairs (\code{sampling.type = "SS"}),
- or the original sampling scheme of Meinshausen & Buehlmann (2010).}
- \item{p}{ number of possible predictors (including intercept if
- applicable).}
- \item{B}{ number of subsampling replicates. Per default, we use 50
- complementary pairs for the error bounds of Shah & Samworth (2013)
- and 100 for the error bound derived in Meinshausen & Buehlmann
- (2010). As we use \eqn{B} complementray pairs in the former case
- this leads to \eqn{2B} subsamples.}
- \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 of \code{\link{cvrisk}} for some details).}
- \item{verbose}{ logical (default: \code{TRUE}) that determines wether
- \code{warnings} should be issued. }
- \item{FWER}{ deprecated. Only for compatibility with older versions,
- use PFER instead.}
- \item{eval}{ logical. Determines whether stability selection is
- evaluated (\code{eval = TRUE}; default) or if only the parameter
- combination is returned.}
- \item{x}{object of class \code{stabsel}.}
- \item{main}{main title for the plot.}
- \item{type}{plot type; either stability paths (\code{"paths"}) or a
- plot of the maximum selection frequency (\code{"maxsel"}).}
- \item{col}{a vector of colors; Typically, one can specify a single
- color or one color for each variable. Per default, colors depend on
- the maximal selection frequency of the variable and range from grey
- to red.}
- \item{ymargin}{(temporarily) specifies the y margin of of the plot in
- lines (see argument \code{"mar"} of function \code{\link{par}}).
- This only affects the right margin for \code{type = "paths"} and
- the left margin for \code{type = "maxsel"}. Explicit user specified
- margins are kept and are not overwritten.}
- \item{np}{number of variables to plot for the maximum selection
- frequency plot (\code{type = "maxsel"}); the first \code{np}
- variables with highest selection frequency are plotted.}
- \item{labels}{variable labels for the plot; one label per base-learner
- must be specified. Per default, names of base-learners are used.}
- \item{\dots}{additional arguments to \code{\link{cvrisk}} and further
- arguments to parallel apply methods such as \code{\link{mclapply}}
- or additional arguments to plot functions.}
-
-}
-\details{
-
- This function implements the stability selection procedure
- by Meinshausen and Buehlmann (2010) and the improved error bounds by
- Shah and Samworth (2013).
-
- Two of the three arguments \code{cutoff}, \code{q} and \code{PFER}
- \emph{must} be specified. The per-family error rate (PFER), i.e., the
- expected number of false positives \eqn{E(V)}, where \eqn{V} is the
- number of false positives, is bounded by the argument \code{PFER}.
-
- As controlling the PFER is more conservative as controlling the
- family-wise error rate (FWER), the procedure also controlls the FWER,
- i.e., the probability of selecting at least one non-influential
- variable (or model component) is less than \code{PFER}.
-
-}
-\value{
- An object of class \code{stabsel} with a special \code{print} method.
- The object has the following elements:
- \item{phat}{selection probabilities.}
- \item{selected}{elements with maximal selection probability greater
- \code{cutoff}.}
- \item{max}{maximum of selection probabilities.}
- \item{cutoff}{cutoff used.}
- \item{q}{average number of selected variables used.}
- \item{PFER}{per-family error rate.}
- \item{sampling.type}{the sampling type used for stability selection.}
- \item{assumption}{the assumptions made on the selection
- probabilities.}
- \item{call}{the call.}
-}
-\references{
-
- N. Meinshausen and P. Buehlmann (2010), Stability selection.
- \emph{Journal of the Royal Statistical Society, Series B},
- \bold{72}, 417--473.
-
- R.D. Shah and R.J. Samworth (2013), Variable selection with error
- control: another look at stability selection. \emph{Journal of the Royal
- Statistical Society, Series B}, \bold{75}, 55--80.
-
-}
-\examples{
-
- data("bodyfat", package = "TH.data")
-
- ### low-dimensional example
- mod <- glmboost(DEXfat ~ ., data = bodyfat)
-
- ## 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,
- sampling.type = "MB")
- ## the same:
- stabsel(mod, q = 3, PFER = 1, sampling.type = "MB", eval = FALSE)
-
- ## now run stability selection; to make results reproducible
- set.seed(1234)
- (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{nonparametric}
Modified: pkg/mboostDevel/tests/regtest-inference.R
===================================================================
--- pkg/mboostDevel/tests/regtest-inference.R 2014-08-20 14:56:51 UTC (rev 790)
+++ pkg/mboostDevel/tests/regtest-inference.R 2014-09-18 15:36:36 UTC (rev 791)
@@ -3,225 +3,6 @@
set.seed(1907)
-### (Slightly) modified version of the code accompanying the paper:
-### Shah, R. D. and Samworth, R. J. (2013), Variable selection with error
-### control: Another look at Stability Selection, J. Roy. Statist. Soc., Ser.
-### B, 75, 55-80. DOI: 10.1111/j.1467-9868.2011.01034.x
-###
-### Original code available from
-### http://www.statslab.cam.ac.uk/~rds37/papers/r_concave_tail.R
-### or
-### http://www.statslab.cam.ac.uk/~rjs57/r_concave_tail.R
-r.TailProbs <- function(eta, 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
-
- MAXa <- 100000
- MINa <- 0.00001
-
- s <- -1/r
- etaB <- eta * B
- k_start <- (ceiling(2 * etaB) + 1)
- output <- rep(1, B)
- if (k_start > B)
- return(output)
-
- a_vec <- rep(MAXa,B)
-
- Find.a <- function(prev_a)
- uniroot(Calc.a, lower = MINa, upper = prev_a,
- tol = .Machine$double.eps^0.75)$root
-
- Calc.a <- function(a) {
- denom <- sum((a + 0:k)^(-s))
- num <- sum((0:k) * (a + 0:k)^(-s))
- num / denom - etaB
- }
-
- 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
- }
-
- prev_k <- k_start
- for(t in k_start:B) {
- cur_optim <- rep(0, B)
- 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)
- }
- return(output)
-}
-
-pminD <- function(theta, B, r = c(-1/2, -1/4)) {
- pmin(c(rep(1, B), r.TailProbs(theta^2, B, r[1])),
- r.TailProbs(theta, 2*B, r[2]))
-}
-
-## test r-concave bound
-B <- 50
-x <- (1:(2 * B))/(2 * B)
-p <- 1000
-q <- 50
-theta <- q/p
-if (FALSE) {
- plot(x, log(pminD(theta, B)), xlab = "pi")
- abline(v = ceiling(2 * theta * 100) / 100)
- Ds <- cbind(c(rep(1, B), r.TailProbs(theta^2, B, -1/2)),
- r.TailProbs(theta, 2*B, -1/4))
- round(log(Ds), 2)
- lines(x, log(Ds[,1]), col = "red", lwd = 2)
- lines(x, log(Ds[,2]), col = "blue", lwd = 2)
-}
-
-## r-concave bound of Shah & Samworth (2013)
-bound_ss <- (pminD(theta, B) * p)[40:100]
-plot(x[40:100], bound_ss, xlab = "pi", ylim = c(0, 50))
-## 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)))
-
-## 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
-(q <- optimal_q(p = p, cutoff = cutoff, PFER = PFER, B = B,
- assumption = "r-concave"))
-# check:
-(a <- round(minD(q, p, cutoff, B) * p, 3))
-(b <- round(minD(q + 1, p, cutoff, B) * p, 3))
-stopifnot(a < PFER && b > PFER)
-
-## same for unimodal bound
-(q <- optimal_q(p = p, cutoff = cutoff, PFER = PFER, B = B,
- assumption = "unimodal"))
-
-### computation of cutoff from other values
-PFER <- 0.2
-B <- 50
-p <- 200
-q <- 7
-(cutoff <- optimal_cutoff(p = p, q = q, PFER = PFER, B = B,
- assumption = "r-concave"))
-# check:
-(a <- round(minD(q, p, cutoff, B) * p, 3))
-(b <- round(minD(q, p, cutoff - 1e-2, B) * p, 3))
-stopifnot(a < PFER && b > PFER)
-
-## same for unimodal bound
-(cutoff <- optimal_cutoff(p = p, q = q, PFER = PFER, B = B,
- assumption = "unimodal"))
-
-### check stabsel interface
-data("bodyfat", package = "TH.data")
-mod <- glmboost(DEXfat ~ ., data = bodyfat)
-(sbody <- stabsel(mod, q = 3, PFER = 0.2, sampling.type = "MB"))
-dim(sbody$phat)
-(sbody <- stabsel(mod, q = 3, PFER = 0.2, sampling.type = "SS"))
-dim(sbody$phat)
-
-## check interface of stabsel_parameters
-stabsel(mod, q = 3, PFER = 0.2, sampling.type = "SS", eval = FALSE)
-
-
-## check stabsel_parameters and (theoretical) error control
-cutoff <- 0.6
-for (i in 1:10) {
- print(stabsel_parameters(cutoff = cutoff, q = i, p = 100, sampling.type = "MB"))
-}
-for (i in 1:10) {
- print(stabsel_parameters(cutoff = cutoff, q = i, p = 100, sampling.type = "SS",
- assumption = "unimodal"))
- print(stabsel_parameters(cutoff = cutoff, q = i, p = 100, sampling.type = "SS",
- assumption = "r-concave"))
-}
-
-## 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,
- sampling.type = "SS", assumption = "r-concave"))
-stabsel_parameters(p = p, cutoff = cutoff, q = res$q, B = B,
- sampling.type = "SS", assumption = "r-concave")
-# high PFER
-PFER <- 50
-(res <- stabsel_parameters(p = p, cutoff = cutoff, PFER = PFER, B = B,
- sampling.type = "SS", assumption = "r-concave"))
-stabsel_parameters(p = p, cutoff = cutoff, q = res$q, B = B,
- sampling.type = "SS", assumption = "r-concave")
-# medium PFER
-PFER <- 1
-(res <- stabsel_parameters(p = p, cutoff = cutoff, PFER = PFER, B = B,
- sampling.type = "SS", assumption = "r-concave"))
-stabsel_parameters(p = p, cutoff = cutoff, q = res$q, B = B,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/mboost -r 791
More information about the Mboost-commits
mailing list