Function Linear.R

  ###################################################
### Refit the model with the selected variables ###
###################################################

### Inputs:
        # x: design matrix
        # y: response vector
        # beta: original coeffecients
### Output:
        # updated coefficients

update.beta <- function(x, y, beta) {
        pos <- which(abs(beta[-1]) > 1e-10)
        beta.refit <- rep(0, length(beta))
        if (length(pos) >= 1) {
                refit <- lm(y ~ x[, pos])
                beta.refit[c(1, pos + 1)] <- refit$coef
                return(beta.refit)
        }
        else {
                # refit <- lm(y ~ 1)
                beta.refit[1] <- mean(y)
                return(beta.refit)
    }
}

########################################
### Calculate mean squared residuals ###
########################################

### Inputs:
        # as in update.beta
### Output:
        # mean squared residuals

MSR.function <- function(x, y, beta) {
        eta <- cbind(rep(1, nrow(x)), x) %*% beta
        res <- y - eta
        return(sum(res * res)/nrow(x))
}

#######################################################
### pLASSO in linear regression by cross validation ###
#######################################################

### This is the main function for linear regression with cross validation.
### It can find each unweighted penalization estimators: LASSO, p, pLASSO.
### It can find each weighted penalization estimators: LASSO-A, p-A, pLASSO-A.
### For details of the above terms, refer to the paper.

### Inputs:
        # x: design matrix
        # y: response vector
        # eta.seq: the sequence of eta, the parameter tuning the prior information [(6) in the paper]
                # for LASSO/p/pLASSO, use eta sequence
                # for LASSO-A/p-A/pLASSO-A, use 0
        # tau.seq: the sequence of tau, the parameter tuning the L1 penalty weights [(35) in the paper]
                # for LASSO/p/pLASSO, use 0
                # for LASSO-A/p-A/pLASSO-A, use tau sequence
        # beta.prior: either the prior estimator or the initial estimator
                # for LASSO/p: use a vector with all 0
                # for pLASSO: use prior estimator
                # for LASSO-A: use LASSO estimator
                # for p-A: use prior estimator
                # for pLASSO-A: use pLASSO estimator
        # penalty.factor: a vector of indicators of which predictors are penalized (yes/no: 1/0)
        # cv.group: the index of cross validation groups
        # is.refit: is the estimator refitted (yes/no: 1/0)?
### Output:
        # the optimal estimator with cross validation
        
find.beta <- function(x, y, eta.seq, tau.seq, beta.prior, penalty.factor, cv.group, is.refit)
{
        n.group <- length(unique(cv.group))
        y.prior <- cbind(1, x) %*% beta.prior
        log.liklhd.eta.tau <- matrix(0, length(eta.seq), length(tau.seq))
        betas.eta.tau <- array(0, dim = c(ncol(x) + 1, length(eta.seq), length(tau.seq)))
        
        for(k in 1 : length(eta.seq)) 
        {
                eta <- eta.seq[k]
                y.tilde <- (y + eta * y.prior)/(1 + eta)
                for(l in 1 : length(tau.seq))
                {
                        tau <- tau.seq[l]
                        Winv <- 1 + tau * abs(beta.prior[-1])
                        x.Winv <- x * (matrix(1, nrow = nrow(x)) %*% Winv)
                        fits <- glmnet(x.Winv, y.tilde, standardize = FALSE, penalty.factor = penalty.factor)
                        
                        lambda.seq <- fits$lambda
                        log.liklhd <- matrix(0, length(lambda.seq), n.group)
                        a0s <- fits$a0
                        betas <- fits$beta
                
                        for(j in 1 : n.group)
                        {
                                sel <- (cv.group != j)
                                x.Winv.train <- x.Winv[sel, ]
                                y.tilde.train <- y.tilde[sel]
                                x.tune <- x[!sel, ]
                                y.tune <- y[!sel]
                                fits <- glmnet(x.Winv.train, y.tilde.train, standardize = FALSE, penalty.factor = penalty.factor)
                                
                                for(i in 1 : length(lambda.seq))
                                {
                                        beta <- coef(fits, s = lambda.seq[i])
                                        beta.refit <- rep(0, length(beta))
                                        if(is.refit) {
                                                temp <- update.beta(x.Winv.train, y.tilde.train, beta)
                                        }
                                        else {
                                                temp <- beta
                                        }
                                        beta.refit[1] <- temp[1]
                                        beta.refit[-1] <- Winv * temp[-1]
                                        log.liklhd[i, j] <- MSR.function(x.tune, y.tune, beta.refit) * (-1)
                                }
                        }
                        
                        mean.log.liklhd <- apply(log.liklhd, 1, mean)
                        se.log.liklhd <- apply(log.liklhd, 1, sd)/sqrt(n.group)
                        pos.max <- which.max(mean.log.liklhd)
                        # cat("pos.max = ", pos.max, "\n")
                        # one standard error rule
                        temp <- which(mean.log.liklhd >= mean.log.liklhd[pos.max] - se.log.liklhd[pos.max])
                        pos.1se <- temp[1]
                        # cat("pos.1se = ", pos.1se, "\n")

                        log.liklhd.eta.tau[k, l] <- mean.log.liklhd[pos.1se]
                        betas.eta.tau[, k, l] <- c(a0s[pos.1se], betas[, pos.1se])
                }
        }
        
        pos.max.eta.tau <- which(log.liklhd.eta.tau == max(log.liklhd.eta.tau), arr.ind = TRUE)
        if(nrow(pos.max.eta.tau) > 1) {
                cat("pos.max.eta.tau has more than one row", "\n")
                pos.max.eta.tau <- pos.max.eta.tau[nrow(pos.max.eta.tau), ]
        }
        eta <- eta.seq[pos.max.eta.tau[1]]
        tau <- tau.seq[pos.max.eta.tau[2]]
        
        if(length(eta.seq) > 1 && length(tau.seq) == 1) {
                cat("plasso: selected eta is ", eta, "\n")
        }
        if(length(eta.seq) == 1 && length(tau.seq) > 1) {
                cat("alasso: selected tau is ", tau, "\n")
        }
        if(length(eta.seq) > 1 && length(tau.seq) > 1) {
                cat("plasso + alasso: selected eta is ", eta, " and selected tau is ", tau, "\n")
        }
        
        y.tilde <- (y + eta * y.prior)/(1 + eta)
        Winv <- 1 + tau * abs(beta.prior[-1])
        x.Winv <- x * (matrix(1, nrow = nrow(x)) %*% Winv)
        beta <- betas.eta.tau[, pos.max.eta.tau[1], pos.max.eta.tau[2]]
        beta.refit <- rep(0, length(beta))
        if(is.refit) {
                temp <- update.beta(x.Winv, y.tilde, beta)
        }
        else {
                temp <- beta
        }
        beta.refit[1] <- temp[1]
        beta.refit[-1] <- Winv * temp[-1]
        
        return(beta.refit = beta.refit)
}