Conditional independence testing

This page provides a step-by-step guide on how to conduct the Graphical conditional randomization test (G-CRT) ( Algorithm 4 in the accompanying paper). Our methodology is based on a combination of the conditional randomization test CRT and our algorithm for generating exchangeable copies in Section 2.

library(glmnet)
library(rags2ridges)
library(glasso)
library(Matrix)
library(MASS)
library(mvtnorm)
library(reshape)
library(ggplot2)
source("CRT-Simulate.R")
source("CRT-Stat.R")
source("DBL_GLMs_functions.R")
source("OneExp.R")
source("Run-Testing.R")
source("dCRT_function.R")
source("GoF-Settings.R")
source("CSS-GGM.R")
source("GGMTesting.R")
source("CIT_function_for_webpage.R")

Looking into the CRT function

Take the Gaussian linear regression for example

setting = "l_l_w"             # gaussian linear regression (low-dimensional setting)
n = c(50)             
p = c(20)
setT=1:8

K = 6
s = 0.2
grid_signal= seq(0,2.5,by=0.25)   # control signal strength

# get graph G
true_graph=graph_band(p,K,s)
perm.var=sample(p,p)
sigma=true_graph$sigma[perm.var,perm.var]
inv.sigma=true_graph$inv.sigma[perm.var,perm.var]
diag_s=sqrt(diag(sigma))
sigma=sweep(sweep(sigma,1,1/diag_s,"*"),2,1/diag_s, "*")
inv.sigma=sweep(sweep(inv.sigma,1,diag_s,"*"),2,diag_s, "*")
G = prec_to_adj(inv.sigma)       # graph G

# Simulate the data X
X <- simulate_data(sigma, n)
beta=random_beta(p)              # regression coefficient
noise_unif=runif(n)

Sampling M exchangeable copies of X, L is the number of iterations in each Markov chain (Algorithm 2)

CSS_sample = exchangeable_sampling(X, G, I = setT, bReduced = T, M=100, L=3)
X_copies  = CSS_sample$X_copies

We choose statistics LM-SST and LM-SSR by comparison to the classical F-test for the null hypothesis \(H_0: \beta_{\mathcal{T}}=0\) and the dCRT using a Gaussian Lasso model (Bonferroni adjusted).

# statistics
List_Stat_Comb=c("Stat_Comb_Linear_SST","Stat_Comb_Linear_Dev")

CRT test (Algorithm 4)

test_result = GRT_test_webpage(X, CSS_sample, setting, setT, beta, niose_unif, grid_signal, List_Stat_Comb,IncludeBaseline = T)

Output the p-value

for(ind.signal in 1:length(grid_signal)){
  print(paste("theta=",test_result[[ind.signal]]$theta))
  output_pvalue = test_result[[ind.signal]]$pvalues[-which(names(test_result[[ind.signal]]$pvalues)=="dCRT_RF")]
  names(output_pvalue) = c("LM-SST","LM-SSR","F-test","dCRT")
  print(output_pvalue)
}
## [1] "theta= 0"
##    LM-SST    LM-SSR    F-test      dCRT 
## 0.9900990 0.9900990 0.9833458 0.9944809 
## [1] "theta= 0.25"
##    LM-SST    LM-SSR    F-test      dCRT 
## 0.9108911 0.9405941 0.9285549 1.0000000 
## [1] "theta= 0.5"
##    LM-SST    LM-SSR    F-test      dCRT 
## 0.6534653 0.6534653 0.6953275 1.0000000 
## [1] "theta= 0.75"
##    LM-SST    LM-SSR    F-test      dCRT 
## 0.2574257 0.3465347 0.3362194 0.5788404 
## [1] "theta= 1"
##     LM-SST     LM-SSR     F-test       dCRT 
## 0.05940594 0.12871287 0.10276437 0.35787418 
## [1] "theta= 1.25"
##     LM-SST     LM-SSR     F-test       dCRT 
## 0.02970297 0.02970297 0.02213120 0.11063836 
## [1] "theta= 1.5"
##      LM-SST      LM-SSR      F-test        dCRT 
## 0.019801980 0.029702970 0.003815606 0.038739981 
## [1] "theta= 1.75"
##       LM-SST       LM-SSR       F-test         dCRT 
## 0.0099009901 0.0198019802 0.0005840345 0.0160886540 
## [1] "theta= 2"
##       LM-SST       LM-SSR       F-test         dCRT 
## 9.900990e-03 9.900990e-03 8.538935e-05 2.388902e-03 
## [1] "theta= 2.25"
##       LM-SST       LM-SSR       F-test         dCRT 
## 9.900990e-03 9.900990e-03 1.251375e-05 2.313675e-03 
## [1] "theta= 2.5"
##       LM-SST       LM-SSR       F-test         dCRT 
## 9.900990e-03 9.900990e-03 1.894425e-06 6.591402e-04

Take theta = 0, statistics: LM-SST for example, show \(T(Y, X), T(Y, \tilde{X}^{(1)}), \cdots, T(Y, \tilde{X}^{(M)})\)

hist(test_result$`theta= 0`$Ttildes$Stat_Comb_Linear_SST,xlab = "LM-SST_Ttilde",main="Histogram of LM-SST_Ttilde")
abline(v=test_result$`theta= 0`$T0s$Stat_Comb_Linear_SST,col="red")
text(test_result$`theta= 0`$T0s$Stat_Comb_Linear_SST, 0,labels = "T0", pos = 3, col = "red")