Skip to content

Commit ce15b68

Browse files
jonathan-taylorjonathan-taylor
committed
checking comparison example again -- looks good
1 parent c2460a1 commit ce15b68

File tree

7 files changed

+1138
-7
lines changed

7 files changed

+1138
-7
lines changed

‎selectiveInference/R/funs.fixed.R

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ fixedLassoInf <- function(x, y, beta,
88
sigma=NULL, alpha=0.1,
99
type=c("partial", "full"), tol.beta=1e-5, tol.kkt=0.1,
1010
gridrange=c(-100,100), bits=NULL, verbose=FALSE,
11-
linesearch.try=10) {
11+
linesearch.try=10, offset_correction=TRUE) {
1212

1313
family = match.arg(family)
1414
this.call = match.call()
@@ -197,6 +197,9 @@ fixedLassoInf <- function(x, y, beta,
197197
M = M[-1,] # remove intercept row
198198
null_value = null_value[-1] # remove intercept element
199199
}
200+
if (!offset_correction) {
201+
null_value = 0 * null_value
202+
}
200203
} else if (type=="partial" || p > n) {
201204
xa = x[,vars,drop=F]
202205
M = pinv(crossprod(xa)) %*% t(xa)
@@ -325,13 +328,13 @@ debiasingMatrix = function(Xinfo, # could be X or t(X) %*% X / n d
325328
max_active=NULL, # how big can active set get?
326329
max_try=10, # how many steps in linesearch?
327330
warn_kkt=FALSE, # warn if KKT does not seem to be satisfied?
328-
max_iter=100, # how many iterations for each optimization problem
331+
max_iter=50, # how many iterations for each optimization problem
329332
kkt_stop=TRUE, # stop based on KKT conditions?
330333
parameter_stop=TRUE, # stop based on relative convergence of parameter?
331334
objective_stop=TRUE, # stop based on relative decrease in objective?
332335
kkt_tol=1.e-4, # tolerance for the KKT conditions
333336
parameter_tol=1.e-4, # tolerance for relative convergence of parameter
334-
objective_tol=1.e-8 # tolerance for relative decrease in objective
337+
objective_tol=1.e-4 # tolerance for relative decrease in objective
335338
) {
336339

337340

@@ -399,13 +402,13 @@ debiasingRow = function (Xinfo, # could be X or t(X) %*% X / n dep
399402
max_active=NULL, # how big can active set get?
400403
max_try=10, # how many steps in linesearch?
401404
warn_kkt=FALSE, # warn if KKT does not seem to be satisfied?
402-
max_iter=100, # how many iterations for each optimization problem
405+
max_iter=50, # how many iterations for each optimization problem
403406
kkt_stop=TRUE, # stop based on KKT conditions?
404407
parameter_stop=TRUE, # stop based on relative convergence of parameter?
405408
objective_stop=TRUE, # stop based on relative decrease in objective?
406409
kkt_tol=1.e-4, # tolerance for the KKT conditions
407410
parameter_tol=1.e-4, # tolerance for relative convergence of parameter
408-
objective_tol=1.e-8 # tolerance for relative decrease in objective
411+
objective_tol=1.e-4 # tolerance for relative decrease in objective
409412
) {
410413

411414
p = ncol(Xinfo)

‎selectiveInference/man/debiasingMatrix.Rd

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -22,13 +22,13 @@ debiasingMatrix(Xinfo,
2222
max_active=NULL,
2323
max_try=10,
2424
warn_kkt=FALSE,
25-
max_iter=100,
25+
max_iter=50,
2626
kkt_stop=TRUE,
2727
parameter_stop=TRUE,
2828
objective_stop=TRUE,
2929
kkt_tol=1.e-4,
3030
parameter_tol=1.e-4,
31-
objective_tol=1.e-8)
31+
objective_tol=1.e-4)
3232
}
3333
\arguments{
3434
\item{Xinfo}{

‎tests/debiased_lasso/comparison_scaled.R

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
source('javanmard_montanari.R')
2+
3+
##############################################
4+
5+
# Runs nsims simulations under the global null, computing p-values
6+
# using both the old code (slow one using Adel's code) and the new
7+
# code (faster using Jon's code), and produces qq-plots for both.
8+
# Runing 50 sims takes about 10-15 mins because old code is slow, so
9+
# feel free to lower nsims if you want
10+
11+
12+
library(selectiveInference)
13+
library(glmnet)
14+
15+
# set.seed(424)
16+
17+
n=100
18+
p=200
19+
20+
sigma=.5
21+
22+
theor_lambda = sigma * sqrt(2 * log(p))
23+
lambda=c(0.25, 0.5, 1, 0.8 * theor_lambda, theor_lambda)
24+
25+
for (j in c(3,4,5,1,2)) {
26+
27+
thresh = 1e-10
28+
29+
beta=rep(0,p)
30+
type="full"
31+
32+
nsim = 20
33+
34+
scaling = sqrt(n)
35+
pvs_old = c()
36+
pvs_new <- c()
37+
pvs_old_0 = c() # don't add the offset correction
38+
pvs_new_0 = c() # don't add the offset correction
39+
for (i in 1:nsim) {
40+
cat(i,fill=T)
41+
x = matrix(rnorm(n*p),n,p)
42+
x = scale(x,T,T) / scaling
43+
mu = x%*%beta
44+
y=mu+sigma*rnorm(n)
45+
46+
# first run glmnet
47+
gfit=glmnet(x,y,intercept=F,standardize=F,thresh=thresh)
48+
49+
bhat = coef(gfit, s=lambda[j]/(sqrt(n) * scaling), exact=TRUE,x=x,y=y)[-1]
50+
51+
if(sum(bhat != 0) > 0) {
52+
53+
# compute fixed lambda p-values and selection intervals
54+
55+
aa = fixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type)
56+
bb = oldFixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type)
57+
cc = fixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type, offset_correction=FALSE)
58+
dd = oldFixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type, offset_correction=FALSE)
59+
pvs_new <- c(pvs_new, aa$pv, recursive=TRUE)
60+
pvs_old <- c(pvs_old, bb$pv,recursive=TRUE)
61+
pvs_new_0 <- c(pvs_new_0, cc$pv, recursive=TRUE)
62+
pvs_old_0 <- c(pvs_old_0, dd$pv, recursive=TRUE)
63+
64+
cat()
65+
}
66+
}
67+
68+
#check uniformity
69+
70+
png(paste('comparison_scaled', j, '.png', sep=''))
71+
plot(ecdf(pvs_old), pch=23, col='green', xlim=c(0,1), ylim=c(0,1), main='ECDF of p-values')
72+
plot(ecdf(pvs_new), pch=24, col='purple', add=TRUE)
73+
plot(ecdf(pvs_old_0), pch=23, col='red', add=TRUE)
74+
plot(ecdf(pvs_new_0), pch=24, col='black', add=TRUE)
75+
abline(0,1)
76+
legend("bottomright", legend=c("Old","New", "Old 0", "New 0"), pch=c(23,24,23,24), pt.bg=c("green","purple","red","black"))
77+
dev.off()
78+
}

