diff options
Diffstat (limited to 'man/rsc_cv.Rd')
| -rwxr-xr-x | man/rsc_cv.Rd | 228 |
1 files changed, 228 insertions, 0 deletions
diff --git a/man/rsc_cv.Rd b/man/rsc_cv.Rd new file mode 100755 index 0000000..3cd2da6 --- /dev/null +++ b/man/rsc_cv.Rd @@ -0,0 +1,228 @@ +\name{rsc_cv} + +\alias{rsc_cv} + +\title{Optimal threshold selection for the RSC estimator} + +\description{ + Perform cross-validation to select an adaptive optimal threshold for + the RSC estimator proposed in Serra et al. (2018). +} + + +\usage{ + rsc_cv(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) +} + + +\arguments{ + \item{x}{ + A matrix or a data.frame. Rows of \code{x} correspond to sample units + and columns correspond to variables. Categorical variables are not + allowed. + } + \item{cv.type}{ + A character string indicating the cross-validation algorithm. Possible + values are \code{"kfold"} for repeated K-fold cross-validation, and + \code{"random"} for random cross-validation (see \emph{Details}). + } + \item{R}{ + An integer corresponding to the number of repeated foldings when + \code{cv.type = "kfold"}. When \code{cv.type = "random"} \code{R} + defines the number of random splits (see \emph{Details}). + } + \item{K}{ + An integer corresponding to the number of \emph{folds} in K-fold + cross-validation. Therefore this argument is not relevant when + \code{cv.type = "random"}. + } + \item{threshold}{ + A sequence of reals taken onto the interval (0,1) defining the + threshold values at which the loss is estimated. + } + \item{even.correction}{ + A logical value. It sets the parameter \code{even.correction} in + each of the underlying RMAD computations (see \emph{Details} in + \code{\link{rmad}}). + } + \item{na.rm}{ + A logical value, it defines the treatment of missing values in + each of the underlying RMAD computations (see \emph{Details}). + } + \item{ncores}{ + An integer value defining the number of cores used for parallel + computing. When \code{ncores=NULL} (default), the number \code{r} of + available cores is detected, and \code{(r-1)} of them are used + (see \emph{Details}). + } + \item{monitor}{ + A logical value. If \code{TRUE} progress messages are + printed on screen. + } +} + + +\details{ + The \code{rsc_cv} function performs cross-validation to estimate the + expected Frobenius loss proposed in Bickel and Levina (2008). The + original contribution of Bickel and Levina (2008), and its extension + in Serra et al. (2018), is based on a random + cross-validation algorithm where the training/test size depends on + the sample size \emph{n}. The latter is implemented selecting + \code{cv.type = "ramdom"}, and fixing an appropriate number \code{R} of random + train/test splits. \code{R} should be as large as possible, but + in practice this impacts the computing time strongly for + high-dimensional data sets. + + Although Serra et al. (2018) showed that the random cross-validation + of Bickel and Levina (2008) works well for the RSC estimator, + subsequent experiments suggested that repeated K-fold cross-validation + on average produces better results. Repeated K-fold cross-validation + is implemented with the default \code{cv.type = "kfold"}. In this case + \code{K} defines the number of \emph{folds}, while \code{R} defines + the number of times that the K-fold cross-validation is repeated with + \code{R} independent shuffles of the original data. Selecting + \code{R=1} and \code{K=10} one performs the standard 10-fold + cross-validation. Ten replicates (\code{R=10}) of the K-fold + cross-validation are generally sufficient to obtain reasonable + estimates of the underlying loss, but for extremely high-dimensional + data \code{R} may be varied to speed up calculations. + + On multi-core hardware the cross-validation is executed in parallel + setting \code{ncores}. The parallelism is implemented on the + total number of data splits, that is \code{R} for the random + cross-validation, and \code{R*K} for the repeated K-fold + cross-validation. The software is optimized so that generally the + total computing time scales almost linearly with the number of + available computer cores (\code{ncores}). + + For both the random and the K-fold cross-validation it is computed the + normalized version of the expected squared Frobenius loss proposed in + Bickel and Levina (2008). The normalization is such + that the squared Frobenius norm of the identity matrix equals to 1 + whatever is its dimension. + + Two optimal threshold selection types are reported with flags (see + \emph{Value} section below): \code{"minimum"} and + \code{"minimum1se"}. The flag \code{"minimum"} denotes the threshold + value that minimizes the average loss. The flag \code{"minimum1se"} + implements the so called + \emph{1-SE rule}: this is the maximum threshold value such that the + corresponding average loss is within \emph{1-standard-error} with + respect to the threshold that minimizes the average loss + (that is the one corresponding to the \code{"minimum"} flag). + + Since unbiased standard errors for the K-fold cross-validation are + impossible to compute (see Bengio and Grandvalet, 2004), when + \code{cv.type="kfold"} the reported standard errors have to be + considered as a downward biased approximation. +} + + + +\value{ + An S3 object of class \code{'cv_rsc'} with the following components: + \item{rmadvec}{ + A vector containing the lower triangle of the underlying RMAD + matrix. + } + \item{varnames}{ + A character vector if variable names are available for the input + data set \code{x}. Otherwise this is \code{NULL}. + } + \item{loss}{ + A data.frame reporting cross-validation estimates. Columns of + \code{loss} are as follows: \code{loss$Threshold} is the threshold value; + \code{loss$Average} is averaged loss; \code{loss$SE} is the standard error + for the average loss; \code{loss$Flag="minimum"} denotes the threshold + achieving the minimum average loss; \code{loss$Flag = "*"} denotes threshold + values such that the average loss is within \emph{1-standard-error} + with respect to the \code{"minimum"} solution. + } + \item{minimum}{ + A numeric value. This is the minimum of the average loss. This + corresponds to the flag \code{"minimum"} in the loss component + above (see \emph{Details}). + } + \item{minimum1se}{ + A numeric value. This is the largest threshold such that the + corresponding \code{flag = "*"}. In practice this selects the + optimal threshold based on the \emph{1-SE rule} discussed in the + \emph{Details} Section above. + } +} + + + + +\section{References}{ + Bengio, Y., and Grandvalet, Y. (2004). No unbiased estimator of the + variance of k-fold cross-validation. \emph{Journal of Machine Learning + Research}, 5(Sep), 1089-1105. + + Bickel, P. J., and Levina, E. (2008). Covariance regularization by + thresholding. The \emph{Annals of Statistics}, 36(6), 2577-2604. + doi:10.1214/08-AOS600 + + Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018). + Robust and sparsecorrelation matrix estimation for the analysis of + high-dimensional genomics data. \emph{Bioinformatics}, 34(4), + 625-634. doi:10.1093/bioinformatics/btx642 +} + + + +\seealso{ + \code{rsc}, \code{plot.rsc_cv} +} + + + + + + + +\examples{ +\donttest{ +## simulate a random sample from a multivariate Cauchy distribution +## note: example in high-dimension are obtained increasing p +set.seed(1) +n <- 100 # sample size +p <- 10 # dimension +dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p) +colnames(dat) <- paste0("Var", 1:p) + + +## perform 10-fold cross-validation repeated R=10 times +## note: for multi-core machines experiment with 'ncores' +set.seed(2) +a <- rsc_cv(x = dat, R = 10, K = 10) +a + + +## threshold selection: note that here, knowing the sampling designs, +## we would like to threshold any correlation larger than zero in +## absolute value +## +a$minimum ## "minimum" flagged solution +a$minimum1se ## "minimum1se" flagged solution + +## plot the cross-validation estimates +plot(a) + +## to obtain the RSC matrix we pass 'a' to the rsc() function +b <- rsc(cv = a, threshold = "minimum") +b + +d <- rsc(cv = a, threshold = "minimum1se") +d + +## since the object 'a' stores the RMAD underlying estimator, we can +## apply thresholding at any level without re-estimating the RMAD +## matrix +e <- rsc(cv = a, threshold = 0.5) +e +} +} + |