version 1.0.1

DESCRIPTION

@@ -1,14 +1,14 @@
 Package: qut
 Type: Package
 Title: Quantile Universal Threshold
-Version: 1.0
+Version: 1.0.1
 Author: Jairo Diaz, Sylvain Sardy, Caroline Giacobino, Nick Hengartner.
 Maintainer: Jairo Diaz <jairo.diaz@unige.ch>
 Description: Selection of a threshold parameter for GLM-lasso to obtain a sparse model 
  with a good compromise between high true positive rate and low false discovery rate.
 License: GPL-2
 Depends: Matrix, glmnet, lars
 NeedsCompilation: no
-Packaged: 2015-09-20 07:54:10 UTC; Jairo
+Packaged: 2015-11-22 20:28:55 UTC; Jairo
 Repository: CRAN
-Date/Publication: 2015-09-20 15:41:39
+Date/Publication: 2015-11-23 12:20:39

MD5

@@ -1,18 +1,19 @@
-796386f5b46abf28de97ed5e6bfa948a *DESCRIPTION
-7e939e6767941731e8bc9e537d7ebf75 *NAMESPACE
-b35d2c3fe626c987a1236bf500ff89a1 *R/lambdaqut.R
+e4c9210192d1c05af27334a477b198f5 *DESCRIPTION
+4f9cf21e088bc66df0acd208ededb1df *NAMESPACE
+6ec7a07d1de47c081c07680ff189c26c *R/coef.qut.R
+788f694e8b29a08e21e1b634a51c61c2 *R/lambdaqut.R
 c5c0e9083a2d2266484c69cec2892cfe *R/predict.qut.R
-60b8c469d26b0a0a513bad3854bf3aca *R/qut.R
+dacf6caafe796d2c4ab4f614ab8ad9e1 *R/qut.R
 07d830a8b07ef1d53e36e1e72b24b789 *R/sigmaqut.R
 83c9097fad099cce1c20298145c00225 *R/sigmarcv.R
 6366e619fc2df2a00b86d383beaf2eaf *data/datalist
 fb4808ca9012f6abbb606d7d1b21a039 *data/leukemia.rda
 3f8a1ef02e84a88f448620114f857d20 *data/riboflavin.rda
-82c882af4ff01af36d89fc694e28510f *man/lambdaqut.Rd
+a8447c6ed4b177381723ab76d60774c8 *man/lambdaqut.Rd
 ed49b029b582232aab02fc804364a362 *man/leukemia.Rd
-2fc547bfe76a865d91d34c6954aff496 *man/predict.qut.Rd
+a55c86328f997036f9dd03d5fe418434 *man/predict.qut.Rd
 b52ed035629cdab2584ed1bd53bd9541 *man/qut-package.Rd
-7701f3ae6406752095a020ef17f50a08 *man/qut.Rd
+75face335c37e6495e57e559df198f78 *man/qut.Rd
 e1e2b3e8a62002755cc157be3174b951 *man/riboflavin.Rd
 5cb8388ddf7ecb6984c3ac7afb91319a *man/sigmaqut.Rd
 504c41ff245b7e6c27cd3a1eb75b6a43 *man/sigmarcv.Rd

NAMESPACE

@@ -1,7 +1,8 @@
-export(qut,lambdaqut,sigmaqut,sigmarcv,predict.qut)
+export(qut,lambdaqut,sigmaqut,sigmarcv,predict.qut,coef.qut)
 import(Matrix, glmnet, lars)
 importFrom("stats", "coef", "gaussian", "glm", "glm.fit", "lm",
              "median", "predict", "quantile", "rbinom", "rnorm", "rpois",
              "var")
 
-S3method(predict, qut)
+S3method(predict, qut)
+S3method(coef, qut)

R/coef.qut.R

@@ -0,0 +1,11 @@
+coef.qut <- 
+function(object, mode='glm',...){
+	if(mode=='lasso') return(object$beta)
+	else{
+		if(!attr(object$betaglm,'converged')){
+			warning('GLM did not converged, returning LASSO coefficients')
+			return(object$lasso)
+		}
+		else return(object$betaglm)
+	}
+}

R/lambdaqut.R

@@ -1,5 +1,5 @@
 lambdaqut <-
