1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
|
\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{
## 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, ncores = 1)
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
}
|
