In LuZhangH/floodgate: Floodgate: inference for model-free variable importance

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Setup

set.seed(1995)
## load required packages
library(floodgate)
library(methods)
library(conformalInference)
library(glmnet)
library(lars)
library(randomForest)
library(SAM)
library(ggplot2)
## load utility functions: related to model fitting
source("../utils/algo_utils.R")

#### problem setup
n = 1000 # sample size
p = 1000 # covariate dimension

Xmodel = "gaussian" # covariate distribution
rho = 0.3 # auto-correlation coefficient

Ydist = "gaussian" # conditional model of response
s = 20 # number of non-nulls
amplitude = 5 # signal amplitude value

split.prop = 0.5 # splitting proportion
K = 500 # number of null replicates
alevel = 0.05 # confidence level

Prepare data

## load model parameters of the covariate distribution
load(paste0("../inst/rho", rho, "_Sigma.RData"))
load(paste0("../inst/rho", rho, "_X_paras_gaussian.RData"))

## choose non-null varaibles randomly 
S_star = sort(sample(1:p,s))
beta = rep(0,p)
beta[S_star] = sample(c(-1,1), s, replace = TRUE) * amplitude/sqrt(n)

## generate the covaraites X 
X = matrix(rnorm(n*p),n,p)%*% Sigma.chol
## Generate the response Y from a linear model
Y = gen.Y(X, beta, Ydist = Ydist)

Sample null covariates and compute variable importance measures

## sample null covariates 
nulls.list = sample.gaussian.nulls(X = X, S = as.list(1:p), K = K, gamma_X.list_S = gamma_X.list, 
                                   sigma_X.list_S = sigma_X.list)

## compute the mMSEgap values
mMSEgap = compute.movi(beta = beta, Xmodel = Xmodel, Ydist = Ydist,
                       sigma_X.list = sigma_X.list, nulls.list = nulls.list)

Run floodgate

## sample splitting
i1 = sample(1:n, floor(n*split.prop))
i2 = (1:n)[-i1]
n1 = length(i1)
n2 = length(i2)

## use LASSO to estimate the conditional mean
algo = "lasso"
funs = funs.list[[algo]]
## run floodgate to obtain LCBs
fg.out = floodgate(X, Y, i1, i2, nulls.list = nulls.list,
                   gamma_X.list = gamma_X.list, sigma_X.list = sigma_X.list,
                   Xmodel = Xmodel, funs, algo = algo,
                   alevel = alevel, verbose = TRUE)

Plot the results

The following plot shows the floodgate lower confidence bound (LCB): the horizontal bar with a black color and the mMSEgap: the star-shaped point with a red color. The step of model fitting on the training data also outputs a selected subset $S$, in additional to a regression function estimator $\mu$. We only plot the LCBs for covariates in $S$ here.

## extract output and produce plots
inf.out = as.data.frame(fg.out$inf.out)
S = unlist(fg.out$S)
inf.out$mMSEgap = mMSEgap[S]
ggplot(data = inf.out, aes(x = S, y = mMSEgap)) + 
      ylim(0, max(inf.out$mMSEgap, inf.out$LCB) + 0.05) + 
      ggtitle(paste0("algo = ", algo )) +
      ylab("mMSEgap and LCB") + xlab("Selected variables") + 
      geom_point(color = "red", shape = 8, size = 3.5) + 
      geom_errorbar(aes(ymin=LCB, ymax=LCB), width = 15, color = "black") +
      geom_segment(aes(x = S, y = LCB, xend = S, yend = mMSEgap),
                    arrow = arrow(length = unit(0.15, "cm"), type = "closed"))

The arrow for a given covaraite starts from the LCB and ends at the mMSEgap. The arrow provides a good illustration of LCB's performance, with its length being the half-width and its direction indicating coverage/miscoverage (upward: coverage; downward: miscovergae; leftward: coverage with $\text{LCB}=\mathcal{I} =0$).

Try with a different fitting algorithm

## use SAM to estimate the conditional mean
algo = "sam"
funs = funs.list[[algo]]
## run floodgate to obtain LCBs
fg.out = floodgate(X, Y, i1, i2, nulls.list = nulls.list,
                   gamma_X.list = gamma_X.list, sigma_X.list = sigma_X.list,
                   Xmodel = Xmodel, funs, algo = algo,
                   alevel = alevel, verbose = TRUE)
## extract output and produce plots
inf.out = as.data.frame(fg.out$inf.out)
S = unlist(fg.out$S)
inf.out$mMSEgap = mMSEgap[S]
ggplot(data = inf.out, aes(x = S, y = mMSEgap)) + 
      ylim(0, max(inf.out$mMSEgap, inf.out$LCB) + 0.05) + 
      ggtitle(paste0("algo = ", algo )) +
      ylab("mMSEgap and LCB") + xlab("Selected variables") + 
      geom_point(color = "red", shape = 8, size = 3.5) + 
      geom_errorbar(aes(ymin=LCB, ymax=LCB), width = 15, color = "black") +
      geom_segment(aes(x = S, y = LCB, xend = S, yend = mMSEgap),
                    arrow = arrow(length = unit(0.15, "cm"), type = "closed"))