diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml
new file mode 100644
index 0000000..d3e57f1
--- /dev/null
+++ b/.github/workflows/lint.yaml
@@ -0,0 +1,78 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+
+name: lint
+
+jobs:
+  lint:
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    steps:
+      - uses: actions/checkout@v3
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::lintr, local::.
+          needs: lint
+
+      - name: Lint
+        run: |
+          library(lintr)
+          critical_rules <- c(
+            absolute_path_linter(), any_duplicated_linter(),
+            any_is_na_linter(), assignment_linter(), backport_linter(),
+            brace_linter(), class_equals_linter(), commas_linter(),
+            commented_code_linter(), condition_message_linter(),
+            conjunct_test_linter(), consecutive_assertion_linter(),
+            duplicate_argument_linter(), empty_assignment_linter(),
+            equals_na_linter(), expect_comparison_linter(),
+            expect_identical_linter(), expect_not_linter(),
+            expect_null_linter(), expect_s3_class_linter(),
+            expect_s4_class_linter(), expect_true_false_linter(),
+            expect_type_linter(), fixed_regex_linter(),
+            for_loop_index_linter(), function_left_parentheses_linter(),
+            function_return_linter(), if_not_else_linter(),
+            ifelse_censor_linter(), implicit_assignment_linter(),
+            indentation_linter(), infix_spaces_linter(),
+            inner_combine_linter(), is_numeric_linter(),
+            keyword_quote_linter(), length_levels_linter(),
+            length_test_linter(), lengths_linter(),
+            library_call_linter(), literal_coercion_linter(),
+            missing_argument_linter(), missing_package_linter(),
+            namespace_linter(), nested_ifelse_linter(),
+            nonportable_path_linter(), numeric_leading_zero_linter(),
+            object_length_linter(), object_usage_linter(),
+            outer_negation_linter(), package_hooks_linter(),
+            paren_body_linter(), paste_linter(), pipe_call_linter(),
+            pipe_consistency_linter(), pipe_continuation_linter(),
+            quotes_linter(), redundant_equals_linter(),
+            redundant_ifelse_linter(), regex_subset_linter(),
+            repeat_linter(), routine_registration_linter(),
+            scalar_in_linter(), semicolon_linter(), sort_linter(),
+            spaces_inside_linter(), spaces_left_parentheses_linter(),
+            sprintf_linter(), string_boundary_linter(),
+            strings_as_factors_linter(), system_file_linter(),
+            T_and_F_symbol_linter(), undesirable_function_linter(),
+            undesirable_operator_linter(),
+            unnecessary_concatenation_linter(), unnecessary_lambda_linter(),
+            unnecessary_nested_if_linter(),
+            unnecessary_placeholder_linter(), unreachable_code_linter(),
+            unused_import_linter(), vector_logic_linter(),
+            whitespace_linter(), yoda_test_linter()
+          )
+
+          lint_package(linters = critical_rules, show_progress = FALSE)
+
+        shell: Rscript {0}
+        env:
+          LINTR_ERROR_ON_LINT: true
diff --git a/DESCRIPTION b/DESCRIPTION
index ee1ff72..911a4df 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -46,7 +46,8 @@ Imports:
 LinkingTo: Rcpp, RcppArmadillo
 Suggests:
     testthat (>= 3.0.0),
-    pracma
+    pracma,
+    lintr
 License: GPL-3
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
diff --git a/R/MADMMplasso.R b/R/MADMMplasso.R
index e90b2b2..9167ba5 100644
--- a/R/MADMMplasso.R
+++ b/R/MADMMplasso.R
@@ -50,7 +50,7 @@
 
 #' @example inst/examples/MADMMplasso_example.R
 #' @export
-MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max_it = 50000, e.abs = 1E-3, e.rel = 1E-3, maxgrid, nlambda, rho = 5, my_print = F, alph = 1.8, tree, parallel = T, pal = 0, gg = NULL, tol = 1E-4, cl = 4, legacy = FALSE) {
+MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = 0.001, max_it = 50000, e.abs = 1E-3, e.rel = 1E-3, maxgrid, nlambda, rho = 5, my_print = FALSE, alph = 1.8, tree, parallel = TRUE, pal = 0, gg = NULL, tol = 1E-4, cl = 4, legacy = FALSE) {
   N <- nrow(X)
 
   p <- ncol(X)
@@ -93,7 +93,6 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max
   rat <- lambda_min
 
   if (is.null(my_lambda)) {
-    lamda_new <- matrix(0, dim(y)[2])
     r <- y
 
     lammax <- lapply(
@@ -126,13 +125,9 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max
   }
 
   lam_list <- list()
-  beta_0_list <- list()
-  theta_0_list <- list()
-  beta_list <- list()
-  theta_list <- list()
-  obj <- c()
-  n_main_terms <- c()
-  non_zero_theta <- c()
+  obj <- NULL
+  n_main_terms <- NULL
+  non_zero_theta <- NULL
   my_obj <- list()
 
   my_W_hat <- generate_my_w(X = X, Z = Z)
@@ -152,7 +147,7 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max
   I <- matrix(0, nrow = nrow(C) * dim(y)[2], ncol = dim(y)[2])
 
   II <- input[multiple_of_D]
-  diag(I[c(1:dim(y)[2]), ]) <- C[1, ] * (CW[1])
+  diag(I[1:dim(y)[2], ]) <- C[1, ] * (CW[1])
 
   c_count <- 2
   for (e in II[-length(II)]) {
@@ -161,10 +156,10 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max
   }
   new_I <- diag(t(I) %*% I)
   new_G <- matrix(0, (p + p * K))
-  new_G[c(1:p)] <- 1
-  new_G[-c(1:p)] <- 2
-  new_G[c(1:p)] <- rho * (1 + new_G[c(1:p)])
-  new_G[-c(1:p)] <- rho * (1 + new_G[-c(1:p)])
+  new_G[1:p] <- 1
+  new_G[-1:-p] <- 2
+  new_G[1:p] <- rho * (1 + new_G[1:p])
+  new_G[-1:-p] <- rho * (1 + new_G[-1:-p])
 
   invmat <- list() # denominator of the beta estimates
   for (rr in 1:D) {
@@ -211,29 +206,25 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max
       admm_MADMMplasso(
         beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat, XtY,
         y, N, e.abs, e.rel, alpha, lam[i, ], alph, svd.w, tree, my_print,
-        invmat, gg[i, ],legacy
+        invmat, gg[i, ], legacy
       )
     }
     parallel::stopCluster(cl)
-  } else if (parallel == F & pal == 0) {
+  } else if (parallel && pal == 0) {
     my_values <- lapply(
       seq_len(nlambda),
       function(g) {
         admm_MADMMplasso(
           beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat,
           XtY, y, N, e.abs, e.rel, alpha, lam[g, ], alph, svd.w, tree, my_print,
-          invmat, gg[g, ],legacy
+          invmat, gg[g, ], legacy
         )
       }
     )
   }
 
-  repeat_loop <- 0
   hh <- 1
   while (hh <= nlambda) {
-    res_dual <- 0 # dual residual
-    res_pri <- 0 # primal residual
-
     lambda <- lam[hh, ]
 
     start_time <- Sys.time()
@@ -241,12 +232,11 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max
       my_values <- admm_MADMMplasso(
         beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat, XtY,
         y, N, e.abs, e.rel, alpha, lambda, alph, svd.w, tree, my_print, invmat,
-        gg[hh, ],legacy
+        gg[hh, ], legacy
       )
 
       beta <- my_values$beta
       theta <- my_values$theta
-      converge <- my_values$converge
       my_obj[[hh]] <- list(my_values$obj)
       beta0 <- my_values$beta0
       theta0 <- my_values$theta0 ### iteration
@@ -255,19 +245,17 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = .001, max
     }
     cost_time <- Sys.time() - start_time
     print(cost_time)
-    if (parallel == T & pal == 0) {
+    if (parallel && pal == 0) {
       beta <- my_values[hh, ]$beta
       theta <- my_values[hh, ]$theta
-      converge <- my_values[hh, ]$converge
       my_obj[[hh]] <- list(my_values[hh, ]$obj)
       beta0 <- my_values[hh, ]$beta0
       theta0 <- my_values[hh, ]$theta0 ### iteration
       beta_hat <- my_values[hh, ]$beta_hat
       y_hat <- my_values[hh, ]$y_hat
-    } else if (parallel == F & pal == 0) {
+    } else if (parallel && pal == 0) {
       beta <- my_values[[hh]]$beta
       theta <- my_values[[hh]]$theta
-      converge <- my_values[[hh]]$converge
       my_obj[[hh]] <- list(my_values[[hh]]$obj)
       beta0 <- my_values[[hh]]$beta0
       theta0 <- my_values[[hh]]$theta0 ### iteration
diff --git a/R/admm_MADMMplasso.R b/R/admm_MADMMplasso.R
index d09c37a..6bb0775 100644
--- a/R/admm_MADMMplasso.R
+++ b/R/admm_MADMMplasso.R
@@ -48,7 +48,7 @@
 
 
 #' @export
-admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, W_hat, XtY, y, N, e.abs, e.rel, alpha, lambda, alph, svd.w, tree, my_print = T, invmat, gg = 0.2, legacy = FALSE) {
+admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, W_hat, XtY, y, N, e.abs, e.rel, alpha, lambda, alph, svd.w, tree, my_print = TRUE, invmat, gg = 0.2, legacy = FALSE) {
   if (!legacy) {
     out <- admm_MADMMplasso_cpp(
       beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, W_hat, XtY, y,
@@ -88,7 +88,7 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
 
   I <- matrix(0, nrow = nrow(C) * dim(y)[2], ncol = dim(y)[2])
   II <- input[multiple_of_D]
-  diag(I[c(1:dim(y)[2]), ]) <- C[1, ] * (CW[1])
+  diag(I[1:dim(y)[2], ]) <- C[1, ] * (CW[1])
 
   c_count <- 2
   for (e in II[-length(II)]) {
@@ -122,7 +122,7 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
   EE_old <- EE
   res_pri <- 0
   res_dual <- 0
-  obj <- c()
+  obj <- NULL
 
   SVD_D <- Diagonal(x = svd.w$d)
   R_svd <- (svd.w$u %*% SVD_D) / N
@@ -139,8 +139,6 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
 
     XtY <- crossprod((W_hat), (new_y))
 
-    main_beta <- array(0, c(p, K + 1, D))
-
     res_val <- rho * (t(I) %*% (E) - (t(I) %*% (H)))
 
     v.diff1 <- matrix(0, D)
@@ -148,12 +146,12 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
     ee.diff1 <- matrix(0, D)
 
     new_G <- matrix(0, (p + p * K))
-    new_G[c(1:p)] <- 1
-    new_G[-c(1:p)] <- 2
-    new_G[c(1:p)] <- rho * (1 + new_G[c(1:p)])
-    new_G[-c(1:p)] <- rho * (1 + new_G[-c(1:p)])
+    new_G[1:p] <- 1
+    new_G[-1:-p] <- 2
+    new_G[1:p] <- rho * (1 + new_G[1:p])
+    new_G[-1:-p] <- rho * (1 + new_G[-1:-p])
 
-    invmat <- lapply(seq_len(D), function(j) {(new_G + rho * (new_I[j] + 1))})
+    invmat <- lapply(seq_len(D), function(j) return(new_G + rho * (new_I[j] + 1)))
 
     for (jj in 1:D) {
       group <- (rho) * (t(G) %*% t(V[, , jj]) - t(G) %*% t(O[, , jj]))
@@ -189,7 +187,7 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
       Q[, -1, jj] <- sign(new.mat) * pmax(abs(new.mat) - ((alpha * lambda[jj]) / (rho)), 0)
       b_hat <- alph * beta_hat1 + (1 - alph) * EE[, , jj]
       new.mat <- b_hat + HH[, , jj]
-      row.norm1 <- sqrt(apply(new.mat^2, 1, sum, na.rm = T))
+      row.norm1 <- sqrt(apply(new.mat^2, 1, sum, na.rm = TRUE))
       coef.term1 <- pmax(1 - (gg[2]) / rho / (row.norm1), 0)
       ee1 <- scale(t(as.matrix(new.mat)), center = FALSE, scale = 1 / coef.term1)
       EE[, , jj] <- t(ee1)
@@ -200,10 +198,10 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
 
       # Now we have the main part.
       new.mat <- Big_beta1 + O[, , jj]
-      new.mat1 <- new.mat[, c(1:(K + 1))]
-      new.mat2 <- new.mat[, -c(1:(K + 1))]
-      row.norm1 <- sqrt(apply(new.mat1^2, 1, sum, na.rm = T))
-      row.norm2 <- sqrt(apply(new.mat2^2, 1, sum, na.rm = T))
+      new.mat1 <- new.mat[, 1:(K + 1)]
+      new.mat2 <- new.mat[, -1:-(K + 1)]
+      row.norm1 <- sqrt(apply(new.mat1^2, 1, sum, na.rm = TRUE))
+      row.norm2 <- sqrt(apply(new.mat2^2, 1, sum, na.rm = TRUE))
 
       coef.term1 <- pmax(1 - ((1 - alpha) * lambda[jj]) / (rho) / (row.norm1), 0)
       coef.term2 <- pmax(1 - ((1 - alpha) * lambda[jj]) / (rho) / (row.norm2), 0)
@@ -230,11 +228,11 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
     beta.group <- (array(NA, c(p + p * K, dim(y)[2], dim(C)[1])))
     N_E <- list()
     II <- input[multiple_of_D]
-    new.mat_group[, , 1] <- t((new.mat[c(1:dim(y)[2]), ]))
-    beta.group[, , 1] <- t((Big_beta_respone[c(1:dim(y)[2]), ]))
+    new.mat_group[, , 1] <- t((new.mat[1:dim(y)[2], ]))
+    beta.group[, , 1] <- t((Big_beta_respone[1:dim(y)[2], ]))
 
     beta_transform <- matrix(0, p, (K + 1) * dim(y)[2])
-    beta_transform[, c(1:(1 + K))] <- matrix(new.mat_group[, 1, 1], ncol = (K + 1), nrow = p)
+    beta_transform[, 1:(1 + K)] <- matrix(new.mat_group[, 1, 1], ncol = (K + 1), nrow = p)
     input2 <- 1:(dim(y)[2] * (1 + K))
     multiple_of_K <- (input2 %% (K + 1)) == 0
     II2 <- input2[multiple_of_K]
@@ -245,14 +243,14 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
       e2 <- II2[c_count2]
     }
 
-    norm_res <- ((apply(beta_transform, c(1), twonorm)))
+    norm_res <- ((apply(beta_transform, 1, twonorm)))
     coef.term1 <- pmax(1 - (gg[1]) / rho / (norm_res), 0)
 
     N_E1 <- scale(t(beta_transform), center = FALSE, scale = 1 / coef.term1)
 
     N_E1 <- t(N_E1)
     beta_transform1 <- matrix(0, p + p * K, dim(y)[2])
-    beta_transform1[, 1] <- as.vector(N_E1[, c(1:(K + 1))])
+    beta_transform1[, 1] <- as.vector(N_E1[, 1:(K + 1)])
 
     input3 <- 1:(dim(y)[2] * (1 + K))
     multiple_of_K <- (input3 %% (K + 1)) == 0
@@ -272,7 +270,7 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
       beta.group[, , c_count] <- t(Big_beta_respone[c((e + 1):(c_count * dim(y)[2])), ])
 
       beta_transform <- matrix(0, p, (K + 1) * dim(y)[2])
-      beta_transform[, c(1:(1 + K))] <- matrix(new.mat_group[, 1, c_count], ncol = (K + 1), nrow = p)
+      beta_transform[, 1:(1 + K)] <- matrix(new.mat_group[, 1, c_count], ncol = (K + 1), nrow = p)
       input2 <- 1:(dim(y)[2] * (1 + K))
       multiple_of_K <- (input2 %% (K + 1)) == 0
       II2 <- input2[multiple_of_K]
@@ -283,14 +281,14 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
         e2 <- II2[c_count2]
       }
 
-      norm_res <- ((apply(beta_transform, c(1), twonorm)))
+      norm_res <- ((apply(beta_transform, 1, twonorm)))
       coef.term1 <- pmax(1 - (gg[1]) / rho / (norm_res), 0)
 
       N_E1 <- scale(t(beta_transform), center = FALSE, scale = 1 / coef.term1)
 
       N_E1 <- t(N_E1)
       beta_transform1 <- matrix(0, p + p * K, dim(y)[2])
-      beta_transform1[, 1] <- as.vector(N_E1[, c(1:(K + 1))])
+      beta_transform1[, 1] <- as.vector(N_E1[, 1:(K + 1)])
 
       input3 <- 1:(dim(y)[2] * (1 + K))
       multiple_of_K <- (input3 %% (K + 1)) == 0
@@ -307,9 +305,7 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
       e <- II[c_count]
     }
 
-    N_beta.group <- apply(beta.group, 3, twonorm)
-
-    E[c(1:dim(C)[2]), ] <- N_E[[1]]
+    E[1:dim(C)[2], ] <- N_E[[1]]
 
     c_count <- 2
     e <- II[-length(II)][1]
@@ -354,7 +350,7 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
       rho <- rho / 2
     }
 
-    if (my_print == T) {
+    if (my_print) {
       print(c(res_dual, e.dual, res_pri, e.primal))
     }
     if (res_pri <= e.primal && res_dual <= e.dual) {
@@ -362,10 +358,10 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
 
       # Update convergence message
       message("Convergence reached after ", i, " iterations")
-      converge <- T
+      converge <- TRUE
       break
     }
-    converge <- F
+    converge <- FALSE
   } ### iteration
 
   res_val <- t(I) %*% (E)
@@ -379,8 +375,8 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
     new_group[, -1] <- group2
     new_g_theta <- as.vector(new_group)
 
-    finB1 <- as.vector(beta_hat[c(1:p), jj]) * (new_g_theta[c(1:p)] != 0) * (as.vector((Q[, 1, jj])) != 0)
-    finB2 <- as.vector(beta_hat[-c(1:p), jj]) * (new_g_theta[-c(1:p)] != 0) * (as.vector((Q[, -1, jj])) != 0)
+    finB1 <- as.vector(beta_hat[1:p, jj]) * (new_g_theta[1:p] != 0) * (as.vector((Q[, 1, jj])) != 0)
+    finB2 <- as.vector(beta_hat[-1:-p, jj]) * (new_g_theta[-1:-p] != 0) * (as.vector((Q[, -1, jj])) != 0)
 
     beta_hat1 <- matrix(c(finB1, finB2), ncol = (K + 1), nrow = p)
     beta[, jj] <- beta_hat1[, 1]
diff --git a/R/compute_pliable.R b/R/compute_pliable.R
index 0b557b0..177b6d2 100644
--- a/R/compute_pliable.R
+++ b/R/compute_pliable.R
@@ -3,8 +3,8 @@
 #' @param X N by p matrix of predictors
 #' @param Z N by K matrix of modifying variables. The elements of Z  may represent quantitative or categorical variables, or a mixture of the two.
 #' Categorical variables should be coded by 0-1 dummy variables: for a k-level variable, one can use either k or k-1  dummy variables.
-#' @param theta    theta coefficients for a single response ncol(X) by ncol(Z)  
-#' @return a vector of length N of the calculated interaction term for a single response 
+#' @param theta    theta coefficients for a single response ncol(X) by ncol(Z)
+#' @return a vector of length N of the calculated interaction term for a single response
 #' @export
 compute_pliable <- function(X, Z, theta) {
   p <- ncol(X)
diff --git a/R/cv_MADMMplasso.R b/R/cv_MADMMplasso.R
index 7295f4c..b87da9a 100644
--- a/R/cv_MADMMplasso.R
+++ b/R/cv_MADMMplasso.R
@@ -29,32 +29,26 @@
 #' @return  results containing the CV values
 #' @example inst/examples/cv_MADMMplasso_example.R
 #' @export
-cv_MADMMplasso <- function(fit, nfolds, X, Z, y, alpha = 0.5, lambda = fit$Lambdas, max_it = 50000, e.abs = 1E-3, e.rel = 1E-3, nlambda, rho = 5, my_print = F, alph = 1, foldid = NULL, parallel = T, pal = 0, gg = c(7, 0.5), TT, tol = 1E-4, cl = 2) {
+cv_MADMMplasso <- function(fit, nfolds, X, Z, y, alpha = 0.5, lambda = fit$Lambdas, max_it = 50000, e.abs = 1E-3, e.rel = 1E-3, nlambda, rho = 5, my_print = FALSE, alph = 1, foldid = NULL, parallel = TRUE, pal = 0, gg = c(7, 0.5), TT, tol = 1E-4, cl = 2) {
   BIG <- 10e9
   no <- nrow(X)
-  ni <- ncol(X)
-  nz <- ncol(Z)
   ggg <- vector("list", nfolds)
 
   yhat <- array(NA, c(no, dim(y)[2], length(lambda[, 1])))
 
-  my_nzero <- matrix(0, nfolds, length(lambda[, 1]))
-
   if (is.null(foldid)) {
     foldid <- sample(rep(1:nfolds, ceiling(no / nfolds)), no, replace = FALSE)
   }
 
   nfolds <- length(table(foldid))
 
-  status.in <- NULL
-
   for (ii in 1:nfolds) {
     print(c("fold,", ii))
     oo <- foldid == ii
 
-    ggg[[ii]] <- MADMMplasso(X = X[!oo, , drop = F], Z = Z[!oo, , drop = F], y = y[!oo, , drop = F], alpha = alpha, my_lambda = lambda, lambda_min = .01, max_it = max_it, e.abs = e.abs, e.rel = e.rel, nlambda = length(lambda[, 1]), rho = rho, tree = TT, my_print = my_print, alph = alph, parallel = parallel, pal = pal, gg = gg, tol = tol, cl = cl)
+    ggg[[ii]] <- MADMMplasso(X = X[!oo, , drop = FALSE], Z = Z[!oo, , drop = FALSE], y = y[!oo, , drop = FALSE], alpha = alpha, my_lambda = lambda, lambda_min = 0.01, max_it = max_it, e.abs = e.abs, e.rel = e.rel, nlambda = length(lambda[, 1]), rho = rho, tree = TT, my_print = my_print, alph = alph, parallel = parallel, pal = pal, gg = gg, tol = tol, cl = cl)
 
-    cv_p <- predict.MADMMplasso(ggg[[ii]], X = X[oo, , drop = F], Z = Z[oo, ], y = y[oo, ])
+    cv_p <- predict.MADMMplasso(ggg[[ii]], X = X[oo, , drop = FALSE], Z = Z[oo, ], y = y[oo, ])
     ggg[[ii]] <- 0
     yhat[oo, , seq(nlambda)] <- cv_p$y_hat[, , seq(nlambda)]
   }
@@ -64,7 +58,6 @@ cv_MADMMplasso <- function(fit, nfolds, X, Z, y, alpha = 0.5, lambda = fit$Lambd
 
   mat <- err
   outmat <- matrix(NA, nfolds, ncol(mat))
-  good <- matrix(0, nfolds, ncol(mat))
   mat[is.infinite(mat)] <- NA
   for (i in seq(nfolds)) {
     mati <- mat[foldid == i, , drop = FALSE]
@@ -73,17 +66,20 @@ cv_MADMMplasso <- function(fit, nfolds, X, Z, y, alpha = 0.5, lambda = fit$Lambd
 
   non_zero <- c(fit$path$nzero)
 
-  cvm <- (apply(err, 2, mean, na.rm = T)) / dim(y)[2]
-  nn <- apply(!is.na(err), 2, sum, na.rm = T)
-  cvsd <- sqrt(apply(err, 2, var, na.rm = T) / (dim(y)[2] * nn))
+  cvm <- (apply(err, 2, mean, na.rm = TRUE)) / dim(y)[2]
+  nn <- apply(!is.na(err), 2, sum, na.rm = TRUE)
+  cvsd <- sqrt(apply(err, 2, var, na.rm = TRUE) / (dim(y)[2] * nn))
 
   cvm.nz <- cvm
   cvm.nz[non_zero == 0] <- BIG
   imin <- which.min(cvm.nz)
   imin.1se <- which(cvm < cvm[imin] + cvsd[imin])[1]
 
-  out <- list(lambda = fit$Lambdas, cvm = cvm, cvsd = cvsd, cvup = cvm +
-    cvsd, cvlo = cvm - cvsd, nz = c(fit$path$nzero), lambda.min = fit$Lambdas[imin, 1], lambda.1se = fit$Lambdas[imin.1se, 1])
+  out <- list(
+    lambda = fit$Lambdas, cvm = cvm, cvsd = cvsd, cvup = cvm + cvsd,
+    cvlo = cvm - cvsd, nz = c(fit$path$nzero),
+    lambda.min = fit$Lambdas[imin, 1], lambda.1se = fit$Lambdas[imin.1se, 1]
+  )
   class(out) <- "cv_MADMMplasso"
 
   return(out)
diff --git a/R/generate_my_w.R b/R/generate_my_w.R
index 701fe2e..10afad3 100644
--- a/R/generate_my_w.R
+++ b/R/generate_my_w.R
@@ -11,11 +11,9 @@
 #' @export
 generate_my_w <- function(X = matrix(), Z = matrix()) {
   p1 <- ncol(X)
-  p2 <- ncol(Z)
 
   # Just in case we have only one oberservation? Not sure why I did this
   if (is.vector(X)) p1 <- length(X)
-  if (is.vector(Z)) p2 <- length(Z)
 
   # Add the intercept
   x <- X
diff --git a/R/plot_coeff.R b/R/plot_coeff.R
index b24c6f8..25a60c3 100644
--- a/R/plot_coeff.R
+++ b/R/plot_coeff.R
@@ -15,14 +15,14 @@ plot_coeff <- function(beta, theta, error, nz, p, K, D, nlambda, Lambda) {
     n <- dim(my_beta)[2]
     matplot(b, t(my_beta), type = "n", col = "red", ylim = range(my_beta), xlab = "Log Lambda", ylab = (paste("coefficient", ii)))
     axis(
-      side = 3, at = (as.matrix(b)), labels = paste(as.matrix(nz[c(1:gg)])),
+      side = 3, at = (as.matrix(b)), labels = paste(as.matrix(nz[1:gg])),
       tick = FALSE, line = 0
     )
 
     for (i in 1:(p)) {
       lines(b, (my_beta[i, ]), col = i + 1, lty = 1)
 
-      text((min(b - .1)), my_beta[i, n], labels = i, cex = .7)
+      text((min(b - 0.1)), my_beta[i, n], labels = i, cex = 0.7)
     }
 
     my_beta <- (my_beta)
@@ -37,7 +37,7 @@ plot_coeff <- function(beta, theta, error, nz, p, K, D, nlambda, Lambda) {
     sbeta <- (my_beta)
     for (j in act) {
       for (i in 1:length(index)) {
-        if (ntheta[j, i] > 0) text(index[i], sbeta[j, i], label = "x", cex = .7)
+        if (ntheta[j, i] > 0) text(index[i], sbeta[j, i], label = "x", cex = 0.7)
       }
     }
   }
diff --git a/R/predict.MADMMplasso.R b/R/predict.MADMMplasso.R
index 64d7ac8..858ba11 100644
--- a/R/predict.MADMMplasso.R
+++ b/R/predict.MADMMplasso.R
@@ -33,7 +33,6 @@ predict.MADMMplasso <- function(object, X, Z, y, lambda = NULL, ...) {
 
   yh <- array(0, c(N, D, length(isel)))
   DEV <- matrix(NA, length(isel))
-  my_theta <- array(0, c(ncol(X), ncol(Z), length(isel)))
   pBETA0 <- lapply(
     seq_len(length(isel)),
     function(j) (matrix(0, nrow = (D)))
@@ -57,11 +56,6 @@ predict.MADMMplasso <- function(object, X, Z, y, lambda = NULL, ...) {
     function(j) (array(0, c(p, K, (D))))
   )
 
-  pY_HAT <- lapply(
-    seq_len(length(isel)),
-    function(j) (matrix(0, nrow = N, ncol = (D)))
-  )
-
   ii <- 0
   for (m in isel) {
     ii <- ii + 1
@@ -70,10 +64,6 @@ predict.MADMMplasso <- function(object, X, Z, y, lambda = NULL, ...) {
       seq_len(max(1)),
       function(j) (matrix(0, nrow = N))
     )
-    pr <- lapply(
-      seq_len(max(1)),
-      function(j) (matrix(0, nrow = N))
-    )
     beta0 <- object$beta0[[z]]
     beta <- as.matrix(object$beta[[z]])
     beta_hat <- as.matrix(object$BETA_hat[[z]])
diff --git a/R/quick_func.R b/R/quick_func.R
index 71d01a5..baa1e81 100644
--- a/R/quick_func.R
+++ b/R/quick_func.R
@@ -1,3 +1,3 @@
-quick_func <- function(xz = c(), xn) {
+quick_func <- function(xz = NULL, xn) {
   as.vector(xz[1:xn] %o% xz[-(1:xn)])
 }
diff --git a/R/sim2.R b/R/sim2.R
index dd44193..0eae753 100644
--- a/R/sim2.R
+++ b/R/sim2.R
@@ -14,7 +14,7 @@
 #'  Z: a N by K matrix of modifiers
 
 #' @export
-sim2 <- function(p = 500, n = 100, m = 24, nz = 4, rho = .4, B.elem = 0.5) {
+sim2 <- function(p = 500, n = 100, m = 24, nz = 4, rho = 0.4, B.elem = 0.5) {
   b <- 10
   if (!is.na(p)) {
     # generate covariance matrix
diff --git a/R/tree_parms.R b/R/tree_parms.R
index 5ac764b..d09e571 100644
--- a/R/tree_parms.R
+++ b/R/tree_parms.R
@@ -9,7 +9,7 @@
 #'  y.colnames: names of the response
 
 #' @export
-tree_parms <- function(y = y, h = .7) {
+tree_parms <- function(y = y, h = 0.7) {
   m <- dim(y)[2]
   myDist0 <- 1 - abs(fast_corr(y))
   myDist <- myDist0[lower.tri(myDist0)]
diff --git a/inst/examples/MADMMplasso_example.R b/inst/examples/MADMMplasso_example.R
index 46a2c2f..15c6a70 100644
--- a/inst/examples/MADMMplasso_example.R
+++ b/inst/examples/MADMMplasso_example.R
@@ -29,7 +29,7 @@ beta_5[11:15] <- c(-2, -2, -2, -2, -2)
 beta_6[11:15] <- c(-2, -2, -2, -2, -2)
 
 Beta <- cbind(beta_1, beta_2, beta_3, beta_4, beta_5, beta_6)
-colnames(Beta) <- c(1:6)
+colnames(Beta) <- 1:6
 
 theta <- array(0, c(p, K, 6))
 theta[1, 1, 1] <- 2
@@ -57,7 +57,6 @@ theta[13, 2, 6] <- 2
 theta[14, 3, 6] <- -2
 theta[15, 4, 6] <- -2
 
-library(MASS)
 pliable <- matrix(0, N, 6)
 for (e in 1:6) {
   pliable[, e] <- compute_pliable(X, Z, theta[, , e])
@@ -68,7 +67,7 @@ e <- MASS::mvrnorm(N, mu = rep(0, 6), Sigma = esd)
 y_train <- X %*% Beta + pliable + e
 y <- y_train
 
-colnames(y) <- c(paste("y", 1:(ncol(y)), sep = ""))
+colnames(y) <- c(paste0("y", 1:(ncol(y))))
 TT <- tree_parms(y)
 plot(TT$h_clust)
 gg1 <- matrix(0, 2, 2)
@@ -78,7 +77,7 @@ gg1[2, ] <- c(0.02, 0.02)
 nlambda <- 1
 e.abs <- 1E-4
 e.rel <- 1E-2
-alpha <- .2
+alpha <- 0.2
 tol <- 1E-3
 fit <- MADMMplasso(
   X, Z, y,
diff --git a/inst/examples/cv_MADMMplasso_example.R b/inst/examples/cv_MADMMplasso_example.R
index 821b3c2..b18bf91 100644
--- a/inst/examples/cv_MADMMplasso_example.R
+++ b/inst/examples/cv_MADMMplasso_example.R
@@ -1,3 +1,4 @@
+ # nolint start: indentation_linter
 \dontrun{
   # Train the model
   # generate some data
@@ -93,3 +94,4 @@
   )
   plot(cv_admp)
 }
+# nolint end: indentation_linter
diff --git a/man/MADMMplasso.Rd b/man/MADMMplasso.Rd
index 1e49518..adc7e13 100644
--- a/man/MADMMplasso.Rd
+++ b/man/MADMMplasso.Rd
@@ -18,10 +18,10 @@ MADMMplasso(
   maxgrid,
   nlambda,
   rho = 5,
-  my_print = F,
+  my_print = FALSE,
   alph = 1.8,
   tree,
-  parallel = T,
+  parallel = TRUE,
   pal = 0,
   gg = NULL,
   tol = 1e-04,
@@ -133,7 +133,7 @@ beta_5[11:15] <- c(-2, -2, -2, -2, -2)
 beta_6[11:15] <- c(-2, -2, -2, -2, -2)
 
 Beta <- cbind(beta_1, beta_2, beta_3, beta_4, beta_5, beta_6)
-colnames(Beta) <- c(1:6)
+colnames(Beta) <- 1:6
 
 theta <- array(0, c(p, K, 6))
 theta[1, 1, 1] <- 2
@@ -161,7 +161,6 @@ theta[13, 2, 6] <- 2
 theta[14, 3, 6] <- -2
 theta[15, 4, 6] <- -2
 
-library(MASS)
 pliable <- matrix(0, N, 6)
 for (e in 1:6) {
   pliable[, e] <- compute_pliable(X, Z, theta[, , e])
@@ -172,7 +171,7 @@ e <- MASS::mvrnorm(N, mu = rep(0, 6), Sigma = esd)
 y_train <- X \%*\% Beta + pliable + e
 y <- y_train
 
-colnames(y) <- c(paste("y", 1:(ncol(y)), sep = ""))
+colnames(y) <- c(paste0("y", 1:(ncol(y))))
 TT <- tree_parms(y)
 plot(TT$h_clust)
 gg1 <- matrix(0, 2, 2)
@@ -182,7 +181,7 @@ gg1[2, ] <- c(0.02, 0.02)
 nlambda <- 1
 e.abs <- 1E-4
 e.rel <- 1E-2
-alpha <- .2
+alpha <- 0.2
 tol <- 1E-3
 fit <- MADMMplasso(
   X, Z, y,
diff --git a/man/admm_MADMMplasso.Rd b/man/admm_MADMMplasso.Rd
index 68f0b77..18996a8 100644
--- a/man/admm_MADMMplasso.Rd
+++ b/man/admm_MADMMplasso.Rd
@@ -25,7 +25,7 @@ admm_MADMMplasso(
   alph,
   svd.w,
   tree,
-  my_print = T,
+  my_print = TRUE,
   invmat,
   gg = 0.2,
   legacy = FALSE
diff --git a/man/cv_MADMMplasso.Rd b/man/cv_MADMMplasso.Rd
index 295756f..0c92fff 100644
--- a/man/cv_MADMMplasso.Rd
+++ b/man/cv_MADMMplasso.Rd
@@ -18,10 +18,10 @@ cv_MADMMplasso(
   e.rel = 0.001,
   nlambda,
   rho = 5,
-  my_print = F,
+  my_print = FALSE,
   alph = 1,
   foldid = NULL,
-  parallel = T,
+  parallel = TRUE,
   pal = 0,
   gg = c(7, 0.5),
   TT,
@@ -85,6 +85,7 @@ Carries out cross-validation for  a  MADMMplasso model over a path of regulariza
 @description Carries out cross-validation for  a  MADMMplasso model over a path of regularization values
 }
 \examples{
+ # nolint start: indentation_linter
 \dontrun{
   # Train the model
   # generate some data
@@ -180,4 +181,5 @@ Carries out cross-validation for  a  MADMMplasso model over a path of regulariza
   )
   plot(cv_admp)
 }
+# nolint end: indentation_linter
 }
diff --git a/tests/testthat/test-admm_MADMMplasso_cpp.R b/tests/testthat/test-admm_MADMMplasso_cpp.R
index b2fe762..b686c4d 100644
--- a/tests/testthat/test-admm_MADMMplasso_cpp.R
+++ b/tests/testthat/test-admm_MADMMplasso_cpp.R
@@ -1,13 +1,9 @@
 # Auxiliary funcions ===========================================================
 # TODO: ask why these are not in the package?
 model <- function(beta0, theta0, beta, theta, X, Z) {
-  p <- ncol(X)
   N <- nrow(X)
-  K <- ncol(Z)
-  D <- dim(beta0)[2]
   intercepts <- (matrix(1, N)) %*% beta0 + Z %*% ((theta0))
   shared_model <- X %*% (beta)
-  pliable <- matrix(0, N, D)
   return(intercepts + shared_model)
 }
 
@@ -22,7 +18,7 @@ reg_temp <- function(r, Z) {
     my_inv <- pracma::pinv(t(my_w) %*% my_w)
     my_res <- my_inv %*% (t(my_w) %*% r[, e])
     beta01[e] <- matrix(my_res[(K + 1)])
-    theta01[, e] <- matrix(my_res[c(1:(K))])
+    theta01[, e] <- matrix(my_res[1:K])
   }
   return(list(beta0 = beta01, theta0 = theta01))
 }
@@ -57,7 +53,7 @@ beta_4[6:10] <- c(2, 2, 2, -2, -2)
 beta_5[11:15] <- c(-2, -2, -2, -2, -2)
 beta_6[11:15] <- c(-2, -2, -2, -2, -2)
 Beta <- cbind(beta_1, beta_2, beta_3, beta_4, beta_5, beta_6)
-colnames(Beta) <- c(1:6)
+colnames(Beta) <- 1:6
 theta <- array(0, c(p, K, 6))
 theta[1, 1, 1] <- 2
 theta[3, 2, 1] <- 2
@@ -91,8 +87,8 @@ esd <- diag(6)
 e <- mvrnorm(N, mu = rep(0, 6), Sigma = esd)
 y_train <- X %*% Beta + pliable + e
 y <- y_train
-colnames(y) <- c(1:6)
-colnames(y) <- c(paste("y", 1:(ncol(y)), sep = ""))
+colnames(y) <- 1:6
+colnames(y) <- c(paste0("y", 1:(ncol(y))))
 TT <- tree_parms(y)
 C <- TT$Tree
 CW <- TT$Tw
@@ -104,7 +100,7 @@ lambda <- rep(0.5, 6)
 alpha <- 0.5
 e.abs <- 1E-4
 e.rel <- 1E-2
-alpha <- .2
+alpha <- 0.2
 tol <- 1E-3
 alph <- 1
 rho <- 5
@@ -118,7 +114,7 @@ input <- 1:(dim(y)[2] * nrow(C))
 multiple_of_D <- (input %% dim(y)[2]) == 0
 I <- matrix(0, nrow = nrow(C) * dim(y)[2], ncol = dim(y)[2])
 II <- input[multiple_of_D]
-diag(I[c(1:dim(y)[2]), ]) <- C[1, ] * (CW[1])
+diag(I[1:dim(y)[2], ]) <- C[1, ] * (CW[1])
 c_count <- 2
 for (e in II[-length(II)]) {
   diag(I[c((e + 1):(c_count * dim(y)[2])), ]) <- C[c_count, ] * (CW[c_count])
@@ -126,15 +122,15 @@ for (e in II[-length(II)]) {
 }
 new_I <- diag(t(I) %*% I)
 new_G <- matrix(0, (p + p * K))
-new_G[c(1:p)] <- 1
-new_G[-c(1:p)] <- 2
-new_G[c(1:p)] <- rho * (1 + new_G[c(1:p)])
-new_G[-c(1:p)] <- rho * (1 + new_G[-c(1:p)])
+new_G[1:p] <- 1
+new_G[-1:-p] <- 2
+new_G[1:p] <- rho * (1 + new_G[1:p])
+new_G[-1:-p] <- rho * (1 + new_G[-1:-p])
 invmat <- list() # denominator of the beta estimates
 for (rr in 1:D) {
   DD1 <- rho * (new_I[rr] + 1)
   DD2 <- new_G + DD1
-  invmat[[rr]] <- DD2 # Matrix::chol2inv( Matrix::chol(new_sparse) )
+  invmat[[rr]] <- DD2
 }
 beta0 <- matrix(0, 1, D) # estimates$Beta0
 theta0 <- matrix(0, K, D)
@@ -151,7 +147,7 @@ Q <- (array(0, c(p, (1 + K), D)))
 P <- (array(0, c(p, (1 + K), D)))
 H <- (matrix(0, dim(y)[2] * nrow(C), (p + p * K))) # response multiplier
 HH <- (array(0, c(p, (1 + K), D)))
-r_current <- y #-model(beta0,theta0,beta=beta_hat, theta, X=W_hat, Z)
+r_current <- y
 b <- reg_temp(r_current, Z) # Analytic solution how no sample lower bound (Z.T @ Z + cI)^-1 @ (Z.T @ r)
 beta0 <- b$beta0
 theta0 <- b$theta0
@@ -175,13 +171,13 @@ beta_hat <- my_values$beta_hat
 y_hat <- my_values$y_hat
 
 test_that("final objects have correct dimensions", {
-  expect_equal(dim(beta0), c(1, 6))
-  expect_equal(dim(theta0), c(4, 6))
-  expect_equal(dim(beta), c(50, 6))
-  expect_equal(dim(theta), c(50, 4, 6))
-  expect_equal(length(converge), 1)
-  expect_equal(dim(beta_hat), c(250, 6))
-  expect_equal(dim(y_hat), c(100, 6))
+  expect_identical(dim(beta0), c(1L, 6L))
+  expect_identical(dim(theta0), c(4L, 6L))
+  expect_identical(dim(beta), c(50L, 6L))
+  expect_identical(dim(theta), c(50L, 4L, 6L))
+  expect_length(converge, 1L)
+  expect_identical(dim(beta_hat), c(250L, 6L))
+  expect_identical(dim(y_hat), c(100L, 6L))
 })
 
 test_that("mean values of final objects are expected", {
@@ -190,7 +186,7 @@ test_that("mean values of final objects are expected", {
   expect_equal(mean(theta0), 5.123034e-02, tolerance = tole)
   expect_equal(mean(beta), 2.104393e-02, tolerance = tole)
   expect_equal(mean(theta), 2.841666e-04, tolerance = tole)
-  expect_equal(converge, TRUE)
+  expect_true(converge)
   expect_equal(mean(beta_hat), 4.436118e-03, tolerance = tole)
   expect_equal(mean(y_hat), -8.380419e-02, tolerance = tole)
 })
@@ -204,8 +200,8 @@ my_values_cpp <- admm_MADMMplasso(
 )
 
 test_that("C++ function output structure", {
-  expect_equal(length(my_values_cpp), length(my_values))
-  expect_equal(names(my_values_cpp), names(my_values))
+  expect_identical(length(my_values_cpp), length(my_values))
+  expect_identical(names(my_values_cpp), names(my_values))
 })
 
 test_that("Values are the same", {
@@ -214,7 +210,7 @@ test_that("Values are the same", {
   expect_equal(my_values$theta0, my_values_cpp$theta0, tolerance = tl)
   expect_equal(my_values$beta, my_values_cpp$beta, tolerance = tl)
   expect_equal(my_values$theta, my_values_cpp$theta, tolerance = tl)
-  expect_equal(my_values$converge, my_values_cpp$converge)
+  expect_identical(my_values$converge, my_values_cpp$converge)
   expect_equal(my_values$beta_hat, my_values_cpp$beta_hat, tolerance = tl)
   expect_equal(my_values$y_hat, my_values_cpp$y_hat, tolerance = tl)
 })
diff --git a/tests/testthat/test-aux_cpp_functions.R b/tests/testthat/test-aux_cpp_functions.R
index 54ba17c..b80e5ca 100644
--- a/tests/testthat/test-aux_cpp_functions.R
+++ b/tests/testthat/test-aux_cpp_functions.R
@@ -3,24 +3,24 @@ set.seed(1282897)
 y <- rpois(n = 10, lambda = 10)
 test_that("multiples_of() works", {
   expect_error(multiples_of(x, 0), "divisor cannot be 0")
-  expect_equal(multiples_of(x, 1), matrix(1:12))
-  expect_equal(multiples_of(x, 2), matrix(c(2, 4, 6, 8, 10, 12)))
-  expect_equal(multiples_of(x, 3), matrix(c(3, 6, 9, 12)))
-  expect_equal(multiples_of(x, 4), matrix(c(4, 8, 12)))
-  expect_equal(multiples_of(x, 5), matrix(c(5, 10)))
-  expect_equal(multiples_of(x, 6), matrix(c(6, 12)))
+  expect_identical(multiples_of(x, 1), matrix(1:12))
+  expect_identical(multiples_of(x, 2), matrix(c(2L, 4L, 6L, 8L, 10L, 12L)))
+  expect_identical(multiples_of(x, 3), matrix(c(3L, 6L, 9L, 12L)))
+  expect_identical(multiples_of(x, 4), matrix(c(4L, 8L, 12L)))
+  expect_identical(multiples_of(x, 5), matrix(c(5L, 10L)))
+  expect_identical(multiples_of(x, 6), matrix(c(6L, 12L)))
 
-  expect_equal(multiples_of(y, 4), matrix(c(2, 3, 4, 6)))
-  expect_equal(multiples_of(y, 5), matrix(c(1, 5, 9)))
-  expect_equal(multiples_of(y, 7), matrix(c(8, 10)))
-  expect_equal(multiples_of(y, 15), matrix(c(1, 5, 9)))
-  expect_equal(multiples_of(y, 10), matrix(as.integer(NULL)))
+  expect_identical(multiples_of(y, 4), matrix(c(2L, 3L, 4L, 6L)))
+  expect_identical(multiples_of(y, 5), matrix(c(1L, 5L, 9L)))
+  expect_identical(multiples_of(y, 7), matrix(c(8L, 10L)))
+  expect_identical(multiples_of(y, 15), matrix(c(1L, 5L, 9L)))
+  expect_identical(multiples_of(y, 10), matrix(as.integer(NULL)))
 
-  expect_equal(multiples_of(y, 4, TRUE), matrix(c(y[2], y[3], y[4], y[6])))
-  expect_equal(multiples_of(y, 5, TRUE), matrix(c(y[1], y[5], y[9])))
-  expect_equal(multiples_of(y, 7, TRUE), matrix(c(y[8], y[10])))
-  expect_equal(multiples_of(y, 15, TRUE), matrix(c(y[1], y[5], y[9])))
-  expect_equal(multiples_of(y, 10, TRUE), matrix(as.integer(NULL)))
+  expect_identical(multiples_of(y, 4, TRUE), matrix(c(y[2], y[3], y[4], y[6])))
+  expect_identical(multiples_of(y, 5, TRUE), matrix(c(y[1], y[5], y[9])))
+  expect_identical(multiples_of(y, 7, TRUE), matrix(c(y[8], y[10])))
+  expect_identical(multiples_of(y, 15, TRUE), matrix(c(y[1], y[5], y[9])))
+  expect_identical(multiples_of(y, 10, TRUE), matrix(as.integer(NULL)))
 })
 
 N <- rpois(1, 10)
@@ -44,6 +44,6 @@ test_that("modulo() works", {
   top <- rpois(1, 20)
   for (i in seq_len(top)) {
     x <- matrix(rpois(top, i))
-    expect_equal(x %% i, modulo(x, i))
+    expect_equal(x %% i, modulo(x, i), ignore_attr = TRUE)
   }
 })