‎tests/debiased_lasso/comparison_unscaled.R

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
source('javanmard_montanari.R')
2+
3+
##############################################
4+
5+
# Runs nsims simulations under the global null, computing p-values
6+
# using both the old code (slow one using Adel's code) and the new
7+
# code (faster using Jon's code), and produces qq-plots for both.
8+
# Runing 50 sims takes about 10-15 mins because old code is slow, so
9+
# feel free to lower nsims if you want
10+
11+
12+
library(selectiveInference)
13+
library(glmnet)
14+
15+
# set.seed(424)
16+
17+
n=100
18+
p=200
19+
20+
sigma=.5
21+
22+
theor_lambda = sigma * sqrt(2 * log(p))
23+
lambda=c(0.25, 0.5, 1, 0.8 * theor_lambda, theor_lambda)
24+
25+
for (j in c(3,4,5,1,2)) {
26+
27+
thresh = 1e-10
28+
29+
beta=rep(0,p)
30+
type="full"
31+
32+
nsim = 20
33+
34+
scaling = sqrt(n)
35+
pvs_old = c()
36+
pvs_new <- c()
37+
pvs_old_0 = c() # don't add the offset correction
38+
pvs_new_0 = c() # don't add the offset correction
39+
for (i in 1:nsim) {
40+
cat(i,fill=T)
41+
x = matrix(rnorm(n*p),n,p)
42+
x = scale(x,T,T) / scaling
43+
mu = x%*%beta
44+
y=mu+sigma*rnorm(n)
45+
46+
# first run glmnet
47+
gfit=glmnet(x,y,intercept=F,standardize=F,thresh=thresh)
48+
49+
bhat = coef(gfit, s=lambda[j]/(sqrt(n) * scaling), exact=TRUE,x=x,y=y)[-1]
50+
51+
if(sum(bhat != 0) > 0) {
52+
53+
# compute fixed lambda p-values and selection intervals
54+
55+
aa = fixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type)
56+
bb = oldFixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type)
57+
cc = fixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type, offset_correction=FALSE)
58+
dd = oldFixedLassoInf(x,y,bhat,lambda[j]*sqrt(n) / scaling,intercept=F,sigma=sigma,type=type, offset_correction=FALSE)
59+
pvs_new <- c(pvs_new, aa$pv, recursive=TRUE)
60+
pvs_old <- c(pvs_old, bb$pv,recursive=TRUE)
61+
pvs_new_0 <- c(pvs_new_0, cc$pv, recursive=TRUE)
62+
pvs_old_0 <- c(pvs_old_0, dd$pv, recursive=TRUE)
63+
64+
cat()
65+
}
66+
}
67+
68+
#check uniformity
69+
70+
png(paste('comparison_unscaled', j, '.png', sep=''))
71+
plot(ecdf(pvs_old), pch=23, col='green', xlim=c(0,1), ylim=c(0,1), main='ECDF of p-values')
72+
plot(ecdf(pvs_new), pch=24, col='purple', add=TRUE)
73+
plot(ecdf(pvs_old_0), pch=23, col='red', add=TRUE)
74+
plot(ecdf(pvs_new_0), pch=24, col='black', add=TRUE)
75+
abline(0,1)
76+
legend("bottomright", legend=c("Old","New", "Old 0", "New 0"), pch=c(23,24,23,24), pt.bg=c("green","purple","red","black"))
77+
dev.off()
78+
}

0 commit comments

Comments
 (0)