Skip to content

Commit 1535c7c

Browse files
kevinbfrykevinbfry
committed
testing travisci
1 parent 7e16174 commit 1535c7c

File tree

1 file changed

+137
-137
lines changed

1 file changed

+137
-137
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+
# }

0 commit comments

Comments
 (0)