require(lars) require(quantreg) ################################################################### # Regularized Rank Regression. # Reference: Leng, C. (2010). Variable Selection and Coefficient # Estimation via Regularized Rank Regression. # Statistica Sinica, 20, 167-181. # input: observations x and y # output: # beta: the whole solution path of beta # beta.sic: the beta chosen by sic ################################################################### rrr <- function(x,y) { n <- nrow(x) X <- do.call('rbind',sapply(2:n, ff, x)) Y <- do.call('c', sapply(2:n, ff1, y)) fit <- rq.fit.fn(X,Y) beta <- lars.ls(t(x)%*%x, coef(fit)) res <- drop(y)-x%*%t(beta) tmp <- apply(res,2,rank)-(n+1)/2 II <- x%*%solve(t(x)%*%x,t(x)) T.l <- 3*diag(t(tmp)%*%II%*%tmp)/n/n SIC <- T.l+ apply(abs(beta)>1e-10,1,sum)*log(log(dim(x)[1])) s.sic <- sort(SIC, index=T) beta.sic <- beta[s.sic$ix[1],] list(beta=beta, beta.sic=beta.sic) } ## Obj: (b-b0)' Sigma0 (b-b0)+ adaptive lasso penalty ## Input: Sigma0 and b0 ## Output: beta, the whole solution path lars.ls <- function (Sigma0, b0, gamma=1, type = c("lasso", "lar"), eps = .Machine$double.eps, max.steps) { type <- match.arg(type) TYPE <- switch(type, lasso = "LASSO", lar = "LAR") n1 <- dim(Sigma0)[1] Sigma <- Sigma0 b <- b0 Sigma <- diag(abs(b)^gamma)%*%Sigma%*%diag(abs(b)^gamma) b <- b/abs(b)^gamma nm <- dim(Sigma) m <- nm[2] im <- inactive <- seq(m) Cvec <- drop(t(b)%*%Sigma) ssy <- sum(Cvec*b) if (missing(max.steps)) max.steps <- 8 * m beta <- matrix(0, max.steps + 1, m) Gamrat <- NULL arc.length <- NULL R2 <- 1 RSS <- ssy first.in <- integer(m) active <- NULL actions <- as.list(seq(max.steps)) drops <- FALSE Sign <- NULL R <- NULL k <- 0 ignores <- NULL while ((k < max.steps) & (length(active) < m)) { action <- NULL k <- k + 1 C <- Cvec[inactive] Cmax <- max(abs(C)) if (!any(drops)) { new <- abs(C) >= Cmax - eps C <- C[!new] new <- inactive[new] for (inew in new) { R <- updateR(Sigma[inew, inew], R, drop(Sigma[inew, active]), Gram = TRUE,eps=eps) if(attr(R, "rank") == length(active)) { ##singularity; back out nR <- seq(length(active)) R <- R[nR, nR, drop = FALSE] attr(R, "rank") <- length(active) ignores <- c(ignores, inew) action <- c(action, - inew) } else { if(first.in[inew] == 0) first.in[inew] <- k active <- c(active, inew) Sign <- c(Sign, sign(Cvec[inew])) action <- c(action, inew) } } } else action <- -dropid Gi1 <- backsolve(R, backsolvet(R, Sign)) dropouts <- NULL A <- 1/sqrt(sum(Gi1 * Sign)) w <- A * Gi1 if (length(active) >= m) { gamhat <- Cmax/A } else { a <- drop(w %*% Sigma[active, -c(active,ignores), drop = FALSE]) gam <- c((Cmax - C)/(A - a), (Cmax + C)/(A + a)) gamhat <- min(gam[gam > eps], Cmax/A) } if (type == "lasso") { dropid <- NULL b1 <- beta[k, active] z1 <- -b1/w zmin <- min(z1[z1 > eps], gamhat) # cat('zmin ',zmin, ' gamhat ',gamhat,'\n') if (zmin < gamhat) { gamhat <- zmin drops <- z1 == zmin } else drops <- FALSE } beta[k + 1, ] <- beta[k, ] beta[k + 1, active] <- beta[k + 1, active] + gamhat * w Cvec <- Cvec - gamhat * Sigma[, active, drop = FALSE] %*% w Gamrat <- c(Gamrat, gamhat/(Cmax/A)) arc.length <- c(arc.length, gamhat) if (type == "lasso" && any(drops)) { dropid <- seq(drops)[drops] for (id in rev(dropid)) { R <- downdateR(R,id) } dropid <- active[drops] beta[k + 1, dropid] <- 0 active <- active[!drops] Sign <- Sign[!drops] } actions[[k]] <- action inactive <- im[-c(active)] } beta <- beta[seq(k + 1), ] #beta <- t(abs(b0)^gamma*t(beta)) beta <- beta%*%diag(abs(b0)) beta } ## pairwise difference operator ff <- function(i, mat) { mmt <- t(mat) z <- mmt[,i-1] - mmt[,i:nrow(mat),drop=F] t(z) } ff1 <- function(i, vec) vec[i-1]-vec[i:length(vec)]