通过Model-X Knockoffs控制variable Selection

翻译于https://web.stanford.edu/group/candes/knockoffs/software/knockoffs/tutorial-1-r.html

使用R中的knockoffs包

这个小插图说明了使用Model-X knockoffs的knockoff包的基础使用方法。在这种情况下，我们假设预测变量分布是已知的（或者它可以很好地近似），但我们不对响应的条件分布做出假设。简单起见，我们将使用通过线性模型构建的合成数据，使得响应仅取决于有小部分变量。

set.seed(1234)

# Problem parameters
n = 1000				# number of observations
p = 1000				# number of variables
k = 60					# number of variables with nonzero coefficients
amplitude = 4.5				# signal amplitude (for noise level = 1)

# Generate the variables from a multivariate normal distribution
mu = rep(0, p)
rho = 0.25
Sigma = toeplitz(rho^(0:(p-1)))
X = matrix(rnorm(n*p),n) %*% chol(Sigma)

# Generate the response from a linear model
nonzero = sample(p,k)
beta = amplitude * (1:p %in% nonzero) / sqrt(n)
y.sample = function(X) X %*% beta + rnorm(n)
y = y.sample(X)

第一个例子

首先，我们使用所有默认设置调用knockoff.filter。

library(knockoff)
result = knockoff.filter(X,y)

我们可以显示结果通过：

print(result)

## Call:
## knockoff.filter(X = x,y = y)
##
## Selected variables:
## [1]   3   9  40  44  46  61  67  78  85 108 148 153 172 173 177 210 223
## [18] 238 248 281 295 301 302 317 319 326 334 343 360 364 378 384 389 421
## [35] 426 428 451 494 506 510 528 534 557 559 595 596 617 668 676 682 708
## [52] 770 775 787 836 844 875 893 906 913 931 937 953 959

目标错误发现率的默认值为0.1.在这个实验中，错误发现比例为：

fdp = function(selected) sum(beta[selected] == 0) / max(1,length(selected))
fdp(result$selected)

## [1] 0.171875

默认情况下，knockoff filter创建二阶高斯model-X knockoffs。这种结构从数据中预估$X$行的均值$mu$和协方差$sigma$，而不是使用真实参数抽取变量。
knockoff包还包括其他knockoff构造，所有这些方法的名称都以knockoff.create为前缀。在下一个片段中，我们使用真实模型参数生成knodkoffs。

gaussian_knockoffs = function(X) create.gaussian(X, mu, Sigma)
result = knockoff.filter(X, y, knockoffs=guassian_knockoffs)
print(result)

## Call:
## knockoff.filter(X = X, y = y, knockoffs = gaussian_knockoffs)
## 
## Selected variables:
##  [1]   3   9  40  44  46  61  85 108 148 153 172 173 177 210 223 238 248
## [18] 295 301 319 326 334 343 360 364 378 384 389 421 426 428 451 506 557
## [35] 559 595 668 708 770 775 844 893 906 913 931 937 953 959

现在错误发现比例为：

fdp = function(selected) sum(beta[selectecd] == 0) / max(1, length(selected))
fdp(result$selected)

## [1] 0.0625

默认情况下，knockoff filter使用基于lasso的测试统计。具体来说，它使用stat.glmnet_coefdiff统计，它计$W_j=|Z_j|-|\widehat{Z_j}|$，其中$Z_j$和$\widehat{Z_j}$分别是第j个变量及其knockoff的lasso系数估计值。正则化参数lambda的值通过交叉验证并使用glmnet计算。
还有其他几个内置统计函数，所有统计数据的名称都以sta为前缀。例如，我们可以使用基于随机森林的统计数据。除了选择不同的统计数据外，我们还可以改变目标FDR级别（例如我们将它增加到0.2）。

result = knockoff.filter(X, y = y, knockoffs = gaussian_knockoffs,statistic = stat.random_forest,fdr=0.2)
print(result)

## Call:
## knockoff.filter(X = X, y = y, knockoffs = gaussian_knockoffs, 
##     statistic = stat.random_forest, fdr = 0.2)
## 
## Selected variables:
##  [1]   9  85 108 148 172 173 210 223 238 248 301 334 343 384 426 428 557
## [18] 595 668 708 770 785 906 931 953

fdp(result$selected)

## [1] 0.08

用户定义的测试统计信息

