diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/check_inputs.R | 64 | ||||
| -rw-r--r-- | R/cv_loss.R | 17 | ||||
| -rw-r--r-- | R/plot_print_methods.R | 26 | ||||
| -rw-r--r-- | R/rmad.R | 31 | ||||
| -rw-r--r-- | R/rsc.R | 42 | ||||
| -rw-r--r-- | R/rsc_cv.R | 183 | ||||
| -rwxr-xr-x | R/zzz.R | 4 |
7 files changed, 367 insertions, 0 deletions
diff --git a/R/check_inputs.R b/R/check_inputs.R new file mode 100644 index 0000000..4cbcf0a --- /dev/null +++ b/R/check_inputs.R @@ -0,0 +1,64 @@ +.check_input_data_matrix <- function(x, y, na.rm) { + if (is.null(y)) { + if (!is.array(x) & !is.data.frame(x) & !is.matrix(x)) { + stop("\"x\" must be a numeric matrix or any other array-type object that can be converted or a \"matrix\" object with with ncol>=2.\n\n") + } + if (is.vector(x)) { + stop("\"x\" must be a numeric matrix or any other array-type object that can be converted or a \"matrix\" object. with with ncol>=2.\n\n") + } + if (!is.matrix(x)) { + x <- data.matrix(x) + } + if (!is.numeric(x)) { + stop("\"x\" must be numeric.") + } + if (nrow(x) < 2 | ncol(x) < 2) { + stop("nrow(x)>=2 and ncol(xa)>=2 are required\n\n") + } + is_na_data <- is.na(x) + if (any(is_na_data)) { + if (na.rm == FALSE) { + stop("\"x\" contains NA records. You may want to filter NAs by setting \"na.rm=TRUE\" (see documentation for more details).\n\n") + } + else { + idx_na <- which(rowSums(is_na_data) >= 1) + x <- x[-idx_na, , drop = FALSE] + if (nrow(x) < 2) { + stop("nrow(x)<2 after NA removal.\n\n") + } + } + } + if (any(!is.finite(x))) { + stop("\"x\" contains Inf values\n\n") + } + } + else { + if (!is.vector(x) | !is.vector(y)) { + stop("If \"y\" is given, \"x\" and \"y\" must be both numeric.\n\n") + } + if (!is.numeric(x) | !is.numeric(y)) { + stop("\"x\" and \"y\" must be numeric.\n\n") + } + if (length(x) != length(y)) { + stop("\"x\" and \"y\" have different length.\n\n") + } + x <- cbind(x, y, deparse.level = 0) + is_na_data <- is.na(x) + if (any(is_na_data)) { + if (na.rm == FALSE) { + stop("\"x\" or \"y\" contains NA records. You may want to filter NAs by setting \"na.rm=TRUE\" (see documentation for more details).\n\n") + } + else { + idx_na <- which(rowSums(is_na_data) >= 1) + x <- x[-idx_na, , drop = FALSE] + if (nrow(x) < 2) { + stop("length(x)<2 and/or length(y)<2 after NA removal.\n\n") + } + } + } + if (any(!is.finite(x))) { + stop("\"x\" and/or \"y\" contains Inf values\n\n") + } + } + return(x) +} diff --git a/R/cv_loss.R b/R/cv_loss.R new file mode 100644 index 0000000..ea46007 --- /dev/null +++ b/R/cv_loss.R @@ -0,0 +1,17 @@ +.cv_loss <- function(idx, dat, evencorrection, threshold, grid.length, p, nc) { + res <- numeric(nc) + n1 <- as.integer(sum(idx)) + C1 <- .Fortran("cormadvecdp", matrix = dat[idx, ], nrow = n1, ncol = p, res = res, + ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res + n2 <- as.integer(sum(!idx)) + C2 <- .Fortran("cormadvecdp", matrix = dat[!idx, ], nrow = n2, ncol = p, res = res, + ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res + ans <- rep(0, times = grid.length) + for (h in 1:grid.length) { + C1[abs(C1) < threshold[h]] <- 0 + ans[h] <- sum(2 * { + C1 - C2 + }^2)/p + } + return(ans) +} diff --git a/R/plot_print_methods.R b/R/plot_print_methods.R new file mode 100644 index 0000000..ca163cb --- /dev/null +++ b/R/plot_print_methods.R @@ -0,0 +1,26 @@ +print.rsc_cv <- function(x, ...) { + cat("\n") + cat("================================================\n") + cat(" Expected Normalized Squared Frobenius Loss \n") + cat("================================================\n") + print(x$loss) + cat("================================================\n") + cat("\n") +} +plot.rsc_cv <- function(x, ...) { + tstar <- which(x$loss$Flag == "minimum") + hstar <- x$loss$Threshold[tstar] + inf_loss <- x$loss$Average - x$loss$SE + sup_loss <- x$loss$Average + x$loss$SE + a <- inf_loss[tstar] + b <- sup_loss[tstar] + hstar1se <- max(x$loss$Threshold[which(x$loss$Flag == "*")]) + Ylim <- range(c(inf_loss, sup_loss)) + plot(x$loss$Threshold, x$loss$Average, t = "b", ylim = Ylim, pch = 20, col = "#0052A5", + lwd = 2, main = "RSC Optimal Threshold Selection", xlab = "Threshold", ylab = "Average loss", + ...) + arrows(x$loss$Threshold, inf_loss, x$loss$Threshold, sup_loss, length = 0.05, + angle = 90, code = 3, col = "#0052A5") + abline(v = hstar1se, col = "#31A853", lty = 2, lwd = 2) + abline(v = hstar, col = "#E0162B", lty = 2, lwd = 2) +} diff --git a/R/rmad.R b/R/rmad.R new file mode 100644 index 0000000..b8e9d41 --- /dev/null +++ b/R/rmad.R @@ -0,0 +1,31 @@ +rmad <- function(x, y = NULL, na.rm = FALSE, even.correction = FALSE) { + dat <- .check_input_data_matrix(x = x, y = y, na.rm = na.rm) + colnames_original <- colnames(dat) + storage.mode(dat) <- "double" + n <- as.integer(nrow(dat)) + p <- as.integer(ncol(dat)) + nc <- as.integer({ + p^2 - p + }/2) + if (even.correction) { + evencorrection <- 1L + } + else { + evencorrection <- 0L + } + u <- .Fortran("cormadvecdp", matrix = dat, nrow = n, ncol = p, res = numeric(nc), + ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res + if (!is.null(y)) { + return(u) + } + else { + R <- Matrix(1, nrow = p, ncol = p, sparse = FALSE) + R[lower.tri(R, diag = FALSE)] <- u + R <- forceSymmetric(R, uplo = "L") + R <- as(R, "dspMatrix") + if (!is.null(colnames_original)) { + dimnames(R)[[1]] <- dimnames(R)[[2]] <- colnames_original + } + return(R) + } +} @@ -0,0 +1,42 @@ +rsc <- function(cv, threshold = "minimum") { + if (class(cv) == "rsc_cv") { + if (is.numeric(threshold)) { + if (length(threshold) > 1) { + stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)") + } + else if (threshold <= 0 | threshold >= 1) { + stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)") + } + } + else { + if ({ + threshold != "minimum" + } & { + threshold != "minimum1se" + }) { + stop("'threshold' must be one of the following: 'minimum', 'minimum1se', a numeric value in (0,1).") + } + if (threshold == "minimum") { + threshold <- cv$minimum + } + else if (threshold == "minimum1se") { + threshold <- cv$minimum1se + } + } + cv$rmadvec[abs(cv$rmadvec) < threshold] <- 0 + nc <- length(cv$rmadvec) + p <- { + 1 + sqrt(1 + 8 * nc) + }/2 + R <- Matrix(1, nrow = p, ncol = p, sparse = TRUE) + R[lower.tri(R, diag = FALSE)] <- cv$rmadvec + R <- forceSymmetric(R, uplo = "L") + if (!is.null(cv$varnames)) { + dimnames(R)[[1]] <- dimnames(R)[[2]] <- cv$varnames + } + } + else { + stop("'cv' must be a an object of class 'rsc_cv' obtained from 'rsc::rsc_cv'") + } + return(R) +} diff --git a/R/rsc_cv.R b/R/rsc_cv.R new file mode 100644 index 0000000..4153322 --- /dev/null +++ b/R/rsc_cv.R @@ -0,0 +1,183 @@ +rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0.95, + by = 0.025), even.correction = FALSE, na.rm = FALSE, ncores = NULL, monitor = TRUE) { + dat <- .check_input_data_matrix(x = x, y = NULL, na.rm = na.rm) + colnames_original <- colnames(dat) + storage.mode(dat) <- "double" + n <- as.integer(nrow(dat)) + p <- as.integer(ncol(dat)) + nc <- as.integer({ + p^2 - p + }/2) + if ({ + cv.type != "random" + } & { + cv.type != "kfold" + }) { + stop("\"cv.type\" must be either \"random\" (default) or \"kfold\"") + } + if (!is.numeric(R)) { + stop("\"R\" must be an integer > 1") + } + else if (R < 1) { + stop("\"R\" must be an integer > 1") + } + if (!is.numeric(K)) { + stop("\"K\" must be an integer > 1") + } + else if (R < 1) { + stop("\"K\" must be an integer > 1") + } + if (length(threshold) == 1) { + if (threshold <= 0 | threshold >= 1) { + stop("\"threshold\" value does not belong to the interval (0,1).") + } + } + else if (length(threshold) > 1) { + if (any(threshold < 0) | any(threshold >= 1)) { + stop("Some of the \"threshold\" values do not belong to the interval (0,1).") + } + grid.length <- length(threshold) + } + if (even.correction) { + evencorrection <- 1L + } + else { + evencorrection <- 0L + } + if (is.null(ncores)) { + DetectedCores <- detectCores() + if (DetectedCores <= 2) { + ncores <- 1 + } + else { + ncores <- { + DetectedCores - 1 + } + } + } + else { + ncores <- as.integer(ncores) + if (ncores <= 0) { + stop("\"ncores\" must be an integer larger or equal to 1.") + } + } + if (monitor) { + cat("\n") + message("Computing the RMAD matrix") + t0 <- Sys.time() + } + rmad_vec <- .Fortran("cormadvecdp", matrix = dat, nrow = n, ncol = p, res = numeric(nc), + ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res + if (monitor) { + t1 <- Sys.time() + dt01 <- difftime(t1, t0, units = "auto") + message("* RMAD computing time:...... ", round(dt01, 2), " [", attributes(dt01)$units, + "]") + } + if (cv.type == "random") { + if (monitor) { + cat("\n") + message("Performing cross-validation") + t_hat <- round(1.2 * { + { + dt01 * R * 2 + }/ncores + }, 2) + message("* predicted end time (worst case):...... ", Sys.time() + t_hat) + } + n1 <- n - floor(n/log(n)) + IDX <- array(FALSE, dim = c(R, n)) + for (r in 1:R) { + IDX[r, ][sample(1:n, size = n1, replace = FALSE)] <- TRUE + } + registerDoParallel(ncores) + U <- foreach(r = 1:R) %dopar% { + .cv_loss(idx = IDX[r, ], dat = dat, evencorrection = evencorrection, + threshold = threshold, grid.length = grid.length, p = p, nc = nc) + } + stopImplicitCluster() + LOSS <- array(0, dim = c(R, grid.length)) + for (r in 1:R) { + LOSS[r, ] <- U[[r]] + } + avg_loss <- apply(LOSS, 2, mean) + se_loss <- apply(LOSS, 2, sd)/sqrt(R) + } + if (cv.type == "kfold") { + if (monitor) { + cat("\n") + message("Performing cross-validation") + t_hat <- round(1.2 * { + { + dt01 * R * K * 2 + }/ncores + }, 2) + message("* predicted end time (worst case):...... ", Sys.time() + t_hat) + } + idx_fold <- cut(1:n, breaks = K, labels = FALSE) + IDX <- array(TRUE, dim = c(R * K, n)) + row_count <- 1L + for (r in 1:R) { + idx_fold_shuffle <- sample(idx_fold, size = n, replace = FALSE) + for (k in 1:K) { + IDX[row_count, ][idx_fold_shuffle == k] <- FALSE + row_count <- 1L + row_count + } + } + registerDoParallel(ncores) + U <- foreach(r = 1:{ + R * K + }) %dopar% { + .cv_loss(idx = IDX[r, ], dat = dat, evencorrection = evencorrection, + threshold = threshold, grid.length = grid.length, p = p, nc = nc) + } + stopImplicitCluster() + if (R == 1) { + LOSS <- array(0, dim = c(K, grid.length)) + for (k in 1:K) { + LOSS[k, ] <- U[[k]] + } + } + else { + LOSS <- array(0, dim = c(K, grid.length, R)) + dimnames(LOSS)[[3]] <- paste0("r", 1:R) + dimnames(LOSS)[[1]] <- paste0("k", 1:K) + row_count <- 1L + for (r in 1:R) { + for (k in 1:K) { + LOSS[k, , r] <- U[[row_count]] + row_count <- 1L + row_count + } + } + } + avg_loss_r <- sd_loss_r <- matrix(0, nrow = R, ncol = grid.length) + for (r in 1:R) { + avg_loss_r[r, ] <- apply(LOSS[, , r], 2, mean) + sd_loss_r[r, ] <- apply(LOSS[, , r], 2, sd)/sqrt(K) + } + avg_loss <- apply(avg_loss_r, 2, mean) + se_loss <- apply(sd_loss_r, 2, mean) + if (monitor) { + t2 <- Sys.time() + dt02 <- difftime(t2, t0, units = "auto") + } + } + if (monitor) { + message("* finished on:.......................... ", Sys.time()) + } + tstar <- which.min(avg_loss) + flags <- rep("", grid.length) + a <- avg_loss[tstar] - se_loss[tstar] + b <- avg_loss[tstar] + se_loss[tstar] + flags[{ + avg_loss >= a + } & { + avg_loss <= b + }] <- "*" + flags[tstar] <- "minimum" + res <- data.frame(Threshold = threshold, Average = avg_loss, SE = se_loss, Flag = flags) + ans <- list(rmadvec = rmad_vec, varnames = colnames_original, loss = res, minimum = threshold[tstar], + minimum1se = max(threshold[avg_loss >= a & avg_loss <= b])) + class(ans) <- "rsc_cv" + return(ans) +} @@ -0,0 +1,4 @@ +.onAttach <- function(lib, pkg) { + packageStartupMessage("\nRSC: robust and sparse correlation matrix estimation\n Type 'citation(\"RSC\")' for citing this package\n") + invisible() +} |