joint estimation of multiple regression coefficient matrices

Share:

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 <a.caballe@sms.ed.ac.uk>, 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.

See Also

plot.wfrl for graphical representation.
wfgl for joint partial correlation estimation.

Examples

 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)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.