[Mboost-commits] r715 - in pkg/mboostPatch: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 27 11:48:24 CEST 2013
Author: hofner
Date: 2013-06-27 11:48:24 +0200 (Thu, 27 Jun 2013)
New Revision: 715
Modified:
pkg/mboostPatch/R/family.R
pkg/mboostPatch/R/inference.R
pkg/mboostPatch/inst/CHANGES
pkg/mboostPatch/man/cvrisk.Rd
pkg/mboostPatch/man/stabsel.Rd
Log:
o stabsel was recoded and now uses different terminology and a better
tested code base
o fixed bugs in survival families:
- offset in all survival families was based on max(survtime) instead
of max(log(survtime));
- offset in CoxPH can't be computed from Cox Partial LH as constants
are canceled out; Use fixed offset instead;
Modified: pkg/mboostPatch/R/family.R
===================================================================
--- pkg/mboostPatch/R/family.R 2013-04-23 12:34:23 UTC (rev 714)
+++ pkg/mboostPatch/R/family.R 2013-06-27 09:48:24 UTC (rev 715)
@@ -281,9 +281,9 @@
.Call("ngradientCoxPLik", time, event, f, w, package = "mboost")
},
risk = risk <- function(y, f, w = 1) -sum(plloss(y, f, w), na.rm = TRUE),
- offset = function(y, w)
- optimize(risk, interval = c(0, max(y[,1], na.rm = TRUE)),
- y = y, w = w)$minimum,
+ offset = function(y, w = 1) 0, ## perhaps use something different
+ ## Note: offset cannot be computed from Cox Partial LH as
+ ## PLH doesn't depend on constant
check_y = function(y) {
if (!inherits(y, "Surv"))
stop("response is not an object of class ", sQuote("Surv"),
@@ -460,7 +460,7 @@
Family(ngradient = ngradient, risk = risk,
offset = function(y, w)
- optimize(risk, interval = c(0, max(y[,1], na.rm = TRUE)),
+ optimize(risk, interval = c(0, max(log(y[,1]), na.rm = TRUE)),
y = y, w = w)$minimum,
check_y = function(y) {
if (!inherits(y, "Surv"))
@@ -513,7 +513,7 @@
Family(ngradient = ngradient, risk = risk,
offset = function(y, w)
- optimize(risk, interval = c(0, max(y[,1], na.rm = TRUE)),
+ optimize(risk, interval = c(0, max(log(y[,1]), na.rm = TRUE)),
y = y, w = w)$minimum,
check_y = function(y) {
if (!inherits(y, "Surv"))
@@ -564,7 +564,7 @@
Family(ngradient = ngradient, risk = risk,
offset = function(y, w)
- optimize(risk, interval = c(0, max(y[,1], na.rm = TRUE)),
+ optimize(risk, interval = c(0, max(log(y[,1]), na.rm = TRUE)),
y = y, w = w)$minimum,
check_y = function(y) {
if (!inherits(y, "Surv"))
@@ -816,7 +816,7 @@
weights = "case",
offset = function(y,w){
optimize(risk,
- interval = c(0, max(y[,1], na.rm=TRUE)), y = y, w = w)$minimum
+ interval = c(0, max(log(y[,1]), na.rm=TRUE)), y = y, w = w)$minimum
},
check_y = function(y) {
if (!inherits(y,"Surv"))
Modified: pkg/mboostPatch/R/inference.R
===================================================================
--- pkg/mboostPatch/R/inference.R 2013-04-23 12:34:23 UTC (rev 714)
+++ pkg/mboostPatch/R/inference.R 2013-06-27 09:48:24 UTC (rev 715)
@@ -1,42 +1,74 @@
-stabsel <- function(object, FWER = 0.05, cutoff, q,
+stabsel <- function(object, cutoff, q, PFER,
folds = cv(model.weights(object), type = "subsampling", B = 100),
- papply = mclapply, verbose = TRUE, ...) {
+ papply = mclapply, verbose = TRUE, FWER, ...) {
p <- length(variable.names(object))
ibase <- 1:p
- if (!missing(q) && 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"))
+ ## 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 (!(FWER > 0 && FWER < 0.5))
- stop(sQuote("FWER"), " must be between 0 and 0.5")
+ 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 (! xor(missing(cutoff), missing(q)))
- stop(" Either ", sQuote("cutoff"), " or ", sQuote("q"),
- "must be specified (but not both).")
+ 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)) {
- cutoff <- min(0.9, tmp <- (q^2 / (FWER * p) + 1) / 2)
+ cutoff <- min(0.9, tmp <- (q^2 / (PFER * p) + 1) / 2)
upperbound <- q^2 / p / (2 * cutoff - 1)
- if (verbose && tmp > 0.9 && upperbound - FWER > FWER/2) {
- warning("Upper bound for FWER >> ", FWER,
+ 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 = ", min(1, round(upperbound, 2)), ")")
+ " (true upper bound = ", round(upperbound, 2), ")")
}
}
- if (missing(q)){
- stopifnot(cutoff >= 0.5)
- q <- ceiling(sqrt(FWER * (2 * cutoff - 1) * p))
+
+ if (missing(q)) {
+ q <- ceiling(sqrt(PFER * (2 * cutoff - 1) * p))
upperbound <- q^2 / p / (2 * cutoff - 1)
- if (verbose && upperbound - FWER > FWER/2)
- warning("Upper bound for FWER >> ", FWER,
+ 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)) {
+ upperbound <- PFER <- q^2 / p / (2 * cutoff - 1)
+ }
+ 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])))
@@ -77,12 +109,12 @@
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)
+ max = mm, cutoff = cutoff, q = q, PFER = upperbound)
class(ret) <- "stabsel"
ret
}
-print.stabsel <- function(x, ...) {
+print.stabsel <- function(x, decreasing = FALSE, ...) {
cat("\tStability Selection\n")
if (length(x$selected) > 0) {
@@ -92,9 +124,10 @@
cat("\nNo base-learner selected\n")
}
cat("\nSelection probabilities:\n")
- print(x$max[x$max > 0])
+ print(sort(x$max[x$max > 0], decreasing = decreasing))
cat("\nCutoff: ", x$cutoff, "; ", sep = "")
- cat("q: ", x$q, "\n\n")
+ cat("q: ", x$q, "; ", sep = "")
+ cat("PFER: ", x$PFER, "\n\n")
invisible(x)
}
Modified: pkg/mboostPatch/inst/CHANGES
===================================================================
--- pkg/mboostPatch/inst/CHANGES 2013-04-23 12:34:23 UTC (rev 714)
+++ pkg/mboostPatch/inst/CHANGES 2013-06-27 09:48:24 UTC (rev 715)
@@ -1,6 +1,15 @@
- CHANGES in `mboost' VERSION 2.2-2 (2013-XX-XX, rYYY)
+ CHANGES in `mboost' VERSION 2.2-3 (2013-XX-XX, rYYY)
- o speed up checking of manual by changing some computaions (e.g. reduce
+ o stabsel was recoded and now uses different terminology and a better
+ tested code base
+
+ o fixed bugs in survival families:
+ - offset in all survival families was based on max(survtime) instead
+ of max(log(survtime));
+ - offset in CoxPH can't be computed from Cox Partial LH as constants
+ are canceled out; Use fixed offset instead;
+
+ o speed up checking of manual by changing some computations (e.g. reduce
mstop) or exclude code from checking via \dontrun{}
o small improvements in manual
@@ -25,7 +34,7 @@
o fixed bug in brandom: now really use contrasts.arg = "contr.dummy" per
default.
- o removed tests/ folder and .Rout.save files for vignettes from the CRAN
+ o removed tests/ folder and .Rout.save files for vignettes from the CRAN
release
o small improvements in manual
Modified: pkg/mboostPatch/man/cvrisk.Rd
===================================================================
--- pkg/mboostPatch/man/cvrisk.Rd 2013-04-23 12:34:23 UTC (rev 714)
+++ pkg/mboostPatch/man/cvrisk.Rd 2013-06-27 09:48:24 UTC (rev 715)
@@ -22,9 +22,13 @@
using function \code{cv} and defaults to 25 bootstrap samples.}
\item{grid}{ a vector of stopping parameters the empirical risk
is to be evaluated for. }
- \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 details).}
+ \item{papply}{
+ (parallel) apply function, defaults to \code{\link[parallel]{mclapply}}.
+ Alternatively, \code{\link[parallel]{parLapply}} can be used. In the
+ latter case, usually more setup is needed (see example for some
+ details). To run \code{cvrisk} sequentially (i.e. not in parallel),
+ one can use \code{\link{lapply}}.
+ }
\item{fun}{ if \code{fun} is NULL, the out-of-sample risk is returned. \code{fun},
as a function of \code{object}, may extract any other characteristic
of the cross-validated models. These are returned as is.}
@@ -37,8 +41,7 @@
\item{prob}{ percentage of observations to be included in the learning samples
for subsampling.}
\item{strata}{ a factor of the same length as \code{weights} for stratification.}
- \item{...}{additional arguments passed to \code{\link[parallel]{mclapply}}
- eventually.}
+ \item{...}{additional arguments passed to \code{\link[parallel]{mclapply}}.}
}
\details{
Modified: pkg/mboostPatch/man/stabsel.Rd
===================================================================
--- pkg/mboostPatch/man/stabsel.Rd 2013-04-23 12:34:23 UTC (rev 714)
+++ pkg/mboostPatch/man/stabsel.Rd 2013-06-27 09:48:24 UTC (rev 715)
@@ -7,15 +7,18 @@
Selection of influential variables or model components with error control.
}
\usage{
-stabsel(object, FWER = 0.05, cutoff, q,
+stabsel(object, cutoff, q, PFER,
folds = cv(model.weights(object), type = "subsampling", B = 100),
- papply = mclapply, verbose = TRUE, ...)
+ papply = mclapply, verbose = TRUE, FWER, ...)
}
\arguments{
\item{object}{an \code{mboost} object.}
- \item{FWER}{family-wise error rate to be controlled by the selection procedure.}
- \item{cutoff}{cutoff between 0.5 and 1.}
- \item{q}{average number of selected base-learners.}
+ \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{papply}{ (parallel) apply function, defaults to \code{\link[parallel]{mclapply}}.
@@ -24,6 +27,8 @@
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{\dots}{additional arguments to \code{\link{cvrisk}}.}
}
\details{
@@ -31,10 +36,15 @@
This function implements the "stability selection" procedure
by Meinshausen and Buehlmann (2010).
- Either \code{cutoff} or \code{q} must be specified. The probability
- of selecting at least one non-influential variable (or model component)
- is less than \code{FWER}.
+ Two of the three arguments \code{cutoff}, \code{q} and \code{PFER}
+ \emph{must} be specified. The expected number of false positives E(V), where
+ V is the number of false positives, is controlled by \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 elements
@@ -43,7 +53,9 @@
\item{max }{maximum of selection probabilities.}
\item{cutoff }{cutoff used.}
\item{q }{average number of selected variables used.}
- \item{FWER }{family-wise error rate.}
+ \item{PFER }{per-family error rate.}
+
+ A special print method for objects of class exists.
}
\references{
@@ -56,10 +68,10 @@
data(bodyfat)
- ### (too) low-dimensional example
+ ### low-dimensional example
mod <- glmboost(DEXfat ~ ., data = bodyfat)
- (sbody <- stabsel(mod, q = 3,
- folds = cv(model.weights(mod), type = "subsampling", B = 25)))
+ (sbody <- stabsel(mod, q = 3, PFER = 1,
+ folds = cv(model.weights(mod), type = "subsampling", B = 100)))
opar <- par(mai = par("mai") * c(1, 1, 1, 2.7))
plot(sbody)
par(opar)
More information about the Mboost-commits
mailing list