CHG: refine the r code by enhancing readbility, reformatting, and red…

…ucing redundant arguments
egpivo · Jan 12, 2024 · ac0a6e6 · ac0a6e6
1 parent e91be26
commit ac0a6e6
Show file tree

Hide file tree

Showing 7 changed files with 102 additions and 175 deletions.
diff --git a/R/helper.R b/R/helper.R
@@ -1,12 +1,13 @@
 #' Orthogonalized B-splines
+#'
 #' @param knots Array. The knots that define the spline.
 #' @param boundary_knots Array. The breakpoints that define the spline.
 #' @param degree Integer. The degree of the piecewise polynomial.
 #' @param predictors Array. The predictor variables with size p.
 #' @param is_approx Boolean. The default is `FALSE`.
 #' @return A list containing:
-#'   \item{bsplines}{Matrix of orthogonalized B-splines with dimension (p, length(knots) + degree + 1)}
-#'   \item{z}{Predictors used in generation}
+#'   \item{\code{bsplines}}{Matrix of orthogonalized B-splines with dimensions \eqn{(p, \text{length}(knots) + \text{degree} + 1)}}
+#'   \item{\code{z}}{Predictors used in generation}
 #' @export
 #' @examples
 #' p <- 30
@@ -83,39 +84,38 @@ orthogonize_bspline <- function(
   ))
 }
 
-
-#' Internal function: Validate new locations for a qrglasso_object object
+#' Internal function: Validate Parameters for Predicting with a `qrglasso` Object
 #'
 #' @keywords internal
-#' @param qrglasso_object An `qrglasso` class object.
-#' @param metric_type Character. A metric type for gamma selection. e.g., `BIC`, `BIC-log`. Default is `BIC`.
-#' @param top_k Integer. A matrix of the top K estimated functions.
+#' @param qrglasso_object A `qrglasso` class object.
+#' @param metric_type Character. Metric type for gamma selection, e.g., `BIC`, `BIC-log`. Default is `BIC`.
+#' @param top_k Integer. Top K estimated functions.
 #' @param degree Integer. Degree of the piecewise polynomial.
 #' @param boundaries Array. Two boundary points.
 #' @return `NULL`.
 #'
 check_predict_parameters <- function(qrglasso_object, metric_type, top_k, degree, boundaries) {
   if (!inherits(qrglasso_object, "qrglasso")) {
-    stop("Invalid object! Please enter a `qrglasso` object")
+    stop("Invalid object! Please enter a `qrglasso` object.")
   }
   if (!(metric_type %in% c("BIC", "BIC-log"))) {
-    stop("Only accept types: `BIC` and `BIC-log`")
+    stop("Only accept types: `BIC` and `BIC-log`.")
   }
   if (top_k <= 0) {
-    stop("Please enter a positive top k")
+    stop("Please enter a positive top k.")
   }
   if (degree <= 0) {
-    stop("Please enter a positive degree")
+    stop("Please enter a positive degree.")
   }
   if (length(boundaries) != 2) {
     stop("Please enter a size 2 boundaries.")
   }
   if (boundaries[1] >= boundaries[2]) {
     stop("Please input valid boundaries consisting of two elements in ascending order.")
   }
-  total_knots = qrglasso_object$L - degree + 1
+  total_knots <- qrglasso_object$L - degree + 1
   if (total_knots <= 0) {
-    stop("Please enter a smaller degree")
+    stop("Please enter a smaller degree.")
   }
 }
 
