Introduction

missoNet is a package that fits penalized multi-task regression – that is, with multiple correlated tasks or response variables – to simultaneously estimate the coefficients of a set of predictor variables for all tasks and the conditional response network structure given all predictors, via penalized maximum likelihood in an undirected conditional Gaussian graphical model. In contrast to most penalized multi-task regression (conditional graphical lasso) methods, missoNet has the capability of obtaining estimates even when the response data is corrupted by missing values. The method automatically enjoys the theoretical and computational benefits of convexity, and returns solutions that are comparable/close to the estimates without any missing values.

missoNet assumes the model $$ \mathbf{Y} = \mathbf{X}\mathbf{B}^* + \mathbf{E},\ \ \mathbf{E}_{i\cdot} \stackrel{iid}{\sim} \mathcal{MVN}(0_q,(\mathbf{\Theta}^*)^{-1})\ \ \forall i = 1,...,n, $$ where Y ∈ ℝ^n × q and X ∈ ℝ^n × p represent the centered response data matrix and predictor data matrix, respectively (thus the intercepts are omitted). E_i⋅ ∈ ℝ^q is the ith row (ith sample) of the error matrix and is drawn from a multivariate Gaussian distribution. n, p, q denote the sample size, the number of predictors and the number of responses, respectively. The regression coefficient matrix B^* ∈ ℝ^p × q and the precision (inverse covariance) matrix Θ^* ∈ ℝ^q × q are parameters to be estimated in this model. Note that the entries of Θ^* have a one-to-one correspondence with partial correlations. That is, the random variables Y_j and Y_k are conditionally independent (i.e. Θ_jk^* = 0) given X and all other remaining variables in Y iff the corresponding partial correlation coefficient is zero.

To estimate the parameters, missoNet solves a penalized MLE problem: $$ (\hat{\mathbf{\Theta}},\hat{\mathbf{B}}) = \underset{\mathbf{\Theta} \succeq 0,\ \mathbf{B}}{\text{argmin}}\ \mathrm{tr}\left[\frac{1}{n}(\mathbf{Y} - \mathbf{XB})^\top(\mathbf{Y} - \mathbf{XB}) \mathbf{\Theta}\right] - \mathrm{log}|\mathbf{\Theta}| + \lambda_{\Theta}(\|\mathbf{\Theta}\|_{1,\mathrm{off}} + 1_{n\leq \mathrm{max}(p,q)} \|\mathbf{\Theta}\|_{1,\mathrm{diag}}) + \lambda_{B}\|\mathbf{B}\|_1. $$ B^* that mapping predictors to responses and Θ^* that revealing the responses’ conditional dependencies are estimated for the lasso penalties – λ_B and λ_Θ, over a grid of values (on the log scale) covering the entire range of possible solutions. To learn more about sparse multi-task regression with inverse covariance estimation using datasets without missing values, see MRCE.

However, missingness in real data is inevitable. Most penalized multi-task regression / conditional graphical lasso methods either assume the data is fully-observed (eg. MRCE), or deal with missing data through a likelihood-based method such as the EM algorithm (e.g, cglasso). An important challenge in this context is the possible non-convexity of the associated optimization problem when there is missing data, as well as the high computational cost from iteratively updating the estimations for parameters.

missoNet aims to handle a specific class of datasets where the response matrix Y contains data which is missing at random (MAR). To use missoNet, users do not need to possess additional knowledge for pre-processing missing data (e.g. imputation) nor for selection of regularization parameters (λ_B and λ_Θ), the method provides a unified framework for automatically solving a convex modification of the multi-task learning problem defined above, using corrupted datasets.

The package provides an integrated set of core routines for data simulation, model fitting and cross-validation, visualization of results, as well as predictions in new data. The function arguments are in the same style as those of glmnet, making it easy for experienced users to get started.

Installation

To install the package missoNet from CRAN, type the following command in the R console:

install.packages("missoNet")

Or install the development version of missoNet from GitHub:

