[Mboost-commits] r753 - in pkg/mboostDevel: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 28 17:20:41 CET 2013
Author: hofner
Date: 2013-10-28 17:20:41 +0100 (Mon, 28 Oct 2013)
New Revision: 753
Modified:
pkg/mboostDevel/R/bmono.R
pkg/mboostDevel/R/helpers.R
pkg/mboostDevel/R/inference.R
pkg/mboostDevel/man/baselearners.Rd
pkg/mboostDevel/man/stabsel.Rd
pkg/mboostDevel/tests/regtest-inference.R
Log:
- bugfix in bmono with type = "quad.prog" when some directions are unpenalized
- implemented a less stringent error.bound [see Shah, R. D. and Samworth, R. J. (2013),
Variable selection with error control: Another look at Stability Selection, JRSSB]
- updated manuals
Modified: pkg/mboostDevel/R/bmono.R
===================================================================
--- pkg/mboostDevel/R/bmono.R 2013-10-25 12:33:29 UTC (rev 752)
+++ pkg/mboostDevel/R/bmono.R 2013-10-28 16:20:41 UTC (rev 753)
@@ -232,6 +232,8 @@
## set lambda2 = 0 if no constraint is used
if (any(none <- args$constraint == "none"))
lambda2[none] <- 0
+ if (any(none <- lambda2 == 0))
+ args$constraint[none] <- "none"
## <FIXME> Boundary constraints for bivariate smooths are currently not
## implemented
if (args$boundary.constraints)
@@ -297,8 +299,7 @@
}
} else { ## i.e. type == "quad.prog"
if (lambda2[[2]] == 0) {
- coef <- solveLSEI(XtX, crossprod(X, y),
- constraint = args$constraint)
+ coef <- solveLSEI(XtX, crossprod(X, y), D = D[[1]])
} else {
coef <- solveLSEI(XtX, crossprod(X, y), D = D)
}
Modified: pkg/mboostDevel/R/helpers.R
===================================================================
--- pkg/mboostDevel/R/helpers.R 2013-10-25 12:33:29 UTC (rev 752)
+++ pkg/mboostDevel/R/helpers.R 2013-10-28 16:20:41 UTC (rev 753)
@@ -224,26 +224,13 @@
## least squares with equalities and inequalities
-solveLSEI <- function(XtX, Xty, D = NULL, constraint = "none") {
- if (constraint == "none" && is.null(D))
+solveLSEI <- function(XtX, Xty, D = NULL) {
+ if (is.null(D) || all(sapply(D, is.null)))
return(solve(XtX, Xty, LINPACK = FALSE))
- if (is.null(D)) {
- ## i.e., if constraint != "none"
- D <- differences(constraint, ncol = ncol(XtX))
- }
-
if (!isMATRIX(D)) ## i.e. D is a list with 2 entries
D <- rbind(D[[1]], D[[2]])
## NOTE: Currently both constraints get the same weight
-
-## first we used package limSolve
-# cf <- lsei(A= XtX, B = Xty, G = D, H = rep(0, nrow(D)), type = 1,
-# E = matrix(0, nrow(XtX), nrow(XtX)), F = rep(0, nrow(XtX)),
-# tol = .Machine$double.eps,
-# tolrank = c(.Machine$double.eps, .Machine$double.eps),
-# fulloutput = TRUE)$X
-## but to reduce computational overhead we directly use quadprog
cf <- solve.QP(Dmat = XtX, dvec = as.vector(Xty), Amat = t(D),
bvec = rep(0, nrow(D)))$solution
cf
Modified: pkg/mboostDevel/R/inference.R
===================================================================
--- pkg/mboostDevel/R/inference.R 2013-10-25 12:33:29 UTC (rev 752)
+++ pkg/mboostDevel/R/inference.R 2013-10-28 16:20:41 UTC (rev 753)
@@ -1,14 +1,15 @@
stabsel <- function(object, cutoff, q, PFER,
folds = cv(model.weights(object), type = "subsampling",
- B = ifelse(error.bound = "MB", 100, 50)),
+ B = ifelse(error.bound == "MB", 100, 50)),
papply = mclapply, verbose = TRUE, FWER,
error.bound = c("MB", "SS"), ...) {
p <- length(variable.names(object))
ibase <- 1:p
- error.bound <- match.args(error.bound)
+ error.bound <- match.arg(error.bound)
+ B <- ncol(folds)
## only two of the four arguments can be specified
if ((nmiss <- sum(missing(PFER), missing(cutoff),
@@ -53,11 +54,11 @@
cutoff <- min(0.9, tmp <- (q^2 / (PFER * p) + 1) / 2)
upperbound <- q^2 / p / (2 * cutoff - 1)
} else {
- objective <- function(cutoff) {
+ objective_cf <- 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)
+ root <- uniroot(objective_cf, lower = 0.5, upper = 0.9)$root
+ cutoff <- min(0.9, tmp <- floor(root * 2 * B) / (2* B))
upperbound <- minD(q, p, cutoff, B) * p
}
upperbound <- signif(upperbound, 3)
@@ -73,10 +74,10 @@
q <- ceiling(sqrt(PFER * (2 * cutoff - 1) * p))
upperbound <- q^2 / p / (2 * cutoff - 1)
} else {
- objective <- function(q) {
+ objective_q <- function(q) {
PFER / p - minD(q, p, cutoff, B)
}
- root <- uniroot(objective, lower = 1,
+ root <- uniroot(objective_q, lower = 1,
upper = min(sqrt((B - 1) / (2 * B) * p^2),
(B - 1) / (2 * B) * p))$root
q <- ceiling(root)
@@ -112,7 +113,7 @@
}
ss <- cvrisk(object, fun = fun,
folds = folds,
- papply = papply, ....)
+ papply = papply, ...)
if (verbose){
qq <- sapply(ss, function(x) length(unique(x)))
Modified: pkg/mboostDevel/man/baselearners.Rd
===================================================================
--- pkg/mboostDevel/man/baselearners.Rd 2013-10-25 12:33:29 UTC (rev 752)
+++ pkg/mboostDevel/man/baselearners.Rd 2013-10-28 16:20:41 UTC (rev 753)
@@ -50,7 +50,7 @@
bmono(...,
constraint = c("increasing", "decreasing", "convex", "concave",
"none", "positive", "negative"),
- type = c("iterative", "quad.prog"),
+ type = c("quad.prog", "iterative"),
by = NULL, index = NULL, knots = 20, boundary.knots = NULL,
degree = 3, differences = 2, df = 4, lambda = NULL,
lambda2 = 1e6, niter=10, intercept = TRUE,
Modified: pkg/mboostDevel/man/stabsel.Rd
===================================================================
--- pkg/mboostDevel/man/stabsel.Rd 2013-10-25 12:33:29 UTC (rev 752)
+++ pkg/mboostDevel/man/stabsel.Rd 2013-10-28 16:20:41 UTC (rev 753)
@@ -8,8 +8,10 @@
}
\usage{
stabsel(object, cutoff, q, PFER,
- folds = cv(model.weights(object), type = "subsampling", B = 100),
- papply = mclapply, verbose = TRUE, FWER, ...)
+ folds = cv(model.weights(object), type = "subsampling",
+ B = ifelse(error.bound == "MB", 100, 50)),
+ papply = mclapply, verbose = TRUE, FWER,
+ error.bound = c("MB", "SS"), ...)
}
\arguments{
\item{object}{an \code{mboost} object.}
@@ -29,6 +31,10 @@
\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{\dots}{additional arguments to \code{\link{cvrisk}}.}
}
\details{
@@ -61,8 +67,12 @@
N. Meinshausen and P. Buehlmann (2010), Stability selection.
\emph{Journal of the Royal Statistical Society, Series B},
- \bold{72}(4).
+ \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{
Modified: pkg/mboostDevel/tests/regtest-inference.R
===================================================================
--- pkg/mboostDevel/tests/regtest-inference.R 2013-10-25 12:33:29 UTC (rev 752)
+++ pkg/mboostDevel/tests/regtest-inference.R 2013-10-28 16:20:41 UTC (rev 753)
@@ -90,7 +90,7 @@
bound[i - 39] <- minD(q, p, i/100, B) * p
}
points((40:100)/100, bound, col = "green")
-stopifnot(bound == bound_ss)
+stopifnot(all((bound - bound_ss) < sqrt(.Machine$double.eps)))
### computation of q from other values
cutoff <- 0.6
@@ -124,8 +124,10 @@
round(minD(q, p, cutoff + 1e-5, B) * p, 3)
-weights <- c(2, 3, 2, rep(1, 90), rep(0, 10))
-folds <- cv(weights, type = "subsampling", B = 10)
-head(folds)
-tail(folds)
-rowSums(cbind(folds, weights - folds))
+### check stabsel interface
+data("bodyfat", package = "TH.data")
+mod <- glmboost(DEXfat ~ ., data = bodyfat)
+(sbody <- stabsel(mod, q = 3, PFER = 0.2))
+dim(sbody$phat)
+(sbody <- stabsel(mod, q = 3, PFER = 0.2, error.bound = "SS"))
+dim(sbody$phat)
More information about the Mboost-commits
mailing list