bugfixes, documentation improvements, and generic function agreement

9 files changed, 145 insertions, 154 deletions

diff --git a/R/coef.sparsestep.R b/R/coef.sparsestep.R
index 3a6286c..713bcb2 100644
--- a/R/coef.sparsestep.R
+++ b/R/coef.sparsestep.R
@@ -2,7 +2,8 @@
 #'
 #' @description Returns the coefficients of the SparseStep model.
 #'
-#' @param obj a "sparsestep" object
+#' @param object a "sparsestep" object
+#' @param ... further argument are ignored
 #'
 #' @return The coefficients of the SparseStep model (i.e. the betas). If the
 #' model was fitted with an intercept this will be the first value in the
@@ -17,12 +18,12 @@
 #' fit <- sparsestep(x, y)
 #' coef(fit)
 #'
-coef.sparsestep <- function(obj, ...)
+coef.sparsestep <- function(object, ...)
 {
-  if (obj$intercept) {
-    beta <- rbind(obj$a0, obj$beta)
+  if (object$intercept) {
+    beta <- rbind(object$a0, object$beta)
   } else {
-    beta <- obj$beta
+    beta <- object$beta
   }
   nbeta <- drop0(Matrix(as.matrix(beta)))
 
diff --git a/R/fit.sparsestep.R b/R/fit.sparsestep.R
deleted file mode 100644
index 020f29b..0000000
--- a/R/fit.sparsestep.R
+++ /dev/null
@@ -1,92 +0,0 @@
-#' @title Fits the SparseStep model
-#'
-#' @description Fits the SparseStep model for a single value of the 
-#' regularization parameter.
-#'
-#' @param x matrix of predictors
-#' @param y response
-#' @param lambda regularization parameter
-#' @param gamma0 starting value of the gamma parameter
-#' @param gammastop stopping value of the gamma parameter
-#' @param IMsteps number of steps of the majorization algorithm to perform for 
-#' each value of gamma
-#' @param gammastep factor to decrease gamma with at each step
-#' @param normalize if TRUE, each variable is standardized to have unit L2 
-#' norm, otherwise it is left alone.
-#' @param intercept if TRUE, an intercept is included in the model (and not 
-#' penalized), otherwise no intercept is included
-#' @param force.zero if TRUE, absolute coefficients smaller than the provided 
-#' threshold value are set to absolute zero as a post-processing step, 
-#' otherwise no thresholding is performed
-#' @param threshold threshold value to use for setting coefficients to 
-#' absolute zero
-#' @param XX The X'X matrix; useful for repeated runs where X'X stays the same
-#' @param Xy The X'y matrix; useful for repeated runs where X'y stays the same
-#' @param use.XX whether or not to compute X'X and return it
-#' @param use.Xy whether or not to compute X'y and return it
-#'
-#' @return A "sparsestep" S3 object is returned, for which print, predict, 
-#' coef, and plot methods exist. It has the following items:
-#' \item{call}{The call that was used to construct the model.}
-#' \item{lambda}{The value(s) of lambda used to construct the model.}
-#' \item{gamma0}{The gamma0 value of the model.}
-#' \item{gammastop}{The gammastop value of the model}
-#' \item{IMsteps}{The IMsteps value of the model}
-#' \item{gammastep}{The gammastep value of the model}
-#' \item{intercept}{Boolean indicating if an intercept was fitted in the 
-#' model}
-#' \item{force.zero}{Boolean indicating if a force zero-setting was 
-#' performed.}
-#' \item{threshold}{The threshold used for a forced zero-setting}
-#' \item{beta}{The resulting coefficients stored in a sparse matrix format 
-#' (dgCMatrix). This matrix has dimensions nvar x nlambda}
-#' \item{a0}{The intercept vector for each value of gamma of length nlambda}
-#' \item{normx}{Vector used to normalize the columns of x}
-#' \item{meanx}{Vector of column means of x}
-#' \item{XX}{The matrix X'X if use.XX was set to TRUE}
-#' \item{Xy}{The matrix X'y if use.Xy was set to TRUE}
-#'
-#' @author
-#' Gertjan van den Burg (author and maintainer).
-#'
-#' @seealso
-#' \code{\link{coef}}, \code{\link{print}}, \code{\link{predict}}, 
-#' \code{\link{plot}}, and \code{\link{sparsestep.path}}.
-#'
-#' @export
-#'
-#' @examples
-#' x <- matrix(rnorm(100*20), 100, 20)
-#' y <- rnorm(100)
-#' fit <- sparsestep(x, y)
-#'
-sparsestep <- function(x, y, lambda=c(0.1, 0.5, 1.0, 5, 10), gamma0=1e3,
-		       gammastop=1e-4, IMsteps=2, gammastep=2.0,
-		       normalize=TRUE, intercept=TRUE, force.zero=TRUE,
-		       threshold=1e-7, XX=NULL, Xy=NULL, use.XX = TRUE,
-		       use.Xy = TRUE)
-{
-	call <- match.call()
-	prep <- preprocess(x, y, normalize, intercept, XX, Xy, use.XX, use.Xy)
-
-	betas <- NULL
-	for (lmd in lambda) {
-		beta <- run.sparsestep(prep$XX, prep$Xy, prep$nvars, lmd,
-				       gamma0, gammastep, gammastop, IMsteps,
-				       force.zero, threshold)
-		betas <- cbind(beta, betas)
-	}
-
-	post <- postprocess(t(betas), prep$a0, x, prep$normx, prep$nvars,
-			    length(lambda))
-
-	object <- list(call = call, lambda = lambda, gamma0 = gamma0,
-		       gammastop = gammastop, IMsteps = IMsteps,
-		       gammastep = gammastep, intercept = intercept,
-		       force.zero = force.zero, threshold = threshold,
-		       beta = post$beta, a0 = post$a0, normx = prep$normx,
-		       meanx = prep$meanx, XX = if(use.XX) prep$XX else NULL,
-		       Xy = if(use.Xy) prep$Xy else NULL)
-	class(object) <- "sparsestep"
-	return(object)
-}
diff --git a/R/path.sparsestep.R b/R/path.sparsestep.R
index 2081346..6f013e0 100644
--- a/R/path.sparsestep.R
+++ b/R/path.sparsestep.R
@@ -57,14 +57,11 @@
 #' \item{XX}{The matrix X'X if use.XX was set to TRUE}
 #' \item{Xy}{The matrix X'y if use.Xy was set to TRUE}
 #'
-#' @author
-#' Gertjan van den Burg (author and maintainer).
-#'
 #' @export
 #'
 #' @seealso
 #' \code{\link{coef}}, \code{\link{print}}, \code{\link{predict}}, 
-#' \code{\link{plot}}, and \code{\link{sparsestep.path}}.
+#' \code{\link{plot}}, and \code{\link{sparsestep}}.
 #'
 #' @examples
 #' x <- matrix(rnorm(100*20), 100, 20)
@@ -90,9 +87,10 @@ path.sparsestep <- function(x, y, max.depth=10, gamma0=1e3, gammastop=1e-4,
 	beta <- NULL
 	while (1) {
 		last.beta <- beta
-		beta <- run.sparsestep(XX, Xy, nvars, lambda.max, gamma0,
-				      gammastep, gammastop, IMsteps,
-				      force.zero, threshold)
+		beta <- run.sparsestep(prep$x, prep$y, XX, Xy, nvars,
+				       lambda.max, gamma0, gammastep,
+				       gammastop, IMsteps, force.zero,
+				       threshold)
 		iter <- iter + 1
 		if (all(beta == 0)) {
 			lambda.max <- lambda.max / 2
@@ -114,9 +112,10 @@ path.sparsestep <- function(x, y, max.depth=10, gamma0=1e3, gammastop=1e-4,
 	beta <- NULL
 	while (1) {
 		last.beta <- beta
-		beta <- run.sparsestep(XX, Xy, nvars, lambda.min, gamma0,
-				      gammastep, gammastop, IMsteps,
-				      force.zero, threshold)
+		beta <- run.sparsestep(prep$x, prep$y, XX, Xy, nvars,
+				       lambda.min, gamma0, gammastep,
+				       gammastop, IMsteps, force.zero,
+				       threshold)
 		iter <- iter + 1
 		if (all(beta != 0)) {
 			lambda.min <- lambda.min * 2
@@ -141,9 +140,9 @@ path.sparsestep <- function(x, y, max.depth=10, gamma0=1e3, gammastop=1e-4,
 	left <- log(lambda.min)/log(2)
 	right <- log(lambda.max)/log(2)
 
-	l <- lambda.search(0, max.depth, have.zeros, left, right, 1, 
-			   nvars+1, XX, Xy, gamma0, gammastep, gammastop,
-			   IMsteps, force.zero, threshold)
+	l <- lambda.search(prep$x, prep$y, 0, max.depth, have.zeros, left,
+			   right, 1, nvars+1, XX, Xy, gamma0, gammastep,
+			   gammastop, IMsteps, force.zero, threshold)
 	have.zeros <- have.zeros | l$have.zeros
 	lambdas <- c(lambda.min, l$lambdas, lambda.max)
 	betas <- t(cbind(betas.min, l$betas, betas.max))
@@ -153,7 +152,7 @@ path.sparsestep <- function(x, y, max.depth=10, gamma0=1e3, gammastop=1e-4,
 	lambdas <- lambdas[ord]
 	betas <- betas[ord, ]
 
-	post <- postprocess(betas, prep$a0, x, prep$normx, nvars,
+	post <- postprocess(betas, prep$a0, prep$x, prep$normx, nvars,
 			    length(lambdas))
 
 	object <- list(call = call, lambda = lambdas, gamma0 = gamma0,
@@ -167,7 +166,7 @@ path.sparsestep <- function(x, y, max.depth=10, gamma0=1e3, gammastop=1e-4,
 	return(object)
 }
 
-lambda.search <- function(depth, max.depth, have.zeros, left, right,
+lambda.search <- function(x, y, depth, max.depth, have.zeros, left, right,
 			  lidx, ridx, XX, Xy, gamma0, gammastep, gammastop,
 			   IMsteps, force.zero, threshold)
 {
@@ -179,7 +178,7 @@ lambda.search <- function(depth, max.depth, have.zeros, left, right,
 
 	middle <- left + (right - left)/2
 	lambda <- 2^middle
-	beta <- run.sparsestep(XX, Xy, nvars, lambda, gamma0, gammastep,
+	beta <- run.sparsestep(x, y, XX, Xy, nvars, lambda, gamma0, gammastep,
 			       gammastop, IMsteps, force.zero, threshold)
 	iter <- 1
 
@@ -199,10 +198,10 @@ lambda.search <- function(depth, max.depth, have.zeros, left, right,
 		b1 <- bnd[r, 1]
 		b2 <- bnd[r, 2]
 		if (depth < max.depth && any(have.zeros[i1:i2] == F)) {
-			ds <- lambda.search(depth+1, max.depth, have.zeros,
-					    b1, b2, i1, i2, XX, Xy, gamma0,
-					    gammastep, gammastop, IMsteps,
-					    force.zero, threshold)
+			ds <- lambda.search(x, y, depth+1, max.depth,
+					    have.zeros, b1, b2, i1, i2, XX,
+					    Xy, gamma0, gammastep, gammastop,
+					    IMsteps, force.zero, threshold)
 			have.zeros <- have.zeros | ds$have.zeros
 			lambdas <- c(lambdas, ds$lambdas)
 			betas <- cbind(betas, ds$betas)
diff --git a/R/plot.sparsestep.R b/R/plot.sparsestep.R
index 87ab661..336a36b 100644
--- a/R/plot.sparsestep.R
+++ b/R/plot.sparsestep.R
@@ -2,29 +2,27 @@
 #'
 #' @description Plot the coefficients of the SparseStep path
 #'
-#' @param obj a \code{sparsestep} object
-#' @param type the plot type, default: "s".
-#' @param lty line type, default: 1
+#' @param x a \code{sparsestep} object
 #' @param ... further argument to matplot
 #'
-#' @author
-#' Gertjan van den Burg (author and maintainer).
-#'
 #' @export
 #' @aliases plot
 #'
 #' @examples
-#' data(diabetes)
-#' attach(diabetes)
+#' x <- matrix(rnorm(100*20), 100, 20)
+#' y <- rnorm(100)
 #' fit <- sparsestep(x, y)
 #' plot(fit)
-#' pth <- sparsestep.path(x, y)
+#' pth <- path.sparsestep(x, y)
 #' plot(pth)
-plot.sparsestep <- function(obj, type="s", lty=1, ...)
+plot.sparsestep <- function(x, ...)
 {
-	index <- log(obj$lambda)
-	coefs <- t(as.matrix(obj$beta))
+	index <- log(x$lambda)
+	coefs <- t(as.matrix(x$beta))
 
+	dotlist <- list(...)
+	type <- ifelse(is.null(dotlist$type), "s", dotlist$type)
+	lty <- ifelse(is.null(dotlist$lty), 1, dotlist$lty)
 	matplot(index, coefs, xlab="Log lambda", ylab="Coefficients",
-	       	type=type, lty=lty, ...)
+		type=type, lty=lty, ...)
 }
diff --git a/R/predict.sparsestep.R b/R/predict.sparsestep.R
index 8a4be1a..b2ebd77 100644
--- a/R/predict.sparsestep.R
+++ b/R/predict.sparsestep.R
@@ -5,7 +5,8 @@
 #'
 #' @param object Fitted "sparsestep" object
 #' @param newx Matrix of new values for `x` at which predictions are to be made.
-#' 
+#' @param ... further argument are ignored
+#'
 #' @return a matrix of numerical predictions of size nobs x nlambda
 #'
 #' @export
@@ -15,13 +16,13 @@
 #' x <- matrix(rnorm(100*20), 100, 20)
 #' y <- rnorm(100)
 #' fit <- sparsestep(x, y)
-#' yhat <- predict(fit)
+#' yhat <- predict(fit, x)
 #'
-predict.sparsestep <- function(obj, newx)
+predict.sparsestep <- function(object, newx, ...)
 {
-  yhat <- newx %*% as.matrix(obj$beta)
-  if (obj$intercept) {
-	  yhat <- yhat + rep(1, nrow(yhat)) %*% obj$a0
+  yhat <- newx %*% as.matrix(object$beta)
+  if (object$intercept) {
+	  yhat <- yhat + rep(1, nrow(yhat)) %*% object$a0
   }
   return(yhat)
 }
diff --git a/R/preprocess.R b/R/preprocess.R
index 253394e..f72c4f4 100644
--- a/R/preprocess.R
+++ b/R/preprocess.R
@@ -20,7 +20,6 @@ preprocess <- function(x, y, normalize, intercept, XX, Xy, use.XX, use.Xy)
 		normx <- sqrt(drop(one %*% (x^2)))
 		names(normx) <- NULL
 		x <- scale(x, FALSE, normx)
-		cat("Normalizing in sparsestep\n")
 	} else {
 		normx <- rep(1, nvars)
 	}
@@ -32,7 +31,7 @@ preprocess <- function(x, y, normalize, intercept, XX, Xy, use.XX, use.Xy)
 		Xy <- t(x) %*% y
 	}
 
-	out <- list(nvars = nvars, normx = normx, a0 = a0,
+	out <- list(x = x, y = y, nvars = nvars, normx = normx, a0 = a0,
 		    XX = XX, Xy = Xy, meanx = meanx)
 	return(out)
 }
diff --git a/R/print.sparsestep.R b/R/print.sparsestep.R
index 336d705..549815b 100644
--- a/R/print.sparsestep.R
+++ b/R/print.sparsestep.R
@@ -2,23 +2,23 @@
 #'
 #' @description Prints a short text of a fitted SparseStep model
 #'
-#' @param obj a "sparsestep" object to print
+#' @param x a "sparsestep" object to print
+#' @param ... further argument are ignored
 #'
 #' @export
 #'
 #' @aliases print
 #'
 #' @examples
-#' data(diabetes)
-#' attach(diabetes)
-#' object <- sparsestep(x, y)
-#' print(object)
-#' detach(diabetes)
+#' x <- matrix(rnorm(100*20), 100, 20)
+#' y <- rnorm(100)
+#' fit <- sparsestep(x, y)
+#' print(fit)
 #'
-print.sparsestep <- function(obj, ...)
+print.sparsestep <- function(x, ...)
 {
 	cat("\nCall:\n")
-	dput(obj$call)
-	cat("Lambda:", obj$lambda, "\n")
-	invisible(obj)
+	dput(x$call)
+	cat("Lambda:", x$lambda, "\n")
+	invisible(x)
 }
diff --git a/R/run.sparsestep.R b/R/run.sparsestep.R
index 43017c9..0b32fd4 100644
--- a/R/run.sparsestep.R
+++ b/R/run.sparsestep.R
@@ -1,8 +1,8 @@
 # Core SparseStep routine. This could be implemented in a low-level language 
 # in the future, if necessary.
 #
-run.sparsestep <- function(XX, Xy, nvars, lambda, gamma0, gammastep, gammastop,
-			   IMsteps, force.zero, threshold) {
+run.sparsestep <- function(x, y, XX, Xy, nvars, lambda, gamma0, gammastep,
+			   gammastop, IMsteps, force.zero, threshold) {
 	# Start solving SparseStep
 	gamma <- gamma0
 	beta <- matrix(0.0, nvars, 1)
diff --git a/R/sparsestep.R b/R/sparsestep.R
index 05cf9ff..b1f9a58 100644
--- a/R/sparsestep.R
+++ b/R/sparsestep.R
@@ -1,5 +1,90 @@
-#' sparsestep.
+#' @title Fits the SparseStep model
 #'
-#' @name sparsestep
-#' @docType package
-NULL
+#' @description Fits the SparseStep model for a single value of the 
+#' regularization parameter.
+#'
+#' @param x matrix of predictors
+#' @param y response
+#' @param lambda regularization parameter
+#' @param gamma0 starting value of the gamma parameter
+#' @param gammastop stopping value of the gamma parameter
+#' @param IMsteps number of steps of the majorization algorithm to perform for 
+#' each value of gamma
+#' @param gammastep factor to decrease gamma with at each step
+#' @param normalize if TRUE, each variable is standardized to have unit L2 
+#' norm, otherwise it is left alone.
+#' @param intercept if TRUE, an intercept is included in the model (and not 
+#' penalized), otherwise no intercept is included
+#' @param force.zero if TRUE, absolute coefficients smaller than the provided 
+#' threshold value are set to absolute zero as a post-processing step, 
+#' otherwise no thresholding is performed
+#' @param threshold threshold value to use for setting coefficients to 
+#' absolute zero
+#' @param XX The X'X matrix; useful for repeated runs where X'X stays the same
+#' @param Xy The X'y matrix; useful for repeated runs where X'y stays the same
+#' @param use.XX whether or not to compute X'X and return it
+#' @param use.Xy whether or not to compute X'y and return it
+#'
+#' @return A "sparsestep" S3 object is returned, for which print, predict, 
+#' coef, and plot methods exist. It has the following items:
+#' \item{call}{The call that was used to construct the model.}
+#' \item{lambda}{The value(s) of lambda used to construct the model.}
+#' \item{gamma0}{The gamma0 value of the model.}
+#' \item{gammastop}{The gammastop value of the model}
+#' \item{IMsteps}{The IMsteps value of the model}
+#' \item{gammastep}{The gammastep value of the model}
+#' \item{intercept}{Boolean indicating if an intercept was fitted in the 
+#' model}
+#' \item{force.zero}{Boolean indicating if a force zero-setting was 
+#' performed.}
+#' \item{threshold}{The threshold used for a forced zero-setting}
+#' \item{beta}{The resulting coefficients stored in a sparse matrix format 
+#' (dgCMatrix). This matrix has dimensions nvar x nlambda}
+#' \item{a0}{The intercept vector for each value of gamma of length nlambda}
+#' \item{normx}{Vector used to normalize the columns of x}
+#' \item{meanx}{Vector of column means of x}
+#' \item{XX}{The matrix X'X if use.XX was set to TRUE}
+#' \item{Xy}{The matrix X'y if use.Xy was set to TRUE}
+#'
+#' @seealso
+#' \code{\link{coef}}, \code{\link{print}}, \code{\link{predict}}, 
+#' \code{\link{plot}}, and \code{\link{path.sparsestep}}.
+#'
+#' @export
+#'
+#' @examples
+#' x <- matrix(rnorm(100*20), 100, 20)
+#' y <- rnorm(100)
+#' fit <- sparsestep(x, y)
+#'
+sparsestep <- function(x, y, lambda=c(0.1, 0.5, 1.0, 5, 10), gamma0=1e3,
+		       gammastop=1e-4, IMsteps=2, gammastep=2.0,
+		       normalize=TRUE, intercept=TRUE, force.zero=TRUE,
+		       threshold=1e-7, XX=NULL, Xy=NULL, use.XX = TRUE,
+		       use.Xy = TRUE)
+{
+	call <- match.call()
+	prep <- preprocess(x, y, normalize, intercept, XX, Xy, use.XX, use.Xy)
+
+	betas <- NULL
+	for (lmd in lambda) {
+		beta <- run.sparsestep(prep$x, prep$y, prep$XX, prep$Xy,
+				       prep$nvars, lmd, gamma0, gammastep,
+				       gammastop, IMsteps, force.zero,
+				       threshold)
+		betas <- cbind(betas, beta)
+	}
+
+	post <- postprocess(t(betas), prep$a0, prep$x, prep$normx, prep$nvars,
+			    length(lambda))
+
+	object <- list(call = call, lambda = lambda, gamma0 = gamma0,
+		       gammastop = gammastop, IMsteps = IMsteps,
+		       gammastep = gammastep, intercept = intercept,
+		       force.zero = force.zero, threshold = threshold,
+		       beta = post$beta, a0 = post$a0, normx = prep$normx,
+		       meanx = prep$meanx, XX = if(use.XX) prep$XX else NULL,
+		       Xy = if(use.Xy) prep$Xy else NULL)
+	class(object) <- "sparsestep"
+	return(object)
+}