if(!require("devtools")) {
  install.packages("devtools")
}
devtools::install_github("yixiao-zeng/missoNet")

Quick Start

The purpose of this section is to give users a general overview of the functions and their usages. We will briefly go over the main functions, basic operations and outputs.

First, we load the missoNet package:

library(missoNet)

We will use a synthetic dataset for illustration. To generate a set of data with missing response values, we can simply call the built-in function generateData:

## Specify a random seed for reproducibility.
## The overall missing rate in the response matrix is around 10%.
## The missing mechanism can also be "MCAR" or "MNAR".
sim.dat <- generateData(n = 150, p = 15, q = 12, rho = 0.1, missing.type = "MAR", with.seed = 1512)

tr <- 1:120  # training set indices
tst <- 121:150  # test set indices

This command returns an object containing all necessary components for analysis including:

• 'X', 'Y': a predictor matrix X ∈ ℝ^n × p and a complete (clean) response matrix Y ∈ ℝ^n × q, in which rows correspond to samples and columns correspond to variables. The rows of Y is sampled from ℳ𝒱𝒩(XB^*, (Θ^*)⁻¹).

• 'Beta', 'Theta': the true parameters B^* ∈ ℝ^p × q and Θ^* ∈ ℝ^q × q for the generative model that will be estimated later (see the later section for more details on their structures and sparsity).

• 'Z': a corrupted response matrix Z ∈ ℝ^n × q having elements Z_ij (∀ i = 1, ..., n,  j = 1, ..., q) which is either NA or equal to Y_ij.

The probabilities of missingness for the response variables and the missing mechanism are specified by 'rho' and 'missing.type', respectively. 'rho' can be a scalar or a vector of length q. In the code above, 'rho' = 0.1 indicates an overall missing rate of around 10%. In addition to "MAR", the 'missing.type' can also be "MCAR" or "MNAR", see the later section for more details on how missing values are generated under different mechanisms.

NOTE 1: the program only accepts missing values that are coded as either NAs or NaNs.

NOTE 2: although the program will provide estimates even when n ≤ max(p, q), convergence is likely to be rather slow, and the variance of estimation is likely to be excessively large. Therefore, we suggest that both the predictor columns and the response columns be reduced (filtered) if possible, prior to fitting any models. For example in genomics, where X can very wide (i.e. large p), we often filter features based on variance or coefficient of variation.

We can easily visualize how missing data is distributed in the corrupted response matrix Z using package visdat:

Z <- sim.dat$Z  # corrupted response matrix
visdat::vis_miss(as.data.frame(Z))

A single model can be fitted using the most basic call to missoNet:

## Training set
X.tr <- sim.dat$X[tr, ]  # predictor matrix
Z.tr <- sim.dat$Z[tr, ]  # corrupted response matrix

## Using the training set to fit the model.
## 'lambda.Beta' and 'lambda.Theta' are arbitrarily set to 0.1.
## 'verbose' = 0 suppresses printing of messages.
fit <- missoNet(X = X.tr, Y = Z.tr, lambda.Beta = 0.1, lambda.Theta = 0.1, verbose = 0)

However, missoNet should be more commonly used with a grid of values for λ_B and λ_Θ. The two penalty vectors must have the same length, and the missoNet function will be called with each consecutive pair of values – i.e. with the first elements of the two vectors, then with the second elements, etc.

lambda.Beta.vec <- 10^(seq(from = 0, to = -1, length.out = 10))  # 10 values on the log scale, from 1 to 0.1
lambda.Theta.vec <- rep(0.1, 10)  # for each value of 'lambda.Beta', 'lambda.Theta' remains constant at 0.1
fit_list <- missoNet(X = X.tr, Y = Z.tr, lambda.Beta = lambda.Beta.vec, lambda.Theta = lambda.Theta.vec, verbose = 0)

The command above returns a sequence of models for users to choose from. In many cases, users may prefer that the software selects the best fit among them. Cross-validation is perhaps the simplest and most widely used method for this task. cv.missoNet is the main function to do cross-validation, along with supporting methods such as plot and predict.

