structuredEstimate: Estimation of Partial Correlation Matrix Given Zero...

In GGMridge: Gaussian Graphical Models Using Ridge Penalty Followed by Thresholding and Reestimation

Estimation of Partial Correlation Matrix Given Zero Structure.

Description

Estimation of nonzero entries of the partial correlation matrix given zero structure.

Usage

structuredEstimate(x, E)

Arguments

x	n by p data matrix with the number of variables p and sample size n.
E	The row and column indices of zero entries of the partial correlation matrix.

Value

A list containing

R	The partial correlation matrix.
K	The inverse covariance matrix.
RSS	The residual sum of squares.

Author(s)

Min Jin Ha

References

Ha, M. J. and Sun, W. (2014). Partial correlation matrix estimation using ridge penalty followed by thresholding and re-estimation. Biometrics, 70, 762–770.

Examples

 p <- 100 # number of variables
 n <- 50 # sample size
 
 ###############################
 # Simulate data
 ###############################
 simulation <- simulateData(G = p, etaA = 0.02, n = n, r = 1)
 data <- simulation$data[[1L]]
 stddata <- scale(x = data, center = TRUE, scale = TRUE)
 
 ###############################
 # estimate ridge parameter
 ###############################
 lambda.array <- seq(from = 0.1, to = 20, by = 0.1) * (n-1.0)
 fit <- lambda.cv(x = stddata, lambda = lambda.array, fold = 10L)
 lambda <- fit$lambda[which.min(fit$spe)]/(n-1)
 
 ###############################
 # calculate partial correlation 
 # using ridge inverse
 ###############################
 w.upper <- which(upper.tri(diag(p)))
 
 partial <- solve(lambda * diag(p) + cor(data))
 partial <- (-scaledMat(x = partial))[w.upper]
 
 ###############################
 # get p-values from empirical 
 # null distribution 
 ###############################
 efron.fit <- getEfronp(z = transFisher(x = partial), 
                        bins = 50L, 
                        maxQ = 13)
 
 ###############################
 # estimate the edge set of 
 # partial correlation graph with 
 # FDR control at level 0.01
 ###############################
 w.array <- which(upper.tri(diag(p)),arr.ind=TRUE)
 th <- 0.01
 wsig <- which(p.adjust(efron.fit$correctp, method="BH") < th )
 E <- w.array[wsig,]
 dim(E)
 
 ###############################
 # structured estimation
 ###############################
 fit <- structuredEstimate(x = stddata, E = E)
 th.partial <- fit$R

GGMridge documentation built on Nov. 25, 2023, 1:08 a.m.