library(glmnet) #################################################################### #################################################################### #################################################################### auc.lasso=function(x, y, error0=1e-5, df0=10, ratio=1.08, print=F){ p=length(x[1,]) n0=sum(y==0) n1=sum(y==1) deltax=matrix(0, n1*n0, p) for(j in 1:p) {score0=x[y==0,j] score1=x[y==1,j] deltax[,j]=rep(score1, rep(n0 ,n1))-rep(score0, n1) } n=n1*n0 x=deltax beta=rep(0, p) scorehat=rep(0, n) prob=1/(1+exp(scorehat)) slope.dir=rep(0, p) for(j in 1:p) {slope.dir[j]=abs(mean(x[,j]*prob))} lambda.max=max(slope.dir) temp=ceiling(8/log(ratio)) lambda.grd=exp(log(lambda.max)-(1:temp)*log(ratio)) df=llike=gammar=betar=NULL max.df=0 nlambda=length(lambda.grd) count=0 while(count < nlambda && max.df1) {prob=pmin(1/(1+exp(scorehat)), 1-1e-8) slope.dir=rep(0, p) for(j in 1:p) {slope.dir[j]=abs(mean(x[,j]*prob))} lambda.bound=2*lambda0-lambda.grd[count-1] } id.sel=(1:p)[slope.dir>=lambda.bound] p.sel=length(id.sel) x.sel=x[,id.sel,drop=F] beta.old=beta[id.sel] beta.old=rnorm(p) L.old=-mean(log(1-prob))+lambda0*sum(abs(beta.old)) error=1 while(error>error0) {prob=pmin(1/(1+exp(scorehat)), 1-1e-8) ynew=scorehat+1/(1-prob) weightnew=prob*(1-prob) fit=cd.naiveupdate(x.sel, ynew, weight=weightnew, lambda=lambda0, beta=beta.old, error.tol=1e-8, trace=F) beta.old=fit$beta scorehat=as.vector(x.sel%*%beta.old) prob=pmin(1/(1+exp(scorehat)), 1-1e-8) L.new=-mean(log(1-prob))+lambda0*sum(abs(beta.old)) error=(L.old-L.new)/abs(L.old) if(print==T) {cat("\n") print(error) cat("\n") } L.old=L.new } beta=rep(0, p) beta[id.sel]=beta.old betar=rbind(betar, beta) max.df=sum(beta!=0) df=c(df, max.df) llike=c(llike, -mean(log(1-prob))) if(print==T) {print(c(count, p.sel, max.df))} } return(list(beta=betar, loglike=llike, lambda=lambda.grd, df=df)) } ############################################################################################## sh=function(x, lambda) {sign(x)*(abs(x)-lambda)*(abs(x)>lambda)} cd.naiveupdate=function(x, y, weight, lambda, beta, error.tol=1e-12, trace=F) {p=length(x[1,]) n=length(x[,1]) x=x*sqrt(weight) y=y*sqrt(weight) sigmax=rep(0, p) for(j in 1:p) sigmax[j]=mean(x[,j]^2) beta.old=beta=rnorm(p) yhat=as.vector(x%*%beta) res=y-yhat L.old=0.5*mean(res^2)+lambda*sum(abs(beta)) error=1 while(error>error.tol) {for(j in 1:p) {betanew=(sh(mean(res*x[,j])+sigmax[j]*beta[j], lambda))/sigmax[j] if(betanew!= beta[j]) {res=res-(betanew-beta[j])*x[,j] beta[j]=betanew } } L.new=0.5*mean(res^2)+lambda*sum(abs(beta)) error=(L.old-L.new)/L.old L.old=L.new if(trace==T) {print(error)} } return(list(beta=beta)) }