Here we use cv.missoNet to perform a five-fold cross-validation, samples will be permuted before splitting into multiple folds. For reproducibility, we assign a random seed to the permutation.

## If 'permute' = FALSE, the samples will be split into k-folds in their original orders,
## i.e. the first (n/'kfold') samples will be assigned to the first fold an so on.

cvfit <- cv.missoNet(X = X.tr, Y = Z.tr, kfold = 5,
                     lambda.Beta = lambda.Beta.vec, lambda.Theta = lambda.Theta.vec,
                     permute = TRUE, with.seed = 433, verbose = 0)
#> Warning in boundaryCheck(lambda.Theta = lambda.Theta, lambda.Beta = lambda.Beta, : 
#> The optimal `lambda.Theta` is close to the upper boundary of the search scope, 
#> try to provide a new sequence covering larger values for the following argument:
#> 
#>     1. lambda.Theta

The program has warned us that the range of λ_Θ is inappropriate so that we need to supply a sequence of λ_Θ covering larger values. However, picking a suitable range for such hyperparameters may require experience or a series of trials and errors. Therefore, cv.missoNet provides a smarter way to solve this problem.

Let’s fit the cross-validation model again, this time running all folds in parallel on two CPU cores. The parallelization of missoNet models relies on package parallel, so make sure the parallel clusters have been registered beforehand.

## '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 = Z.tr, kfold = 5,
                     fit.1se = TRUE, parallel = TRUE, cl = cl,
                     permute = TRUE, with.seed = 433, verbose = 1)
parallel::stopCluster(cl)

Note that we have not explicitly specified the regularization parameters λ_B and λ_Θ in the command above. In this case, a grid of λ_B and λ_Θ values in a (hopefully) reasonable range will be computed and used by the program. Users can even provide just one of λ_B and λ_Θ vectors such that the program will compute the other. Generally, it is difficult at the beginning to choose suitable λ sequences, so we suggest users start analysis using cv.missoNet, and let the program compute the appropriate λ_B and λ_Θ values itself, and then automatically pick the optimal combination of them at which the minimum cross-validated error is achieved.

cv.missoNet returns a 'cv.missoNet' object, a list with all the ingredients of the cross-validated fit. We can execute the plot method

## The values of mean cross-validated errors along with upper and lower standard deviation bounds 
## can be accessed via 'cvfit$cvm', 'cvfit$cvup' and 'cvfit$cvlo', respectively.

plot(cvfit)

to visualize the mean cross-validated error in a heatmap style. The white solid box marks out the position of the minimum mean cross-validated error with corresponding λ_B and λ_Θ ("lambda.min" := (λ_Bmin, λ_Θmin)), and the white dashed boxes indicate the largest λ_B and λ_Θ at which the mean cross-validated error is within one standard error of the minimum, by fixing the other one at "lambda.min" (i.e. "lambda.1se.Beta" := (λ_B1se, λ_Θmin) and "lambda.1se.Theta" := (λ_Bmin, λ_Θ1se)). We often use this “one-standard-error” rule when selecting the most parsimonious model; this acknowledges the fact that the cross-validation surface is estimated with error, so error on the side of parsimony is preferable.

We can also plot the errors in a 3D scatter plot (without surfaces):

plot(cvfit, type = "cv.scatter", plt.surf = FALSE)

After cross-validation, B^* and Θ^* can be estimated at three special λ pairs discussed above – "lambda.min", "lambda.1se.Beta" and "lambda.1se.Theta". The corresponding results, along with the λ values, are stored in three separate lists: '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).

Let’s extract the estimates of the model parameters $\hat{\mathbf{B}}$ and $\hat{\mathbf{\Theta}}$ at "lambda.min" and "lambda.1se.Beta" then plot them beside the ground truth B^* and Θ^*:

