Merge pull request #33 from jonathan-taylor/master

Fixing one-sided pvalues in "full"

Author: tibshirani

selectiveInference/R/funs.fixed.R

@@ -82,17 +82,20 @@ fixedLassoInf <- function(x, y, beta,
 
     tol.coef = tol.beta * sqrt(n / colSums(x^2))
     vars = which(abs(beta) > tol.coef)
+    sign_vars = sign(beta[vars])
 
     if(sum(vars)==0){
       cat("Empty model",fill=T)
       return()
     }
-    if (any(sign(g[vars]) != sign(beta[vars])))
+
+    if (any(sign(g[vars]) != sign_vars)) {
       warning(paste("Solution beta does not satisfy the KKT conditions",
                     "(to within specified tolerances). You might try rerunning",
                     "glmnet with a lower setting of the",
                     "'thresh' parameter, for a more accurate convergence."))
-    
+    }
+
     # Get lasso polyhedral region, of form Gy >= u
 
     logical.vars=rep(FALSE,p)
@@ -132,13 +135,19 @@ fixedLassoInf <- function(x, y, beta,
     }
 
     # add additional targets for inference if provided
-    if (!is.null(add.targets)) vars = sort(unique(c(vars,add.targets,recursive=T)))
-    
-      k = length(vars)
+    if (!is.null(add.targets)) {
+       # vars is boolean...
+       old_vars = vars & TRUE
+       vars[add.targets] = TRUE
+       sign_vars = sign(beta[vars]) 
+       sign_vars[!old_vars] = NA
+       stop("`add.targets` not fully implemented yet")
+    }
+
+    k = length(vars)
     pv = vlo = vup = numeric(k)
     vmat = matrix(0,k,n)
     ci = tailarea = matrix(0,k,2)
-    sign = numeric(k)
 
     if (type=="full" & p > n) {
       if (intercept == T) {
@@ -202,28 +211,36 @@ fixedLassoInf <- function(x, y, beta,
     vj = M[j,]
     mj = sqrt(sum(vj^2))
     vj = vj / mj        # Standardize (divide by norm of vj)
-    sign[j] = sign(sum(vj*y))
-    vj = sign[j] * vj
+
+    if (!is.na(sign_vars[j])) {
+        vj = sign_vars[j] * vj
+    }
 
     limits.info = TG.limits(y, A, b, vj, Sigma=diag(rep(sigma^2, n)))
     a = TG.pvalue.base(limits.info, null_value=null_value[j], bits=bits)
     pv[j] = a$pv
+    if (is.na(sign_vars[j])) { # for variables not in the active set, report 2-sided pvalue
+       pv[j] = 2 * min(pv[j], 1 - pv[j])
+    }
     vlo[j] = a$vlo * mj # Unstandardize (mult by norm of vj)
     vup[j] = a$vup * mj # Unstandardize (mult by norm of vj)
-    vmat[j,] = vj * mj * sign[j]  # Unstandardize (mult by norm of vj)
-
+    if (!is.na(sign_vars[j])) { 
+        vmat[j,] = vj * mj * sign_vars[j]  # Unstandardize (mult by norm of vj) and fix sign
+    } else {
+        vmat[j,] = vj * mj # Unstandardize (mult by norm of vj)
+    }
     a = TG.interval.base(limits.info, 
                          alpha=alpha,
                          gridrange=gridrange,
-			 flip=(sign[j]==-1),
+			 flip=(sign_vars[j]==-1),
                          bits=bits)
     ci[j,] = (a$int-null_value[j]) * mj # Unstandardize (mult by norm of vj)
     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,sigma=sigma,alpha=alpha,
+    vars=vars,sign=sign_vars,sigma=sigma,alpha=alpha,
     sd=sigma*sqrt(rowSums(vmat^2)),
     coef0=vmat%*%y,
     call=this.call)

selectiveInference/man/fixedLassoInf.Rd

@@ -257,29 +257,54 @@ set.seed(43)
      out = fixedLassoInf(x, tim, beta_hat, lambda, status=status, family="cox")
      out
 
-# Debiased lasso or "full"
+     # Debiased lasso or "full"
 
-n = 50
-p = 100
-sigma = 1
+     n = 50
+     p = 100
+     sigma = 1
 
-x = matrix(rnorm(n*p),n,p)
-x = scale(x,TRUE,TRUE)
+     x = matrix(rnorm(n*p),n,p)
+     x = scale(x,TRUE,TRUE)
 
-beta = c(3,2,rep(0,p-2))
-y = x\%*\%beta + sigma*rnorm(n)
+     beta = c(3,2,rep(0,p-2))
+     y = x\%*\%beta + sigma*rnorm(n)
 
-# first run glmnet
-gfit = glmnet(x, y, standardize=FALSE, intercept=FALSE)
+     # first run glmnet
+     gfit = glmnet(x, y, standardize=FALSE, intercept=FALSE)
 
-# extract coef for a given lambda; note the 1/n factor!
-# (and we don't save the intercept term)
-lambda = 2.8
-beta = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)[-1]
+     # extract coef for a given lambda; note the 1/n factor!
+     # (and we don't save the intercept term)
+     lambda = 2.8
+     beta = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)[-1]
 
-# compute fixed lambda p-values and selection intervals
-out = fixedLassoInf(x, y, beta, lambda, sigma=sigma, type='full', intercept=FALSE)
-out
+     # compute fixed lambda p-values and selection intervals
+     out = fixedLassoInf(x, y, beta, lambda, sigma=sigma, type='full', intercept=FALSE)
+     out
+
+     # When n > p and "full" we use the full inverse
+     # instead of Javanmard and Montanari's approximate inverse
+
+     n = 200
+     p = 50
+     sigma = 1
+
+     x = matrix(rnorm(n*p),n,p)
+     x = scale(x,TRUE,TRUE)
+
+     beta = c(3,2,rep(0,p-2))
+     y = x\%*\%beta + sigma*rnorm(n)
+
+     # first run glmnet
+     gfit = glmnet(x, y, standardize=FALSE, intercept=FALSE)
+
+     # extract coef for a given lambda; note the 1/n factor!
+     # (and we don't save the intercept term)
+     lambda = 2.8
+     beta = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)[-1]
+
+     # compute fixed lambda p-values and selection intervals
+     out = fixedLassoInf(x, y, beta, lambda, sigma=sigma, type='full', intercept=FALSE)
+     out
 
 }