cv.missoNet {missoNet}

Cross-validation for missoNet

Description

This function performs k-fold cross-validation for ‘missoNet’. The regularization path is computed for all possible combinations of values given in the two regularization parameter sequences, namely \lambda_B and \lambda_\Theta. ‘cv.missoNet’ will select the most suitable model among all cross-validated fits along the path. See the details of ‘missoNet’ for the model definition. To help users, the ‘cv.missoNet’ function is designed to automatically determine the likely ranges of the regularization parameters over which the cross-validation searches.

Usage

cv.missoNet(
  X,
  Y,
  kfold = 5,
  rho = NULL,
  lambda.Beta = NULL,
  lambda.Theta = NULL,
  lamBeta.min.ratio = NULL,
  lamTheta.min.ratio = NULL,
  n.lamBeta = NULL,
  n.lamTheta = NULL,
  lamBeta.scale.factor = 1,
  lamTheta.scale.factor = 1,
  Beta.maxit = 1000,
  Beta.thr = 1e-04,
  eta = 0.8,
  Theta.maxit = 1000,
  Theta.thr = 1e-04,
  eps = 1e-08,
  penalize.diagonal = TRUE,
  diag.penalty.factor = NULL,
  standardize = TRUE,
  standardize.response = TRUE,
  fit.1se = FALSE,
  fit.relax = FALSE,
  permute = TRUE,
  with.seed = NULL,
  parallel = FALSE,
  cl = NULL,
  verbose = 1
)

Arguments

X

Numeric predictor matrix (n\times p): columns correspond to predictor variables and rows correspond to samples. Missing values are not allowed. There is no need for centering or scaling of the variables. 'X' should not include a column of ones for an intercept.

Y

Numeric response matrix (n\times q): columns correspond to response variables and rows correspond to samples. Missing values should be coded as either 'NA's or 'NaN's. There is no need for centering or scaling of the variables.

kfold

Number of folds for cross-validation – the default is '5'.

rho

(Optional) A scalar or a numeric vector of length q: the elements are user-supplied probabilities of missingness for the response variables. The default is 'rho = NULL' and the program will compute the empirical missing rates for each of the columns of 'Y' and use them as the working missing probabilities. The default setting should suffice in most cases; misspecified missing probabilities would introduce biases into the model.

lambda.Beta

(Optional) Numeric vector: a user-supplied sequence of non-negative values for {\lambda_B} penalizing the elements of the coefficient matrix \mathbf{B} among which the cross-validation procedure searches. The default is 'lambda.Beta = NULL', in which case the program computes an appropriate range of \lambda_B values using 'n.lamBeta' and 'lamBeta.min.ratio'. Supplying a vector overrides this default. Note that the supplied sequence will be automatically arranged, internally, in a descending order.

lambda.Theta

(Optional) Numeric vector: a user-supplied sequence of non-negative values for {\lambda_\Theta} penalizing the (off-diagonal) elements of the precision matrix \mathbf{\Theta} among which the cross-validation procedure searches. The default is 'lambda.Theta = NULL', in which case the program computes an appropriate range of \lambda_\Theta values using 'n.lamTheta' and 'lamTheta.min.ratio'. Supplying a vector overrides this default. Note that the supplied sequence will be automatically arranged, internally, in a descending order.

lamBeta.min.ratio

The smallest value of \lambda_B is calculated as the data-derived \mathrm{max}(\lambda_B) multiplied by 'lamBeta.min.ratio'. The default depends on the sample size, n, relative to the number of predictors, p. If n > p, the default is '1.0E-4', otherwise it is '1.0E-3'. A very small value of 'lamBeta.min.ratio' may significantly increase runtime and lead to a saturated fit in the n \leq p case. This is only needed when 'lambda.Beta = NULL'.

lamTheta.min.ratio

The smallest value of \lambda_\Theta is calculated as the data-derived \mathrm{max}(\lambda_\Theta) multiplied by 'lamTheta.min.ratio'. The default depends on the sample size, n, relative to the number of responses, q. If n > q, the default is '1.0E-4', otherwise it is '1.0E-3'. A very small value of 'lamTheta.min.ratio' may significantly increase runtime and lead to a saturated fit in the n \leq q case. This is only needed when 'lambda.Theta = NULL'.

