[Mboost-commits] r774 - / pkg/mboostPatch pkg/mboostPatch/R pkg/mboostPatch/inst pkg/mboostPatch/man pkg/mboostPatch/tests pkg/mboostPatch/tests/Examples pkg/mboostPatch/vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 27 19:06:32 CEST 2014
Author: hofner
Date: 2014-06-27 19:06:32 +0200 (Fri, 27 Jun 2014)
New Revision: 774
Added:
pkg/mboostPatch/tests/regtest-inference.R
pkg/mboostPatch/tests/regtest-inference.Rout.save
Removed:
pkg/mboostPatch/man/Westbc.Rd
pkg/mboostPatch/man/birds.Rd
pkg/mboostPatch/man/bodyfat.Rd
pkg/mboostPatch/man/wpbc.Rd
Modified:
pkg/mboostPatch/.Rbuildignore
pkg/mboostPatch/.RbuildignoreCRAN
pkg/mboostPatch/DESCRIPTION
pkg/mboostPatch/NAMESPACE
pkg/mboostPatch/R/AAA.R
pkg/mboostPatch/R/bkronecker.R
pkg/mboostPatch/R/bl.R
pkg/mboostPatch/R/bmono.R
pkg/mboostPatch/R/control.R
pkg/mboostPatch/R/crossvalidation.R
pkg/mboostPatch/R/family.R
pkg/mboostPatch/R/helpers.R
pkg/mboostPatch/R/inference.R
pkg/mboostPatch/R/mboost.R
pkg/mboostPatch/R/methods.R
pkg/mboostPatch/inst/CHANGES
pkg/mboostPatch/inst/birds_Biometrics.R
pkg/mboostPatch/man/FP.Rd
pkg/mboostPatch/man/Family.Rd
pkg/mboostPatch/man/baselearners.Rd
pkg/mboostPatch/man/control.Rd
pkg/mboostPatch/man/cvrisk.Rd
pkg/mboostPatch/man/gamboost.Rd
pkg/mboostPatch/man/glmboost.Rd
pkg/mboostPatch/man/mboost.Rd
pkg/mboostPatch/man/mboost_package.Rd
pkg/mboostPatch/man/methods.Rd
pkg/mboostPatch/man/stabsel.Rd
pkg/mboostPatch/tests/Examples/mboost-Ex.Rout.save
pkg/mboostPatch/tests/birds_Biometrics.Rout.save
pkg/mboostPatch/tests/bugfixes.R
pkg/mboostPatch/tests/bugfixes.Rout.save
pkg/mboostPatch/tests/regtest-baselearner.R
pkg/mboostPatch/tests/regtest-baselearner.Rout.save
pkg/mboostPatch/tests/regtest-blackboost.Rout.save
pkg/mboostPatch/tests/regtest-family.R
pkg/mboostPatch/tests/regtest-family.Rout.save
pkg/mboostPatch/tests/regtest-gamboost.R
pkg/mboostPatch/tests/regtest-gamboost.Rout.save
pkg/mboostPatch/tests/regtest-glmboost.Rout.save
pkg/mboostPatch/tests/regtest-hatmatrix.Rout.save
pkg/mboostPatch/vignettes/SurvivalEnsembles.Rout.save
pkg/mboostPatch/vignettes/mboost.Rnw
pkg/mboostPatch/vignettes/mboost.Rout.save
pkg/mboostPatch/vignettes/mboost_illustrations.Rnw
pkg/mboostPatch/vignettes/mboost_illustrations.Rout.save
pkg/mboostPatch/vignettes/mboost_tutorial.Rnw
pkg/mboostPatch/vignettes/mboost_tutorial.Rout.save
pkg/mboostPatch/vignettes/setup.R
svn_release.txt
Log:
merge mboostDevel to mboostPatch
Modified: pkg/mboostPatch/.Rbuildignore
===================================================================
--- pkg/mboostPatch/.Rbuildignore 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/.Rbuildignore 2014-06-27 17:06:32 UTC (rev 774)
@@ -1,3 +1,4 @@
demo
to_do_list.txt
-^\..*
\ No newline at end of file
+^\..*
+.*/auto
Modified: pkg/mboostPatch/.RbuildignoreCRAN
===================================================================
--- pkg/mboostPatch/.RbuildignoreCRAN 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/.RbuildignoreCRAN 2014-06-27 17:06:32 UTC (rev 774)
@@ -2,4 +2,5 @@
to_do_list.txt
test
^\..*
-vignettes/.*\.Rout\.save$
\ No newline at end of file
+.*/auto
+vignettes/.*\.Rout\.save$
Modified: pkg/mboostPatch/DESCRIPTION
===================================================================
--- pkg/mboostPatch/DESCRIPTION 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/DESCRIPTION 2014-06-27 17:06:32 UTC (rev 774)
@@ -1,7 +1,7 @@
Package: mboost
Title: Model-Based Boosting
-Version: 2.2-4
-Date: 2014-04-15
+Version: 2.3-0
+Date: 2014-06-26
Authors at R: c(person("Torsten", "Hothorn", role = c("aut", "cre"),
email = "Torsten.Hothorn at R-project.org"),
person("Peter", "Buehlmann", role = "aut"),
@@ -15,10 +15,10 @@
component-wise (penalised) least squares estimates or regression
trees as base-learners for fitting generalized linear, additive
and interaction models to potentially high-dimensional data.
-Depends: R (>= 2.14.0), methods, stats, parallel, survival
-Imports: Matrix, splines, lattice, utils
-Suggests: party (>= 1.0-3), TH.data, MASS, fields,
- BayesX, gbm, mlbench, RColorBrewer, rpart (>= 4.0-3)
+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)
LazyData: yes
License: GPL-2
-URL: http://r-forge.r-project.org/projects/mboost/
\ No newline at end of file
+URL: http://r-forge.r-project.org/projects/mboost/
Modified: pkg/mboostPatch/NAMESPACE
===================================================================
--- pkg/mboostPatch/NAMESPACE 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/NAMESPACE 2014-06-27 17:06:32 UTC (rev 774)
@@ -6,6 +6,8 @@
importFrom(survival, Surv, survfit)
importFrom(splines, bs, splineDesign)
importFrom(lattice, levelplot)
+importFrom(nnls, nnls)
+importFrom(quadprog, solve.QP)
importFrom(utils, packageDescription)
export(glmboost,
@@ -15,9 +17,10 @@
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, 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)
+ nuisance, "%+%", "%X%", "%O%", extract, risk, "mstop<-")
###, basesel, fitsel)
exportClasses("boost_family")
exportMethods("show")
@@ -68,7 +71,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)
@@ -76,5 +81,6 @@
S3method(extract, bl_lin)
S3method(extract, bl_tree)
S3method(residuals, mboost)
+S3method(risk, mboost)
useDynLib(mboost)
Modified: pkg/mboostPatch/R/AAA.R
===================================================================
--- pkg/mboostPatch/R/AAA.R 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/R/AAA.R 2014-06-27 17:06:32 UTC (rev 774)
@@ -31,11 +31,7 @@
packageStartupMessage("This is mboost ", vers, ". ", "See ",
sQuote("package?mboost"), " and the NEWS file\n",
"for a complete list of changes.\n",
- "Note: The default for the computation",
- " of the degrees of freedom has changed.\n",
- " For details see section ",
- sQuote("Global Options"), " of ",
- sQuote("?bols"), ".", appendLF = TRUE)
+ appendLF = TRUE)
return(TRUE)
}
Modified: pkg/mboostPatch/R/bkronecker.R
===================================================================
--- pkg/mboostPatch/R/bkronecker.R 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/R/bkronecker.R 2014-06-27 17:06:32 UTC (rev 774)
@@ -34,6 +34,9 @@
dpp <- function(weights) {
+ if (!is.null(attr(X$X1, "deriv")) || !is.null(attr(X$X2, "deriv")))
+ stop("fitting of derivatives of B-splines not implemented")
+
W <- matrix(weights, nrow = n1, ncol = n2)
### X = kronecker(X2, X1)
@@ -43,10 +46,11 @@
XtX <- array(XtX, c(c1, c1, c2, c2))
XtX <- mymatrix(aperm(XtX, c(1, 3, 2, 4)), nrow = c1 * c2)
- ### <FIXME> This does not happen in bl_lin / df2lambda.
- ### For one base learner only, it makes sense to allow
- ### for a direct choice of lambda (regardless of df)
- ### </FIXME>
+ ### If lambda was given in both baselearners, we
+ ### directly multiply the marginal penalty matrices by lambda
+ ### and then compute the total penalty as the kronecker sum.
+ ### args$lambda is NA in this case and we don't compute
+ ### the corresponding df's (unlike bl_lin)
if (is.null(args$lambda)) {
### <FIXME>: is there a better way to feed XtX into lambdadf?
@@ -55,17 +59,29 @@
dmat = K, weights = weights, XtX = XtX)
### </FIXME>
lambda <- lambdadf["lambda"]
+ K <- lambda * K
} else {
- lambda <- args$lambda
+ lambdadf <- args[c("lambda", "df")]
}
- XtX <- XtX + lambda * K
+ ### note: K already contains the lambda penalty parameter(s)
+ XtX <- XtX + K
+ ### nnls
+ constr <- (!is.null(attr(X$X1, "constraint"))) +
+ (!is.null(attr(X$X2, "constraint")))
+
+ if (constr == 2)
+ stop("only one dimension may be subject to constraints")
+ constr <- constr > 0
+
## matrizes of class dgeMatrix are dense generic matrices; they should
## be coerced to class matrix and handled in the standard way
if (is(XtX, "Matrix") && !extends(class(XtX), "dgeMatrix")) {
XtXC <- Cholesky(forceSymmetric(XtX))
mysolve <- function(y) {
Y <- matrix(y, nrow = n1) * W
+ if (constr)
+ return(nnls2D(X, as(XtXC, "matrix"), Y))
XWY <- as.vector(crossprod(X$X1, Y) %*% X$X2)
solve(XtXC, XWY) ## special solve routine from
## package Matrix
@@ -77,6 +93,8 @@
}
mysolve <- function(y) {
Y <- matrix(y, nrow = n1) * W
+ if (constr)
+ return(nnls2D(X, as(XtX, "matrix"), Y))
XWY <- crossprod(X$X1, Y) %*% X$X2
solve(XtX, matrix(as(XWY, "matrix"), ncol = 1),
LINPACK = FALSE)
@@ -208,8 +226,12 @@
args2 <- environment(bl2$dpp)$args
l1 <- args1$lambda
l2 <- args2$lambda
+ if (xor(is.null(l1), is.null(l2)))
+ stop("lambda needs to be given in both baselearners combined with ",
+ sQuote("%O%"))
if (!is.null(l1) && !is.null(l2)) {
- args <- list(lambda = l1 + l2, df = NULL)
+ ### there is no common lambda!
+ args <- list(lambda = NA, df = NA)
} else {
args <- list(lambda = NULL,
df = ifelse(is.null(args1$df), 1, args1$df) *
Modified: pkg/mboostPatch/R/bl.R
===================================================================
--- pkg/mboostPatch/R/bl.R 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/R/bl.R 2014-06-27 17:06:32 UTC (rev 774)
@@ -167,7 +167,7 @@
### hyper parameters for P-splines baselearner (including tensor product P-splines)
hyper_bbs <- function(mf, vary, knots = 20, boundary.knots = NULL, degree = 3,
differences = 2, df = 4, lambda = NULL, center = FALSE,
- cyclic = FALSE) {
+ cyclic = FALSE, constraint = "none", deriv = 0L) {
knotf <- function(x, knots, boundary.knots) {
if (is.null(boundary.knots))
@@ -189,7 +189,7 @@
stop("variable names and knot names must be the same")
if (is.list(boundary.knots)) if(!all(names(boundary.knots) %in% nm))
stop("variable names and boundary.knot names must be the same")
- if (isTRUE(center) && cyclic)
+ if (!identical(center, FALSE) && cyclic)
stop("centering of cyclic covariates not yet implemented")
ret <- vector(mode = "list", length = length(nm))
names(ret) <- nm
@@ -197,8 +197,12 @@
ret[[n]] <- knotf(mf[[n]], if (is.list(knots)) knots[[n]] else knots,
if (is.list(boundary.knots)) boundary.knots[[n]]
else boundary.knots)
+ if (cyclic & constraint != "none")
+ stop("constraints not implemented for cyclic B-splines")
+ stopifnot(is.numeric(deriv) & length(deriv) == 1)
list(knots = ret, degree = degree, differences = differences,
- df = df, lambda = lambda, center = center, cyclic = cyclic)
+ df = df, lambda = lambda, center = center, cyclic = cyclic,
+ Ts_constraint = constraint, deriv = deriv)
}
### model.matrix for P-splines baselearner (including tensor product P-splines)
@@ -209,12 +213,15 @@
X <- bsplines(mf[[i]],
knots = args$knots[[i]]$knots,
boundary.knots = args$knots[[i]]$boundary.knots,
- degree = args$degree)
+ degree = args$degree,
+ Ts_constraint = args$Ts_constraint,
+ deriv = args$deriv)
if (args$cyclic) {
X <- cbs(mf[[i]],
knots = args$knots[[i]]$knots,
boundary.knots = args$knots[[i]]$boundary.knots,
- degree = args$degree)
+ degree = args$degree,
+ deriv = args$deriv)
}
class(X) <- "matrix"
return(X)
@@ -269,7 +276,7 @@
if (vary != "" && ncol(by) > 1){ # build block diagonal penalty
suppressMessages(K <- kronecker(diag(ncol(by)), K))
}
- if (isTRUE(args$center)) {
+ if (!identical(args$center, FALSE)) {
tmp <- attributes(X)[c("degree", "knots", "Boundary.knots")]
center <- match.arg(as.character(args$center),
choices = c("TRUE", "differenceMatrix", "spectralDecomp"))
@@ -291,6 +298,10 @@
} else {
K <- crossprod(K)
}
+ if (!is.null(attr(X, "Ts_constraint"))) {
+ D <- attr(X, "D")
+ K <- crossprod(D, K) %*% D
+ }
}
if (length(mm) == 2) {
suppressMessages(
@@ -346,7 +357,7 @@
if (vary != "" && ncol(by) > 1){ # build block diagonal penalty
suppressMessages(K <- kronecker(diag(ncol(by)), K))
}
- if (isTRUE(args$center)) {
+ if (!identical(args$center, FALSE)) {
### L = \Gamma \Omega^1/2 in Section 2.3. of Fahrmeir et al.
### (2004, Stat Sinica), always
L <- eigen(K, symmetric = TRUE, EISPACK = FALSE)
@@ -467,7 +478,8 @@
### P-spline (and tensor-product spline) baselearner
bbs <- function(..., by = NULL, index = NULL, knots = 20, boundary.knots = NULL,
degree = 3, differences = 2, df = 4, lambda = NULL, center = FALSE,
- cyclic = FALSE) {
+ cyclic = FALSE, constraint = c("none", "increasing", "decreasing"),
+ deriv = 0) {
if (!is.null(lambda)) df <- NULL
@@ -541,13 +553,14 @@
ret$dpp <- bl_lin(ret, Xfun = X_bbs,
args = hyper_bbs(mf, vary, knots = knots, boundary.knots =
boundary.knots, degree = degree, differences = differences,
- df = df, lambda = lambda, center = center, cyclic = cyclic))
+ df = df, lambda = lambda, center = center, cyclic = cyclic,
+ constraint = match.arg(constraint), deriv = deriv))
return(ret)
}
### cyclic B-splines
### adapted version of mgcv:cSplineDes from S.N. Wood
-cbs <- function (x, knots, boundary.knots, degree = 3) {
+cbs <- function (x, knots, boundary.knots, degree = 3, deriv = 0L) {
# require(splines)
nx <- names(x)
x <- as.vector(x)
@@ -564,10 +577,11 @@
(boundary.knots[2] - knots[(nKnots - ord + 1):(nKnots - 1)]),
knots)
ind <- x > xc
- X <- splineDesign(knots, x, ord, outer.ok = TRUE)
+ X <- splineDesign(knots, x, ord, derivs = rep(deriv, length(x)), outer.ok = TRUE)
x[ind] <- x[ind] - boundary.knots[2] + boundary.knots[1]
if (sum(ind)) {
- Xtmp <- splineDesign(knots, x[ind], ord, outer.ok = TRUE)
+ Xtmp <- splineDesign(knots, x[ind], ord, derivs = rep(deriv, length(x[ind])),
+ outer.ok = TRUE)
X[ind, ] <- X[ind, ] + Xtmp
}
## handling of NAs
@@ -580,11 +594,14 @@
attr(X, "degree") <- degree
attr(X,"knots") <- knots
attr(X,"boundary.knots") <- boundary.knots
+ if (deriv != 0)
+ attr(X, "deriv") <- deriv
dimnames(X) <- list(nx, 1L:ncol(X))
return(X)
}
-bsplines <- function(x, knots, boundary.knots, degree){
+bsplines <- function(x, knots, boundary.knots, degree,
+ Ts_constraint = "none", deriv = 0L){
nx <- names(x)
x <- as.vector(x)
## handling of NAs
@@ -600,17 +617,30 @@
## complete knot mesh
k <- c(bk_lower, knots, bk_upper)
## construct design matrix
- X <- splineDesign(k, x, degree + 1, outer.ok = TRUE)
+ X <- splineDesign(k, x, degree + 1, derivs = rep(deriv, length(x)),
+ outer.ok = TRUE)
## handling of NAs
if (nas) {
tmp <- matrix(NA, length(nax), ncol(X))
tmp[!nax, ] <- X
X <- tmp
}
+ ### constraints; experimental
+ D <- diag(ncol(X))
+ D[lower.tri(D)] <- 1
+ X <- switch(Ts_constraint, "none" = X,
+ "increasing" = X %*% D,
+ "decreasing" = -X %*% D)
## add attributes
attr(X, "degree") <- degree
- attr(X,"knots") <- knots
- attr(X,"boundary.knots") <- list(lower = bk_lower, upper = bk_upper)
+ attr(X, "knots") <- knots
+ attr(X, "boundary.knots") <- list(lower = bk_lower, upper = bk_upper)
+ if (Ts_constraint != "none")
+ attr(X, "Ts_constraint") <- Ts_constraint
+ if (Ts_constraint != "none")
+ attr(X, "D") <- D
+ if (deriv != 0)
+ attr(X, "deriv") <- deriv
dimnames(X) <- list(nx, 1L:ncol(X))
return(X)
}
@@ -639,6 +669,9 @@
dpp <- function(weights) {
+ if (!is.null(attr(X, "deriv")))
+ stop("fitting of derivatives of B-splines not implemented")
+
weights[!Complete.cases(mf)] <- 0
w <- weights
if (!is.null(index))
@@ -653,17 +686,25 @@
## be coerced to class matrix and handled in the standard way
if (is(X, "Matrix") && !extends(class(XtX), "dgeMatrix")) {
XtXC <- Cholesky(forceSymmetric(XtX))
- mysolve <- function(y)
- solve(XtXC, crossprod(X, y)) ## special solve routine from
- ## package Matrix
+ mysolve <- function(y) {
+ if (is.null(attr(X, "Ts_constraint")))
+ return(solve(XtXC, crossprod(X, y))) ## special solve routine from
+ ## package Matrix
+ ### non-negative LS only at the moment
+ return(nnls1D(as(XtX, "matrix"), as(X, "matrix"), y))
+ }
} else {
if (is(X, "Matrix")) {
## coerce Matrix to matrix
X <- as(X, "matrix")
XtX <- as(XtX, "matrix")
}
- mysolve <- function(y)
- solve(XtX, crossprod(X, y), LINPACK = FALSE)
+ mysolve <- function(y) {
+ if (is.null(attr(X, "Ts_constraint")))
+ return(solve(XtX, crossprod(X, y), LINPACK = FALSE))
+ ### non-negative LS only at the moment
+ return(nnls1D(XtX, X, y))
+ }
}
fit <- function(y) {
@@ -747,9 +788,21 @@
}
### random-effects (Ridge-penalized ANOVA) baselearner
-brandom <- function (..., contrasts.arg = "contr.dummy", df = 4) {
+brandom <- function(..., by = NULL, index = NULL, df = 4, lambda = NULL,
+ contrasts.arg = "contr.dummy") {
cl <- cltmp <- match.call()
- if (is.null(cl$df))
+ x <- list(...)
+ ## drop further arguments to be passed to bols
+ if (!is.null(names(x)))
+ x <- x[names(x) == ""]
+
+ if (!all(sapply(x, is.factor) |
+ sapply(x, is.matrix) |
+ sapply(x, is.data.frame)))
+ stop(sQuote("..."), " must be a factor or design matrix in ",
+ sQuote("brandom"))
+
+ if (is.null(cl$df) && is.null(cl$lambda))
cl$df <- df
if (is.null(cl$contrasts.arg))
cl$contrasts.arg <- contrasts.arg
Modified: pkg/mboostPatch/R/bmono.R
===================================================================
--- pkg/mboostPatch/R/bmono.R 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/R/bmono.R 2014-06-27 17:06:32 UTC (rev 774)
@@ -1,6 +1,8 @@
### P-spline base-learner with (monotonicity) constraints
bmono <- function(..., constraint = c("increasing", "decreasing",
- "convex", "concave", "none"),
+ "convex", "concave", "none",
+ "positive", "negative"),
+ 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,
@@ -20,7 +22,7 @@
if (!is.list(constraint)) {
constraint <- match.arg(constraint)
} else {
- c.args <- eval(formals(sys.function(sys.parent()))[["constraint"]])
+ c.args <- eval(formals(sys.function())[["constraint"]])
constraint <- lapply(constraint, match.arg, choices = c.args)
}
if (length(mf) > 1){
@@ -32,6 +34,8 @@
}
##
+ type = match.arg(type)
+
if (length(mf) == 1 && (is.matrix(mf[[1]]) || is.data.frame(mf[[1]]))) {
mf <- as.data.frame(mf[[1]])
} else {
@@ -104,10 +108,11 @@
degree = degree, differences = differences,
df = df, lambda = lambda, center = FALSE)
args$constraint <- constraint
+ args$type <- type
args$lambda2 <- lambda2
args$niter <- niter
args$boundary.constraints <- boundary.constraints
- if(boundary.constraints){
+ if(boundary.constraints) {
if (is.null(cons.arg$n)){
## use 10% of the knots on each side per default
cons.arg$n <- sapply(args$knots,
@@ -118,14 +123,11 @@
}
## <fixme> was passiert bei bivariatem bmono? </fixme>
}
+ ## diff_order for boundary constraints
if(is.null(cons.arg$diff_order)){
## use same difference order as defined by "constraint":
## <FIXME> args$constraint may be a list of length 2 for spatial effects
- if (args$constraint %in% c("increasing", "decreasing")){
- cons.arg$diff_order <- 1
- } else { # i.e. args$constraint %in% c("convex", "concave")
- cons.arg$diff_order <- 2
- }
+ cons.arg$diff_order <- differences(args$constraint)
}
if(is.null(cons.arg$lambda)){
cons.arg$lambda <- 1e6
@@ -139,9 +141,10 @@
intercept = intercept,
contrasts.arg = contrasts.arg)
args$constraint <- constraint
+ args$type <- type
args$lambda2 <- lambda2
args$niter <- niter
- ## <FIXME> Was machen wir bei cat. Effekten? Da müsste das doch auch gehen!
+ ## <FIXME> Was machen wir bei kateg. Effekten? Da müsste das doch auch gehen!
args$boundary.constraints <- boundary.constraints
args$cons.arg$n <- cons.arg$n
ret$dpp <- bl_mono(ret, Xfun = X_ols,
@@ -168,45 +171,35 @@
X <- X$X
if (length(args$constraint) == 1) {
+ D <- V <- lambda2 <- vector(mode = "list", length = 2)
+ ## set up difference matrix
+ D[[1]] <- differences(args$constraint, ncol(X))
- if (args$constraint %in% c("increasing", "decreasing")){
- diff_order <- 1
- } else { # i.e. args$constraint %in% c("convex", "concave")
- diff_order <- 2
- }
-
- D <- V <- lambda2 <- vector(mode = "list", length =2)
-
if (is.factor(mf[[1]]) && args$intercept) {
- D[[1]] <- diff(diag(ncol(X)), differences = diff_order)
D[[1]][1,1] <- 0
- } else {
- D[[1]] <- diff(diag(ncol(X)), differences = diff_order)
+ }
+ if (!is.factor(mf[[1]]) && args$boundary.constraints) {
## set up boundary constraints
- if (args$boundary.constraints){
- cons.arg <- args$cons.arg
- idx <- rep(0, ncol(X) - cons.arg$diff_order)
- if (cons.arg$n[1] == 0) {
- lower <- 0
- } else {
- lower <- 1:cons.arg$n[1]
- }
- if (cons.arg$n[2] == length(idx)) {
- upper <- 0
- } else {
- upper <- length(idx) - 1:cons.arg$n[2] + 1
- }
- idx[c(lower, upper)] <- 1
- V3 <- diag(idx)
- D3 <- V3 %*% diff(diag(ncol(X)),
- differences = cons.arg$diff_order)
+ cons.arg <- args$cons.arg
+ idx <- rep(0, ncol(X) - cons.arg$diff_order)
+ if (cons.arg$n[1] == 0) {
+ lower <- 0
+ } else {
+ lower <- 1:cons.arg$n[1]
}
+ if (cons.arg$n[2] == length(idx)) {
+ upper <- 0
+ } else {
+ upper <- length(idx) - 1:cons.arg$n[2] + 1
+ }
+ idx[c(lower, upper)] <- 1
+ V3 <- diag(idx)
+ D3 <- V3 %*% diff(diag(ncol(X)), differences = cons.arg$diff_order)
}
- V[[1]] <- matrix(0, ncol = nrow(D[[1]]),
- nrow = nrow(D[[1]]))
- lambda2[[1]] <- args$lambda2
+ V[[1]] <- matrix(0, ncol = nrow(D[[1]]), nrow = nrow(D[[1]]))
+ lambda2[[1]] <- ifelse(args$constraint == "none", 0, args$lambda2)
lambda2[[2]] <- 0
if (args$boundary.constraints) {
lambda3 <- cons.arg$lambda
@@ -215,37 +208,32 @@
}
}
if (length(args$constraint) == 2) {
- diff_order <- lapply(args$constraint, function(x){
- ifelse( x %in% c("increasing", "decreasing"), 1, 2) } )
-
if (is.factor(mf[[1]]))
stop(paste("Bivariate monotonic effects currently not",
"implemented for ordered factors"))
## ncol1 = length(knots[[1]]) + degree + 1
- ## ncol2 = length(knots[[1]]) + degree + 1
+ ## ncol2 = length(knots[[2]]) + degree + 1
## ncol(X) = ncol1 * ncol2
ncoli <- lapply(args$knots, function(x)
length(x$knots) + args$degree + 1)
stopifnot(ncoli[[1]] * ncoli[[2]] == ncol(X))
D <- V <- lambda2 <- vector(mode = "list", length =2)
- suppressMessages(
- D[[1]] <- kronecker(diff(diag(ncoli[[1]]),
- differences = diff_order[[1]]),
- diag(ncoli[[2]]))
- )
- suppressMessages(
- D[[2]] <- kronecker(diag(ncoli[[1]]),
- diff(diag(ncoli[[2]]),
- differences = diff_order[[2]]))
- )
- V[[1]] <- matrix(0, ncol = nrow(D[[1]]), nrow = nrow(D[[1]]))
- V[[2]] <- matrix(0, ncol = nrow(D[[2]]), nrow = nrow(D[[2]]))
+ ## set up difference matrices
+ D <- differences(args$constraint, ncoli)
+ idx <- !sapply(D, is.null)
+ V[idx] <- lapply(D[idx], function(m) matrix(0, nrow(m), nrow(m)))
+
if (length(args$lambda2) == 1) {
lambda2[[1]] <- lambda2[[2]] <- args$lambda2
} else {
lambda2 <- args$lambda2
}
+ ## 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)
@@ -266,30 +254,11 @@
XtX <- crossprod(X * w, X)
XtX <- XtX + lambda * K
- ## Define solvers:
- ## define function as text and eval(parse()) later.
- l2txt <- "+ lambda2[[2]] * crossprod(D[[2]], V[[2]] %*% D[[2]])"
- l3txt <- "+ lambda3 * crossprod(D3, V3 %*% D3)"
- fct <- c("function(y, V) {",
- " XtXC <- Cholesky(forceSymmetric(XtX +",
- " lambda2[[1]] * crossprod(D[[1]], V[[1]] %*% D[[1]])",
- ## add if lambda2[[2]] != 0
- ifelse(lambda2[[2]] != 0, l2txt,""),
- ## add if lambda3 != 0
- ifelse(lambda3 != 0, l3txt,""),
- " ))",
- " solve(XtXC, crossprod(X, y), LINPACK = FALSE)",
- "}"
- )
-
- if (!is(X, "Matrix")) {
- ## some lines must be replaced in order to solve directly
- fct[2] <- ' solve(XtX +'
- fct[6] <- " , crossprod(X, y),"
- fct[7] <- ' LINPACK = FALSE)'
+ if (args$type == "iterative") {
+ fct <- define_solver(lambda2, lambda3, X)
+ ## deparsing and parsing again needed to tidy-up code.
+ mysolve <- eval(parse(text = deparse(eval(parse(text = fct)))))
}
- ## deparsing and parsing again needed to tidy-up code.
- mysolve <- eval(parse(text = deparse(eval(parse(text = fct)))))
fit <- function(y) {
if (!is.null(index)) {
@@ -299,32 +268,41 @@
y <- y * weights
}
- for (i in 1:args$niter){
- coef <- mysolve(y, V)
- ## compare old and new V
- tmp1 <- do.call(args$constraint[[1]],
- args=list(D[[1]] %*% coef))
- if (lambda2[[2]] != 0)
- tmp2 <- do.call(args$constraint[[2]],
- args=list(D[[2]] %*% coef))
+ if (args$type == "iterative") {
+ for (i in 1:args$niter){
+ coef <- mysolve(y, V)
+ if (args$constraint[[1]] == "none")
+ break ## as there is no need to iterate
- if ( all( V[[1]] == tmp1 ) &&
- ( lambda2[[2]] == 0 || all( V[[2]] == tmp2 ) ) )
- break # if both are equal: done!
- #if (args$boundary.constraints &&
- # all( V[[1]][-idxB, -idxB] == tmp1[-idxB, -idxB]) )
- # break # if both are equal (without V for boundary
- # # constraints): done!
- V[[1]] <- tmp1
- #if (args$boundary.constraints) {
- # V[[1]][idxFlat, idxFlat] <- diag(rep(1, length(idxFlat)))
- #}
- if (lambda2[[2]] != 0)
- V[[2]] <- tmp2
- if (i == args$niter)
- warning("no convergence of coef in bmono\n",
- "You could try increasing ", sQuote("niter"),
- " or ", sQuote("lambda2"))
+ ## compare old and new V
+ tmp1 <- violations(D[[1]] %*% coef)
+ if (lambda2[[2]] != 0)
+ tmp2 <- violations(D[[2]] %*% coef)
+
+ if ( all( V[[1]] == tmp1 ) &&
+ ( lambda2[[2]] == 0 || all( V[[2]] == tmp2 ) ) )
+ break # if both are equal: done!
+ #if (args$boundary.constraints &&
+ # all( V[[1]][-idxB, -idxB] == tmp1[-idxB, -idxB]) )
+ # break # if both are equal (without V for boundary
+ # # constraints): done!
+ V[[1]] <- tmp1
+ #if (args$boundary.constraints) {
+ # V[[1]][idxFlat, idxFlat] <- diag(rep(1, length(idxFlat)))
+ #}
+ if (lambda2[[2]] != 0)
+ V[[2]] <- tmp2
+ if (i == args$niter)
+ warning("no convergence of coef in bmono\n",
+ "You could try increasing ", sQuote("niter"),
+ " or ", sQuote("lambda2"))
+ }
+ } else { ## i.e. type == "quad.prog"
+ if (lambda2[[2]] == 0) {
+ coef <- solveLSEI(XtX, crossprod(X, y), D = D[[1]])
+ } else {
+ coef <- solveLSEI(XtX, crossprod(X, y), D = D)
+ }
}
ret <- list(model = coef,
@@ -392,17 +370,30 @@
return(dpp)
}
-none <- function(diffs)
- diag(rep(0,length(diffs)))
-
-increasing <- function(diffs)
+violations <- function(diffs)
diag(c(as.numeric(diffs)) <= 0)
-decreasing <- function(diffs)
- diag(c(as.numeric(diffs)) >= 0)
+define_solver <- function(lambda2, lambda3, X) {
+ ## define function as text and eval(parse()) later.
+ l2txt <- "+ lambda2[[2]] * crossprod(D[[2]], V[[2]] %*% D[[2]])"
+ l3txt <- "+ lambda3 * crossprod(D3, V3 %*% D3)"
+ fct <- c("function(y, V) {",
+ " XtXC <- Cholesky(forceSymmetric(XtX +",
+ " lambda2[[1]] * crossprod(D[[1]], V[[1]] %*% D[[1]])",
+ ## add if lambda2[[2]] != 0
+ ifelse(lambda2[[2]] != 0, l2txt,""),
+ ## add if lambda3 != 0
+ ifelse(lambda3 != 0, l3txt,""),
+ " ))",
+ " solve(XtXC, crossprod(X, y), LINPACK = FALSE)",
+ "}"
+ )
-convex <- function(diffs)
- diag(c(as.numeric(diffs)) <= 0)
-
-concave <- function(diffs)
- diag(c(as.numeric(diffs)) >= 0)
+ if (!is(X, "Matrix")) {
+ ## some lines must be replaced in order to solve directly
+ fct[2] <- ' solve(XtX +'
+ fct[6] <- " , crossprod(X, y),"
+ fct[7] <- ' LINPACK = FALSE)'
+ }
+ fct
+}
Modified: pkg/mboostPatch/R/control.R
===================================================================
--- pkg/mboostPatch/R/control.R 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/R/control.R 2014-06-27 17:06:32 UTC (rev 774)
@@ -1,12 +1,14 @@
boost_control <- function(mstop = 100, nu = 0.1,
risk = c("inbag", "oobag", "none"),
+ stopintern = FALSE,
center = TRUE, trace = FALSE) {
risk <- match.arg(risk)
+ stopintern <- stopintern & (risk == "oobag")
RET <- list(mstop = mstop, nu = nu,
- risk = risk, center = center,
- trace = trace)
+ risk = risk, stopintern = stopintern,
+ center = center, trace = trace)
class(RET) <- c("boost_control")
RET
}
Modified: pkg/mboostPatch/R/crossvalidation.R
===================================================================
--- pkg/mboostPatch/R/crossvalidation.R 2014-06-25 19:42:07 UTC (rev 773)
+++ pkg/mboostPatch/R/crossvalidation.R 2014-06-27 17:06:32 UTC (rev 774)
@@ -39,6 +39,7 @@
fun(mod)
}
}
+ ## use case weights as out-of-bag weights (but set inbag to 0)
OOBweights <- matrix(rep(weights, ncol(folds)), ncol = ncol(folds))
OOBweights[folds > 0] <- 0
oobrisk <- papply(1:ncol(folds),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/mboost -r 774
More information about the Mboost-commits
mailing list