通过Model-XKnockoffs控制variableSelection
通过Model-X Knockoffs控制variable Selection
翻译于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)
variable used in lambda
##
## 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。这种结构从数据中预估⾏的均值和协⽅差，⽽不是使⽤真实参数抽取变量。knockoff包还包括其他knockoff构造，所有这些⽅法的名称都以ate 为前缀。在下⼀个⽚段中，我们使⽤真实模型参数⽣成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 统计，它计，其中和分别是第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 论⽂中最简单的测试统计之⼀，即。
X mu sigma W =j ∣Z ∣−j ∣∣Z j Z j Z j W =j X ∗j y −∗X j 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简介。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。