n.lamBeta

The number of \lambda_B values. If n > p, the default is '40', otherwise it is '30'. Avoid supplying an excessively large number since the program will fit ('n.lamBeta' * 'n.lamTheta') models in total for each fold of the cross-validation. Typically we suggest 'n.lamBeta' = -log10('lamBeta.min.ratio') * c, where c \in [10, 20]. This is only needed when 'lambda.Beta = NULL'.

n.lamTheta

The number of \lambda_\Theta values. If n > q, the default is '40', otherwise it is '30'. Avoid supplying an excessively large number since the program will fit ('n.lamBeta' * 'n.lamTheta') models in total for each fold of the cross-validation. Typically we suggest 'n.lamTheta' = -log10('lamTheta.min.ratio') * c, where c \in [10, 20]. This is only needed when 'lambda.Theta = NULL'.

lamBeta.scale.factor

A positive multiplication factor for scaling the entire \lambda_B sequence; the default is '1'. A typical usage is when the magnitudes of the auto-computed \lambda_B values are inappropriate. For example, this factor would be needed if the optimal value of \lambda_B selected by the cross-validation (i.e. {\lambda_B}_\mathrm{min} with the minimum cross-validated error) approaches either boundary of the search range. This is only needed when 'lambda.Beta = NULL'.

lamTheta.scale.factor

A positive multiplication factor for scaling the entire \lambda_\Theta sequence; the default is '1'. A typical usage is when the magnitudes of the auto-computed \lambda_\Theta values are inappropriate. For example, this factor would be needed if the optimal value of \lambda_\Theta selected by the cross-validation (i.e. {\lambda_\Theta}_\mathrm{min} with the minimum cross-validated error) approaches either boundary of the search range. This is only needed when 'lambda.Theta = NULL'.

Beta.maxit

The maximum number of iterations of the fast iterative shrinkage-thresholding algorithm (FISTA) for updating \hat{\mathbf{B}}. The default is 'Beta.maxit = 1000'.

Beta.thr

The convergence threshold of the FISTA algorithm for updating \hat{\mathbf{B}}; the default is 'Beta.thr = 1.0E-4'. Iterations stop when the absolute parameter change is less than ('Beta.thr' * sum(abs(\hat{\mathbf{B}}))).

eta

The backtracking line search shrinkage factor; the default is 'eta = 0.8'. Most users will be able to use the default value, some experienced users may want to tune 'eta' according to the properties of a specific dataset for a faster convergence of the FISTA algorithm. Note that 'eta' must be in (0, 1).

Theta.maxit

The maximum number of iterations of the ‘glasso’ algorithm for updating \hat{\mathbf{\Theta}}. The default is 'Theta.maxit = 1000'.

Theta.thr

The convergence threshold of the ‘glasso’ algorithm for updating \hat{\mathbf{\Theta}}; the default is 'Theta.thr = 1.0E-4'. Iterations stop when the average absolute parameter change is less than ('Theta.thr' * ave(abs(offdiag(\hat{\mathbf{\Sigma}})))), where \hat{\mathbf{\Sigma}} denotes the empirical working covariance matrix.

eps

A numeric tolerance level for the L1 projection of the empirical covariance matrix; the default is 'eps = 1.0E-8'. The empirical covariance matrix will be projected onto a L1 ball to have min(eigen(\hat{\mathbf{\Sigma}})$value) == 'eps', if any of the eigenvalues is less than the specified tolerance. Most users will be able to use the default value.

penalize.diagonal

Logical: should the diagonal elements of \mathbf{\Theta} be penalized? The default is 'TRUE'.

diag.penalty.factor

Numeric: a separate penalty multiplication factor for the diagonal elements of \mathbf{\Theta} when 'penalize.diagonal = TRUE'. \lambda_\Theta is multiplied by this number to allow a differential shrinkage of the diagonal elements. The default is 'NULL' and the program will guess a value based on an initial estimate of \mathbf{\Theta}. This factor could be '0' for no shrinkage (equivalent to 'penalize.diagonal = FALSE') or '1' for an equal shrinkage.

