now debias lasso only computes theta for the active set

Author: kevinbfry

selectiveInference/R/funs.fixed.R

@@ -148,12 +148,15 @@ sigma=NULL, alpha=0.1,
       hsigma <- 1/n*(t(Xordered)%*%Xordered)
       hsigmaS <- 1/n*(t(XS)%*%XS) # hsigma[S,S]
       hsigmaSinv <- solve(hsigmaS) # pinv(hsigmaS)
-      
+
       # Approximate inverse covariance matrix for when (n < p) from lasso_Inference.R
-      htheta <- InverseLinfty(hsigma, n, verbose=FALSE)
+      htheta <- InverseLinfty(hsigma, n, length(S), verbose=FALSE)
+      # htheta <- InverseLinfty(hsigma, n, verbose=FALSE)
 
       FS = rbind(diag(length(S)),matrix(0,pp-length(S),length(S)))
-      ithetasigma = (diag(pp)-(htheta%*%hsigma))
+      GS = cbind(diag(length(S)),matrix(0,length(S),pp-length(S)))
+      ithetasigma = (GS-(htheta%*%hsigma))
+      # ithetasigma = (diag(pp) - (htheta%*%hsigma))
 
       M <- (((htheta%*%t(Xordered))+ithetasigma%*%FS%*%hsigmaSinv%*%t(XS))/n)
       # vector which is offset for testing debiased beta's
@@ -254,20 +257,21 @@ fixedLasso.poly=
 
 ##############################
 
-### Functions borrowed from lasso_inference.R
+### Functions borrowed and slightly modified from lasso_inference.R
 
 ## Approximates inverse covariance matrix theta
-InverseLinfty <- function(sigma, n, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE) {
-  isgiven <- 1;
+InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE) {
+  # InverseLinfty <- function(sigma, n, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE) {
+    isgiven <- 1;
   if (is.null(mu)){
     isgiven <- 0;
   }
 
   p <- nrow(sigma);
-  M <- matrix(0, p, p);
+  M <- matrix(0, e, p);
   xperc = 0;
   xp = round(p/10);
-  for (i in 1:p) {
+  for (i in 1:e) {
     if ((i %% xp)==0){
       xperc = xperc+10;
       if (verbose) {

selectiveInference/R/linear.tests.R

@@ -15,7 +15,7 @@ robs.test <- function() {
   hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=T))
 
 
-  return(fixedLassoInf(X,y,hbeta[-1],lambda,family="gaussian",type="partial",intercept=T,sigma=sigma)$pv)
+  return(fixedLassoInf(X,y,hbeta[-1],lambda,family="gaussian",type="partial",intercept=T,sigma=sigma))
 }
 
 
@@ -134,39 +134,4 @@ power.partial.pval.dist <- function(n,p,intercept=T,lambda=1) {
   }
   qqplot(x=runif(length(pvs)),y=pvs,xlab="Expected",ylab="Observed",main="Partial Coef. 10 Corr. X")
   abline(0,1)
-}
-
-
-
-
-## Tests pop inf for X and y randomly generated
-nullPopTest <- function(n,p,intercept=T,lambda=1) {
-  y <- matrix(rnorm(n),ncol=1) # rand N(0,1) response
-  X <- matrix(rnorm(p*n),ncol = p) # p rand N(0,1) predictors
-  
-  # lambda <- 1
-  X=scale(X,T,T)/sqrt(n-1)
-  
-  # lambda <- 1
-  las <- glmnet(X,y,family="gaussian",alpha=1,standardize=F,intercept=intercept)
-  hbeta <- as.numeric(coef(las,x=X,y=y,s=lambda/n,exact=TRUE,intercept=intercept))
-  
-  ### perform post selection inference
-  
-  sigma = estimateSigma(X,y)$sigmahat
-  
-  if (intercept) return(fixedLassoInf(X,y,hbeta,lambda,family="gaussian",type="full",intercept=intercept,sigma=sigma))
-  else return(fixedLassoInf(X,y,hbeta[-1],lambda,family="gaussian",type="full",intercept=intercept,sigma=sigma))
-}
-
-## QQ plot of p-values for all null data now that bug fix is implemented
-null.pop.pval.dist <- function(n,p,intercept=T,lambda=1) {
-  pvs <- c()
-  for(i in 1:10) {
-    a <- nullPopTest(n,p,intercept,lambda)
-    ps <- a$pv
-    pvs <- c(pvs,ps,recursive=T)
-  }
-  qqplot(x=runif(length(pvs)),y=pvs,xlab="Expected",ylab="Observed",main="Pop Coef. Null X")
-  abline(0,1)
 }