Merge pull request #36 from jonathan-taylor/master

fix of debiasingMatrix path, beginning of randomized LASSO

Author: tibshirani

selectiveInference/DESCRIPTION

@@ -9,7 +9,7 @@ Maintainer: Rob Tibshirani <[email protected]>
 Depends:
     glmnet,
     intervals,
-    survival
+    survival,
 Suggests:
     Rmpfr
 Description: New tools for post-selection inference, for use with forward

selectiveInference/NAMESPACE

@@ -43,5 +43,6 @@ importFrom("stats", dnorm, lsfit, pexp, pnorm, predict,
            qnorm, rnorm, sd, uniroot, dchisq, model.matrix, pchisq)
 importFrom("stats", "coef", "df", "lm", "pf")
 importFrom("stats", "glm", "residuals", "vcov")
+importFrom("stats", "rbinom", "rexp")
 importFrom("Rcpp", "sourceCpp")
 

selectiveInference/R/RcppExports.R

@@ -6,7 +6,7 @@ fixedLassoInf <- function(x, y, beta,
                           lambda, family=c("gaussian","binomial","cox"),
                           intercept=TRUE, add.targets=NULL, status=NULL,
                           sigma=NULL, alpha=0.1,
-                          type=c("partial","full"), tol.beta=1e-5, tol.kkt=0.1,
+                          type=c("partial", "full"), tol.beta=1e-5, tol.kkt=0.1,
                           gridrange=c(-100,100), bits=NULL, verbose=FALSE, 
                           linesearch.try=10) {
 
@@ -150,7 +150,7 @@ fixedLassoInf <- function(x, y, beta,
     ci = tailarea = matrix(0,k,2)
 
     if (type=="full" & p > n) {
-      if (intercept == T) {
+      if (intercept == TRUE) {
         pp=p+1
         Xint <- cbind(rep(1,n),x)
         # indices of selected predictors
@@ -189,8 +189,10 @@ fixedLassoInf <- function(x, y, beta,
       }
 
       M <- (((htheta%*%t(Xordered))+ithetasigma%*%FS%*%hsigmaSinv%*%t(XS))/n)
+
       # vector which is offset for testing debiased beta's
       null_value <- (((ithetasigma%*%FS%*%hsigmaSinv)%*%sign(hbetaS))*lambda/n)
+
       if (intercept == T) {
         M = M[-1,] # remove intercept row
         null_value = null_value[-1] # remove intercept element
@@ -238,12 +240,23 @@ fixedLassoInf <- function(x, y, beta,
     tailarea[j,] = a$tailarea
   }
 
-  out = list(type=type,lambda=lambda,pv=pv,ci=ci,
-    tailarea=tailarea,vlo=vlo,vup=vup,vmat=vmat,y=y,
-    vars=vars,sign=sign_vars,sigma=sigma,alpha=alpha,
-    sd=sigma*sqrt(rowSums(vmat^2)),
-    coef0=vmat%*%y,
-    call=this.call)
+  out = list(type=type,
+             lambda=lambda,
+             pv=pv,
+             ci=ci,
+             tailarea=tailarea,
+             vlo=vlo,
+             vup=vup,
+             vmat=vmat,
+             y=y,
+             vars=vars,
+             sign=sign_vars,
+             sigma=sigma,
+             alpha=alpha,
+             sd=sigma*sqrt(rowSums(vmat^2)),
+             coef0=vmat%*%y,
+             call=this.call)
+
   class(out) = "fixedLassoInf"
   return(out)
 }
@@ -306,15 +319,19 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
                            nsample, 
                            rows, 
  		           verbose=FALSE, 
-		           mu=NULL,             # starting value of mu
+		           bound=NULL,             # starting value of bound
    			   linesearch=TRUE,     # do a linesearch?
    		           scaling_factor=1.5,  # multiplicative factor for linesearch
 			   max_active=NULL,     # how big can active set get?
 			   max_try=10,          # how many steps in linesearch?
 			   warn_kkt=FALSE,      # warn if KKT does not seem to be satisfied?
-			   max_iter=100,        # how many iterations for each optimization problem
+			   max_iter=50,         # how many iterations for each optimization problem
+                           kkt_stop=TRUE,       # stop based on KKT conditions?
+                           parameter_stop=TRUE, # stop based on relative convergence of parameter?
+			   objective_stop=TRUE, # stop based on relative decrease in objective?
                            kkt_tol=1.e-4,       # tolerance for the KKT conditions
-			   objective_tol=1.e-8  # tolerance for relative decrease in objective
+                           parameter_tol=1.e-4, # tolerance for relative convergence of parameter
+			   objective_tol=1.e-4  # tolerance for relative decrease in objective
                            ) {
 
 
@@ -325,8 +342,8 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
   p = ncol(Xinfo);
   M = matrix(0, length(rows), p);
 
-  if (is.null(mu)) {
-      mu = (1/sqrt(nsample)) * qnorm(1-(0.1/(p^2)))
+  if (is.null(bound)) {
+      bound = (1/sqrt(nsample)) * qnorm(1-(0.1/(p^2)))
   }
 
   xperc = 0;
@@ -342,14 +359,18 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
     output = debiasingRow(Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
                           is_wide,
                           row,
-                          mu,
+                          bound,
                           linesearch=linesearch,
                           scaling_factor=scaling_factor,
                           max_active=max_active,
 			  max_try=max_try,
 			  warn_kkt=FALSE,
 			  max_iter=max_iter,
+			  kkt_stop=kkt_stop,
+			  parameter_stop=parameter_stop,
+			  objective_stop=objective_stop,
 			  kkt_tol=kkt_tol,
+			  parameter_tol=parameter_tol,
 			  objective_tol=objective_tol)
 
     if (warn_kkt && (!output$kkt_check)) {
@@ -372,15 +393,19 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
 debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
                          is_wide, 
                          row, 
-                         mu, 
+                         bound, 
 	                 linesearch=TRUE,     # do a linesearch?
-		         scaling_factor=1.2,  # multiplicative factor for linesearch
+		         scaling_factor=1.5,  # multiplicative factor for linesearch
 		         max_active=NULL,     # how big can active set get?
 			 max_try=10,          # how many steps in linesearch?
 			 warn_kkt=FALSE,      # warn if KKT does not seem to be satisfied?
-			 max_iter=100,        # how many iterations for each optimization problem
+			 max_iter=50,         # how many iterations for each optimization problem
+                         kkt_stop=TRUE,       # stop based on KKT conditions?
+                         parameter_stop=TRUE, # stop based on relative convergence of parameter?
+                         objective_stop=TRUE, # stop based on relative decrease in objective?
                          kkt_tol=1.e-4,       # tolerance for the KKT conditions
-			 objective_tol=1.e-8  # tolerance for relative decrease in objective
+			 parameter_tol=1.e-4, # tolerance for relative convergence of parameter
+			 objective_tol=1.e-4  # tolerance for relative decrease in objective
                          ) {
 
   p = ncol(Xinfo)
@@ -389,9 +414,11 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
       max_active = min(nrow(Xinfo), ncol(Xinfo))
   }
 
+   
   # Initialize variables 
 
   soln = rep(0, p)
+  soln = as.numeric(soln)
   ever_active = rep(0, p)
   ever_active[1] = row      # 1-based
   ever_active = as.integer(ever_active)
@@ -407,11 +434,15 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
 
   last_output = NULL
 
+  if (is_wide) {
+     Xsoln = as.numeric(rep(0, nrow(Xinfo)))
+  }
+
   while (counter_idx < max_try) {
 
       if (!is_wide) {
           result = solve_QP(Xinfo, # this is non-neg-def matrix
-                            mu, 
+                            bound, 
                             max_iter, 
                             soln, 
                             linear_func, 
@@ -420,11 +451,15 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
                             nactive, 
                             kkt_tol, 
                             objective_tol, 
-                            max_active) 
+			    parameter_tol,
+                            max_active,
+			    kkt_stop,
+			    objective_stop,
+			    parameter_stop)
       } else {
-          Xsoln = rep(0, nrow(Xinfo))
-          result = solve_QP_wide(Xinfo, # this is a design matrix
-                                 mu, 
+          result = solve_QP_wide(Xinfo,                      # this is a design matrix
+                                 as.numeric(rep(bound, p)),  # vector of Lagrange multipliers
+				 0,                          # ridge_term 
                                  max_iter, 
                                  soln, 
                                  linear_func, 
@@ -434,7 +469,11 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
                                  nactive, 
                                  kkt_tol, 
                                  objective_tol, 
-                                 max_active) 
+				 parameter_tol,
+                                 max_active,
+				 kkt_stop,
+				 objective_stop,	
+				 parameter_stop)
 
       }
 
@@ -458,13 +497,13 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
          if ((iter < (max_iter+1)) && (counter_idx > 1)) { 
            break;      # we've found a feasible point and solved the problem            
          }
-         mu = mu * scaling_factor;
+         bound = bound * scaling_factor;
       } else {         # trying to drop the bound parameter further
          if ((iter == (max_iter + 1)) && (counter_idx > 1)) {
             result = last_output; # problem seems infeasible because we didn't solve it
    	    break;                # so we revert to previously found solution
          }
-         mu = mu / scaling_factor;
+         bound = bound / scaling_factor;
       }
 
       # If the active set has grown to a certain size
@@ -490,7 +529,8 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
   } 
 
   return(list(soln=result$soln,
-              kkt_check=result$kkt_check))
+              kkt_check=result$kkt_check,
+	      gradient=result$gradient))
 
 }
 

selectiveInference/R/funs.fixed.R

@@ -29,7 +29,7 @@ if( sum(status==0)+sum(status==1)!=length(y)) stop("status vector must have valu
 vars=which(m)
 if(sum(m)>0){
     bhat=beta[beta!=0] #penalized coefs just for active variables
-    s2=sign(bhat)
+    sign_bhat=sign(bhat)
 
  #check KKT
 
@@ -40,7 +40,7 @@ if(sum(m)>0){
     res=residuals(aaa,type="score")
 if(!is.matrix(res)) res=matrix(res,ncol=1)
 scor=colSums(res)
-    g=(scor+lambda*s2)/(2*lambda)
+    g=(scor+lambda*sign_bhat)/(2*lambda)
 #    cat(c(g,lambda,tol.kkt),fill=T)
      if (any(abs(g) > 1+tol.kkt) )
     warning(paste("Solution beta does not satisfy the KKT conditions",
@@ -49,9 +49,9 @@ scor=colSums(res)
 # Hessian of partial likelihood at the LASSO solution    
 MM=vcov(aaa)
 
-bbar=(bhat+lambda*MM%*%s2)
-A1=-(mydiag(s2))
-b1= -(mydiag(s2)%*%MM)%*%s2*lambda
+bbar=(bhat+lambda*MM%*%sign_bhat)
+A1=-(mydiag(sign_bhat))
+b1= -(mydiag(sign_bhat)%*%MM)%*%sign_bhat*lambda
 
    temp=max(A1%*%bbar-b1)
 
@@ -63,7 +63,7 @@ b1= -(mydiag(s2)%*%MM)%*%s2*lambda
 # the one sided p-values are a bit off
 
     for(jj in 1:length(bbar)){
-      vj=rep(0,length(bbar));vj[jj]=s2[jj]
+      vj=rep(0,length(bbar));vj[jj]=sign_bhat[jj]
 
 
       junk=TG.pvalue(bbar, A1, b1, vj,MM)
@@ -73,7 +73,7 @@ b1= -(mydiag(s2)%*%MM)%*%s2*lambda
        vup[jj]=junk$vup
        sd[jj]=junk$sd
 
-      junk2=TG.interval(bbar, A1, b1, vj, MM, alpha, flip=(s2[jj]==-1))
+      junk2=TG.interval(bbar, A1, b1, vj, MM, alpha, flip=(sign_bhat[jj]==-1))
        ci[jj,]=junk2$int
        tailarea[jj,] = junk2$tailarea
 

selectiveInference/R/funs.fixedCox.R

@@ -32,7 +32,7 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
  m=beta[-1]!=0  #active set
 
     bhat=c(beta[1],beta[-1][beta[-1]!=0]) # intcpt plus active vars
-     s2=sign(bhat)
+     sign_bhat=sign(bhat)
      lam2m=diag(c(0,rep(lambda,sum(m))))
 
 
@@ -66,14 +66,14 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
  # MM=solve(t(xxm)%*%w%*%xxm)
    MM=solve(scale(t(xxm),F,1/ww)%*%xxm)
   gm = c(0,-g[vars]*lambda) # gradient at LASSO solution, first entry is 0 because intercept is unpenalized
-                            # at exact LASSO solution it should be s2[-1]
+                            # at exact LASSO solution it should be sign_bhat[-1]
   dbeta = MM %*% gm
 
-  # bbar=(bhat+lam2m%*%MM%*%s2)  # JT: this is wrong, shouldn't use sign of intercept anywhere...
+  # bbar=(bhat+lam2m%*%MM%*%sign_bhat)  # JT: this is wrong, shouldn't use sign of intercept anywhere...
   bbar = bhat - dbeta
 
-  A1=-(mydiag(s2))[-1,]
-  b1= (s2 * dbeta)[-1]
+  A1=-(mydiag(sign_bhat))[-1,]
+  b1= (sign_bhat * dbeta)[-1]
 
   tol.poly = 0.01 
   if (max((A1 %*% bbar) - b1) > tol.poly)
@@ -87,7 +87,7 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
 
 
     for(jj in 1:sum(m)){
-       vj=c(rep(0,sum(m)+1));vj[jj+1]=s2[jj+1]
+       vj=c(rep(0,sum(m)+1));vj[jj+1]=sign_bhat[jj+1]
       # compute p-values
       junk=TG.pvalue(bbar, A1, b1, vj, MM)
       pv[jj] = junk$pv
@@ -96,7 +96,7 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
    vup[jj]=junk$vup
        sd[jj]=junk$sd
 
-     junk2=TG.interval(bbar, A1, b1, vj, MM,alpha=alpha, flip=(s2[jj+1]==-1))
+     junk2=TG.interval(bbar, A1, b1, vj, MM,alpha=alpha, flip=(sign_bhat[jj+1]==-1))
 
      ci[jj,]=junk2$int
      tailarea[jj,] = junk2$tailarea