Merge pull request #30 from jonathan-taylor/master

Updating to new glmnet `coef`

Author: jonathan-taylor

.travis.yml

@@ -11,6 +11,6 @@ addons:
 warnings_are_errors: true
 before_install:
   - tlmgr install index # for texlive and vignette?
-  - R -e 'install.packages("Rcpp", repos="http://cloud.r-project.org")'
+  - R -e 'install.packages(c("Rcpp", "intervals"), repos="http://cloud.r-project.org")'
   - make Rcpp
   - cd selectiveInference

selectiveInference/R/funs.fixed.R

@@ -2,10 +2,13 @@
 # for the solution of
 # min 1/2 || y - \beta_0 - X \beta ||_2^2 + \lambda || \beta ||_1
 
-fixedLassoInf <- function(x, y, beta, lambda, family=c("gaussian","binomial","cox"),intercept=TRUE, add.targets=NULL, status=NULL,
+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,
-                          gridrange=c(-100,100), bits=NULL, verbose=FALSE, linesearch.try=10) {
+                          gridrange=c(-100,100), bits=NULL, verbose=FALSE, 
+                          linesearch.try=10) {
 
   family = match.arg(family)
   this.call = match.call()

selectiveInference/R/funs.fixedCox.R

@@ -1,6 +1,10 @@
-fixedCoxLassoInf=function(x,y,status,beta,lambda,alpha=.1, type=c("partial"),tol.beta=1e-5, tol.kkt=0.1,
-                     gridrange=c(-100,100), bits=NULL, verbose=FALSE,this.call=NULL){
-
+fixedCoxLassoInf=function(x, y, status,
+                          beta, lambda,
+                          alpha=.1, type=c("partial"),
+                          tol.beta=1e-5, tol.kkt=0.1,
+                          gridrange=c(-100,100), 
+                          bits=NULL, verbose=FALSE,
+                          this.call=NULL){
 
  checkargs.xy(x,y)
  if(is.null(status)) stop("Must supply `status' argument")

selectiveInference/man/fixedLassoInf.Rd

@@ -160,7 +160,7 @@ gfit = glmnet(x,y,standardize=FALSE)
 # extract coef for a given lambda; note the 1/n factor!
 # (and we don't save the intercept term)
 lambda = .8
-beta = coef(gfit, s=lambda/n, exact=TRUE)[-1]
+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)
@@ -170,8 +170,8 @@ out
 ## as above, but use lar function instead to get initial 
 ## lasso fit (should get same results)
 lfit = lar(x,y,normalize=FALSE)
-beta = coef(lfit,s=lambda,mode="lambda")
-out2 = fixedLassoInf(x,y,beta,lambda,sigma=sigma)
+beta = coef(lfit, s=lambda, mode="lambda")
+out2 = fixedLassoInf(x, y, beta, lambda, sigma=sigma)
 out2
 
 ## mimic different penalty factors by first scaling x
@@ -189,12 +189,12 @@ pf=c(rep(1,7),rep(.1,3))  #define penalty factors
 pf=p*pf/sum(pf)   # penalty factors should be rescaled so they sum to p
 xs=scale(x,FALSE,pf) #scale cols of x by penalty factors
 # first run glmnet
-gfit = glmnet(xs,y,standardize=FALSE)
+gfit = glmnet(xs, y, standardize=FALSE)
 
 # extract coef for a given lambda; note the 1/n factor!
 # (and we don't save the intercept term)
 lambda = .8
-beta_hat = coef(gfit, s=lambda/n, exact=TRUE)[-1]
+beta_hat = coef(gfit, x=xs, y=y, s=lambda/n, exact=TRUE)[-1]
 
 # compute fixed lambda p-values and selection intervals
 out = fixedLassoInf(xs,y,beta_hat,lambda,sigma=sigma)
@@ -221,37 +221,40 @@ set.seed(43)
      # extract coef for a given lambda; note the 1/n factor!
      # (and here  we DO  include the intercept term)
      lambda = .8
-     beta_hat = coef(gfit, s=lambda/n, exact=TRUE)
+     beta_hat = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)
 
      # compute fixed lambda p-values and selection intervals
      out = fixedLassoInf(x,y,beta_hat,lambda,family="binomial")
      out
 
-#Cox model
-set.seed(43)
+
+     # Cox model
+     
+     set.seed(43)
      n = 50
      p = 10
      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))
      tim = as.vector(x\%*\%beta + sigma*rnorm(n))
      tim= tim-min(tim)+1
-status=sample(c(0,1),size=n,replace=TRUE)
+     status=sample(c(0,1),size=n,replace=TRUE)
      # first run glmnet
 
 
-     gfit = glmnet(x,Surv(tim,status),standardize=FALSE,family="cox")
+     y = Surv(tim,status)
+     gfit = glmnet(x, y, standardize=FALSE, family="cox")
 
      # extract coef for a given lambda; note the 1/n factor!
 
      lambda = 1.5
-     beta_hat = as.numeric(coef(gfit, s=lambda/n, exact=TRUE))
-     
+     beta_hat = as.numeric(coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE))
+
      # compute fixed lambda p-values and selection intervals
-     out = fixedLassoInf(x,tim,beta_hat,lambda,status=status,family="cox")
+     out = fixedLassoInf(x, tim, beta_hat, lambda, status=status, family="cox")
      out
 
 # Debiased lasso or "full"
@@ -272,7 +275,7 @@ 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, s=lambda/n, exact=TRUE)[-1]
+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)