diff options
Diffstat (limited to 'man')
| -rwxr-xr-x | man/plot.cv_rsc.Rd | 84 | ||||
| -rwxr-xr-x | man/rmad.Rd | 136 | ||||
| -rwxr-xr-x | man/rsc.Rd | 118 | ||||
| -rwxr-xr-x | man/rsc_cv.Rd | 228 |
4 files changed, 566 insertions, 0 deletions
diff --git a/man/plot.cv_rsc.Rd b/man/plot.cv_rsc.Rd new file mode 100755 index 0000000..0d7b2c2 --- /dev/null +++ b/man/plot.cv_rsc.Rd @@ -0,0 +1,84 @@ +\name{plot.rsc_cv} +\alias{plot.rsc_cv} + +\title{ + Plot method for rsc_cv objects +} + +\description{ + Plot the cross-validation estimates of the Frobenius loss. +} + + +\usage{ + \method{plot}{rsc_cv}(x, \dots) +} + + + + +\arguments{ + \item{x}{ + Output from \code{\link{rsc_cv}}, that is an S3 object of class \code{"rsc_cv"}. + } + \item{\dots}{ + additional arguments passed to \code{\link[graphics]{plot.default}}. + } +} + + +\value{ + Plot the Frobenius loss estimated via cross-validation (y-axis) vs + threshold values (x-axis). The dotted blue line represents the average + expected normalized Frobenius loss, while the vertical segments + around the average are \emph{1-standard-error} error bars + (see \emph{Details} in \code{\link{rsc_cv}}. + + The vertical dashed red line identifies the minimum of the average + loss, that is the optimal threshold flagged as \code{"minimum"}. The + vertical dashed green line identifies the optimal selection flagged + as \code{"minimum1se"} in the output of \code{\link{rsc_cv}} (see + \emph{Details} in \code{\link{rsc_cv}}). +} + + + + +\section{References}{ + 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{\link{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 + +## plot the cross-validation estimates +plot(a) + +## pass additional parameters to graphics::plot +plot(a , cex = 2) +} +} diff --git a/man/rmad.Rd b/man/rmad.Rd new file mode 100755 index 0000000..01247cf --- /dev/null +++ b/man/rmad.Rd @@ -0,0 +1,136 @@ +\name{rmad} + +\alias{rmad} + +\title{RMAD correlation matrix} + +\description{ + Compute the RMAD robust correlation matrix proposed in Serra et + al. (2018) based on the robust correlation coefficient proposed in + Pasman and Shevlyakov (1987). +} + + +\usage{ + rmad(x , y = NULL, na.rm = FALSE , even.correction = FALSE) +} + + +\arguments{ + \item{x}{ + A numeric vector, a matrix or a data.frame. If \code{x} is a matrix + or a data.frame, rows of \code{x} correspond to sample units + and columns correspond to variables. If \code{x} is a numerical + vector, and \code{y} is not \code{NULL}, the RMAD correlation + coefficient between \code{x} and \code{y} is computed. Categorical + variables are not allowed. + } + \item{y}{ + A numerical vector if not \code{NULL}. If both \code{x} and \code{y} + are numerical vectors, the RMAD correlation coefficient between + \code{x} and \code{y} is computed. + } + \item{na.rm}{ + A logical value, if \code{TRUE} sample observation + containing \code{NA} values are excluded (see \emph{Details}). + } + \item{even.correction}{ + A logical value, if \code{TRUE} a correction + for the calculation of the medians is applied to reduce the bias + when the number of samples even (see \emph{Details}). + } +} + + +\details{ + The \code{rmad} function computes the correlation matrix based on the + pairwise robust correlation coefficient of Pasman and Shevlyakov + (1987). This correlation coefficient is based on repeated median + calculations for all pairs of variables. This is a computational + intensive task when the number of variables (that is \code{ncol(x)}) + is large. + + The software is optimized for large dimensional data sets, the median + is approximated as the central observation obtained based on the + \emph{introselect} sorting algorithm of Musser (1997) implemented in + Fortran 95 language. For small samples this may be a crude + approximation, however, it makes the computational cost feasible for + high-dimensional data sets. With the option \code{even.correction + = TRUE} a correction is applied to reduce the bias for data sets with + an even number of samples. Although \code{even.correction = TRUE} + has a small computational cost for each pair of variables, it is + suggested to use the default \code{even.correction = FALSE} for large + dimensional data sets. + + The function can handle a data matrix with missing values (\code{NA} + records). If \code{na.rm = TRUE} then missing values are handled by + casewise deletion (and if there are no complete cases, an error is + returned). In practice, if \code{na.rm = TRUE} all rows of + \code{x} that contain at least an \code{NA} are removed. + + Since the software is optimized to work with high-dimensional data sets, + the output RMAD matrix is packed into a storage efficient format + using the \code{"dspMatrix"} S4 class from the \code{\link{Matrix}} + package. The latter is specifically designed for dense real symmetric + matrices. A sparse correlation matrix can be obtained applying + thresholding using the \code{\link{rsc_cv}} and \code{\link{rsc}}. +} + + + +\value{ + \item{If \code{x} is a matrix or a data.frame}{ + Returns a correlation matrix of class \code{"dspMatrix"} (S4 class object) + as defined in the \code{\link{Matrix}} package. + } + \item{If \code{x} and \code{y} are numerical vectors}{ + Returns a numerical value, that is the RMAD correlation coefficient + between \code{x} and \code{y}. + } +} + + + +\section{References}{ + Musser, D. R. (1997). Introspective sorting and selection algorithms. + \emph{Software: Practice and Experience}, 27(8), 983-993. + + Pasman,V. and Shevlyakov,G. (1987). Robust methods of estimation of + correlation coefficient. \emph{Automation Remote Control}, 48, 332-340. + + 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_cv}, \code{rsc} +} + + + + + + + +\examples{ +## simulate a random sample from a multivariate Cauchy distribution +set.seed(1) +n <- 100 # sample size +p <- 7 # dimension +dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p) +colnames(dat) <- paste0("Var", 1:p) + + +## compute the rmad correlation coefficient between dat[,1] and dat[,2] +a <- rmad(x = dat[,1], y = dat[,2]) + + +## compute the RMAD correlaiton matrix +b <- rmad(x = dat) +b +} + diff --git a/man/rsc.Rd b/man/rsc.Rd new file mode 100755 index 0000000..cd252a6 --- /dev/null +++ b/man/rsc.Rd @@ -0,0 +1,118 @@ +\name{rsc} + +\alias{rsc} + +\title{Robust and Sparse Correlation Matrix Estimator} + +\description{ + Compute the Robust and Sparse Correlation Matrix (RSC) estimator + proposed in Serra et al. (2018). +} + + +\usage{ + rsc(cv, threshold = "minimum") +} + + +\arguments{ + \item{cv}{ + An S3 object of class \code{"rsc_cv"} (see \code{\link{rsc_cv}}). + } + \item{threshold}{ + Threshold parameter to compute the RSC estimate. This + is a numeric value taken onto the interval (0,1), or it is + equal to \code{"minimum"} or \code{"minimum1se"} for selecting the + optimal threshold according to the selection performed in + \code{\link{rsc_cv}}. + } +} + + +\details{ + The setting \code{threshold = "minimum"} or \code{threshold = + "minimum1se"} applies thresholding according to the criteria + discussed in the \emph{Details} section in \code{\link{rsc_cv}}. + When \code{cv} is obtained using \code{\link{rsc_cv}} with + \code{cv.type = "random"}, the default settings for \code{\link{rsc}} + implements exactly the RSC estimator proposed in Serra et al., + (2018). + + Although \code{threshold = "minimum"} is the default choice, in + high-dimensional situations \code{threshold = "minimum1se"} usually + provides a more parsimonious representation of the correlation + structure. Since the underlying RMAD matrix is passed through the + \code{cv} input, any other hand-tuned threshold to the RMAD matrix + can be applied without significant additional computational + costs. The latter can be done setting \code{threshold} to any value + onto the (0,1) interval. + + The software is optimized to handle high-dimensional data sets, + therefore, the output RSC matrix is packed into a storage efficient + sparse format using the \code{"dsCMatrix"} S4 class from the + \code{\link{Matrix}} package. The latter is specifically designed for + sparse real symmetric matrices. +} + + + +\value{ + Returns a sparse correlaiton matrix of class \code{"dsCMatrix"} + (S4 class object) as defined in the \code{\link{Matrix}} package. +} + + + + +\section{References}{ + 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{\link{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 + +## obtain the RSC matrix with "minimum" flagged solution +b <- rsc(cv = a, threshold = "minimum") +b + +## obtain the RSC matrix with "minimum1se" flagged solution +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 +} +} + 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 +} +} + |