From 9d68ff63ad9752452ecaccaf9a215cff2bcabac3 Mon Sep 17 00:00:00 2001 From: Patrick Breheny Date: Sun, 21 Apr 2024 12:01:42 -0500 Subject: [PATCH] Version 1.6.0 --- .version.json | 2 +- DESCRIPTION | 9 +- NEWS.md | 21 +++- R/biglasso-package.R | 47 ++++---- R/biglasso.R | 103 +++++++++-------- R/biglasso_fit.R | 78 ++++++------- R/biglasso_path.R | 80 +++++++------ R/{cv.biglasso.R => cv-biglasso.R} | 106 +++++++++--------- R/loss.R | 20 ++-- R/{plot.biglasso.R => plot-biglasso.R} | 20 ++-- R/plot-cv-biglasso.R | 42 +++++++ R/{plot.mbiglasso.R => plot-mbiglasso.R} | 17 ++- R/plot.cv.biglasso.R | 44 -------- R/{predict.cv.R => predict-cv.R} | 51 ++++----- R/predict.R | 56 ++++----- R/setupX.R | 48 ++++---- ...ry.cv.biglasso.R => summary-cv-biglasso.R} | 45 ++++---- README.md | 43 +++---- inst/CITATION | 17 ++- man/biglasso-package.Rd | 43 ++++--- man/biglasso.Rd | 49 ++++---- man/biglasso_fit.Rd | 93 ++++++++------- man/biglasso_path.Rd | 99 ++++++++-------- man/cv.biglasso.Rd | 47 ++++---- man/loss.biglasso.Rd | 4 +- man/plot.biglasso.Rd | 17 +-- man/plot.cv.biglasso.Rd | 14 +-- man/plot.mbiglasso.Rd | 14 +-- man/predict.biglasso.Rd | 31 +++-- man/predict.cv.biglasso.Rd | 32 +++--- man/setupX.Rd | 17 +-- man/summary.cv.biglasso.Rd | 35 +++--- src/Makevars | 2 - vignettes/biglasso.Rmd | 30 +++-- 34 files changed, 651 insertions(+), 725 deletions(-) rename R/{cv.biglasso.R => cv-biglasso.R} (66%) rename R/{plot.biglasso.R => plot-biglasso.R} (80%) create mode 100644 R/plot-cv-biglasso.R rename R/{plot.mbiglasso.R => plot-mbiglasso.R} (90%) delete mode 100644 R/plot.cv.biglasso.R rename R/{predict.cv.R => predict-cv.R} (54%) rename R/{summary.cv.biglasso.R => summary-cv-biglasso.R} (64%) diff --git a/.version.json b/.version.json index 87b449d..323e140 100644 --- a/.version.json +++ b/.version.json @@ -1,6 +1,6 @@ { "schemaVersion": 1, "label": "GitHub", - "message": "1.5.2.1", + "message": "1.6.0", "color": "blue" } diff --git a/DESCRIPTION b/DESCRIPTION index 4863b37..5b7cc03 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: biglasso -Version: 1.5.2.1 -Date: 2024-03-19 +Version: 1.6.0 +Date: 2024-04-21 Title: Extending Lasso Model Fitting to Big Data Authors@R: c( person("Yaohui", "Zeng", role = c("aut")), @@ -15,8 +15,8 @@ Description: Extend lasso and elastic-net model fitting for ultra lasso-fitting packages like 'glmnet' and 'ncvreg', thus allowing the user to analyze big data analysis even on an ordinary laptop. License: GPL-3 -URL: https://yaohuizeng.github.io/biglasso/index.html, https://github.com/YaohuiZeng/biglasso, https://arxiv.org/abs/1701.05936 -BugReports: https://github.com/YaohuiZeng/biglasso/issues +URL: https://pbreheny.github.io/biglasso/index.html, https://github.com/pbreheny/biglasso, https://arxiv.org/abs/1701.05936 +BugReports: https://github.com/pbreheny/biglasso/issues Depends: R (>= 3.2.0), bigmemory (>= 4.5.0), Matrix, ncvreg Imports: Rcpp (>= 0.12.1), methods LinkingTo: Rcpp, RcppArmadillo (>= 0.8.600), bigmemory, BH @@ -28,5 +28,6 @@ Suggests: survival, knitr, rmarkdown +Roxygen: list(markdown = TRUE) RoxygenNote: 7.3.1 Encoding: UTF-8 diff --git a/NEWS.md b/NEWS.md index 1b9c6b8..863a151 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,10 +1,15 @@ +# biglasso 1.6.0 + * New: functions biglasso_fit() and biglasso_path(), which allow users to turn + off standardization and intercept + # biglasso 1.5.2 * Update coercion for compatibility with Matrix 1.5 * Now using GitHub Actions instead of Travis for CI # biglasso 1.5.1 * Internal Cpp changes: initialize Xty, remove unused cutoff variable (#48) - * Eliminate CV test against ncvreg (the two packages no longer use the same approach (#47) + * Eliminate CV test against ncvreg (the two packages no longer use the same + approach (#47) # biglasso 1.5.0 * Update headers to maintain compatibility with new version of Rcpp (#40) @@ -13,14 +18,17 @@ * changed R package maintainer to Chuyi Wang (wwaa0208@gmail.com) * fixed bugs * Add 'auc', 'class' options to cv.biglasso eval.metric - * predict.cv now predicts standard error over CV folds by default; set 'grouped' argument to FALSE for old behaviour. - * predict.cv.biglasso accepts 'lambda.min', 'lambda.1se' argument, similar to predict.cv.glmnet() + * predict.cv now predicts standard error over CV folds by default; set + 'grouped' argument to FALSE for old behaviour. + * predict.cv.biglasso accepts 'lambda.min', 'lambda.1se' argument, similar to + predict.cv.glmnet() # biglasso 1.4-0 * adaptive screening methods were implemented and set as default when applicable * added sparse Cox regression - * removed uncompetitive screening methods and combined naming of screening methods - * version 1.4-0 for CRAN submission + * removed uncompetitive screening methods and combined naming of screening + methods + * version 1.4-0 for CRAN submission # biglasso 1.3-7 * update email to personal email @@ -30,7 +38,8 @@ # biglasso 1.3-6 * optimized the code for computing the slores rule. - * added Slores screening without active cycling (-NAC) for logistic regression, research usage only. + * added Slores screening without active cycling (-NAC) for logistic + regression, research usage only. * corrected BEDPP for elastic net. * fixed a bug related to "exporting SSR-BEDPP". diff --git a/R/biglasso-package.R b/R/biglasso-package.R index aa41abd..12612b3 100644 --- a/R/biglasso-package.R +++ b/R/biglasso-package.R @@ -26,7 +26,7 @@ #' Data in R. Version >= 1.2-3 represents a major redesign where the source #' code is converted into C++ (previously in C), and new feature screening #' rules, as well as OpenMP parallel computing, are implemented. Some key -#' features of \code{biglasso} are summarized as below: \enumerate{ \item it +#' features of `biglasso` are summarized as below: \enumerate{ \item it #' utilizes memory-mapped files to store the massive data on the disk, only #' loading data into memory when necessary during model fitting. Consequently, #' it's able to seamlessly data-larger-than-RAM cases. \item it is built upon @@ -38,57 +38,54 @@ #' additional 1.5x to 4x speedup. \item the implementation is designed to be as #' memory-efficient as possible by eliminating extra copies of the data created #' by other R packages, making it at least 2x more memory-efficient than -#' \code{glmnet}. \item the underlying computation is implemented in C++, and +#' `glmnet`. \item the underlying computation is implemented in C++, and #' parallel computing with OpenMP is also supported. } #' -#' \strong{For more information:} \itemize{ \item Benchmarking results: -#' \url{https://github.com/YaohuiZeng/biglasso}. -#' \item Tutorial: -#' \url{http://yaohuizeng.github.io/biglasso/articles/biglasso.html} -#' \item Technical paper: -#' \url{https://arxiv.org/abs/1701.05936} } +#' **For more information:** +#' * Benchmarking results: \url{https://github.com/pbreheny/biglasso} +#' * Tutorial: \url{https://pbreheny.github.io/biglasso/articles/biglasso.html} +#' * Technical paper: \url{https://arxiv.org/abs/1701.05936} #' #' @name biglasso-package #' -#' @note The input design matrix X must be a \code{\link[bigmemory]{big.matrix}} object. -#' This can be created by the function \code{as.big.matrix} in the R package +#' @note The input design matrix X must be a [bigmemory::big.matrix()] object. +#' This can be created by the function `as.big.matrix` in the R package #' \href{https://CRAN.R-project.org//package=bigmemory}{bigmemory}. #' If the data (design matrix) is very large (e.g. 10 GB) and stored in an external #' file, which is often the case for big data, X can be created by calling the -#' function \code{\link{setupX}}. +#' function [setupX()]. #' \strong{In this case, there are several restrictions about the data file:} #' \enumerate{ \item the data file must be a well-formated ASCII-file, with #' each row corresponding to an observation and each column a variable; \item #' the data file must contain only one single type. Current version only -#' supports \code{double} type; \item the data file must contain only numeric +#' supports `double` type; \item the data file must contain only numeric #' variables. If there are categorical variables, the user needs to create #' dummy variables for each categorical varable (by adding additional columns).} #' Future versions will try to address these restrictions. #' -#' Denote the number of observations and variables be, respectively, \code{n} -#' and \code{p}. It's worth noting that the package is more suitable for wide -#' data (ultrahigh-dimensional, \code{p >> n}) as compared to long data -#' (\code{n >> p}). This is because the model fitting algorithm takes advantage +#' Denote the number of observations and variables be, respectively, `n` +#' and `p`. It's worth noting that the package is more suitable for wide +#' data (ultrahigh-dimensional, `p >> n`) as compared to long data +#' (`n >> p`). This is because the model fitting algorithm takes advantage #' of sparsity assumption of high-dimensional data. To just give the user some #' ideas, below are some benchmarking results of the total computing time (in #' seconds) for solving lasso-penalized linear regression along a sequence of #' 100 values of the tuning parameter. In all cases, assume 20 non-zero #' coefficients equal +/- 2 in the true model. (Based on Version 1.2-3, #' screening rule "SSR-BEDPP" is used) -#' \itemize{ \item For wide data case (\code{p > n}), \code{n = 1,000}: -#' \tabular{ccccc}{ \code{p} \tab 1,000 \tab 10,000 \tab 100,000 \tab 1,000,000 -#' \cr Size of \code{X} \tab 9.5 MB \tab 95 MB \tab 950 MB \tab 9.5 GB \cr +#' \itemize{ \item For wide data case (`p > n`), `n = 1,000`: +#' \tabular{ccccc}{ `p` \tab 1,000 \tab 10,000 \tab 100,000 \tab 1,000,000 +#' \cr Size of `X` \tab 9.5 MB \tab 95 MB \tab 950 MB \tab 9.5 GB \cr #' Elapsed time (s) \tab 0.11 \tab 0.83 \tab 8.47 \tab 85.50 \cr } -#' %\item For long data case (\code{n >> p}), \code{p = 1,000}: +#' %\item For long data case (`n >> p`), `p = 1,000`: #' %\tabular{ccccc}{ -#' %\code{n} \tab 1,000 \tab 10,000 \tab 100,000 \tab 1,000,000 \cr -#' %Size of \code{X} \tab 9.5 MB \tab 95 MB \tab 950 MB \tab 9.5 GB \cr +#' %`n` \tab 1,000 \tab 10,000 \tab 100,000 \tab 1,000,000 \cr +#' %Size of `X` \tab 9.5 MB \tab 95 MB \tab 950 MB \tab 9.5 GB \cr #' %Elapsed time (s) \tab 2.50 \tab 11.43 \tab 83.69 \tab 1090.62 \cr %} #' } #' -#' @author Yaohui Zeng, Chuyi Wang and Patrick Breheny +#' @author Yaohui Zeng, Chuyi Wang, Tabitha Peter, and Patrick Breheny #' -#' Maintainer: Yaohui Zeng and Chuyi Wang #' @references \itemize{ \item Zeng, Y., and Breheny, P. (2017). The biglasso #' Package: A Memory- and Computation-Efficient Solver for Lasso Model Fitting #' with Big Data in R. \url{https://arxiv.org/abs/1701.05936}. \item @@ -104,7 +101,7 @@ #' 2137-2140). IEEE. \item Wang, J., Zhou, J., Liu, J., Wonka, P., and Ye, J. #' (2014). A safe screening rule for sparse logistic regression. \emph{In #' Advances in Neural Information Processing Systems}, pp. 1053-1061. } -#' @keywords package +#' #' @examples #' \dontrun{ #' ## Example of reading data from external big data file, fit lasso model, diff --git a/R/biglasso.R b/R/biglasso.R index af0e89d..c63022e 100644 --- a/R/biglasso.R +++ b/R/biglasso.R @@ -6,52 +6,52 @@ #' parameter lambda. #' #' The objective function for linear regression or multiple responses linear regression -#' (\code{family = "gaussian"} or \code{family = "mgaussian"}) is +#' (`family = "gaussian"` or `family = "mgaussian"`) is #' \deqn{\frac{1}{2n}\textrm{RSS} + \lambda*\textrm{penalty},}{(1/(2n))*RSS+ #' \lambda*penalty,} -#' where for \code{family = "mgaussian"}), a group-lasso type penalty is applied. +#' where for `family = "mgaussian"`), a group-lasso type penalty is applied. #' For logistic regression -#' (\code{family = "binomial"}) it is \deqn{-\frac{1}{n} loglike + +#' (`family = "binomial"`) it is \deqn{-\frac{1}{n} loglike + #' \lambda*\textrm{penalty},}{-(1/n)*loglike+\lambda*penalty}, for cox regression, #' breslow approximation for ties is applied. #' #' Several advanced feature screening rules are implemented. For -#' lasso-penalized linear regression, all the options of \code{screen} are -#' applicable. Our proposal adaptive rule - \code{"Adaptive"} - achieves highest speedup +#' lasso-penalized linear regression, all the options of `screen` are +#' applicable. Our proposal adaptive rule - `"Adaptive"` - achieves highest speedup #' so it's the recommended one, especially for ultrahigh-dimensional large-scale #' data sets. For cox regression and/or the elastic net penalty, only -#' \code{"SSR"} is applicable for now. More efficient rules are under development. +#' `"SSR"` is applicable for now. More efficient rules are under development. #' #' @param X The design matrix, without an intercept. It must be a -#' double type \code{\link[bigmemory]{big.matrix}} object. The function +#' double type [bigmemory::big.matrix()] object. The function #' standardizes the data and includes an intercept internally by default during #' the model fitting. -#' @param y The response vector for \code{family="gaussian"} or \code{family="binomial"}. -#' For \code{family="cox"}, \code{y} should be a two-column matrix with columns +#' @param y The response vector for `family="gaussian"` or `family="binomial"`. +#' For `family="cox"`, `y` should be a two-column matrix with columns #' 'time' and 'status'. The latter is a binary variable, with '1' indicating death, -#' and '0' indicating right censored. For \code{family="mgaussin"}, \code{y} +#' and '0' indicating right censored. For `family="mgaussin"`, `y` #' should be a n*m matrix where n is the sample size and m is the number of #' responses. -#' @param row.idx The integer vector of row indices of \code{X} that used for -#' fitting the model. \code{1:nrow(X)} by default. -#' @param penalty The penalty to be applied to the model. Either \code{"lasso"} -#' (the default), \code{"ridge"}, or \code{"enet"} (elastic net). -#' @param family Either \code{"gaussian"}, \code{"binomial"}, \code{"cox"} or -#' \code{"mgaussian"} depending on the response. +#' @param row.idx The integer vector of row indices of `X` that used for +#' fitting the model. `1:nrow(X)` by default. +#' @param penalty The penalty to be applied to the model. Either `"lasso"` +#' (the default), `"ridge"`, or `"enet"` (elastic net). +#' @param family Either `"gaussian"`, `"binomial"`, `"cox"` or +#' `"mgaussian"` depending on the response. #' @param alg.logistic The algorithm used in logistic regression. If "Newton" #' then the exact hessian is used (default); if "MM" then a #' majorization-minimization algorithm is used to set an upper-bound on the #' hessian matrix. This can be faster, particularly in data-larger-than-RAM #' case. -#' @param screen The feature screening rule used at each \code{lambda} that -#' discards features to speed up computation: \code{"SSR"} (default if -#' \code{penalty="ridge"} or \code{penalty="enet"} )is the sequential strong rule; -#' \code{"Hybrid"} is our newly proposed hybrid screening rules which combine the -#' strong rule with a safe rule. \code{"Adaptive"} (default for \code{penalty="lasso"} -#' without \code{penalty.factor}) is our newly proposed adaptive rules which +#' @param screen The feature screening rule used at each `lambda` that +#' discards features to speed up computation: `"SSR"` (default if +#' `penalty="ridge"` or `penalty="enet"` )is the sequential strong rule; +#' `"Hybrid"` is our newly proposed hybrid screening rules which combine the +#' strong rule with a safe rule. `"Adaptive"` (default for `penalty="lasso"` +#' without `penalty.factor`) is our newly proposed adaptive rules which #' reuse screening reference for multiple lambda values. \strong{Note that:} -#' (1) for linear regression with elastic net penalty, both \code{"SSR"} and -#' \code{"Hybrid"} are applicable since version 1.3-0; (2) only \code{"SSR"} is +#' (1) for linear regression with elastic net penalty, both `"SSR"` and +#' `"Hybrid"` are applicable since version 1.3-0; (2) only `"SSR"` is #' applicable to elastic-net-penalized logistic regression or cox regression; #' (3) active set cycling strategy is incorporated with these screening rules. #' @param safe.thresh the threshold value between 0 and 1 that controls when to @@ -67,8 +67,8 @@ #' @param alpha The elastic-net mixing parameter that controls the relative #' contribution from the lasso (l1) and the ridge (l2) penalty. The penalty is #' defined as \deqn{ \alpha||\beta||_1 + (1-\alpha)/2||\beta||_2^2.} -#' \code{alpha=1} is the lasso penalty, \code{alpha=0} the ridge penalty, -#' \code{alpha} in between 0 and 1 is the elastic-net ("enet") penalty. +#' `alpha=1` is the lasso penalty, `alpha=0` the ridge penalty, +#' `alpha` in between 0 and 1 is the elastic-net ("enet") penalty. #' @param lambda.min The smallest value for lambda, as a fraction of #' lambda.max. Default is .001 if the number of observations is larger than #' the number of covariates and .05 otherwise. @@ -76,23 +76,23 @@ #' @param lambda.log.scale Whether compute the grid values of lambda on log #' scale (default) or linear scale. #' @param lambda A user-specified sequence of lambda values. By default, a -#' sequence of values of length \code{nlambda} is computed, equally spaced on +#' sequence of values of length `nlambda` is computed, equally spaced on #' the log scale. #' @param eps Convergence threshold for inner coordinate descent. The #' algorithm iterates until the maximum change in the objective after any -#' coefficient update is less than \code{eps} times the null deviance. Default -#' value is \code{1e-7}. +#' coefficient update is less than `eps` times the null deviance. Default +#' value is `1e-7`. #' @param max.iter Maximum number of iterations. Default is 1000. #' @param dfmax Upper bound for the number of nonzero coefficients. Default is #' no upper bound. However, for large data sets, computational burden may be #' heavy for models with a large number of nonzero coefficients. #' @param penalty.factor A multiplicative factor for the penalty applied to -#' each coefficient. If supplied, \code{penalty.factor} must be a numeric -#' vector of length equal to the number of columns of \code{X}. The purpose of -#' \code{penalty.factor} is to apply differential penalization if some +#' each coefficient. If supplied, `penalty.factor` must be a numeric +#' vector of length equal to the number of columns of `X`. The purpose of +#' `penalty.factor` is to apply differential penalization if some #' coefficients are thought to be more likely than others to be in the model. #' Current package doesn't allow unpenalized coefficients. That -#' is\code{penalty.factor} cannot be 0. \code{penalty.factor} is only supported +#' is`penalty.factor` cannot be 0. `penalty.factor` is only supported #' for "SSR" screen. #' @param warn Return warning messages for failures to converge and model #' saturation? Default is TRUE. @@ -102,45 +102,44 @@ #' fitting. Default is TRUE. #' @param verbose Whether to output the timing of each lambda iteration. #' Default is FALSE. -#' @return An object with S3 class \code{"biglasso"} for -#' \code{"gaussian", "binomial", "cox"} families, or an object with S3 class -#' \code{"mbiglasso"} for \code{"mgaussian"} family, with following variables. +#' +#' @returns An object with S3 class `"biglasso"` for +#' `"gaussian", "binomial", "cox"` families, or an object with S3 class +#' `"mbiglasso"` for `"mgaussian"` family, with following variables. #' \item{beta}{The fitted matrix of coefficients, store in sparse matrix #' representation. The number of rows is equal to the number of coefficients, -#' whereas the number of columns is equal to \code{nlambda}. For \code{"mgaussian"} +#' whereas the number of columns is equal to `nlambda`. For `"mgaussian"` #' family with m responses, it is a list of m such matrices.} -#' \item{iter}{A vector of length \code{nlambda} containing the number of -#' iterations until convergence at each value of \code{lambda}.} +#' \item{iter}{A vector of length `nlambda` containing the number of +#' iterations until convergence at each value of `lambda`.} #' \item{lambda}{The sequence of regularization parameter values in the path.} #' \item{penalty}{Same as above.} #' \item{family}{Same as above.} #' \item{alpha}{Same as above.} #' \item{loss}{A vector containing either the residual sum of squares -#' (for \code{"gaussian", "mgaussian"}) or negative log-likelihood -#' (for \code{"binomial", "cox"}) of the fitted model at each value of \code{lambda}.} +#' (for `"gaussian", "mgaussian"`) or negative log-likelihood +#' (for `"binomial", "cox"`) of the fitted model at each value of `lambda`.} #' \item{penalty.factor}{Same as above.} #' \item{n}{The number of observations used in the model fitting. It's equal to -#' \code{length(row.idx)}.} +#' `length(row.idx)`.} #' \item{center}{The sample mean vector of the variables, i.e., column mean of -#' the sub-matrix of \code{X} used for model fitting.} +#' the sub-matrix of `X` used for model fitting.} #' \item{scale}{The sample standard deviation of the variables, i.e., column -#' standard deviation of the sub-matrix of \code{X} used for model fitting.} +#' standard deviation of the sub-matrix of `X` used for model fitting.} #' \item{y}{The response vector used in the model fitting. Depending on -#' \code{row.idx}, it could be a subset of the raw input of the response vector y.} +#' `row.idx`, it could be a subset of the raw input of the response vector y.} #' \item{screen}{Same as above.} #' \item{col.idx}{The indices of features that have 'scale' value greater than #' 1e-6. Features with 'scale' less than 1e-6 are removed from model fitting.} -#' \item{rejections}{The number of features rejected at each value of \code{lambda}.} +#' \item{rejections}{The number of features rejected at each value of `lambda`.} #' \item{safe_rejections}{The number of features rejected by safe rules at each -#' value of \code{lambda}.} +#' value of `lambda`.} +#' #' @author Yaohui Zeng, Chuyi Wang and Patrick Breheny #' -#' Maintainer: Yaohui Zeng and Chuyi Wang -#' @seealso \code{\link{biglasso-package}}, \code{\link{setupX}}, -#' \code{\link{cv.biglasso}}, \code{\link{plot.biglasso}}, -#' \code{\link[ncvreg]{ncvreg}} -#' @examples +#' @seealso [biglasso-package], [setupX()], [cv.biglasso()], [plot.biglasso()], [ncvreg::ncvreg()] #' +#' @examples #' ## Linear regression #' data(colon) #' X <- colon$X diff --git a/R/biglasso_fit.R b/R/biglasso_fit.R index 4ecd239..344c2f7 100644 --- a/R/biglasso_fit.R +++ b/R/biglasso_fit.R @@ -1,30 +1,27 @@ -#' Simplified call to biglasso: a gaussian model fit with no 'bells and whistles' (e.g., no SSR) +#' Direct interface to biglasso fitting, no preprocessing #' -#' NOTE: this function is designed for users who have a strong understanding of -#' statistics and know exactly what they are doing. This is a simplification of -#' the main `biglasso()` function with more flexible settings. +#' This function is intended for users who know exactly what they're doing and +#' want complete control over the fitting process. It +#' * does NOT add an intercept +#' * does NOT standardize the design matrix +#' * does NOT set up a path for lambda (the lasso tuning parameter) +#' all of the above are critical steps in data analysis. However, a direct API +#' has been provided for use in situations where the lasso fitting process is +#' an internal component of a more complicated algorithm and standardization +#' must be handled externally. #' -#' Of note, this function: -#' -#' * does NOT add an intercept -#' * does NOT standardize the design matrix -#' * does NOT set up a path for lambda (the lasso tuning parameter) -#' -#' all of the above are among the best practices for data analysis. This function -#' is made for use in situations where these steps have already been addressed prior -#' to model fitting. -#' -#' In other words, `biglasso_fit()` is to `biglasso()` what `ncvreg::ncvfit()` -#' is to `ncvreg::ncvreg()`. -#' -#' For now, this function only works with linear regression (`family = 'gaussian'`) +#' Note: +#' * Hybrid safe-strong rules are turned off for `biglasso_fit()`, as these rely +#' on standardization +#' * Currently, the function only works with linear regression +#' (`family = 'gaussian'`). #' #' @param X The design matrix, without an intercept. It must be a -#' double type \code{\link[bigmemory]{big.matrix}} object. +#' double type [bigmemory::big.matrix()] object. #' @param y The response vector #' @param r Residuals (length n vector) corresponding to `init`. -#' WARNING: If you supply an incorrect value of `r`, the -#' solution will be incorrect. +#' WARNING: If you supply an incorrect value of `r`, the +#' solution will be incorrect. #' @param init Initial values for beta. Default: zero (length p vector) #' @param xtx X scales: the jth element should equal `crossprod(X[,j])/n`. #' In particular, if X is standardized, one should pass @@ -37,23 +34,23 @@ #' contribution from the lasso (l1) and the ridge (l2) penalty. #' The penalty is defined as: #' \deqn{ \alpha||\beta||_1 + (1-\alpha)/2||\beta||_2^2.} -#' \code{alpha=1} is the lasso penalty, \code{alpha=0} the ridge penalty, -#' \code{alpha} in between 0 and 1 is the elastic-net ("enet") penalty. +#' `alpha=1` is the lasso penalty, `alpha=0` the ridge penalty, +#' `alpha` in between 0 and 1 is the elastic-net ("enet") penalty. #' @param gamma Tuning parameter value for nonconvex penalty. Defaults are #' 3.7 for `penalty = 'SCAD'` and 3 for `penalty = 'MCP'` #' @param ncores The number of OpenMP threads used for parallel computing. #' @param max.iter Maximum number of iterations. Default is 1000. #' @param eps Convergence threshold for inner coordinate descent. The #' algorithm iterates until the maximum change in the objective -#' after any coefficient update is less than \code{eps} times -#' the null deviance. Default value is \code{1e-7}. +#' after any coefficient update is less than `eps` times +#' the null deviance. Default value is `1e-7`. #' @param dfmax Upper bound for the number of nonzero coefficients. Default is #' no upper bound. However, for large data sets, #' computational burden may be heavy for models with a large #' number of nonzero coefficients. #' @param penalty.factor A multiplicative factor for the penalty applied to -#' each coefficient. If supplied, \code{penalty.factor} must be a numeric -#' vector of length equal to the number of columns of \code{X}. +#' each coefficient. If supplied, `penalty.factor` must be a numeric +#' vector of length equal to the number of columns of `X`. #' @param warn Return warning messages for failures to converge and model #' saturation? Default is TRUE. #' @param output.time Whether to print out the start and end time of the model @@ -61,9 +58,9 @@ #' @param return.time Whether to return the computing time of the model #' fitting. Default is TRUE. #' -#' @return An object with S3 class \code{"biglasso"} with following variables. +#' @returns An object with S3 class `"biglasso"` with following variables. #' \item{beta}{The vector of estimated coefficients} -#' \item{iter}{A vector of length \code{nlambda} containing the number of +#' \item{iter}{A vector of length `nlambda` containing the number of #' iterations until convergence} #' \item{resid}{Vector of residuals calculated from estimated coefficients.} #' \item{lambda}{The sequence of regularization parameter values in the path.} @@ -72,25 +69,25 @@ #' \item{penalty.factor}{Same as in `biglasso()`.} #' \item{n}{The number of observations used in the model fitting.} #' \item{y}{The response vector used in the model fitting.} +#' #' @author Tabitha Peter and Patrick Breheny #' #' @examples -#' #' data(Prostate) -#' X <- cbind(1, Prostate$X) |> ncvreg::std() # standardizing -> xtx is all 1s +#' X <- cbind(1, Prostate$X) +#' xtx <- apply(X, 2, crossprod)/nrow(X) #' y <- Prostate$y #' X.bm <- as.big.matrix(X) -#' init <- rep(0, ncol(X)) # using cold starts - will need more iterations -#' r <- y - X%*%init -#' fit <- biglasso_fit(X = X.bm, y = y, r = r, init = init, -#' xtx = rep(1, ncol(X)),lambda = 0.1, penalty.factor=c(0, rep(1, ncol(X)-1)), -#' max.iter = 10000) -#' -#' fit <- biglasso_fit(X = X.bm, y = y, r = r, init = init, penalty = 'MCP', -#' xtx = rep(1, ncol(X)), lambda = 0.005, penalty.factor=c(0, rep(1, ncol(X)-1)), -#' max.iter = 10000) +#' init <- rep(0, ncol(X)) +#' fit <- biglasso_fit(X = X.bm, y = y, r=y, init = init, xtx = xtx, +#' lambda = 0.1, penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter = 10000) +#' fit$beta #' +#' fit <- biglasso_fit(X = X.bm, y = y, r=y, init = init, xtx = xtx, penalty='MCP', +#' lambda = 0.1, penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter = 10000) +#' fit$beta #' @export biglasso_fit + biglasso_fit <- function(X, y, r, @@ -203,4 +200,3 @@ biglasso_fit <- function(X, val <- structure(return.val, class = c("biglasso", 'ncvreg')) } - diff --git a/R/biglasso_path.R b/R/biglasso_path.R index d5b7cc6..bcc26f2 100644 --- a/R/biglasso_path.R +++ b/R/biglasso_path.R @@ -1,26 +1,25 @@ -#' Simplified call to biglasso: a gaussian model fit with no 'bells and whistles' (e.g., no SSR) +#' Direct interface to biglasso fitting, no preprocessing, path version #' -#' NOTE: this function is designed for users who have a strong understanding of -#' statistics and know exactly what they are doing. This is a simplification of -#' the main `biglasso()` function with more flexible settings. +#' This function is intended for users who know exactly what they're doing and +#' want complete control over the fitting process. It +#' * does NOT add an intercept +#' * does NOT standardize the design matrix +#' both of the above are critical steps in data analysis. However, a direct API +#' has been provided for use in situations where the lasso fitting process is +#' an internal component of a more complicated algorithm and standardization +#' must be handled externally. #' -#' Of note, this function: -#' -#' * does NOT add an intercept -#' * does NOT standardize the design matrix -#' * does NOT set up a path for lambda (the lasso tuning parameter) automatically; -#' This vector of values must be user-supplied. -#' -#' This function is made for use in situations where these steps have already been addressed prior -#' to model fitting. +#' `biglasso_path()` works identically to [biglasso_fit()] except it offers the +#' additional option of fitting models across a path of tuning parameter values. #' -#' In other words, `biglasso_path()` is doing the same thing as `biglasso_fit()`, -#' with the additional option to fit models across a path of tuning parameter values. -#' -#' For now, this function only works with linear regression (`family = 'gaussian'`) +#' Note: +#' * Hybrid safe-strong rules are turned off for `biglasso_fit()`, as these rely +#' on standardization +#' * Currently, the function only works with linear regression +#' (`family = 'gaussian'`). #' #' @param X The design matrix, without an intercept. It must be a -#' double type \code{\link[bigmemory]{big.matrix}} object. +#' double type [bigmemory::big.matrix()] object. #' @param y The response vector #' @param r Residuals (length n vector) corresponding to `init`. #' WARNING: If you supply an incorrect value of `r`, the @@ -28,7 +27,7 @@ #' @param init Initial values for beta. Default: zero (length p vector) #' @param xtx X scales: the jth element should equal `crossprod(X[,j])/n`. #' In particular, if X is standardized, one should pass -#' `xtx = rep(1, p)`. WARNING: If you supply an incorrect value of +#' `xtx = rep(1, p)`. WARNING: If you supply an incorrect value of #' `xtx`, the solution will be incorrect. (length p vector) #' @param penalty String specifying which penalty to use. Default is 'lasso', #' Other options are 'SCAD' and 'MCP' (the latter are non-convex) @@ -37,23 +36,23 @@ #' contribution from the lasso (l1) and the ridge (l2) penalty. #' The penalty is defined as: #' \deqn{ \alpha||\beta||_1 + (1-\alpha)/2||\beta||_2^2.} -#' \code{alpha=1} is the lasso penalty, \code{alpha=0} the ridge penalty, -#' \code{alpha} in between 0 and 1 is the elastic-net ("enet") penalty. +#' `alpha=1` is the lasso penalty, `alpha=0` the ridge penalty, +#' `alpha` in between 0 and 1 is the elastic-net ("enet") penalty. #' @param gamma Tuning parameter value for nonconvex penalty. Defaults are #' 3.7 for `penalty = 'SCAD'` and 3 for `penalty = 'MCP'` #' @param ncores The number of OpenMP threads used for parallel computing. #' @param max.iter Maximum number of iterations. Default is 1000. #' @param eps Convergence threshold for inner coordinate descent. The #' algorithm iterates until the maximum change in the objective -#' after any coefficient update is less than \code{eps} times -#' the null deviance. Default value is \code{1e-7}. +#' after any coefficient update is less than `eps` times +#' the null deviance. Default value is `1e-7`. #' @param dfmax Upper bound for the number of nonzero coefficients. Default is #' no upper bound. However, for large data sets, #' computational burden may be heavy for models with a large #' number of nonzero coefficients. #' @param penalty.factor A multiplicative factor for the penalty applied to -#' each coefficient. If supplied, \code{penalty.factor} must be a numeric -#' vector of length equal to the number of columns of \code{X}. +#' each coefficient. If supplied, `penalty.factor` must be a numeric +#' vector of length equal to the number of columns of `X`. #' @param warn Return warning messages for failures to converge and model #' saturation? Default is TRUE. #' @param output.time Whether to print out the start and end time of the model @@ -61,9 +60,9 @@ #' @param return.time Whether to return the computing time of the model #' fitting. Default is TRUE. #' -#' @return An object with S3 class \code{"biglasso"} with following variables. +#' @returns An object with S3 class `"biglasso"` with following variables. #' \item{beta}{A sparse matrix where rows are estimates a given coefficient across all values of lambda} -#' \item{iter}{A vector of length \code{nlambda} containing the number of +#' \item{iter}{A vector of length `nlambda` containing the number of #' iterations until convergence} #' \item{resid}{Vector of residuals calculated from estimated coefficients.} #' \item{lambda}{The sequence of regularization parameter values in the path.} @@ -72,27 +71,27 @@ #' \item{penalty.factor}{Same as in `biglasso()`.} #' \item{n}{The number of observations used in the model fitting.} #' \item{y}{The response vector used in the model fitting.} +#' #' @author Tabitha Peter and Patrick Breheny #' #' @examples -#' #' data(Prostate) -#' X <- cbind(1, Prostate$X) |> ncvreg::std() # standardizing -> xtx is all 1s +#' X <- cbind(1, Prostate$X) +#' xtx <- apply(X, 2, crossprod)/nrow(X) #' y <- Prostate$y #' X.bm <- as.big.matrix(X) -#' init <- rep(0, ncol(X)) # using cold starts - will need more iterations -#' r <- y - X%*%init -#' fit_lasso <- biglasso_path(X = X.bm, y = y, r = r, init = init, -#' xtx = rep(1, ncol(X)), lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), -#' penalty.factor=c(0, rep(1, ncol(X)-1)), -#' max.iter = 10000) -#' -#' fit_mcp <- biglasso_path(X = X.bm, y = y, r = r, init = init, -#' xtx = rep(1, ncol(X)), lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), -#' penalty.factor=c(0, rep(1, ncol(X)-1)), -#' max.iter = 10000, penalty= 'MCP') +#' init <- rep(0, ncol(X)) +#' fit <- biglasso_path(X = X.bm, y = y, r = y, init = init, xtx = xtx, +#' lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), +#' penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter=2000) +#' fit$beta #' +#' fit <- biglasso_path(X = X.bm, y = y, r = y, init = init, xtx = xtx, +#' lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), penalty='MCP', +#' penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter = 2000) +#' fit$beta #' @export biglasso_path + biglasso_path <- function(X, y, r, @@ -206,4 +205,3 @@ biglasso_path <- function(X, val <- structure(return.val, class = c("biglasso", 'ncvreg')) } - diff --git a/R/cv.biglasso.R b/R/cv-biglasso.R similarity index 66% rename from R/cv.biglasso.R rename to R/cv-biglasso.R index be8a62a..75c6b03 100644 --- a/R/cv.biglasso.R +++ b/R/cv-biglasso.R @@ -3,69 +3,66 @@ #' Perform k-fold cross validation for penalized regression models over a grid #' of values for the regularization parameter lambda. #' -#' The function calls \code{biglasso} \code{nfolds} times, each time leaving -#' out 1/\code{nfolds} of the data. The cross-validation error is based on the -#' residual sum of squares when \code{family="gaussian"} and the binomial -#' deviance when \code{family="binomial"}.\cr \cr The S3 class object -#' \code{cv.biglasso} inherits class \code{\link[ncvreg]{cv.ncvreg}}. So S3 -#' functions such as \code{"summary", "plot"} can be directly applied to the -#' \code{cv.biglasso} object. +#' The function calls `biglasso` `nfolds` times, each time leaving +#' out 1/`nfolds` of the data. The cross-validation error is based on the +#' residual sum of squares when `family="gaussian"` and the binomial +#' deviance when `family="binomial"`. #' -#' @param X The design matrix, without an intercept, as in -#' \code{\link{biglasso}}. -#' @param y The response vector, as in \code{biglasso}. -#' @param row.idx The integer vector of row indices of \code{X} that used for -#' fitting the model. as in \code{biglasso}. -#' @param family Either \code{"gaussian"}, \code{"binomial"}, \code{"cox"} or -#' \code{"mgaussian"} depending on the response. \code{"cox"} and \code{"mgaussian"} -#' are not supported yet. +#' The S3 class object `cv.biglasso` inherits class [ncvreg::cv.ncvreg()]. So S3 +#' functions such as `"summary", "plot"` can be directly applied to the +#' `cv.biglasso` object. +#' +#' @param X The design matrix, without an intercept, as in [biglasso()]. +#' @param y The response vector, as in `biglasso`. +#' @param row.idx The integer vector of row indices of `X` that used for +#' fitting the model. as in `biglasso`. +#' @param family Either `"gaussian"`, `"binomial"`, `"cox"` or +#' `"mgaussian"` depending on the response. `"cox"` and `"mgaussian"` +#' are not supported yet. #' @param eval.metric The evaluation metric for the cross-validated error and -#' for choosing optimal \code{lambda}. "default" for linear regression is MSE -#' (mean squared error), for logistic regression is binomial deviance. -#' "MAPE", for linear regression only, is the Mean Absolute Percentage Error. -#' "auc", for binary classification, is the area under the receiver operating -#' characteristic curve (ROC). -#' "class", for binary classification, gives the misclassification error. +#' for choosing optimal `lambda`. "default" for linear regression is MSE +#' (mean squared error), for logistic regression is binomial deviance. +#' "MAPE", for linear regression only, is the Mean Absolute Percentage Error. +#' "auc", for binary classification, is the area under the receiver operating +#' characteristic curve (ROC). +#' "class", for binary classification, gives the misclassification error. #' @param ncores The number of cores to use for parallel execution of the -#' cross-validation folds, run on a cluster created by the \code{parallel} -#' package. (This is also supplied to the \code{ncores} argument in -#' \code{\link{biglasso}}, which is the number of OpenMP threads, but only for -#' the first call of \code{\link{biglasso}} that is run on the entire data. The -#' individual calls of \code{\link{biglasso}} for the CV folds are run without -#' the \code{ncores} argument.) -#' @param ... Additional arguments to \code{biglasso}. +#' cross-validation folds, run on a cluster created by the `parallel` +#' package. (This is also supplied to the `ncores` argument in +#' [biglasso()], which is the number of OpenMP threads, but only for +#' the first call of [biglasso()] that is run on the entire data. The +#' individual calls of [biglasso()] for the CV folds are run without +#' the `ncores` argument.) +#' @param ... Additional arguments to `biglasso`. #' @param nfolds The number of cross-validation folds. Default is 5. #' @param seed The seed of the random number generator in order to obtain -#' reproducible results. +#' reproducible results. #' @param cv.ind Which fold each observation belongs to. By default the -#' observations are randomly assigned by \code{cv.biglasso}. +#' observations are randomly assigned by `cv.biglasso`. #' @param trace If set to TRUE, cv.biglasso will inform the user of its -#' progress by announcing the beginning of each CV fold. Default is FALSE. -#' @param grouped Whether to calculate CV standard error (\code{cvse}) over -#' CV folds (\code{TRUE}), or over all cross-validated predictions. Ignored -#' when \code{eval.metric} is 'auc'. +#' progress by announcing the beginning of each CV fold. Default is FALSE. +#' @param grouped Whether to calculate CV standard error (`cvse`) over +#' CV folds (`TRUE`), or over all cross-validated predictions. Ignored +#' when `eval.metric` is 'auc'. #' -#' @return An object with S3 class \code{"cv.biglasso"} which inherits from -#' class \code{"cv.ncvreg"}. The following variables are contained in the -#' class (adopted from \code{\link[ncvreg]{cv.ncvreg}}). \item{cve}{The error -#' for each value of \code{lambda}, averaged across the cross-validation -#' folds.} \item{cvse}{The estimated standard error associated with each value -#' of for \code{cve}.} \item{lambda}{The sequence of regularization parameter -#' values along which the cross-validation error was calculated.} -#' \item{fit}{The fitted \code{biglasso} object for the whole data.} -#' \item{min}{The index of \code{lambda} corresponding to \code{lambda.min}.} -#' \item{lambda.min}{The value of \code{lambda} with the minimum -#' cross-validation error.} \item{lambda.1se}{The largest value of \code{lambda} -#' for which the cross-validation error is at most one standard error larger -#' than the minimum cross-validation error.} \item{null.dev}{The deviance for -#' the intercept-only model.} \item{pe}{If \code{family="binomial"}, the -#' cross-validation prediction error for each value of \code{lambda}.} +#' @returns An object with S3 class `"cv.biglasso"` which inherits from +#' class `"cv.ncvreg"`. The following variables are contained in the +#' class (adopted from [ncvreg::cv.ncvreg()]). +#' \item{cve}{The error for each value of `lambda`, averaged across the cross-validation folds.} +#' \item{cvse}{The estimated standard error associated with each value of for `cve`.} +#' \item{lambda}{The sequence of regularization parameter values along which the cross-validation error was calculated.} +#' \item{fit}{The fitted `biglasso` object for the whole data.} +#' \item{min}{The index of `lambda` corresponding to `lambda.min`.} +#' \item{lambda.min}{The value of `lambda` with the minimum cross-validation error.} +#' \item{lambda.1se}{The largest value of `lambda` for which the cross-validation error is at most one standard error larger than the minimum cross-validation error.} +#' \item{null.dev}{The deviance for the intercept-only model.} +#' \item{pe}{If `family="binomial"`, the cross-validation prediction error for each value of `lambda`.} #' \item{cv.ind}{Same as above.} +#' #' @author Yaohui Zeng and Patrick Breheny #' -#' Maintainer: Yaohui Zeng -#' @seealso \code{\link{biglasso}}, \code{\link{plot.cv.biglasso}}, -#' \code{\link{summary.cv.biglasso}}, \code{\link{setupX}} +#' @seealso [biglasso()], [plot.cv.biglasso()], [summary.cv.biglasso()], [setupX()] +#' #' @examples #' \dontrun{ #' ## cv.biglasso @@ -80,9 +77,8 @@ #' plot(cvfit, type = 'all') #' summary(cvfit) #' } -#' -#' @export cv.biglasso -#' +#' @export + cv.biglasso <- function(X, y, row.idx = 1:nrow(X), family = c("gaussian", "binomial", "cox", "mgaussian"), eval.metric = c("default", "MAPE", "auc", "class"), diff --git a/R/loss.R b/R/loss.R index e381565..2dfef78 100644 --- a/R/loss.R +++ b/R/loss.R @@ -2,26 +2,26 @@ #' #' Internal biglasso functions #' -#' These are not intended for use by users.\code{loss.biglasso} calculates the +#' These are not intended for use by users. `loss.biglasso` calculates the #' value of the loss function for the given predictions (used for cross-validation). #' #' @param y The observed response vector. #' @param yhat The predicted response vector. #' @param family Either "gaussian" or "binomial", depending on the response. #' @param eval.metric The evaluation metric for the cross-validated error and -#' for choosing optimal \code{lambda}. "default" for linear regression is MSE -#' (mean squared error), for logistic regression is misclassification error. -#' "MAPE", for linear regression only, is the Mean Absolute Percentage Error. -#' "auc", for logistic regression, is the area under the receiver operating -#' characteristic curve (ROC). +#' for choosing optimal \code{lambda}. "default" for linear regression is MSE +#' (mean squared error), for logistic regression is misclassification error. +#' "MAPE", for linear regression only, is the Mean Absolute Percentage Error. +#' "auc", for logistic regression, is the area under the receiver operating +#' characteristic curve (ROC). #' @param grouped Whether to calculate loss for the entire CV fold -#' (\code{TRUE}), or for predictions individually. Must be \code{TRUE} when -#' \code{eval.metric} is 'auc'. +#' (`TRUE`), or for predictions individually. Must be `TRUE` when +#' `eval.metric` is 'auc'. +#' #' @author Yaohui Zeng and Patrick Breheny #' -#' Maintainer: Yaohui Zeng #' @keywords internal -#' + loss.biglasso <- function(y, yhat, family, eval.metric, grouped = TRUE) { n <- length(y) if (!is.matrix(yhat)) { diff --git a/R/plot.biglasso.R b/R/plot-biglasso.R similarity index 80% rename from R/plot.biglasso.R rename to R/plot-biglasso.R index 4f3bc4b..bd34cfc 100644 --- a/R/plot.biglasso.R +++ b/R/plot-biglasso.R @@ -1,25 +1,21 @@ #' Plot coefficients from a "biglasso" object #' -#' Produce a plot of the coefficient paths for a fitted \code{\link{biglasso}} -#' object. +#' Produce a plot of the coefficient paths for a fitted [biglasso()] object. #' -#' -#' @param x Fitted \code{"biglasso"} model. +#' @param x Fitted [biglasso()] model. #' @param alpha Controls alpha-blending, helpful when the number of covariates -#' is large. Default is alpha=1. +#' is large. Default is alpha=1. #' @param log.l Should horizontal axis be on the log scale? Default is TRUE. -#' @param \dots Other graphical parameters to \code{plot} +#' @param \dots Other graphical parameters to [plot()] +#' #' @author Yaohui Zeng and Patrick Breheny #' -#' Maintainer: Yaohui Zeng -#' @seealso \code{\link{biglasso}}, \code{\link{cv.biglasso}} -#' @keywords models regression -#' @examples +#' @seealso [biglasso()], [cv.biglasso()] #' +#' @examples #' ## See examples in "biglasso" -#' #' @export -#' + plot.biglasso <- function(x, alpha = 1, log.l = TRUE, ...) { YY <- if (length(x$penalty.factor)==nrow(x$beta)) coef(x) else coef(x)[-1,,drop=FALSE] diff --git a/R/plot-cv-biglasso.R b/R/plot-cv-biglasso.R new file mode 100644 index 0000000..de81ef2 --- /dev/null +++ b/R/plot-cv-biglasso.R @@ -0,0 +1,42 @@ +#' Plots the cross-validation curve from a "cv.biglasso" object +#' +#' Plot the cross-validation curve from a [cv.biglasso()] object, +#' along with standard error bars. +#' +#' Error bars representing approximate 68\% confidence intervals are plotted +#' along with the estimates at value of `lambda`. For `rsq` and +#' `snr`, these confidence intervals are quite crude, especially near. +#' +#' @param x A `"cv.biglasso"` object. +#' @param log.l Should horizontal axis be on the log scale? Default is TRUE. +#' @param type What to plot on the vertical axis. `cve` plots the +#' cross-validation error (deviance); `rsq` plots an estimate of the +#' fraction of the deviance explained by the model (R-squared); `snr` +#' plots an estimate of the signal-to-noise ratio; `scale` plots, for +#' `family="gaussian"`, an estimate of the scale parameter (standard +#' deviation); `pred` plots, for `family="binomial"`, the estimated +#' prediction error; `all` produces all of the above. +#' @param selected If `TRUE` (the default), places an axis on top of the +#' plot denoting the number of variables in the model (i.e., that have a +#' nonzero regression coefficient) at that value of `lambda`. +#' @param vertical.line If `TRUE` (the default), draws a vertical line at +#' the value where cross-validaton error is minimized. +#' @param col Controls the color of the dots (CV estimates). +#' @param \dots Other graphical parameters to `plot` +#' +#' @author Yaohui Zeng and Patrick Breheny +#' +#' @seealso [biglasso()], [cv.biglasso()] +#' +#' @examples +#' ## See examples in "cv.biglasso" +#' @export + +plot.cv.biglasso <- function(x, log.l = TRUE, type = c("cve", "rsq", "scale", + "snr", "pred", "all"), + selected = TRUE, vertical.line = TRUE, col = "red", ...) { + # inherits cv.ncvreg + class(x) <- 'cv.ncvreg' + plot(x = x, log.l = log.l, type = type, selected = selected, + vertical.line = vertical.line, col = col, ...) +} diff --git a/R/plot.mbiglasso.R b/R/plot-mbiglasso.R similarity index 90% rename from R/plot.mbiglasso.R rename to R/plot-mbiglasso.R index 9baffaf..b6c03a4 100644 --- a/R/plot.mbiglasso.R +++ b/R/plot-mbiglasso.R @@ -1,27 +1,24 @@ #' Plot coefficients from a "mbiglasso" object #' #' Produce a plot of the coefficient paths for a fitted multiple responses -#' \code{mbiglasso} object. +#' `mbiglasso` object. #' -#' -#' @param x Fitted \code{"mbiglasso"} model. +#' @param x Fitted `mbiglasso` model. #' @param alpha Controls alpha-blending, helpful when the number of covariates #' is large. Default is alpha=1. #' @param log.l Should horizontal axis be on the log scale? Default is TRUE. #' @param norm.beta Should the vertical axis be the l2 norm of coefficients for each variable? #' Default is TRUE. If False, the vertical axis is the coefficients. -#' @param \dots Other graphical parameters to \code{plot} +#' @param \dots Other graphical parameters to [plot()] +#' #' @author Chuyi Wang #' -#' Maintainer: Chuyi Wang -#' @seealso \code{\link{biglasso}}, -#' @keywords models regression -#' @examples +#' @seealso [biglasso()] #' +#' @examples #' ## See examples in "biglasso" -#' #' @export -#' + plot.mbiglasso <- function(x, alpha = 1, log.l = TRUE, norm.beta = TRUE, ...) { YY <- coef(x, intercept = FALSE) ## currently not support unpenalized coefficients. NOT USED diff --git a/R/plot.cv.biglasso.R b/R/plot.cv.biglasso.R deleted file mode 100644 index c1a88f8..0000000 --- a/R/plot.cv.biglasso.R +++ /dev/null @@ -1,44 +0,0 @@ -#' Plots the cross-validation curve from a "cv.biglasso" object -#' -#' Plot the cross-validation curve from a \code{\link{cv.biglasso}} object, -#' along with standard error bars. -#' -#' Error bars representing approximate 68\% confidence intervals are plotted -#' along with the estimates at value of \code{lambda}. For \code{rsq} and -#' \code{snr}, these confidence intervals are quite crude, especially near. -#' -#' @param x A \code{"cv.biglasso"} object. -#' @param log.l Should horizontal axis be on the log scale? Default is TRUE. -#' @param type What to plot on the vertical axis. \code{cve} plots the -#' cross-validation error (deviance); \code{rsq} plots an estimate of the -#' fraction of the deviance explained by the model (R-squared); \code{snr} -#' plots an estimate of the signal-to-noise ratio; \code{scale} plots, for -#' \code{family="gaussian"}, an estimate of the scale parameter (standard -#' deviation); \code{pred} plots, for \code{family="binomial"}, the estimated -#' prediction error; \code{all} produces all of the above. -#' @param selected If \code{TRUE} (the default), places an axis on top of the -#' plot denoting the number of variables in the model (i.e., that have a -#' nonzero regression coefficient) at that value of \code{lambda}. -#' @param vertical.line If \code{TRUE} (the default), draws a vertical line at -#' the value where cross-validaton error is minimized. -#' @param col Controls the color of the dots (CV estimates). -#' @param \dots Other graphical parameters to \code{plot} -#' @author Yaohui Zeng and Patrick Breheny -#' -#' Maintainer: Yaohui Zeng -#' @seealso \code{\link{biglasso}}, \code{\link{cv.biglasso}} -#' @keywords models regression -#' @examples -#' -#' ## See examples in "cv.biglasso" -#' -#' @export -#' -plot.cv.biglasso <- function(x, log.l = TRUE, type = c("cve", "rsq", "scale", - "snr", "pred", "all"), - selected = TRUE, vertical.line = TRUE, col = "red", ...) { - # inherits cv.ncvreg - class(x) <- 'cv.ncvreg' - plot(x = x, log.l = log.l, type = type, selected = selected, - vertical.line = vertical.line, col = col, ...) -} diff --git a/R/predict.cv.R b/R/predict-cv.R similarity index 54% rename from R/predict.cv.R rename to R/predict-cv.R index e8f974f..dabb271 100644 --- a/R/predict.cv.R +++ b/R/predict-cv.R @@ -1,42 +1,43 @@ -#' Model predictions based on a fitted \code{\link{cv.biglasso}} object +#' Model predictions based on a fitted [cv.biglasso()] object #' -#' Extract predictions from a fitted \code{\link{cv.biglasso}} object. +#' Extract predictions from a fitted [cv.biglasso()] object. #' #' @name predict.cv.biglasso #' @rdname predict.cv.biglasso #' @method predict cv.biglasso #' -#' @param object A fitted \code{"cv.biglasso"} model object. +#' @param object A fitted `"cv.biglasso"` model object. #' @param X Matrix of values at which predictions are to be made. It must be a -#' \code{\link[bigmemory]{big.matrix}} object. Not used for -#' \code{type="coefficients"}. -#' @param row.idx Similar to that in \code{\link[biglasso]{biglasso}}, it's a -#' vector of the row indices of \code{X} that used for the prediction. -#' \code{1:nrow(X)} by default. -#' @param type Type of prediction: \code{"link"} returns the linear predictors; -#' \code{"response"} gives the fitted values; \code{"class"} returns the -#' binomial outcome with the highest probability; \code{"coefficients"} returns -#' the coefficients; \code{"vars"} returns a list containing the indices and -#' names of the nonzero variables at each value of \code{lambda}; -#' \code{"nvars"} returns the number of nonzero coefficients at each value of -#' \code{lambda}. -#' @param lambda Values of the regularization parameter \code{lambda} at which +#' [bigmemory::big.matrix()] object. Not used for +#' `type="coefficients"`. +#' @param row.idx Similar to that in [biglasso()], it's a +#' vector of the row indices of `X` that used for the prediction. +#' `1:nrow(X)` by default. +#' @param type Type of prediction: +#' * `"link"` returns the linear predictors +#' * `"response"` gives the fitted values +#' * `"class"` returns the binomial outcome with the highest probability +#' * `"coefficients"` returns the coefficients +#' * `"vars"` returns a list containing the indices and names of the nonzero variables at each value of `lambda` +#' * `"nvars"` returns the number of nonzero coefficients at each value of `lambda` +#' @param lambda Values of the regularization parameter `lambda` at which #' predictions are requested. The default value is the one corresponding to #' the minimum cross-validation error. Accepted values are also the strings -#' "lambda.min" (\code{lambda} of minimum cross-validation error) and -#' "lambda.1se" (Largest value of \code{lambda} for which the cross-validation +#' "lambda.min" (`lambda` of minimum cross-validation error) and +#' "lambda.1se" (Largest value of `lambda` for which the cross-validation #' error was at most one standard error larger than the minimum.). -#' @param which Indices of the penalty parameter \code{lambda} at which +#' @param which Indices of the penalty parameter `lambda` at which #' predictions are requested. The default value is the index of lambda -#' corresponding to lambda.min. Note: this is overridden if \code{lambda} is +#' corresponding to lambda.min. Note: this is overridden if `lambda` is #' specified. #' @param \dots Not used. -#' @return The object returned depends on \code{type}. +#' +#' @returns The object returned depends on `type`. +#' #' @author Yaohui Zeng and Patrick Breheny #' -#' Maintainer: Yaohui Zeng -#' @seealso \code{\link{biglasso}}, \code{\link{cv.biglasso}} -#' @keywords models regression +#' @seealso [biglasso()], [cv.biglasso()] +#' #' @examples #' \dontrun{ #' ## predict.cv.biglasso @@ -54,7 +55,7 @@ #' predict(cvfit, X.bm, lambda = "lambda.1se") #' } #' @export -#' + predict.cv.biglasso <- function(object, X, row.idx = 1:nrow(X), type = c("link","response","class", "coefficients","vars","nvars"), diff --git a/R/predict.R b/R/predict.R index 18795e9..3045cca 100644 --- a/R/predict.R +++ b/R/predict.R @@ -1,46 +1,46 @@ -#' Model predictions based on a fitted \code{biglasso} object +#' Model predictions based on a fitted `biglasso` object #' #' Extract predictions (fitted reponse, coefficients, etc.) from a -#' fitted \code{\link{biglasso}} object. +#' fitted [biglasso()] object. #' #' @name predict.biglasso #' @rdname predict.biglasso #' @method predict biglasso #' -#' @param object A fitted \code{"biglasso"} model object. +#' @param object A fitted `"biglasso"` model object. #' @param X Matrix of values at which predictions are to be made. It must be a -#' \code{\link[bigmemory]{big.matrix}} object. Not used for -#' \code{type="coefficients"}. -#' @param row.idx Similar to that in \code{\link[biglasso]{biglasso}}, it's a -#' vector of the row indices of \code{X} that used for the prediction. -#' \code{1:nrow(X)} by default. -#' @param type Type of prediction: \code{"link"} returns the linear predictors; -#' \code{"response"} gives the fitted values; \code{"class"} returns the -#' binomial outcome with the highest probability; \code{"coefficients"} returns -#' the coefficients; \code{"vars"} returns a list containing the indices and -#' names of the nonzero variables at each value of \code{lambda}; -#' \code{"nvars"} returns the number of nonzero coefficients at each value of -#' \code{lambda}. -#' @param lambda Values of the regularization parameter \code{lambda} at which +#' [bigmemory::big.matrix()] object. Not used for `type="coefficients"`. +#' @param row.idx Similar to that in [biglasso()], it's a +#' vector of the row indices of `X` that used for the prediction. +#' `1:nrow(X)` by default. +#' @param type Type of prediction: +#' * `"link"` returns the linear predictors +#' * `"response"` gives the fitted values +#' * `"class"` returns the binomial outcome with the highest probability +#' * `"coefficients"` returns the coefficients +#' * `"vars"` returns a list containing the indices and names of the nonzero variables at each value of `lambda` +#' * `"nvars"` returns the number of nonzero coefficients at each value of `lambda` +#' @param lambda Values of the regularization parameter `lambda` at which #' predictions are requested. Linear interpolation is used for values of -#' \code{lambda} not in the sequence of lambda values in the fitted models. +#' `lambda` not in the sequence of lambda values in the fitted models. #' @param k Index of the response to predict in multiple responses regression ( -#' \code{family="mgaussian"}). -#' @param which Indices of the penalty parameter \code{lambda} at which +#' `family="mgaussian"`). +#' @param which Indices of the penalty parameter `lambda` at which #' predictions are required. By default, all indices are returned. If -#' \code{lambda} is specified, this will override \code{which}. +#' `lambda` is specified, this will override `which`. #' @param intercept Whether the intercept should be included in the returned -#' coefficients. For \code{family="mgaussian"} only. -#' @param drop If coefficients for a single value of \code{lambda} are to be -#' returned, reduce dimensions to a vector? Setting \code{drop=FALSE} returns +#' coefficients. For `family="mgaussian"` only. +#' @param drop If coefficients for a single value of `lambda` are to be +#' returned, reduce dimensions to a vector? Setting `drop=FALSE` returns #' a 1-column matrix. #' @param \dots Not used. -#' @return The object returned depends on \code{type}. +#' +#' @returns The object returned depends on `type`. +#' #' @author Yaohui Zeng and Patrick Breheny #' -#' Maintainer: Yaohui Zeng -#' @seealso \code{\link{biglasso}}, \code{\link{cv.biglasso}} -#' @keywords models regression +#' @seealso [biglasso()], [cv.biglasso()] +#' #' @examples #' ## Logistic regression #' data(colon) @@ -173,4 +173,4 @@ coef.mbiglasso <- function(object, lambda, which = 1:length(object$lambda), inte } } return(beta) -} \ No newline at end of file +} diff --git a/R/setupX.R b/R/setupX.R index 242a3a5..6d9e3c4 100644 --- a/R/setupX.R +++ b/R/setupX.R @@ -1,56 +1,52 @@ -## read and set up design matrix X from external ASCII-file - - #' Set up design matrix X by reading data from big data file #' -#' Set up the design matrix X as a \code{big.matrix} object based on external +#' Set up the design matrix X as a `big.matrix` object based on external #' massive data file stored on disk that cannot be fullly loaded into memory. #' The data file must be a well-formated ASCII-file, and contains only one -#' single type. Current version only supports \code{double} type. Other +#' single type. Current version only supports `double` type. Other #' restrictions about the data file are described in -#' \code{\link{biglasso-package}}. This function reads the massive data, and -#' creates a \code{big.matrix} object. By default, the resulting -#' \code{big.matrix} is file-backed, and can be shared across processors or +#' [biglasso-package]. This function reads the massive data, and +#' creates a `big.matrix` object. By default, the resulting +#' `big.matrix` is file-backed, and can be shared across processors or #' nodes of a cluster. #' #' For a data set, this function needs to be called only one time to set up the -#' \code{big.matrix} object with two backing files (.bin, .desc) created in +#' `big.matrix` object with two backing files (.bin, .desc) created in #' current working directory. Once set up, the data can be "loaded" into any -#' (new) R session by calling \code{attach.big.matrix(discriptorfile)}. +#' (new) R session by calling `attach.big.matrix(discriptorfile)`. #' -#' This function is a simple wrapper of -#' \code{\link[bigmemory]{read.big.matrix}}. See -#' \code{\link[bigmemory]{read.big.matrix}} and the package +#' This function is a simple wrapper of [bigmemory::read.big.matrix()]. See #' \href{https://CRAN.R-project.org/package=bigmemory}{bigmemory} for more #' details. #' #' @param filename The name of the data file. For example, "dat.txt". #' @param dir The directory used to store the binary and descriptor files -#' associated with the \code{big.matrix}. The default is current working +#' associated with the `big.matrix`. The default is current working #' directory. #' @param sep The field separator character. For example, "," for #' comma-delimited files (the default); "\\t" for tab-delimited files. #' @param backingfile The binary file associated with the file-backed -#' \code{big.matrix}. By default, its name is the same as \code{filename} with +#' `big.matrix`. By default, its name is the same as `filename` with #' the extension replaced by ".bin". #' @param descriptorfile The descriptor file used for the description of the -#' file-backed \code{big.matrix}. By default, its name is the same as -#' \code{filename} with the extension replaced by ".desc". +#' file-backed `big.matrix`. By default, its name is the same as +#' `filename` with the extension replaced by ".desc". #' @param type The data type. Only "double" is supported for now. #' @param ... Additional arguments that can be passed into function -#' \code{\link[bigmemory]{read.big.matrix}}. -#' @return A \code{big.matrix} object corresponding to a file-backed -#' \code{big.matrix}. It's ready to be used as the design matrix \code{X} in -#' \code{\link{biglasso}} and \code{\link{cv.biglasso}}. +#' [bigmemory::read.big.matrix()]. +#' +#' @returns A `big.matrix` object corresponding to a file-backed +#' [bigmemory::big.matrix()]. It's ready to be used as the design matrix `X` in +#' [biglasso()] and [cv.biglasso()]. +#' #' @author Yaohui Zeng and Patrick Breheny #' -#' Maintainer: Yaohui Zeng -#' @seealso \code{\link{biglasso}}, \code{\link{cv.ncvreg}} +#' @seealso [biglasso()], [cv.ncvreg()], [biglasso-package] +#' #' @examples #' ## see the example in "biglasso-package" -#' -#' @export setupX -#' +#' @export + setupX <- function(filename, dir = getwd(), sep = ",", backingfile = paste0(unlist(strsplit(filename, split = "\\."))[1], diff --git a/R/summary.cv.biglasso.R b/R/summary-cv-biglasso.R similarity index 64% rename from R/summary.cv.biglasso.R rename to R/summary-cv-biglasso.R index 63f6522..54aac46 100644 --- a/R/summary.cv.biglasso.R +++ b/R/summary-cv-biglasso.R @@ -1,42 +1,41 @@ #' Summarizing inferences based on cross-validation #' -#' Summary method for \code{cv.biglasso} objects. +#' Summary method for `cv.biglasso` objects. #' #' @name summary.cv.biglasso #' @rdname summary.cv.biglasso #' @method summary cv.biglasso #' -#' @param object A \code{cv.biglasso} object. -#' @param x A \code{"summary.cv.biglasso"} object. +#' @param object A `cv.biglasso` object. +#' @param x A `"summary.cv.biglasso"` object. #' @param digits Number of digits past the decimal point to print out. Can be #' a vector specifying different display digits for each of the five #' non-integer printed values. #' @param \dots Further arguments passed to or from other methods. -#' @return \code{summary.cv.biglasso} produces an object with S3 class -#' \code{"summary.cv.biglasso"}. The class has its own print method and contains -#' the following list elements: \item{penalty}{The penalty used by -#' \code{biglasso}.} \item{model}{Either \code{"linear"} or \code{"logistic"}, -#' depending on the \code{family} option in \code{biglasso}.} \item{n}{Number -#' of observations} \item{p}{Number of regression coefficients (not including -#' the intercept).} \item{min}{The index of \code{lambda} with the smallest -#' cross-validation error.} \item{lambda}{The sequence of \code{lambda} values -#' used by \code{cv.biglasso}.} \item{cve}{Cross-validation error (deviance).} -#' \item{r.squared}{Proportion of variance explained by the model, as estimated -#' by cross-validation.} \item{snr}{Signal to noise ratio, as estimated by -#' cross-validation.} \item{sigma}{For linear regression models, the scale -#' parameter estimate.} \item{pe}{For logistic regression models, the -#' prediction error (misclassification error).} +#' +#' @returns `summary.cv.biglasso` produces an object with S3 class +#' `"summary.cv.biglasso"`. The class has its own print method and contains +#' the following list elements: +#' \item{penalty}{The penalty used by `biglasso`.} +#' \item{model}{Either `"linear"` or `"logistic"`, depending on the `family` option in `biglasso`.} +#' \item{n}{Number of observations} +#' \item{p}{Number of regression coefficients (not including the intercept).} +#' \item{min}{The index of `lambda` with the smallest cross-validation error.} +#' \item{lambda}{The sequence of `lambda` values used by `cv.biglasso`.} +#' \item{cve}{Cross-validation error (deviance).} +#' \item{r.squared}{Proportion of variance explained by the model, as estimated by cross-validation.} +#' \item{snr}{Signal to noise ratio, as estimated by cross-validation.} +#' \item{sigma}{For linear regression models, the scale parameter estimate.} +#' \item{pe}{For logistic regression models, the prediction error (misclassification error).} +#' #' @author Yaohui Zeng and Patrick Breheny #' -#' Maintainer: Yaohui Zeng -#' @seealso \code{\link{biglasso}}, \code{\link{cv.biglasso}}, -#' \code{\link{plot.cv.biglasso}} -#' @keywords models regression +#' @seealso [biglasso()], [cv.biglasso()], [plot.cv.biglasso()], [biglasso-package] +#' #' @examples #' ## See examples in "cv.biglasso" and "biglasso-package" #' @export -#' -#' + summary.cv.biglasso <- function(object, ...) { S <- pmax(object$null.dev - object$cve, 0) if (!inherits(object, 'cv.ncvsurv') && object$fit$family=="gaussian") { diff --git a/README.md b/README.md index dba167e..e302796 100644 --- a/README.md +++ b/README.md @@ -1,15 +1,15 @@ -[![GitHub version](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/YaohuiZeng/biglasso/master/.version.json&style=flat&logo=github)](https://github.com/YaohuiZeng/biglasso) +[![GitHub version](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/pbreheny/biglasso/master/.version.json&style=flat&logo=github)](https://github.com/pbreheny/biglasso) [![CRAN version](https://img.shields.io/cran/v/biglasso?logo=R)](https://cran.r-project.org/package=biglasso) [![downloads](https://cranlogs.r-pkg.org/badges/biglasso)](https://cran.r-project.org/package=biglasso) -[![R-CMD-check](https://github.com/YaohuiZeng/biglasso/workflows/R-CMD-check/badge.svg)](https://github.com/YaohuiZeng/biglasso/actions) +[![R-CMD-check](https://github.com/pbreheny/biglasso/workflows/R-CMD-check/badge.svg)](https://github.com/pbreheny/biglasso/actions) -# [biglasso: Extend Lasso Model Fitting to Big Data in R](https://yaohuizeng.github.io/biglasso/index.html) +# [biglasso: Extend Lasso Model Fitting to Big Data in R](https://pbreheny.github.io/biglasso/index.html) `biglasso` extends lasso and elastic-net linear and logistic regression models for ultrahigh-dimensional, multi-gigabyte data sets that cannot be loaded into memory. It utilizes memory-mapped files to store the massive data on the disk and only read those into memory whenever necessary during model fitting. Moreover, some advanced feature screening rules are proposed and implemented to accelerate the model fitting. **As a result, this package is much more memory- and computation-efficient and highly scalable as compared to existing lasso-fitting packages such as [glmnet](https://CRAN.R-project.org/package=glmnet) and [ncvreg](https://CRAN.R-project.org/package=ncvreg)**. Bechmarking experiments using both simulated and real data sets show that `biglasso` is not only 1.5x to 4x times faster than existing packages, but also at least 2x more memory-efficient. More importantly, to the best of our knowledge, `biglasso` is the first R package that enables users to fit lasso models with data sets that are larger than available RAM, thus allowing for powerful big data analysis on an ordinary laptop. -## Installation: +## Installation To install the latest stable release version from CRAN: @@ -20,29 +20,30 @@ install.packages("biglasso") To install the latest development version from GitHub: ```r -remotes::install_github("YaohuiZeng/biglasso") +remotes::install_github("pbreheny/biglasso") ``` -## News: -* See NEWS.md for latest news. -* This package was ranked top 3 for [2017 ASA Chambers Statistical Software Award](http://stat-computing.org/awards/jmc/). -* The technical paper of this package was selected as a Winner of [2017 ASA Student Paper Competiton from Section on Statistical Computing](http://stat-computing.org/awards/student/winners.html). +## News +* See [NEWS.md](https://pbreheny.github.io/biglasso/news/index.html) for latest news. +* The technical paper of this package was selected as a Winner of [2017 ASA Student Paper Competiton from Section on Statistical Computing](https://community.amstat.org/jointscsg-section/awards/student-paper-competition). +* This package finished in the top 3 for [2017 ASA Chambers Statistical Software Award](https://community.amstat.org/jointscsg-section/awards/john-m-chambers). -## Documentation: -* Here are the [R Reference manual](https://CRAN.R-project.org/package=biglasso/biglasso.pdf) and [Package Website](https://yaohuizeng.github.io/biglasso/index.html) +## Documentation + +* Here are the [R Reference manual](https://CRAN.R-project.org/package=biglasso/biglasso.pdf) and [Package Website](https://pbreheny.github.io/biglasso/index.html) * Here are the technical papers of the package: i) [The software paper](https://arxiv.org/abs/1701.05936); and ii) [the paper of hybrid safe-strong rules](https://arxiv.org/abs/1704.08742) -## Features: +## Features + 1. It utilizes memory-mapped files to store the massive data on the disk, only loading data into memory when necessary during model fitting. Consequently, it's able to seamlessly handle out-of-core computation. 2. It is built upon pathwise coordinate descent algorithm with *warm start, active set cycling, and feature screening* strategies, which has been proven to be one of fastest lasso solvers. 3. We develop new, adaptive feature screening rules that outperform state-of-the-art screening rules such as the sequential strong rule (SSR) and the sequential EDPP rule (SEDPP) with additional 1.5x to 4x speedup. 4. The implementation is designed to be as memory-efficient as possible by eliminating extra copies of the data created by other R packages, making `biglasso` at least 2x more memory-efficient than `glmnet`. 5. The underlying computation is implemented in C++, and parallel computing with OpenMP is also supported. - ## Benchmarks: ### Simulated data: @@ -63,7 +64,7 @@ remotes::install_github("YaohuiZeng/biglasso") ![Alt text](/vignettes/2020-12-18_vary_n_pkgs.png) --> - + In all the settings, `biglasso` (1 core) is uniformly faster than `picasso`, `glmnet` and `ncvreg`. When the data gets bigger, `biglasso` achieves 6-9x speed-up compared to other packages. @@ -93,8 +94,8 @@ The maximum RSS (in **GB**) used by a single fit and 10-fold cross validation is ### Real data: The performance of the packages are also tested using diverse real data sets: -* [Breast cancer gene expression data](http://myweb.uiowa.edu/pbreheny/data/bcTCGA.html) (GENE); -* [MNIST handwritten image data](http://yann.lecun.com/exdb/mnist/) (MNIST); +* [Breast cancer gene expression data](https://iowabiostat.github.io/data-sets/brca1/brca1.html) (GENE); +* [MNIST handwritten image data](https://www.kaggle.com/datasets/hojjatk/mnist-dataset/data) (MNIST); * [Cardiac fibrosis genome-wide association study data](https://arxiv.org/abs/1607.05636) (GWAS); * [Subset of New York Times bag-of-words data](https://archive.ics.uci.edu/ml/datasets/Bag+of+Words) (NYT). @@ -132,13 +133,15 @@ Since other three packages cannot handle this data-larger-than-RAM case, we comp ## Reference: -* Zeng Y and Breheny P (2021). The biglasso Package: A Memory- and Computation-Efficient Solver for Lasso Model Fitting with Big Data in R. R Journal, 12: 6-19. URL [https://doi.org/10.32614/RJ-2021-001](https://doi.org/10.32614/RJ-2021-001). -* Zeng Y, Yang T, and Breheny P (2021). Hybrid safe-strong rules for efficient optimization in lasso-type problems. Computational Statistics and Data Analysis, 153: 107063. URL [http://www.sciencedirect.com/science/article/pii/S0167947320301547](http://www.sciencedirect.com/science/article/pii/S0167947320301547) -* Wang C and Breheny P (2022). Adaptive hybrid screening for efficient lasso optimization. Journal of Statistical Computation and Simulation, 92: 2233–2256. URL [https://doi.org/10.1080/00949655.2021.2025376](https://doi.org/10.1080/00949655.2021.2025376) + +* Zeng Y and Breheny P (2021). The biglasso Package: A Memory- and Computation-Efficient Solver for Lasso Model Fitting with Big Data in R. R Journal, 12: 6-19. URL +* Zeng Y, Yang T, and Breheny P (2021). Hybrid safe-strong rules for efficient optimization in lasso-type problems. Computational Statistics and Data Analysis, 153: 107063. URL +* Wang C and Breheny P (2022). Adaptive hybrid screening for efficient lasso optimization. Journal of Statistical Computation and Simulation, 92: 2233–2256. URL * Tibshirani, R., Bien, J., Friedman, J., Hastie, T., Simon, N., Taylor, J., and Tibshirani, R. J. (2012). Strong rules for discarding predictors in lasso-type problems. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74 (2), 245-266. * Wang, J., Zhou, J., Wonka, P., and Ye, J. (2013). Lasso screening rules via dual polytope projection. In Advances in Neural Information Processing Systems, pp. 1070-1078. * Xiang, Z. J., and Ramadge, P. J. (2012, March). Fast lasso screening tests based on correlations. In Acoustics, Speech and Signal Processing (ICASSP), 2012 IEEE International Conference on (pp. 2137-2140). IEEE. * Wang, J., Zhou, J., Liu, J., Wonka, P., and Ye, J. (2014). A safe screening rule for sparse logistic regression. In Advances in Neural Information Processing Systems, pp. 1053-1061. ## Report bugs： -* open an [issue](https://github.com/YaohuiZeng/biglasso/issues) or send an email to Patrick Breheny at + +* open an [issue](https://github.com/pbreheny/biglasso/issues) or send an email to Patrick Breheny at diff --git a/inst/CITATION b/inst/CITATION index 26a657d..a4a715d 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -1,16 +1,13 @@ -citHeader("To cite biglasso in publications use:") - -citEntry( - entry ="Article", - author = personList(as.person("Yaohui Zeng"), as.person("Patrick Breheny")), +bibentry( + bibtype="Article", + author = c(person(given="Yaohui", family="Zeng"), + person(given="Patrick", family="Breheny")), title = "The biglasso Package: A Memory- and Computation-Efficient Solver for Lasso Model Fitting with Big Data in R", journal = "R Journal", eprint = "1701.05936", year = "2021", number = "2", pages = "6--19", - volume = "12", - url = "https://doi.org/10.32614/RJ-2021-001", - - textVersion = "Zeng Y and Breheny P (2017). The biglasso Package: A Memory- and Computation-Efficient Solver for Lasso Model Fitting with Big Data in R. R Journal 12: 6-19." -) + volume = "12", + doi = "10.32614/RJ-2021-001", + url = "https://doi.org/10.32614/RJ-2021-001") diff --git a/man/biglasso-package.Rd b/man/biglasso-package.Rd index f8356b7..1022823 100644 --- a/man/biglasso-package.Rd +++ b/man/biglasso-package.Rd @@ -46,20 +46,20 @@ by other R packages, making it at least 2x more memory-efficient than \code{glmnet}. \item the underlying computation is implemented in C++, and parallel computing with OpenMP is also supported. } -\strong{For more information:} \itemize{ \item Benchmarking results: -\url{https://github.com/YaohuiZeng/biglasso}. -\item Tutorial: -\url{http://yaohuizeng.github.io/biglasso/articles/biglasso.html} -\item Technical paper: -\url{https://arxiv.org/abs/1701.05936} } +\strong{For more information:} +\itemize{ +\item Benchmarking results: \url{https://github.com/pbreheny/biglasso} +\item Tutorial: \url{https://pbreheny.github.io/biglasso/articles/biglasso.html} +\item Technical paper: \url{https://arxiv.org/abs/1701.05936} +} } \note{ -The input design matrix X must be a \code{\link[bigmemory]{big.matrix}} object. -This can be created by the function \code{as.big.matrix} in the R package -\href{https://CRAN.R-project.org//package=bigmemory}{bigmemory}. -If the data (design matrix) is very large (e.g. 10 GB) and stored in an external +The input design matrix X must be a \code{\link[bigmemory:big.matrix]{bigmemory::big.matrix()}} object. +This can be created by the function \code{as.big.matrix} in the R package +\href{https://CRAN.R-project.org//package=bigmemory}{bigmemory}. +If the data (design matrix) is very large (e.g. 10 GB) and stored in an external file, which is often the case for big data, X can be created by calling the -function \code{\link{setupX}}. +function \code{\link[=setupX]{setupX()}}. \strong{In this case, there are several restrictions about the data file:} \enumerate{ \item the data file must be a well-formated ASCII-file, with each row corresponding to an observation and each column a variable; \item @@ -71,23 +71,23 @@ Future versions will try to address these restrictions. Denote the number of observations and variables be, respectively, \code{n} and \code{p}. It's worth noting that the package is more suitable for wide -data (ultrahigh-dimensional, \code{p >> n}) as compared to long data -(\code{n >> p}). This is because the model fitting algorithm takes advantage +data (ultrahigh-dimensional, \verb{p >> n}) as compared to long data +(\verb{n >> p}). This is because the model fitting algorithm takes advantage of sparsity assumption of high-dimensional data. To just give the user some ideas, below are some benchmarking results of the total computing time (in seconds) for solving lasso-penalized linear regression along a sequence of 100 values of the tuning parameter. In all cases, assume 20 non-zero coefficients equal +/- 2 in the true model. (Based on Version 1.2-3, screening rule "SSR-BEDPP" is used) -\itemize{ \item For wide data case (\code{p > n}), \code{n = 1,000}: +\itemize{ \item For wide data case (\code{p > n}), \verb{n = 1,000}: \tabular{ccccc}{ \code{p} \tab 1,000 \tab 10,000 \tab 100,000 \tab 1,000,000 \cr Size of \code{X} \tab 9.5 MB \tab 95 MB \tab 950 MB \tab 9.5 GB \cr Elapsed time (s) \tab 0.11 \tab 0.83 \tab 8.47 \tab 85.50 \cr } -%\item For long data case (\code{n >> p}), \code{p = 1,000}: -%\tabular{ccccc}{ -%\code{n} \tab 1,000 \tab 10,000 \tab 100,000 \tab 1,000,000 \cr -%Size of \code{X} \tab 9.5 MB \tab 95 MB \tab 950 MB \tab 9.5 GB \cr -%Elapsed time (s) \tab 2.50 \tab 11.43 \tab 83.69 \tab 1090.62 \cr %} +\%\item For long data case (\verb{n >> p}), \verb{p = 1,000}: +\%\tabular{ccccc}{ +\%\code{n} \tab 1,000 \tab 10,000 \tab 100,000 \tab 1,000,000 \cr +\%Size of \code{X} \tab 9.5 MB \tab 95 MB \tab 950 MB \tab 9.5 GB \cr +\%Elapsed time (s) \tab 2.50 \tab 11.43 \tab 83.69 \tab 1090.62 \cr \%} } } \examples{ @@ -147,8 +147,5 @@ and Signal Processing (ICASSP), 2012 IEEE International Conference on} (pp. Advances in Neural Information Processing Systems}, pp. 1053-1061. } } \author{ -Yaohui Zeng, Chuyi Wang and Patrick Breheny - -Maintainer: Yaohui Zeng and Chuyi Wang +Yaohui Zeng, Chuyi Wang, Tabitha Peter, and Patrick Breheny } -\keyword{package} diff --git a/man/biglasso.Rd b/man/biglasso.Rd index 0143190..a44163b 100644 --- a/man/biglasso.Rd +++ b/man/biglasso.Rd @@ -32,16 +32,16 @@ biglasso( } \arguments{ \item{X}{The design matrix, without an intercept. It must be a -double type \code{\link[bigmemory]{big.matrix}} object. The function +double type \code{\link[bigmemory:big.matrix]{bigmemory::big.matrix()}} object. The function standardizes the data and includes an intercept internally by default during the model fitting.} \item{y}{The response vector for \code{family="gaussian"} or \code{family="binomial"}. For \code{family="cox"}, \code{y} should be a two-column matrix with columns 'time' and 'status'. The latter is a binary variable, with '1' indicating death, - and '0' indicating right censored. For \code{family="mgaussin"}, \code{y} - should be a n*m matrix where n is the sample size and m is the number of - responses.} +and '0' indicating right censored. For \code{family="mgaussin"}, \code{y} +should be a n*m matrix where n is the sample size and m is the number of +responses.} \item{row.idx}{The integer vector of row indices of \code{X} that used for fitting the model. \code{1:nrow(X)} by default.} @@ -71,9 +71,9 @@ applicable to elastic-net-penalized logistic regression or cox regression; (3) active set cycling strategy is incorporated with these screening rules.} \item{safe.thresh}{the threshold value between 0 and 1 that controls when to -stop safe test. For example, 0.01 means to stop safe test at next lambda +stop safe test. For example, 0.01 means to stop safe test at next lambda iteration if the number of features rejected by safe test at current lambda -iteration is not larger than 1\% of the total number of features. So 1 means +iteration is not larger than 1\\% of the total number of features. So 1 means to always turn off safe test, whereas 0 (default) means to turn off safe test if the number of features rejected by safe test is 0 at current lambda.} @@ -136,33 +136,33 @@ Default is FALSE.} } \value{ An object with S3 class \code{"biglasso"} for -\code{"gaussian", "binomial", "cox"} families, or an object with S3 class +\verb{"gaussian", "binomial", "cox"} families, or an object with S3 class \code{"mbiglasso"} for \code{"mgaussian"} family, with following variables. \item{beta}{The fitted matrix of coefficients, store in sparse matrix representation. The number of rows is equal to the number of coefficients, whereas the number of columns is equal to \code{nlambda}. For \code{"mgaussian"} -family with m responses, it is a list of m such matrices.} -\item{iter}{A vector of length \code{nlambda} containing the number of -iterations until convergence at each value of \code{lambda}.} +family with m responses, it is a list of m such matrices.} +\item{iter}{A vector of length \code{nlambda} containing the number of +iterations until convergence at each value of \code{lambda}.} \item{lambda}{The sequence of regularization parameter values in the path.} \item{penalty}{Same as above.} \item{family}{Same as above.} -\item{alpha}{Same as above.} -\item{loss}{A vector containing either the residual sum of squares -(for \code{"gaussian", "mgaussian"}) or negative log-likelihood -(for \code{"binomial", "cox"}) of the fitted model at each value of \code{lambda}.} +\item{alpha}{Same as above.} +\item{loss}{A vector containing either the residual sum of squares +(for \verb{"gaussian", "mgaussian"}) or negative log-likelihood +(for \verb{"binomial", "cox"}) of the fitted model at each value of \code{lambda}.} \item{penalty.factor}{Same as above.} \item{n}{The number of observations used in the model fitting. It's equal to -\code{length(row.idx)}.} +\code{length(row.idx)}.} \item{center}{The sample mean vector of the variables, i.e., column mean of -the sub-matrix of \code{X} used for model fitting.} +the sub-matrix of \code{X} used for model fitting.} \item{scale}{The sample standard deviation of the variables, i.e., column -standard deviation of the sub-matrix of \code{X} used for model fitting.} +standard deviation of the sub-matrix of \code{X} used for model fitting.} \item{y}{The response vector used in the model fitting. Depending on \code{row.idx}, it could be a subset of the raw input of the response vector y.} -\item{screen}{Same as above.} +\item{screen}{Same as above.} \item{col.idx}{The indices of features that have 'scale' value greater than -1e-6. Features with 'scale' less than 1e-6 are removed from model fitting.} +1e-6. Features with 'scale' less than 1e-6 are removed from model fitting.} \item{rejections}{The number of features rejected at each value of \code{lambda}.} \item{safe_rejections}{The number of features rejected by safe rules at each value of \code{lambda}.} @@ -174,7 +174,7 @@ lasso, ridge, or elastic-net over a grid of values for the regularization parameter lambda. } \details{ -The objective function for linear regression or multiple responses linear regression +The objective function for linear regression or multiple responses linear regression (\code{family = "gaussian"} or \code{family = "mgaussian"}) is \deqn{\frac{1}{2n}\textrm{RSS} + \lambda*\textrm{penalty},}{(1/(2n))*RSS+ \lambda*penalty,} @@ -182,7 +182,7 @@ where for \code{family = "mgaussian"}), a group-lasso type penalty is applied. For logistic regression (\code{family = "binomial"}) it is \deqn{-\frac{1}{n} loglike + \lambda*\textrm{penalty},}{-(1/n)*loglike+\lambda*penalty}, for cox regression, - breslow approximation for ties is applied. +breslow approximation for ties is applied. Several advanced feature screening rules are implemented. For lasso-penalized linear regression, all the options of \code{screen} are @@ -192,7 +192,6 @@ data sets. For cox regression and/or the elastic net penalty, only \code{"SSR"} is applicable for now. More efficient rules are under development. } \examples{ - ## Linear regression data(colon) X <- colon$X @@ -245,12 +244,8 @@ plot(fit, main = "mgaussian") } \seealso{ -\code{\link{biglasso-package}}, \code{\link{setupX}}, -\code{\link{cv.biglasso}}, \code{\link{plot.biglasso}}, -\code{\link[ncvreg]{ncvreg}} +\link{biglasso-package}, \code{\link[=setupX]{setupX()}}, \code{\link[=cv.biglasso]{cv.biglasso()}}, \code{\link[=plot.biglasso]{plot.biglasso()}}, \code{\link[ncvreg:ncvreg]{ncvreg::ncvreg()}} } \author{ Yaohui Zeng, Chuyi Wang and Patrick Breheny - -Maintainer: Yaohui Zeng and Chuyi Wang } diff --git a/man/biglasso_fit.Rd b/man/biglasso_fit.Rd index df664b5..9f601a2 100644 --- a/man/biglasso_fit.Rd +++ b/man/biglasso_fit.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/biglasso_fit.R \name{biglasso_fit} \alias{biglasso_fit} -\title{Simplified call to biglasso: a gaussian model fit with no 'bells and whistles' (e.g., no SSR)} +\title{Direct interface to biglasso fitting, no preprocessing} \usage{ biglasso_fit( X, @@ -26,48 +26,48 @@ biglasso_fit( } \arguments{ \item{X}{The design matrix, without an intercept. It must be a -double type \code{\link[bigmemory]{big.matrix}} object.} +double type \code{\link[bigmemory:big.matrix]{bigmemory::big.matrix()}} object.} \item{y}{The response vector} -\item{r}{Residuals (length n vector) corresponding to `init`. -WARNING: If you supply an incorrect value of `r`, the +\item{r}{Residuals (length n vector) corresponding to \code{init}. +WARNING: If you supply an incorrect value of \code{r}, the solution will be incorrect.} \item{init}{Initial values for beta. Default: zero (length p vector)} -\item{xtx}{X scales: the jth element should equal `crossprod(X[,j])/n`. +\item{xtx}{X scales: the jth element should equal \code{crossprod(X[,j])/n}. In particular, if X is standardized, one should pass -`xtx = rep(1, p)`. WARNING: If you supply an incorrect value of -`xtx`, the solution will be incorrect. (length p vector)} +\code{xtx = rep(1, p)}. WARNING: If you supply an incorrect value of +\code{xtx}, the solution will be incorrect. (length p vector)} -\item{penalty}{String specifying which penalty to use. Default is 'lasso', +\item{penalty}{String specifying which penalty to use. Default is 'lasso', Other options are 'SCAD' and 'MCP' (the latter are non-convex)} \item{lambda}{A single value for the lasso tuning parameter.} \item{alpha}{The elastic-net mixing parameter that controls the relative -contribution from the lasso (l1) and the ridge (l2) penalty. +contribution from the lasso (l1) and the ridge (l2) penalty. The penalty is defined as: \deqn{ \alpha||\beta||_1 + (1-\alpha)/2||\beta||_2^2.} \code{alpha=1} is the lasso penalty, \code{alpha=0} the ridge penalty, \code{alpha} in between 0 and 1 is the elastic-net ("enet") penalty.} \item{gamma}{Tuning parameter value for nonconvex penalty. Defaults are -3.7 for `penalty = 'SCAD'` and 3 for `penalty = 'MCP'`} +3.7 for \code{penalty = 'SCAD'} and 3 for \code{penalty = 'MCP'}} \item{ncores}{The number of OpenMP threads used for parallel computing.} \item{max.iter}{Maximum number of iterations. Default is 1000.} \item{eps}{Convergence threshold for inner coordinate descent. The -algorithm iterates until the maximum change in the objective -after any coefficient update is less than \code{eps} times +algorithm iterates until the maximum change in the objective +after any coefficient update is less than \code{eps} times the null deviance. Default value is \code{1e-7}.} \item{dfmax}{Upper bound for the number of nonzero coefficients. Default is -no upper bound. However, for large data sets, -computational burden may be heavy for models with a large +no upper bound. However, for large data sets, +computational burden may be heavy for models with a large number of nonzero coefficients.} \item{penalty.factor}{A multiplicative factor for the penalty applied to @@ -85,54 +85,53 @@ fitting. Default is TRUE.} } \value{ An object with S3 class \code{"biglasso"} with following variables. -\item{beta}{The vector of estimated coefficients} -\item{iter}{A vector of length \code{nlambda} containing the number of -iterations until convergence} +\item{beta}{The vector of estimated coefficients} +\item{iter}{A vector of length \code{nlambda} containing the number of +iterations until convergence} \item{resid}{Vector of residuals calculated from estimated coefficients.} \item{lambda}{The sequence of regularization parameter values in the path.} -\item{alpha}{Same as in `biglasso()`} +\item{alpha}{Same as in \code{biglasso()}} \item{loss}{A vector containing either the residual sum of squares of the fitted model at each value of lambda.} -\item{penalty.factor}{Same as in `biglasso()`.} +\item{penalty.factor}{Same as in \code{biglasso()}.} \item{n}{The number of observations used in the model fitting.} \item{y}{The response vector used in the model fitting.} } \description{ -NOTE: this function is designed for users who have a strong understanding of -statistics and know exactly what they are doing. This is a simplification of -the main `biglasso()` function with more flexible settings. +This function is intended for users who know exactly what they're doing and +want complete control over the fitting process. It +\itemize{ +\item does NOT add an intercept +\item does NOT standardize the design matrix +\item does NOT set up a path for lambda (the lasso tuning parameter) +all of the above are critical steps in data analysis. However, a direct API +has been provided for use in situations where the lasso fitting process is +an internal component of a more complicated algorithm and standardization +must be handled externally. +} } \details{ -Of note, this function: - - * does NOT add an intercept - * does NOT standardize the design matrix - * does NOT set up a path for lambda (the lasso tuning parameter) - - all of the above are among the best practices for data analysis. This function - is made for use in situations where these steps have already been addressed prior - to model fitting. - - In other words, `biglasso_fit()` is to `biglasso()` what `ncvreg::ncvfit()` - is to `ncvreg::ncvreg()`. - - For now, this function only works with linear regression (`family = 'gaussian'`) +Note: +\itemize{ +\item Hybrid safe-strong rules are turned off for \code{biglasso_fit()}, as these rely +on standardization +\item Currently, the function only works with linear regression +(\code{family = 'gaussian'}). +} } \examples{ - data(Prostate) -X <- cbind(1, Prostate$X) |> ncvreg::std() # standardizing -> xtx is all 1s +X <- cbind(1, Prostate$X) +xtx <- apply(X, 2, crossprod)/nrow(X) y <- Prostate$y X.bm <- as.big.matrix(X) -init <- rep(0, ncol(X)) # using cold starts - will need more iterations -r <- y - X\%*\%init -fit <- biglasso_fit(X = X.bm, y = y, r = r, init = init, - xtx = rep(1, ncol(X)),lambda = 0.1, penalty.factor=c(0, rep(1, ncol(X)-1)), - max.iter = 10000) - -fit <- biglasso_fit(X = X.bm, y = y, r = r, init = init, penalty = 'MCP', - xtx = rep(1, ncol(X)), lambda = 0.005, penalty.factor=c(0, rep(1, ncol(X)-1)), - max.iter = 10000) +init <- rep(0, ncol(X)) +fit <- biglasso_fit(X = X.bm, y = y, r=y, init = init, xtx = xtx, + lambda = 0.1, penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter = 10000) +fit$beta +fit <- biglasso_fit(X = X.bm, y = y, r=y, init = init, xtx = xtx, penalty='MCP', + lambda = 0.1, penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter = 10000) +fit$beta } \author{ Tabitha Peter and Patrick Breheny diff --git a/man/biglasso_path.Rd b/man/biglasso_path.Rd index 973d66c..d60b4e4 100644 --- a/man/biglasso_path.Rd +++ b/man/biglasso_path.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/biglasso_path.R \name{biglasso_path} \alias{biglasso_path} -\title{Simplified call to biglasso: a gaussian model fit with no 'bells and whistles' (e.g., no SSR)} +\title{Direct interface to biglasso fitting, no preprocessing, path version} \usage{ biglasso_path( X, @@ -26,48 +26,48 @@ biglasso_path( } \arguments{ \item{X}{The design matrix, without an intercept. It must be a -double type \code{\link[bigmemory]{big.matrix}} object.} +double type \code{\link[bigmemory:big.matrix]{bigmemory::big.matrix()}} object.} \item{y}{The response vector} -\item{r}{Residuals (length n vector) corresponding to `init`. -WARNING: If you supply an incorrect value of `r`, the +\item{r}{Residuals (length n vector) corresponding to \code{init}. +WARNING: If you supply an incorrect value of \code{r}, the solution will be incorrect.} \item{init}{Initial values for beta. Default: zero (length p vector)} -\item{xtx}{X scales: the jth element should equal `crossprod(X[,j])/n`. +\item{xtx}{X scales: the jth element should equal \code{crossprod(X[,j])/n}. In particular, if X is standardized, one should pass -`xtx = rep(1, p)`. WARNING: If you supply an incorrect value of -`xtx`, the solution will be incorrect. (length p vector)} +\code{xtx = rep(1, p)}. WARNING: If you supply an incorrect value of +\code{xtx}, the solution will be incorrect. (length p vector)} -\item{penalty}{String specifying which penalty to use. Default is 'lasso', +\item{penalty}{String specifying which penalty to use. Default is 'lasso', Other options are 'SCAD' and 'MCP' (the latter are non-convex)} \item{lambda}{A vector of numeric values the lasso tuning parameter.} \item{alpha}{The elastic-net mixing parameter that controls the relative -contribution from the lasso (l1) and the ridge (l2) penalty. +contribution from the lasso (l1) and the ridge (l2) penalty. The penalty is defined as: \deqn{ \alpha||\beta||_1 + (1-\alpha)/2||\beta||_2^2.} \code{alpha=1} is the lasso penalty, \code{alpha=0} the ridge penalty, \code{alpha} in between 0 and 1 is the elastic-net ("enet") penalty.} \item{gamma}{Tuning parameter value for nonconvex penalty. Defaults are -3.7 for `penalty = 'SCAD'` and 3 for `penalty = 'MCP'`} +3.7 for \code{penalty = 'SCAD'} and 3 for \code{penalty = 'MCP'}} \item{ncores}{The number of OpenMP threads used for parallel computing.} \item{max.iter}{Maximum number of iterations. Default is 1000.} \item{eps}{Convergence threshold for inner coordinate descent. The -algorithm iterates until the maximum change in the objective -after any coefficient update is less than \code{eps} times +algorithm iterates until the maximum change in the objective +after any coefficient update is less than \code{eps} times the null deviance. Default value is \code{1e-7}.} \item{dfmax}{Upper bound for the number of nonzero coefficients. Default is -no upper bound. However, for large data sets, -computational burden may be heavy for models with a large +no upper bound. However, for large data sets, +computational burden may be heavy for models with a large number of nonzero coefficients.} \item{penalty.factor}{A multiplicative factor for the penalty applied to @@ -85,56 +85,57 @@ fitting. Default is TRUE.} } \value{ An object with S3 class \code{"biglasso"} with following variables. -\item{beta}{A sparse matrix where rows are estimates a given coefficient across all values of lambda} -\item{iter}{A vector of length \code{nlambda} containing the number of -iterations until convergence} +\item{beta}{A sparse matrix where rows are estimates a given coefficient across all values of lambda} +\item{iter}{A vector of length \code{nlambda} containing the number of +iterations until convergence} \item{resid}{Vector of residuals calculated from estimated coefficients.} \item{lambda}{The sequence of regularization parameter values in the path.} -\item{alpha}{Same as in `biglasso()`} +\item{alpha}{Same as in \code{biglasso()}} \item{loss}{A vector containing either the residual sum of squares of the fitted model at each value of lambda.} -\item{penalty.factor}{Same as in `biglasso()`.} +\item{penalty.factor}{Same as in \code{biglasso()}.} \item{n}{The number of observations used in the model fitting.} \item{y}{The response vector used in the model fitting.} } \description{ -NOTE: this function is designed for users who have a strong understanding of -statistics and know exactly what they are doing. This is a simplification of -the main `biglasso()` function with more flexible settings. +This function is intended for users who know exactly what they're doing and +want complete control over the fitting process. It +\itemize{ +\item does NOT add an intercept +\item does NOT standardize the design matrix +both of the above are critical steps in data analysis. However, a direct API +has been provided for use in situations where the lasso fitting process is +an internal component of a more complicated algorithm and standardization +must be handled externally. +} } \details{ -Of note, this function: - - * does NOT add an intercept - * does NOT standardize the design matrix - * does NOT set up a path for lambda (the lasso tuning parameter) automatically; - This vector of values must be user-supplied. - - This function is made for use in situations where these steps have already been addressed prior - to model fitting. - - In other words, `biglasso_path()` is doing the same thing as `biglasso_fit()`, - with the additional option to fit models across a path of tuning parameter values. - - For now, this function only works with linear regression (`family = 'gaussian'`) +\code{biglasso_path()} works identically to \code{\link[=biglasso_fit]{biglasso_fit()}} except it offers the +additional option of fitting models across a path of tuning parameter values. + +Note: +\itemize{ +\item Hybrid safe-strong rules are turned off for \code{biglasso_fit()}, as these rely +on standardization +\item Currently, the function only works with linear regression +(\code{family = 'gaussian'}). +} } \examples{ - data(Prostate) -X <- cbind(1, Prostate$X) |> ncvreg::std() # standardizing -> xtx is all 1s +X <- cbind(1, Prostate$X) +xtx <- apply(X, 2, crossprod)/nrow(X) y <- Prostate$y X.bm <- as.big.matrix(X) -init <- rep(0, ncol(X)) # using cold starts - will need more iterations -r <- y - X\%*\%init -fit_lasso <- biglasso_path(X = X.bm, y = y, r = r, init = init, - xtx = rep(1, ncol(X)), lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), - penalty.factor=c(0, rep(1, ncol(X)-1)), - max.iter = 10000) - -fit_mcp <- biglasso_path(X = X.bm, y = y, r = r, init = init, - xtx = rep(1, ncol(X)), lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), - penalty.factor=c(0, rep(1, ncol(X)-1)), - max.iter = 10000, penalty= 'MCP') +init <- rep(0, ncol(X)) +fit <- biglasso_path(X = X.bm, y = y, r = y, init = init, xtx = xtx, + lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), + penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter=2000) +fit$beta +fit <- biglasso_path(X = X.bm, y = y, r = y, init = init, xtx = xtx, + lambda = c(0.5, 0.1, 0.05, 0.01, 0.001), penalty='MCP', + penalty.factor=c(0, rep(1, ncol(X)-1)), max.iter = 2000) +fit$beta } \author{ Tabitha Peter and Patrick Breheny diff --git a/man/cv.biglasso.Rd b/man/cv.biglasso.Rd index 3be57f9..fd9b7f3 100644 --- a/man/cv.biglasso.Rd +++ b/man/cv.biglasso.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/cv.biglasso.R +% Please edit documentation in R/cv-biglasso.R \name{cv.biglasso} \alias{cv.biglasso} \title{Cross-validation for biglasso} @@ -20,8 +20,7 @@ cv.biglasso( ) } \arguments{ -\item{X}{The design matrix, without an intercept, as in -\code{\link{biglasso}}.} +\item{X}{The design matrix, without an intercept, as in \code{\link[=biglasso]{biglasso()}}.} \item{y}{The response vector, as in \code{biglasso}.} @@ -43,9 +42,9 @@ characteristic curve (ROC). \item{ncores}{The number of cores to use for parallel execution of the cross-validation folds, run on a cluster created by the \code{parallel} package. (This is also supplied to the \code{ncores} argument in -\code{\link{biglasso}}, which is the number of OpenMP threads, but only for -the first call of \code{\link{biglasso}} that is run on the entire data. The -individual calls of \code{\link{biglasso}} for the CV folds are run without +\code{\link[=biglasso]{biglasso()}}, which is the number of OpenMP threads, but only for +the first call of \code{\link[=biglasso]{biglasso()}} that is run on the entire data. The +individual calls of \code{\link[=biglasso]{biglasso()}} for the CV folds are run without the \code{ncores} argument.)} \item{...}{Additional arguments to \code{biglasso}.} @@ -67,20 +66,17 @@ when \code{eval.metric} is 'auc'.} } \value{ An object with S3 class \code{"cv.biglasso"} which inherits from -class \code{"cv.ncvreg"}. The following variables are contained in the -class (adopted from \code{\link[ncvreg]{cv.ncvreg}}). \item{cve}{The error -for each value of \code{lambda}, averaged across the cross-validation -folds.} \item{cvse}{The estimated standard error associated with each value -of for \code{cve}.} \item{lambda}{The sequence of regularization parameter -values along which the cross-validation error was calculated.} +class \code{"cv.ncvreg"}. The following variables are contained in the +class (adopted from \code{\link[ncvreg:cv.ncvreg]{ncvreg::cv.ncvreg()}}). +\item{cve}{The error for each value of \code{lambda}, averaged across the cross-validation folds.} +\item{cvse}{The estimated standard error associated with each value of for \code{cve}.} +\item{lambda}{The sequence of regularization parameter values along which the cross-validation error was calculated.} \item{fit}{The fitted \code{biglasso} object for the whole data.} \item{min}{The index of \code{lambda} corresponding to \code{lambda.min}.} -\item{lambda.min}{The value of \code{lambda} with the minimum -cross-validation error.} \item{lambda.1se}{The largest value of \code{lambda} -for which the cross-validation error is at most one standard error larger -than the minimum cross-validation error.} \item{null.dev}{The deviance for -the intercept-only model.} \item{pe}{If \code{family="binomial"}, the -cross-validation prediction error for each value of \code{lambda}.} +\item{lambda.min}{The value of \code{lambda} with the minimum cross-validation error.} +\item{lambda.1se}{The largest value of \code{lambda} for which the cross-validation error is at most one standard error larger than the minimum cross-validation error.} +\item{null.dev}{The deviance for the intercept-only model.} +\item{pe}{If \code{family="binomial"}, the cross-validation prediction error for each value of \code{lambda}.} \item{cv.ind}{Same as above.} } \description{ @@ -89,11 +85,12 @@ of values for the regularization parameter lambda. } \details{ The function calls \code{biglasso} \code{nfolds} times, each time leaving -out 1/\code{nfolds} of the data. The cross-validation error is based on the +out 1/\code{nfolds} of the data. The cross-validation error is based on the residual sum of squares when \code{family="gaussian"} and the binomial -deviance when \code{family="binomial"}.\cr \cr The S3 class object -\code{cv.biglasso} inherits class \code{\link[ncvreg]{cv.ncvreg}}. So S3 -functions such as \code{"summary", "plot"} can be directly applied to the +deviance when \code{family="binomial"}. + +The S3 class object \code{cv.biglasso} inherits class \code{\link[ncvreg:cv.ncvreg]{ncvreg::cv.ncvreg()}}. So S3 +functions such as \verb{"summary", "plot"} can be directly applied to the \code{cv.biglasso} object. } \examples{ @@ -110,14 +107,10 @@ par(mfrow = c(2, 2)) plot(cvfit, type = 'all') summary(cvfit) } - } \seealso{ -\code{\link{biglasso}}, \code{\link{plot.cv.biglasso}}, -\code{\link{summary.cv.biglasso}}, \code{\link{setupX}} +\code{\link[=biglasso]{biglasso()}}, \code{\link[=plot.cv.biglasso]{plot.cv.biglasso()}}, \code{\link[=summary.cv.biglasso]{summary.cv.biglasso()}}, \code{\link[=setupX]{setupX()}} } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } diff --git a/man/loss.biglasso.Rd b/man/loss.biglasso.Rd index 46b3c92..41f8959 100644 --- a/man/loss.biglasso.Rd +++ b/man/loss.biglasso.Rd @@ -28,12 +28,10 @@ characteristic curve (ROC).} Internal biglasso functions } \details{ -These are not intended for use by users.\code{loss.biglasso} calculates the +These are not intended for use by users. \code{loss.biglasso} calculates the value of the loss function for the given predictions (used for cross-validation). } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } \keyword{internal} diff --git a/man/plot.biglasso.Rd b/man/plot.biglasso.Rd index 318bf68..1173078 100644 --- a/man/plot.biglasso.Rd +++ b/man/plot.biglasso.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot.biglasso.R +% Please edit documentation in R/plot-biglasso.R \name{plot.biglasso} \alias{plot.biglasso} \title{Plot coefficients from a "biglasso" object} @@ -7,31 +7,24 @@ \method{plot}{biglasso}(x, alpha = 1, log.l = TRUE, ...) } \arguments{ -\item{x}{Fitted \code{"biglasso"} model.} +\item{x}{Fitted \code{\link[=biglasso]{biglasso()}} model.} \item{alpha}{Controls alpha-blending, helpful when the number of covariates is large. Default is alpha=1.} \item{log.l}{Should horizontal axis be on the log scale? Default is TRUE.} -\item{\dots}{Other graphical parameters to \code{plot}} +\item{\dots}{Other graphical parameters to \code{\link[=plot]{plot()}}} } \description{ -Produce a plot of the coefficient paths for a fitted \code{\link{biglasso}} -object. +Produce a plot of the coefficient paths for a fitted \code{\link[=biglasso]{biglasso()}} object. } \examples{ - ## See examples in "biglasso" - } \seealso{ -\code{\link{biglasso}}, \code{\link{cv.biglasso}} +\code{\link[=biglasso]{biglasso()}}, \code{\link[=cv.biglasso]{cv.biglasso()}} } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } -\keyword{models} -\keyword{regression} diff --git a/man/plot.cv.biglasso.Rd b/man/plot.cv.biglasso.Rd index a679a30..1100483 100644 --- a/man/plot.cv.biglasso.Rd +++ b/man/plot.cv.biglasso.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot.cv.biglasso.R +% Please edit documentation in R/plot-cv-biglasso.R \name{plot.cv.biglasso} \alias{plot.cv.biglasso} \title{Plots the cross-validation curve from a "cv.biglasso" object} @@ -39,26 +39,20 @@ the value where cross-validaton error is minimized.} \item{\dots}{Other graphical parameters to \code{plot}} } \description{ -Plot the cross-validation curve from a \code{\link{cv.biglasso}} object, +Plot the cross-validation curve from a \code{\link[=cv.biglasso]{cv.biglasso()}} object, along with standard error bars. } \details{ -Error bars representing approximate 68\% confidence intervals are plotted +Error bars representing approximate 68\\% confidence intervals are plotted along with the estimates at value of \code{lambda}. For \code{rsq} and \code{snr}, these confidence intervals are quite crude, especially near. } \examples{ - ## See examples in "cv.biglasso" - } \seealso{ -\code{\link{biglasso}}, \code{\link{cv.biglasso}} +\code{\link[=biglasso]{biglasso()}}, \code{\link[=cv.biglasso]{cv.biglasso()}} } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } -\keyword{models} -\keyword{regression} diff --git a/man/plot.mbiglasso.Rd b/man/plot.mbiglasso.Rd index 69cdc98..5d36e76 100644 --- a/man/plot.mbiglasso.Rd +++ b/man/plot.mbiglasso.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot.mbiglasso.R +% Please edit documentation in R/plot-mbiglasso.R \name{plot.mbiglasso} \alias{plot.mbiglasso} \title{Plot coefficients from a "mbiglasso" object} @@ -7,7 +7,7 @@ \method{plot}{mbiglasso}(x, alpha = 1, log.l = TRUE, norm.beta = TRUE, ...) } \arguments{ -\item{x}{Fitted \code{"mbiglasso"} model.} +\item{x}{Fitted \code{mbiglasso} model.} \item{alpha}{Controls alpha-blending, helpful when the number of covariates is large. Default is alpha=1.} @@ -17,24 +17,18 @@ is large. Default is alpha=1.} \item{norm.beta}{Should the vertical axis be the l2 norm of coefficients for each variable? Default is TRUE. If False, the vertical axis is the coefficients.} -\item{\dots}{Other graphical parameters to \code{plot}} +\item{\dots}{Other graphical parameters to \code{\link[=plot]{plot()}}} } \description{ Produce a plot of the coefficient paths for a fitted multiple responses \code{mbiglasso} object. } \examples{ - ## See examples in "biglasso" - } \seealso{ -\code{\link{biglasso}}, +\code{\link[=biglasso]{biglasso()}} } \author{ Chuyi Wang - -Maintainer: Chuyi Wang } -\keyword{models} -\keyword{regression} diff --git a/man/predict.biglasso.Rd b/man/predict.biglasso.Rd index d9c5449..86cda8a 100644 --- a/man/predict.biglasso.Rd +++ b/man/predict.biglasso.Rd @@ -36,20 +36,21 @@ \item{object}{A fitted \code{"biglasso"} model object.} \item{X}{Matrix of values at which predictions are to be made. It must be a -\code{\link[bigmemory]{big.matrix}} object. Not used for -\code{type="coefficients"}.} +\code{\link[bigmemory:big.matrix]{bigmemory::big.matrix()}} object. Not used for \code{type="coefficients"}.} -\item{row.idx}{Similar to that in \code{\link[biglasso]{biglasso}}, it's a +\item{row.idx}{Similar to that in \code{\link[=biglasso]{biglasso()}}, it's a vector of the row indices of \code{X} that used for the prediction. \code{1:nrow(X)} by default.} -\item{type}{Type of prediction: \code{"link"} returns the linear predictors; -\code{"response"} gives the fitted values; \code{"class"} returns the -binomial outcome with the highest probability; \code{"coefficients"} returns -the coefficients; \code{"vars"} returns a list containing the indices and -names of the nonzero variables at each value of \code{lambda}; -\code{"nvars"} returns the number of nonzero coefficients at each value of -\code{lambda}.} +\item{type}{Type of prediction: +\itemize{ +\item \code{"link"} returns the linear predictors +\item \code{"response"} gives the fitted values +\item \code{"class"} returns the binomial outcome with the highest probability +\item \code{"coefficients"} returns the coefficients +\item \code{"vars"} returns a list containing the indices and names of the nonzero variables at each value of \code{lambda} +\item \code{"nvars"} returns the number of nonzero coefficients at each value of \code{lambda} +}} \item{lambda}{Values of the regularization parameter \code{lambda} at which predictions are requested. Linear interpolation is used for values of @@ -75,8 +76,8 @@ coefficients. For \code{family="mgaussian"} only.} The object returned depends on \code{type}. } \description{ -Extract predictions (fitted reponse, coefficients, etc.) from a -fitted \code{\link{biglasso}} object. +Extract predictions (fitted reponse, coefficients, etc.) from a +fitted \code{\link[=biglasso]{biglasso()}} object. } \examples{ ## Logistic regression @@ -94,12 +95,8 @@ predict(fit, type="vars", lambda=c(0.05, 0.1)) predict(fit, type="nvars", lambda=c(0.05, 0.1)) } \seealso{ -\code{\link{biglasso}}, \code{\link{cv.biglasso}} +\code{\link[=biglasso]{biglasso()}}, \code{\link[=cv.biglasso]{cv.biglasso()}} } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } -\keyword{models} -\keyword{regression} diff --git a/man/predict.cv.biglasso.Rd b/man/predict.cv.biglasso.Rd index ea4cca0..705b9f7 100644 --- a/man/predict.cv.biglasso.Rd +++ b/man/predict.cv.biglasso.Rd @@ -1,9 +1,9 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/predict.cv.R +% Please edit documentation in R/predict-cv.R \name{predict.cv.biglasso} \alias{predict.cv.biglasso} \alias{coef.cv.biglasso} -\title{Model predictions based on a fitted \code{\link{cv.biglasso}} object} +\title{Model predictions based on a fitted \code{\link[=cv.biglasso]{cv.biglasso()}} object} \usage{ \method{predict}{cv.biglasso}( object, @@ -21,20 +21,22 @@ \item{object}{A fitted \code{"cv.biglasso"} model object.} \item{X}{Matrix of values at which predictions are to be made. It must be a -\code{\link[bigmemory]{big.matrix}} object. Not used for +\code{\link[bigmemory:big.matrix]{bigmemory::big.matrix()}} object. Not used for \code{type="coefficients"}.} -\item{row.idx}{Similar to that in \code{\link[biglasso]{biglasso}}, it's a +\item{row.idx}{Similar to that in \code{\link[=biglasso]{biglasso()}}, it's a vector of the row indices of \code{X} that used for the prediction. \code{1:nrow(X)} by default.} -\item{type}{Type of prediction: \code{"link"} returns the linear predictors; -\code{"response"} gives the fitted values; \code{"class"} returns the -binomial outcome with the highest probability; \code{"coefficients"} returns -the coefficients; \code{"vars"} returns a list containing the indices and -names of the nonzero variables at each value of \code{lambda}; -\code{"nvars"} returns the number of nonzero coefficients at each value of -\code{lambda}.} +\item{type}{Type of prediction: +\itemize{ +\item \code{"link"} returns the linear predictors +\item \code{"response"} gives the fitted values +\item \code{"class"} returns the binomial outcome with the highest probability +\item \code{"coefficients"} returns the coefficients +\item \code{"vars"} returns a list containing the indices and names of the nonzero variables at each value of \code{lambda} +\item \code{"nvars"} returns the number of nonzero coefficients at each value of \code{lambda} +}} \item{lambda}{Values of the regularization parameter \code{lambda} at which predictions are requested. The default value is the one corresponding to @@ -54,7 +56,7 @@ specified.} The object returned depends on \code{type}. } \description{ -Extract predictions from a fitted \code{\link{cv.biglasso}} object. +Extract predictions from a fitted \code{\link[=cv.biglasso]{cv.biglasso()}} object. } \examples{ \dontrun{ @@ -74,12 +76,8 @@ predict(cvfit, X.bm, lambda = "lambda.1se") } } \seealso{ -\code{\link{biglasso}}, \code{\link{cv.biglasso}} +\code{\link[=biglasso]{biglasso()}}, \code{\link[=cv.biglasso]{cv.biglasso()}} } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } -\keyword{models} -\keyword{regression} diff --git a/man/setupX.Rd b/man/setupX.Rd index 55ba902..89e15e3 100644 --- a/man/setupX.Rd +++ b/man/setupX.Rd @@ -35,12 +35,12 @@ file-backed \code{big.matrix}. By default, its name is the same as \item{type}{The data type. Only "double" is supported for now.} \item{...}{Additional arguments that can be passed into function -\code{\link[bigmemory]{read.big.matrix}}.} +\code{\link[bigmemory:write.big.matrix]{bigmemory::read.big.matrix()}}.} } \value{ A \code{big.matrix} object corresponding to a file-backed -\code{big.matrix}. It's ready to be used as the design matrix \code{X} in -\code{\link{biglasso}} and \code{\link{cv.biglasso}}. +\code{\link[bigmemory:big.matrix]{bigmemory::big.matrix()}}. It's ready to be used as the design matrix \code{X} in +\code{\link[=biglasso]{biglasso()}} and \code{\link[=cv.biglasso]{cv.biglasso()}}. } \description{ Set up the design matrix X as a \code{big.matrix} object based on external @@ -48,7 +48,7 @@ massive data file stored on disk that cannot be fullly loaded into memory. The data file must be a well-formated ASCII-file, and contains only one single type. Current version only supports \code{double} type. Other restrictions about the data file are described in -\code{\link{biglasso-package}}. This function reads the massive data, and +\link{biglasso-package}. This function reads the massive data, and creates a \code{big.matrix} object. By default, the resulting \code{big.matrix} is file-backed, and can be shared across processors or nodes of a cluster. @@ -59,21 +59,16 @@ For a data set, this function needs to be called only one time to set up the current working directory. Once set up, the data can be "loaded" into any (new) R session by calling \code{attach.big.matrix(discriptorfile)}. -This function is a simple wrapper of -\code{\link[bigmemory]{read.big.matrix}}. See -\code{\link[bigmemory]{read.big.matrix}} and the package +This function is a simple wrapper of \code{\link[bigmemory:write.big.matrix]{bigmemory::read.big.matrix()}}. See \href{https://CRAN.R-project.org/package=bigmemory}{bigmemory} for more details. } \examples{ ## see the example in "biglasso-package" - } \seealso{ -\code{\link{biglasso}}, \code{\link{cv.ncvreg}} +\code{\link[=biglasso]{biglasso()}}, \code{\link[=cv.ncvreg]{cv.ncvreg()}}, \link{biglasso-package} } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } diff --git a/man/summary.cv.biglasso.Rd b/man/summary.cv.biglasso.Rd index 26bca24..b1a220b 100644 --- a/man/summary.cv.biglasso.Rd +++ b/man/summary.cv.biglasso.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/summary.cv.biglasso.R +% Please edit documentation in R/summary-cv-biglasso.R \name{summary.cv.biglasso} \alias{summary.cv.biglasso} \alias{print.summary.cv.biglasso} @@ -22,19 +22,19 @@ non-integer printed values.} } \value{ \code{summary.cv.biglasso} produces an object with S3 class -\code{"summary.cv.biglasso"}. The class has its own print method and contains -the following list elements: \item{penalty}{The penalty used by -\code{biglasso}.} \item{model}{Either \code{"linear"} or \code{"logistic"}, -depending on the \code{family} option in \code{biglasso}.} \item{n}{Number -of observations} \item{p}{Number of regression coefficients (not including -the intercept).} \item{min}{The index of \code{lambda} with the smallest -cross-validation error.} \item{lambda}{The sequence of \code{lambda} values -used by \code{cv.biglasso}.} \item{cve}{Cross-validation error (deviance).} -\item{r.squared}{Proportion of variance explained by the model, as estimated -by cross-validation.} \item{snr}{Signal to noise ratio, as estimated by -cross-validation.} \item{sigma}{For linear regression models, the scale -parameter estimate.} \item{pe}{For logistic regression models, the -prediction error (misclassification error).} +\code{"summary.cv.biglasso"}. The class has its own print method and contains +the following list elements: +\item{penalty}{The penalty used by \code{biglasso}.} +\item{model}{Either \code{"linear"} or \code{"logistic"}, depending on the \code{family} option in \code{biglasso}.} +\item{n}{Number of observations} +\item{p}{Number of regression coefficients (not including the intercept).} +\item{min}{The index of \code{lambda} with the smallest cross-validation error.} +\item{lambda}{The sequence of \code{lambda} values used by \code{cv.biglasso}.} +\item{cve}{Cross-validation error (deviance).} +\item{r.squared}{Proportion of variance explained by the model, as estimated by cross-validation.} +\item{snr}{Signal to noise ratio, as estimated by cross-validation.} +\item{sigma}{For linear regression models, the scale parameter estimate.} +\item{pe}{For logistic regression models, the prediction error (misclassification error).} } \description{ Summary method for \code{cv.biglasso} objects. @@ -43,13 +43,8 @@ Summary method for \code{cv.biglasso} objects. ## See examples in "cv.biglasso" and "biglasso-package" } \seealso{ -\code{\link{biglasso}}, \code{\link{cv.biglasso}}, -\code{\link{plot.cv.biglasso}} +\code{\link[=biglasso]{biglasso()}}, \code{\link[=cv.biglasso]{cv.biglasso()}}, \code{\link[=plot.cv.biglasso]{plot.cv.biglasso()}}, \link{biglasso-package} } \author{ Yaohui Zeng and Patrick Breheny - -Maintainer: Yaohui Zeng } -\keyword{models} -\keyword{regression} diff --git a/src/Makevars b/src/Makevars index aea00f0..3f8523a 100644 --- a/src/Makevars +++ b/src/Makevars @@ -5,11 +5,9 @@ #PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) #PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) -#CXX_STD = CXX11 #PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) #PKG_LIBS = $(SHLIB_OPENMP_CFLAGS) -# CXX_STD = CXX11 # CXX11FLAGS=-O0 -g # Uncomment this FOR DEBUGGING ONLY PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/vignettes/biglasso.Rmd b/vignettes/biglasso.Rmd index ab42f77..e570254 100644 --- a/vignettes/biglasso.Rmd +++ b/vignettes/biglasso.Rmd @@ -16,11 +16,9 @@ knitr::opts_chunk$set( ) ``` -# 1 User Guide +# Small Data -## 1.1 Small Data - -### 1.1.1 Standar Lasso +## Linear regression ```{r} library(biglasso) @@ -42,7 +40,6 @@ X.bm <- as.big.matrix(X) str(X.bm) ``` - ```{r} dim(X.bm) ``` @@ -52,15 +49,14 @@ X.bm[1:5, 1:5] ## same results as X[1:5, 1:5] ``` - - ```{r} ## fit entire solution path, using our newly proposed "Adaptive" screening rule (default) fit <- biglasso(X.bm, y) plot(fit) ``` -### 1.1.2 Cross-Validation +## Cross-Validation + ```{r} ## 10-fold cross-valiation in parallel cvfit <- tryCatch( @@ -95,7 +91,7 @@ summary(cvfit) coef(cvfit)[which(coef(cvfit) != 0),] ``` -### 1.1.3 Logistic Regression +## Logistic Regression ```{r} data(Heart) @@ -106,7 +102,7 @@ fit <- biglasso(X.bm, y, family = "binomial") plot(fit) ``` -### 1.1.4 Cox Regression +## Cox Regression ```{r} library(survival) @@ -117,7 +113,7 @@ fit <- biglasso(X.bm, y, family = "cox") plot(fit) ``` -### 1.1.5 Multiple responses Linear Regression +## Multiple response Linear Regression (multi-task learning) ```{r} set.seed(10101) @@ -130,7 +126,7 @@ fit = biglasso(x.bm, y, family = "mgaussian") plot(fit) ``` -## 1.2 Big Data +# Big Data When the raw data file is very large, it's better to convert the raw data file into a file-backed `big.matrix` by using a file cache. We can call function @@ -191,9 +187,9 @@ par(mfrow = c(2, 2), mar = c(3.5, 3.5, 3, 1), mgp = c(2.5, 0.5, 0)) plot(cvfit, type = "all") ``` -# 2 Useful Reference +# Links -* biglasso R manual: https://cran.r-project.org/package=biglasso/biglasso.pdf -* biglasso on GitHub: https://github.com/YaohuiZeng/biglasso -* biglasso website: https://yaohuizeng.github.io/biglasso/index.html -* big.matrix manipulation: https://cran.r-project.org/package=bigmemory/index.html +* [biglasso on CRAN](https://cran.r-project.org/package=biglasso) +* [biglasso on GitHub](https://github.com/pbreheny/biglasso) +* [biglasso website](https://pbreheny.github.io/biglasso/index.html) +* [big.matrix manipulation](https://cran.r-project.org/package=bigmemory/index.html)