[Mboost-commits] r851 - / pkg/mboostDevel/R pkg/mboostDevel/inst pkg/mboostDevel/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 22 19:29:07 CEST 2015
Author: hofner
Date: 2015-04-22 19:29:06 +0200 (Wed, 22 Apr 2015)
New Revision: 851
Modified:
README.md
pkg/mboostDevel/R/bkronecker.R
pkg/mboostDevel/R/bl.R
pkg/mboostDevel/R/bmono.R
pkg/mboostDevel/R/bmrf.R
pkg/mboostDevel/inst/NEWS.Rd
pkg/mboostDevel/tests/regtest-baselearner.R
Log:
(bbs) added linear extrapolation for prediction
and propagate changes of bl_lin() to bl_mono() and bl_lin_matrix()
Modified: README.md
===================================================================
--- README.md 2015-04-22 12:35:22 UTC (rev 850)
+++ README.md 2015-04-22 17:29:06 UTC (rev 851)
@@ -1,7 +1,7 @@
mboost
======
-[](https://travis-ci.org/hofnerb/mboost)
+[](https://travis-ci.org/hofnerb/mboost)
`mboost` implements boosting algorithms for fitting generalized linear, additive and interaction models
to potentially high-dimensional data.
Modified: pkg/mboostDevel/R/bkronecker.R
===================================================================
--- pkg/mboostDevel/R/bkronecker.R 2015-04-22 12:35:22 UTC (rev 850)
+++ pkg/mboostDevel/R/bkronecker.R 2015-04-22 17:29:06 UTC (rev 851)
@@ -8,12 +8,19 @@
index <- blg$get_index()
vary <- blg$get_vary()
- newX <- function(newdata = NULL) {
+ newX <- function(newdata = NULL, prediction = FALSE) {
if (!is.null(newdata)) {
- stopifnot(all(names(newdata) == names(blg)))
+ if (!all(names(blg) %in% names(newdata)))
+ stop("Variable(s) missing in ", sQuote("newdata"), ":\n\t",
+ names(blg)[!names(blg) %in% names(newdata)])
# stopifnot(all(class(newdata) == class(mf)))
- mf <- newdata[names(blg)]
+ nm <- names(blg)
+ if (any(duplicated(nm))) ## removes duplicates
+ nm <- unique(nm)
+ mf <- newdata[nm]
}
+ ## this argument is currently only used in X_bbs --> bsplines
+ args$prediction <- prediction
return(Xfun(mf, vary, args))
}
X <- newX()
@@ -34,7 +41,7 @@
dpp <- function(weights) {
- if (!is.null(attr(X$X1, "deriv")) || !is.null(attr(X$X2, "deriv")))
+ 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)
@@ -46,7 +53,7 @@
XtX <- array(XtX, c(c1, c1, c2, c2))
XtX <- mymatrix(aperm(XtX, c(1, 3, 2, 4)), nrow = c1 * c2)
- ### If lambda was given in both baselearners, we
+ ### 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
@@ -67,10 +74,10 @@
XtX <- XtX + K
### nnls
- constr <- (!is.null(attr(X$X1, "constraint"))) +
+ constr <- (!is.null(attr(X$X1, "constraint"))) +
(!is.null(attr(X$X2, "constraint")))
- if (constr == 2)
+ if (constr == 2)
stop("only one dimension may be subject to constraints")
constr <- constr > 0
@@ -137,7 +144,7 @@
index <- NULL
nm <- names(blg)
newdata <- newdata[nm]
- X <- newX(newdata)$X
+ X <- newX(newdata, prediction = TRUE)$X
}
ncfprod <- function(b)
as.vector(as(tcrossprod(X$X1 %*% b, X$X2), "matrix"))
@@ -227,7 +234,7 @@
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 ",
+ stop("lambda needs to be given in both baselearners combined with ",
sQuote("%O%"))
if (!is.null(l1) && !is.null(l2)) {
### there is no common lambda!
Modified: pkg/mboostDevel/R/bl.R
===================================================================
--- pkg/mboostDevel/R/bl.R 2015-04-22 12:35:22 UTC (rev 850)
+++ pkg/mboostDevel/R/bl.R 2015-04-22 17:29:06 UTC (rev 851)
@@ -206,9 +206,10 @@
if (cyclic & constraint != "none")
stop("constraints not implemented for cyclic B-splines")
stopifnot(is.numeric(deriv) & length(deriv) == 1)
+ ## prediction is usually set in/by newX()
list(knots = ret, degree = degree, differences = differences,
df = df, lambda = lambda, center = center, cyclic = cyclic,
- Ts_constraint = constraint, deriv = deriv)
+ Ts_constraint = constraint, deriv = deriv, prediction = FALSE)
}
### model.matrix for P-splines baselearner (including tensor product P-splines)
@@ -222,7 +223,7 @@
boundary.knots = args$knots[[i]]$boundary.knots,
degree = args$degree,
Ts_constraint = args$Ts_constraint,
- deriv = args$deriv)
+ deriv = args$deriv, extrapolation = args$prediction)
} else { ## if cyclic spline
X <- cbs(mf[[i]],
knots = args$knots[[i]]$knots,
@@ -569,7 +570,8 @@
### adapted version of mgcv::cSplineDes from S.N. Wood
cbs <- function (x, knots, boundary.knots, degree = 3, deriv = 0L) {
- if (any(x < boundary.knots[1]) | any(x > boundary.knots[2]))
+ if (any(x < boundary.knots[1], na.rm = TRUE) |
+ any(x > boundary.knots[2], na.rm = TRUE))
stop("some ", sQuote("x"), " values are beyond ",
sQuote("boundary.knots"))
@@ -605,19 +607,30 @@
attr(X, "degree") <- degree
attr(X,"knots") <- knots
attr(X,"boundary.knots") <- boundary.knots
- if (deriv != 0)
+ if (length(deriv) > 1 || deriv != 0)
attr(X, "deriv") <- deriv
dimnames(X) <- list(nx, 1L:ncol(X))
return(X)
}
bsplines <- function(x, knots, boundary.knots, degree,
- Ts_constraint = "none", deriv = 0L){
+ Ts_constraint = "none", deriv = 0L,
+ extrapolation = FALSE) {
- if (any(x < boundary.knots[1]) | any(x > boundary.knots[2]))
- warning("some ", sQuote("x"), " values are beyond ",
- sQuote("boundary.knots"))
+ ## do not allow data beyond boundary knots while fitting
+ if (!extrapolation && (any(x < boundary.knots[1], na.rm = TRUE) |
+ any(x > boundary.knots[2], na.rm = TRUE)))
+ stop("some ", sQuote("x"), " values are beyond ",
+ sQuote("boundary.knots"))
+ ## allow extrapolation when predicting
+ if (extrapolation <- extrapolation &&
+ (any(x < boundary.knots[1], na.rm = TRUE) |
+ any(x > boundary.knots[2], na.rm = TRUE))) {
+ warning("Some ", sQuote("x"), " values are beyond ",
+ sQuote("boundary.knots"), "; Linear extrapolation used.")
+ }
+
nx <- names(x)
x <- as.vector(x)
## handling of NAs
@@ -635,6 +648,27 @@
## construct design matrix
X <- splineDesign(k, x, degree + 1, derivs = rep(deriv, length(x)),
outer.ok = TRUE)
+
+ ## code along the lines of mgcv::Predict.matrix.pspline.smooth
+ if (extrapolation) {
+ ## Build matrix to map coeficients to value (deriv = 0) and
+ ## slope (deriv = 1) at end points.
+ if (deriv != 0L) {
+ warning("deriv != 0L; Linear extrapolation overwritten")
+ } else {
+ deriv <- c(0, 1, 0, 1)
+ }
+ D <- splineDesign(knots = k, x = rep(boundary.knots, each = 2),
+ ord = degree + 1, deriv)
+ ## Add rows for linear extrapolation
+ idx <- x < boundary.knots[1]
+ if (any(idx, na.rm = TRUE))
+ X[idx,] <- cbind(1, x[idx] - boundary.knots[1]) %*% D[1:2, ]
+ idx <- x > boundary.knots[2]
+ if (any(idx, na.rm = TRUE))
+ X[idx,] <- cbind(1, x[idx] - boundary.knots[2]) %*% D[3:4, ]
+ }
+
## handling of NAs
if (nas) {
tmp <- matrix(NA, length(nax), ncol(X))
@@ -655,7 +689,7 @@
attr(X, "Ts_constraint") <- Ts_constraint
if (Ts_constraint != "none")
attr(X, "D") <- D
- if (deriv != 0)
+ if (length(deriv) > 1 || deriv != 0)
attr(X, "deriv") <- deriv
dimnames(X) <- list(nx, 1L:ncol(X))
return(X)
@@ -668,7 +702,7 @@
index <- blg$get_index()
vary <- blg$get_vary()
- newX <- function(newdata = NULL) {
+ newX <- function(newdata = NULL, prediction = FALSE) {
if (!is.null(newdata)) {
if (!all(names(blg) %in% names(newdata)))
stop("Variable(s) missing in ", sQuote("newdata"), ":\n\t",
@@ -682,6 +716,8 @@
nm <- unique(nm)
mf <- newdata[, nm, drop = FALSE]
}
+ ## this argument is currently only used in X_bbs --> bsplines
+ args$prediction <- prediction
return(Xfun(mf, vary, args))
}
X <- newX()
@@ -773,7 +809,7 @@
newdata <- newdata[index[[1]],,drop = FALSE]
index <- index[[2]]
}
- X <- newX(newdata)$X
+ X <- newX(newdata, prediction = TRUE)$X
}
aggregate <- match.arg(aggregate)
pr <- switch(aggregate, "sum" =
Modified: pkg/mboostDevel/R/bmono.R
===================================================================
--- pkg/mboostDevel/R/bmono.R 2015-04-22 12:35:22 UTC (rev 850)
+++ pkg/mboostDevel/R/bmono.R 2015-04-22 17:29:06 UTC (rev 851)
@@ -134,8 +134,7 @@
}
}
args$cons.arg <- cons.arg
- ret$dpp <- bl_mono(ret, Xfun = X_bbs,
- args = args)
+ ret$dpp <- bl_mono(ret, Xfun = X_bbs, args = args)
} else {
args <- hyper_ols(df = df, lambda = lambda,
intercept = intercept,
@@ -147,8 +146,7 @@
## <FIXME> Was machen wir bei kateg. Effekten? Da muesste das doch auch gehen!
args$boundary.constraints <- boundary.constraints
args$cons.arg$n <- cons.arg$n
- ret$dpp <- bl_mono(ret, Xfun = X_ols,
- args = args)
+ ret$dpp <- bl_mono(ret, Xfun = X_ols, args = args)
}
return(ret)
}
@@ -158,12 +156,22 @@
index <- blg$get_index()
vary <- blg$get_vary()
- newX <- function(newdata = NULL) {
+ newX <- function(newdata = NULL, prediction = FALSE) {
if (!is.null(newdata)) {
- stopifnot(all(names(newdata) == names(blg)))
- stopifnot(all(class(newdata) == class(mf)))
- mf <- newdata[,names(blg),drop = FALSE]
+ if (!all(names(blg) %in% names(newdata)))
+ stop("Variable(s) missing in ", sQuote("newdata"), ":\n\t",
+ names(blg)[!names(blg) %in% names(newdata)])
+ if (!all(class(newdata) == class(mf)))
+ stop(sQuote("newdata"),
+ " must have the same class as the original data:\n\t",
+ class(mf))
+ nm <- names(blg)
+ if (any(duplicated(nm))) ## removes duplicates
+ nm <- unique(nm)
+ mf <- newdata[, nm, drop = FALSE]
}
+ ## this argument is currently only used in X_bbs --> bsplines
+ args$prediction <- prediction
return(Xfun(mf, vary, args))
}
X <- newX()
@@ -346,7 +354,7 @@
newdata <- newdata[index[[1]],,drop = FALSE]
index <- index[[2]]
}
- X <- newX(newdata)$X
+ X <- newX(newdata, prediction = TRUE)$X
}
aggregate <- match.arg(aggregate)
pr <- switch(aggregate, "sum" =
Modified: pkg/mboostDevel/R/bmrf.R
===================================================================
--- pkg/mboostDevel/R/bmrf.R 2015-04-22 12:35:22 UTC (rev 850)
+++ pkg/mboostDevel/R/bmrf.R 2015-04-22 17:29:06 UTC (rev 851)
@@ -1,7 +1,5 @@
-bmrf <-
-function (..., by = NULL, index = NULL, bnd = NULL, df = 4, lambda = NULL,
- center = FALSE)
-{
+bmrf <- function (..., by = NULL, index = NULL, bnd = NULL, df = 4,
+ lambda = NULL, center = FALSE) {
if (!requireNamespace("BayesX"))
stop("cannot load ", sQuote("BayesX"))
@@ -71,7 +69,7 @@
K <- Matrix(unclass(K))
nm <- colnames(mf)[colnames(mf) != vary]
list(K = K, bnd = bnd, pen = TRUE, df = df, lambda = lambda,
- center = center)
+ center = center)
}
X_bmrf <- function (mf, vary, args) {
Modified: pkg/mboostDevel/inst/NEWS.Rd
===================================================================
--- pkg/mboostDevel/inst/NEWS.Rd 2015-04-22 12:35:22 UTC (rev 850)
+++ pkg/mboostDevel/inst/NEWS.Rd 2015-04-22 17:29:06 UTC (rev 851)
@@ -10,6 +10,11 @@
\item Better handling of errors in (single) folds of \code{cvrisk}:
results of folds without errors are used and a \code{warning} is
issued.
+ \item \code{bbs} and \code{bmono} no longer allow data outside of
+ the \code{boundary.knots} during model fitting.
+ \item Predictions for \code{bbs} and \code{bmono} now use linear
+ extrapolation (user request inspired by
+ \code{mgcv::Predict.matrix.pspline.smooth}).
}
}
\subsection{Miscellaneous}{
Modified: pkg/mboostDevel/tests/regtest-baselearner.R
===================================================================
--- pkg/mboostDevel/tests/regtest-baselearner.R 2015-04-22 12:35:22 UTC (rev 850)
+++ pkg/mboostDevel/tests/regtest-baselearner.R 2015-04-22 17:29:06 UTC (rev 851)
@@ -477,3 +477,25 @@
round(extract(brandom(z1, df = 3)$dpp(rep(1, 100)), what = "df"), 2)
round(extract(brandom(z1, lambda = 50.39)$dpp(rep(1, 100)), what = "lambda"), 2)
round(extract(brandom(z1, lambda = 50.39)$dpp(rep(1, 100)), what = "df"), 2)
+
+
+### check if data beyond boundary knots is permitted
+set.seed(1234)
+x <- rnorm(100)
+y <- sin(x) + rnorm(100, sd = 0.1)
+plot(x, y, xlim = c(-3, 5))
+## should not work:
+try(mod <- mboost(y ~ bbs(x, boundary.knots = c(-1, 1))))
+try(mod <- mboost(y ~ bbs(x, cyclic = TRUE, boundary.knots = c(-1, 1))))
+## now fit models and check linear extrapolation
+mod <- mboost(y ~ bbs(x))
+tail(pr <- predict(mod, newdata = data.frame(x = seq(-3, 5, by = 0.1))))
+lines(seq(-3, 5, by = 0.1), pr)
+## now with bmono
+mod <- mboost(y ~ bmono(x))
+tail(pr2 <- predict(mod, newdata = data.frame(x = seq(-3, 5, by = 0.1))))
+lines(seq(-3, 5, by = 0.1), pr2, col = "red")
+## check same with cyclic splines
+mod <- mboost(y ~ bbs(x, cyclic = TRUE))
+try(predict(mod, newdata = data.frame(x = seq(-3, 5, by = 0.1))))
+
More information about the Mboost-commits
mailing list