wfrl: joint estimation of multiple regression coefficient matrices In ldstatsHD: Linear Dependence Statistics for High-Dimensional Data

Description

wfrl estimates jointly two regression coefficient matrices from multivariate normal distributed datasets using an ADMM based algorithm.

Usage

 1 2 3 wfrl(D1, D2, lambda1, lambda2, automLambdas = TRUE, paired = TRUE, sigmaEstimate = "CRmad", maxiter=30, tol=1e-05, nsubset = 10000, rho = 1, rho.increment = 1, notOnlyLambda2 = TRUE) 

Arguments

 D1 list with the response variables. Two matrices in the list corresponding to the response variables of the two populations. D2 list with the explanatory variables. Two matrices in the list corresponding to the explanatory variables of the two populations. lambda1 tuning parameter for sparsity in the regression coefficients. lambda2 tuning parameter for similarity between the regression coefficients in the two populations. automLambdas if TRUE the lambda's are estimated automatically with lambda1 and lambda2 being expected false positive rate levels. paired if TRUE, observations in D1 and D2 are assumed to be matched (n_1 must be equal to n_2). sigmaEstimate robust method used to estimate the variance of estimated partial correlations: name that uniquely identifies "mad", "IQR" or "CRmad" (default). This measure is used to automatically select the tuning parameter (when automLambdas = TRUE). maxiter maximum number of iterations for the ADMM algorithm. tol convergence tolerance. nsubset maximum number of estimated partial correlation coefficients (chosen randomly) used to select lambda1 and lambda2 automatically. rho regularization parameter used to compute matrix inverse by eigen value decomposition (default of 1). rho.increment default of 1. notOnlyLambda2 if FALSE only lambda2 is found automatically.

Details

wfrl uses a weighted-fused least squares lasso maximum likelihood estimator by solving:

[\hat{β}_H,\hat{β}_T] = \arg\min\limits_{β_H,β_T} ≤ft[ \frac{1}{2n}||Y-β_HX||^2_2 + \frac{1}{2n}||Q-β_TW||^2_2 +P_{λ_1,λ_2,V}(β)\right]

with

P_{λ_1,λ_2,V}(β) = λ_1||β_H||_1 + λ_1||β_T||_1 + λ_2||V \circ (β_T-β_H)||_1.

where λ_1 is the sparsity tuning parameter, λ_2 is the similarity tuning parameter, and V = [v_{ij}] is a p\times p matrix to weight λ_2 for each coefficient of the differential precision matrix. If datasets are independent (paired = "FALSE"), then it is assumed that v_{ij} = 1 for all pairs (i,j). Otherwise (paired = "TRUE"), weights are estimated in order to account for the dependence structure between datasets in the differential network estimation. An ADMM-type recursive algorithm is used to solve the optimization problem.

See details in wfgl for transforming the selection problem of the tuning parameters λ_1 and λ_2.

Value

An object of class wfrl containing the following components:

 regCoef  regression coefficients. path  non-zero structure for the regression coefficients. diff_value  convergence control. iters  number of iterations used.

Author(s)

Caballe, Adria <[email protected]>, Natalia Bochkina and Claus Mayer.

References

Danaher, P., P. Wang, and D. Witten (2014). The joint graphical lasso for inverse covariance estimation across multiple classes. Journal of the Royal Statistical Society: Series B (Statistical Methodology) (2006), 1-20.

Boyd, S. (2010). Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends in Machine Learning 3(1), 1-122.

plot.wfrl for graphical representation.
wfgl for joint partial correlation estimation.
  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  # example to use of wfrl N <- 200 EX2 <- pcorSimulatorJoint(nobs = N, nclusters = 3, nnodesxcluster = c(30, 30,30), pattern = "pow", diffType = "cluster", dataDepend = "diag", low.strength = 0.5, sup.strength = 0.9, pdiff = 0.5, nhubs = 5, degree.hubs = 20, nOtherEdges = 30, alpha = 2.3, plus = 0, prob = 0.05, perturb.clust = 0.2, mu = 0, diagCCtype = "dicot", diagNZ.strength = 0.6, mixProb = 0.5, probSign = 0.7, exactZeroTh = 0.05) P <- EX2$P q <- 50 BETA1 <- array(0, dim = c(P, q)) diag(BETA1) <- rep(0.35,q) BETA2 <- BETA1 diag(BETA2)[c(1:floor(q/2))] <- 0 sigma2 <- 1.3 Q <- scale(EX2$D1) W <- scale(EX2\$D2) X <- Q%*%BETA1 + mvrnorm(N,rep(0,q),diag(rep(sigma2,q))) Y <- W%*%BETA2 + mvrnorm(N,rep(0,q),diag(rep(sigma2,q))) D1 <- list(scale(X), scale(Y)) D2 <- list(scale(Q), scale(W)) ## not run #wfrl1 <- wfrl(D1, D2, lambda1 = 0.05, lambda2 = 0.05, automLambdas = TRUE, paired = FALSE, # sigmaEstimate = "CRmad", maxiter = 30, tol = 1e-05, nsubset = 10000, rho = 1, # rho.increment = 1, notOnlyLambda2 = TRUE) #print(wfrl1)