[Mboost-commits] r778 - in pkg/mboostDevel: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 21 19:22:26 CEST 2014
Author: hofner
Date: 2014-07-21 19:22:26 +0200 (Mon, 21 Jul 2014)
New Revision: 778
Added:
pkg/mboostDevel/man/confint.Rd
Modified:
pkg/mboostDevel/R/confint.R
pkg/mboostDevel/man/methods.Rd
pkg/mboostDevel/tests/regtest-inference.R
Log:
- changes in the confint functions
- added manual
Modified: pkg/mboostDevel/R/confint.R
===================================================================
--- pkg/mboostDevel/R/confint.R 2014-07-03 15:47:27 UTC (rev 777)
+++ pkg/mboostDevel/R/confint.R 2014-07-21 17:22:26 UTC (rev 778)
@@ -1,8 +1,11 @@
-confint.mboost <- function(object, B = 1000, newdata = NULL,
- B.mstop = 25, which = NULL, ...) {
+confint.mboost <- function(object, parm = NULL, level = 0.95,
+ B = 1000, B.mstop = 25, newdata = NULL,
+ which = parm, ...) {
which <- object$which(which, usedonly = FALSE)
+ if (!all(which %in% object$which(NULL, usedonly = FALSE)))
+ stop(sQuote("which"), " is wrongly specified")
## create new data and/or restructure data
newdata <- .create_newdata(object, newdata, which)
@@ -22,16 +25,23 @@
predictions[[i]] <- .predict_confint(mod, newdata = newdata,
which = which)
}
- res <- list(boot_pred = predictions, data = newdata, model = object)
+
+ ## prepare returned object
+ res <- list(level = level, boot_pred = predictions, data = newdata,
+ model = object)
+ attr(res, "which") <- which
class(res) <- "mboost.ci"
return(res)
}
-confint.glmboost <- function(object, B = 1000, B.mstop = 25,
- which = NULL, ...) {
+confint.glmboost <- function(object, parm = NULL, level = 0.95,
+ B = 1000, B.mstop = 25,
+ which = parm, ...) {
outer.folds <- cv(model.weights(object), B = B)
which <- object$which(which, usedonly = FALSE)
+ if (!all(which %in% object$which(NULL, usedonly = FALSE)))
+ stop(sQuote("which"), " is wrongly specified")
coefficients <- matrix(NA, ncol = length(which), nrow = B)
colnames(coefficients) <- names(coef(object, which = which))
@@ -47,18 +57,53 @@
}
coefficients[i, ] <- unlist(coef(mod, which = which, off2int = TRUE))
}
- res <- list(boot_coefs = coefficients, model = object)
+
+ ## prepare returned object
+ res <- list(confint = .ci_glmboost(coefficients, level = level, which = which),
+ level = level, boot_coefs = coefficients, model = object)
+ attr(res, "which") <- which
class(res) <- "glmboost.ci"
return(res)
}
-print.glmboost.ci <- function(x, which = NULL, level = 0.95) {
+.ci_glmboost <- function(coefficients, level, which = NULL) {
quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
- which <- x$model$which(which, usedonly = FALSE)
- tmp <- apply(x$boot_coefs[, which], 2, FUN = quantile, probs = quantiles)
- print(t(tmp))
+
+ tmp <- apply(coefficients, 2, FUN = quantile, probs = quantiles)
+ CI <- as.data.frame(t(tmp))[which, ]
+ return(CI)
}
+## pe = add point estimte
+print.glmboost.ci <- function(x, which = NULL, level = x$level, pe = FALSE, ...) {
+
+ if (is.null(which)) {
+ which <- attr(x, "which")
+ } else {
+ which <- x$model$which(which, usedonly = FALSE)
+ if (!all(which %in% attr(x, "which")))
+ stop(sQuote("which"), " is wrongly specified")
+ }
+
+ if (!is.null(level) && level != x$level) {
+ CI <- .ci_glmboost(x$boot_coefs, level = level, which = which)
+ } else {
+ CI <- x$confint[which, ]
+ }
+
+ if (pe) {
+ tmp <- data.frame(beta = coef(x$model, which))
+ CI <- cbind(tmp, CI)
+ }
+ if (length(which) > 1) {
+ cat("\tBootstrap Confidence Intervals\n")
+ } else {
+ cat("\tBootstrap Confidence Interval\n")
+ }
+ print(CI, ...)
+ return(invisible(x))
+}
+
## ## check for varing...
## data <- model.frame(x, which = w)[[1]]
## get_vary <- x$baselearner[[w]]$get_vary
@@ -74,8 +119,9 @@
## ## Aditionally needed: Check for multivariate base-learners (except bols)
-## check for by variable and bivariate base-learners which both need a different
-## data set for prediction
+## FIXME: check for by variable and bivariate base-learners which both need a
+## different data set for prediction
+## FIXME: what about factor variables? do we get the correct levels?
.create_newdata <- function(object, newdata = NULL, which, ...) {
if (is.null(newdata)) {
data <- newdata <- model.frame(object, which = which)
@@ -104,25 +150,6 @@
return(newdata)
}
-## .create_newdata.glmboost <- function(object, newdata, ...) {
-## if (is.null(newdata)) {
-## data <- model.frame(object)
-## ## make grid!
-## tmp <- data[rep(1, 100), ]
-## grid <- function(x) {
-## if (is.numeric(x)) {
-## return(seq(min(x), max(x), length = 100))
-## } else {
-## return(rep(levels(x), length.out = 100))
-## }
-## }
-## for (j in 1:ncol(data))
-## tmp[, colnames(data)[j]] <- grid(data[,j])
-## newdata <- tmp
-## }
-## return(newdata)
-## }
-
## special prediction function for the construction of confidence intervals:
.predict_confint <- function(object, newdata = NULL, which, ...) {
predictions <- matrix(NA, ncol = length(which), nrow = nrow(newdata[[1]]))
@@ -132,28 +159,28 @@
return(predictions)
}
-# .predict_confint.glmboost <- function(object, newdata, which, ...) {
-# warning("shouldn't we return confints for coef?")
-# predict(object, newdata = newdata, which = which)
-# }
-plot.mboost.ci <- function(x, which, level = 0.95,
+### plot functions
+plot.mboost.ci <- function(x, which, level = x$level,
ylim = NULL, type = "l", col = "black",
- ci.col = "grey", raw = FALSE, ...) {
+ ci.col = rgb(170, 170, 170, alpha = 85,
+ maxColorValue = 255),
+ raw = FALSE, ...) {
- which <- x$model$which(which, usedonly = FALSE)
+ if (missing(which)) {
+ which <- attr(x, "which")
+ } else {
+ which <- x$model$which(which, usedonly = FALSE)
+ if (!all(which %in% attr(x, "which")))
+ stop(sQuote("which"), " is wrongly specified")
+ }
+ if (length(which) > 1)
+ stop("Specify a single base-learner using ", sQuote("which"))
+ CI <- .ci_mboost(x$boot_pred, level = level, which = which, raw = raw)
+
if (is.null(ylim)) {
- preds <- sapply(x$boot_pred, function(p) p[, which])
- if (!raw) {
- quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
- tmp <- apply(preds, 1, FUN = quantile, probs = quantiles)
- if (is.null(ylim))
- ylim <- range(tmp)
- } else {
- if (is.null(ylim))
- ylim <- range(preds)
- }
+ ylim <- range(CI)
}
plot(x$model, which = which, type = "n", ylim = ylim,
@@ -162,22 +189,47 @@
lines(x$model, which = which, type = "l", col = col, ...)
}
-lines.mboost.ci <- function(x, which, level = 0.95, col = "grey",
+lines.mboost.ci <- function(x, which, level = x$level,
+ col = rgb(170, 170, 170, alpha = 85,
+ maxColorValue = 255),
raw = FALSE, ...) {
- preds <- sapply(x$boot_pred, function(p) p[, which])
+
+ if (missing(which)) {
+ which <- attr(x, "which")
+ } else {
+ which <- x$model$which(which, usedonly = FALSE)
+ if (!all(which %in% attr(x, "which")))
+ stop(sQuote("which"), " is wrongly specified")
+ }
+ if (length(which) > 1)
+ stop("Specify a single base-learner using ", sQuote("which"))
+
+
+ CI <- .ci_mboost(x$boot_pred, level = level, which = which, raw = raw)
+
x.data <- x$data[[which]]
if (ncol(x.data) > 1) {
stop("Cannot plot lines for more than 1 dimenstion")
} else {
x.data <- x.data[, 1]
}
+
if (!raw) {
- quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
- tmp <- apply(preds, 1, FUN = quantile, probs = quantiles)
polygon(c(x.data, rev(x.data)),
- c(tmp[1, ], rev(tmp[2,])),
+ c(CI[1, ], rev(CI[2,])),
col = col, border = col)
} else {
- matlines(x$data[[which]], preds, col = col, lty = "solid", ...)
+ matlines(x$data[[which]], CI, col = col, lty = "solid", ...)
}
}
+
+.ci_mboost <- function(predictions, level, which = NULL, raw = FALSE) {
+
+ preds <- sapply(predictions, function(p) p[, which])
+ if (!raw) {
+ quantiles <- c((1 - level)/2, 1 - (1 - level)/2)
+ preds <- apply(preds, 1, FUN = quantile, probs = quantiles)
+ }
+
+ return(preds)
+}
Added: pkg/mboostDevel/man/confint.Rd
===================================================================
--- pkg/mboostDevel/man/confint.Rd (rev 0)
+++ pkg/mboostDevel/man/confint.Rd 2014-07-21 17:22:26 UTC (rev 778)
@@ -0,0 +1,143 @@
+\name{confint.mboost}
+
+\alias{confint.mboost}
+\alias{confint.glmboost}
+
+\alias{plot.mboost.ci}
+\alias{lines.mboost.ci}
+\alias{print.glmboost.ci}
+
+\title{
+ Pointwise Bootstrap Confidence Intervals
+}
+\description{
+ Compute and display pointwise confidence intervals
+}
+\usage{
+\method{confint}{mboost}(object, parm = NULL, level = 0.95, B = 1000,
+ B.mstop = 25, newdata = NULL, which = parm, ...)
+\method{plot}{mboost.ci}(x, which, level = x$level, ylim = NULL, type = "l", col = "black",
+ ci.col = rgb(170, 170, 170, alpha = 85, maxColorValue = 255),
+ raw = FALSE, ...)
+\method{lines}{mboost.ci}(x, which, level = x$level,
+ col = rgb(170, 170, 170, alpha = 85, maxColorValue = 255),
+ raw = FALSE, ...)
+
+
+\method{confint}{glmboost}(object, parm = NULL, level = 0.95,
+ B = 1000, B.mstop = 25, which = parm, ...)
+\method{print}{glmboost.ci}(x, which = NULL, level = x$level, pe = FALSE, ...)
+}
+
+\arguments{
+ \item{object}{
+ a fitted model object of class \code{glmboost}, \code{gamboost} or
+ \code{mboost} for which the confidence intervals should be computed.
+ }
+ \item{parm, which}{
+ a subset of base-learners to take into account for computing
+ confidence intervals. See \code{\link{mboost_methods}} for details.
+ \code{parm} is just a synonyme for \code{which} to be in line with
+ the generic \code{confint} function. Preferably use \code{which}.
+ }
+ \item{level}{
+ the confidence level required.
+ }
+ \item{B}{
+ number of outer bootstrap replicates used to compute the empirical
+ bootstrap confidence intervals.
+ }
+ \item{newdata}{
+ optionally, a data frame on which to compute the predictions for the
+ confidence intervals.
+ }
+ \item{B.mstop}{
+ number of inner bootstrap replicates used to determine the optimal
+ mstop on each of the \code{B} bootstrap samples.
+ }
+ \item{x}{
+ a confidence interval object.
+ }
+ \item{ylim}{
+ limits of the y scale. Per default computed from the data to plot.
+ }
+ \item{type}{
+ type of graphic for the point estimate, i.e., the predicted function.
+methods }
+ \item{col}{
+ color of the point estimate, i.e., the predicted function.
+ }
+ \item{ci.col}{
+ color of the confidence interval.
+ }
+ \item{raw}{
+ logical, should the raw function estimates or the derived confidence
+ estimates be plotted?
+ }
+ \item{pe}{
+ logical, should the point estimtate (PE) be also returned?
+ }
+ \item{\dots}{
+ additional arguments
+ }
+}
+\details{
+ Use a nested boostrap approach to compute pointwise confidence
+ intervals for the predicted partial functions or regression
+ parameters.
+}
+\value{
+ An object of class \code{glmboost.ci} or \code{mboost.ci} with special
+ \code{print} and/or \code{plot} functions.
+}
+\references{
+ %% ~put references to the literature/web site here ~
+}
+\author{
+ Benjamin Hofner <benjamin.hofner at fau.de>
+}
+\seealso{
+ \code{\link{cvrisk}} for crossvalidation approaches and
+ \code{\link{mboost_methods}} for other methods.
+}
+\examples{
+### a simple linear example
+set.seed(1907)
+data <- data.frame(x1 = rnorm(100), x2 = rnorm(100),
+ z = factor(sample(1:3, 100, replace = TRUE)))
+data$y <- rnorm(100, mean = data$x1 - data$x2 - 1 * (data$z == 2) +
+ 1 * (data$z == 3), sd = 0.1)
+linmod <- glmboost(y ~ x1 + x2 + z, data = data,
+ control = boost_control(mstop = 200))
+
+## compute confidence interval from 10 samples. Usually one should use
+## at least 1000 samples.
+CI <- confint(linmod, B = 10, level = 0.9)
+CI
+
+## to compute a confidence interval for another level simply change the
+## level in the print function:
+print(CI, level = 0.8)
+## or print a subset (with point estimates):
+print(CI, level = 0.8, pe = TRUE, which = "z")
+
+
+### a simple smooth example
+set.seed(1907)
+data <- data.frame(x1 = rnorm(100), x2 = rnorm(100))
+data$y <- rnorm(100, mean = data$x1^2 - sin(data$x2), sd = 0.1)
+gam <- gamboost(y ~ x1 + x2, data = data,
+ control = boost_control(mstop = 200))
+
+## compute confidence interval from 10 samples. Usually one should use
+## at least 1000 samples.
+CI_gam <- confint(gam, B = 10, level = 0.9)
+
+par(mfrow = c(1, 2))
+plot(CI_gam, which = 1)
+plot(CI_gam, which = 2)
+## to compute a confidence interval for another level simply change the
+## level in the plot or lines function:
+lines(CI_gam, which = 2, level = 0.8)
+}
+\keyword{methods}
Modified: pkg/mboostDevel/man/methods.Rd
===================================================================
--- pkg/mboostDevel/man/methods.Rd 2014-07-03 15:47:27 UTC (rev 777)
+++ pkg/mboostDevel/man/methods.Rd 2014-07-21 17:22:26 UTC (rev 778)
@@ -1,4 +1,5 @@
\name{methods}
+\alias{mboost_methods}
\alias{print.glmboost}
\alias{print.mboost}
Modified: pkg/mboostDevel/tests/regtest-inference.R
===================================================================
--- pkg/mboostDevel/tests/regtest-inference.R 2014-07-03 15:47:27 UTC (rev 777)
+++ pkg/mboostDevel/tests/regtest-inference.R 2014-07-21 17:22:26 UTC (rev 778)
@@ -231,10 +231,11 @@
refit <- glm$update(weights = model.weights(glm), risk = "inbag")
stopifnot(all.equal(coef(refit), coef(glm)))
+glm[200]
confint.glm <- confint(glm, B = 100, B.mstop = 2)
confint.glm
-confint.gam <- confint(gam, B = 100, B.mstop = 2)
+confint.gam <- confint(gam, B = 100, B.mstop = 1)
plot(confint.gam, which = 1)
plot(confint.gam, which = 2)
plot(confint.gam, which = 3)
More information about the Mboost-commits
mailing list