## Define a plotting function
plot_heatmap <- function(est, col, legend = FALSE, lgd_name, title) {
  return(ComplexHeatmap::Heatmap(est, cluster_rows = FALSE, cluster_columns = FALSE,
                                 col = col, show_heatmap_legend = legend, name = lgd_name,
                                 column_names_gp = grid::gpar(fontsize = 8),
                                 row_names_gp = grid::gpar(fontsize = 8), row_names_side = "left", 
                                 border = TRUE, column_title = title))
}

## Color space
col <- circlize::colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))

## Beta*
Beta_star <- sim.dat$Beta
Beta_star_ht <- plot_heatmap(Beta_star, col, title = expression(paste(bold(Beta), "*")))

## Beta_hat at "lambda.min"
Beta_hat_min <- cvfit$est.min$Beta
Beta_hat_min_ht <- plot_heatmap(Beta_hat_min, col, title = expression(paste(hat(bold(Beta)), " at `lambda.min`")))

## Beta_hat at "lambda.1se.Beta"
Beta_hat_1se <- cvfit$est.1se.B$Beta
Beta_hat_1se_ht <- plot_heatmap(Beta_hat_1se, col, legend = TRUE, lgd_name = "values",
                                title = expression(paste(hat(bold(Beta)), " at `lambda.Beta.1se`")))

## Theta*
Theta_star <- sim.dat$Theta
Theta_star_ht <- plot_heatmap(Theta_star, col, title = expression(paste(bold(Theta), "*")))

## Theta_hat at "lambda.min"
Theta_hat_min <- cvfit$est.min$Theta
Theta_hat_min_ht <- plot_heatmap(Theta_hat_min, col, title = expression(paste(hat(bold(Theta)), " at `lambda.min`")))
  
## Theta_hat at "lambda.1se.Beta"
Theta_hat_1se <- cvfit$est.1se.B$Theta
Theta_hat_1se_ht <- plot_heatmap(Theta_hat_1se, col, legend = TRUE, lgd_name = "values",
                                 title = expression(paste(hat(bold(Theta)), " at `lambda.Beta.1se`")))

## Plot
Beta_star_ht + Beta_hat_min_ht + Beta_hat_1se_ht
Theta_star_ht + Theta_hat_min_ht + Theta_hat_1se_ht

Like other regression methods, predictions can be made based on the fitted 'cv.missoNet' object as well. The code below gives predictions for a new input matrix 'newx' at "lambda.min" ('s' = "lambda.1se.Beta" and 's' = "lambda.1se.Theta" is available only when 'fit.1se' = TRUE):

## Predictions on the test set: newy = mu + newx %*% Beta.

## 's' can also be "lambda.1se.Beta" or "lambda.1se.Theta"
## (why 's' but not 'lambda'? Because we follow the naming rules of glmnet).

newy <- predict(cvfit, newx = sim.dat$X[tst, ], s = "lambda.min")

cat("dim(newy):", dim(newy))

#> dim(newy): 30 12

Users now should be able to use missoNet package with the functions introduced so far. There are many more arguments in the functions that give users a great deal of flexibility. To learn more, move on to the next section. In the section on a real data application, we demonstrate a complete process of analyzing a neuroimaging and genetic dataset using missoNet, and we highly recommend interested users to browse this content.

Other Commonly Used Function Arguments

missoNet includes a variety of functions for data simulation, goodness-of-fit evaluation, regularization parameter tuning, and visualization of results, as well as predictions in new data. In addition to the basic function arguments introduced so far, there are some other commonly used arguments available for users to customize the fit.

An Application to Structural Neuroimaging and Genetic Data

In this section, we will apply missoNet to a relatively high-dimensional dataset – bgsmtr_example_data, as a part of the bgsmtr package, obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI-1) database.

This example dataset consists of 15 structural neuroimaging measures and 486 single nucleotide polymorphisms (SNPs) determined from a sample of 632 subjects. Importantly, the 486 SNPs cover 33 genes deemed associated with Alzheimer’s disease. In this case, X is a 632-by-486 matrix containing minor allele counts (i.e. 0, 1 or 2) for 486 SNPs on the 632 subjects. Y is a 632-by-15 matrix containing volumetric and cortical thickness measures for 15 regions of interest. For more details, see Greenlaw et al. (2017) <doi:10.1093/bioinformatics/btx215> and bgsmtr.