@@ -125,8 +125,8 @@ check_predict_parameters <- function(qrglasso_object, metric_type, top_k, degree
 #' @return `NULL`
 #' 
 plot_sequentially <- function(objs) {
-  originalPar <- par(no.readonly = TRUE)
-  on.exit(par(originalPar))
+  original_par <- par(no.readonly = TRUE)
+  on.exit(par(original_par))
   par(ask = TRUE)
   suppressWarnings({
     for (obj in objs) {
@@ -138,8 +138,8 @@ plot_sequentially <- function(objs) {
 
 #' Internal function: Plot Coefficient Function
 #' @keywords internal
-#' @param data A dataframe contains columns ``z``, ``coefficient``
-#' @param variate A character represent the title
+#' @param data A dataframe containing columns ``z``, ``coefficient``
+#' @param variate A character representing the title
 #' @return A ggplot object
 plot_coefficient_function <- function(data, variate) {
   default_theme <- theme_classic() +
@@ -148,17 +148,16 @@ plot_coefficient_function <- function(data, variate) {
       plot.title = element_text(hjust = 0.5)
     )
   result <- ggplot(data, aes(x = z, y = coefficient)) +
-    geom_point(col="#4634eb") +
+    geom_point(col = "#4634eb") +
     ggtitle(variate) +
     default_theme
   return(result)
 }
 
-
 #' Internal function: Plot BIC Results w.r.t. lambda
 #' @keywords internal
-#' @param data A dataframe contains columns ``lambda``, ``bic``
-#' @param variate A character represent the title
+#' @param data A dataframe containing columns ``lambda``, ``bic``
+#' @param variate A character representing the title
 #' @return A ggplot object
 plot_bic_result <- function(data, variate) {
   default_theme <- theme_classic() +
@@ -167,12 +166,11 @@ plot_bic_result <- function(data, variate) {
       plot.title = element_text(hjust = 0.5)
     )
   result <- ggplot(data, aes(x = lambda, y = bic)) +
-    geom_point(col="#4634eb") +
+    geom_point(col = "#4634eb") +
     geom_line() +
     ggtitle(variate) +
     xlab(expression(lambda)) +
     ylab("BIC") +
     default_theme
   return(result)
 }
-
diff --git a/R/qrglasso.R b/R/qrglasso.R
@@ -1,11 +1,10 @@
-#' Adaptively weighted group Lasso
+#' @title Adaptively Weighted Group Lasso
 #'
-#' @param Y A  \eqn{n \times 1} data matrix
-#' @param W A  \eqn{n \times pL}B-spline matrix.
-#' @param L The number of groups
-#' @param omega A $p x 1$ weight matrix. Default value is NULL.
+#' @param Y A \eqn{n \times 1} data matrix where \eqn{n} is the sample size.
+#' @param W A \eqn{n \times (p \times L) } B-spline matrix where \eqn{L} is the number of groups and \eqn{p} is the number of covariates.
+#' @param omega A \eqn{p \times 1} weight matrix. Default value is NULL.
 #' @param tau A numeric quantile of interest. Default value is 0.5.
-#' @param qn  A nuneric bound parameter for HDIC.  Default value is 1.
+#' @param qn A numeric bound parameter for HDIC. Default value is 1.
 #' @param lambda A sequence of tuning parameters. Default value is NULL.
 #' @param maxit The maximum number of iterations. Default value is 1000.
 #' @param thr Threshold for convergence. Default value is \eqn{10^{-4}}.
@@ -15,6 +14,7 @@
 #'   \item{\code{phi}}{An auxiliary estimate in the ADMM algorithm.}
 #'   \item{\code{BIC}}{A sequence of BIC values with respect to different lambdas.}
 #'   \item{\code{lambda}}{A sequence of tuning parameters used in the algorithm.}
+#'   \item{\code{L}}{The number of groups.}
 #'   \item{\code{omega}}{A \eqn{p \times 1} weight matrix used in the algorithm.}
 #' @author Wen-Ting Wang
 #' @references Toshio Honda, Ching-Kang Ing, Wei-Ying Wu (2019). Adaptively weighted group Lasso for semiparametric quantile regression models. \emph{Bernoulli} \bold{225} 4B.
@@ -56,104 +56,69 @@
 #' }
 #' 
 #' # Perform quantile regression with group Lasso
-#' result <- qrglasso(as.matrix(y), W, L - 1)
-#' 
-#' 
-#' # Extract relevant information from the result
-#' approx_bsplines <- orthogonize_bspline(knots, boundaries, degree)
-#' bsplines <- approx_bsplines$bsplines[, -1]
-#' z <- approx_bsplines$z
-#' L_star <- L - 1
-#' gamma_hat <- result$gamma[, which.min(result$BIC[, 1])]
-#' g1_BIC <- bsplines %*% gamma_hat[1:L_star]
-#' g2_BIC <- bsplines %*% gamma_hat[(L_star + 1):(2 * L_star)]
-#' 
-#' g10 <- bsplines %*% result$gamma[, 1][1:L_star]
-#' g20 <- bsplines %*% result$gamma[, 1][(L_star + 1):(2 * L_star)]
-#' 
-#' # Plotting
-#' original_par <- par(no.readonly = TRUE)
-#' par(mfrow = c(1, 2))
-#' 
-#' # Plot for g1(z)
-#' plot(z, g1(z), type = 'l', lwd = 4, cex.lab = 2, lty = 2, cex.axis = 1.5, ylim = c(-3, 3))
-#' lines(z, g10, col = 2, lwd = 2)
-#' lines(z, g1_BIC, col = 3, lwd = 2)
-#' legend("topleft", c("True", expression(lambda == 0), expression(lambda[BIC])), col = c(1, 2, 3), 
-#'        lty = c(2, 1, 1), lwd = c(3, 3, 3), horiz = TRUE, bty = "n")
-#'
-#' # Plot for g2(z)
-#' plot(z, rep(0, length(z)), type = 'l', lwd = 4, cex.lab = 2, lty = 2, cex.axis = 1.5, 
-#'      ylab = "g2(z)", ylim = c(-0.1, 0.1))
-#' lines(z, g20, col = 2, lwd = 2)
-#' lines(z, g2_BIC, col = 3, lwd = 2)
-#' legend("topleft", c("True", expression(lambda == 0), expression(lambda[BIC])), col = c(1, 2, 3, 4), 
-#'        lty = c(2, 1, 1, 1), lwd = c(3, 3, 3, 3), horiz = TRUE, bty = "n")
-#'        
-#' par(original_par)
+#' result <- qrglasso(as.matrix(y), W)
+#' # BIC Results
+#' plot(result)
+#' # Prediction
+#' estimate = predict(result)
+#' plot(estimate)
 #' 
-qrglasso <-
-  function(Y,
-           W,
-           L,
-           omega = NULL,
-           tau = 0.5,
-           qn = 1,
-           lambda = NULL,
-           maxit = 1000,
-           thr = 1e-04) {
-    if (is.null(lambda)) {
-      nlambda <- 51
-      max.lambda <- 10
-      lambda <-
-        c(0, exp(seq(
-          log(max.lambda / 1e4), log(max.lambda), length = (nlambda - 1)
-        )))
-    } else {
-      nlambda <- length(lambda)
-    }
-
-    zeta <- 10
-    zetaincre <- 1
-
-    if (is.null(omega))
-      result <-
-      awgl(Y, W, lambda, tau, L, qn, zeta, zetaincre, maxit, thr)
-    else
-      result <-
-      awgl_omega(Y, W, omega, lambda, tau, qn, zeta, zetaincre, maxit, thr)
-
-    result$phi[, 1] <- result$gamma[, 1]
-
-    obj.cv <- list(
-      gamma = result$gamma,
-      xi = result$xi,
-      phi = result$phi,
-      BIC = result$BIC,
-      lambda = lambda,
-      L = L,
-      omega = result$omega
-    )
-
-    class(obj.cv) <- "qrglasso"
-    return(obj.cv)
+#' @export
+qrglasso <- function(Y,
+                     W,
+                     omega = NULL,
+                     tau = 0.5,
+                     qn = 1,
+                     lambda = NULL,
+                     maxit = 1000,
+                     thr = 1e-04) {
+  if (is.null(lambda)) {
+    nlambda <- 51
+    max.lambda <- 10
+    lambda <- c(0, exp(seq(log(max.lambda / 1e4), log(max.lambda), length = (nlambda - 1))))
+  } else {
+    nlambda <- length(lambda)
   }
+
+  zeta <- 10
+  zetaincre <- 1
+  L_star <- dim(W)[2] / dim(W)[1]
+
+  if (is.null(omega))
+    result <- awgl(Y, W, lambda, tau, L_star, qn, zeta, zetaincre, maxit, thr)
+  else
+    result <- awgl_omega(Y, W, omega, lambda, tau, qn, zeta, zetaincre, maxit, thr)
+
+  result$phi[, 1] <- result$gamma[, 1]
+
+  obj.cv <- list(
+    gamma = result$gamma,
+    xi = result$xi,
+    phi = result$phi,
+    BIC = result$BIC,
+    lambda = lambda,
+    L = L_star + 1,
+    omega = result$omega
+  )
+
+  class(obj.cv) <- "qrglasso"
+  return(obj.cv)
+}
 
 
-
-#' @title Predict the coefficient functions
+#' @title Predict Coefficient Functions
 #'
 #' @description Predict the top-k coefficient functions
 #'
-#' @param qrglasso_object An  \code{qrglasso} class object.
-#' @param metric_type Character. A metric type for gamma selection. e.g., `BIC`, `BIC-log`. Default is `BIC`.
+#' @param qrglasso_object A \code{qrglasso} class object.
+#' @param metric_type Character. Metric type for gamma selection, e.g., `BIC`, `BIC-log`. Default is `BIC`.
 #' @param top_k Integer. A matrix of the top K estimated functions. Default is 5.
 #' @param degree Integer. Degree of the piecewise polynomial. Default is 2.
 #' @param boundaries Array. Two boundary points. Default is c(0, 1).
 #' @param is_approx Logical. If TRUE, the size of covariate indexes will be 1e6; otherwise, 1e4. Default is FALSE.
 #' @seealso \code{\link{qrglasso}}
 #' @return A list containing:
-#'   \item{coef_functions}{Matrix. Top-k coefficient function estimates with dimenstion (\eqn{m \times k}) where $m$ is size of `z`.}
+#'   \item{coef_functions}{Matrix. Top-k coefficient function estimates with dimension (\eqn{m \times k}) where $m$ is the size of `z`.}
 #'   \item{z}{Array. Index predictors used in generation}
 #' @examples
 #' set.seed(123)
@@ -164,7 +129,7 @@ qrglasso <-
 #' W <- matrix(rnorm(n * p * (L - 1)), n, p * (L - 1))
 #'
 #' # Call qrglasso with default parameters
-#' result <- qrglasso(Y = Y, W = W, L = 5)
+#' result <- qrglasso(Y = Y, W = W)
 #' estimate <- predict(result) 
 #' print(dim(estimate$coef_functions))
 #' 
@@ -215,29 +180,28 @@ predict <- function(qrglasso_object,
 #' Y <- matrix(rnorm(n), n, 1)
 #' W <- matrix(rnorm(n * p * (L - 1)), n, p * (L - 1))
 #'
-#' result <- qrglasso(Y = Y, W = W, L = 5)
+#' result <- qrglasso(Y = Y, W = W)
 #' plot(result)
 #'
 plot.qrglasso <- function(x, ...) {
   if (!inherits(x, "qrglasso")) {
-    stop("Invalid object! Please enter a `qrglasso` object")
+    stop("Invalid object! Please enter a `qrglasso` object.")
   }
   originalPar <- par(no.readonly = TRUE)
   result <- list()
   variates <- c("BIC", "BIC-log")
   for (i in 1:2) {
     variate <- variates[i]
-    data <- data.frame(lambda = x$lambda, bic = x$BIC[,i])
+    data <- data.frame(lambda = x$lambda, bic = x$BIC[, i])
     result[[variate]] <- plot_bic_result(data, variate)
   }
   plot_sequentially(result)
   par(originalPar)
 }
 
-
-#' @title  Display the estimated coefficient functions
+#' @title  Display Predicted Coefficient Functions
 #'
-#' @description Display the estimated coefficient functions by BIC
+#' @description Display the predicted coefficient functions by BIC
 #'
 #' @param x An object of class \code{qrglasso.predict} for the \code{plot} method
 #' @param ... Not used directly
@@ -250,27 +214,25 @@ plot.qrglasso <- function(x, ...) {
 #' set.seed(123)
 #' n <- 100
 #' p <- 5
-#' L <- 5
 #' Y <- matrix(rnorm(n), n, 1)
 #' W <- matrix(rnorm(n * p * (L - 1)), n, p * (L - 1))
 #'
-#' result <- qrglasso(Y = Y, W = W, L = 5)
+#' result <- qrglasso(Y = Y, W = W)
 #' estimate <- predict(result, top_k = 2)
 #' plot(estimate)
 #'
 plot.qrglasso.predict <- function(x, ...) {
   if (!inherits(x, "qrglasso.predict")) {
-    stop("Invalid object! Please enter a `qrglasso.predict` object")
+    stop("Invalid object! Please enter a `qrglasso.predict` object.")
   }
   originalPar <- par(no.readonly = TRUE)
   k <- dim(x$coef_functions)[2]
   result <- list()
   for (i in 1:k) {
     variate <- paste0("Coefficient function - g", i)
-    data <- data.frame(z = x$z, coefficient = x$coef_functions[,i])
+    data <- data.frame(z = x$z, coefficient = x$coef_functions[, i])
     result[[variate]] <- plot_coefficient_function(data, variate)
   }
   plot_sequentially(result)
   par(originalPar)
 }
-
diff --git a/man/plot.qrglasso.Rd b/man/plot.qrglasso.Rd
diff --git a/man/plot.qrglasso.predict.Rd b/man/plot.qrglasso.predict.Rd