standardize

Logical: should the columns of 'X' be standardized so each has unit variance? The default is 'TRUE'. The estimated results will always be returned on the original scale. ‘cv.missoNet’ computes appropriate \lambda sequences relying on standardization, if 'X' has been standardized prior to fitting the model, you might not wish to standardize it inside the algorithm.

standardize.response

Logical: should the columns of 'Y' be standardized so each has unit variance? The default is 'TRUE'. The estimated results will always be returned on the original scale. ‘cv.missoNet’ computes appropriate \lambda sequences relying on standardization, if 'Y' has been standardized prior to fitting the model, you might not wish to standardize it inside the algorithm.

fit.1se

Logical: the default is 'FALSE'. If 'TRUE', two additional models will be fitted with the largest values of \lambda_B and \lambda_\Theta respectively at which the cross-validated error is within one standard error of the minimum.

fit.relax

Logical: the default is 'FALSE'. If 'TRUE', the program will re-estimate the edges in the active set (i.e. nonzero off-diagonal elements) of the network structure \hat{\mathbf{\Theta}} without penalization (\lambda_\Theta=0). This debiased estimate of \mathbf{\Theta} could be useful for some interdependency analyses. WARNING: there may be convergence issues if the empirical covariance matrix is not of full rank (e.g. n < q)).

permute

Logical: should the subject indices for the cross-validation be permuted? The default is 'TRUE'.

with.seed

A random number seed for the permutation.

parallel

Logical: the default is 'FALSE'. If 'TRUE', the program uses clusters to compute the cross-validation folds in parallel. Must register parallel clusters beforehand, see examples below.

cl

A cluster object created by ‘parallel::makeCluster’ for parallel evaluations. This is only needed when 'parallel = TRUE'.

verbose

Value of '0', '1' or '2'. 'verbose = 0' – silent; 'verbose = 1' (the default) – limited tracing with progress bars; 'verbose = 2' – detailed tracing. Note that displaying the progress bars slightly increases the computation overhead compared to the silent mode. The detailed tracing should be used for monitoring progress only when the program runs extremely slowly, and it is not supported under 'parallel = TRUE'.

Details

The ‘cv.missoNet’ function fits ‘missoNet’ models ('kfold' * 'n.lamBeta' * 'n.lamTheta') times in the whole cross-validation process. That is, for the kth-fold (k=1,...,K) computation, the models are fitted at each of the all ('n.lamBeta' * 'n.lamTheta') possible combinations of the regularization parameters (\lambda_B, \lambda_\Theta), with the kth fold of the training data omitted. The errors are accumulated, and the averaged errors as well as the standard deviations are computed over all folds. Note that the results of ‘cv.missoNet’ are random, since the samples are randomly split into k-folds. Users can eliminate this randomness by setting 'permute = FALSE', or by explicitly assigning a seed to the permutation of samples.

A user-supplied sequence for {\lambda_B} and/or {\lambda_\Theta} is permitted, otherwise the program computes an appropriate range of values for the regularization parameters using other control arguments. Note that ‘cv.missoNet’ standardizes 'X' and 'Y' to have unit variances before computing its \lambda sequences (and then unstandardizes the resulting coefficients); if you wish to reproduce/compare results with those of other softwares, it is best to supply pre-standardized 'X' and 'Y'. If the algorithm is running slowly, track its progress with 'verbose = 2'. In most cases, choosing a sparser grid for the regularization parameters (e.g. smaller 'n.lamBeta' and/or 'n.lamTheta') or setting 'Beta.thr = 1.0E-3' (or even bigger) allows the algorithm to make faster progress.

After cross-validation, the regression coefficient matrix \mathbf{B} and the precision matrix \mathbf{\Theta} can be estimated at three special \lambda pairs, by reapplying ‘missoNet’ to the entire training dataset:

1. "lambda.min" := ({\lambda_B}_\mathrm{min}, {\lambda_\Theta}_\mathrm{min}), at which the minimum mean cross-validated error is achieved;

2. "lambda.1se.Beta" := ({\lambda_B}_\mathrm{1se}, {\lambda_\Theta}_\mathrm{min}), where {\lambda_B}_\mathrm{1se} is the largest \lambda_B at which the mean cross-validated error is within one standard error of the minimum;