除了使用预定义的测试统计信息之外，还可以使用您自己的自定义测试统计信息。为了说明这个功能，我们实现了原始knockoff filter论文中最简单的测试统计之一，即$W_j=X_j\ast y-\widehat{X_j}\ast y$。

my_knockoff_stat = function(X, X_k, y){
	abs(t(X) %*% y) - abs(t(X_k) %*% y)
}
result = knockoff.filter(X, y, knockoffs = gaussian_knockoffs, statistic = my_knockoff_stat)
print(result)

## Call:
## knockoff.filter(X = X, y = y, knockoffs = gaussian_knockoffs, 
##     statistic = my_knockoff_stat)
## 
## Selected variables:
##  [1] 108 148 173 223 238 248 274 301 421 426 668 708 906 931 937 953 959

fdp(result$selected)

## [1] 0.2394366

参数nlambda由stat.glmnet_coefdiff传递给glment，gkmnet用于计算lasso path。有关此参数和其他参数的详细信息，请参阅stat.glmnet_coefdiff或glmnet.glmnet的文档。

用户定义的knockoff生成函数

除了使用构造knockoff变量的预定义程序之外，还可以创建自己的knockoffs。为了说明这个功能，我们实现了一个简单的包装器来构造二阶Model-X knockoffs。

create_knockoffs = function(X){
	create.second_order(X, shrink=T)
}
result = knockoff.filter(X, y, knockoffs=create_knockoffs)
print(result)

## Call:
## knockoff.filter(X = X, y = y, knockoffs = create_knockoffs)
## 
## Selected variables:
##  [1]   9  40  61  85 108 148 153 172 173 177 210 223 248 295 301 319 326
## [18] 334 343 360 364 378 384 389 421 426 428 451 506 559 595 668 708 770
## [35] 844 893 906 913 931 937 953 959

fdp(result$selected)

## [1] 0

近似与完全SDP knockoffs

knockoff包支持knockoff变量的两种主要类型，半定规划（SDP）knockoffs（默认）和等相关knockoffs。虽然计算成本很高，但是SDPknockoffs在统计上句有更高的功率。为了创建SDP knockoffs，这个包依赖于R库Rdsdp来有效的解决半定规划。在高维设置中，该程序在计算上变得难以处理。然后近似SDP提供了解决方案，解决这个问题通过解决一个基于协方差矩阵的块对角矩阵的更简单的轻松问题。默认情况下，如果p < 500，knockoff filter则使用SDP knockoffs，否则使用ASDP knockoffs。
在这个案例中，我们使用估计的模型参数和完整的SDP结构生成二阶高斯knockoffs。然后，我们像往常一样运行knockoff filter。

gaussian_knockoffs = function(X) create.second_order(X, method='sdp',shrink = T)
result = knockoff.filter(X, y, knockoffs = gaussian_knockoffs)
print(result)

## Call:
## knockoff.filter(X = X, y = y, knockoffs = gaussian_knockoffs)
## 
## Selected variables:
##  [1]   9  40  44  46  61  78  85 108 141 146 148 153 172 173 177 210 238
## [18] 248 274 295 301 302 319 326 334 343 360 364 378 384 389 421 426 428
## [35] 451 494 506 510 528 559 595 617 651 668 676 682 702 708 718 770 775
## [52] 836 844 875 893 906 913 931 937 953 959

fdp(result$selected)

## [1] 0.1967213

等相关knockoffs

等效相关knockoffs提供了一种计算上更廉价的SDP knockoffs替代方案，代价是降低了统计功效。在这个例子中使用估计的模型参数和等相关结构生成二阶高斯knockoffs。然后我们运行knockoff filter。

guassian_knockoffs = function(X) create.second_order(X, method='equi', shrink=T)
result = knockoff.filter(X, y, knockoffs = gaussian_knockoffs)
print(result)

## Call:
## knockoff.filter(X = X, y = y, knockoffs = gaussian_knockoffs)
## 
## Selected variables:
##  [1]   3   9  40  46  61  78  85 108 148 153 172 173 177 210 223 238 248
## [18] 281 295 301 326 334 343 360 364 378 384 389 421 426 428 451 494 506
## [35] 557 559 595 651 668 682 702 708 770 775 787 844 893 906 913 931 937
## [52] 953 959

fdp(result$selected)

## [1] 0.0754717

另见

如果想查看knockoff filter，请参阅高级简介。如果想查看如何对Fixed-X变量使用knockoffs，请参阅Fixed-X简介。