Skip to content

Commit be89769

Browse files
jonathan-taylorjonathan-taylor
committed
Merge branch 'master' of github.com:selective-inference/R-software into debias_lasso_linear
2 parents 411143e + 1672cee commit be89769

File tree

2 files changed

+139
-138
lines changed

2 files changed

+139
-138
lines changed

‎selectiveInference/R/linear.tests.R

Lines changed: 137 additions & 137 deletions
Original file line numberDiff line numberDiff line change
@@ -1,137 +1,137 @@
1-
robs.test <- function() {
2-
n <- 100
3-
p <- 200
4-
5-
set.seed(11332)
6-
7-
y <- matrix(rnorm(n),ncol=1) # rand N(0,1) response
8-
X <- matrix(rnorm(p*n),ncol = p) # p rand N(0,1) predictors
9-
10-
X=scale(X,T,T)/sqrt(n-1)
11-
lambda=1
12-
sigma = estimateSigma(X,y)$sigmahat
13-
14-
las <- glmnet(X,y,family="gaussian",alpha=1,standardize=F,intercept=T)
15-
hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=T))
16-
17-
18-
return(fixedLassoInf(X,y,hbeta[-1],lambda,family="gaussian",type="partial",intercept=T,sigma=sigma))
19-
}
20-
21-
22-
## Tests partial inf for X and y randomly generated from N(0,1)
23-
nullTest <- function(X,y,lambda,intercept,type=c("full","partial")) {
24-
n=nrow(X)
25-
X=scale(X,T,T)/sqrt(n-1)
26-
27-
sigma = estimateSigma(X,y)$sigmahat
28-
29-
las <- glmnet(X,y,family="gaussian",alpha=1,standardize=F,intercept=intercept)
30-
hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=intercept))
31-
32-
if (type=="partial" || intercept==F) hbeta = hbeta[-1]
33-
34-
return(fixedLassoInf(X,y,hbeta,lambda,family="gaussian",type=type,intercept=intercept,sigma=sigma))
35-
}
36-
37-
## Test partial inf for X and y where 10 variables are y with random additive N(0,0.5) noise
38-
corrTest <- function(X,y,lambda,intercept,type=c("full","partial")) {
39-
n=nrow(X)
40-
corr.X = rep(y,10) + matrix(rnorm(n*10,0,0.5),ncol = 10)
41-
X = cbind(corr.X,X)
42-
X=scale(X,T,T)/sqrt(n-1)
43-
44-
las <- glmnet(X,y,family="gaussian",alpha=1,standardize=F,intercept=intercept)
45-
hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=intercept))
46-
47-
sigma = estimateSigma(X,y)$sigmahat
48-
49-
if (type=="partial" || intercept==F) hbeta = hbeta[-1]
50-
51-
return(fixedLassoInf(X,y,hbeta,lambda,family="gaussian",type=type,intercept=intercept,sigma=sigma))
52-
}
53-
54-
## QQ plot of p-values for all null data now that bug fix is implemented
55-
partial.qq.test <- function() {
56-
n <- 100
57-
p <- 200
58-
59-
lambda=1
60-
61-
null.int.pvs <- c()
62-
corr.int.pvs <- c()
63-
null.pvs <- c()
64-
corr.pvs <- c()
65-
for(i in 1:25) {
66-
y <- matrix(rnorm(n),ncol=1) # rand N(0,1) response
67-
X <- matrix(rnorm(p*n),ncol=p) # p rand N(0,1) predictors
68-
69-
null <- nullTest(X,y,lambda,F,type="partial")
70-
corr <- corrTest(X,y,lambda,F,type="partial")
71-
null.pvs <- c(null.pvs,null$pv,recursive=T)
72-
corr.pvs <- c(corr.pvs,corr$pv,recursive=T)
73-
null.int <- nullTest(X,y,lambda,T,type="partial")
74-
corr.int <- corrTest(X,y,lambda,T,type="partial")
75-
null.int.pvs <- c(null.int.pvs,null.int$pv,recursive=T)
76-
corr.int.pvs <- c(corr.int.pvs,corr.int$pv,recursive=T)
77-
}
78-
79-
qqplot(x=runif(length(null.pvs)),y=null.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. Null X w/o Intercept")
80-
abline(0,1)
81-
qqplot(x=runif(length(corr.pvs)),y=corr.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. 10 Corr. X w/o Intercept")
82-
abline(0,1)
83-
qqplot(x=runif(length(null.int.pvs)),y=null.int.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. Null X w/ Intercept")
84-
abline(0,1)
85-
qqplot(x=runif(length(corr.int.pvs)),y=corr.int.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. 10 Corr. X w/ Intercept")
86-
abline(0,1)
87-
}
88-
89-
## QQ plot of p-values for all null data now that bug fix is implemented
90-
pop.qq.test <- function() {
91-
n <- 100
92-
p <- 200
93-
94-
lambda=1
95-
96-
null.int.pvs <- c()
97-
corr.int.pvs <- c()
98-
null.pvs <- c()
99-
corr.pvs <- c()
100-
for(i in 1:25) {
101-
y <- matrix(rnorm(n),ncol=1) # rand N(0,1) response
102-
X <- matrix(rnorm(p*n),ncol=p) # p rand N(0,1) predictors
103-
104-
null <- nullTest(X,y,lambda,F,type="full")
105-
corr <- corrTest(X,y,lambda,F,type="full")
106-
null.pvs <- c(null.pvs,null$pv,recursive=T)
107-
corr.pvs <- c(corr.pvs,corr$pv,recursive=T)
108-
null.int <- nullTest(X,y,lambda,T,type="full")
109-
corr.int <- corrTest(X,y,lambda,T,type="full")
110-
null.int.pvs <- c(null.int.pvs,null.int$pv,recursive=T)
111-
corr.int.pvs <- c(corr.int.pvs,corr.int$pv,recursive=T)
112-
}
113-
114-
qqplot(x=runif(length(null.pvs)),y=null.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. Null X w/o Intercept")
115-
abline(0,1)
116-
qqplot(x=runif(length(corr.pvs)),y=corr.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. 10 Corr. X w/o Intercept")
117-
abline(0,1)
118-
qqplot(x=runif(length(null.int.pvs)),y=null.int.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. Null X w/ Intercept")
119-
abline(0,1)
120-
qqplot(x=runif(length(corr.int.pvs)),y=corr.int.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. 10 Corr. X w/ Intercept")
121-
abline(0,1)
122-
}
123-
124-
125-
126-
127-
## QQ plot of p-values for data with correlated x now that bug fix implemented
128-
power.partial.pval.dist <- function(n,p,intercept=T,lambda=1) {
129-
pvs <- c()
130-
for(i in 1:10) {
131-
a <- powerPartialTest(n,p,intercept,lambda)
132-
ps <- a$pv
133-
pvs <- c(pvs,ps,recursive=T)
134-
}
135-
qqplot(x=runif(length(pvs)),y=pvs,xlab="Expected",ylab="Observed",main="Partial Coef. 10 Corr. X")
136-
abline(0,1)
137-
}
1+
# robs.test <- function() {
2+
# n <- 100
3+
# p <- 200
4+
#
5+
# set.seed(11332)
6+
#
7+
# y <- matrix(rnorm(n),ncol=1) # rand N(0,1) response
8+
# X <- matrix(rnorm(p*n),ncol = p) # p rand N(0,1) predictors
9+
#
10+
# X=scale(X,T,T)/sqrt(n-1)
11+
# lambda=1
12+
# sigma = estimateSigma(X,y)$sigmahat
13+
#
14+
# las <- glmnet(X,y,family="gaussian",alpha=1,standardize=F,intercept=T)
15+
# hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=T))
16+
#
17+
#
18+
# return(fixedLassoInf(X,y,hbeta[-1],lambda,family="gaussian",type="partial",intercept=T,sigma=sigma))
19+
# }
20+
#
21+
#
22+
# ## Tests partial inf for X and y randomly generated from N(0,1)
23+
# nullTest <- function(X,y,lambda,intercept,type=c("full","partial")) {
24+
# n=nrow(X)
25+
# X=scale(X,T,T)/sqrt(n-1)
26+
#
27+
# sigma = estimateSigma(X,y)$sigmahat
28+
#
29+
# las <- glmnet(X,y,family="gaussian",alpha=1,standardize=F,intercept=intercept)
30+
# hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=intercept))
31+
#
32+
# if (type=="partial" || intercept==F) hbeta = hbeta[-1]
33+
#
34+
# return(fixedLassoInf(X,y,hbeta,lambda,family="gaussian",type=type,intercept=intercept,sigma=sigma))
35+
# }
36+
#
37+
# ## Test partial inf for X and y where 10 variables are y with random additive N(0,0.5) noise
38+
# corrTest <- function(X,y,lambda,intercept,type=c("full","partial")) {
39+
# n=nrow(X)
40+
# corr.X = rep(y,10) + matrix(rnorm(n*10,0,0.5),ncol = 10)
41+
# X = cbind(corr.X,X)
42+
# X=scale(X,T,T)/sqrt(n-1)
43+
#
44+
# las <- glmnet(X,y,family="gaussian",alpha=1,standardize=F,intercept=intercept)
45+
# hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=intercept))
46+
#
47+
# sigma = estimateSigma(X,y)$sigmahat
48+
#
49+
# if (type=="partial" || intercept==F) hbeta = hbeta[-1]
50+
#
51+
# return(fixedLassoInf(X,y,hbeta,lambda,family="gaussian",type=type,intercept=intercept,sigma=sigma))
52+
# }
53+
#
54+
# ## QQ plot of p-values for all null data now that bug fix is implemented
55+
# partial.qq.test <- function() {
56+
# n <- 100
57+
# p <- 200
58+
#
59+
# lambda=1
60+
#
61+
# null.int.pvs <- c()
62+
# corr.int.pvs <- c()
63+
# null.pvs <- c()
64+
# corr.pvs <- c()
65+
# for(i in 1:25) {
66+
# y <- matrix(rnorm(n),ncol=1) # rand N(0,1) response
67+
# X <- matrix(rnorm(p*n),ncol=p) # p rand N(0,1) predictors
68+
#
69+
# null <- nullTest(X,y,lambda,F,type="partial")
70+
# corr <- corrTest(X,y,lambda,F,type="partial")
71+
# null.pvs <- c(null.pvs,null$pv,recursive=T)
72+
# corr.pvs <- c(corr.pvs,corr$pv,recursive=T)
73+
# null.int <- nullTest(X,y,lambda,T,type="partial")
74+
# corr.int <- corrTest(X,y,lambda,T,type="partial")
75+
# null.int.pvs <- c(null.int.pvs,null.int$pv,recursive=T)
76+
# corr.int.pvs <- c(corr.int.pvs,corr.int$pv,recursive=T)
77+
# }
78+
#
79+
# qqplot(x=runif(length(null.pvs)),y=null.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. Null X w/o Intercept")
80+
# abline(0,1)
81+
# qqplot(x=runif(length(corr.pvs)),y=corr.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. 10 Corr. X w/o Intercept")
82+
# abline(0,1)
83+
# qqplot(x=runif(length(null.int.pvs)),y=null.int.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. Null X w/ Intercept")
84+
# abline(0,1)
85+
# qqplot(x=runif(length(corr.int.pvs)),y=corr.int.pvs,xlab="Expected",ylab="Observed",main="Partial Coef. 10 Corr. X w/ Intercept")
86+
# abline(0,1)
87+
# }
88+
#
89+
# ## QQ plot of p-values for all null data now that bug fix is implemented
90+
# pop.qq.test <- function() {
91+
# n <- 100
92+
# p <- 200
93+
#
94+
# lambda=1
95+
#
96+
# null.int.pvs <- c()
97+
# corr.int.pvs <- c()
98+
# null.pvs <- c()
99+
# corr.pvs <- c()
100+
# for(i in 1:25) {
101+
# y <- matrix(rnorm(n),ncol=1) # rand N(0,1) response
102+
# X <- matrix(rnorm(p*n),ncol=p) # p rand N(0,1) predictors
103+
#
104+
# null <- nullTest(X,y,lambda,F,type="full")
105+
# corr <- corrTest(X,y,lambda,F,type="full")
106+
# null.pvs <- c(null.pvs,null$pv,recursive=T)
107+
# corr.pvs <- c(corr.pvs,corr$pv,recursive=T)
108+
# null.int <- nullTest(X,y,lambda,T,type="full")
109+
# corr.int <- corrTest(X,y,lambda,T,type="full")
110+
# null.int.pvs <- c(null.int.pvs,null.int$pv,recursive=T)
111+
# corr.int.pvs <- c(corr.int.pvs,corr.int$pv,recursive=T)
112+
# }
113+
#
114+
# qqplot(x=runif(length(null.pvs)),y=null.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. Null X w/o Intercept")
115+
# abline(0,1)
116+
# qqplot(x=runif(length(corr.pvs)),y=corr.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. 10 Corr. X w/o Intercept")
117+
# abline(0,1)
118+
# qqplot(x=runif(length(null.int.pvs)),y=null.int.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. Null X w/ Intercept")
119+
# abline(0,1)
120+
# qqplot(x=runif(length(corr.int.pvs)),y=corr.int.pvs,xlab="Expected",ylab="Observed",main="Pop Coef. 10 Corr. X w/ Intercept")
121+
# abline(0,1)
122+
# }
123+
#
124+
#
125+
#
126+
#
127+
# ## QQ plot of p-values for data with correlated x now that bug fix implemented
128+
# power.partial.pval.dist <- function(n,p,intercept=T,lambda=1) {
129+
# pvs <- c()
130+
# for(i in 1:10) {
131+
# a <- powerPartialTest(n,p,intercept,lambda)
132+
# ps <- a$pv
133+
# pvs <- c(pvs,ps,recursive=T)
134+
# }
135+
# qqplot(x=runif(length(pvs)),y=pvs,xlab="Expected",ylab="Observed",main="Partial Coef. 10 Corr. X")
136+
# abline(0,1)
137+
# }

‎selectiveInference/man/fixedLassoInf.Rd

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ fixed value of the tuning parameter lambda
1010
}
1111
\usage{
1212
fixedLassoInf(x, y, beta, lambda, family = c("gaussian", "binomial",
13-
"cox"),intercept=TRUE, status=NULL, sigma=NULL, alpha=0.1,
13+
"cox"),intercept=TRUE, add.targets=NULL, status=NULL, sigma=NULL, alpha=0.1,
1414
type=c("partial","full"), tol.beta=1e-5, tol.kkt=0.1,
1515
gridrange=c(-100,100), bits=NULL, verbose=FALSE)
1616
}
@@ -56,6 +56,7 @@ Significance level for confidence intervals (target is miscoverage alpha/2 in ea
5656
Was the lasso problem solved (e.g., by glmnet) with an intercept in the model?
5757
Default is TRUE. Must be TRUE for "binomial" family. Not used for 'cox" family, where no intercept is assumed.
5858
}
59+
\item{add.targets}{Optional vector of predictors to be included as targets of inference, regardless of whether or not they are selected by the lasso. Default is NULL.}
5960
\item{status}{Censoring status for Cox model; 1=failurem 0=censored}
6061
\item{type}{Contrast type for p-values and confidence intervals: default is
6162
"partial"---meaning that the contrasts tested are the partial population

0 commit comments

Comments
 (0)