require(MASS)
require(glmnet)
require(arm)
options(warn=-1)

plasso_logistic = function(x,y) {

# This function implements the (adaptive) pLasso for logistic regression 
# using the Gelman et al. (2008) prior and the plug-in method 
# for estimating the predictions \bar y_i^f.
# Inputs:
# 	x - the desgin matrix (without the column of 1s)
# 	y - the binary response vector
# Outputs:
#	a vector of pLasso estimate of coefficients
# Reference: "The predictive Lasso" by Tran, Nott and Leng

    
n = dim(x)[1]
p = dim(x)[2]
M2 <- bayesglm (y ~ x, family=binomial(link="logit"))
beta = M2$coef
xx = cbind(as.matrix(rep(1,n)),x)
yf = exp(xx%*%beta)/(1+exp(xx%*%beta))
beta.plasso = aglm(x,as.numeric(yf))
return(beta.plasso)

}

nplasso_logistic = function(x,y) {

# This function implements the non-adaptive pLasso for logistic regression 
    
n = dim(x)[1]
p = dim(x)[2]
M2 <- bayesglm (y ~ x, family=binomial(link="logit"))
beta = M2$coef
xx = cbind(as.matrix(rep(1,n)),x)
yf = exp(xx%*%beta)/(1+exp(xx%*%beta))
beta.nplasso = lasso_plasso(x,as.numeric(yf))
return(beta.nplasso)

}

lasso_plasso = function(x,y)
# helper function for non-adaptive plasso for logistic regression
{
n = dim(x)[1]
p = dim(x)[2]
Y = cbind(1-y,y)
fit = glmnet(x, Y, weights = rep(1,length(y)),family='binomial',standardize = FALSE)
cv = cv.glmnet(x,y)
coef = coef(fit,cv$lambda.min)
return(as.numeric(coef))
}



standardize1 = function(x)
{
xx = matrix(rep(0,prod(dim(x))),dim(x)[1])
for (i in 1:dim(x)[2])
{
	xx[,i] = 0.5*(x[,i]-mean(x[,i]))/sqrt(var(x[,i]))
}
return(xx)
}

aglm = function(x,y)
# helper function for adaptive pLasso for logistic regression
{
n = dim(x)[1]
p = dim(x)[2]
Y = cbind(1-y,y)
if (n>p)
	{ fit  = glm(Y ~ x, family = 'binomial',weights = rep(1,length(y)))
	  wei = abs(fit$coefficients)
} else
	{ 
	  fit = glmnet(x, Y, weights = rep(1,length(y)),family='binomial',standardize = FALSE)
	  cv = cv.glmnet(x,y)
        wei = abs(coef(fit,cv$lambda.min))
}
	Wei = diag(wei[2:length(wei)])
	X = x%*%Wei
	fit = glmnet(X,Y,family='binomial',standardize = FALSE,weights = rep(1,length(y)))
	cv = cv.glmnet(X,y)
	coef = coef(fit,cv$lambda.min)
	coef = c(coef[1],Wei%*%coef[2:length(coef)])
}

alasso_logistic = function(x,y)
# adaptive Lasso for logistic regression
{
n = dim(x)[1]
p = dim(x)[2]
if (n>p)
	{ fit  = glm(y ~ x, family = 'binomial')
	  wei = abs(fit$coefficients)
} else
	{ 
	  fit = glmnet(x, y,family='binomial',standardize = FALSE)
	  cv = cv.glmnet(x,y)
        wei = abs(coef(fit,cv$lambda.min))
}
	Wei = diag(wei[2:length(wei)])
	X = x%*%Wei
	fit = glmnet(X,y,family='binomial',standardize = FALSE)
	cv = cv.glmnet(X,y)
	coef = coef(fit,cv$lambda.min)
	coef = c(coef[1],Wei%*%coef[2:length(coef)])
}


lasso_logistic = function(x,y)
# Lasso for logistic regression
{
n = dim(x)[1]
p = dim(x)[2]

fit = glmnet(x,y,family='binomial',standardize = FALSE)
cv = cv.glmnet(x,y)
coef = coef(fit,cv$lambda.min)
coef = c(coef[1],coef[2:length(coef)])
}

