Sample_linear.R

  #############################################################################
### Sample code to get LASSO, p, pLASS, LASSO-A, p-A, pLASSO-A estimators ###
#############################################################################

library(glmnet)
library(MASS)
source("function_linear.R")

###====== Read x and y ======###
x <- read.table("x_linear.txt", header = FALSE)
y <- read.table("y_linear.txt", header = FALSE)
x <- as.matrix(x)
y <- as.matrix(y)
n <- nrow(x)
p <- ncol(x)
        
###====== Settings ======###    
# prior set
pos.prior <- c(1 : 10)  

# refit or not?
is.refit <- 1

# sequence of tuning parameters
eta.seq <- c(seq(10, 1, by = -2), seq(1, 0.2, by = -0.2))
tau.seq <- c(2, 1, 0.5)

# 3-fold cross validation grouping
cv.group <- c(rep(1 : 3, each = floor(n/3)), rep(3, n - 3 * floor(n/3)))

###====== LASSO ======###
# penalty.factor: set the nonpenalized terms to 0
penalty.factor <- rep(1, p)
beta.1 <- find.beta(x, y, 0, 0, rep(0, p + 1), penalty.factor, cv.group, is.refit)

###====== prior ======###
penalty.factor <- rep(1, p)
penalty.factor[pos.prior] <- 0
beta.2 <- find.beta(x, y, 0, 0, rep(0, p + 1), penalty.factor, cv.group, is.refit)
        
###====== pLASSO ======###
penalty.factor <- rep(1, p)
beta.3 <- find.beta(x, y, eta.seq, 0, beta.2, penalty.factor, cv.group, is.refit)

###====== LASSO-A ======###
penalty.factor <- rep(1, p)
beta.4 <- find.beta(x, y, 0, tau.seq, beta.1, penalty.factor, cv.group, is.refit)

###====== p-A ======###
penalty.factor <- rep(1, p)
beta.5 <- find.beta(x, y, 0, tau.seq, beta.2, penalty.factor, cv.group, is.refit)

###====== pLASSO-A ======###
penalty.factor <- rep(1, p)
beta.6 <- find.beta(x, y, 0, tau.seq, beta.3, penalty.factor, cv.group, is.refit)

cbind(beta.1, beta.2, beta.3, beta.4, beta.5, beta.6)