3. "lambda.1se.Theta" := ({\lambda_B}_\mathrm{min}, {\lambda_\Theta}_\mathrm{1se}), where {\lambda_\Theta}_\mathrm{1se} is the largest \lambda_\Theta at which the mean cross-validated error is within one standard error of the minimum.

The corresponding estimates, along with the \lambda values, are stored in three separate lists inside the returned object: 'est.min', 'est.1se.B' and 'est.1se.Tht' ('est.1se.B' and 'est.1se.Tht' are 'NULL' unless the argument 'fit.1se = TRUE' when calling ‘cv.missoNet’).

The ‘cv.missoNet’ function returns an R object of S3 class 'cv.missoNet' for which there are a set of accessory functions available. The plotting function ‘plot.cv.missoNet’ can be used to graphically identify the optimal pair of the regularization parameters, and the prediction function ‘predict.cv.missoNet’ can be used to make predictions of response values given new input 'X'. See the vignette for examples.

Value

This function returns a 'cv.missoNet' object containing a named 'list' with all the ingredients of the cross-validated fit:

est.min

A 'list' of results estimated at "lambda.min" := ({\lambda_B}_\mathrm{min}, {\lambda_\Theta}_\mathrm{min}) that gives the minimum mean cross-validated error. It consists of the following components:

• Beta: the penalized estimate of the regression coefficient matrix \hat{\mathbf{B}} (p\times q).

• Theta: the penalized estimate of the precision matrix \hat{\mathbf{\Theta}} (q\times q).

• mu: a vector of length q storing the model intercept \hat{\mu}.

• lambda.Beta: the value of \lambda_B (i.e. {\lambda_B}_\mathrm{min}) used to fit the model.

• lambda.Theta: the value of \lambda_\Theta (i.e. {\lambda_\Theta}_\mathrm{min}) used to fit the model.

• relax.net: the relaxed (debiased) estimate of the conditional network structure \hat{\mathbf{\Theta}}_\mathrm{rlx} (q\times q) if 'fit.relax = TRUE' when calling ‘cv.missoNet’.

est.1se.B

A 'list' of results estimated at "lambda.1se.Beta" := ({\lambda_B}_\mathrm{1se}, {\lambda_\Theta}_\mathrm{min}) if 'fit.1se = TRUE' when calling ‘cv.missoNet’. "lambda.1se.Beta" refers to the largest \lambda_B at which the mean cross-validated error is within one standard error of the minimum, by fixing \lambda_\Theta at {\lambda_\Theta}_\mathrm{min}. This 'list' consists of the same components as 'est.min'.

est.1se.Tht

A 'list' of results estimated at "lambda.1se.Theta" := ({\lambda_B}_\mathrm{min}, {\lambda_\Theta}_\mathrm{1se}) if 'fit.1se = TRUE' when calling ‘cv.missoNet’. "lambda.1se.Theta" refers to the largest \lambda_\Theta at which the mean cross-validated error is within one standard error of the minimum, by fixing \lambda_B at {\lambda_B}_\mathrm{min}. This 'list' consists of the same components as 'est.min'.

rho

A vector of length q storing the working missing probabilities for the q response variables.

fold.index

The subject indices identifying which fold each observation is in.

lambda.Beta.vec

A flattened vector of length ('n.lamBeta' * 'n.lamTheta') storing the \lambda_B values along the regularization path. More specifically, 'lambda.Beta.vec' = rep('lambda.Beta', each = 'n.lamTheta').

lambda.Theta.vec

A flattened vector of length ('n.lamBeta' * 'n.lamTheta') storing the \lambda_\Theta values along the regularization path. More specifically, 'lambda.Theta.vec' = rep('lambda.Theta', times = 'n.lamBeta').

cvm

A flattened vector of length ('n.lamBeta' * 'n.lamTheta') storing the (standardized) mean cross-validated errors along the regularization path.

cvup

Upper cross-validated errors.

cvlo

Lower cross-validated errors.

penalize.diagonal

Logical: whether the diagonal elements of \mathbf{\Theta} were penalized.

