rob

Author: tibshirani

selectiveInference/DESCRIPTION

@@ -1,8 +1,8 @@
 Package: selectiveInference
 Type: Package
 Title: Tools for Post-Selection Inference
-Version: 1.1.3
-Date: 2016-02-8
+Version: 1.2.1
+Date: 2016-07-3
 Author: Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor,
     Joshua Loftus, Stephen Reid
 Maintainer: Rob Tibshirani <[email protected]>
@@ -14,5 +14,5 @@ Suggests:
     Rmpfr
 Description: New tools for post-selection inference, for use
     with forward stepwise regression, least angle regression, the
-    lasso, and the many means problem.
+    lasso, and the many means problem. The lasso function implements  Gaussian,  logistic and Cox survival models.
 License: GPL-2

selectiveInference/NAMESPACE

@@ -40,4 +40,6 @@ importFrom("graphics", abline, axis, matplot)
 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")
+
 

selectiveInference/R/funs.fixed.R

@@ -41,7 +41,7 @@ else{
   n = nrow(x)
   p = ncol(x)
   beta = as.numeric(beta)
-  if (length(beta) != p) stop("beta must have length equal to ncol(x)")
+  if (length(beta) != p) stop("Since family='gaussian', beta must have length equal to ncol(x)")
 
   # If glmnet was run with an intercept term, center x and y
   if (intercept==TRUE) {

selectiveInference/R/funs.fixedLogit.R

@@ -18,7 +18,7 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
      n=length(y)
      p=ncol(x)
  # I assume that intcpt was used
-     if(length(beta)!=p+1) stop("beta must be of length ncol(x)+1, that is, it should include an intercept")
+     if(length(beta)!=p+1) stop("Since family='binomial', beta must be of length ncol(x)+1, that is, it should include an intercept")
      nvar=sum(beta[-1]!=0)
      pv=vlo=vup=sd=rep(NA, nvar)
      ci=tailarea=matrix(NA,nvar,2)
@@ -41,11 +41,12 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
    etahat = xxm %*% bhat
    prhat = as.vector(exp(etahat) / (1 + exp(etahat)))
    ww=prhat*(1-prhat)
-   w=diag(ww)
+ #  w=diag(ww)
 
 #check KKT
    z=etahat+(y-prhat)/ww
-       g=  t(x)%*%w%*%(z-etahat)/lambda # negative gradient scaled by lambda
+ #      g=  t(x)%*%w%*%(z-etahat)/lambda # negative gradient scaled by lambda
+  gg=scale(t(x),FALSE,1/ww)%*%(z-etahat)/lambda # negative gradient scaled by lambda
      if (any(abs(g) > 1+tol.kkt) )
     warning(paste("Solution beta does not satisfy the KKT conditions",
                   "(to within specified tolerances)"))
@@ -62,7 +63,8 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
                   "'thresh' parameter, for a more accurate convergence."))
 
  #constraints for active variables             
-  MM=solve(t(xxm)%*%w%*%xxm)
+ # 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]
   dbeta = MM %*% gm

tests/test.fixed.R

@@ -510,3 +510,47 @@ g=t(X)%*%(Y-X%*%b)/lam2
 out = fixedLassoInf(X,Y,beta,lam*n)
 
 
+
+
+
+#
+
+#gaussian
+n=50
+p=10
+sigma=.7
+beta=c(0,0,0,0,rep(0,p-4))
+set.seed(43)
+nsim = 1000
+pvals <- matrix(NA, nrow=nsim, ncol=p)
+x = matrix(rnorm(n*p),n,p)
+x = scale(x,T,T)/sqrt(n-1)
+mu = x%*%beta
+for (i in 1:nsim) {
+    cat(i)
+y=mu+sigma*rnorm(n)
+#y=y-mean(y)
+# first run  glmnet
+    pf=c(rep(.001,4),rep(1,p-4))
+     xs=scale(x,FALSE,pf) #scale cols of x by penalty factors
+     # first run glmnet
+     gfit = glmnet(xs,y,standardize=FALSE)
+     
+    
+     lambda = .8
+     beta = coef(gfit, s=lambda/n, exact=TRUE)[-1]
+     
+     # compute fixed lambda p-values and selection intervals
+     aa = fixedLassoInf(xs,y,beta,lambda,sigma=sigma)
+     
+pvals[i, which(beta != 0)] <- aa$pv
+}
+nulls = 1:nsim
+np = pvals[nulls,-(1:4)]
+mean(np[!is.na(np)] < 0.1)
+o=!is.na(np)
+plot((1:sum(o))/sum(o),sort(np))
+abline(0,1)
+#####
+
+