Skip to content

Commit 749c136

Browse files
tibshiranitibshirani
committed
rob added unifTest.R
1 parent c5bf9f6 commit 749c136

File tree

3 files changed

+149
-26
lines changed

3 files changed

+149
-26
lines changed

‎selectiveInference/R/RcppExports.R

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,19 @@
1-
# This file was generated by Rcpp::compileAttributes
1+
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
22
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
33

44
solve_QP <- function(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active) {
5-
.Call('selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)
5+
.Call('_selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)
66
}
77

88
solve_QP_wide <- function(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active) {
9-
.Call('selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active)
9+
.Call('_selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active)
1010
}
1111

1212
update1_ <- function(Q2, w, m, k) {
13-
.Call('selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)
13+
.Call('_selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)
1414
}
1515

1616
downdate1_ <- function(Q1, R, j0, m, n) {
17-
.Call('selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)
17+
.Call('_selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)
1818
}
1919

‎selectiveInference/src/RcppExports.cpp

Lines changed: 34 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// This file was generated by Rcpp::compileAttributes
1+
// Generated by using Rcpp::compileAttributes() -> do not edit by hand
22
// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
33

44
#include <Rcpp.h>
@@ -7,10 +7,10 @@ using namespace Rcpp;
77

88
// solve_QP
99
Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma, double bound, int maxiter, Rcpp::NumericVector theta, Rcpp::NumericVector linear_func, Rcpp::NumericVector gradient, Rcpp::IntegerVector ever_active, Rcpp::IntegerVector nactive, double kkt_tol, double objective_tol, int max_active);
10-
RcppExport SEXP selectiveInference_solve_QP(SEXP SigmaSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
10+
RcppExport SEXP _selectiveInference_solve_QP(SEXP SigmaSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
1111
BEGIN_RCPP
12-
Rcpp::RObject __result;
13-
Rcpp::RNGScope __rngScope;
12+
Rcpp::RObject rcpp_result_gen;
13+
Rcpp::RNGScope rcpp_rngScope_gen;
1414
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Sigma(SigmaSEXP);
1515
Rcpp::traits::input_parameter< double >::type bound(boundSEXP);
1616
Rcpp::traits::input_parameter< int >::type maxiter(maxiterSEXP);
@@ -22,16 +22,16 @@ BEGIN_RCPP
2222
Rcpp::traits::input_parameter< double >::type kkt_tol(kkt_tolSEXP);
2323
Rcpp::traits::input_parameter< double >::type objective_tol(objective_tolSEXP);
2424
Rcpp::traits::input_parameter< int >::type max_active(max_activeSEXP);
25-
__result = Rcpp::wrap(solve_QP(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active));
26-
return __result;
25+
rcpp_result_gen = Rcpp::wrap(solve_QP(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active));
26+
return rcpp_result_gen;
2727
END_RCPP
2828
}
2929
// solve_QP_wide
3030
Rcpp::List solve_QP_wide(Rcpp::NumericMatrix X, double bound, int maxiter, Rcpp::NumericVector theta, Rcpp::NumericVector linear_func, Rcpp::NumericVector gradient, Rcpp::NumericVector X_theta, Rcpp::IntegerVector ever_active, Rcpp::IntegerVector nactive, double kkt_tol, double objective_tol, int max_active);
31-
RcppExport SEXP selectiveInference_solve_QP_wide(SEXP XSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP X_thetaSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
31+
RcppExport SEXP _selectiveInference_solve_QP_wide(SEXP XSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP X_thetaSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
3232
BEGIN_RCPP
33-
Rcpp::RObject __result;
34-
Rcpp::RNGScope __rngScope;
33+
Rcpp::RObject rcpp_result_gen;
34+
Rcpp::RNGScope rcpp_rngScope_gen;
3535
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type X(XSEXP);
3636
Rcpp::traits::input_parameter< double >::type bound(boundSEXP);
3737
Rcpp::traits::input_parameter< int >::type maxiter(maxiterSEXP);
@@ -44,36 +44,49 @@ BEGIN_RCPP
4444
Rcpp::traits::input_parameter< double >::type kkt_tol(kkt_tolSEXP);
4545
Rcpp::traits::input_parameter< double >::type objective_tol(objective_tolSEXP);
4646
Rcpp::traits::input_parameter< int >::type max_active(max_activeSEXP);
47-
__result = Rcpp::wrap(solve_QP_wide(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active));
48-
return __result;
47+
rcpp_result_gen = Rcpp::wrap(solve_QP_wide(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active));
48+
return rcpp_result_gen;
4949
END_RCPP
5050
}
5151
// update1_
5252
Rcpp::List update1_(Rcpp::NumericMatrix Q2, Rcpp::NumericVector w, int m, int k);
53-
RcppExport SEXP selectiveInference_update1_(SEXP Q2SEXP, SEXP wSEXP, SEXP mSEXP, SEXP kSEXP) {
53+
RcppExport SEXP _selectiveInference_update1_(SEXP Q2SEXP, SEXP wSEXP, SEXP mSEXP, SEXP kSEXP) {
5454
BEGIN_RCPP
55-
Rcpp::RObject __result;
56-
Rcpp::RNGScope __rngScope;
55+
Rcpp::RObject rcpp_result_gen;
56+
Rcpp::RNGScope rcpp_rngScope_gen;
5757
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Q2(Q2SEXP);
5858
Rcpp::traits::input_parameter< Rcpp::NumericVector >::type w(wSEXP);
5959
Rcpp::traits::input_parameter< int >::type m(mSEXP);
6060
Rcpp::traits::input_parameter< int >::type k(kSEXP);
61-
__result = Rcpp::wrap(update1_(Q2, w, m, k));
62-
return __result;
61+
rcpp_result_gen = Rcpp::wrap(update1_(Q2, w, m, k));
62+
return rcpp_result_gen;
6363
END_RCPP
6464
}
6565
// downdate1_
6666
Rcpp::List downdate1_(Rcpp::NumericMatrix Q1, Rcpp::NumericMatrix R, int j0, int m, int n);
67-
RcppExport SEXP selectiveInference_downdate1_(SEXP Q1SEXP, SEXP RSEXP, SEXP j0SEXP, SEXP mSEXP, SEXP nSEXP) {
67+
RcppExport SEXP _selectiveInference_downdate1_(SEXP Q1SEXP, SEXP RSEXP, SEXP j0SEXP, SEXP mSEXP, SEXP nSEXP) {
6868
BEGIN_RCPP
69-
Rcpp::RObject __result;
70-
Rcpp::RNGScope __rngScope;
69+
Rcpp::RObject rcpp_result_gen;
70+
Rcpp::RNGScope rcpp_rngScope_gen;
7171
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Q1(Q1SEXP);
7272
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type R(RSEXP);
7373
Rcpp::traits::input_parameter< int >::type j0(j0SEXP);
7474
Rcpp::traits::input_parameter< int >::type m(mSEXP);
7575
Rcpp::traits::input_parameter< int >::type n(nSEXP);
76-
__result = Rcpp::wrap(downdate1_(Q1, R, j0, m, n));
77-
return __result;
76+
rcpp_result_gen = Rcpp::wrap(downdate1_(Q1, R, j0, m, n));
77+
return rcpp_result_gen;
7878
END_RCPP
7979
}
80+
81+
static const R_CallMethodDef CallEntries[] = {
82+
{"_selectiveInference_solve_QP", (DL_FUNC) &_selectiveInference_solve_QP, 11},
83+
{"_selectiveInference_solve_QP_wide", (DL_FUNC) &_selectiveInference_solve_QP_wide, 12},
84+
{"_selectiveInference_update1_", (DL_FUNC) &_selectiveInference_update1_, 4},
85+
{"_selectiveInference_downdate1_", (DL_FUNC) &_selectiveInference_downdate1_, 5},
86+
{NULL, NULL, 0}
87+
};
88+
89+
RcppExport void R_init_selectiveInference(DllInfo *dll) {
90+
R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
91+
R_useDynamicSymbols(dll, FALSE);
92+
}

‎tests/unifTest.R

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
2+
library(selectiveInference)
3+
4+
library(glmnet)
5+
6+
set.seed(424)
7+
8+
#n=100
9+
#p=30
10+
11+
n=20
12+
p=40
13+
sigma=.4
14+
beta=c(3,2,-1,4,-2,2,rep(0,p-6))
15+
#beta=rep(0,p)
16+
17+
tr=beta!=0
18+
19+
#type="full"
20+
type="part"
21+
22+
nsim = 1000
23+
lambda=.3
24+
nzb=0
25+
pvals <- matrix(NA, nrow=nsim, ncol=p)
26+
x = matrix(rnorm(n*p),n,p)
27+
x = scale(x,T,T)/sqrt(n-1)
28+
mu = x%*%beta
29+
30+
for (i in 1:nsim) {
31+
cat(i)
32+
y=mu+sigma*rnorm(n)
33+
y=y-mean(y)
34+
# first run glmnet
35+
gfit=glmnet(x,y,intercept=F,standardize=F,thresh=1e-8)
36+
37+
#extract coef for a given lambda; Note the 1/n factor!
38+
bhat = coef(gfit, s=lambda/n, exact=TRUE,x=x,y=y)[-1]
39+
nzb=nzb+sum(bhat!=0)
40+
# compute fixed lambda p-values and selection intervals
41+
aa = fixedLassoInf(x,y,bhat,lambda,intercept=F,sigma=sigma,type=type)
42+
pvals[i, aa$vars] <- aa$pv
43+
}
44+
45+
# summarize results
46+
47+
if(type=="partial"){
48+
nulls=rowSums(is.na(pvals[,tr]))==0 # for type=partial, nonnull setting
49+
np = pvals[nulls,-(1:sum(beta!=0))]
50+
}
51+
52+
if(type=="full"){
53+
nulls=1:nrow(pvals) # for type=full non null setting
54+
np = pvals[nulls,-(1:sum(beta!=0))]
55+
}
56+
57+
58+
59+
#np=pvals #for null setting
60+
61+
o=!is.na(np)
62+
63+
#check uniformity
64+
65+
plot((1:sum(o))/sum(o),sort(np[o]),xlab="Expected pvalue",ylab="Observed pvalue")
66+
abline(0,1)
67+
68+
69+
# estimate and plot FDR
70+
71+
pvadj=pvadj.by=pvadj.holm=matrix(NA,nsim,p)
72+
for(ii in 1:nsim){
73+
o=!is.na(pvals[ii,])
74+
pvadj[ii,o]=p.adjust(pvals[ii,o],method="BH")
75+
pvadj.by[ii,o]=p.adjust(pvals[ii,o],method="BY")
76+
pvadj.holm[ii,o]=p.adjust(pvals[ii,o],method="holm")
77+
}
78+
qqlist=fdr=se=fdr.by=se.by=fdr.holm=se.holm=c(.05, .1,.15,.2,.25,.3)
79+
jj=0
80+
for(qq in qqlist){
81+
jj=jj+1
82+
83+
r=v=r.by=v.by=r.holm=v.holm=rep(NA,nsim)
84+
for(ii in 1:nsim){
85+
v[ii]=sum( (pvadj[ii,]<qq & !tr), na.rm=T)
86+
r[ii]=sum( (pvadj[ii,]<qq), na.rm=T)
87+
v.by[ii]=sum( (pvadj.by[ii,]<qq & !tr), na.rm=T)
88+
r.by[ii]=sum( (pvadj.by[ii,]<qq), na.rm=T)
89+
v.holm[ii]=sum( (pvadj.holm[ii,]<qq & !tr), na.rm=T)
90+
r.holm[ii]=sum( (pvadj.holm[ii,]<qq), na.rm=T)
91+
92+
}
93+
oo=r!=0
94+
fdr[jj]=mean((v/r)[oo])
95+
se[jj]=sqrt(var((v/r)[oo])/sum(oo))
96+
oo=r.by!=0
97+
fdr.by[jj]=mean((v.by/r.by)[oo])
98+
se.by[jj]=sqrt(var((v.by/r.by)[oo])/sum(oo))
99+
oo=r.by!=0
100+
fdr.holm[jj]=mean((v.holm/r.holm)[oo])
101+
se.holm[jj]=sqrt(var((v.holm/r.holm)[oo])/sum(oo))
102+
}
103+
104+
105+
plot(qqlist,fdr,type="b",xlab="target FDR",ylab="observed FDR",ylim=c(0,.6),xlim=c(0,.6))
106+
lines(qqlist,fdr.by,type="b",col=3)
107+
lines(qqlist,fdr.holm,type="b",col=4)
108+
abline(0,1,lty=2)
109+
title(paste("n=",as.character(n)," p=",as.character(p)," ",as.character(type)))
110+
legend("bottomright",c("BH","BY","Holm"),col=c(1,3,4),lty=1)

0 commit comments

Comments
 (0)