High-Dimensional Gaussian Graphical Models

This page shows an example of how to generate conditional knockoffs for a multivariate gaussian graphical model as in Section 3.2 of the accompanying paper, which allows the dimension of covariates to be much larger than the sample size.

library(knockoff)
source('src/knockoff_measure.R')
source('src/util.R')
source('src/cknock_ldg.R')
source('src/cknock_ggm.R')

In this simulation we assume that the distribution of the covariates belongs to a Gaussian graphical model and a superset of its edges is known. Examples include autoregressive models and spatial-lattice models. Here we use AR(10) as an example.

First, we set up the distributional parameters. Specifically, the covariance matrix is obtained by renormalizing a banded precision matrix \(\mathbf{\Omega}_0\), whose entries satisfy \[-\mathbf{\Omega}_0 (i,j) / \sqrt{\mathbf{\Omega}_0 (i,i) \mathbf{\Omega}_0 (j,j) }= \left\{ \begin{array}{c l} 0.05 \quad &\text{, if } |i-j|\le 10, \\ 0 \quad &\text{, if } |i-j|>10. \end{array} \right.\]

set.seed(2019)
# Problem parameters
p = 2000    # number of covariates
n = 400     # number of observations
k = 60      # number of non-zero regression coefficients
A = 60      # signal amptitude 
nonzero = sample(p, k) 
beta = 1 * (1:p %in% nonzero) * sample(c(-1,1),p,rep = T) 


# setup the banded precison matrix
nn=10            # bandwidth
rho0=-0.05       # the non-zero partial correlation
Omega0 = (1 - rho0) * diag(p) + rho0
Omega0[abs(outer(1:p, 1:p, '-')) > nn] = 0
S=solve(Omega0)
# renormalize so that the marginal variance = 1
ds=diag(S)
Omega=sqrt(ds)%diag*%Omega0%*diag%sqrt(ds)
Sigma=sqrt(1/ds)%diag*%S%*diag%sqrt(1/ds)
Sigma.chol = chol(Sigma)

Then we generate the covariates and response.

# Generate the covariate 
X = matrix(rnorm(n * p), n) %*% Sigma.chol
# Known graphical information: adjacent matrix
Adj = (Omega!=0 )

# Generate the response from a linear model
Y = X%*%beta*A/sqrt(n) + rnorm(n)

The conditional knockoffs can be generated by blocking and two-fold data spliting, where the input threshold \(n'=50\), the number of folds \(K=2\) and the method for computing \(\mathbf{s}\) is mix (as discussed in Appendix B.1.1 of the paper). This function runs Algorithm 12 in Appendix B.2.2 to search for blocking sets and then runs Algorithm 4 to generate knockoffs for each fold.

Xk.cond = cknockoff.ggm(X,Adj,threshold=50,K=2,method='mix')

Once we have generated conditional knockoffs, all the usual knockoffs machinery applies.

fdr.level = 0.2
knockoff.stat.fun =  mod.stat.lasso_coefdiff
filter.cond = knockoff.filter(X,Y,
    knockoffs = function(x) {
     Xk.cond
    },
    statistic = knockoff.stat.fun,
    fdr = fdr.level,offset = 1)
# false positive and false negative
c(fp(filter.cond$selected, beta),
 fn(filter.cond$selected, beta))

## [1] 0.1800000 0.6833333