diag.penalty.factor

The additional penalty multiplication factor for the diagonal elements of \mathbf{\Theta} when 'penalize.diagonal' was returned as 'TRUE'.

Author(s)

Yixiao Zeng [e-mail], Celia M.T. Greenwood and Archer Yi Yang.

Examples

## Simulate a dataset with response values missing completely at random (MCAR), 
## the overall missing rate is around 10%.
set.seed(123)  # reproducibility
sim.dat <- generateData(n = 300, p = 50, q = 20, rho = 0.1, missing.type = "MCAR")
tr <- 1:240  # training set indices
tst <- 241:300  # test set indices
X.tr <- sim.dat$X[tr, ]  # predictor matrix
Y.tr <- sim.dat$Z[tr, ]  # corrupted response matrix


## Perform a five-fold cross-validation WITH specified 'lambda.Beta' and 'lambda.Theta'.
## 'standardize' and 'standardize.response' are 'TRUE' by default.
lamB.vec <- 10^(seq(from = 0, to = -1, length.out = 5))
lamTht.vec <- 10^(seq(from = 0, to = -1, length.out = 5))
cvfit <- cv.missoNet(X = X.tr, Y = Y.tr, kfold = 5,
                     lambda.Beta = lamB.vec, lambda.Theta = lamTht.vec)


## Perform a five-fold cross-validation WITHOUT specified 'lambda.Beta' and 'lambda.Theta'.
## In this case, a grid of 'lambda.Beta' and 'lambda.Theta' values in a (hopefully) reasonable 
## range will be computed and used by the program.
## 
## < This simplest command should be a good start for most analyses. >
cvfit <- cv.missoNet(X = X.tr, Y = Y.tr, kfold = 5)


## Alternatively, compute the cross-validation folds in parallel on a cluster with 2 cores.
## 
## 'fit.1se = TRUE' tells the program to make additional estimations of the parameters at the 
## largest value of 'lambda.Beta' / 'lambda.Theta' that gives the most regularized model such 
## that the cross-validated error is within one standard error of the minimum.
cl <- parallel::makeCluster(min(parallel::detectCores()-1, 2))
cvfit <- cv.missoNet(X = X.tr, Y = Y.tr, kfold = 5, fit.1se = TRUE,
                     parallel = TRUE, cl = cl,
                     permute = TRUE, with.seed = 486)  # permute with seed for reproducibility
parallel::stopCluster(cl)


## Use PRE-STANDARDIZED training data if you wish to compare the results with other softwares. 
X.tr.std <- scale(X.tr, center = TRUE, scale = TRUE)
Y.tr.std <- scale(Y.tr, center = TRUE, scale = TRUE)
cvfit.std <- cv.missoNet(X = X.tr.std, Y = Y.tr.std, kfold = 5,
                         standardize = FALSE, standardize.response = FALSE)


## Plot the (standardized) mean cross-validated errors in a heatmap.
plot(cvfit, type = "cv.heatmap")

## Plot the (standardized) mean cross-validated errors in a 3D scatterplot.
plot(cvfit, type = "cv.scatter", plt.surf = TRUE)


## Extract the estimates at "lambda.min".
Beta.hat1 <- cvfit$est.min$Beta
Theta.hat1 <- cvfit$est.min$Theta

## Extract the estimates at "lambda.1se.Beta" (if 'fit.1se' = TRUE).
Beta.hat2 <- cvfit$est.1se.B$Beta
Theta.hat2 <- cvfit$est.1se.B$Theta

## Extract the estimates at "lambda.1se.Theta" (if 'fit.1se' = TRUE).
Beta.hat3 <- cvfit$est.1se.Tht$Beta
Theta.hat3 <- cvfit$est.1se.Tht$Theta


## Make predictions of response values on the test set.
newy1 <- predict(cvfit, newx = sim.dat$X[tst, ], s = "lambda.min")
newy2 <- predict(cvfit, newx = sim.dat$X[tst, ], s = "lambda.1se.Beta")  # 'fit.1se' = TRUE
newy3 <- predict(cvfit, newx = sim.dat$X[tst, ], s = "lambda.1se.Theta")  # 'fit.1se' = TRUE