First let’s load the dataset and randomly split the data into training and test sets. As a demonstration, we only use the SNPs with large variance (i.e. variance in the top 50% among all SNPs, p = 243) as the predictor variables. In most practical cases, this kind of unsupervised filtering of variables for dimension reduction before analyses can help the algorithm make faster progress and get more robust results.

## Load data
library(bgsmtr, quietly = TRUE)
data(bgsmtr_example_data)

## Transpose data matrix to make rows correspond to 
## samples and columns correspond to variables.
SNP <- t(bgsmtr_example_data$SNP_data)  # predictor matrix
BM <- t(bgsmtr_example_data$BrainMeasures)  # complete response matrix

## Unsupervised filtering of the top 50% variables.
SNP_var <- apply(SNP, 2, var)
SNP_subset <- SNP_var > quantile(SNP_var, 0.5)
SNP <- SNP[ ,SNP_subset]

set.seed(123)  # a random seed for reproducibility
tr <- sample(1:632, 550, replace = FALSE)  # training set indices
tst <- c(1:632)[-tr]  # test set indices

cat("dim(SNP):", dim(SNP[tr, ]), "
dim(BM):", dim(BM[tr, ]))

#> dim(SNP): 550 243
#> dim(BM): 550 15

Since this dataset does not contain missing values, we manually drop around 5% of values in the matrix of brain measures (outcomes) following the MCAR mechanism:

## Generate the indicator matrix of missingness.
M <- matrix(1, nrow(BM), ncol(BM))
NA.pos <- do.call("cbind", lapply(1:ncol(BM), function(x) {
  rbinom(nrow(BM), size = 1, prob = 0.05) == 1
}))
## All missing values should be coded as 'NA's or 'NaN's.
M[NA.pos] <- NA
BM.mis = BM * M

cat(BM.mis[1:6, 5:8])

#>     Left_InfLatVent.adj Left_LatVent.adj Left_EntCtx.adj Left_Fusiform.adj
#> V1           -473.42061       -4595.2672      0.60678817                NA
#> V2            100.12989       -9163.7594      0.29868885       0.196951758
#> V3            373.36131      -10659.8137      0.18835349       0.105863582
#> V4            876.03547               NA     -1.19832801      -0.233082858
#> V5           -887.40929       -6789.0606      0.09203347       0.098447465
#> V6            354.12157               NA     -0.32757347       0.310694508

We remark that this is a challenging task because the relatively small sample size and the existences of non-Gaussian variables make the estimation results have large variances, and heavily depend on how the data is split for the cross-validation, which is very likely to result in difficulty in convergence and abnormal cross-validation behaviors (e.g. the optimal area of λ := (λ_B, λ_Θ) is heavily skewed towards the boundary of the search domain). Therefore, we adopt a dynamic strategy of “coarse + fine” search here to avoid the unnecessary over-searching in those areas with obviously excessive cross-validated errors.

We first perform a large-scope coarse search for the optimal lambda pair (λ_Bmin^coarse, λ_Θmin^coarse) that achieves the minimum prediction error using cv.missoNet. We set 'lamBeta.scale.factor' = 3 and 'lamTheta.scale.factor' = 3 (the default is 1) to raise the upper boundaries of the search ranges for λ_B and λ_Θ (i.e. max(λ_B) and max(λ_Θ)) a little bit, meanwhile 'lamBeta.min.ratio' and 'lamTheta.min.ratio' are given a small value 1E-4 (the default) close to zero to ensure that the entire ranges of all possible solutions are covered. As a rough search, the number of lambda values ('n.lamBeta' and 'n.lamTheta') usually does not need to be very large, because we only need approximate ranges of λ_B and λ_Θ where (λ_Bmin^coarse, λ_Θmin^coarse) is most likely to exist.

## Model I (the coarse-grid fit)

## To be more in line with real-world usage, the model is trained without specified 
## missing probabilities for the response variables ('rho' = NULL) and lambda sequences 
## for the lasso penalties ('lambda.Beta' = NULL, 'lambda.Theta' = NULL), in which case 
## the program will automatically compute and use reasonable values.

cl <- parallel::makeCluster(min(parallel::detectCores()-1, 2))
cvfit.BM <- cv.missoNet(X = SNP[tr, ], Y = BM.mis[tr, ], kfold = 5, 
                        rho = NULL, lambda.Beta = NULL, lambda.Theta = NULL,
                        lamBeta.min.ratio = 1e-4, lamTheta.min.ratio = 1e-4,
                        n.lamBeta = 20, n.lamTheta = 20,
                        lamBeta.scale.factor = 3, lamTheta.scale.factor = 3,
                        standardize = TRUE, standardize.response = TRUE,
                        permute = TRUE, with.seed = 433,
                        parallel = TRUE, cl = cl, verbose = 1)
parallel::stopCluster(cl)

Then we plot the heatmap of the standardized mean cross-validated errors:

plot(cvfit.BM)

Now we can pick more reasonable ranges for λ_B and λ_Θ for a further fine-grid search according to the figure above. We observed that the current optimal λ = (λ_Bmin^coarse, λ_Θmin^coarse) achieving the minimum cross-validated error (white box) is about (0.24, 0.62), while max(λ_B) and max(λ_Θ) are around 2.66 (appropriate) and 29.94 (too large), respectively. Therefore, we decide to shrink 'lamTheta.scale.factor' from 3 to 0.3 (= 3/10) as we believe 2.99 (= 29.94/10) is a more appropriate upper boundary for λ_Θ. 'lamBeta.scale.factor' can remain at 3 so that max(λ_B) will still be 2.66. For the same reason, the new min(λ_B) and min(λ_Θ) are determined to be 0.027 (= 2.66 * 0.01) and 0.030 (= 2.99 * 0.01) respectively by setting 'lamBeta.min.ratio' = 0.01 and 'lamTheta.min.ratio' = 0.01.

TIPS: Sometimes the lambda sequences automatically computed by the program may have inappropriate ranges, if the optimal lambda pair "lambda.min" := (λ_Bmin, λ_Θmin) selected by the cross-validation is located at the boundary of the search domain, users can use the arguments 'x.scale.factor' and/or 'x.min.ratio' to shift and/or zoom the search ranges accordingly.

Because we are going to do a fine-grid search, the number of lambda values can be slightly increased (here we increase 'n.lamBeta' and 'n.lamTheta' from 20 to 40). Typically we suggest 'n.x' = -log10('x.min.ratio') * c, where c ∈ [10, 20] (x is a placeholder for lamBeta or lamTheta), to achieve a balance between the computational fidelity and speed.

## Model II (the fine-grid fit)

## For brevity, arguments taking the defaults are not specified explicitly.
cl <- parallel::makeCluster(min(parallel::detectCores()-1, 2))
cvfit2.BM <- cv.missoNet(X = SNP[tr, ], Y = BM.mis[tr, ], kfold = 5,
                         lamBeta.min.ratio = 0.01, lamTheta.min.ratio = 0.01,
                         n.lamBeta = 40, n.lamTheta = 40,
                         lamBeta.scale.factor = 3, lamTheta.scale.factor = 0.3,
                         fit.relax = TRUE, permute = TRUE, with.seed = 433,
                         parallel = TRUE, cl = cl, verbose = 1)
parallel::stopCluster(cl)

## Certainly, users can specify the lambda values (linear on the log scale) by themselves
## for the cross-validation search rather than use the control arguments. The following 
## commands should return almost the same results as the above ones (subtle differences 
## come from rounding the float numbers).

# cl <- parallel::makeCluster(min(parallel::detectCores()-1, 2))
# cvfit2.BM <- cv.missoNet(X = SNP[tr, ], Y = BM.mis[tr, ], kfold = 5,
#                          lambda.Beta = 10^(seq(from = log10(2.66), to = log10(2.66*0.01), length.out = 40)),
#                          lambda.Theta = 10^(seq(from = log10(2.99), to = log10(2.99*0.01), length.out = 40)),
#                          fit.relax = TRUE, permute = TRUE, with.seed = 433,
#                          parallel = TRUE, cl = cl, verbose = 0)
# parallel::stopCluster(cl)

Plot the standardized mean cross-validated errors of this refitted model, 'detailed.axes' = FALSE can prevent axes to get squeezed together:

cat("\"lambda.min\":
  - lambda.Beta:", cvfit2.BM$est.min$lambda.Beta, "
  - lambda.Theta:", cvfit2.BM$est.min$lambda.Theta)

plot(cvfit2.BM, detailed.axes = FALSE)

#> "lambda.min":
#>   - lambda.Beta: 0.198283
#>   - lambda.Theta: 0.573265

We found the new "lambda.min" = (λ_Bmin^fine, λ_Θmin^fine) is (0.20, 0.57), which is actually very close to the one from the coarse search (λ_Bmin^coarse, λ_Θmin^coarse) = (0.24, 0.62), showing that one single search is able to provide reasonable accuracy most of the time as long as sufficiently large search ranges are given.

As we did before, the estimated parameters $\hat{\mathbf{B}}$ and $\hat{\mathbf{\Theta}}$ at the new "lambda.min" can be visualized by heatmaps. For the sake of brevity, we shall skip the discussion of "lambda.1se.Beta" and "lambda.1se.Theta" here.

## Theta_hat at "lambda.min"
Theta_hat <- cvfit2.BM$est.min$Theta
rownames(Theta_hat) <- colnames(Theta_hat) <- colnames(BM)
col <- circlize::colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
ht1 <- plot_heatmap(Theta_hat, col = col, legend = TRUE, lgd_name = "values",
                    title = expression(paste(hat(bold(Theta)), " at `lambda.min` (unclustered)")))
plot(ht1)

## Beta_hat at "lambda.min"
Beta_hat <- cvfit2.BM$est.min$Beta
colnames(Beta_hat) <- colnames(BM)
col <- circlize::colorRamp2(c(-0.05, 0, 0.05), c("blue", "white", "red"))
ht2 <- plot_heatmap(Beta_hat, col = col, legend = TRUE, lgd_name = "values",
                    title = expression(paste(hat(bold(Beta)), " at `lambda.min` (unclustered)")))
plot(ht2)

## Relaxed network
relax.net <- cvfit2.BM$est.min$relax.net
rownames(relax.net) <- colnames(relax.net) <- colnames(BM)
col <- circlize::colorRamp2(c(-4, 0, 4), c("blue", "white", "red"))
ht3 <- plot_heatmap(relax.net, col = col, legend = TRUE, lgd_name = "values",
                    title = "Relaxed fit of network (unclustered)")
plot(ht3)

## Plot the partial correlation matrix of the 15 neuroimaging measures.
## Must supply a dataset without missing values, so we use the original dataset here.
pcor.mat <- ppcor::pcor(BM)  # not BM.mis
corrplot::corrplot(pcor.mat$estimate, p.mat = pcor.mat$p.value, 
                   sig.level = 0.01, insig = "blank", 
                   tl.cex = 0.5, tl.col = "black", mar = c(.5, .1, 2, .1),
                   title = "Partial correlation matrix (unclustered)")

The conditional network structure $\hat{\mathbf{\Theta}}$ can be represented by an undirected graph 𝒢, consisting of a set of vertices V and a set of edges E, where an edge connects a pair of variables (Y_j and Y_k) if and only if they are conditionally dependent ($\hat{\mathbf{\Theta}}_{jk} \neq 0$). One important feature of missoNet is that the program supports fitting a relaxed conditional graphical model upon response variables (by setting 'fit.relax' = TRUE). The active edges (non-zero off-diagonal elements) of $\hat{\mathbf{\Theta}}$ will be re-estimated without penalization ('lambda.Theta' = 0), such a debiased estimate of the network structure would be useful for exploring the conditional interdependencies among traits of interest.

The two figures in the second row above show the relaxed conditional network and the partial correlation matrix of the 15 neuroimaging measures. As will be readily seen that missoNet sharply simplifies dependencies compared with the naive partial correlation analysis below right (insignificant correlations with p > 0.01 were not plotted), by regressing on SNPs and adjusting for the genetic effects, which helps us better focus on those truly important, and non-genetic induced connections in brain. Moreover, we have to emphasize that unlike the traditional partial correlation analysis which requires complete datasets (although pair-wise complete cases can be used, the results tend to be biased), missoNet substantially reduces the difficulty of data acquisition and processing by allowing for missing values in the tasks.

If we cluster the estimated coefficient matrix $\hat{\mathbf{B}}$ by rows (SNPs) in an unsupervised manner:

colnames(Beta_hat) <- colnames(BM)  # brain measure names
rownames(Beta_hat) <- colnames(SNP)  # SNP names
Gene_groups <- bgsmtr_example_data$SNP_groups[SNP_subset]  # gene names

ComplexHeatmap::pheatmap(Beta_hat,
                         annotation_row = data.frame("Gene" = Gene_groups, row.names = rownames(Beta_hat)),
                         color = circlize::colorRamp2(c(-0.05,0,0.05), c("blue","white","red")), show_rownames = FALSE,
                         border_color = NA, border = TRUE, name = "values", cluster_rows = TRUE, cluster_cols = FALSE,
                         column_title = expression(paste(hat(bold(Beta)), " at `lambda.min` (clustered)")))

The figure above depicts a very large proportion of coefficients with zero values, and there is an obvious clustering of the contributing SNPs with respect to the genes (see row-wise annotation), implying potential pathways between certain SNPs and changes in the local structures of brain.

Finally, we end this section by refitting the model using the original dataset without missing values:

## Model III (the complete-data fit)

## All arguments are kept unchanged compared to model II, except BM.mis -> BM.
cl <- parallel::makeCluster(min(parallel::detectCores()-1, 2))
cvfit3.BM <- cv.missoNet(X = SNP[tr, ], Y = BM[tr, ], kfold = 5,
                         lamBeta.min.ratio = 0.01, lamTheta.min.ratio = 0.01,
                         n.lamBeta = 40, n.lamTheta = 40,
                         lamBeta.scale.factor = 3, lamTheta.scale.factor = 0.3,
                         fit.relax = TRUE, permute = TRUE, with.seed = 433,
                         parallel = TRUE, cl = cl, verbose = 1)
parallel::stopCluster(cl)

Then we can assess the prediction performance of the models trained using the corrupted vs. complete dataset.

## Predictions on the test set.
newy.model2 <- predict(cvfit2.BM, newx = SNP[tst, ], s = "lambda.min")
newy.model3 <- predict(cvfit3.BM, newx = SNP[tst, ], s = "lambda.min")

cat("MAE on the test set:
  - Model II (corrupted):", mean(abs(BM[tst, ] - newy.model2)), "
  - Model III (complete):", mean(abs(BM[tst, ] - newy.model3)))

#> MAE on the test set:
#>   - Model II (corrupted): 2543.588
#>   - Model III (complete): 2543.591

As expected, there is no obvious difference between the two models with regard to the mean absolute prediction error (MAE) on the test set, since we generated the missing data simply following the MCAR mechanism. Real missing data usually falls on a continuum between one extreme – where the systematic missingness pattern depends entirely on the observed data (pure MAR), and the other extreme – where the systematic missingness pattern depends entirely on the missing data (pure MNAR). Even though, missoNet can still provide good enough estimates compared with those algorithms that cannot tolerate missing values.