-function(y,X,family=gaussian,alpha.level='default',M=1000,qut.standardize=TRUE,intercept=TRUE,no.penalty=NULL,offset=NULL,bootstrap=TRUE){
+function(y,X,family=gaussian,alpha.level='default',M=1000,qut.standardize=TRUE,intercept=TRUE,no.penalty=NULL,offset=NULL,bootstrap=TRUE,beta0=NA){
 #FUNCTION TO OBTAIN THE QUANTILE UNIVERSAL THRESHOLD
 
 	#Hidden function for vectorized bootstrapping
@@ -60,32 +60,34 @@ function(y,X,family=gaussian,alpha.level='default',M=1000,qut.standardize=TRUE,i
 		if(intercept) A=cBind(rep(1,n),A)  #if there is an intercept (column of ones)
 
 		#Estimate beta0 as glm(y~A)
-		beta0=glm.fit(y=y,x=as.matrix(A),intercept=FALSE,family=family,offset=offset)$coef	
+		if(family$family=='gaussian') beta0=rep(0,ncol(A))
+		else if(!is.numeric(beta0)) beta0=glm.fit(y=y,x=as.matrix(A),intercept=FALSE,family=family,offset=offset)$coef
+		else if(length(beta0)!=(length(no.penalty)+intercept)) stop("length of beta0 is different from the number of unpenalized covariates or the intercept has not been included")
 		muhat=family$linkinv(as.matrix(A)%*%beta0+O)
 	}
 
 	#Set default alpha.level
 	if(alpha.level=='default')	alpha.level=1/(sqrt(pi*log(p)))
 
-	#adjustment gamma to the quantile (existence of glm)
-	gamma=1
-	if(intercept&is.null(no.penalty)){	
-		if(family$family=='poisson') gamma=1-exp(-n*muhat[1])
-		else if(family$family=='binomial') gamma=1-muhat[1]^n+(1-muhat[1])^n
-	}
-	else warning("Adjustment gamma for glm existence is not considered")
+	#Check if A=1 and there is an explicit characterization of D
+	if(!(intercept&is.null(no.penalty)&sum(O==0)==n)) warning("Explicit characterization of D is not defined")
 
 	#Monte Carlo Simulation of the null model
+	znotinD=0
 	if(family$family=='gaussian') z=matrix(rnorm(n*M,mean=as.vector(muhat),sd=1),n,M)
 	else if(family$family=='poisson'){
 		z=matrix(rpois(n*M,lambda=as.vector(muhat)),n,M)
+		znotinD=sum(apply(z,2,sum)==0)
 		z=z[,apply(z,2,sum)!=0]
 	}	
 	else if(family$family=='binomial'){
 		z=matrix(rbinom(n*M,size=1,prob=as.vector(muhat)),n,M)
+		znotinD=sum((apply(z,2,sum)==0)|(apply(z,2,sum)==n))
 		z=z[,(apply(z,2,sum)!=0)&(apply(z,2,sum)!=n)]
 	}
 	else stop("Not available family")
+
+	if(znotinD>(M-2)) stop("Can't generate valid simulations under the null hypothesis, try increasing M")
 
 	#obtain estimate of beta0 for each simulation
 	if(!intercept&is.null(no.penalty)){
@@ -108,7 +110,7 @@ function(y,X,family=gaussian,alpha.level='default',M=1000,qut.standardize=TRUE,i
 
 		#quantile-based standardization
 		if(qut.standardize){
-			divX=apply(abs(bp),1,quantile,prob=(1-alpha.level)/gamma)
+			divX=apply(abs(bp),1,quantile,prob=(1-alpha.level))
 			divX[which(divX==0)]=1
 
 			#Alignment/scaling
@@ -125,7 +127,7 @@ function(y,X,family=gaussian,alpha.level='default',M=1000,qut.standardize=TRUE,i
 		if(qut.standardize){
 			bp=abs(t(X)%*%(z-muhatZ))
 			scale.factor=rep(1,p)
-			divX=apply(abs(bp),1,quantile,prob=(1-alpha.level)/gamma)
+			divX=apply(abs(bp),1,quantile,prob=(1-alpha.level))
 			divX[which(divX==0)]=1
 
 			#Alignment/scaling
@@ -134,15 +136,13 @@ function(y,X,family=gaussian,alpha.level='default',M=1000,qut.standardize=TRUE,i
 
 		#Bootstrapping
 		bp=apply(rbind(muhatZ,z),2,qutbootstrap)
-		resultsMC=apply(bp,2,max,na.rm=TRUE)
-		lambda=quantile(resultsMC, prob=(1-alpha.level)/gamma)
 	}
 
 	#Obtain distribution of Lambda=max(X'z)
 	resultsMC=apply(bp,2,max,na.rm=TRUE)
-
+	resultsMC=c(resultsMC,rep(Inf,znotinD))
 	#obtain lambda.qut for the given quantile
-	lambda=quantile(resultsMC, prob=(1-alpha.level)/gamma)
+	lambda=quantile(resultsMC, prob=(1-alpha.level))
 
 	lambdamax=max(abs(t(X)%*%(y-muhat)))